Completed ocn_init_mct

Author: alperaltuntas

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,52 +201,67 @@ 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)
 
   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)
 
+  ! 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
 
 
 !***********************************************************************