From 249a0788d8a041cd7a4cef497aa1c5b2cf75931d Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 11 Jul 2023 17:40:32 -0400 Subject: [PATCH 1/3] +*Revise units of arguments to vert_fill_TS Pass dt_kappa_smooth to calc_isoneutral_slopes and vert_fill_TS in units of [H Z ~> m2 or kg m-1] instead of [Z2 ~> m2] for consistency with the units of other diffusivities in the code and to reduce the depenency on the Boussinesq reference density in non-Boussinesq configurations. In addition to the changes to the units of these two arguments, there is a new unit_scale_type argument to vert_fill_TS and MOM_calc_varT and a new verticalGrid_type argument to MOM_stoch_eos_init. The units of 4 vertical diffusivities in the control structures in 4 different modules are also changed accordingly. All answers are bitwise identical in Boussinesq mode, but they can change for some non-Boussinesq configurations. There are new mandatory arguments to three publicly visible routines. --- src/core/MOM.F90 | 4 +-- src/core/MOM_isopycnal_slopes.F90 | 34 +++++++++++-------- src/core/MOM_stoch_eos.F90 | 12 ++++--- src/parameterizations/lateral/MOM_MEKE.F90 | 3 +- .../lateral/MOM_lateral_mixing_coeffs.F90 | 4 +-- .../lateral/MOM_thickness_diffuse.F90 | 10 +++--- .../vertical/MOM_internal_tide_input.F90 | 6 ++-- .../vertical/MOM_set_diffusivity.F90 | 12 +++---- 8 files changed, 46 insertions(+), 39 deletions(-) diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index d2df26524d..d7e2d74735 100644 --- a/src/core/MOM.F90 +++ b/src/core/MOM.F90 @@ -1090,7 +1090,7 @@ subroutine step_MOM_dynamics(forces, p_surf_begin, p_surf_end, dt, dt_thermo, & call cpu_clock_end(id_clock_stoch) call cpu_clock_begin(id_clock_varT) if (CS%use_stochastic_EOS) then - call MOM_calc_varT(G, GV, h, CS%tv, CS%stoch_eos_CS, dt) + call MOM_calc_varT(G, GV, US, h, CS%tv, CS%stoch_eos_CS, dt) if (associated(CS%tv%varT)) call pass_var(CS%tv%varT, G%Domain, clock=id_clock_pass, halo=1) endif call cpu_clock_end(id_clock_varT) @@ -3022,7 +3022,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, & new_sim = is_new_run(restart_CSp) if (use_temperature) then - CS%use_stochastic_EOS = MOM_stoch_eos_init(Time, G, US, param_file, diag, CS%stoch_eos_CS, restart_CSp) + CS%use_stochastic_EOS = MOM_stoch_eos_init(Time, G, GV, US, param_file, diag, CS%stoch_eos_CS, restart_CSp) else CS%use_stochastic_EOS = .false. endif diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 07dd19b0a6..73ae8a9816 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -36,8 +36,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: e !< Interface heights [Z ~> m] type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various !! thermodynamic variables - real, intent(in) :: dt_kappa_smooth !< A smoothing vertical diffusivity - !! times a smoothing timescale [Z2 ~> m2]. + real, intent(in) :: dt_kappa_smooth !< A smoothing vertical + !! diffusivity times a smoothing + !! timescale [H Z ~> m2 or kg m-1] logical, intent(in) :: use_stanley !< turn on stanley param in slope real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: slope_x !< Isopycnal slope in i-dir [Z L-1 ~> nondim] real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(inout) :: slope_y !< Isopycnal slope in j-dir [Z L-1 ~> nondim] @@ -142,7 +143,7 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan h_neglect = GV%H_subroundoff ; h_neglect2 = h_neglect**2 - dz_neglect = GV%H_subroundoff * GV%H_to_Z + dz_neglect = GV%dZ_subroundoff local_open_u_BC = .false. local_open_v_BC = .false. @@ -195,9 +196,9 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan if (use_EOS) then if (present(halo)) then - call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, halo+1) + call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, US, halo+1) else - call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, 1) + call vert_fill_TS(h, tv%T, tv%S, dt_kappa_smooth, T, S, G, GV, US, 1) endif endif @@ -341,9 +342,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan slope = slope * max(g%mask2dT(i,j),g%mask2dT(i+1,j)) endif slope_x(I,j,K) = slope - if (present(dzSxN)) dzSxN(I,j,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) & - + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N - * abs(slope) * G%mask2dCu(I,j) ! x-direction contribution to S^2 + if (present(dzSxN)) & + dzSxN(I,j,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) & + + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N + * abs(slope) * G%mask2dCu(I,j) ! x-direction contribution to S^2 enddo ! I enddo ; enddo ! end of j-loop @@ -477,9 +479,10 @@ subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, use_stan slope = slope * max(g%mask2dT(i,j),g%mask2dT(i,j+1)) endif slope_y(i,J,K) = slope - if (present(dzSyN)) dzSyN(i,J,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) & - + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N - * abs(slope) * G%mask2dCv(i,J) ! x-direction contribution to S^2 + if (present(dzSyN)) & + dzSyN(i,J,K) = sqrt( G_Rho0 * max(0., wtL * ( dzaL * drdkL ) & + + wtR * ( dzaR * drdkR )) / (wtL + wtR) ) & ! dz * N + * abs(slope) * G%mask2dCv(i,J) ! x-direction contribution to S^2 enddo ! i enddo ; enddo ! end of j-loop @@ -488,14 +491,15 @@ end subroutine calc_isoneutral_slopes !> Returns tracer arrays (nominally T and S) with massless layers filled with !! sensible values, by diffusing vertically with a small but constant diffusivity. -subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, larger_h_denom) +subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, US, halo_here, larger_h_denom) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: T_in !< Input temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: S_in !< Input salinity [S ~> ppt] real, intent(in) :: kappa_dt !< A vertical diffusivity to use for smoothing - !! times a smoothing timescale [Z2 ~> m2]. + !! times a smoothing timescale [H Z ~> m2 or kg m-1] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T_f !< Filled temperature [C ~> degC] real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S_f !< Filled salinity [S ~> ppt] integer, optional, intent(in) :: halo_here !< Number of halo points to work on, @@ -525,10 +529,10 @@ subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, lar is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke h_neglect = GV%H_subroundoff - kap_dt_x2 = (2.0*kappa_dt)*GV%Z_to_H**2 + kap_dt_x2 = (2.0*kappa_dt) * (US%Z_to_m*GV%m_to_H) ! Usually the latter term is GV%Z_to_H. h0 = h_neglect if (present(larger_h_denom)) then - if (larger_h_denom) h0 = 1.0e-16*sqrt(kappa_dt)*GV%Z_to_H + if (larger_h_denom) h0 = 1.0e-16*sqrt(0.5*kap_dt_x2) endif if (kap_dt_x2 <= 0.0) then diff --git a/src/core/MOM_stoch_eos.F90 b/src/core/MOM_stoch_eos.F90 index deb878e99c..2bd742be6d 100644 --- a/src/core/MOM_stoch_eos.F90 +++ b/src/core/MOM_stoch_eos.F90 @@ -40,7 +40,7 @@ module MOM_stoch_eos real :: stanley_coeff !< Coefficient correlating the temperature gradient !! and SGS T variance [nondim]; if <0, turn off scheme in all codes real :: stanley_a !< a in exp(aX) in stochastic coefficient [nondim] - real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1] + real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] !>@{ Diagnostic IDs integer :: id_stoch_eos = -1, id_stoch_phi = -1, id_tvar_sgs = -1 @@ -51,9 +51,10 @@ module MOM_stoch_eos contains !> Initializes MOM_stoch_eos module, returning a logical indicating whether this module will be used. -logical function MOM_stoch_eos_init(Time, G, US, param_file, diag, CS, restart_CS) +logical function MOM_stoch_eos_init(Time, G, GV, US, param_file, diag, CS, restart_CS) type(time_type), intent(in) :: Time !< Time for stochastic process type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type type(param_file_type), intent(in) :: param_file !< structure indicating parameter file to parse type(diag_ctrl), target, intent(inout) :: diag !< Structure used to control diagnostics @@ -80,7 +81,7 @@ logical function MOM_stoch_eos_init(Time, G, US, param_file, diag, CS, restart_C call get_param(param_file, "MOM_stoch_eos", "KD_SMOOTH", CS%kappa_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m_to_Z**2*US%T_to_s, & + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T, & do_not_log=(CS%stanley_coeff<0.0)) ! Don't run anything if STANLEY_COEFF < 0 @@ -193,9 +194,10 @@ subroutine post_stoch_EOS_diags(CS, tv, diag) end subroutine post_stoch_EOS_diags !> Computes a parameterization of the SGS temperature variance -subroutine MOM_calc_varT(G, GV, h, tv, CS, dt) +subroutine MOM_calc_varT(G, GV, US, h, tv, CS, dt) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type real, dimension(SZI_(G),SZJ_(G),SZK_(G)), & intent(in) :: h !< Layer thickness [H ~> m] type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure @@ -219,7 +221,7 @@ subroutine MOM_calc_varT(G, GV, h, tv, CS, dt) ! extreme gradients along layers which are vanished against topography. It is ! still a poor approximation in the interior when coordinates are strongly tilted. if (.not. associated(tv%varT)) allocate(tv%varT(G%isd:G%ied, G%jsd:G%jed, GV%ke), source=0.0) - call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, halo_here=1, larger_h_denom=.true.) + call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, US, halo_here=1, larger_h_denom=.true.) do k=1,G%ke do j=G%jsc,G%jec diff --git a/src/parameterizations/lateral/MOM_MEKE.F90 b/src/parameterizations/lateral/MOM_MEKE.F90 index 02338fab96..0ef261a956 100644 --- a/src/parameterizations/lateral/MOM_MEKE.F90 +++ b/src/parameterizations/lateral/MOM_MEKE.F90 @@ -1578,7 +1578,8 @@ subroutine ML_MEKE_calculate_features(G, GV, US, CS, Rd_dx_h, u, v, tv, h, dt, f h_v(i,J,k) = 0.5*(h(i,j,k)*G%mask2dT(i,j) + h(i,j+1,k)*G%mask2dT(i,j+1)) + GV%Angstrom_H enddo; enddo; enddo; call find_eta(h, tv, G, GV, US, e, halo_size=2) - call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*1.e-7, .false., slope_x, slope_y) + ! Note the hard-coded dimenisional constant in the following line. + call calc_isoneutral_slopes(G, GV, US, h, e, tv, dt*1.e-7*GV%m2_s_to_HZ_T, .false., slope_x, slope_y) call pass_vector(slope_x, slope_y, G%Domain) do j=js-1,je+1; do i=is-1,ie+1 slope_x_vert_avg(I,j) = vertical_average_interface(slope_x(i,j,:), h_u(i,j,:), GV%H_subroundoff) diff --git a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 index 8f0aa02b12..e74e48055a 100644 --- a/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 +++ b/src/parameterizations/lateral/MOM_lateral_mixing_coeffs.F90 @@ -135,7 +135,7 @@ module MOM_lateral_mixing_coeffs !! F = 1 / (1 + (Res_coef_visc*Ld/dx)^Res_fn_power) real :: depth_scaled_khth_h0 !< The depth above which KHTH is linearly scaled away [Z ~> m] real :: depth_scaled_khth_exp !< The exponent used in the depth dependent scaling function for KHTH [nondim] - real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [Z2 T-1 ~> m2 s-1] + real :: kappa_smooth !< A diffusivity for smoothing T/S in vanished layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] integer :: Res_fn_power_khth !< The power of dx/Ld in the KhTh resolution function. Any !! positive integer power may be used, but even powers !! and especially 2 are coded to be more efficient. @@ -1265,7 +1265,7 @@ subroutine VarMix_init(Time, G, GV, US, param_file, diag, CS) call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m_to_Z**2*US%T_to_s) + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T) endif if (CS%calculate_Eady_growth_rate) then diff --git a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 index 8617795e16..f24f790d06 100644 --- a/src/parameterizations/lateral/MOM_thickness_diffuse.F90 +++ b/src/parameterizations/lateral/MOM_thickness_diffuse.F90 @@ -44,9 +44,9 @@ module MOM_thickness_diffuse real :: Kh_eta_bg !< Background isopycnal height diffusivity [L2 T-1 ~> m2 s-1] real :: Kh_eta_vel !< Velocity scale that is multiplied by the grid spacing to give !! the isopycnal height diffusivity [L T-1 ~> m s-1] - real :: slope_max !< Slopes steeper than slope_max are limited in some way [Z L-1 ~> nondim]. - real :: kappa_smooth !< Vertical diffusivity used to interpolate more - !! sensible values of T & S into thin layers [Z2 T-1 ~> m2 s-1]. + real :: slope_max !< Slopes steeper than slope_max are limited in some way [Z L-1 ~> nondim] + real :: kappa_smooth !< Vertical diffusivity used to interpolate more sensible values + !! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] logical :: thickness_diffuse !< If true, interfaces heights are diffused. logical :: use_FGNV_streamfn !< If true, use the streamfunction formulation of !! Ferrari et al., 2010, which effectively emphasizes @@ -798,7 +798,7 @@ subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV if (use_EOS) then halo = 1 ! Default halo to fill is 1 - call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, halo, larger_h_denom=.true.) + call vert_fill_TS(h, tv%T, tv%S, CS%kappa_smooth*dt, T, S, G, GV, US, halo, larger_h_denom=.true.) endif ! Rescale the thicknesses, perhaps using the specific volume. @@ -2191,7 +2191,7 @@ subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS) call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m_to_Z**2*US%T_to_s) + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T) call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", CS%use_FGNV_streamfn, & "If true, use the streamfunction formulation of "//& "Ferrari et al., 2010, which effectively emphasizes "//& diff --git a/src/parameterizations/vertical/MOM_internal_tide_input.F90 b/src/parameterizations/vertical/MOM_internal_tide_input.F90 index 7ec612f141..95e33929df 100644 --- a/src/parameterizations/vertical/MOM_internal_tide_input.F90 +++ b/src/parameterizations/vertical/MOM_internal_tide_input.F90 @@ -41,7 +41,7 @@ module MOM_int_tide_input real :: TKE_itide_max !< Maximum Internal tide conversion !! available to mix above the BBL [R Z3 T-3 ~> W m-2] real :: kappa_fill !< Vertical diffusivity used to interpolate sensible values - !! of T & S into thin layers [Z2 T-1 ~> m2 s-1]. + !! of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] real, allocatable, dimension(:,:) :: TKE_itidal_coef !< The time-invariant field that enters the TKE_itidal input calculation [R Z3 T-2 ~> J m-2]. @@ -118,7 +118,7 @@ subroutine set_int_tide_input(u, v, h, tv, fluxes, itide, dt, G, GV, US, CS) ! Smooth the properties through massless layers. if (use_EOS) then - call vert_fill_TS(h, tv%T, tv%S, CS%kappa_fill*dt, T_f, S_f, G, GV, larger_h_denom=.true.) + call vert_fill_TS(h, tv%T, tv%S, CS%kappa_fill*dt, T_f, S_f, G, GV, US, larger_h_denom=.true.) endif call find_N2_bottom(h, tv, T_f, S_f, itide%h2, fluxes, G, GV, US, N2_bot) @@ -352,7 +352,7 @@ subroutine int_tide_input_init(Time, G, GV, US, param_file, diag, CS, itide) call get_param(param_file, mdl, "KD_SMOOTH", CS%kappa_fill, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m2_s_to_Z2_T) + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T) call get_param(param_file, mdl, "UTIDE", utide, & "The constant tidal amplitude used with INT_TIDE_DISSIPATION.", & diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index 9dc7b81c46..dfd264c92a 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -87,7 +87,7 @@ module MOM_set_diffusivity real :: Kd_add !< uniform diffusivity added everywhere without !! filtering or scaling [Z2 T-1 ~> m2 s-1]. real :: Kd_smooth !< Vertical diffusivity used to interpolate more - !! sensible values of T & S into thin layers [Z2 T-1 ~> m2 s-1]. + !! sensible values of T & S into thin layers [H Z T-1 ~> m2 s-1 or kg m-1 s-1] type(diag_ctrl), pointer :: diag => NULL() !< structure to regulate diagnostic output timing logical :: limit_dissipation !< If enabled, dissipation is limited to be larger @@ -274,7 +274,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i integer :: i, j, k, is, ie, js, je, nz, isd, ied, jsd, jed - real :: kappa_dt_fill ! diffusivity times a timestep used to fill massless layers [Z2 ~> m2] + real :: kappa_dt_fill ! diffusivity times a timestep used to fill massless layers [H Z ~> m2 or kg m-1] is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -289,7 +289,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i if (CS%answer_date < 20190101) then ! These hard-coded dimensional parameters are being replaced. - kappa_dt_fill = US%m_to_Z**2 * 1.e-3 * 7200. + kappa_dt_fill = 1.e-3*GV%m2_s_to_HZ_T * 7200.*US%s_to_T else kappa_dt_fill = CS%Kd_smooth * dt endif @@ -340,7 +340,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i call cpu_clock_begin(id_clock_kappaShear) if (CS%Vertex_shear) then call full_convection(G, GV, US, h, tv, T_f, S_f, fluxes%p_surf, & - (GV%Z_to_H**2)*kappa_dt_fill, halo=1) + GV%Z_to_H*kappa_dt_fill, halo=1) call calc_kappa_shear_vertex(u, v, h, T_f, S_f, tv, fluxes%p_surf, visc%Kd_shear, & visc%TKE_turb, visc%Kv_shear_Bu, dt, G, GV, US, CS%kappaShear_CSp) @@ -380,7 +380,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i call hchksum(tv%S, "before vert_fill_TS tv%S", G%HI, scale=US%S_to_ppt) call hchksum(h, "before vert_fill_TS h",G%HI, scale=GV%H_to_m) endif - call vert_fill_TS(h, tv%T, tv%S, kappa_dt_fill, T_f, S_f, G, GV, larger_h_denom=.true.) + call vert_fill_TS(h, tv%T, tv%S, kappa_dt_fill, T_f, S_f, G, GV, US, larger_h_denom=.true.) if (CS%debug) then call hchksum(tv%T, "after vert_fill_TS tv%T", G%HI, scale=US%C_to_degC) call hchksum(tv%S, "after vert_fill_TS tv%S", G%HI, scale=US%S_to_ppt) @@ -2212,7 +2212,7 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ call get_param(param_file, mdl, "KD_SMOOTH", CS%Kd_smooth, & "A diapycnal diffusivity that is used to interpolate "//& "more sensible values of T & S into thin layers.", & - units="m2 s-1", default=1.0e-6, scale=US%m2_s_to_Z2_T) + units="m2 s-1", default=1.0e-6, scale=GV%m2_s_to_HZ_T) call get_param(param_file, mdl, "DEBUG", CS%debug, & "If true, write out verbose debugging data.", & From e465b1f6f913c97905b9d8aa9d9aec1110154a01 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Wed, 26 Jul 2023 04:55:19 -0400 Subject: [PATCH 2/3] Add comment justifying rescaling in vert_fill_TS Added a comment justifying the use of a fixed rescaling factor for the diffusivity used in vert_fill_TS. All answers and output are identical. --- src/core/MOM_isopycnal_slopes.F90 | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/core/MOM_isopycnal_slopes.F90 b/src/core/MOM_isopycnal_slopes.F90 index 73ae8a9816..29c547148d 100644 --- a/src/core/MOM_isopycnal_slopes.F90 +++ b/src/core/MOM_isopycnal_slopes.F90 @@ -529,6 +529,11 @@ subroutine vert_fill_TS(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, US, halo_here, is = G%isc-halo ; ie = G%iec+halo ; js = G%jsc-halo ; je = G%jec+halo ; nz = GV%ke h_neglect = GV%H_subroundoff + ! The use of the fixed rescaling factor in the next line avoids an extra call to thickness_to_dz() + ! and the use of an extra 3-d array of vertical distnaces across layers (dz). This would be more + ! physically consistent, but it would also be more expensive, and given that this routine applies + ! a small (but arbitrary) amount of mixing to clean up the properties of nearly massless layers, + ! the added expense is hard to justify. kap_dt_x2 = (2.0*kappa_dt) * (US%Z_to_m*GV%m_to_H) ! Usually the latter term is GV%Z_to_H. h0 = h_neglect if (present(larger_h_denom)) then From 636d6109b9c14f7c5479dc5514913159d7cd97e3 Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Tue, 18 Jul 2023 05:27:01 -0400 Subject: [PATCH 3/3] +*Add and use find_ustar Added the new public interface find_ustar to extract the friction velocity from either a forcing type argument, or a mech_forcing_type argument, either directly or from tau_mag, and in non-Boussinesq mode by using the time-evolving surface specific volume. Find_ustar is an overloaded interface to find_ustar_fluxes or find_ustar_mech_forcing, which are the same but for the type of one of their arguments. For now, the subroutines bulkmixedlayer, mixedlayer_restrajt_OM4, mixedlayer_restrat_Bodner and mixedlayer_restrat_BML are calling find_ustar to avoid code duplication during the transition to work in fully non-Boussinesq mode, but it will eventually be used in about another half dozen other places. All Boussinesq answers are bitwise identical, but non-Boussinesq answers will change and become less dependent on the Boussinesq reference density, and there is a new publicly visible interface wrapping two subroutines. --- src/core/MOM_forcing_type.F90 | 141 +++++++++++++++++- .../lateral/MOM_mixed_layer_restrat.F90 | 36 +++-- .../vertical/MOM_bulk_mixed_layer.F90 | 19 ++- 3 files changed, 182 insertions(+), 14 deletions(-) diff --git a/src/core/MOM_forcing_type.F90 b/src/core/MOM_forcing_type.F90 index 9623256b88..a6d35903ee 100644 --- a/src/core/MOM_forcing_type.F90 +++ b/src/core/MOM_forcing_type.F90 @@ -29,7 +29,7 @@ module MOM_forcing_type public extractFluxes1d, extractFluxes2d, optics_type public MOM_forcing_chksum, MOM_mech_forcing_chksum -public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d +public calculateBuoyancyFlux1d, calculateBuoyancyFlux2d, find_ustar public forcing_accumulate, fluxes_accumulate public forcing_SinglePointPrint, mech_forcing_diags, forcing_diagnostics public register_forcing_type_diags, allocate_forcing_type, deallocate_forcing_type @@ -53,6 +53,12 @@ module MOM_forcing_type module procedure allocate_mech_forcing_from_ref end interface allocate_mech_forcing +!> Determine the friction velocity from a forcing type or a mechanical forcing type. +interface find_ustar + module procedure find_ustar_fluxes + module procedure find_ustar_mech_forcing +end interface find_ustar + ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units @@ -1077,6 +1083,139 @@ subroutine calculateBuoyancyFlux2d(G, GV, US, fluxes, optics, h, Temp, Salt, tv, end subroutine calculateBuoyancyFlux2d +!> Determine the friction velocity from the contenxts of a forcing type, perhaps +!! using the evolving surface density. +subroutine find_ustar_fluxes(fluxes, tv, U_star, G, GV, US, halo, H_T_units) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(forcing), intent(in) :: fluxes !< Surface fluxes container + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any + !! available thermodynamic fields. + real, dimension(SZI_(G),SZJ_(G)), & + intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1] or + !! [H T-1 ~> m s-1 or kg m-2 s-1], depending on H_T_units. + integer, optional, intent(in) :: halo !< The extra halo size to fill in, 0 by default + logical, optional, intent(in) :: H_T_units !< If present and true, return U_star in units + !! of [H T-1 ~> m s-1 or kg m-2 s-1] + + ! Local variables + real :: I_rho ! The inverse of the reference density times a ratio of scaling + ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases + ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3] + logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is + ! returned in [H T-1 ~> m s-1 or kg m-2 s-1] + integer :: i, j, k, is, ie, js, je, hs + + hs = 0 ; if (present(halo)) hs = max(halo, 0) + is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs + + Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units + + if (.not.(associated(fluxes%ustar) .or. associated(fluxes%tau_mag))) & + call MOM_error(FATAL, "find_ustar_fluxes requires that either ustar or tau_mag be associated.") + + if (associated(fluxes%ustar) .and. (GV%Boussinesq .or. .not.associated(fluxes%tau_mag))) then + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = fluxes%ustar(i,j) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%Z_to_H * fluxes%ustar(i,j) + enddo ; enddo + endif + elseif (allocated(tv%SpV_avg)) then + if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, & + "find_ustar_fluxes called in non-Boussinesq mode with invalid values of SpV_avg.") + if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, & + "find_ustar_fluxes called in non-Boussinesq mode with insufficient valid values of SpV_avg.") + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(US%L_to_Z*fluxes%tau_mag(i,j) * tv%SpV_avg(i,j,1)) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*fluxes%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + enddo ; enddo + endif + else + I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H + if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(fluxes%tau_mag(i,j) * I_rho) + enddo ; enddo + endif + +end subroutine find_ustar_fluxes + + +!> Determine the friction velocity from the contenxts of a forcing type, perhaps +!! using the evolving surface density. +subroutine find_ustar_mech_forcing(forces, tv, U_star, G, GV, US, halo, H_T_units) + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(mech_forcing), intent(in) :: forces !< Surface forces container + type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any + !! available thermodynamic fields. + real, dimension(SZI_(G),SZJ_(G)), & + intent(out) :: U_star !< The surface friction velocity [Z T-1 ~> m s-1] + integer, optional, intent(in) :: halo !< The extra halo size to fill in, 0 by default + logical, optional, intent(in) :: H_T_units !< If present and true, return U_star in units + !! of [H T-1 ~> m s-1 or kg m-2 s-1] + + ! Local variables + real :: I_rho ! The inverse of the reference density times a ratio of scaling + ! factors [Z L-1 R-1 ~> m3 kg-1] or in some semi-Boussinesq cases + ! the rescaled reference density [H2 Z-1 L-1 R-1 ~> m3 kg-1 or kg m-3] + logical :: Z_T_units ! If true, U_star is returned in units of [Z T-1 ~> m s-1], otherwise it is + ! returned in [H T-1 ~> m s-1 or kg m-2 s-1] + integer :: i, j, k, is, ie, js, je, hs + + hs = 0 ; if (present(halo)) hs = max(halo, 0) + is = G%isc - hs ; ie = G%iec + hs ; js = G%jsc - hs ; je = G%jec + hs + + Z_T_units = .true. ; if (present(H_T_units)) Z_T_units = .not.H_T_units + + if (.not.(associated(forces%ustar) .or. associated(forces%tau_mag))) & + call MOM_error(FATAL, "find_ustar_mech requires that either ustar or tau_mag be associated.") + + if (associated(forces%ustar) .and. (GV%Boussinesq .or. .not.associated(forces%tau_mag))) then + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = forces%ustar(i,j) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%Z_to_H * forces%ustar(i,j) + enddo ; enddo + endif + elseif (allocated(tv%SpV_avg)) then + if (tv%valid_SpV_halo < 0) call MOM_error(FATAL, & + "find_ustar_mech called in non-Boussinesq mode with invalid values of SpV_avg.") + if (tv%valid_SpV_halo < hs) call MOM_error(FATAL, & + "find_ustar_mech called in non-Boussinesq mode with insufficient valid values of SpV_avg.") + if (Z_T_units) then + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(US%L_to_Z*forces%tau_mag(i,j) * tv%SpV_avg(i,j,1)) + enddo ; enddo + else + do j=js,je ; do i=is,ie + U_star(i,j) = GV%RZ_to_H * sqrt(US%L_to_Z*forces%tau_mag(i,j) / tv%SpV_avg(i,j,1)) + enddo ; enddo + endif + else + I_rho = US%L_to_Z * GV%Z_to_H * GV%RZ_to_H + if (Z_T_units) I_rho = US%L_to_Z * GV%H_to_Z * GV%RZ_to_H ! == US%L_to_Z / GV%Rho0 + do j=js,je ; do i=is,ie + U_star(i,j) = sqrt(forces%tau_mag(i,j) * I_rho) + enddo ; enddo + endif + +end subroutine find_ustar_mech_forcing + + !> Write out chksums for thermodynamic fluxes. subroutine MOM_forcing_chksum(mesg, fluxes, G, US, haloshift) character(len=*), intent(in) :: mesg !< message diff --git a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 index 206773ecb0..5b7ec60dee 100644 --- a/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 +++ b/src/parameterizations/lateral/MOM_mixed_layer_restrat.F90 @@ -11,7 +11,7 @@ module MOM_mixed_layer_restrat use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_version, param_file_type use MOM_file_parser, only : openParameterBlock, closeParameterBlock -use MOM_forcing_type, only : mech_forcing +use MOM_forcing_type, only : mech_forcing, find_ustar use MOM_grid, only : ocean_grid_type use MOM_hor_index, only : hor_index_type use MOM_lateral_mixing_coeffs, only : VarMix_CS @@ -184,6 +184,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, h_avail ! The volume available for diffusion out of each face of each ! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZI_(G),SZJ_(G)) :: & + U_star_2d, & ! The wind friction velocity, calculated using + ! the Boussinesq reference density or the time-evolving surface density + ! in non-Boussinesq mode [Z T-1 ~> m s-1] MLD_fast, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2] htot_fast, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2] Rml_av_fast, & ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2] @@ -254,6 +257,9 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, call MOM_error(FATAL, "mixedlayer_restrat_OM4: "// & "The resolution argument, Rd/dx, was not associated.") + ! Extract the friction velocity from the forcing type. + call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1) + if (CS%MLE_density_diff > 0.) then ! We need to calculate a mixed layer depth, MLD. !! TODO: use derivatives and mid-MLD pressure. Currently this is sigma-0. -AJA pRef_MLD(:) = 0. @@ -408,7 +414,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, if (CS%debug) then call hchksum(h,'mixed_layer_restrat: h', G%HI, haloshift=1, scale=GV%H_to_m) - call hchksum(forces%ustar,'mixed_layer_restrat: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) + call hchksum(U_star_2d, 'mixed_layer_restrat: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) call hchksum(MLD_fast,'mixed_layer_restrat: MLD', G%HI, haloshift=1, scale=GV%H_to_m) call hchksum(Rml_av_fast,'mixed_layer_restrat: rml', G%HI, haloshift=1, & scale=US%m_to_Z*US%L_T_to_m_s**2) @@ -421,7 +427,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, ! U - Component !$OMP do do j=js,je ; do I=is-1,ie - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j))) absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) ! If needed, res_scaling_fac = min( ds, L_d ) / l_f @@ -508,7 +514,7 @@ subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, MLD_in, VarMix, ! V- component !$OMP do do J=js-1,je ; do i=is,ie - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1))) absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) ! If needed, res_scaling_fac = min( ds, L_d ) / l_f @@ -711,6 +717,9 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d wpup ! Turbulent vertical momentum [ ????? ~> m2 s-2] real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1] real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1] + real :: U_star_2d(SZI_(G),SZJ_(G)) ! The wind friction velocity, calculated using the Boussinesq + ! reference density or the time-evolving surface density in non-Boussinesq + ! mode [Z T-1 ~> m s-1] real :: covTS(SZI_(G)) ! SGS TS covariance in Stanley param; currently 0 [degC ppt] real :: varS(SZI_(G)) ! SGS S variance in Stanley param; currently 0 [ppt2] real :: dmu(SZK_(GV)) ! Change in mu(z) across layer k [nondim] @@ -762,12 +771,15 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d call pass_var(bflux, G%domain, halo=1) + ! Extract the friction velocity from the forcing type. + call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1) + if (CS%debug) then call hchksum(h,'mixed_Bodner: h', G%HI, haloshift=1, scale=GV%H_to_m) call hchksum(BLD, 'mle_Bodner: BLD in', G%HI, haloshift=1, scale=US%Z_to_m) if (associated(bflux)) & call hchksum(bflux, 'mle_Bodner: bflux', G%HI, haloshift=1, scale=US%Z_to_m**2*US%s_to_T**3) - call hchksum(forces%ustar,'mle_Bodner: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) + call hchksum(U_star_2d, 'mle_Bodner: u*', G%HI, haloshift=1, scale=US%Z_to_m*US%s_to_T) call hchksum(CS%MLD_filtered, 'mle_Bodner: MLD_filtered 1', & G%HI, haloshift=1, scale=US%Z_to_m) call hchksum(CS%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 1', & @@ -793,7 +805,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d do j = js-1, je+1 ; do i = is-1, ie+1 w_star3 = max(0., -bflux(i,j)) * BLD(i,j) & ! (this line in Z3 T-3 ~> m3 s-3) * ( ( US%Z_to_m * US%s_to_T )**3 ) ! m3 s-3 - u_star3 = ( US%Z_to_m * US%s_to_T * forces%ustar(i,j) )**3 ! m3 s-3 + u_star3 = ( US%Z_to_m * US%s_to_T * U_star_2d(i,j) )**3 ! m3 s-3 wpup(i,j) = max( CS%min_wstar2, & ! The max() avoids division by zero later ( CS%mstar * u_star3 + CS%nstar * w_star3 )**two_thirds ) & ! (this line m2 s-2) * ( ( US%m_to_Z * US%T_to_s )**2 ) ! Z2 T-2 ~> m2 s-2 @@ -965,7 +977,7 @@ subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, d ! Offer diagnostic fields for averaging. if (query_averaging_enabled(CS%diag)) then - if (CS%id_ustar > 0) call post_data(CS%id_ustar, forces%ustar, CS%diag) + if (CS%id_ustar > 0) call post_data(CS%id_ustar, U_star_2d, CS%diag) if (CS%id_bflux > 0) call post_data(CS%id_bflux, bflux, CS%diag) if (CS%id_wpup > 0) call post_data(CS%id_wpup, wpup, CS%diag) if (CS%id_Rml > 0) call post_data(CS%id_Rml, buoy_av, CS%diag) @@ -1053,6 +1065,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) h_avail ! The volume available for diffusion out of each face of each ! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1]. real, dimension(SZI_(G),SZJ_(G)) :: & + U_star_2d, & ! The wind friction velocity, calculated using + ! the Boussinesq reference density or the time-evolving surface density + ! in non-Boussinesq mode [Z T-1 ~> m s-1] htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2] Rml_av ! g_Rho0 times the average mixed layer density [L2 Z-1 T-2 ~> m s-2] real :: g_Rho0 ! G_Earth/Rho0 [L2 Z-1 T-2 R-1 ~> m4 s-2 kg-1] @@ -1110,6 +1125,9 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) if (CS%use_Stanley_ML) call MOM_error(FATAL, "mixedlayer_restrat_BML: "// & "The Stanley parameterization is not available with the BML.") + ! Extract the friction velocity from the forcing type. + call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1) + ! Fix this later for nkml >= 3. p0(:) = 0.0 @@ -1145,7 +1163,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) do j=js,je ; do I=is-1,ie h_vel = 0.5*(htot(i,j) + htot(i+1,j)) * GV%H_to_Z - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i+1,j))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j))) absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) @@ -1196,7 +1214,7 @@ subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS) do J=js-1,je ; do i=is,ie h_vel = 0.5*(htot(i,j) + htot(i,j+1)) * GV%H_to_Z - u_star = max(CS%ustar_min, 0.5*(forces%ustar(i,j) + forces%ustar(i,j+1))) + u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1))) absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) diff --git a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 index 66e2dfa6b2..ceba8dad1a 100644 --- a/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 +++ b/src/parameterizations/vertical/MOM_bulk_mixed_layer.F90 @@ -9,7 +9,7 @@ module MOM_bulk_mixed_layer use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type use MOM_error_handler, only : MOM_error, FATAL, WARNING use MOM_file_parser, only : get_param, log_param, log_version, param_file_type -use MOM_forcing_type, only : extractFluxes1d, forcing +use MOM_forcing_type, only : extractFluxes1d, forcing, find_ustar use MOM_grid, only : ocean_grid_type use MOM_opacity, only : absorbRemainingSW, optics_type, extract_optics_slice use MOM_unit_scaling, only : unit_scale_type @@ -235,6 +235,9 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ksort ! The sorted k-index that each original layer goes to. real, dimension(SZI_(G),SZJ_(G)) :: & h_miss ! The summed absolute mismatch [H ~> m or kg m-2]. + real, dimension(SZI_(G),SZJ_(G)) :: & + U_star_2d ! The wind friction velocity, calculated using the Boussinesq reference density or + ! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1] real, dimension(SZI_(G)) :: & TKE, & ! The turbulent kinetic energy available for mixing over a ! time step [Z L2 T-2 ~> m3 s-2]. @@ -412,6 +415,9 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C max_BL_det(:) = -1 EOSdom(:) = EOS_domain(G%HI) + ! Extract the friction velocity from the forcing type. + call find_ustar(fluxes, tv, U_star_2d, G, GV, US) + !$OMP parallel default(shared) firstprivate(dKE_CA,cTKE,h_CA,max_BL_det,p_ref,p_ref_cv) & !$OMP private(h,u,v,h_orig,eps,T,S,opacity_band,d_ea,d_eb,R0,Rcv,ksort, & !$OMP dR0_dT,dR0_dS,dRcv_dT,dRcv_dS,htot,Ttot,Stot,TKE,Conv_en, & @@ -513,7 +519,7 @@ subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, US, C ! First the TKE at the depth of free convection that is available ! to drive mixing is calculated. - call find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, & + call find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_FC, dKE_CA, & TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, & j, ksort, G, GV, US, CS) @@ -1252,7 +1258,7 @@ end subroutine mixedlayer_convection !> This subroutine determines the TKE available at the depth of free !! convection to drive mechanical entrainment. -subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, & +subroutine find_starting_TKE(htot, h_CA, fluxes, U_star_2d, Conv_En, cTKE, dKE_FC, dKE_CA, & TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, & j, ksort, G, GV, US, CS) type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. @@ -1265,6 +1271,10 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, type(forcing), intent(in) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields !! have NULL pointers. + real, dimension(SZI_(G),SZJ_(G)), intent(in) :: U_star_2d !< The wind friction velocity, calculated + !! using the Boussinesq reference density or + !! the time-evolving surface density in + !! non-Boussinesq mode [Z T-1 ~> m s-1] real, dimension(SZI_(G)), intent(inout) :: Conv_En !< The buoyant turbulent kinetic energy source !! due to free convection [Z L2 T-2 ~> m3 s-2]. real, dimension(SZI_(G)), intent(in) :: dKE_FC !< The vertically integrated change in @@ -1325,7 +1335,8 @@ subroutine find_starting_TKE(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, if (CS%omega_frac >= 1.0) absf = 2.0*CS%omega do i=is,ie - U_star = fluxes%ustar(i,j) + U_star = U_star_2d(i,j) + if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then if (fluxes%frac_shelf_h(i,j) > 0.0) & U_star = (1.0 - fluxes%frac_shelf_h(i,j)) * U_star + &