diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 5c7764c613..4facca958b 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] @@ -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. @@ -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) @@ -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 @@ -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. @@ -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 "//& @@ -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 "//& @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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 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) diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 5389de3e1d..4ba224ecc2 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 @@ -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] @@ -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) :: & @@ -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) @@ -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 @@ -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)) :: & @@ -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) @@ -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 @@ -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] @@ -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 "//&