From 67e27cffb7900818d7badeed77a03d2108ae6543 Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Tue, 8 Sep 2026 14:09:13 +1000 Subject: [PATCH 1/5] Scale diabatic forcing differently under ice shelf The boundary layer used for the melt parameterisation must match the depth over which the forcing is applied. The only controls for these are global: HMIX_SFC_PROP/HMIX_UV_SFC_PROP control the boundary layer for computing the parameterisation, and MINIMUM_FORCING_DEPTH controls the depth over which the melt is applied. We introduce new optional parameters to applyBoundaryFluxesInOut() to modify IforcingDepthScale to use a different `shelf_forcing_depth` according to the fraction of the cell covered by ice shelf. --- .../vertical/MOM_diabatic_aux.F90 | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_aux.F90 b/src/parameterizations/vertical/MOM_diabatic_aux.F90 index 25c5fdef14..dd820994a7 100644 --- a/src/parameterizations/vertical/MOM_diabatic_aux.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_aux.F90 @@ -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 @@ -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 @@ -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, & @@ -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) From c4a00009d7e5758067853432e9f01820b31d9986 Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Tue, 8 Sep 2026 14:30:36 +1000 Subject: [PATCH 2/5] Read the SHELF_FORCING_DEPTH parameter This is the analogue to MINIMUM_FORCING_DEPTH, applied under ice shelves. It is only read when an ice shelf coverage field is provided to the diabatic driver initialisation routine. --- .../vertical/MOM_diabatic_driver.F90 | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 5389de3e1d..e656f7e16a 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -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 @@ -3214,7 +3216,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 @@ -3233,6 +3236,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] @@ -3380,6 +3384,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 "//& From 65c4e34f0656a87a3852965438773e9e2f668a4d Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Tue, 8 Sep 2026 14:41:22 +1000 Subject: [PATCH 3/5] Pass ice shelf fraction into diabatic routine This only affects the path when using the non-legacy ALE diabatic driver. In this case, the frac_shelf_h array is passed through, and used to apply a different forcing depth as required. --- src/core/MOM.F90 | 26 +++++++++++++++---- .../vertical/MOM_diabatic_driver.F90 | 17 +++++++++--- 2 files changed, 34 insertions(+), 9 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 5c7764c613..2bdf817716 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1611,6 +1611,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. @@ -1621,6 +1622,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) @@ -1677,8 +1680,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 @@ -3609,9 +3618,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) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e656f7e16a..4ba224ecc2 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -279,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] @@ -303,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) :: & @@ -412,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) @@ -1245,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 @@ -1268,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)) :: & @@ -1545,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) From 42334c215f87802ac2cce72fcf582a122570c413 Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Tue, 8 Sep 2026 15:12:40 +1000 Subject: [PATCH 4/5] Add shelf variants of HMIX This adds HMIX_SHELF_PROP and HMIX_UV_SHELF_PROP, which are analogues to the `SFC` parameters. These determine the boundary layer depth for near-shelf temperature, salinity and velocities, which are used for calculation of melt parameterisation. --- src/core/MOM.F90 | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 2bdf817716..62b3815597 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -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] @@ -2293,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. @@ -2597,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 "//& @@ -2608,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 "//& @@ -2903,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 From 8a93ac83036eb6e155bc68d53884e57b83bcf626 Mon Sep 17 00:00:00 2001 From: Angus Gibson Date: Tue, 8 Sep 2026 15:43:29 +1000 Subject: [PATCH 5/5] Split surface calculation between shelf and open ocean This now takes frac_shelf_h into account when computing depth_ml used for calculating the surface state. --- src/core/MOM.F90 | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 62b3815597..4facca958b 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -4034,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 @@ -4118,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 @@ -4182,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 @@ -4206,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