Skip to content

Commit

Permalink
Merge branch 'dev/emc' into feature/optmesh
Browse files Browse the repository at this point in the history
  • Loading branch information
@DeniseWorthen
DeniseWorthen committed Jan 29, 2021
2 parents 13a5a2e + d531a32 commit 12b3895
Show file tree
Hide file tree
Showing 11 changed files with 123 additions and 165 deletions.
25 changes: 1 addition & 24 deletions config_src/coupled_driver/MOM_surface_forcing_gfdl.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -189,20 +189,13 @@ module MOM_surface_forcing_gfdl
!! ice-shelves, expressed as a coefficient
!! for divergence damping, as determined
!! outside of the ocean model [m3 s-1]
real, pointer, dimension(:,:) :: ustk0 => NULL() !<
real, pointer, dimension(:,:) :: vstk0 => NULL() !<
real, pointer, dimension(:) :: stk_wavenumbers => NULL() !<
real, pointer, dimension(:,:,:) :: ustkb => NULL() !<
real, pointer, dimension(:,:,:) :: vstkb => NULL() !<

integer :: xtype !< The type of the exchange - REGRID, REDIST or DIRECT
type(coupler_2d_bc_type) :: fluxes !< A structure that may contain an array of named fields
!! used for passive tracer fluxes.
integer :: wind_stagger = -999 !< A flag indicating the spatial discretization of wind stresses.
!! This flag may be set by the flux-exchange code, based on what
!! the sea-ice model is providing. Otherwise, the value from
!! the surface_forcing_CS is used.
integer :: num_stk_bands !< Number of Stokes drift bands passed through the coupler
end type ice_ocean_boundary_type

integer :: id_clock_forcing !< A CPU time clock
Expand Down Expand Up @@ -684,7 +677,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
real :: mass_eff ! effective mass of sea ice for rigidity [R Z ~> kg m-2]
real :: wt1, wt2 ! Relative weights of previous and current values of ustar, ND.

integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0, istk
integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq, i0, j0
integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, isr, ier, jsr, jer
integer :: isc_bnd, iec_bnd, jsc_bnd, jec_bnd

Expand Down Expand Up @@ -725,9 +718,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
(associated(IOB%mass_berg) .and. (.not. associated(forces%mass_berg))) ) &
call allocate_mech_forcing(G, forces, iceberg=.true.)

if ( associated(IOB%ustk0) ) &
call allocate_mech_forcing(G, forces, waves=.true., num_stk_bands=IOB%num_stk_bands)

if (associated(IOB%ice_rigidity)) then
rigidity_at_h(:,:) = 0.0
call safe_alloc_ptr(forces%rigidity_ice_u,IsdB,IedB,jsd,jed)
Expand Down Expand Up @@ -787,19 +777,6 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS, dt_
forces%ustar(i,j) = wt1*forces%ustar(i,j) + wt2*ustar_tmp(i,j)
enddo ; enddo
endif
forces%stk_wavenumbers(:) = IOB%stk_wavenumbers
do j=js,je; do i=is,ie
forces%ustk0(i,j) = IOB%ustk0(i-I0,j-J0) ! How to be careful here that the domains are right?
forces%vstk0(i,j) = IOB%vstk0(i-I0,j-J0)
enddo ; enddo
call pass_vector(forces%ustk0,forces%vstk0, G%domain )
do istk = 1,IOB%num_stk_bands
do j=js,je; do i=is,ie
forces%ustkb(i,j,istk) = IOB%ustkb(i-I0,j-J0,istk)
forces%vstkb(i,j,istk) = IOB%vstkb(i-I0,j-J0,istk)
enddo; enddo
call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
enddo

! Find the net mass source in the input forcing without other adjustments.
if (CS%approx_net_mass_src .and. associated(forces%net_mass_src)) then
Expand Down
2 changes: 1 addition & 1 deletion config_src/coupled_driver/ocean_model_MOM.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -565,7 +565,7 @@ subroutine update_ocean_model(Ice_ocean_boundary, OS, Ocean_sfc, time_start_upda
! For now, the waves are only updated on the thermodynamics steps, because that is where
! the wave intensities are actually used to drive mixing. At some point, the wave updates
! might also need to become a part of the ocean dynamics, according to B. Reichl.
call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves, OS%forces)
call Update_Surface_Waves(OS%grid, OS%GV, OS%US, OS%time, ocean_coupling_time_step, OS%waves)
endif

