From 3a923a224ceb0a3bbe03034c45628b71b03e8dd5 Mon Sep 17 00:00:00 2001 From: alperaltuntas Date: Thu, 20 Jul 2017 10:49:33 -0600 Subject: [PATCH] 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