From 1082f8fb45b40469f08e3c54ecb87c3eb25dfc7f Mon Sep 17 00:00:00 2001 From: Denise Worthen Date: Wed, 28 Jan 2026 11:26:04 -0500 Subject: [PATCH] update cap to receive u_surf,v_surf on C-grid --- config_src/drivers/nuopc_cap/mom_cap.F90 | 28 +++- .../drivers/nuopc_cap/mom_cap_methods.F90 | 121 +++++++++++++----- .../nuopc_cap/mom_ocean_model_nuopc.F90 | 6 +- 3 files changed, 116 insertions(+), 39 deletions(-) diff --git a/config_src/drivers/nuopc_cap/mom_cap.F90 b/config_src/drivers/nuopc_cap/mom_cap.F90 index ad8cbf3dce..314e057c8e 100644 --- a/config_src/drivers/nuopc_cap/mom_cap.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap.F90 @@ -109,13 +109,13 @@ module MOM_cap_mod type(ocean_public_type), pointer :: ocean_public_type_ptr type(ocean_state_type), pointer :: ocean_state_type_ptr type(ice_ocean_boundary_type), pointer :: ice_ocean_boundary_type_ptr -end type +end type ocean_internalstate_type !> Wrapper-derived type required to associate an internal state instance !! with the ESMF/NUOPC component type ocean_internalstate_wrapper type(ocean_internalstate_type), pointer :: ptr -end type +end type ocean_internalstate_wrapper !> Contains field information type fld_list_type @@ -167,6 +167,7 @@ module MOM_cap_mod logical :: pointer_date = .true. ! append date to rpointer real(8) :: timere integer :: mype = -1 +character(len=1) :: ocean_surface_stagger = '' contains @@ -749,10 +750,21 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) ocean_public%is_ocean_pe = .true. if (cesm_coupled .and. len_trim(inst_suffix)>0) then - call ocean_model_init(ocean_public, ocean_state, time0, time_start, & + call ocean_model_init(ocean_public, ocean_state, time0, time_start, ocean_surface_stagger, & input_restart_file=trim(adjustl(restartfiles)), inst_index=inst_index) else - call ocean_model_init(ocean_public, ocean_state, time0, time_start, input_restart_file=trim(adjustl(restartfiles))) + call ocean_model_init(ocean_public, ocean_state, time0, time_start, ocean_surface_stagger, & + input_restart_file=trim(adjustl(restartfiles))) + endif + if (ocean_surface_stagger /= 'A' .and. ocean_surface_stagger /= 'C') then + call MOM_error(FATAL,'OCEAN_SURFACE_STAGGER must be A or C for NUOPC cap ') + end if + if (is_root_pe()) then + if (ocean_surface_stagger == 'A') write(stdout,*) 'ocean_surface_stagger is A: ', & + 'all exports from NUOPC cap are on the A grid' + if (ocean_surface_stagger == 'C') write(stdout,*) 'ocean_surface_stagger is C: ', & + 'surface slopes and velocities will be exported on C grid; A-grid velocites ', & + 'will also be exported from the NUOPC cap' endif ! GMM, this call is not needed in CESM. Check with EMC if it can be deleted. @@ -935,6 +947,10 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (cesm_coupled .and. use_MARBL) then call fld_list_add(fldsFrOcn_num, fldsFrOcn, "Faoo_fco2_ocn", "will provide") endif + if (ocean_surface_stagger == 'C') then + call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_uc" , "will provide") + call fld_list_add(fldsFrOcn_num, fldsFrOcn, "So_vc" , "will provide") + end if do n = 1,fldsToOcn_num call NUOPC_Advertise(importState, standardName=fldsToOcn(n)%stdname, name=fldsToOcn(n)%shortname, rc=rc) @@ -1689,7 +1705,7 @@ subroutine DataInitialize(gcomp, rc) ocean_state => ocean_internalstate%ptr%ocean_state_type_ptr call get_ocean_grid(ocean_state, ocean_grid) - call mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, rc=rc) + call mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, ocean_surface_stagger, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call ESMF_StateGet(exportState, itemCount=fieldCount, rc=rc) @@ -1928,7 +1944,7 @@ subroutine ModelAdvance(gcomp, rc) ! Export Data !--------------- - call mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, rc=rc) + call mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, ocean_surface_stagger, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return if (dbug > 0) then diff --git a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 index 8a7a1b3942..17eba11a57 100644 --- a/config_src/drivers/nuopc_cap/mom_cap_methods.F90 +++ b/config_src/drivers/nuopc_cap/mom_cap_methods.F90 @@ -572,12 +572,13 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary, end subroutine mom_import !> Maps outgoing ocean data to ESMF State -subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, rc) +subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, ocean_surface_stagger, rc) type(ocean_public_type) , intent(in) :: ocean_public !< Ocean surface state type(ocean_grid_type) , intent(in) :: ocean_grid !< Ocean model grid type(ocean_state_type) , pointer :: ocean_state !< Ocean state type(ESMF_State) , intent(inout) :: exportState !< outgoing data type(ESMF_Clock) , intent(in) :: clock !< ESMF clock + character(len=*) , intent(in) :: ocean_surface_stagger !< stagger setting for u_surf,v_surf integer , intent(inout) :: rc !< Return code ! Local variables @@ -585,6 +586,7 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, integer :: isc, iec, jsc, jec ! indices integer :: iloc, jloc ! indices integer :: iglob, jglob ! indices + integer :: i0,j0 integer :: n integer :: icount real :: slp_L, slp_R, slp_C @@ -601,6 +603,10 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, real(ESMF_KIND_R8), allocatable :: ssh(:,:) real(ESMF_KIND_R8), allocatable :: dhdx(:,:), dhdy(:,:) real(ESMF_KIND_R8), allocatable :: dhdx_rot(:,:), dhdy_rot(:,:) + real(ESMF_KIND_R8), allocatable :: uc(:,:) + real(ESMF_KIND_R8), allocatable :: vc(:,:) + real(ESMF_KIND_R8), allocatable :: workx(:,:) + real(ESMF_KIND_R8), allocatable :: worky(:,:) character(len=*) , parameter :: subname = '(mom_export)' rc = ESMF_SUCCESS @@ -657,12 +663,41 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, ! ------- ! zonal and meridional currents ! ------- + allocate(ocz(isc:iec, jsc:jec)) + allocate(ocm(isc:iec, jsc:jec)) + if (ocean_surface_stagger == 'C') then + call State_SetExport(exportState, 'So_uc', isc, iec, jsc, jec, ocean_public%u_surf, ocean_grid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call State_SetExport(exportState, 'So_vc', isc, iec, jsc, jec, ocean_public%v_surf, ocean_grid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! Local copy of C-grid velocities with halos to calculate A-grid velocities + allocate(uc(ocean_grid%isd:ocean_grid%ied, ocean_grid%jsd:ocean_grid%jed), source=0.0_ESMF_KIND_R8) + allocate(vc(ocean_grid%isd:ocean_grid%ied, ocean_grid%jsd:ocean_grid%jed), source=0.0_ESMF_KIND_R8) + + uc(ocean_grid%isc:ocean_grid%iec, ocean_grid%jsc:ocean_grid%jec) = ocean_public%u_surf(isc:iec, jsc:jec) + vc(ocean_grid%isc:ocean_grid%iec, ocean_grid%jsc:ocean_grid%jec) = ocean_public%v_surf(isc:iec, jsc:jec) + + call pass_var(uc, ocean_grid%Domain) + call pass_var(vc, ocean_grid%Domain) + i0 = ocean_grid%isc-isc + j0 = ocean_grid%jsc-jsc + do j=jsc,jec ; do i=isc,iec + ocz(i,j) = ocean_grid%mask2dT(i+i0,j+j0) * 0.5*(uc(i+i0,j+j0) + uc(i+i0-1,j+j0)) + ocm(i,j) = ocean_grid%mask2dT(i+i0,j+j0) * 0.5*(vc(i+i0,j+j0) + vc(i+i0,j+j0-1)) + enddo; enddo + else + do j = jsc, jec + jg = j + ocean_grid%jsc - jsc + do i = isc, iec + ig = i + ocean_grid%isc - isc + ocz(i,j) = ocean_public%u_surf(i,j) + ocm(i,j) = ocean_public%v_surf(i,j) + enddo + enddo + endif ! rotate ocn current from tripolar grid back to lat/lon grid x,y => latlon (CCW) ! "ocean_grid" has halos and uses local indexing. - - allocate(ocz(isc:iec, jsc:jec)) - allocate(ocm(isc:iec, jsc:jec)) allocate(ocz_rot(isc:iec, jsc:jec)) allocate(ocm_rot(isc:iec, jsc:jec)) @@ -670,8 +705,6 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, jg = j + ocean_grid%jsc - jsc do i = isc, iec ig = i + ocean_grid%isc - isc - ocz(i,j) = ocean_public%u_surf(i,j) - ocm(i,j) = ocean_public%v_surf(i,j) ocz_rot(i,j) = ocean_grid%cos_rot(ig,jg)*ocz(i,j) + ocean_grid%sin_rot(ig,jg)*ocm(i,j) ocm_rot(i,j) = ocean_grid%cos_rot(ig,jg)*ocm(i,j) - ocean_grid%sin_rot(ig,jg)*ocz(i,j) enddo @@ -684,6 +717,8 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, if (ChkErr(rc,__LINE__,u_FILE_u)) return deallocate(ocz, ocm, ocz_rot, ocm_rot) + if(allocated(uc))deallocate(uc) + if(allocated(vc))deallocate(vc) ! ------- ! Boundary layer depth @@ -735,11 +770,9 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, !---------------- allocate(ssh(ocean_grid%isd:ocean_grid%ied,ocean_grid%jsd:ocean_grid%jed), & ! local indices with halos - dhdx(isc:iec, jsc:jec), & !global indices without halos - dhdy(isc:iec, jsc:jec), & !global indices without halos - source=0.0_ESMF_KIND_R8) - allocate(dhdx_rot(isc:iec, jsc:jec)) !global indices without halos - allocate(dhdy_rot(isc:iec, jsc:jec)) !global indices without halos + dhdx(ocean_grid%isd:ocean_grid%ied,ocean_grid%jsd:ocean_grid%jed), & !global indices with halos + dhdy(ocean_grid%isd:ocean_grid%ied,ocean_grid%jsd:ocean_grid%jed), & !global indices with halos + source=0.0_ESMF_KIND_R8) ! Make a copy of ssh in order to do a halo update (ssh has local indexing with halos) do j = ocean_grid%jsc, ocean_grid%jec @@ -778,8 +811,8 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, ! larger extreme values. slope = 0.0 endif - dhdx(iglob,jglob) = slope * ocean_grid%US%m_to_L*ocean_grid%IdxT(i,j) * ocean_grid%mask2dT(i,j) - if (ocean_grid%mask2dT(i,j)==0.) dhdx(iglob,jglob) = 0.0 + dhdx(i,j) = slope * ocean_grid%US%m_to_L*ocean_grid%IdxT(i,j) * ocean_grid%mask2dT(i,j) + if (ocean_grid%mask2dT(i,j)==0.) dhdx(i,j) = 0.0 enddo enddo @@ -808,28 +841,54 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, ! larger extreme values. slope = 0.0 endif - dhdy(iglob,jglob) = slope * ocean_grid%US%m_to_L*ocean_grid%IdyT(i,j) * ocean_grid%mask2dT(i,j) - if (ocean_grid%mask2dT(i,j)==0.) dhdy(iglob,jglob) = 0.0 + dhdy(i,j) = slope * ocean_grid%US%m_to_L*ocean_grid%IdyT(i,j) * ocean_grid%mask2dT(i,j) + if (ocean_grid%mask2dT(i,j)==0.) dhdy(i,j) = 0.0 enddo enddo - - ! rotate slopes from tripolar grid back to lat/lon grid, x,y => latlon (CCW) - ! "ocean_grid" uses has halos and uses local indexing. - - do j = jsc, jec - jg = j + ocean_grid%jsc - jsc - do i = isc, iec - ig = i + ocean_grid%isc - isc - dhdx_rot(i,j) = ocean_grid%cos_rot(ig,jg)*dhdx(i,j) + ocean_grid%sin_rot(ig,jg)*dhdy(i,j) - dhdy_rot(i,j) = ocean_grid%cos_rot(ig,jg)*dhdy(i,j) - ocean_grid%sin_rot(ig,jg)*dhdx(i,j) + ! Update halo of slopes to calculate gradients on C-grid + call pass_var(dhdx, ocean_grid%domain) + call pass_var(dhdy, ocean_grid%domain) + + if (ocean_surface_stagger == 'C') then + allocate(workx, mold=dhdx) + allocate(worky, mold=dhdy) + workx = 0.0 + worky = 0.0 + do jglob = jsc, jec + j = jglob + ocean_grid%jsc - jsc + do iglob = isc,iec + i = iglob + ocean_grid%isc - isc + workx(i,j) = 0.5*(dhdx(i,j) + dhdx(i+1,j))*ocean_grid%mask2dCu(i,j) + worky(i,j) = 0.5*(dhdy(i,j) + dhdy(i,j+1))*ocean_grid%mask2dCv(i,j) + enddo enddo - enddo - call State_SetExport(exportState, 'So_dhdx', isc, iec, jsc, jec, dhdx_rot, ocean_grid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + call State_SetExport(exportState, 'So_dhdx', isc, iec, jsc, jec, workx, ocean_grid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call State_SetExport(exportState, 'So_dhdy', isc, iec, jsc, jec, worky, ocean_grid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + deallocate(workx, worky) + else + ! rotate slopes from tripolar grid back to lat/lon grid, x,y => latlon (CCW) + ! "ocean_grid" uses has halos and uses local indexing. + allocate(dhdx_rot(isc:iec, jsc:jec), source=0.0_ESMF_KIND_R8) !global indices without halos + allocate(dhdy_rot(isc:iec, jsc:jec), source=0.0_ESMF_KIND_R8) !global indices without halos + + do j = jsc, jec + jg = j + ocean_grid%jsc - jsc + do i = isc, iec + ig = i + ocean_grid%isc - isc + dhdx_rot(i,j) = ocean_grid%cos_rot(ig,jg)*dhdx(ig,jg) + ocean_grid%sin_rot(ig,jg)*dhdy(ig,jg) + dhdy_rot(i,j) = ocean_grid%cos_rot(ig,jg)*dhdy(ig,jg) - ocean_grid%sin_rot(ig,jg)*dhdx(ig,jg) + enddo + enddo - call State_SetExport(exportState, 'So_dhdy', isc, iec, jsc, jec, dhdy_rot, ocean_grid, rc=rc) - if (ChkErr(rc,__LINE__,u_FILE_u)) return + call State_SetExport(exportState, 'So_dhdx', isc, iec, jsc, jec, dhdx_rot, ocean_grid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call State_SetExport(exportState, 'So_dhdy', isc, iec, jsc, jec, dhdy_rot, ocean_grid, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + deallocate(dhdx_rot, dhdy_rot) + endif ! ------- ! CO2 Flux @@ -841,7 +900,7 @@ subroutine mom_export(ocean_public, ocean_grid, ocean_state, exportState, clock, if (ChkErr(rc,__LINE__,u_FILE_u)) return endif - deallocate(ssh, dhdx, dhdy, dhdx_rot, dhdy_rot) + deallocate(ssh, dhdx, dhdy) end subroutine mom_export diff --git a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 index a83576028a..969cae6acb 100644 --- a/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 +++ b/config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90 @@ -230,7 +230,7 @@ module MOM_ocean_model_nuopc !! This subroutine initializes both the ocean state and the ocean surface type. !! Because of the way that indicies and domains are handled, Ocean_sfc must have !! been used in a previous call to initialize_ocean_type. -subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, input_restart_file, inst_index) +subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, ocean_surface_stagger, gas_fields_ocn, input_restart_file, inst_index) type(ocean_public_type), target, & intent(inout) :: Ocean_sfc !< A structure containing various publicly !! visible ocean surface properties after initialization, @@ -240,6 +240,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i !! contain all information about the ocean's interior state. type(time_type), intent(in) :: Time_init !< The start time for the coupled model's calendar type(time_type), intent(in) :: Time_in !< The time at which to initialize the ocean model. + character(len=1), intent(out) :: ocean_surface_stagger !< location of u,v velocities in ocean_public type(coupler_1d_bc_type), & optional, intent(in) :: gas_fields_ocn !< If present, this type describes the !! ocean and surface-ice fields that will participate @@ -377,6 +378,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i else use_melt_pot=.false. endif + ocean_surface_stagger = uppercase(stagger(1:1)) call get_param(param_file, mdl, "USE_WAVES", OS%Use_Waves, & "If true, enables surface wave modules.", default=.false.) @@ -433,7 +435,7 @@ subroutine ocean_model_init(Ocean_sfc, OS, Time_init, Time_in, gas_fields_ocn, i endif call extract_surface_state(OS%MOM_CSp, OS%sfc_state) -! get number of processors and PE list for stocasthci physics initialization +! get number of processors and PE list for stochastic physics initialization call get_param(param_file, mdl, "DO_SPPT", OS%do_sppt, & "If true, then stochastically perturb the thermodynamic "//& "tendencies of T,S, and h. Amplitude and correlations are "//&