Skip to content

Ice shelf-specific forcing depths #78

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

Closed
angus-g wants to merge 5 commits into from
Closed

Ice shelf-specific forcing depths #78

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
67e27cf
Scale diabatic forcing differently under ice shelf
angus-g Sep 8, 2026
c4a0000
Read the SHELF_FORCING_DEPTH parameter
angus-g Sep 8, 2026
65c4e34
Pass ice shelf fraction into diabatic routine
angus-g Sep 8, 2026
42334c2
Add shelf variants of HMIX
angus-g Sep 8, 2026
8a93ac8
Split surface calculation between shelf and open ocean
angus-g Sep 8, 2026
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
65 changes: 60 additions & 5 deletions src/core/MOM.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -377,6 +377,10 @@ module MOM
!! feedback to the coupler/driver [H ~> m or kg m-2] when
!! bulk mixed layer is not used, or a negative value
!! if a bulk mixed layer is being used.
real :: Hmix_shelf !< Diagnostic surface boundary layer thickness over which
!! to average tracer properties under ice shelves [H ~> m or kg m-2].
real :: Hmix_UV_shelf !< Depth scale over which to average sub-ice shelf surface
!! flow [H ~> m or kg m-2].
logical :: check_bad_sfc_vals !< If true, scan surface state for ridiculous values.
real :: bad_val_ssh_max !< Maximum SSH before triggering bad value message [Z ~> m]
real :: bad_val_sst_max !< Maximum SST before triggering bad value message [C ~> degC]
Expand Down Expand Up @@ -1611,6 +1615,7 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &

logical :: debug_redundant ! If true, check redundant values on PE boundaries when debugging.
logical :: showCallTree
logical :: use_ice_shelf
type(group_pass_type) :: pass_T_S, pass_T_S_h, pass_uv_T_S_h
integer :: dynamics_stencil ! The computational stencil for the calculations
! in the dynamic core.
Expand All @@ -1621,6 +1626,8 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
showCallTree = callTree_showQuery()
if (showCallTree) call callTree_enter("step_MOM_thermo(), MOM.F90")
if (CS%debug) call query_debugging_checks(do_redundant=debug_redundant)
use_ice_shelf = .false.
if (associated(CS%frac_shelf_h)) use_ice_shelf = .true.

call enable_averages(dtdia, Time_end_thermo, CS%diag)

Expand Down Expand Up @@ -1677,8 +1684,14 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &

call cpu_clock_begin(id_clock_diabatic)

call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, dtdia, &
Time_end_thermo, G, GV, US, CS%diabatic_CSp, CS%stoch_CS, CS%OBC, Waves)
if (use_ice_shelf) then
call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, dtdia, &
Time_end_thermo, G, GV, US, CS%diabatic_CSp, CS%stoch_CS, CS%OBC, Waves, &
CS%frac_shelf_h)
else
call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, dtdia, &
Time_end_thermo, G, GV, US, CS%diabatic_CSp, CS%stoch_CS, CS%OBC, Waves)
endif
fluxes%fluxes_used = .true.

if (CS%stoch_CS%do_skeb) then
Expand Down Expand Up @@ -2284,6 +2297,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
type(group_pass_type) :: tmp_pass_uv_T_S_h, pass_uv_T_S_h

real :: Hmix_z, Hmix_UV_z ! Temporary variables with averaging depths [Z ~> m]
real :: Hmix_shelf_z, Hmix_UV_shelf_z
real :: HFrz_z ! Temporary variable with the melt potential depth [Z ~> m]
real :: default_val ! The default value for DTBT_RESET_PERIOD [s]
logical :: write_geom_files ! If true, write out the grid geometry files.
Expand Down Expand Up @@ -2588,6 +2602,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &

if (bulkmixedlayer) then
CS%Hmix = -1.0 ; CS%Hmix_UV = -1.0
CS%Hmix_shelf = -1.0 ; CS%Hmix_UV_shelf = -1.0
else
call get_param(param_file, "MOM", "HMIX_SFC_PROP", Hmix_z, &
"If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth "//&
Expand All @@ -2599,6 +2614,17 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
"over which to average to find surface flow properties, "//&
"SSU, SSV. A non-positive value indicates no averaging.", &
units="m", default=0.0, scale=US%m_to_Z)

