Got initialization working

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
 

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

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), &

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