From 37e73567377b742a3e701dc00c81e1d5e8d8adf9 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 13 Jul 2017 08:43:04 -0600 Subject: [PATCH 1/9] Added ocn_import_export.F90 --- config_src/mct_driver/ocn_comp_mct.F90 | 243 ++++++++++++++++---- config_src/mct_driver/ocn_import_export.F90 | 140 +++++++++++ 2 files changed, 341 insertions(+), 42 deletions(-) create mode 100644 config_src/mct_driver/ocn_import_export.F90 diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index b3c760c0fd..8e4dee42c8 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -6,26 +6,32 @@ module ocn_comp_mct ! !INTERFACE: ! !DESCRIPTION: -! This is the main driver for the MOM6 in CIME +! This is the main driver for MOM6 in CIME ! ! !REVISION HISTORY: ! ! !USES: - use esmf - use seq_cdata_mod - use mct_mod - use seq_infodata_mod, only: seq_infodata_type, & - seq_infodata_GetData, & - seq_infodata_start_type_start, & - seq_infodata_start_type_cont, & - seq_infodata_start_type_brnch + use esmf + use seq_cdata_mod + use mct_mod + use seq_flds_mod, only: seq_flds_x2o_fields, & + seq_flds_o2x_fields + use seq_infodata_mod, only: seq_infodata_type, & + seq_infodata_GetData, & + seq_infodata_start_type_start, & + seq_infodata_start_type_cont, & + seq_infodata_start_type_brnch + use seq_comm_mct, only: seq_comm_name, seq_comm_inst, seq_comm_suffix + use perf_mod, only: t_startf, t_stopf + ! From MOM6 - use ocean_model_mod, only: ocean_state_type, ocean_public_type - use ocean_model_mod, only: ocean_model_init - use MOM_time_manager, only: time_type, set_date, set_calendar_type, NOLEAP - use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here - use coupler_indices, only: coupler_indices_init + use ocean_model_mod, only: ocean_state_type, ocean_public_type + use ocean_model_mod, only: ocean_model_init + use MOM_time_manager, only: time_type, set_date, set_calendar_type, NOLEAP + use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here + use coupler_indices, only: coupler_indices_init + use ocn_import_export, only: SBUFF_SUM ! ! !PUBLIC MEMBER FUNCTIONS: @@ -43,13 +49,14 @@ module ocn_comp_mct ! !EOP ! !PRIVATE MODULE FUNCTIONS: + private :: ocn_SetGSMap_mct + private :: ocn_domain_mct -! ! !PRIVATE MODULE VARIABLES - type(ocean_state_type), pointer :: ocn_state => NULL() ! Private state of ocean - type(ocean_public_type), pointer :: ocn_surface => NULL() ! Public surface state of ocean + type(ocean_state_type), pointer :: ocn_state => NULL() ! Private state of ocean + type(ocean_public_type), pointer :: ocn_surface => NULL() ! Public surface state of ocean - type(seq_infodata_type), pointer :: & + type(seq_infodata_type), pointer :: & infodata !======================================================================= @@ -90,19 +97,69 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) character(len=384) :: runid character(len=384) :: runtype character(len=32) :: starttype ! infodata start type - integer :: mpicom + integer :: mpicom_ocn integer :: npes, pe0 integer :: i + integer :: lsize, nsend, nrecv + + ! mct variables (these are local for now) + integer :: MOM_MCT_ID + type(mct_gsMap), pointer :: MOM_MCT_gsMap ! 2d, points to cdata + type(mct_gGrid), pointer :: MOM_MCT_dom ! 2d, points to cdata + type(mct_gsMap) :: MOM_MCT_gsMap3d ! for 3d streams, local + type(mct_gGrid) :: MOM_MCT_dom3d ! for 3d streams, local + + ! instance control vars (these are local for now) + integer(kind=4) :: inst_index + character(len=16) :: inst_name + character(len=16) :: inst_suffix + + !!!DANGER!!!: change the following vars with the corresponding MOM6 vars + integer :: km=62 ! number of vertical levels +!----------------------------------------------------------------------- + + ! set (actually, get from mct) the cdata pointers: + call seq_cdata_setptrs(cdata_o, id=MOM_MCT_ID, mpicom=mpicom_ocn, infodata=infodata) - ! Initialize MOM6 - mpicom = cdata_o%mpicom - call MOM_infra_init(mpicom) + !--------------------------------------------------------------------- + ! Initialize the model run + !--------------------------------------------------------------------- + + call coupler_indices_init() + + call seq_infodata_GetData( infodata, case_name=runid ) + + call seq_infodata_GetData( infodata, start_type=starttype) + + if ( trim(starttype) == trim(seq_infodata_start_type_start)) then + runtype = "initial" + else if (trim(starttype) == trim(seq_infodata_start_type_cont) ) then + runtype = "continue" + else if (trim(starttype) == trim(seq_infodata_start_type_brnch)) then + runtype = "branch" + else + write(*,*) 'ocn_comp_mct ERROR: unknown starttype' + call exit(0) + end if + + ! instance control + + inst_name = seq_comm_name(MOM_MCT_ID) + inst_index = seq_comm_inst(MOM_MCT_ID) + inst_suffix = seq_comm_suffix(MOM_MCT_ID) + + !--------------------------------------------------------------------- + ! Initialize MOM6 + !--------------------------------------------------------------------- + + call t_startf('MOM_init') + + call MOM_infra_init(mpicom_ocn) call ESMF_ClockGet(EClock, currTime=current_time, rc=rc) call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) - ! we need to confirm this: - call set_calendar_type(NOLEAP) + call set_calendar_type(NOLEAP) !TODO: confirm this time_init = set_date(year, month, day, hour, minute, seconds, err_msg=errMsg) time_in = set_date(year, month, day, hour, minute, seconds, err_msg=errMsg) @@ -115,30 +172,59 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) allocate(ocn_surface%pelist(npes)) ocn_surface%pelist(:) = (/(i,i=pe0,pe0+npes)/) - ! initialize the model run + ! initialize the MOM6 model call ocean_model_init(ocn_surface, ocn_state, time_init, time_in) - ! set infodata, a cdata pointer - call seq_cdata_setptrs(cdata_o, infodata=infodata) + call t_stopf('MOM_init') + + !--------------------------------------------------------------------- + ! Initialize MCT attribute vectors and indices + !--------------------------------------------------------------------- + + call t_startf('MOM_mct_init') + + ! Set mct global seg maps: + + call ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, MOM_MCT_GSMap, MOM_MCT_GSMap3d) + lsize = mct_gsMap_lsize(MOM_MCT_gsmap, mpicom_ocn) + + ! Initialize mct ocn domain (needs ocn initialization info) + + call ocn_domain_mct(lsize, MOM_MCT_gsmap, MOM_MCT_dom) + call ocn_domain_mct(lsize*km, MOM_MCT_gsmap3d, MOM_MCT_dom3d) + + ! Inialize mct attribute vectors + + ! Initialize the mct attribute vector x2o_o, given Attribute list and length: + call mct_aVect_init(x2o_o, rList=seq_flds_x2o_fields, lsize=lsize) + ! set the mct attribute vector x2o_o to zero: + call mct_aVect_zero(x2o_o) + + ! Initialize the mct attribute vector o2x_o, given Attribute list and length: + call mct_aVect_init(o2x_o, rList=seq_flds_o2x_fields, lsize=lsize) + ! set the mct attribute vector o2x_o to zero: + call mct_aVect_zero(o2x_o) + + nsend = mct_avect_nRattr(o2x_o) + nrecv = mct_avect_nRattr(x2o_o) + !allocate (SBUFF_SUM(nx_block,ny_block,max_blocks_clinic,nsend)) + + + + + + + + - ! initialize coupler indices - call coupler_indices_init() - ! get runid and starttype: - call seq_infodata_GetData( infodata, case_name=runid ) - call seq_infodata_GetData( infodata, start_type=starttype) - if ( trim(starttype) == trim(seq_infodata_start_type_start)) then - runtype = "initial" - else if (trim(starttype) == trim(seq_infodata_start_type_cont) ) then - runtype = "continue" - else if (trim(starttype) == trim(seq_infodata_start_type_brnch)) then - runtype = "branch" - else - write(*,*) 'ocn_comp_mct ERROR: unknown starttype' - call exit(0) - end if + + + + call t_stopf('MOM_mct_init') + !----------------------------------------------------------------------- @@ -209,6 +295,79 @@ subroutine ocn_final_mct( EClock, cdata_o, x2o_o, o2x_o) end subroutine ocn_final_mct +!*********************************************************************** +!BOP +!IROUTINE: ocn_SetGSMap_mct +! !INTERFACE: + + subroutine ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, gsMap_ocn, gsMap3d_ocn) + +! !DESCRIPTION: +! This routine mct global seg maps for the MOM decomposition +! +! !REVISION HISTORY: +! same as module + +! !INPUT/OUTPUT PARAMETERS: + + implicit none + integer , intent(in) :: mpicom_ocn + integer , intent(in) :: MOM_MCT_ID + type(mct_gsMap), intent(inout) :: gsMap_ocn + type(mct_gsMap), intent(inout) :: gsMap3d_ocn + +!EOP +!BOC +!----------------------------------------------------------------------- +! +! local variables +! +!----------------------------------------------------------------------- + + + +!----------------------------------------------------------------------- +!EOC + + end subroutine ocn_SetGSMap_mct + + +!*********************************************************************** +!BOP +! !IROUTINE: ocn_domain_mct +! !INTERFACE: + + subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn) + +! !DESCRIPTION: +! This routine mct global seg maps for the pop decomposition +! +! !REVISION HISTORY: +! same as module +! +! !INPUT/OUTPUT PARAMETERS: + + implicit none + integer , intent(in) :: lsize + type(mct_gsMap), intent(in) :: gsMap_ocn + type(mct_ggrid), intent(inout) :: dom_ocn + +!EOP +!BOC +!----------------------------------------------------------------------- +! +! local variables +! +!----------------------------------------------------------------------- + + + +!----------------------------------------------------------------------- +!EOC + + end subroutine ocn_domain_mct + + end module ocn_comp_mct !||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| diff --git a/config_src/mct_driver/ocn_import_export.F90 b/config_src/mct_driver/ocn_import_export.F90 new file mode 100644 index 0000000000..6fd16835fb --- /dev/null +++ b/config_src/mct_driver/ocn_import_export.F90 @@ -0,0 +1,140 @@ +module ocn_import_export + + implicit none + public + save + + ! accumulated sum of send buffer quantities for averaging before being sent + !real (r8), dimension(:,:,:,:), allocatable :: SBUFF_SUM + !real (r8) :: tlast_coupled + + !TODO: update the types of following vars + double precision, dimension(:,:,:,:), allocatable :: SBUFF_SUM + double precision :: tlast_coupled +contains + +!*********************************************************************** +!BOP +! !IROUTINE: ocn_import +! !INTERFACE: + + subroutine ocn_import(x2o, ldiag_cpl, errorCode) + +! !DESCRIPTION: +!----------------------------------------------------------------------- +! This routine receives message from cpl7 driver +! +! The following fields are always received from the coupler: +! +! o taux -- zonal wind stress (taux) (W/m2 ) +! o tauy -- meridonal wind stress (tauy) (W/m2 ) +! o snow -- water flux due to snow (kg/m2/s) +! o rain -- water flux due to rain (kg/m2/s) +! o evap -- evaporation flux (kg/m2/s) +! o meltw -- snow melt flux (kg/m2/s) +! o salt -- salt (kg(salt)/m2/s) +! o swnet -- net short-wave heat flux (W/m2 ) +! o sen -- sensible heat flux (W/m2 ) +! o lwup -- longwave radiation (up) (W/m2 ) +! o lwdn -- longwave radiation (down) (W/m2 ) +! o melth -- heat flux from snow&ice melt (W/m2 ) +! o ifrac -- ice fraction +! o rofl -- river runoff flux (kg/m2/s) +! o rofi -- ice runoff flux (kg/m2/s) +! +! The following fields are sometimes received from the coupler, +! depending on model options: +! +! o pslv -- sea-level pressure (Pa) +! o duu10n -- 10m wind speed squared (m^2/s^2) +! o co2prog-- bottom atm level prognostic co2 +! o co2diag-- bottom atm level diagnostic co2 +! +!----------------------------------------------------------------------- +! +! !REVISION HISTORY: +! same as module + +! !INPUT/OUTPUT PARAMETERS: + + !real(r8) , intent(inout) :: x2o(:,:) + !logical (log_kind) , intent(in) :: ldiag_cpl + !integer (POP_i4) , intent(out) :: errorCode ! returned error code + + !TODO: update the types of following params + double precision, intent(inout) :: x2o(:,:) + logical, intent(in) :: ldiag_cpl + integer, intent(out) :: errorCode ! returned error code + +!EOP +!BOC +!----------------------------------------------------------------------- +! +! local variables +! +!----------------------------------------------------------------------- + + + + + + + + + + +!----------------------------------------------------------------------- +!EOC + + end subroutine ocn_import + +!*********************************************************************** +!BOP +! !IROUTINE: ocn_export_mct +! !INTERFACE: + + subroutine ocn_export(o2x, ldiag_cpl, errorCode) + +! !DESCRIPTION: +! This routine calls the routines necessary to send MOM6 fields to +! the CCSM cpl7 driver +! +! !REVISION HISTORY: +! same as module +! +! !INPUT/OUTPUT PARAMETERS: + + !real(r8) , intent(inout) :: o2x(:,:) + !logical (log_kind) , intent(in) :: ldiag_cpl + !integer (POP_i4) , intent(out) :: errorCode ! returned error code + + !TODO: update the types of following params + double precision, intent(inout) :: o2x(:,:) + logical, intent(in) :: ldiag_cpl + integer, intent(out) :: errorCode ! returned error code + +!EOP +!BOC +!----------------------------------------------------------------------- +! +! local variables +! +!----------------------------------------------------------------------- + + + + + + + + +!----------------------------------------------------------------------- +!EOC + + end subroutine ocn_export + +!*********************************************************************** + + +end module ocn_import_export + From 57e60868208767bc4720974ce5b8e30451957698 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 13 Jul 2017 17:35:33 -0600 Subject: [PATCH 2/9] Added methods to return domain/grid shape - get_global_shape() returns niglobal, njglobal from domain type. - get_global_grid_shape() returns niglobal, njglobal by calling get_global_shape(). - This avoids exposing members inside opaque types otherwise needed for initializing under MCT. --- src/core/MOM_grid.F90 | 13 ++++++++++++- src/framework/MOM_domains.F90 | 13 ++++++++++++- 2 files changed, 24 insertions(+), 2 deletions(-) diff --git a/src/core/MOM_grid.F90 b/src/core/MOM_grid.F90 index 6d5597ce4e..0f01700c01 100644 --- a/src/core/MOM_grid.F90 +++ b/src/core/MOM_grid.F90 @@ -5,6 +5,7 @@ module MOM_grid use MOM_hor_index, only : hor_index_type, hor_index_init use MOM_domains, only : MOM_domain_type, get_domain_extent, compute_block_extent +use MOM_domains, only : get_global_shape use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL use MOM_file_parser, only : get_param, log_param, log_version, param_file_type @@ -13,7 +14,7 @@ module MOM_grid #include public MOM_grid_init, MOM_grid_end, set_derived_metrics, set_first_direction -public isPointInCell, hor_index_type +public isPointInCell, hor_index_type, get_global_grid_size !> Ocean grid type. See mom_grid for details. type, public :: ocean_grid_type @@ -443,6 +444,16 @@ subroutine set_first_direction(G, y_first) G%first_direction = y_first end subroutine set_first_direction +!> Return global shape of horizontal grid +subroutine get_global_grid_size(G, niglobal, njglobal) + type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type + integer, intent(out) :: niglobal !< i-index global size of grid + integer, intent(out) :: njglobal !< j-index global size of grid + + call get_global_shape(G%domain, niglobal, njglobal) + +end subroutine get_global_grid_size + !> Allocate memory used by the ocean_grid_type and related structures. subroutine allocate_metrics(G) type(ocean_grid_type), intent(inout) :: G !< The horizontal grid type diff --git a/src/framework/MOM_domains.F90 b/src/framework/MOM_domains.F90 index ca1dc281a2..5e872d0a72 100644 --- a/src/framework/MOM_domains.F90 +++ b/src/framework/MOM_domains.F90 @@ -60,7 +60,7 @@ module MOM_domains public :: To_East, To_West, To_North, To_South, To_All, Omit_Corners public :: create_group_pass, do_group_pass, group_pass_type public :: start_group_pass, complete_group_pass -public :: compute_block_extent +public :: compute_block_extent, get_global_shape interface pass_var module procedure pass_var_3d, pass_var_2d @@ -1599,4 +1599,15 @@ subroutine get_domain_extent(Domain, isc, iec, jsc, jec, isd, ied, jsd, jed, & end subroutine get_domain_extent +!> Returns the global shape of h-point arrays +subroutine get_global_shape(domain, niglobal, njglobal) + type(MOM_domain_type), intent(in) :: domain !< MOM domain + integer, intent(out) :: niglobal !< i-index global size of h-point arrays + integer, intent(out) :: njglobal !< j-index global size of h-point arrays + + niglobal = domain%niglobal + njglobal = domain%njglobal + +end subroutine get_global_shape + end module MOM_domains From 5df3bc3ba707fdf9d1f5b214997479eeecfbba79 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 13 Jul 2017 17:37:09 -0600 Subject: [PATCH 3/9] Added method to point ponters to members of ocean_state_type - To initialize within the MCT coupler we need access to members of ocean_state_type which are private. This method allows us to have local pointers to those members. --- config_src/coupled_driver/ocean_model_MOM.F90 | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 0f0f1f2622..dfcdded6f8 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -89,7 +89,9 @@ module ocean_model_mod public ocean_model_restart public ice_ocn_bnd_type_chksum public ocean_public_type_chksum -public ocean_model_data_get +public ocean_model_data_get +public get_state_pointers + interface ocean_model_data_get module procedure ocean_model_data1D_get module procedure ocean_model_data2D_get @@ -1083,4 +1085,13 @@ subroutine ocean_public_type_chksum(id, timestep, ocn) 100 FORMAT(" CHECKSUM::",A20," = ",Z20) end subroutine ocean_public_type_chksum +!> Returns pointers to objects within ocean_state_type +subroutine get_state_pointers(OS, grid) + type(ocean_state_type), pointer :: OS !< Ocean state type + type(ocean_grid_type), optional, pointer :: grid !< Ocean grid + + if (present(grid)) grid => OS%grid + +end subroutine get_state_pointers + end module ocean_model_mod From aba7dc2a3a1c52233ec5c7382858b0382395ba75 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 13 Jul 2017 17:42:00 -0600 Subject: [PATCH 4/9] Completed ocn_init_mct --- config_src/mct_driver/ocn_comp_mct.F90 | 173 ++++++++++++++++--------- 1 file changed, 111 insertions(+), 62 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 8e4dee42c8..5aa1cd914b 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -20,18 +20,22 @@ module ocn_comp_mct seq_infodata_GetData, & seq_infodata_start_type_start, & seq_infodata_start_type_cont, & - seq_infodata_start_type_brnch + seq_infodata_start_type_brnch, & + seq_infodata_PutData use seq_comm_mct, only: seq_comm_name, seq_comm_inst, seq_comm_suffix + use seq_timemgr_mod, only: seq_timemgr_EClockGetData use perf_mod, only: t_startf, t_stopf - - + + ! From MOM6 use ocean_model_mod, only: ocean_state_type, ocean_public_type - use ocean_model_mod, only: ocean_model_init - use MOM_time_manager, only: time_type, set_date, set_calendar_type, NOLEAP + use ocean_model_mod, only: ocean_model_init, get_state_pointers use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here + use MOM_grid, only: ocean_grid_type, get_global_grid_size + use MOM_error_handler, only: MOM_error, FATAL, is_root_pe + use MOM_time_manager, only: time_type, set_date, set_calendar_type, NOLEAP use coupler_indices, only: coupler_indices_init - use ocn_import_export, only: SBUFF_SUM + use ocn_import_export, only: SBUFF_SUM, ocn_Export, mom_sum_buffer ! ! !PUBLIC MEMBER FUNCTIONS: @@ -72,7 +76,7 @@ module ocn_comp_mct subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! ! !DESCRIPTION: -! Initialize POP +! Initialize POP ! ! !INPUT/OUTPUT PARAMETERS: @@ -99,31 +103,50 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) character(len=32) :: starttype ! infodata start type integer :: mpicom_ocn integer :: npes, pe0 - integer :: i + integer :: i, errorCode integer :: lsize, nsend, nrecv + logical :: ldiag_cpl = .false. ! mct variables (these are local for now) integer :: MOM_MCT_ID - type(mct_gsMap), pointer :: MOM_MCT_gsMap ! 2d, points to cdata - type(mct_gGrid), pointer :: MOM_MCT_dom ! 2d, points to cdata + type(mct_gsMap), pointer :: MOM_MCT_gsMap => NULL() ! 2d, points to cdata + type(mct_gGrid), pointer :: MOM_MCT_dom => NULL() ! 2d, points to cdata type(mct_gsMap) :: MOM_MCT_gsMap3d ! for 3d streams, local type(mct_gGrid) :: MOM_MCT_dom3d ! for 3d streams, local + ! time management + integer :: ocn_cpl_dt + real (kind=8) :: mom_cpl_dt + real (kind=8), parameter :: & + seconds_in_minute = 60.0d0, & + seconds_in_hour = 3600.0d0, & + seconds_in_day = 86400.0d0, & + minutes_in_hour = 60.0d0 + + ! instance control vars (these are local for now) integer(kind=4) :: inst_index character(len=16) :: inst_name character(len=16) :: inst_suffix !!!DANGER!!!: change the following vars with the corresponding MOM6 vars - integer :: km=62 ! number of vertical levels + integer :: km=62 ! number of vertical levels + integer :: nx_block=0, ny_block=0 ! size of block domain in x,y dir including ghost cells + integer :: nx_global, ny_global! size of block domain in x,y dir including ghost cells + integer :: max_blocks_clinic=0 !max number of blocks per processor in each distribution + integer :: ncouple_per_day = 24 + logical :: lsend_precip_fact ! if T,send precip_fact to cpl for use in fw balance + ! (partially-coupled option) + !----------------------------------------------------------------------- ! set (actually, get from mct) the cdata pointers: - call seq_cdata_setptrs(cdata_o, id=MOM_MCT_ID, mpicom=mpicom_ocn, infodata=infodata) + call seq_cdata_setptrs(cdata_o, id=MOM_MCT_ID, mpicom=mpicom_ocn, & + gsMap=MOM_MCT_gsMap, dom=MOM_MCT_dom, infodata=infodata) !--------------------------------------------------------------------- - ! Initialize the model run + ! Initialize the model run !--------------------------------------------------------------------- call coupler_indices_init() @@ -150,14 +173,14 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) inst_suffix = seq_comm_suffix(MOM_MCT_ID) !--------------------------------------------------------------------- - ! Initialize MOM6 + ! Initialize MOM6 !--------------------------------------------------------------------- call t_startf('MOM_init') call MOM_infra_init(mpicom_ocn) - call ESMF_ClockGet(EClock, currTime=current_time, rc=rc) + call ESMF_ClockGet(EClock, currTime=current_time, rc=rc) call ESMF_TimeGet(current_time, yy=year, mm=month, dd=day, h=hour, m=minute, s=seconds, rc=rc) call set_calendar_type(NOLEAP) !TODO: confirm this @@ -178,15 +201,15 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) call t_stopf('MOM_init') !--------------------------------------------------------------------- - ! Initialize MCT attribute vectors and indices + ! Initialize MCT attribute vectors and indices !--------------------------------------------------------------------- call t_startf('MOM_mct_init') ! Set mct global seg maps: - call ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, MOM_MCT_GSMap, MOM_MCT_GSMap3d) - lsize = mct_gsMap_lsize(MOM_MCT_gsmap, mpicom_ocn) + call ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, MOM_MCT_GSMap, MOM_MCT_GSMap3d) + lsize = mct_gsMap_lsize(MOM_MCT_gsmap, mpicom_ocn) ! Initialize mct ocn domain (needs ocn initialization info) @@ -194,36 +217,51 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) call ocn_domain_mct(lsize*km, MOM_MCT_gsmap3d, MOM_MCT_dom3d) ! Inialize mct attribute vectors - + ! Initialize the mct attribute vector x2o_o, given Attribute list and length: call mct_aVect_init(x2o_o, rList=seq_flds_x2o_fields, lsize=lsize) ! set the mct attribute vector x2o_o to zero: call mct_aVect_zero(x2o_o) - + ! Initialize the mct attribute vector o2x_o, given Attribute list and length: call mct_aVect_init(o2x_o, rList=seq_flds_o2x_fields, lsize=lsize) ! set the mct attribute vector o2x_o to zero: call mct_aVect_zero(o2x_o) + ! allocate send buffer nsend = mct_avect_nRattr(o2x_o) - nrecv = mct_avect_nRattr(x2o_o) - !allocate (SBUFF_SUM(nx_block,ny_block,max_blocks_clinic,nsend)) - - - - + nrecv = mct_avect_nRattr(x2o_o) + allocate (SBUFF_SUM(nx_block,ny_block,max_blocks_clinic,nsend)) + ! initialize necessary coupling info + call seq_timemgr_EClockGetData(EClock, dtime=ocn_cpl_dt) + mom_cpl_dt = seconds_in_day / ncouple_per_day + if (mom_cpl_dt /= ocn_cpl_dt) then + write(*,*) 'ERROR pop_cpl_dt and ocn_cpl_dt must be identical' + call exit(0) + end if + ! send initial state to driver + !TODO: + ! if ( lsend_precip_fact ) then + ! call seq_infodata_PutData( infodata, precip_fact=precip_fact) + ! end if + call mom_sum_buffer + call ocn_export(o2x_o%rattr, ldiag_cpl, errorCode) + call t_stopf('MOM_mct_init') - call t_stopf('MOM_mct_init') + call seq_infodata_PutData( infodata, & + ocn_nx = nx_global , ocn_ny = ny_global) + call seq_infodata_PutData( infodata, & + ocn_prognostic=.true., ocnrof_prognostic=.true.) @@ -295,41 +333,52 @@ subroutine ocn_final_mct( EClock, cdata_o, x2o_o, o2x_o) end subroutine ocn_final_mct -!*********************************************************************** -!BOP -!IROUTINE: ocn_SetGSMap_mct -! !INTERFACE: - - subroutine ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, gsMap_ocn, gsMap3d_ocn) - -! !DESCRIPTION: -! This routine mct global seg maps for the MOM decomposition -! -! !REVISION HISTORY: -! same as module - -! !INPUT/OUTPUT PARAMETERS: - - implicit none - integer , intent(in) :: mpicom_ocn - integer , intent(in) :: MOM_MCT_ID - type(mct_gsMap), intent(inout) :: gsMap_ocn - type(mct_gsMap), intent(inout) :: gsMap3d_ocn - -!EOP -!BOC -!----------------------------------------------------------------------- -! -! local variables -! -!----------------------------------------------------------------------- - - - -!----------------------------------------------------------------------- -!EOC - - end subroutine ocn_SetGSMap_mct +!> This routine mct global seg maps for the MOM decomposition +!! +!! \todo Find out if we should only provide indirect indexing for ocean points and not land. +subroutine ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, gsMap_ocn, gsMap3d_ocn) + integer, intent(in) :: mpicom_ocn !< MPI communicator + integer, intent(in) :: MOM_MCT_ID !< MCT component ID + type(mct_gsMap), intent(inout) :: gsMap_ocn !< MCT global segment map for 2d data + type(mct_gsMap), intent(inout) :: gsMap3d_ocn !< MCT global segment map for 3d data + ! Local variables + integer :: lsize ! Local size of indirect indexing array + integer :: i, j, k ! Local indices + integer :: ni, nj ! Declared sizes of h-point arrays + integer :: ig, jg ! Global indices + type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure + integer, allocatable :: gindex(:) ! Indirect indices + + call get_state_pointers(ocn_state, grid=grid) + if (.not. associated(grid)) call MOM_error(FATAL, 'ocn_comp_mct.F90, ocn_SetGSMap_mct():' // & + 'grid returned from get_state_pointers() was not associated!') + + ! Size of computational domain + lsize = ( grid%iec - grid%isc + 1 ) * ( grid%jec - grid%jsc + 1 ) + + ! Size of global domain + call get_global_grid_size(grid, ni, nj) + + ! Create indirect indices for the computational domain + allocate( gindex( lsize ) ) + + ! Set indirect indices in gindex + k = 0 + do j = grid%jsc, grid%jec + jg = j - grid%jdg_offset ! TODO: check this calculation + do i = grid%isc, grid%iec + ig = i - grid%idg_offset ! TODO: check this calculation + k = k + 1 ! Increment position within gindex + gindex(k) = ni * ( jg - 1 ) + ig + enddo + enddo + + ! Tell MCT how to indirectly index into the 2d buffer + call mct_gsMap_init( gsMap_ocn, gindex, mpicom_ocn, MOM_MCT_ID, lsize, ni * nj) + + deallocate( gindex ) + +end subroutine ocn_SetGSMap_mct !*********************************************************************** From ed3ed306863252788464679bad76c852621b3eef Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 13 Jul 2017 17:42:25 -0600 Subject: [PATCH 5/9] Added ocn_import_export module --- config_src/mct_driver/ocn_import_export.F90 | 28 +++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/config_src/mct_driver/ocn_import_export.F90 b/config_src/mct_driver/ocn_import_export.F90 index 6fd16835fb..5cb025a706 100644 --- a/config_src/mct_driver/ocn_import_export.F90 +++ b/config_src/mct_driver/ocn_import_export.F90 @@ -136,5 +136,33 @@ end subroutine ocn_export !*********************************************************************** + subroutine MOM_sum_buffer + +! !DESCRIPTION: +! This routine accumulates sums for averaging fields to +! be sent to the coupler +! +! !REVISION HISTORY: +! same as module +! +!EOP +!BOC + +!----------------------------------------------------------------------- +! +! local variables +! +!----------------------------------------------------------------------- + +!----------------------------------------------------------------------- +!EOC + + end subroutine MOM_sum_buffer + +!*********************************************************************** + + + + end module ocn_import_export From 2540ed80cd9f15be8e5896891e558cecfd2a039c Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 19 Jul 2017 09:51:04 -0600 Subject: [PATCH 6/9] Filled empty subroutines in ocn_comp_mct.F90 --- config_src/mct_driver/ocn_comp_mct.F90 | 90 ++++++++++++++++++++++++-- 1 file changed, 83 insertions(+), 7 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 5aa1cd914b..9b6099fb71 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -15,7 +15,9 @@ module ocn_comp_mct use seq_cdata_mod use mct_mod use seq_flds_mod, only: seq_flds_x2o_fields, & - seq_flds_o2x_fields + seq_flds_o2x_fields, & + SEQ_FLDS_DOM_COORD, & + SEQ_FLDS_DOM_other use seq_infodata_mod, only: seq_infodata_type, & seq_infodata_GetData, & seq_infodata_start_type_start, & @@ -25,6 +27,7 @@ module ocn_comp_mct use seq_comm_mct, only: seq_comm_name, seq_comm_inst, seq_comm_suffix use seq_timemgr_mod, only: seq_timemgr_EClockGetData use perf_mod, only: t_startf, t_stopf + use shr_kind_mod, only: SHR_KIND_R8 ! From MOM6 @@ -134,7 +137,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) integer :: nx_block=0, ny_block=0 ! size of block domain in x,y dir including ghost cells integer :: nx_global, ny_global! size of block domain in x,y dir including ghost cells integer :: max_blocks_clinic=0 !max number of blocks per processor in each distribution - integer :: ncouple_per_day = 24 + integer :: ncouple_per_day = 48 logical :: lsend_precip_fact ! if T,send precip_fact to cpl for use in fw balance ! (partially-coupled option) @@ -386,7 +389,7 @@ end subroutine ocn_SetGSMap_mct ! !IROUTINE: ocn_domain_mct ! !INTERFACE: - subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn) +subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn) ! !DESCRIPTION: ! This routine mct global seg maps for the pop decomposition @@ -396,10 +399,83 @@ subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn) ! ! !INPUT/OUTPUT PARAMETERS: - implicit none - integer , intent(in) :: lsize - type(mct_gsMap), intent(in) :: gsMap_ocn - type(mct_ggrid), intent(inout) :: dom_ocn + implicit none + integer , intent(in) :: lsize + type(mct_gsMap), intent(in) :: gsMap_ocn + type(mct_ggrid), intent(inout) :: dom_ocn + +! Local Variables + integer, parameter :: SHR_REAL_R8 = selected_real_kind (12) + integer, pointer :: idata(:) + integer :: i,j,k + real(kind=SHR_REAL_R8), pointer :: data(:) + real(kind=SHR_REAL_R8) :: m2_to_rad2 + type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure + + call get_state_pointers(ocn_state, grid=grid) + + ! set coords to lat and lon, and areas to rad^2 + call mct_gGrid_init(GGrid=dom_ocn, CoordChars=trim(seq_flds_dom_coord), & + OtherChars=trim(seq_flds_dom_other), lsize=lsize ) + + call mct_avect_zero(dom_ocn%data) + allocate(data(lsize)) + + ! Determine global gridpoint number attribute, GlobGridNum, which is set automatically by MCT + k = pe_here() + call mct_gsMap_orderedPoints(gsMap_ocn, k, idata) + call mct_gGrid_importIAttr(dom_ocn,'GlobGridNum',idata,lsize) + + !initialization + data(:) = -9999.0 + call mct_gGrid_importRAttr(dom_ocn,"lat" ,data,lsize) + call mct_gGrid_importRAttr(dom_ocn,"lon" ,data,lsize) + call mct_gGrid_importRAttr(dom_ocn,"area" ,data,lsize) + call mct_gGrid_importRAttr(dom_ocn,"aream",data,lsize) + data(:) = 0.0 + call mct_gGrid_importRAttr(dom_ocn,"mask",data,lsize) + call mct_gGrid_importRAttr(dom_ocn,"frac",data,lsize) + + k = 0 + do j = grid%jsc, grid%jec + do i = grid%isc, grid%iec + k = k + 1 ! Increment position within gindex + data(k) = grid%geoLonT(i,j) + enddo + enddo + call mct_gGrid_importRattr(dom_ocn,"lon",data,lsize) + + k = 0 + do j = grid%jsc, grid%jec + do i = grid%isc, grid%iec + k = k + 1 ! Increment position within gindex + data(k) = grid%geoLatT(i,j) + enddo + enddo + call mct_gGrid_importRattr(dom_ocn,"lat",data,lsize) + + k = 0 + m2_to_rad2 = 1./grid%Rad_Earth**2 + do j = grid%jsc, grid%jec + do i = grid%isc, grid%iec + k = k + 1 ! Increment position within gindex + data(k) = grid%AreaT(i,j) * m2_to_rad2 + enddo + enddo + call mct_gGrid_importRattr(dom_ocn,"area",data,lsize) + + k = 0 + do j = grid%jsc, grid%jec + do i = grid%isc, grid%iec + k = k + 1 ! Increment position within gindex + data(k) = grid%mask2dT(i,j) + enddo + enddo + call mct_gGrid_importRattr(dom_ocn,"mask",data,lsize) + call mct_gGrid_importRattr(dom_ocn,"frac",data,lsize) + + deallocate(data) + deallocate(idata) !EOP !BOC From 970f8482ace6f6f0a54e803087157ad61432c1f5 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 19 Jul 2017 10:33:22 -0600 Subject: [PATCH 7/9] added debugging --- config_src/mct_driver/ocn_comp_mct.F90 | 29 ++++++++++++++++++++++---- 1 file changed, 25 insertions(+), 4 deletions(-) diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 9b6099fb71..bb67603f06 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -47,6 +47,7 @@ module ocn_comp_mct public :: ocn_run_mct public :: ocn_final_mct private ! By default make data private + logical, parameter :: debug=.true. ! ! ! PUBLIC DATA: @@ -109,6 +110,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) integer :: i, errorCode integer :: lsize, nsend, nrecv logical :: ldiag_cpl = .false. + integer :: ni, nj ! mct variables (these are local for now) integer :: MOM_MCT_ID @@ -133,14 +135,14 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) character(len=16) :: inst_suffix !!!DANGER!!!: change the following vars with the corresponding MOM6 vars - integer :: km=62 ! number of vertical levels + integer :: km=1 ! number of vertical levels integer :: nx_block=0, ny_block=0 ! size of block domain in x,y dir including ghost cells - integer :: nx_global, ny_global! size of block domain in x,y dir including ghost cells integer :: max_blocks_clinic=0 !max number of blocks per processor in each distribution integer :: ncouple_per_day = 48 logical :: lsend_precip_fact ! if T,send precip_fact to cpl for use in fw balance ! (partially-coupled option) + type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure !----------------------------------------------------------------------- @@ -209,6 +211,8 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) call t_startf('MOM_mct_init') + if (debug .and. root_pe().eq.pe_here()) print *, "calling ocn_SetGSMap_mct" + ! Set mct global seg maps: call ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, MOM_MCT_GSMap, MOM_MCT_GSMap3d) @@ -216,16 +220,21 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! Initialize mct ocn domain (needs ocn initialization info) + if (debug .and. root_pe().eq.pe_here()) print *, "calling ocn_domain_mct" call ocn_domain_mct(lsize, MOM_MCT_gsmap, MOM_MCT_dom) - call ocn_domain_mct(lsize*km, MOM_MCT_gsmap3d, MOM_MCT_dom3d) + call ocn_domain_mct(lsize*km, MOM_MCT_gsmap3d, MOM_MCT_dom3d) !TODO: this is not used ! Inialize mct attribute vectors + if (debug .and. root_pe().eq.pe_here()) print *, "calling mct_avect_init a" + ! Initialize the mct attribute vector x2o_o, given Attribute list and length: call mct_aVect_init(x2o_o, rList=seq_flds_x2o_fields, lsize=lsize) ! set the mct attribute vector x2o_o to zero: call mct_aVect_zero(x2o_o) + if (debug .and. root_pe().eq.pe_here()) print *, "calling mct_avect_init b" + ! Initialize the mct attribute vector o2x_o, given Attribute list and length: call mct_aVect_init(o2x_o, rList=seq_flds_o2x_fields, lsize=lsize) ! set the mct attribute vector o2x_o to zero: @@ -238,6 +247,8 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! initialize necessary coupling info + if (debug .and. root_pe().eq.pe_here()) print *, "calling seq_timemgr_eclockgetdata" + call seq_timemgr_EClockGetData(EClock, dtime=ocn_cpl_dt) mom_cpl_dt = seconds_in_day / ncouple_per_day if (mom_cpl_dt /= ocn_cpl_dt) then @@ -252,21 +263,31 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! call seq_infodata_PutData( infodata, precip_fact=precip_fact) ! end if + if (debug .and. root_pe().eq.pe_here()) print *, "calling momo_sum_buffer" call mom_sum_buffer + if (debug .and. root_pe().eq.pe_here()) print *, "calling ocn_export" + call ocn_export(o2x_o%rattr, ldiag_cpl, errorCode) call t_stopf('MOM_mct_init') + if (debug .and. root_pe().eq.pe_here()) print *, "calling get_state_pointers" + + ! Size of global domain + call get_state_pointers(ocn_state, grid=grid) + call get_global_grid_size(grid, ni, nj) + if (debug .and. root_pe().eq.pe_here()) print *, "calling seq_infodata_putdata" call seq_infodata_PutData( infodata, & - ocn_nx = nx_global , ocn_ny = ny_global) + ocn_nx = ni , ocn_ny = nj) call seq_infodata_PutData( infodata, & ocn_prognostic=.true., ocnrof_prognostic=.true.) + if (debug .and. root_pe().eq.pe_here()) print *, "leaving ocean_init_mct" !----------------------------------------------------------------------- !EOC From 8b381416fed31b967f1590ece86c4ad032c28d0f Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Wed, 19 Jul 2017 14:37:30 -0600 Subject: [PATCH 8/9] Wrapped global data --- config_src/mct_driver/coupler_indices.F90 | 9 ++--- config_src/mct_driver/ocn_comp_mct.F90 | 43 +++++++++++++---------- 2 files changed, 29 insertions(+), 23 deletions(-) diff --git a/config_src/mct_driver/coupler_indices.F90 b/config_src/mct_driver/coupler_indices.F90 index bfb42de7f0..e98540084d 100644 --- a/config_src/mct_driver/coupler_indices.F90 +++ b/config_src/mct_driver/coupler_indices.F90 @@ -79,12 +79,13 @@ module coupler_indices end type cpl_indices - ! Module data for storing - type(cpl_indices), public :: ind - contains - subroutine coupler_indices_init( ) + subroutine coupler_indices_init(ind) + + type(cpl_indices), intent(inout) :: ind + + ! Local Variables type(mct_aVect) :: o2x ! temporary type(mct_aVect) :: x2o ! temporary diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index bb67603f06..30e0d1d345 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -37,7 +37,7 @@ module ocn_comp_mct use MOM_grid, only: ocean_grid_type, get_global_grid_size use MOM_error_handler, only: MOM_error, FATAL, is_root_pe use MOM_time_manager, only: time_type, set_date, set_calendar_type, NOLEAP - use coupler_indices, only: coupler_indices_init + use coupler_indices, only: coupler_indices_init, cpl_indices use ocn_import_export, only: SBUFF_SUM, ocn_Export, mom_sum_buffer ! @@ -61,11 +61,16 @@ module ocn_comp_mct private :: ocn_domain_mct ! !PRIVATE MODULE VARIABLES - type(ocean_state_type), pointer :: ocn_state => NULL() ! Private state of ocean - type(ocean_public_type), pointer :: ocn_surface => NULL() ! Public surface state of ocean + type MCT_MOM_Data + type(ocean_state_type), pointer :: ocn_state => NULL() !< Private state of ocean + type(ocean_public_type), pointer :: ocn_surface => NULL() !< Public surface state of ocean - type(seq_infodata_type), pointer :: & - infodata + type(seq_infodata_type), pointer :: infodata + + type(cpl_indices), public :: ind !< Variable IDs + + end type + type(MCT_MOM_Data) :: glb !======================================================================= @@ -148,17 +153,17 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! set (actually, get from mct) the cdata pointers: call seq_cdata_setptrs(cdata_o, id=MOM_MCT_ID, mpicom=mpicom_ocn, & - gsMap=MOM_MCT_gsMap, dom=MOM_MCT_dom, infodata=infodata) + gsMap=MOM_MCT_gsMap, dom=MOM_MCT_dom, infodata=glb%infodata) !--------------------------------------------------------------------- ! Initialize the model run !--------------------------------------------------------------------- - call coupler_indices_init() + call coupler_indices_init(glb%ind) - call seq_infodata_GetData( infodata, case_name=runid ) + call seq_infodata_GetData( glb%infodata, case_name=runid ) - call seq_infodata_GetData( infodata, start_type=starttype) + call seq_infodata_GetData( glb%infodata, start_type=starttype) if ( trim(starttype) == trim(seq_infodata_start_type_start)) then runtype = "initial" @@ -195,13 +200,13 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) npes = num_pes() pe0 = root_pe() - allocate(ocn_surface) - ocn_surface%is_ocean_PE = .true. - allocate(ocn_surface%pelist(npes)) - ocn_surface%pelist(:) = (/(i,i=pe0,pe0+npes)/) + allocate(glb%ocn_surface) + glb%ocn_surface%is_ocean_PE = .true. + allocate(glb%ocn_surface%pelist(npes)) + glb%ocn_surface%pelist(:) = (/(i,i=pe0,pe0+npes)/) ! initialize the MOM6 model - call ocean_model_init(ocn_surface, ocn_state, time_init, time_in) + call ocean_model_init(glb%ocn_surface, glb%ocn_state, time_init, time_in) call t_stopf('MOM_init') @@ -276,14 +281,14 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) if (debug .and. root_pe().eq.pe_here()) print *, "calling get_state_pointers" ! Size of global domain - call get_state_pointers(ocn_state, grid=grid) + call get_state_pointers(glb%ocn_state, grid=grid) call get_global_grid_size(grid, ni, nj) if (debug .and. root_pe().eq.pe_here()) print *, "calling seq_infodata_putdata" - call seq_infodata_PutData( infodata, & + call seq_infodata_PutData( glb%infodata, & ocn_nx = ni , ocn_ny = nj) - call seq_infodata_PutData( infodata, & + call seq_infodata_PutData( glb%infodata, & ocn_prognostic=.true., ocnrof_prognostic=.true.) @@ -373,7 +378,7 @@ subroutine ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, gsMap_ocn, gsMap3d_ocn) type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure integer, allocatable :: gindex(:) ! Indirect indices - call get_state_pointers(ocn_state, grid=grid) + call get_state_pointers(glb%ocn_state, grid=grid) if (.not. associated(grid)) call MOM_error(FATAL, 'ocn_comp_mct.F90, ocn_SetGSMap_mct():' // & 'grid returned from get_state_pointers() was not associated!') @@ -433,7 +438,7 @@ subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn) real(kind=SHR_REAL_R8) :: m2_to_rad2 type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure - call get_state_pointers(ocn_state, grid=grid) + call get_state_pointers(glb%ocn_state, grid=grid) ! set coords to lat and lon, and areas to rad^2 call mct_gGrid_init(GGrid=dom_ocn, CoordChars=trim(seq_flds_dom_coord), & From 3a923a224ceb0a3bbe03034c45628b71b03e8dd5 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 20 Jul 2017 10:49:33 -0600 Subject: [PATCH 9/9] Got initialization working --- config_src/coupled_driver/ocean_model_MOM.F90 | 4 +- config_src/mct_driver/coupler_indices.F90 | 163 ++++++++++++++++++ config_src/mct_driver/ocn_comp_mct.F90 | 53 +++--- config_src/mct_driver/ocn_import_export.F90 | 27 --- 4 files changed, 198 insertions(+), 49 deletions(-) diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index dfcdded6f8..8a5f934d12 100644 --- a/config_src/coupled_driver/ocean_model_MOM.F90 +++ b/config_src/coupled_driver/ocean_model_MOM.F90 @@ -1086,11 +1086,13 @@ subroutine ocean_public_type_chksum(id, timestep, ocn) end subroutine ocean_public_type_chksum !> Returns pointers to objects within ocean_state_type -subroutine get_state_pointers(OS, grid) +subroutine get_state_pointers(OS, grid, surf) type(ocean_state_type), pointer :: OS !< Ocean state type type(ocean_grid_type), optional, pointer :: grid !< Ocean grid + type(surface), optional, pointer :: surf !< Ocean surface state if (present(grid)) grid => OS%grid + if (present(surf)) surf=> OS%state end subroutine get_state_pointers diff --git a/config_src/mct_driver/coupler_indices.F90 b/config_src/mct_driver/coupler_indices.F90 index e98540084d..a2b11a015c 100644 --- a/config_src/mct_driver/coupler_indices.F90 +++ b/config_src/mct_driver/coupler_indices.F90 @@ -1,14 +1,22 @@ module coupler_indices + ! From MCT: use seq_flds_mod, only : seq_flds_x2o_fields, seq_flds_o2x_fields use seq_flds_mod, only : seq_flds_i2o_per_cat, ice_ncat use mct_mod + ! From MOM: + use MOM_grid, only : ocean_grid_type + use MOM_domains, only : pass_var + use MOM_variables, only : surface + implicit none private + public alloc_sbuffer public coupler_indices_init + public time_avg_state type, public :: cpl_indices @@ -77,10 +85,15 @@ module coupler_indices integer, dimension(:), allocatable :: x2o_fracr_col ! fraction of ocean cell used in radiation computations, per column integer, dimension(:), allocatable :: x2o_qsw_fracr_col ! qsw * fracr, per column + real, dimension(:,:,:),allocatable :: time_avg_sbuffer !< time averages of send buffer + real :: accum_time !< time for accumulation + end type cpl_indices contains + + subroutine coupler_indices_init(ind) type(cpl_indices), intent(inout) :: ind @@ -192,4 +205,154 @@ subroutine coupler_indices_init(ind) end subroutine coupler_indices_init + + subroutine alloc_sbuffer(ind, grid, nsend) + + ! Parameters + type(cpl_indices), intent(inout) :: ind + type(ocean_grid_type), intent(in) :: grid + integer, intent(in) :: nsend + + allocate(ind%time_avg_sbuffer(grid%isd:grid%ied,grid%jsd:grid%jed,nsend)) + + end subroutine alloc_sbuffer + + + subroutine time_avg_state(ind, grid, surf_state, dt, reset, last) + + type(cpl_indices), intent(inout) :: ind !< module control structure + type(ocean_grid_type), intent(inout) :: grid !< ocean grid (inout in order to do halo update) + type(surface), intent(in) :: surf_state !< ocean surface state + real, intent(in) :: dt !< time interval to accumulate (seconds) + logical, optional, intent(in) :: reset !< if present and true, reset accumulations + logical, optional, intent(in) :: last !< if present and true, divide by accumulated time + + ! local variables + integer :: i,j,nvar + real :: rtime, slp_L, slp_R, slp_C, u_min, u_max, slope + real, dimension(grid%isd:grid%ied, grid%jsd:grid%jed) :: ssh + + if (present(reset)) then + if (reset) then + ind%time_avg_sbuffer(:,:,:) = 0. + ind%accum_time = 0. + end if + end if + + ! sst + nvar = ind%o2x_So_t + do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + ind%time_avg_sbuffer(i,j,nvar) = ind%time_avg_sbuffer(i,j,nvar)+dt * surf_state%sst(i,j) + end do; end do + + ! sss + nvar = ind%o2x_So_s + do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + ind%time_avg_sbuffer(i,j,nvar) = ind%time_avg_sbuffer(i,j,nvar)+dt * surf_state%sss(i,j) + end do; end do + + + ! u + nvar = ind%o2x_So_u + do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + ind%time_avg_sbuffer(i,j,nvar) = ind%time_avg_sbuffer(i,j,nvar)+dt * & + 0.5*(surf_state%u(I,j)+surf_state%u(I-1,j)) + end do; end do + + ! v + nvar = ind%o2x_So_v + do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + ind%time_avg_sbuffer(i,j,nvar) = ind%time_avg_sbuffer(i,j,nvar)+dt * & + 0.5*(surf_state%v(i,J)+surf_state%v(i,J-1)) + end do; end do + + ! ssh + do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + ssh(i,j) = surf_state%sea_lev(i,j) + end do; end do + call pass_var(ssh, grid%domain) + + ! d/dx ssh + nvar = ind%o2x_So_dhdx + do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + ! This is a simple second-order difference + ! ind%time_avg_sbuffer(i,j,nvar) = ind%time_avg_sbuffer(i,j,nvar) + dt * & + ! 0.5 * (ssh(i+1,j) + ssh(i-1,j)) * grid%IdxT(i,j) * grid%mask2dT(i,j) + ! This is a PLM slope which might be less prone to the A-grid null mode + slp_L = ssh(i,j) - ssh(i-1,j) + slp_R = ssh(i+1,j) - ssh(i,j) + slp_C = 0.5 * (slp_L + slp_R) + if ( (slp_L * slp_R) > 0.0 ) then + ! This limits the slope so that the edge values are bounded by the + ! two cell averages spanning the edge. + u_min = min( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) + u_max = max( ssh(i-1,j), ssh(i,j), ssh(i+1,j) ) + slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) + else + ! Extrema in the mean values require a PCM reconstruction avoid generating + ! larger extreme values. + slope = 0.0 + end if + ind%time_avg_sbuffer(i,j,nvar) = ind%time_avg_sbuffer(i,j,nvar) + dt * slope * grid%IdxT(i,j) * grid%mask2dT(i,j) + end do; end do + + ! d/dy ssh + nvar = ind%o2x_So_dhdy + do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + ! This is a simple second-order difference + ! ind%time_avg_sbuffer(i,j,nvar) = ind%time_avg_sbuffer(i,j,nvar) + dt * & + ! 0.5 * (ssh(i,j+1) + ssh(i,j-1)) * grid%IdyT(i,j) * grid%mask2dT(i,j) + ! This is a PLM slope which might be less prone to the A-grid null mode + slp_L = ssh(i,j) - ssh(i,j-1) + slp_R = ssh(i,j+1) - ssh(i,j) + slp_C = 0.5 * (slp_L + slp_R) + if ( (slp_L * slp_R) > 0.0 ) then + ! This limits the slope so that the edge values are bounded by the + ! two cell averages spanning the edge. + u_min = min( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) + u_max = max( ssh(i,j-1), ssh(i,j), ssh(i,j+1) ) + slope = sign( min( abs(slp_C), 2.*min( ssh(i,j) - u_min, u_max - ssh(i,j) ) ), slp_C ) + else + ! Extrema in the mean values require a PCM reconstruction avoid generating + ! larger extreme values. + slope = 0.0 + end if + ind%time_avg_sbuffer(i,j,nvar) = ind%time_avg_sbuffer(i,j,nvar) + dt * slope * grid%IdyT(i,j) * grid%mask2dT(i,j) + end do; end do + + ! Divide by total accumulated time + ind%accum_time = ind%accum_time + dt + if (present(last)) then + + !! \todo Do dhdx,dhdy need to be rotated before sending to the coupler? + !! \todo Do u,v need to be rotated before sending to the coupler? + + rtime = 1./ind%accum_time + if (last) ind%time_avg_sbuffer(:,:,:) = ind%time_avg_sbuffer(:,:,:) * rtime + end if + + end subroutine time_avg_state + + + subroutine ocn_export(ind, grid, o2x) + + type(cpl_indices), intent(in) :: ind + type(ocean_grid_type), intent(in) :: grid + real(kind=8), intent(inout) :: o2x(:,:) + + integer :: i, j, n + + n = 0 + do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + n = n+1 + o2x(ind%o2x_So_t, n) = ind%time_avg_sbuffer(i,j,ind%o2x_So_t) + o2x(ind%o2x_So_s, n) = ind%time_avg_sbuffer(i,j,ind%o2x_So_s) + o2x(ind%o2x_So_u, n) = ind%time_avg_sbuffer(i,j,ind%o2x_So_u) + o2x(ind%o2x_So_v, n) = ind%time_avg_sbuffer(i,j,ind%o2x_So_v) + o2x(ind%o2x_So_dhdx, n) = ind%time_avg_sbuffer(i,j,ind%o2x_So_dhdx) + o2x(ind%o2x_So_dhdy, n) = ind%time_avg_sbuffer(i,j,ind%o2x_So_dhdy) + end do; end do + + end subroutine ocn_export + end module coupler_indices diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index 30e0d1d345..85a2e5b650 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -31,14 +31,15 @@ module ocn_comp_mct ! From MOM6 - use ocean_model_mod, only: ocean_state_type, ocean_public_type + use ocean_model_mod, only: ocean_state_type, ocean_public_type, ocean_model_init_sfc use ocean_model_mod, only: ocean_model_init, get_state_pointers use MOM_domains, only: MOM_infra_init, num_pes, root_pe, pe_here use MOM_grid, only: ocean_grid_type, get_global_grid_size + use MOM_variables, only: surface use MOM_error_handler, only: MOM_error, FATAL, is_root_pe use MOM_time_manager, only: time_type, set_date, set_calendar_type, NOLEAP - use coupler_indices, only: coupler_indices_init, cpl_indices - use ocn_import_export, only: SBUFF_SUM, ocn_Export, mom_sum_buffer + use coupler_indices, only: coupler_indices_init, cpl_indices, alloc_sbuffer, time_avg_state + use ocn_import_export, only: ocn_Export ! ! !PUBLIC MEMBER FUNCTIONS: @@ -63,7 +64,9 @@ module ocn_comp_mct ! !PRIVATE MODULE VARIABLES type MCT_MOM_Data type(ocean_state_type), pointer :: ocn_state => NULL() !< Private state of ocean - type(ocean_public_type), pointer :: ocn_surface => NULL() !< Public surface state of ocean + type(ocean_public_type), pointer :: ocn_public => NULL() !< Public state of ocean + type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure + type(surface), pointer :: ocn_surface => NULL() !< A pointer to the ocean surface state type(seq_infodata_type), pointer :: infodata @@ -147,8 +150,6 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) logical :: lsend_precip_fact ! if T,send precip_fact to cpl for use in fw balance ! (partially-coupled option) - type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure - !----------------------------------------------------------------------- ! set (actually, get from mct) the cdata pointers: @@ -200,13 +201,19 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) npes = num_pes() pe0 = root_pe() - allocate(glb%ocn_surface) - glb%ocn_surface%is_ocean_PE = .true. - allocate(glb%ocn_surface%pelist(npes)) - glb%ocn_surface%pelist(:) = (/(i,i=pe0,pe0+npes)/) + allocate(glb%ocn_public) + glb%ocn_public%is_ocean_PE = .true. + allocate(glb%ocn_public%pelist(npes)) + glb%ocn_public%pelist(:) = (/(i,i=pe0,pe0+npes)/) + + ! Initialize the MOM6 model + call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in) + + ! Initialize ocn_state%state out of sight + call ocean_model_init_sfc(glb%ocn_state, glb%ocn_public) - ! initialize the MOM6 model - call ocean_model_init(glb%ocn_surface, glb%ocn_state, time_init, time_in) + ! store pointers to components inside MOM + call get_state_pointers(glb%ocn_state, grid=glb%grid, surf=glb%ocn_surface) call t_stopf('MOM_init') @@ -248,7 +255,11 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! allocate send buffer nsend = mct_avect_nRattr(o2x_o) nrecv = mct_avect_nRattr(x2o_o) - allocate (SBUFF_SUM(nx_block,ny_block,max_blocks_clinic,nsend)) + + if (debug .and. root_pe().eq.pe_here()) print *, "calling alloc_sbuffer()", nsend + + call alloc_sbuffer(glb%ind,glb%grid,nsend) + ! initialize necessary coupling info @@ -269,8 +280,9 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) ! end if if (debug .and. root_pe().eq.pe_here()) print *, "calling momo_sum_buffer" - - call mom_sum_buffer + + ! Reset time average of send buffer + call time_avg_state(glb%ind, glb%grid, glb%ocn_surface, 1., reset=.true., last=.true.) if (debug .and. root_pe().eq.pe_here()) print *, "calling ocn_export" @@ -281,8 +293,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) if (debug .and. root_pe().eq.pe_here()) print *, "calling get_state_pointers" ! Size of global domain - call get_state_pointers(glb%ocn_state, grid=grid) - call get_global_grid_size(grid, ni, nj) + call get_global_grid_size(glb%grid, ni, nj) if (debug .and. root_pe().eq.pe_here()) print *, "calling seq_infodata_putdata" @@ -378,7 +389,7 @@ subroutine ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, gsMap_ocn, gsMap3d_ocn) type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure integer, allocatable :: gindex(:) ! Indirect indices - call get_state_pointers(glb%ocn_state, grid=grid) + grid => glb%grid ! for convenience if (.not. associated(grid)) call MOM_error(FATAL, 'ocn_comp_mct.F90, ocn_SetGSMap_mct():' // & 'grid returned from get_state_pointers() was not associated!') @@ -394,9 +405,9 @@ subroutine ocn_SetGSMap_mct(mpicom_ocn, MOM_MCT_ID, gsMap_ocn, gsMap3d_ocn) ! Set indirect indices in gindex k = 0 do j = grid%jsc, grid%jec - jg = j - grid%jdg_offset ! TODO: check this calculation + jg = j + grid%jdg_offset ! TODO: check this calculation do i = grid%isc, grid%iec - ig = i - grid%idg_offset ! TODO: check this calculation + ig = i + grid%idg_offset ! TODO: check this calculation k = k + 1 ! Increment position within gindex gindex(k) = ni * ( jg - 1 ) + ig enddo @@ -438,7 +449,7 @@ subroutine ocn_domain_mct( lsize, gsMap_ocn, dom_ocn) real(kind=SHR_REAL_R8) :: m2_to_rad2 type(ocean_grid_type), pointer :: grid => NULL() ! A pointer to a grid structure - call get_state_pointers(glb%ocn_state, grid=grid) + grid => glb%grid ! for convenience ! set coords to lat and lon, and areas to rad^2 call mct_gGrid_init(GGrid=dom_ocn, CoordChars=trim(seq_flds_dom_coord), & diff --git a/config_src/mct_driver/ocn_import_export.F90 b/config_src/mct_driver/ocn_import_export.F90 index 5cb025a706..a314e60960 100644 --- a/config_src/mct_driver/ocn_import_export.F90 +++ b/config_src/mct_driver/ocn_import_export.F90 @@ -136,33 +136,6 @@ end subroutine ocn_export !*********************************************************************** - subroutine MOM_sum_buffer - -! !DESCRIPTION: -! This routine accumulates sums for averaging fields to -! be sent to the coupler -! -! !REVISION HISTORY: -! same as module -! -!EOP -!BOC - -!----------------------------------------------------------------------- -! -! local variables -! -!----------------------------------------------------------------------- - -!----------------------------------------------------------------------- -!EOC - - end subroutine MOM_sum_buffer - -!*********************************************************************** - - - end module ocn_import_export