call get_param(param_file, "MOM", "HMIX_SHELF_PROP", Hmix_shelf_z, &
"HMIX_SHELF_PROP is the depth over which to average "//&
"to find surface properties like "//&
"SST and SSS or density (but not surface velocities).", &
units="m", default=1.0, scale=US%m_to_Z)
call get_param(param_file, "MOM", "HMIX_UV_SHELF_PROP", Hmix_UV_shelf_z, &
"HMIX_UV_SHELF_PROP is the depth over which to average "//&
"to find surface flow properties, "//&
"SSU, SSV. A non-positive value indicates no averaging.", &
units="m", default=0.0, scale=US%m_to_Z)
endif
call get_param(param_file, "MOM", "HFREEZE", HFrz_z, &
"If HFREEZE > 0, melt potential will be computed. The actual depth "//&
Expand Down Expand Up @@ -2894,6 +2920,8 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
if (.not.bulkmixedlayer) then
CS%Hmix = (US%Z_to_m * GV%m_to_H) * Hmix_z
CS%Hmix_UV = (US%Z_to_m * GV%m_to_H) * Hmix_UV_z
CS%Hmix_shelf = (US%Z_to_m * GV%m_to_H) * Hmix_shelf_z
CS%Hmix_UV_shelf = (US%Z_to_m * GV%m_to_H) * Hmix_UV_shelf_z
endif
CS%HFrz = (US%Z_to_m * GV%m_to_H) * HFrz_z

Expand Down Expand Up @@ -3609,9 +3637,16 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
call adiabatic_driver_init(Time, G, param_file, diag, CS%diabatic_CSp, &
CS%tracer_flow_CSp)
else
call diabatic_driver_init(Time, G, GV, US, param_file, CS%use_ALE_algorithm, diag, &
CS%ADp, CS%CDp, CS%diabatic_CSp, CS%tracer_flow_CSp, &
CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%int_tide_CSp)
if (use_ice_shelf) then
call diabatic_driver_init(Time, G, GV, US, param_file, CS%use_ALE_algorithm, diag, &
CS%ADp, CS%CDp, CS%diabatic_CSp, CS%tracer_flow_CSp, &
CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%int_tide_CSp, &
CS%frac_shelf_h)
else
call diabatic_driver_init(Time, G, GV, US, param_file, CS%use_ALE_algorithm, diag, &
CS%ADp, CS%CDp, CS%diabatic_CSp, CS%tracer_flow_CSp, &
CS%sponge_CSp, CS%ALE_sponge_CSp, CS%oda_incupd_CSp, CS%int_tide_CSp)
endif
endif

CS%vertex_shear = kappa_shear_at_vertex(param_file)
Expand Down Expand Up @@ -3999,6 +4034,7 @@ subroutine extract_surface_state(CS, sfc_state_in)
logical :: use_iceshelves
character(240) :: msg
integer :: turns ! Number of quarter turns
real :: frac_shelf_u, frac_shelf_v

call callTree_enter("extract_surface_state(), MOM.F90")
G => CS%G ; G_in => CS%G_in ; GV => CS%GV ; US => CS%US
Expand Down Expand Up @@ -4083,6 +4119,9 @@ subroutine extract_surface_state(CS, sfc_state_in)
enddo

do k=1,nz ; do i=is,ie
if (use_iceshelves) then
depth_ml = CS%frac_shelf_h(i,j) * CS%Hmix_shelf + (1. - CS%frac_shelf_h(i,j)) * CS%Hmix
endif
if (depth(i) + h(i,j,k)*H_rescale < depth_ml) then
dh = h(i,j,k)*H_rescale
elseif (depth(i) < depth_ml) then
Expand Down Expand Up @@ -4147,6 +4186,14 @@ subroutine extract_surface_state(CS, sfc_state_in)
sfc_state%v(i,J) = 0.0
enddo
do k=1,nz ; do i=is,ie
if (use_iceshelves) then
frac_shelf_v = 0.0
if (G%areaT(i,j) + G%areaT(i,j+1) > 0.0) &
frac_shelf_v = ( CS%frac_shelf_h(i,j) * G%areaT(i,j) &
+ CS%frac_shelf_h(i,j+1) * G%areaT(i,j+1))&
/ (G%areaT(i,j) + G%areaT(i,j+1))
depth_ml = frac_shelf_v * CS%Hmix_UV_shelf + (1. - frac_shelf_v) * CS%Hmix_UV
endif
hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * H_rescale
if (depth(i) + hv < depth_ml) then
dh = hv
Expand All @@ -4171,6 +4218,14 @@ subroutine extract_surface_state(CS, sfc_state_in)
sfc_state%u(I,j) = 0.0
enddo
do k=1,nz ; do I=is-1,ie
if (use_iceshelves) then
frac_shelf_u = 0.0
if (G%areaT(i,j) + G%areaT(i+1,j) > 0.0) &
frac_shelf_u = ( CS%frac_shelf_h(i,j) * G%areaT(i,j) &
+ CS%frac_shelf_h(i+1,j) * G%areaT(i+1,j))&
/ (G%areaT(i+1,j) + G%areaT(i+1,j))
depth_ml = frac_shelf_u * CS%Hmix_UV_shelf + (1. - frac_shelf_u) * CS%Hmix_UV
endif
hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * H_rescale
if (depth(i) + hu < depth_ml) then
dh = hu
Expand Down
16 changes: 13 additions & 3 deletions src/parameterizations/vertical/MOM_diabatic_aux.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -688,7 +688,7 @@ end subroutine set_pen_shortwave
subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, &
aggregate_FW_forcing, evap_CFL_limit, &
minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, &
SkinBuoyFlux, MLD_h)
SkinBuoyFlux, MLD_h, frac_shelf_h, shelf_forcing_depth)
type(diabatic_aux_CS), pointer :: CS !< Control structure for diabatic_aux
type(ocean_grid_type), intent(in) :: G !< Grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
Expand Down Expand Up @@ -720,6 +720,9 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
optional, intent(out) :: SkinBuoyFlux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3].
real, dimension(:,:), &
optional, pointer :: MLD_h !< Mixed layer thickness for brine plumes [H ~> m or kg m-2]
real, dimension(SZI_(G),SZJ_(G)), &
optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim]
real, optional, intent(in) :: shelf_forcing_depth !< Flux depth for shelf regions

