Skip to content

update cap to receive u_surf,v_surf on C-grid #171

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
DeniseWorthen wants to merge 1 commit into
base: dev/emc
Choose a base branch
from
Draft

update cap to receive u_surf,v_surf on C-grid #171

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 22 additions & 6 deletions config_src/drivers/nuopc_cap/mom_cap.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
121 changes: 90 additions & 31 deletions config_src/drivers/nuopc_cap/mom_cap_methods.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -572,19 +572,21 @@ 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
integer :: i, j, ig, jg ! indices
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
Expand All @@ -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
Expand Down Expand Up @@ -657,21 +663,48 @@ 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))

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)
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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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

Expand Down
6 changes: 4 additions & 2 deletions config_src/drivers/nuopc_cap/mom_ocean_model_nuopc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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.)
Expand Down Expand Up @@ -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 "//&
Expand Down
Loading