diff --git a/config_src/coupled_driver/ocean_model_MOM.F90 b/config_src/coupled_driver/ocean_model_MOM.F90 index 0f0f1f2622..8a5f934d12 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,15 @@ 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, 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 + end module ocean_model_mod diff --git a/config_src/mct_driver/coupler_indices.F90 b/config_src/mct_driver/coupler_indices.F90 index bfb42de7f0..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,14 +85,20 @@ 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 - end type cpl_indices + real, dimension(:,:,:),allocatable :: time_avg_sbuffer !< time averages of send buffer + real :: accum_time !< time for accumulation - ! Module data for storing - type(cpl_indices), public :: ind + end type cpl_indices 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 @@ -191,4 +205,154 @@ subroutine coupler_indices_init( ) 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 b3c760c0fd..85a2e5b650 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -6,26 +6,40 @@ 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, & + SEQ_FLDS_DOM_COORD, & + SEQ_FLDS_DOM_other + 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, & + 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 + use shr_kind_mod, only: SHR_KIND_R8 + + ! 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, 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, alloc_sbuffer, time_avg_state + use ocn_import_export, only: ocn_Export ! ! !PUBLIC MEMBER FUNCTIONS: @@ -34,6 +48,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: @@ -43,14 +58,22 @@ 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 MCT_MOM_Data + type(ocean_state_type), pointer :: ocn_state => NULL() !< Private 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 - type(seq_infodata_type), pointer :: & - infodata + type(cpl_indices), public :: ind !< Variable IDs + + end type + type(MCT_MOM_Data) :: glb !======================================================================= @@ -65,7 +88,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: @@ -90,19 +113,87 @@ 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 :: i, errorCode + integer :: lsize, nsend, nrecv + logical :: ldiag_cpl = .false. + integer :: ni, nj + + ! mct variables (these are local for now) + integer :: MOM_MCT_ID + 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=1 ! number of vertical levels + integer :: nx_block=0, ny_block=0 ! 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) + +!----------------------------------------------------------------------- + ! 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=glb%infodata) + !--------------------------------------------------------------------- + ! Initialize the model run + !--------------------------------------------------------------------- + + call coupler_indices_init(glb%ind) + + call seq_infodata_GetData( glb%infodata, case_name=runid ) + + call seq_infodata_GetData( glb%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 - mpicom = cdata_o%mpicom - call MOM_infra_init(mpicom) + !--------------------------------------------------------------------- + + call t_startf('MOM_init') - call ESMF_ClockGet(EClock, currTime=current_time, rc=rc) + 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) @@ -110,36 +201,109 @@ 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_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 model run - call ocean_model_init(ocn_surface, ocn_state, time_init, time_in) + ! Initialize the MOM6 model + call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in) - ! set infodata, a cdata pointer - call seq_cdata_setptrs(cdata_o, infodata=infodata) + ! Initialize ocn_state%state out of sight + call ocean_model_init_sfc(glb%ocn_state, glb%ocn_public) - ! initialize coupler indices - call coupler_indices_init() + ! store pointers to components inside MOM + call get_state_pointers(glb%ocn_state, grid=glb%grid, surf=glb%ocn_surface) - ! get runid and starttype: - call seq_infodata_GetData( infodata, case_name=runid ) - call seq_infodata_GetData( infodata, start_type=starttype) + call t_stopf('MOM_init') + !--------------------------------------------------------------------- + ! Initialize MCT attribute vectors and indices + !--------------------------------------------------------------------- - 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 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) + lsize = mct_gsMap_lsize(MOM_MCT_gsmap, mpicom_ocn) + + ! 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) !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: + call mct_aVect_zero(o2x_o) + + ! allocate send buffer + nsend = mct_avect_nRattr(o2x_o) + nrecv = mct_avect_nRattr(x2o_o) + + 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 + + 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 + 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 + + if (debug .and. root_pe().eq.pe_here()) print *, "calling momo_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" + + 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_global_grid_size(glb%grid, ni, nj) + + if (debug .and. root_pe().eq.pe_here()) print *, "calling seq_infodata_putdata" + + call seq_infodata_PutData( glb%infodata, & + ocn_nx = ni , ocn_ny = nj) + call seq_infodata_PutData( glb%infodata, & + ocn_prognostic=.true., ocnrof_prognostic=.true.) + + + if (debug .and. root_pe().eq.pe_here()) print *, "leaving ocean_init_mct" !----------------------------------------------------------------------- !EOC @@ -209,6 +373,163 @@ subroutine ocn_final_mct( EClock, cdata_o, x2o_o, o2x_o) end subroutine ocn_final_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 + + 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!') + + ! 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 + + +!*********************************************************************** +!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 + +! 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 + + 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), & + 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 +!----------------------------------------------------------------------- +! +! 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..a314e60960 --- /dev/null +++ b/config_src/mct_driver/ocn_import_export.F90 @@ -0,0 +1,141 @@ +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 + 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