! Local variables
integer, parameter :: maxGroundings = 5
Expand Down Expand Up @@ -859,7 +862,8 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
!$OMP minimum_forcing_depth,evap_CFL_limit,dt,EOSdom, &
!$OMP calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho,&
!$OMP calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2, &
!$OMP EnthalpyConst,MLD_h) &
!$OMP EnthalpyConst,MLD_h,frac_shelf_h, &
!$OMP shelf_forcing_depth) &
!$OMP private(opacityBand,h2d,T2d,netMassInOut,netMassOut, &
!$OMP netHeat,netSalt,Pen_SW_bnd,fractionOfForcing, &
!$OMP IforcingDepthScale,g_conv,dSpV_dT,dSpV_dS, &
Expand Down Expand Up @@ -1086,7 +1090,13 @@ subroutine applyBoundaryFluxesInOut(CS, G, GV, US, dt, fluxes, optics, nsw, h, t
! Place forcing into this layer if this layer has nontrivial thickness.
! For layers thin relative to 1/IforcingDepthScale, then distribute
! forcing into deeper layers.
IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) )
if (present(frac_shelf_h)) then
IforcingDepthScale = 1. / max(GV%H_subroundoff, &
frac_shelf_h(i,j) * (shelf_forcing_depth - netMassOut(i)) &
+ (1. - frac_shelf_h(i,j)) * (minimum_forcing_depth - netMassOut(i)))
else
IforcingDepthScale = 1. / max(GV%H_subroundoff, minimum_forcing_depth - netMassOut(i) )
endif
! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
fractionOfForcing = min(1.0, h2d(i,k)*IforcingDepthScale)

Expand Down
32 changes: 27 additions & 5 deletions src/parameterizations/vertical/MOM_diabatic_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -161,6 +161,8 @@ module MOM_diabatic_driver
!! near the bottom [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
real :: minimum_forcing_depth !< The smallest depth over which heat and freshwater
!! fluxes are applied [H ~> m or kg m-2].
real :: shelf_forcing_depth !< The smallest depth over which heat and freshwater
!! fluxes are applied under shelves [H ~> m or kg m-2].
real :: evap_CFL_limit = 0.8 !< The largest fraction of a layer that can be
!! evaporated in one time-step [nondim].
integer :: halo_TS_diff = 0 !< The temperature, salinity and thickness halo size that
Expand Down Expand Up @@ -277,7 +279,7 @@ module MOM_diabatic_driver
!> This subroutine imposes the diapycnal mass fluxes and the
!! accompanying diapycnal advection of momentum and tracers.
subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
G, GV, US, CS, stoch_CS, OBC, Waves)
G, GV, US, CS, stoch_CS, OBC, Waves, frac_shelf_h)
type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< zonal velocity [L T-1 ~> m s-1]
Expand All @@ -301,6 +303,7 @@ subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
type(stochastic_CS), pointer :: stoch_CS !< stochastic control structure
type(ocean_OBC_type), pointer :: OBC !< Open boundaries control structure.
type(Wave_parameters_CS), pointer :: Waves !< Surface gravity waves
real, dimension(:,:), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim]

