From 1f39c84cebed6dbaf0acd6a7df1b4f358349bed0 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 25 Jul 2017 08:47:50 -0600 Subject: [PATCH 1/2] Replaces use of "surface" with "ocean_public_type" - To later facilitate MOM doing the time-averaging within the DT_FORCING time-step, we need to use the ocean_public_type when export surface state to the MCT coupler. - Removes a buffer array inspired by POP interface. - Removes the time-averaging routine we wrote that worked on the POP-like buffer. --- config_src/mct_driver/coupler_indices.F90 | 132 +++++----------------- config_src/mct_driver/ocn_comp_mct.F90 | 20 +--- 2 files changed, 31 insertions(+), 121 deletions(-) diff --git a/config_src/mct_driver/coupler_indices.F90 b/config_src/mct_driver/coupler_indices.F90 index 24e52fa059..dadfdf38e3 100644 --- a/config_src/mct_driver/coupler_indices.F90 +++ b/config_src/mct_driver/coupler_indices.F90 @@ -11,15 +11,14 @@ module coupler_indices use MOM_grid, only : ocean_grid_type ! MOM functions use MOM_domains, only : pass_var - use MOM_variables, only : surface + use ocean_model_mod, only : ocean_public_type implicit none private - public alloc_sbuffer public coupler_indices_init - public time_avg_state + public ocn_export type, public :: cpl_indices @@ -88,9 +87,6 @@ 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 @@ -208,79 +204,38 @@ subroutine coupler_indices_init(ind) end subroutine coupler_indices_init + !> Maps outgoing ocean data to MCT buffer + subroutine ocn_export(ind, ocn_public, grid, o2x) + type(cpl_indices), intent(inout) :: ind !< Index + type(ocean_public_type), intent(in) :: ocn_public !< Ocean surface state + type(ocean_grid_type), intent(in) :: grid !< Ocean model grid + real(kind=8), intent(inout) :: o2x(:,:) !< MCT outgoing bugger + ! Local variables + real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: ssh !< Local copy of sea_lev with updated halo + integer :: i, j, n + real :: slp_L, slp_R, slp_C, slope, u_min, u_max - 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 + n = 0 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)) + n = n+1 + o2x(ind%o2x_So_t, n) = ocn_public%t_surf(i,j) * grid%mask2dT(i,j) + o2x(ind%o2x_So_s, n) = ocn_public%s_surf(i,j) * grid%mask2dT(i,j) + o2x(ind%o2x_So_u, n) = ocn_public%u_surf(i,j) * grid%mask2dT(i,j) + o2x(ind%o2x_So_v, n) = ocn_public%v_surf(i,j) * grid%mask2dT(i,j) end do; end do - ! ssh + ! ssh (make a copy in order to do a halo update) do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec - ssh(i,j) = surf_state%sea_lev(i,j) + ssh(i,j) = ocn_public%sea_lev(i,j) end do; end do call pass_var(ssh, grid%domain) ! d/dx ssh - nvar = ind%o2x_So_dhdx + n = 0 do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec + n = n+1 ! 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) + ! o2x(ind%o2x_So_dhdx, n) = 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) @@ -296,15 +251,13 @@ subroutine time_avg_state(ind, grid, surf_state, dt, reset, last) ! 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) + o2x(ind%o2x_So_dhdx, n) = 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) + ! o2x(ind%o2x_So_dhdy, n) = 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) @@ -320,40 +273,7 @@ subroutine time_avg_state(ind, grid, surf_state, dt, reset, last) ! 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) + o2x(ind%o2x_So_dhdy, n) = slope * grid%IdyT(i,j) * grid%mask2dT(i,j) end do; end do end subroutine ocn_export diff --git a/config_src/mct_driver/ocn_comp_mct.F90 b/config_src/mct_driver/ocn_comp_mct.F90 index ae944a30fc..b98081d914 100644 --- a/config_src/mct_driver/ocn_comp_mct.F90 +++ b/config_src/mct_driver/ocn_comp_mct.F90 @@ -43,8 +43,8 @@ module ocn_comp_mct 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 + use coupler_indices, only: coupler_indices_init, cpl_indices + use coupler_indices, only: ocn_export ! ! !PUBLIC MEMBER FUNCTIONS: @@ -71,7 +71,6 @@ module ocn_comp_mct 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 @@ -229,6 +228,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) glb%ocn_public%is_ocean_PE = .true. allocate(glb%ocn_public%pelist(npes)) glb%ocn_public%pelist(:) = (/(i,i=pe0,pe0+npes)/) + ! \todo Set other bits of glb$ocn_public ! Initialize the MOM6 model call ocean_model_init(glb%ocn_public, glb%ocn_state, time_init, time_in) @@ -237,7 +237,7 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) call ocean_model_init_sfc(glb%ocn_state, glb%ocn_public) ! store pointers to components inside MOM - call get_state_pointers(glb%ocn_state, grid=glb%grid, surf=glb%ocn_surface) + call get_state_pointers(glb%ocn_state, grid=glb%grid) call t_stopf('MOM_init') @@ -280,11 +280,6 @@ subroutine ocn_init_mct( EClock, cdata_o, x2o_o, o2x_o, NLFilename ) 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" @@ -305,14 +300,9 @@ 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" - - ! 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 ocn_export(glb%ind, glb%ocn_public, glb%grid, o2x_o%rattr) call t_stopf('MOM_mct_init') From c83e23833661bc14b34a8b9f38cfd871ef2c4ac9 Mon Sep 17 00:00:00 2001 From: Alistair Adcroft Date: Tue, 25 Jul 2017 09:32:03 -0600 Subject: [PATCH 2/2] Correct indexing when using ocean_public_type - The ocean_public_type uses global indexing without halos. --- config_src/mct_driver/coupler_indices.F90 | 33 ++++++++++++++--------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/config_src/mct_driver/coupler_indices.F90 b/config_src/mct_driver/coupler_indices.F90 index dadfdf38e3..c0ea1eab44 100644 --- a/config_src/mct_driver/coupler_indices.F90 +++ b/config_src/mct_driver/coupler_indices.F90 @@ -212,22 +212,29 @@ subroutine ocn_export(ind, ocn_public, grid, o2x) real(kind=8), intent(inout) :: o2x(:,:) !< MCT outgoing bugger ! Local variables real, dimension(grid%isd:grid%ied,grid%jsd:grid%jed) :: ssh !< Local copy of sea_lev with updated halo - integer :: i, j, n + integer :: i, j, n, ig, jg real :: slp_L, slp_R, slp_C, slope, u_min, u_max + ! Copy from ocn_public to o2x. ocn_public uses global indexing with no halos. + ! The mask comes from "grid" that uses the usual MOM domain that has halos + ! and does not use global indexing. n = 0 - do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec - n = n+1 - o2x(ind%o2x_So_t, n) = ocn_public%t_surf(i,j) * grid%mask2dT(i,j) - o2x(ind%o2x_So_s, n) = ocn_public%s_surf(i,j) * grid%mask2dT(i,j) - o2x(ind%o2x_So_u, n) = ocn_public%u_surf(i,j) * grid%mask2dT(i,j) - o2x(ind%o2x_So_v, n) = ocn_public%v_surf(i,j) * grid%mask2dT(i,j) - end do; end do - - ! ssh (make a copy in order to do a halo update) - do j=grid%jsc, grid%jec ; do i=grid%isc,grid%iec - ssh(i,j) = ocn_public%sea_lev(i,j) - end do; end do + do j=grid%jsc, grid%jec + jg = j + grid%jdg_offset + do i=grid%isc,grid%iec + n = n+1 + ig = i + grid%idg_offset + o2x(ind%o2x_So_t, n) = ocn_public%t_surf(ig,jg) * grid%mask2dT(i,j) + o2x(ind%o2x_So_s, n) = ocn_public%s_surf(ig,jg) * grid%mask2dT(i,j) + o2x(ind%o2x_So_u, n) = ocn_public%u_surf(ig,jg) * grid%mask2dT(i,j) + o2x(ind%o2x_So_v, n) = ocn_public%v_surf(ig,jg) * grid%mask2dT(i,j) + ! Make a copy of ssh in order to do a halo update. We use the usual MOM domain + ! in order to update halos. i.e. does not use global indexing. + ssh(i,j) = ocn_public%sea_lev(ig,jg) + end do + end do + + ! Update halo of ssh so we can calculate gradients call pass_var(ssh, grid%domain) ! d/dx ssh