if ((OS%nstep==0) .and. (OS%nstep_thermo==0)) then ! This is the first call to update_ocean_model.
Expand Down
12 changes: 6 additions & 6 deletions config_src/mct_driver/mom_surface_forcing_mct.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -486,17 +486,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,

! latent heat flux (W/m^2)
fluxes%latent(i,j) = 0.0
! contribution from frozen ppt
! contribution from frozen ppt (notice minus sign since fprec is positive into the ocean)
if (associated(IOB%fprec)) then
fluxes%latent(i,j) = fluxes%latent(i,j) + &
fluxes%latent(i,j) = fluxes%latent(i,j) - &
IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
fluxes%latent_fprec_diag(i,j) = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
fluxes%latent_fprec_diag(i,j) = - G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
endif
! contribution from frozen runoff
! contribution from frozen runoff (notice minus sign since rofi_flux is positive into the ocean)
if (associated(fluxes%frunoff)) then
fluxes%latent(i,j) = fluxes%latent(i,j) + &
fluxes%latent(i,j) = fluxes%latent(i,j) - &
IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * IOB%rofi_flux(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
endif
! contribution from evaporation
if (associated(IOB%q_flux)) then
Expand Down
13 changes: 1 addition & 12 deletions config_src/nuopc_driver/mom_cap.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ module MOM_cap_mod
use MOM_ocean_model_nuopc, only: get_ocean_grid, get_eps_omesh
use MOM_cap_time, only: AlarmInit
use MOM_cap_methods, only: mom_import, mom_export, mom_set_geomtype, state_diagnose
use MOM_cap_methods, only: ChkErr
#ifdef CESMCOUPLED
use shr_file_mod, only: shr_file_setLogUnit, shr_file_getLogUnit
#endif
Expand Down Expand Up @@ -2072,18 +2073,6 @@ subroutine shr_file_getLogUnit(nunit)
end subroutine shr_file_getLogUnit
#endif

logical function chkerr(rc, line, file)
integer, intent(in) :: rc
integer, intent(in) :: line
character(len=*), intent(in) :: file
integer :: lrc
chkerr = .false.
lrc = rc
if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then
chkerr = .true.
endif
end function chkerr

!>
!! @page nuopc_cap NUOPC Cap
!! @author Fei Liu ([email protected])
Expand Down
76 changes: 37 additions & 39 deletions config_src/nuopc_driver/mom_cap_methods.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ module MOM_cap_methods
public :: mom_import
public :: mom_export
public :: state_diagnose
public :: ChkErr

private :: State_getImport
private :: State_setExport
Expand Down Expand Up @@ -251,9 +252,9 @@ subroutine mom_import(ocean_public, ocean_grid, importState, ice_ocean_boundary,


!----
! Partitioned Stokes Drift Components
! Partitioned Stokes Drift Components
!----
if ( associated(ice_ocean_boundary%ustkb) ) then
if ( associated(ice_ocean_boundary%ustkb) ) then
allocate(stkx1(isc:iec,jsc:jec))
allocate(stky1(isc:iec,jsc:jec))
allocate(stkx2(isc:iec,jsc:jec))
Expand Down Expand Up @@ -765,15 +766,18 @@ subroutine State_SetExport(state, fldname, isc, iec, jsc, jec, input, ocean_grid

end subroutine State_SetExport

!> This subroutine writes the minimum, maximum and sum of each field
!! contained within an ESMF state.
subroutine state_diagnose(State, string, rc)

! ----------------------------------------------
! Diagnose status of State
! ----------------------------------------------

type(ESMF_State), intent(in) :: state
character(len=*), intent(in) :: string
integer , intent(out) :: rc
type(ESMF_State), intent(in) :: state !< An ESMF State
character(len=*), intent(in) :: string !< A string indicating whether the State is an
!! import or export State
integer , intent(out) :: rc !< Return code

! local variables
integer :: i,j,n
Expand All @@ -787,19 +791,19 @@ subroutine state_diagnose(State, string, rc)
! ----------------------------------------------

call ESMF_StateGet(state, itemCount=fieldCount, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return
allocate(lfieldnamelist(fieldCount))

call ESMF_StateGet(state, itemNameList=lfieldnamelist, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return

do n = 1, fieldCount

call ESMF_StateGet(state, itemName=lfieldnamelist(n), field=lfield, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return

call field_getfldptr(lfield, fldptr1=dataPtr1d, fldptr2=dataPtr2d, rank=lrank, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return

if (lrank == 0) then
! no local data
Expand Down Expand Up @@ -829,23 +833,16 @@ subroutine state_diagnose(State, string, rc)

end subroutine state_diagnose

!===============================================================================

!> Obtain a pointer to a rank 1 or rank 2 ESMF field
subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)

! ----------------------------------------------
! for a field, determine rank and return fldptr1 or fldptr2
! abort is true by default and will abort if fldptr is not yet allocated in field
! rank returns 0, 1, or 2. 0 means fldptr not allocated and abort=false
! ----------------------------------------------

! input/output variables
type(ESMF_Field) , intent(in) :: field
real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr1(:)
real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr2(:,:)
integer , intent(out) , optional :: rank
logical , intent(in) , optional :: abort
integer , intent(out) , optional :: rc
type(ESMF_Field) , intent(in) :: field !< An ESMF field
real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr1(:) !< A pointer to a rank 1 ESMF field
real(ESMF_KIND_R8), pointer , intent(inout), optional :: fldptr2(:,:) !< A pointer to a rank 2 ESMF field
integer , intent(out) , optional :: rank !< Field rank
logical , intent(in) , optional :: abort !< Abort code
integer , intent(out) , optional :: rc !< Return code

! local variables
type(ESMF_GeomType_Flag) :: geomtype
Expand All @@ -872,7 +869,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
lrank = -99

call ESMF_FieldGet(field, status=status, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return

if (status /= ESMF_FIELDSTATUS_COMPLETE) then
lrank = 0
Expand All @@ -886,20 +883,20 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
else

call ESMF_FieldGet(field, geomtype=geomtype, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return

if (geomtype == ESMF_GEOMTYPE_GRID) then
call ESMF_FieldGet(field, rank=lrank, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return
elseif (geomtype == ESMF_GEOMTYPE_MESH) then
call ESMF_FieldGet(field, rank=lrank, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return
call ESMF_FieldGet(field, mesh=lmesh, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return
call ESMF_MeshGet(lmesh, numOwnedNodes=nnodes, numOwnedElements=nelements, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return
if (nnodes == 0 .and. nelements == 0) lrank = 0
else
else
call ESMF_LogWrite(trim(subname)//": ERROR geomtype not supported ", &
ESMF_LOGMSG_INFO)
rc = ESMF_FAILURE
Expand All @@ -917,7 +914,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
return
endif
call ESMF_FieldGet(field, farrayPtr=fldptr1, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return
elseif (lrank == 2) then
if (.not.present(fldptr2)) then
call ESMF_LogWrite(trim(subname)//": ERROR missing rank=2 array ", &
Expand All @@ -926,7 +923,7 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)
return
endif
call ESMF_FieldGet(field, farrayPtr=fldptr2, rc=rc)
if (chkerr(rc,__LINE__,u_FILE_u)) return
if (ChkErr(rc,__LINE__,u_FILE_u)) return
else
call ESMF_LogWrite(trim(subname)//": ERROR in rank ", &
ESMF_LOGMSG_ERROR, line=__LINE__, file=u_FILE_u)
Expand All @@ -942,16 +939,17 @@ subroutine field_getfldptr(field, fldptr1, fldptr2, rank, abort, rc)

end subroutine field_getfldptr

logical function chkerr(rc, line, file)
integer, intent(in) :: rc
integer, intent(in) :: line
character(len=*), intent(in) :: file
!> Returns true if ESMF_LogFoundError() determines that rc is an error code. Otherwise false.
logical function ChkErr(rc, line, file)
integer, intent(in) :: rc !< return code to check
integer, intent(in) :: line !< Integer source line number
character(len=*), intent(in) :: file !< User-provided source file name
integer :: lrc
chkerr = .false.
ChkErr = .false.
lrc = rc
if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then
chkerr = .true.
ChkErr = .true.
endif
end function chkerr
end function ChkErr

end module MOM_cap_methods
13 changes: 1 addition & 12 deletions config_src/nuopc_driver/mom_cap_time.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ module MOM_cap_time
use ESMF , only : ESMF_RC_ARG_BAD
use ESMF , only : operator(<), operator(/=), operator(+), operator(-), operator(*) , operator(>=)
use ESMF , only : operator(<=), operator(>), operator(==)
use MOM_cap_methods , only : ChkErr

implicit none; private

Expand Down Expand Up @@ -336,16 +337,4 @@ subroutine date2ymd (date, year, month, day)

end subroutine date2ymd

logical function chkerr(rc, line, file)
integer, intent(in) :: rc
integer, intent(in) :: line
character(len=*), intent(in) :: file
integer :: lrc
chkerr = .false.
lrc = rc
if (ESMF_LogFoundError(rcToCheck=lrc, msg=ESMF_LOGERR_PASSTHRU, line=line, file=file)) then
chkerr = .true.
endif
end function chkerr

end module MOM_cap_time
30 changes: 18 additions & 12 deletions config_src/nuopc_driver/mom_surface_forcing_nuopc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -183,11 +183,15 @@ module MOM_surface_forcing_nuopc
!! ice-shelves, expressed as a coefficient
!! for divergence damping, as determined
!! outside of the ocean model in [m3/s]
real, pointer, dimension(:,:) :: ustk0 => NULL() !<
real, pointer, dimension(:,:) :: vstk0 => NULL() !<
real, pointer, dimension(:) :: stk_wavenumbers => NULL() !<
real, pointer, dimension(:,:,:) :: ustkb => NULL() !<
real, pointer, dimension(:,:,:) :: vstkb => NULL() !<
real, pointer, dimension(:,:) :: ustk0 => NULL() !< Surface Stokes drift, zonal [m/s]
real, pointer, dimension(:,:) :: vstk0 => NULL() !< Surface Stokes drift, meridional [m/s]
real, pointer, dimension(:) :: stk_wavenumbers => NULL() !< The central wave number of Stokes bands [rad/m]
real, pointer, dimension(:,:,:) :: ustkb => NULL() !< Stokes Drift spectrum, zonal [m/s]
!! Horizontal - u points
!! 3rd dimension - wavenumber
real, pointer, dimension(:,:,:) :: vstkb => NULL() !< Stokes Drift spectrum, meridional [m/s]
!! Horizontal - v points
!! 3rd dimension - wavenumber
integer :: num_stk_bands !< Number of Stokes drift bands passed through the coupler
integer :: xtype !< The type of the exchange - REGRID, REDIST or DIRECT
type(coupler_2d_bc_type) :: fluxes !< A structure that may contain an array of
Expand Down Expand Up @@ -493,15 +497,17 @@ subroutine convert_IOB_to_fluxes(IOB, fluxes, index_bounds, Time, valid_time, G,
fluxes%seaice_melt(i,j) = kg_m2_s_conversion * G%mask2dT(i,j) * IOB%seaice_melt(i-i0,j-j0)

fluxes%latent(i,j) = 0.0
! notice minus sign since fprec is positive into the ocean
if (associated(IOB%fprec)) then
fluxes%latent(i,j) = fluxes%latent(i,j) + &
fluxes%latent(i,j) = fluxes%latent(i,j) - &
IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
fluxes%latent_fprec_diag(i,j) = G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
fluxes%latent_fprec_diag(i,j) = - G%mask2dT(i,j) * IOB%fprec(i-i0,j-j0)*US%W_m2_to_QRZ_T*CS%latent_heat_fusion
endif
! notice minus sign since frunoff is positive into the ocean
if (associated(IOB%frunoff)) then
fluxes%latent(i,j) = fluxes%latent(i,j) + &
fluxes%latent(i,j) = fluxes%latent(i,j) - &
IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion
fluxes%latent_frunoff_diag(i,j) = G%mask2dT(i,j) * &
fluxes%latent_frunoff_diag(i,j) = - G%mask2dT(i,j) * &
IOB%frunoff(i-i0,j-j0) * US%W_m2_to_QRZ_T * CS%latent_heat_fusion
endif
if (associated(IOB%q_flux)) then
Expand Down Expand Up @@ -790,7 +796,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
endif
forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0*tau_mag)
enddo ; enddo

call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1)
elseif (wind_stagger == AGRID) then
call pass_vector(taux_at_h, tauy_at_h, G%Domain, To_All+Omit_Corners, stagger=AGRID, halo=1)

Expand All @@ -816,7 +822,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
forces%ustar(i,j) = sqrt(gustiness*Irho0 + Irho0 * G%mask2dT(i,j) * &
sqrt(taux_at_h(i,j)**2 + tauy_at_h(i,j)**2))
enddo ; enddo

call pass_vector(forces%taux, forces%tauy, G%Domain, halo=1)
else ! C-grid wind stresses.
if (G%symmetric) &
call fill_symmetric_edges(forces%taux, forces%tauy, G%Domain)
Expand Down Expand Up @@ -858,7 +864,7 @@ subroutine convert_IOB_to_forces(IOB, forces, index_bounds, Time, G, US, CS)
enddo; enddo
call pass_vector(forces%ustkb(:,:,istk),forces%vstkb(:,:,istk), G%domain )
enddo
endif
endif

! sea ice related dynamic fields
if (associated(IOB%ice_rigidity)) then
Expand Down
Loading

0 comments on commit 12b3895

Please sign in to comment.