! local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
Expand Down Expand Up @@ -410,7 +413,7 @@ subroutine diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
G, GV, US, CS, stoch_CS, Waves)
elseif (CS%useALEalgorithm) then
call diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
G, GV, US, CS, stoch_CS, Waves)
G, GV, US, CS, stoch_CS, Waves, frac_shelf_h)
else
call layered_diabatic(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
G, GV, US, CS, Waves)
Expand Down Expand Up @@ -1243,7 +1246,7 @@ end subroutine diabatic_ALE_legacy
!> This subroutine imposes the diapycnal mass fluxes and the
!! accompanying diapycnal advection of momentum and tracers.
subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, &
G, GV, US, CS, stoch_CS, Waves)
G, GV, US, CS, stoch_CS, Waves, frac_shelf_h)
type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
Expand All @@ -1266,6 +1269,7 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
type(diabatic_CS), pointer :: CS !< module control structure
type(stochastic_CS), pointer :: stoch_CS !< stochastic control structure
type(Wave_parameters_CS), pointer :: Waves !< Surface gravity waves
real, dimension(:,:), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim]

! local variables
real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
Expand Down Expand Up @@ -1543,9 +1547,16 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end,
if (CS%use_energetic_PBL) then

skinbuoyflux(:,:) = 0.0
call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, &
if (present(frac_shelf_h)) then
call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, &
optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, &
CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD_h=visc%h_ML, &
frac_shelf_h=frac_shelf_h, shelf_forcing_depth=CS%shelf_forcing_depth)
else
call applyBoundaryFluxesInOut(CS%diabatic_aux_CSp, G, GV, US, dt, fluxes, CS%optics, &
optics_nbands(CS%optics), h, tv, CS%aggregate_FW_forcing, CS%evap_CFL_limit, &
CS%minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, SkinBuoyFlux=SkinBuoyFlux, MLD_h=visc%h_ML)
endif

if (CS%debug) then
call hchksum(ent_t, "after applyBoundaryFluxes ent_t", G%HI, haloshift=0, unscale=GV%H_to_MKS)
Expand Down Expand Up @@ -3214,7 +3225,8 @@ end subroutine adiabatic_driver_init
!> This routine initializes the diabatic driver module.
subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, diag, &
ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, &
ALE_sponge_CSp, oda_incupd_CSp, int_tide_CSp)
ALE_sponge_CSp, oda_incupd_CSp, int_tide_CSp, &
frac_shelf_h)
type(time_type), target :: Time !< model time
type(ocean_grid_type), intent(inout) :: G !< model grid structure
type(verticalGrid_type), intent(in) :: GV !< model vertical grid structure
Expand All @@ -3233,6 +3245,7 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
type(oda_incupd_CS), pointer :: oda_incupd_CSp !< pointer to the ocean data assimilation incremental
!! update module control structure
type(int_tide_CS), pointer :: int_tide_CSp !< pointer to the internal tide structure
real, dimension(:,:), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim]

! Local variables
real :: Kd ! A diffusivity used in the default for other tracer diffusivities [Z2 T-1 ~> m2 s-1]
Expand Down Expand Up @@ -3380,6 +3393,15 @@ subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, di
"relative to this scale, in which case the forcing tendencies "//&
"scaled down by distributing the forcing over this depth scale.", &
units="m", default=0.001, scale=GV%m_to_H)
if (present(frac_shelf_h)) then
call get_param(param_file, mdl, "SHELF_FORCING_DEPTH", CS%shelf_forcing_depth, &
"The smallest depth over which forcing can be applied under an "//&
"ice shelf. This only takes effect under ice shelves, and when "//&
"near-surface layers become thin relative to this scale, in "//&
"which case the forcing tendencies are scaled down by "//&
"distributing the forcing over this depth scale.", &
units="m", default=0.001, scale=GV%m_to_H)
endif
call get_param(param_file, mdl, "EVAP_CFL_LIMIT", CS%evap_CFL_limit, &
"The largest fraction of a layer than can be lost to forcing "//&
"(e.g. evaporation, sea-ice formation) in one time-step. The unused "//&
Expand Down
Loading