diff --git a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 index b2d92c00e7..967667d786 100644 --- a/config_src/drivers/solo_driver/MOM_surface_forcing.F90 +++ b/config_src/drivers/solo_driver/MOM_surface_forcing.F90 @@ -45,9 +45,8 @@ module MOM_surface_forcing use user_surface_forcing, only : USER_surface_forcing_init, user_surface_forcing_CS use user_revise_forcing, only : user_alter_forcing, user_revise_forcing_init use user_revise_forcing, only : user_revise_forcing_CS -use idealized_hurricane, only : idealized_hurricane_wind_init -use idealized_hurricane, only : idealized_hurricane_wind_forcing, SCM_idealized_hurricane_wind_forcing -use idealized_hurricane, only : idealized_hurricane_CS +use idealized_hurricane, only : idealized_hurricane_wind_forcing +use idealized_hurricane, only : idealized_hurricane_wind_init, idealized_hurricane_CS use SCM_CVmix_tests, only : SCM_CVmix_tests_surface_forcing_init use SCM_CVmix_tests, only : SCM_CVmix_tests_wind_forcing use SCM_CVmix_tests, only : SCM_CVmix_tests_buoyancy_forcing @@ -297,7 +296,8 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US elseif (trim(CS%wind_config) == "ideal_hurr") then call idealized_hurricane_wind_forcing(sfc_state, forces, day_center, G, US, CS%idealized_hurricane_CSp) elseif (trim(CS%wind_config) == "SCM_ideal_hurr") then - call SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day_center, G, US, CS%idealized_hurricane_CSp) + call MOM_error(FATAL, "MOM_surface_forcing (set_forcing): "//& + 'WIND_CONFIG = "SCM_ideal_hurr" is a depricated option.') elseif (trim(CS%wind_config) == "SCM_CVmix_tests") then call SCM_CVmix_tests_wind_forcing(sfc_state, forces, day_center, G, US, CS%SCM_CVmix_tests_CSp) elseif (trim(CS%wind_config) == "USER") then @@ -1780,8 +1780,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C call get_param(param_file, mdl, "WIND_CONFIG", CS%wind_config, & "The character string that indicates how wind forcing is specified. Valid "//& "options include (file), (data_override), (2gyre), (1gyre), (gyres), (zero), "//& - "(const), (Neverworld), (scurves), (ideal_hurr), (SCM_ideal_hurr), "//& - "(SCM_CVmix_tests) and (USER).", default="zero") + "(const), (Neverworld), (scurves), (ideal_hurr), (SCM_CVmix_tests) and (USER).", & + default="zero") if (trim(CS%wind_config) == "file") then call get_param(param_file, mdl, "WIND_FILE", CS%wind_file, & "The file in which the wind stresses are found in "//& @@ -1990,9 +1990,13 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C call dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS%dumbbell_forcing_CSp) elseif (trim(CS%wind_config) == "MESO" .or. trim(CS%buoy_config) == "MESO" ) then call MESO_surface_forcing_init(Time, G, US, param_file, diag, CS%MESO_forcing_CSp) - elseif (trim(CS%wind_config) == "ideal_hurr" .or.& - trim(CS%wind_config) == "SCM_ideal_hurr") then + elseif (trim(CS%wind_config) == "ideal_hurr") then call idealized_hurricane_wind_init(Time, G, US, param_file, CS%idealized_hurricane_CSp) + elseif (trim(CS%wind_config) == "SCM_ideal_hurr") then + call MOM_error(FATAL, "MOM_surface_forcing (surface_forcing_init): "//& + 'WIND_CONFIG = "SCM_ideal_hurr" is a depricated option. '//& + 'To obtain mathematically equivalent results set '//& + 'WIND_CONFIG = "ideal_hurr", IDL_HURR_SCM = True and IDL_HURR_X0 = 6.48e+05.') elseif (trim(CS%wind_config) == "const") then call get_param(param_file, mdl, "CONST_WIND_TAUX", CS%tau_x0, & "With wind_config const, this is the constant zonal wind-stress", & diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index e984d5831c..819b06ee50 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/src/parameterizations/vertical/MOM_diabatic_driver.F90 @@ -858,7 +858,7 @@ subroutine diabatic_ALE_legacy(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Tim endif call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) - call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, visc, dt, Kd_ePBL, G, GV, US, & CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) call energetic_PBL_get_MLD(CS%ePBL, BLD(:,:), G, US) @@ -1410,7 +1410,7 @@ subroutine diabatic_ALE(u, v, h, tv, BLD, fluxes, visc, ADp, CDp, dt, Time_end, endif call find_uv_at_h(u, v, h, u_h, v_h, G, GV, US) - call energetic_PBL(h, u_h, v_h, tv, fluxes, dt, Kd_ePBL, G, GV, US, & + call energetic_PBL(h, u_h, v_h, tv, fluxes, visc, dt, Kd_ePBL, G, GV, US, & CS%ePBL, stoch_CS, dSV_dT, dSV_dS, cTKE, SkinBuoyFlux, waves=waves) call energetic_PBL_get_MLD(CS%ePBL, BLD(:,:), G, US) diff --git a/src/parameterizations/vertical/MOM_energetic_PBL.F90 b/src/parameterizations/vertical/MOM_energetic_PBL.F90 index f10e2f445d..d60e08670c 100644 --- a/src/parameterizations/vertical/MOM_energetic_PBL.F90 +++ b/src/parameterizations/vertical/MOM_energetic_PBL.F90 @@ -5,6 +5,7 @@ module MOM_energetic_PBL use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE use MOM_coms, only : EFP_type, real_to_EFP, EFP_to_real, operator(+), assignment(=), EFP_sum_across_PEs +use MOM_debugging, only : hchksum use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_alloc use MOM_diag_mediator, only : time_type, diag_ctrl use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type @@ -16,7 +17,7 @@ module MOM_energetic_PBL use MOM_intrinsic_functions, only : cuberoot use MOM_string_functions, only : uppercase use MOM_unit_scaling, only : unit_scale_type -use MOM_variables, only : thermo_var_ptrs +use MOM_variables, only : thermo_var_ptrs, vertvisc_type use MOM_verticalGrid, only : verticalGrid_type use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number use MOM_stochastics, only : stochastic_CS @@ -59,6 +60,8 @@ module MOM_energetic_PBL !! self-consistent mixed layer depth. Otherwise use the false position !! after a maximum and minimum bound have been evaluated and the !! returned value from the previous guess or bisection before this. + logical :: MLD_iter_bug !< If true use buggy logic that gives the wrong bounds for the next + !! iteration when successive guesses increase by exactly EPBL_MLD_TOLERANCE. integer :: max_MLD_its !< The maximum number of iterations that can be used to find a !! self-consistent mixed layer depth with Use_MLD_iteration. real :: MixLenExponent !< Exponent in the mixing length shape-function [nondim]. @@ -67,6 +70,10 @@ module MOM_energetic_PBL real :: MKE_to_TKE_effic !< The efficiency with which mean kinetic energy released by !! mechanically forced entrainment of the mixed layer is converted to !! TKE [nondim]. + logical :: direct_calc !< If true and there is no conversion from mean kinetic energy to ePBL + !! turbulent kinetic energy, use a direct calculation of the + !! diffusivity that is supported by a given energy input instead of the + !! more general but slower iterative solver. real :: ustar_min !< A minimum value of ustar to avoid numerical problems [Z T-1 ~> m s-1]. !! If the value is small enough, this should not affect the solution. real :: Ekman_scale_coef !< A nondimensional scaling factor controlling the inhibition of the @@ -92,9 +99,8 @@ module MOM_energetic_PBL !! Making this larger increases the diffusivity. real :: vstar_surf_fac !< If (wT_scheme == wT_from_RH18) this is the proportionality coefficient between !! ustar and the surface mechanical contribution to vstar [nondim] - real :: vstar_scale_fac !< An overall nondimensional scaling factor for vstar times a unit - !! conversion factor [Z s T-1 m-1 ~> nondim]. Making this larger increases - !! the diffusivity. + real :: vstar_scale_fac !< An overall nondimensional scaling factor for vstar [nondim]. Making + !! this larger increases the diffusivity. !mstar related options integer :: mstar_scheme !< An encoded integer to determine which formula is used to set mstar @@ -155,6 +161,43 @@ module MOM_energetic_PBL !! the Ekman depth over the Obukhov depth with destabilizing forcing [nondim]. real :: Max_Enhance_M = 5. !< The maximum allowed LT enhancement to the mixing [nondim]. + !/ Bottom boundary layer mixing related options + real :: ePBL_BBL_effic !< The efficiency of bottom boundary layer mixing via ePBL [nondim] + logical :: Use_BBLD_iteration !< If true, use the proximity to the top of the actively turbulent + !! bottom boundary layer to constrain the mixing lengths. + real :: TKE_decay_BBL !< The ratio of the natural Ekman depth to the TKE decay scale for + !! bottom boundary layer mixing [nondim] + real :: min_BBL_mix_len !< The minimum mixing length scale that will be used by ePBL in the bottom + !! boundary layer mixing [Z ~> m]. The default (0) does not set a minimum. + real :: MixLenExponent_BBL !< Exponent in the bottom boundary layer mixing length shape-function [nondim]. + !! 1 is law-of-the-wall at top and bottom, + !! 2 is more KPP like. + real :: BBLD_tol !< The tolerance for the iteratively determined bottom boundary layer depth [Z ~> m]. + !! This is only used with USE_MLD_ITERATION. + integer :: max_BBLD_its !< The maximum number of iterations that can be used to find a self-consistent + !! bottom boundary layer depth. + integer :: wT_scheme_BBL !< An enumerated value indicating the method for finding the bottom boundary + !! layer turbulent velocity scale. There are currently two options: + !! wT_mwT_from_cRoot_TKE is the original (TKE_remaining)^1/3 + !! wT_from_RH18 is the version described by Reichl and Hallberg, 2018 + real :: vstar_scale_fac_BBL !< An overall nondimensional scaling factor for wT in the bottom boundary layer [nondim]. + !! Making this larger increases the bottom boundary layer diffusivity.", & + real :: vstar_surf_fac_BBL !< If (wT_scheme_BBL == wT_from_RH18) this is the proportionality coefficient between + !! ustar and the bottom boundayer layer mechanical contribution to vstar [nondim] + real :: Ekman_scale_coef_BBL !< A nondimensional scaling factor controlling the inhibition of the + !! diffusive length scale by rotation in the bottom boundary layer [nondim]. + !! Making this larger decreases the bottom boundary layer diffusivity. + logical :: decay_adjusted_BBL_TKE !< If true, include an adjustment factor in the bottom boundary layer + !! energetics that accounts for an exponential decay of TKE from a + !! near-bottom source and an assumed piecewise linear linear profile + !! of the buoyancy flux response to a change in a diffusivity. + + !/ Options for documenting differences from parameter choices + integer :: options_diff !< If positive, this is a coded integer indicating a pair of + !! settings whose differences are diagnosed in a passive diagnostic mode + !! via extra calls to ePBL_column. If this is 0 or negative no extra + !! calls occur. + !/ Others type(time_type), pointer :: Time=>NULL() !< A pointer to the ocean model's clock. @@ -170,38 +213,30 @@ module MOM_energetic_PBL !! potential energy change code. Otherwise, it uses a newer version !! that can work with successive increments to the diffusivity in !! upward or downward passes. + logical :: debug !< If true, write verbose checksums for debugging purposes. type(diag_ctrl), pointer :: diag=>NULL() !< A structure that is used to regulate the !! timing of diagnostic output. real, allocatable, dimension(:,:) :: & - ML_depth !< The mixed layer depth determined by active mixing in ePBL [H ~> m or kg m-2] - ! These are terms in the mixed layer TKE budget, all in [R Z3 T-3 ~> W m-2 = kg s-3]. + ML_depth !< The mixed layer depth determined by active mixing in ePBL, which may + !! be used for the first guess in the next time step [H ~> m or kg m-2] real, allocatable, dimension(:,:) :: & - diag_TKE_wind, & !< The wind source of TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_MKE, & !< The resolved KE source of TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_conv, & !< The convective source of TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_forcing, & !< The TKE sink required to mix surface penetrating shortwave heating - !! [R Z3 T-3 ~> W m-2]. - diag_TKE_mech_decay, & !< The decay of mechanical TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_conv_decay, & !< The decay of convective TKE [R Z3 T-3 ~> W m-2]. - diag_TKE_mixing, & !< The work done by TKE to deepen the mixed layer [R Z3 T-3 ~> W m-2]. - ! These additional diagnostics are also 2d. - MSTAR_MIX, & !< Mstar used in EPBL [nondim] - MSTAR_LT, & !< Mstar due to Langmuir turbulence [nondim] - LA, & !< Langmuir number [nondim] - LA_MOD !< Modified Langmuir number [nondim] + BBL_depth !< The bottom boundary layer depth determined by active mixing in ePBL [H ~> m or kg m-2] type(EFP_type), dimension(2) :: sum_its !< The total number of iterations and columns worked on + type(EFP_type), dimension(2) :: sum_its_BBL !< The total number of iterations and columns worked on - real, allocatable, dimension(:,:,:) :: & - Velocity_Scale, & !< The velocity scale used in getting Kd [Z T-1 ~> m s-1] - Mixing_Length !< The length scale used in getting Kd [Z ~> m] !>@{ Diagnostic IDs integer :: id_ML_depth = -1, id_hML_depth = -1, id_TKE_wind = -1, id_TKE_mixing = -1 integer :: id_TKE_MKE = -1, id_TKE_conv = -1, id_TKE_forcing = -1 integer :: id_TKE_mech_decay = -1, id_TKE_conv_decay = -1 integer :: id_Mixing_Length = -1, id_Velocity_Scale = -1 + integer :: id_Kd_BBL = -1, id_BBL_Mix_Length = -1, id_BBL_Vel_Scale = -1 + integer :: id_TKE_BBL = -1, id_TKE_BBL_mixing = -1, id_TKE_BBL_decay = -1 + integer :: id_ustar_BBL = -1, id_BBL_decay_scale = -1, id_BBL_depth = -1 integer :: id_MSTAR_mix = -1, id_LA_mod = -1, id_LA = -1, id_MSTAR_LT = -1 + ! The next options are used when passively diagnosing sensitivities from parameter choices + integer :: id_opt_diff_Kd_ePBL = -1, id_opt_maxdiff_Kd_ePBL = -1, id_opt_diff_hML_depth = -1 !>@} end type energetic_PBL_CS @@ -236,11 +271,14 @@ module MOM_energetic_PBL !>@{ Local column copies of energy change diagnostics, all in [R Z3 T-3 ~> W m-2]. real :: dTKE_conv, dTKE_forcing, dTKE_wind, dTKE_mixing ! Local column diagnostics [R Z3 T-3 ~> W m-2] real :: dTKE_MKE, dTKE_mech_decay, dTKE_conv_decay ! Local column diagnostics [R Z3 T-3 ~> W m-2] + real :: dTKE_BBL, dTKE_BBL_decay, dTKE_BBL_mixing ! Local column diagnostics [R Z3 T-3 ~> W m-2] !>@} real :: LA !< The value of the Langmuir number [nondim] real :: LAmod !< The modified Langmuir number by convection [nondim] real :: mstar !< The value of mstar used in ePBL [nondim] real :: mstar_LT !< The portion of mstar due to Langmuir turbulence [nondim] + integer :: OBL_its !< The number of iterations used to find a self-consistent surface boundary layer depth + integer :: BBL_its !< The number of iterations used to find a self-consistent bottom boundary layer depth end type ePBL_column_diags contains @@ -249,7 +287,7 @@ module MOM_energetic_PBL !! mixed layer model. It assumes that heating, cooling and freshwater fluxes !! have already been applied. All calculations are done implicitly, and there !! is no stability limit on the time step. -subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS, & +subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, visc, dt, Kd_int, G, GV, US, CS, & stoch_CS, dSV_dT, dSV_dS, TKE_forced, buoy_flux, Waves ) type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. @@ -279,6 +317,8 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS type(forcing), intent(inout) :: fluxes !< A structure containing pointers to any !! possible forcing fields. Unused fields have !! NULL ptrs. + type(vertvisc_type), intent(in) :: visc !< Structure with vertical viscosities, + !! BBL properties and related fields real, intent(in) :: dt !< Time increment [T ~> s]. real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), & intent(out) :: Kd_int !< The diagnosed diffusivities at interfaces @@ -287,7 +327,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS real, dimension(SZI_(G),SZJ_(G)), & intent(in) :: buoy_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3]. type(wave_parameters_CS), pointer :: Waves !< Waves control structure for Langmuir turbulence - type(stochastic_CS), pointer :: stoch_CS !< The control structure returned by a previous + type(stochastic_CS), pointer :: stoch_CS !< The control structure returned by a previous ! This subroutine determines the diffusivities from the integrated energetics ! mixed layer model. It assumes that heating, cooling and freshwater fluxes @@ -335,12 +375,17 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS u, & ! The zonal velocity [L T-1 ~> m s-1]. v ! The meridional velocity [L T-1 ~> m s-1]. real, dimension(SZK_(GV)+1) :: & - Kd, & ! The diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + Kd, & ! The diapycnal diffusivity due to ePBL [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. mixvel, & ! A turbulent mixing velocity [Z T-1 ~> m s-1]. mixlen, & ! A turbulent mixing length [Z ~> m]. - SpV_dt ! Specific volume interpolated to interfaces divided by dt or 1.0 / (dt * Rho0) + mixvel_BBL, & ! A bottom boundary layer turbulent mixing velocity [Z T-1 ~> m s-1]. + mixlen_BBL, & ! A bottom boundary layer turbulent mixing length [Z ~> m]. + Kd_BBL, & ! The bottom boundary layer diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + SpV_dt, & ! Specific volume interpolated to interfaces divided by dt or 1.0 / (dt * Rho0), + ! in [R-1 T-1 ~> m3 kg-1 s-1], used to convert local TKE into a turbulence velocity cubed. + SpV_dt_cf ! Specific volume interpolated to interfaces divided by dt or 1.0 / (dt * Rho0) ! times conversion factors for answer dates before 20240101 in - ! [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1] or without the convsersion factors for + ! [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1] or without the conversion factors for ! answer dates of 20240101 and later in [R-1 T-1 ~> m3 kg-1 s-1], used to ! convert local TKE into a turbulence velocity cubed. real :: h_neglect ! A thickness that is so small it is usually lost @@ -351,6 +396,9 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS real :: U_Star_Mean ! The surface friction without gustiness [Z T-1 ~> m s-1]. real :: mech_TKE ! The mechanically generated turbulent kinetic energy available for mixing over a ! timestep before the application of the efficiency in mstar [R Z3 T-2 ~> J m-2] + real :: u_star_BBL ! The bottom boundary layer friction velocity [H T-1 ~> m s-1 or kg m-2 s-1]. + real :: BBL_TKE ! The mechanically generated turbulent kinetic energy available for bottom + ! boundary layer mixing within a timestep [R Z3 T-2 ~> J m-2] real :: I_rho ! The inverse of the Boussinesq reference density times a ratio of scaling ! factors [Z L-1 R-1 ~> m3 kg-1] real :: I_dt ! The Adcroft reciprocal of the timestep [T-1 ~> s-1] @@ -358,9 +406,63 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! step [R-1 T-1 ~> m3 kg-1 s-1] real :: B_Flux ! The surface buoyancy flux [Z2 T-3 ~> m2 s-3] real :: MLD_io ! The mixed layer depth found by ePBL_column [Z ~> m] + real :: BBLD_io ! The bottom boundary layer thickness found by ePBL_BBL_column [Z ~> m] + real :: MLD_in ! The first guess at the mixed layer depth [Z ~> m] + real :: BBLD_in ! The first guess at the bottom boundary layer thickness [Z ~> m] type(ePBL_column_diags) :: eCD ! A container for passing around diagnostics. + ! The following variables are used for diagnostics + real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: & + diag_Velocity_Scale, & ! The velocity scale used in getting Kd [Z T-1 ~> m s-1] + diag_Mixing_Length, & ! The length scale used in getting Kd [Z ~> m] + Kd_BBL_3d, & ! The bottom boundary layer diffusivities [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + BBL_Vel_Scale, & ! The velocity scale used in getting the BBL part of Kd [Z T-1 ~> m s-1] + BBL_Mix_Length ! The length scale used in getting the BBL part of Kd [Z ~> m] + real, dimension(SZI_(G),SZJ_(G)) :: & + ! The next 7 diagnostics are terms in the mixed layer TKE budget, all in [R Z3 T-3 ~> W m-2 = kg s-3]. + diag_TKE_wind, & ! The wind source of TKE [R Z3 T-3 ~> W m-2] + diag_TKE_MKE, & ! The resolved KE source of TKE [R Z3 T-3 ~> W m-2] + diag_TKE_conv, & ! The convective source of TKE [R Z3 T-3 ~> W m-2] + diag_TKE_forcing, & ! The TKE sink required to mix surface penetrating shortwave heating [R Z3 T-3 ~> W m-2] + diag_TKE_mech_decay, & ! The decay of mechanical TKE [R Z3 T-3 ~> W m-2] + diag_TKE_conv_decay, & ! The decay of convective TKE [R Z3 T-3 ~> W m-2] + diag_TKE_mixing, & ! The work done by TKE to deepen the mixed layer [R Z3 T-3 ~> W m-2] + diag_TKE_BBL, & ! The source of TKE to the bottom boundary layer [R Z3 T-3 ~> W m-2]. + diag_TKE_BBL_mixing, & ! The work done by TKE to thicken the bottom boundary layer [R Z3 T-3 ~> W m-2]. + diag_TKE_BBL_decay, & ! The work lost to decy of mechanical TKE in the bottom boundary + ! layer [R Z3 T-3 ~> W m-2]. + diag_ustar_BBL, & ! The bottom boundary layer friction velocity [H T-1 ~> m s-1 or kg m-2 s-1] + diag_BBL_decay_scale, & ! The bottom boundary layer TKE decay length scale [H ~> m] + + diag_mStar_MIX, & ! Mstar used in EPBL [nondim] + diag_mStar_LT, & ! Mstar due to Langmuir turbulence [nondim] + diag_LA, & ! Langmuir number [nondim] + diag_LA_MOD ! Modified Langmuir number [nondim] + + ! The following variables are only used for diagnosing sensitivities to ePBL settings + real, dimension(SZK_(GV)+1) :: & + Kd_1, Kd_2 ! Diapycnal diffusivities found with different ePBL options [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real :: diff_Kd(SZI_(G),SZJ_(G),SZK_(GV)+1) ! The change in diapycnal diffusivities found with different + ! ePBL options [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real :: max_abs_diff_Kd(SZI_(G),SZJ_(G)) ! The column maximum magnitude of the change in diapycnal + ! diffusivities found with different ePBL options [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real :: diff_hML_depth(SZI_(G),SZJ_(G)) ! The change in diagnosed active mixing layer depth with + ! different ePBL options [Z ~> m] + real :: BLD_1, BLD_2 ! Surface or bottom boundary layer depths found with different ePBL_column options [Z ~> m] + real :: SpV_scale1 ! A factor that accounts for the varying scaling of SpV_dt with answer date + ! [nondim] or [T3 m3 Z-3 s-3 ~> 1] + real :: SpV_scale2 ! A factor that accounts for the varying scaling of SpV_dt with answer date + ! [nondim] or [Z3 s3 T-3 m-3 ~> 1] + real :: SpV_dt_tmp(SZK_(GV)+1) ! Specific volume interpolated to interfaces divided by dt or 1.0 / (dt * Rho0) + ! times conversion factors for answer dates before 20240101 in + ! [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1] or without the conversion factors for + ! answer dates of 20240101 and later in [R-1 T-1 ~> m3 kg-1 s-1], used to + ! convert local TKE into a turbulence velocity cubed. + type(ePBL_column_diags) :: eCD_tmp ! A container for not passing around diagnostics. + type(energetic_PBL_CS) :: CS_tmp1, CS_tmp2 ! Copies of the energetic PBL control structure that + ! can be modified to test for sensitivities + integer :: i, j, k, is, ie, js, je, nz is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke @@ -386,16 +488,60 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! Zero out diagnostics before accumulation. if (CS%TKE_diagnostics) then - !!OMP parallel do default(none) shared(is,ie,js,je,CS) + !!OMP parallel do default(shared) do j=js,je ; do i=is,ie - CS%diag_TKE_wind(i,j) = 0.0 ; CS%diag_TKE_MKE(i,j) = 0.0 - CS%diag_TKE_conv(i,j) = 0.0 ; CS%diag_TKE_forcing(i,j) = 0.0 - CS%diag_TKE_mixing(i,j) = 0.0 ; CS%diag_TKE_mech_decay(i,j) = 0.0 - CS%diag_TKE_conv_decay(i,j) = 0.0 !; CS%diag_TKE_unbalanced(i,j) = 0.0 + diag_TKE_wind(i,j) = 0.0 ; diag_TKE_MKE(i,j) = 0.0 + diag_TKE_conv(i,j) = 0.0 ; diag_TKE_forcing(i,j) = 0.0 + diag_TKE_mixing(i,j) = 0.0 ; diag_TKE_mech_decay(i,j) = 0.0 + diag_TKE_conv_decay(i,j) = 0.0 !; diag_TKE_unbalanced(i,j) = 0.0 enddo ; enddo + if (CS%ePBL_BBL_effic > 0.0) then + !!OMP parallel do default(shared) + do j=js,je ; do i=is,ie + diag_TKE_BBL(i,j) = 0.0 ; diag_TKE_BBL_mixing(i,j) = 0.0 + diag_TKE_BBL_decay(i,j) = 0.0 + enddo ; enddo + endif + endif + if (CS%debug .or. (CS%id_Mixing_Length>0)) diag_Mixing_Length(:,:,:) = 0.0 + if (CS%debug .or. (CS%id_Velocity_Scale>0)) diag_Velocity_Scale(:,:,:) = 0.0 + if (CS%ePBL_BBL_effic > 0.0) then + if (CS%debug .or. (CS%id_BBL_Mix_Length>0)) BBL_Mix_Length(:,:,:) = 0.0 + if (CS%debug .or. (CS%id_BBL_Vel_Scale>0)) BBL_Vel_Scale(:,:,:) = 0.0 + if (CS%id_Kd_BBL > 0) Kd_BBL_3d(:,:,:) = 0.0 + if (CS%id_ustar_BBL > 0) diag_ustar_BBL(:,:) = 0.0 + if (CS%id_BBL_decay_scale > 0) diag_BBL_decay_scale(:,:) = 0.0 + endif + + ! CS_tmp is used to test sensitivity to parameter setting changes. + if (CS%options_diff > 0) then + CS_tmp1 = CS ; CS_tmp2 = CS + SpV_scale1 = 1.0 ; SpV_scale2 = 1.0 + + if (CS%options_diff == 1) then + CS_tmp1%orig_PE_calc = .true. ; CS_tmp2%orig_PE_calc = .false. + elseif (CS%options_diff == 2) then + CS_tmp1%answer_date = 20181231 ; CS_tmp2%answer_date = 20240101 + elseif (CS%options_diff == 3) then + CS_tmp1%direct_calc = .true. ; CS_tmp2%direct_calc = .false. + CS_tmp1%MKE_to_TKE_effic = 0.0 ; CS_tmp2%MKE_to_TKE_effic = 0.0 + CS_tmp1%orig_PE_calc = .false. ; CS_tmp2%orig_PE_calc = .false. + elseif (CS%options_diff == 4) then + CS_tmp1%direct_calc = .true. ; CS_tmp2%direct_calc = .false. + CS_tmp1%MKE_to_TKE_effic = 0.0 ; CS_tmp2%MKE_to_TKE_effic = 0.0 + CS_tmp1%ePBL_BBL_effic = 0.2 ; CS_tmp2%ePBL_BBL_effic = 0.2 + elseif (CS%options_diff == 5) then + CS_tmp1%decay_adjusted_BBL_TKE = .true. ; CS_tmp2%decay_adjusted_BBL_TKE = .false. + CS_tmp1%MKE_to_TKE_effic = 0.0 ; CS_tmp2%MKE_to_TKE_effic = 0.0 + CS_tmp1%ePBL_BBL_effic = 0.2 ; CS_tmp2%ePBL_BBL_effic = 0.2 + endif + ! This logic is needed because the scaling of SpV_dt changes with answer date. + if (CS_tmp1%answer_date < 20240101) SpV_scale1 = US%m_to_Z**3 * US%T_to_s**3 + if (CS_tmp2%answer_date < 20240101) SpV_scale2 = US%m_to_Z**3 * US%T_to_s**3 + if (CS%id_opt_diff_Kd_ePBL > 0) diff_Kd(:,:,:) = 0.0 + if (CS%id_opt_maxdiff_Kd_ePBL > 0) max_abs_diff_Kd(:,:) = 0.0 + if (CS%id_opt_diff_hML_depth > 0) diff_hML_depth(:,:) = 0.0 endif - ! if (CS%id_Mixing_Length>0) CS%Mixing_Length(:,:,:) = 0.0 - ! if (CS%id_Velocity_Scale>0) CS%Velocity_Scale(:,:,:) = 0.0 !!OMP parallel do default(private) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt,I_dt, & !!OMP CS,G,GV,US,fluxes,TKE_forced,dSV_dT,dSV_dS,Kd_int) @@ -414,7 +560,7 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS if ((dt > 0.0) .and. GV%Boussinesq .or. .not.allocated(tv%SpV_avg)) then if (CS%answer_date < 20240101) then do K=1,nz+1 - SpV_dt(K) = (US%Z_to_m**3*US%s_to_T**3) / (dt*GV%Rho0) + SpV_dt(K) = 1.0 / (dt*GV%Rho0) enddo else do K=1,nz+1 @@ -457,19 +603,11 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS endif if (allocated(tv%SpV_avg) .and. .not.GV%Boussinesq) then - if (CS%answer_date < 20240101) then - SpV_dt(1) = (US%Z_to_m**3*US%s_to_T**3) * tv%SpV_avg(i,j,1) * I_dt - do K=2,nz - SpV_dt(K) = (US%Z_to_m**3*US%s_to_T**3) * 0.5*(tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k)) * I_dt - enddo - SpV_dt(nz+1) = (US%Z_to_m**3*US%s_to_T**3) * tv%SpV_avg(i,j,nz) * I_dt - else - SpV_dt(1) = tv%SpV_avg(i,j,1) * I_dt - do K=2,nz - SpV_dt(K) = 0.5*(tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k)) * I_dt - enddo - SpV_dt(nz+1) = tv%SpV_avg(i,j,nz) * I_dt - endif + SpV_dt(1) = tv%SpV_avg(i,j,1) * I_dt + do K=2,nz + SpV_dt(K) = 0.5*(tv%SpV_avg(i,j,k-1) + tv%SpV_avg(i,j,k)) * I_dt + enddo + SpV_dt(nz+1) = tv%SpV_avg(i,j,nz) * I_dt endif B_flux = buoy_flux(i,j) @@ -491,77 +629,184 @@ subroutine energetic_PBL(h_3d, u_3d, v_3d, tv, fluxes, dt, Kd_int, G, GV, US, CS ! Perhaps provide a first guess for MLD based on a stored previous value. MLD_io = -1.0 if (CS%MLD_iteration_guess .and. (CS%ML_depth(i,j) > 0.0)) MLD_io = CS%ML_depth(i,j) + BBLD_io = 0.0 + + ! Store the initial guesses at the boundary layer depths for testing sensitivities. + MLD_in = MLD_io + if (CS%answer_date < 20240101) then + do K=1,nz+1 ; SpV_dt_cf(K) = (US%Z_to_m**3*US%s_to_T**3) * SpV_dt(K) ; enddo + else + do K=1,nz+1 ; SpV_dt_cf(K) = SpV_dt(K) ; enddo + endif if (stoch_CS%pert_epbl) then ! stochastics are active - call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, TKE_forcing, B_flux, absf, & + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt_cf, TKE_forcing, B_flux, absf, & u_star, u_star_mean, mech_TKE, dt, MLD_io, Kd, mixvel, mixlen, GV, & US, CS, eCD, Waves, G, i, j, & TKE_gen_stoch=stoch_CS%epbl1_wts(i,j), TKE_diss_stoch=stoch_CS%epbl2_wts(i,j)) else - call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, TKE_forcing, B_flux, absf, & + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt_cf, TKE_forcing, B_flux, absf, & u_star, u_star_mean, mech_TKE, dt, MLD_io, Kd, mixvel, mixlen, GV, & US, CS, eCD, Waves, G, i, j) endif + ! Add the diffusivity due to bottom boundary layer mixing, if there is energy to drive this mixing. + if (CS%ePBL_BBL_effic > 0.0) then + if (CS%MLD_iteration_guess .and. (CS%BBL_depth(i,j) > 0.0)) BBLD_io = CS%BBL_depth(i,j) + BBLD_in = BBLD_io + BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%TKE_BBL(i,j) + u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H) + call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & + u_star_BBL, Kd_BBL, BBLD_io, mixvel_BBL, mixlen_BBL, GV, US, CS, eCD) + + do K=1,nz+1 ; Kd(K) = Kd(K) + Kd_BBL(K) ; enddo + if (CS%id_Kd_BBL > 0) then ; do K=1,nz+1 + Kd_BBL_3d(i,j,K) = Kd_BBL(K) + enddo ; endif + if (CS%id_ustar_BBL > 0) diag_ustar_BBL(i,j) = u_star_BBL + if ((CS%id_BBL_decay_scale > 0) .and. (CS%TKE_decay * absf > 0)) & + diag_BBL_decay_scale(i,j) = u_star_BBL / (CS%TKE_decay * absf) + endif + ! Copy the diffusivities to a 2-d array. do K=1,nz+1 Kd_2d(i,K) = Kd(K) enddo CS%ML_depth(i,j) = MLD_io + CS%BBL_depth(i,j) = BBLD_io if (CS%TKE_diagnostics) then - CS%diag_TKE_MKE(i,j) = CS%diag_TKE_MKE(i,j) + eCD%dTKE_MKE - CS%diag_TKE_conv(i,j) = CS%diag_TKE_conv(i,j) + eCD%dTKE_conv - CS%diag_TKE_forcing(i,j) = CS%diag_TKE_forcing(i,j) + eCD%dTKE_forcing - CS%diag_TKE_wind(i,j) = CS%diag_TKE_wind(i,j) + eCD%dTKE_wind - CS%diag_TKE_mixing(i,j) = CS%diag_TKE_mixing(i,j) + eCD%dTKE_mixing - CS%diag_TKE_mech_decay(i,j) = CS%diag_TKE_mech_decay(i,j) + eCD%dTKE_mech_decay - CS%diag_TKE_conv_decay(i,j) = CS%diag_TKE_conv_decay(i,j) + eCD%dTKE_conv_decay - ! CS%diag_TKE_unbalanced(i,j) = CS%diag_TKE_unbalanced(i,j) + eCD%dTKE_unbalanced + diag_TKE_MKE(i,j) = diag_TKE_MKE(i,j) + eCD%dTKE_MKE + diag_TKE_conv(i,j) = diag_TKE_conv(i,j) + eCD%dTKE_conv + diag_TKE_forcing(i,j) = diag_TKE_forcing(i,j) + eCD%dTKE_forcing + diag_TKE_wind(i,j) = diag_TKE_wind(i,j) + eCD%dTKE_wind + diag_TKE_mixing(i,j) = diag_TKE_mixing(i,j) + eCD%dTKE_mixing + diag_TKE_mech_decay(i,j) = diag_TKE_mech_decay(i,j) + eCD%dTKE_mech_decay + diag_TKE_conv_decay(i,j) = diag_TKE_conv_decay(i,j) + eCD%dTKE_conv_decay + ! diag_TKE_unbalanced(i,j) = diag_TKE_unbalanced(i,j) + eCD%dTKE_unbalanced endif - ! Write to 3-D for outputting Mixing length and velocity scale. - if (CS%id_Mixing_Length>0) then ; do k=1,nz - CS%Mixing_Length(i,j,k) = mixlen(k) + ! Write mixing length and velocity scale to 3-D arrays for diagnostic output + if (CS%debug .or. (CS%id_Mixing_Length > 0)) then ; do K=1,nz+1 + diag_Mixing_Length(i,j,K) = mixlen(K) enddo ; endif - if (CS%id_Velocity_Scale>0) then ; do k=1,nz - CS%Velocity_Scale(i,j,k) = mixvel(k) + if (CS%debug .or. (CS%id_Velocity_Scale > 0)) then ; do K=1,nz+1 + diag_Velocity_Scale(i,j,K) = mixvel(K) enddo ; endif - if (allocated(CS%mstar_mix)) CS%mstar_mix(i,j) = eCD%mstar - if (allocated(CS%mstar_lt)) CS%mstar_lt(i,j) = eCD%mstar_LT - if (allocated(CS%La)) CS%La(i,j) = eCD%LA - if (allocated(CS%La_mod)) CS%La_mod(i,j) = eCD%LAmod + if (CS%ePBL_BBL_effic > 0.0) then + if (CS%debug .or. (CS%id_BBL_Mix_Length>0)) then ; do k=1,nz + BBL_Mix_Length(i,j,k) = mixlen_BBL(k) + enddo ; endif + if (CS%debug .or. (CS%id_BBL_Vel_Scale>0)) then ; do k=1,nz + BBL_Vel_Scale(i,j,k) = mixvel_BBL(k) + enddo ; endif + endif + if (CS%id_MSTAR_MIX > 0) diag_mStar_mix(i,j) = eCD%mstar + if (CS%id_MSTAR_LT > 0) diag_mStar_lt(i,j) = eCD%mstar_LT + if (CS%id_LA > 0) diag_LA(i,j) = eCD%LA + if (CS%id_LA_MOD > 0) diag_LA_mod(i,j) = eCD%LAmod + if (report_avg_its) then + CS%sum_its(1) = CS%sum_its(1) + real_to_EFP(real(eCD%OBL_its)) + CS%sum_its(2) = CS%sum_its(2) + real_to_EFP(1.0) + if (CS%ePBL_BBL_effic > 0.0) then + CS%sum_its_BBL(1) = CS%sum_its_BBL(1) + real_to_EFP(real(eCD%BBL_its)) + CS%sum_its_BBL(2) = CS%sum_its_BBL(2) + real_to_EFP(1.0) + endif + endif + + if (CS%options_diff > 0) then + ! Call ePBL_column of ePBL_BBL_column with different parameter settings to diagnose sensitivities. + ! These do not change the model state, and are only used for diagnostic purposes. + if (CS%options_diff < 4) then + BLD_1 = MLD_in ; BLD_2 = MLD_in + do K=1,nz+1 ; SpV_dt_tmp(K) = SpV_scale1 * SpV_dt(K) ; enddo + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt_tmp, TKE_forcing, & + B_flux, absf, u_star, u_star_mean, mech_TKE, dt, BLD_1, Kd_1, & + mixvel, mixlen, GV, US, CS_tmp1, eCD_tmp, Waves, G, i, j) + do K=1,nz+1 ; SpV_dt_tmp(K) = SpV_scale2 * SpV_dt(K) ; enddo + call ePBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt_tmp, TKE_forcing, & + B_flux, absf, u_star, u_star_mean, mech_TKE, dt, BLD_2, Kd_2, & + mixvel, mixlen, GV, US, CS_tmp2, eCD_tmp, Waves, G, i, j) + else + BLD_1 = BBLD_in ; BLD_2 = BBLD_in + BBL_TKE = CS%ePBL_BBL_effic * GV%H_to_RZ * dt * visc%TKE_BBL(i,j) + u_star_BBL = max(visc%ustar_BBL(i,j), CS%ustar_min*GV%Z_to_H) + call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & + u_star_BBL, Kd_1, BLD_1, mixvel_BBL, mixlen_BBL, GV, US, CS_tmp1, eCD_tmp) + call ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT_1d, dSV_dS_1d, SpV_dt, absf, dt, Kd, BBL_TKE, & + u_star_BBL, Kd_2, BLD_2, mixvel_BBL, mixlen_BBL, GV, US, CS_tmp2, eCD_tmp) + endif + + if (CS%id_opt_diff_Kd_ePBL > 0) then + do K=1,nz+1 ; diff_Kd(i,j,K) = Kd_1(K) - Kd_2(K) ; enddo + endif + if (CS%id_opt_maxdiff_Kd_ePBL > 0) then + max_abs_diff_Kd(i,j) = 0.0 + do K=1,nz+1 ; max_abs_diff_Kd(i,j) = max(max_abs_diff_Kd(i,j), abs(Kd_1(K) - Kd_2(K))) ; enddo + endif + if (CS%id_opt_diff_hML_depth > 0) diff_hML_depth(i,j) = BLD_1 - BLD_2 + endif + else ! End of the ocean-point part of the i-loop ! For masked points, Kd_int must still be set (to 0) because it has intent out. do K=1,nz+1 ; Kd_2d(i,K) = 0. ; enddo CS%ML_depth(i,j) = 0.0 - endif ; enddo ! Close of i-loop - Note unusual loop order! + CS%BBL_depth(i,j) = 0.0 + endif ; enddo ! Close of i-loop - Note the unusual loop order, with k-loops inside i-loops. do K=1,nz+1 ; do i=is,ie ; Kd_int(i,j,K) = Kd_2d(i,K) ; enddo ; enddo enddo ! j-loop + if (CS%debug .and. (CS%ePBL_BBL_effic > 0.0)) then + call hchksum(visc%TKE_BBL, "ePBL visc%TKE_BBL", G%HI, unscale=GV%H_to_MKS*US%Z_to_m**2*US%s_to_T**3) + call hchksum(visc%ustar_BBL, "ePBL visc%ustar_BBL", G%HI, unscale=GV%H_to_MKS*US%s_to_T) + call hchksum(Kd_int, "End of ePBL Kd_int", G%HI, unscale=GV%H_to_MKS*US%Z_to_m*US%s_to_T) + call hchksum(diag_Velocity_Scale, "ePBL Velocity_Scale", G%HI, unscale=US%Z_to_m*US%s_to_T) + call hchksum(diag_Mixing_Length, "ePBL Mixing_Length", G%HI, unscale=US%Z_to_m) + call hchksum(BBL_Vel_Scale, "ePBL BBL_Vel_Scale", G%HI, unscale=US%Z_to_m*US%s_to_T) + call hchksum(BBL_Mix_Length, "ePBL BBL_Mix_Length", G%HI, unscale=US%Z_to_m) + endif + if (CS%id_ML_depth > 0) call post_data(CS%id_ML_depth, CS%ML_depth, CS%diag) if (CS%id_hML_depth > 0) call post_data(CS%id_hML_depth, CS%ML_depth, CS%diag) - if (CS%id_TKE_wind > 0) call post_data(CS%id_TKE_wind, CS%diag_TKE_wind, CS%diag) - if (CS%id_TKE_MKE > 0) call post_data(CS%id_TKE_MKE, CS%diag_TKE_MKE, CS%diag) - if (CS%id_TKE_conv > 0) call post_data(CS%id_TKE_conv, CS%diag_TKE_conv, CS%diag) - if (CS%id_TKE_forcing > 0) call post_data(CS%id_TKE_forcing, CS%diag_TKE_forcing, CS%diag) - if (CS%id_TKE_mixing > 0) call post_data(CS%id_TKE_mixing, CS%diag_TKE_mixing, CS%diag) + if (CS%id_TKE_wind > 0) call post_data(CS%id_TKE_wind, diag_TKE_wind, CS%diag) + if (CS%id_TKE_MKE > 0) call post_data(CS%id_TKE_MKE, diag_TKE_MKE, CS%diag) + if (CS%id_TKE_conv > 0) call post_data(CS%id_TKE_conv, diag_TKE_conv, CS%diag) + if (CS%id_TKE_forcing > 0) call post_data(CS%id_TKE_forcing, diag_TKE_forcing, CS%diag) + if (CS%id_TKE_mixing > 0) call post_data(CS%id_TKE_mixing, diag_TKE_mixing, CS%diag) if (CS%id_TKE_mech_decay > 0) & - call post_data(CS%id_TKE_mech_decay, CS%diag_TKE_mech_decay, CS%diag) + call post_data(CS%id_TKE_mech_decay, diag_TKE_mech_decay, CS%diag) if (CS%id_TKE_conv_decay > 0) & - call post_data(CS%id_TKE_conv_decay, CS%diag_TKE_conv_decay, CS%diag) - if (CS%id_Mixing_Length > 0) call post_data(CS%id_Mixing_Length, CS%Mixing_Length, CS%diag) - if (CS%id_Velocity_Scale >0) call post_data(CS%id_Velocity_Scale, CS%Velocity_Scale, CS%diag) - if (CS%id_MSTAR_MIX > 0) call post_data(CS%id_MSTAR_MIX, CS%MSTAR_MIX, CS%diag) - if (CS%id_LA > 0) call post_data(CS%id_LA, CS%LA, CS%diag) - if (CS%id_LA_MOD > 0) call post_data(CS%id_LA_MOD, CS%LA_MOD, CS%diag) - if (CS%id_MSTAR_LT > 0) call post_data(CS%id_MSTAR_LT, CS%MSTAR_LT, CS%diag) + call post_data(CS%id_TKE_conv_decay, diag_TKE_conv_decay, CS%diag) + if (CS%id_Mixing_Length > 0) call post_data(CS%id_Mixing_Length, diag_Mixing_Length, CS%diag) + if (CS%id_Velocity_Scale >0) call post_data(CS%id_Velocity_Scale, diag_Velocity_Scale, CS%diag) + if (CS%id_MSTAR_MIX > 0) call post_data(CS%id_MSTAR_MIX, diag_mStar_MIX, CS%diag) + if (CS%ePBL_BBL_effic > 0.0) then + if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, Kd_BBL_3d, CS%diag) + if (CS%id_BBL_Mix_Length > 0) call post_data(CS%id_BBL_Mix_Length, BBL_Mix_Length, CS%diag) + if (CS%id_BBL_Vel_Scale > 0) call post_data(CS%id_BBL_Vel_Scale, BBL_Vel_Scale, CS%diag) + if (CS%id_ustar_BBL > 0) call post_data(CS%id_ustar_BBL, diag_ustar_BBL, CS%diag) + if (CS%id_BBL_decay_scale > 0) call post_data(CS%id_BBL_decay_scale, diag_BBL_decay_scale, CS%diag) + if (CS%id_TKE_BBL > 0) call post_data(CS%id_TKE_BBL, diag_TKE_BBL, CS%diag) + if (CS%id_TKE_BBL_mixing > 0) call post_data(CS%id_TKE_BBL_mixing, diag_TKE_BBL_mixing, CS%diag) + if (CS%id_TKE_BBL_decay > 0) call post_data(CS%id_TKE_BBL_decay, diag_TKE_BBL_decay, CS%diag) + if (CS%id_BBL_depth > 0) call post_data(CS%id_BBL_depth, CS%BBL_depth, CS%diag) + endif + if (CS%id_LA > 0) call post_data(CS%id_LA, diag_LA, CS%diag) + if (CS%id_LA_MOD > 0) call post_data(CS%id_LA_MOD, diag_LA_MOD, CS%diag) + if (CS%id_MSTAR_LT > 0) call post_data(CS%id_MSTAR_LT, diag_mStar_LT, CS%diag) if (stoch_CS%pert_epbl) then if (stoch_CS%id_epbl1_wts > 0) call post_data(stoch_CS%id_epbl1_wts, stoch_CS%epbl1_wts, CS%diag) if (stoch_CS%id_epbl2_wts > 0) call post_data(stoch_CS%id_epbl2_wts, stoch_CS%epbl2_wts, CS%diag) endif + if (CS%options_diff > 0) then + ! These diagnostics are only for determining sensitivities to different ePBL settings. + if (CS%id_opt_diff_Kd_ePBL > 0) call post_data(CS%id_opt_diff_Kd_ePBL, diff_Kd, CS%diag) + if (CS%id_opt_maxdiff_Kd_ePBL > 0) call post_data(CS%id_opt_maxdiff_Kd_ePBL, max_abs_diff_Kd, CS%diag) + if (CS%id_opt_diff_hML_depth > 0) call post_data(CS%id_opt_diff_hML_depth, diff_hML_depth, CS%diag) + endif + end subroutine energetic_PBL @@ -591,7 +836,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, !! divided by dt or 1.0 / (dt * Rho0), times conversion !! factors for answer dates before 20240101 in !! [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1] or without - !! the convsersion factors for answer dates of + !! the conversion factors for answer dates of !! 20240101 and later in [R-1 T-1 ~> m3 kg-1 s-1], !! used to convert local TKE into a turbulence !! velocity cubed. @@ -618,14 +863,14 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, !! [Z T-1 ~> m s-1]. real, dimension(SZK_(GV)+1), & intent(out) :: mixlen !< The mixing length scale used in Kd [Z ~> m]. - type(energetic_PBL_CS), intent(inout) :: CS !< Energetic PBL control structure + type(energetic_PBL_CS), intent(in) :: CS !< Energetic PBL control structure type(ePBL_column_diags), intent(inout) :: eCD !< A container for passing around diagnostics. type(wave_parameters_CS), pointer :: Waves !< Waves control structure for Langmuir turbulence - type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure. + type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure. + integer, intent(in) :: i !< The i-index to work on (used for Waves) + integer, intent(in) :: j !< The j-index to work on (used for Waves) real, optional, intent(in) :: TKE_gen_stoch !< random factor used to perturb TKE generation [nondim] real, optional, intent(in) :: TKE_diss_stoch !< random factor used to perturb TKE dissipation [nondim] - integer, intent(in) :: i !< The i-index to work on (used for Waves) - integer, intent(in) :: j !< The i-index to work on (used for Waves) ! This subroutine determines the diffusivities in a single column from the integrated energetics ! planetary boundary layer (ePBL) model. It assumes that heating, cooling and freshwater fluxes @@ -684,6 +929,9 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, Se, & ! Estimated final values of S in the column [S ~> ppt]. dTe, & ! Running (1-way) estimates of temperature change [C ~> degC]. dSe, & ! Running (1-way) estimates of salinity change [S ~> ppt]. + hp_a, & ! An effective pivot thickness of the layer including the effects + ! of coupling with layers above [H ~> m or kg m-2]. This is the first term + ! in the denominator of b1 in a downward-oriented tridiagonal solver. Th_a, & ! An effective temperature times a thickness in the layer above, including implicit ! mixing effects with other yet higher layers [C H ~> degC m or degC kg m-2]. Sh_a, & ! An effective salinity times a thickness in the layer above, including implicit @@ -698,12 +946,9 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, ! the boundary layer [nondim]. h_dz_int, & ! The ratio of the layer thicknesses over the vertical distances ! across the layers surrounding an interface [H Z-1 ~> nondim or kg m-3] - Kddt_h ! The diapycnal diffusivity times a timestep divided by the - ! average thicknesses around a layer [H ~> m or kg m-2]. + Kddt_h ! The total diapycnal diffusivity at an interface times a timestep divided by the + ! average thicknesses around an interface [H ~> m or kg m-2]. real :: b1 ! b1 is inverse of the pivot used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. - real :: hp_a ! An effective pivot thickness of the layer including the effects - ! of coupling with layers above [H ~> m or kg m-2]. This is the first term - ! in the denominator of b1 in a downward-oriented tridiagonal solver. real :: h_neglect ! A thickness that is so small it is usually lost ! in roundoff and can be neglected [H ~> m or kg m-2]. real :: dz_neglect ! A vertical distance that is so small it is usually lost @@ -757,8 +1002,8 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, real :: dMKE_src_dK ! The partial derivative of MKE_src with Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. real :: Kd_guess0 ! A first guess of the diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. real :: PE_chg_g0 ! The potential energy change when Kd is Kd_guess0 [R Z3 T-2 ~> J m-2] - real :: Kddt_h_g0 ! The first guess diapycnal diffusivity times a timestep divided - ! by the average thicknesses around a layer [H ~> m or kg m-2]. + real :: Kddt_h_g0 ! The first guess of the change in diapycnal diffusivity times a timestep + ! divided by the average thicknesses around an interface [H ~> m or kg m-2]. real :: PE_chg_max ! The maximum PE change for very large values of Kddt_h(K) [R Z3 T-2 ~> J m-2]. real :: dPEc_dKd_Kd0 ! The partial derivative of PE change with Kddt_h(K) ! for very small values of Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. @@ -806,12 +1051,18 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, ! can improve this. real :: dMLD_min ! The change in diagnosed mixed layer depth when the guess is min_MLD [Z ~> m] real :: dMLD_max ! The change in diagnosed mixed layer depth when the guess is max_MLD [Z ~> m] - logical :: OBL_converged ! Flag for convergence of MLD integer :: OBL_it ! Iteration counter + real :: TKE_used ! The TKE used to support mixing at an interface [R Z3 T-2 ~> J m-2]. + ! real :: Kd_add ! The additional diffusivity at an interface [H Z T-1 ~> m2 s-1 or kg m-1 s-1] + real :: frac_in_BL ! The fraction of the energy required to support dKd_max that is suppiled by + ! max_PE_chg, used here to determine a fractional layer contribution to the + ! boundary layer thickness [nondim] real :: Surface_Scale ! Surface decay scale for vstar [nondim] logical :: calc_Te ! If true calculate the expected final temperature and salinity values. + logical :: no_MKE_conversion ! If true, there is no conversion from MKE to TKE, so a simpler solver can be used. logical :: debug ! This is used as a hard-coded value for debugging. + logical :: convectively_unstable ! If true, there is convective instability at an interface. ! The following arrays are used only for debugging purposes. real :: dPE_debug ! An estimate of the potential energy change [R Z3 T-2 ~> J m-2] @@ -837,6 +1088,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, debug = .false. ! Change this hard-coded value for debugging. calc_Te = (debug .or. (.not.CS%orig_PE_calc)) + no_MKE_conversion = ((CS%direct_calc) .and. (CS%MKE_to_TKE_effic == 0.0)) h_neglect = GV%H_subroundoff dz_neglect = GV%dZ_subroundoff @@ -905,13 +1157,9 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif ! Iterate to determine a converged EPBL depth. - OBL_converged = .false. do OBL_it=1,CS%Max_MLD_Its - if (.not. OBL_converged) then - ! If not using MLD_Iteration flag loop to only execute once. - if (.not.CS%Use_MLD_iteration) OBL_converged = .true. - + !### The old indenting is being retained for now to simplify the review of some pending changes. if (debug) then ; mech_TKE_k(:) = 0.0 ; conv_PErel_k(:) = 0.0 ; endif ! Reset ML_depth @@ -923,7 +1171,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, call get_Langmuir_Number(LA, G, GV, US, abs(MLD_guess), u_star_mean, i, j, dz, Waves, & U_H=u, V_H=v) call find_mstar(CS, US, B_flux, u_star, MLD_guess, absf, & - MStar_total, Langmuir_Number=La, Convect_Langmuir_Number=LAmod,& + mstar_total, Langmuir_Number=La, Convect_Langmuir_Number=LAmod,& mstar_LT=mstar_LT) else call find_mstar(CS, US, B_flux, u_star, MLD_guess, absf, mstar_total) @@ -931,10 +1179,10 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, !/ Apply MStar to get mech_TKE if ((CS%answer_date < 20190101) .and. (CS%mstar_scheme==Use_Fixed_MStar)) then - mech_TKE = (dt*MSTAR_total*GV%Rho0) * u_star**3 + mech_TKE = (dt*mstar_total*GV%Rho0) * u_star**3 else - mech_TKE = MSTAR_total * mech_TKE_in - ! mech_TKE = MSTAR_total * (dt*GV%Rho0* u_star**3) + mech_TKE = mstar_total * mech_TKE_in + ! mech_TKE = mstar_total * (dt*GV%Rho0* u_star**3) endif ! stochastically perturb mech_TKE in the UFS if (present(TKE_gen_stoch)) mech_TKE = mech_TKE*TKE_gen_stoch @@ -996,7 +1244,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif Kd(1) = 0.0 ; Kddt_h(1) = 0.0 - hp_a = h(1) + hp_a(1) = h(1) dT_to_dPE_a(1) = dT_to_dPE(1) ; dT_to_dColHt_a(1) = dT_to_dColHt(1) dS_to_dPE_a(1) = dS_to_dPE(1) ; dS_to_dColHt_a(1) = dS_to_dColHt(1) @@ -1121,14 +1369,14 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, ! tridiagonal solver for the whole column to be completed for debugging ! purposes, and also allows for something akin to convective adjustment ! in unstable interior regions? - b1 = 1.0 / hp_a + b1 = 1.0 / hp_a(k-1) c1(K) = 0.0 if (CS%orig_PE_calc) then dTe(k-1) = b1 * ( dTe_t2 ) dSe(k-1) = b1 * ( dSe_t2 ) endif - hp_a = h(k) + hp_a(k) = h(k) dT_to_dPE_a(k) = dT_to_dPE(k) dS_to_dPE_a(k) = dS_to_dPE(k) dT_to_dColHt_a(k) = dT_to_dColHt(k) @@ -1145,12 +1393,12 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, dS_km1_t2 = (S0(k)-S0(k-1)) else dT_km1_t2 = (T0(k)-T0(k-1)) - & - (Kddt_h(K-1) / hp_a) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) + (Kddt_h(K-1) / hp_a(k-1)) * ((T0(k-2) - T0(k-1)) + dTe(k-2)) dS_km1_t2 = (S0(k)-S0(k-1)) - & - (Kddt_h(K-1) / hp_a) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) + (Kddt_h(K-1) / hp_a(k-1)) * ((S0(k-2) - S0(k-1)) + dSe(k-2)) endif - dTe_term = dTe_t2 + hp_a * (T0(k-1)-T0(k)) - dSe_term = dSe_t2 + hp_a * (S0(k-1)-S0(k)) + dTe_term = dTe_t2 + hp_a(k-1) * (T0(k-1)-T0(k)) + dSe_term = dSe_t2 + hp_a(k-1) * (S0(k-1)-S0(k)) else if (K<=2) then Th_a(k-1) = h(k-1) * T0(k-1) ; Sh_a(k-1) = h(k-1) * S0(k-1) @@ -1205,7 +1453,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif endif hbs_here = min(hb_hs(K), MixLen_shape(K)) - mixlen(K) = MAX(CS%min_mix_len, ((dz_tt*hbs_here)*vstar) / & + mixlen(K) = max(CS%min_mix_len, ((dz_tt*hbs_here)*vstar) / & ((CS%Ekman_scale_coef * absf) * (dz_tt*hbs_here) + vstar)) !Note setting Kd_guess0 to vstar * CS%vonKar * mixlen(K) here will ! change the answers. Therefore, skipping that. @@ -1221,26 +1469,40 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, mixvel(K) = vstar ! Track vstar Kddt_h_g0 = Kd_guess0 * dt_h - if (CS%orig_PE_calc) then - call find_PE_chg_orig(Kddt_h_g0, h(k), hp_a, dTe_term, dSe_term, & + if (no_MKE_conversion) then + ! Without conversion from MKE to TKE, the updated diffusivity can be determined directly. + ! Replace h(k) with hp_b(k) = h(k), and dT_to_dPE with dT_to_dPE_b, etc., for a 2-direction solver. + call find_Kd_from_PE_chg(0.0, Kd_guess0, dt_h, tot_TKE, hp_a(k-1), h(k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), & + pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & + dT_to_dColHt(k), dS_to_dColHt(k), Kd_add=Kd(K), PE_chg=TKE_used, & + dPE_max=PE_chg_max, frac_dKd_max_PE=frac_in_BL) + convectively_unstable = (PE_chg_max < 0.0) + PE_chg_g0 = TKE_used ! This is only used in the convectively unstable limit. + MKE_src = 0.0 + elseif (CS%orig_PE_calc) then + call find_PE_chg_orig(Kddt_h_g0, h(k), hp_a(k-1), dTe_term, dSe_term, & dT_km1_t2, dS_km1_t2, dT_to_dPE(k), dS_to_dPE(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), & pres_Z(K), dT_to_dColHt(k), dS_to_dColHt(k), & dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & PE_chg=PE_chg_g0, dPE_max=PE_chg_max, dPEc_dKd_0=dPEc_dKd_Kd0 ) + convectively_unstable = (PE_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dPEc_dKd_Kd0 < 0.0)) + MKE_src = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) else - call find_PE_chg(0.0, Kddt_h_g0, hp_a, h(k), & + call find_PE_chg(0.0, Kddt_h_g0, hp_a(k-1), h(k), & Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), & pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & dT_to_dColHt(k), dS_to_dColHt(k), & PE_chg=PE_chg_g0, dPE_max=PE_chg_max, dPEc_dKd_0=dPEc_dKd_Kd0 ) + convectively_unstable = (PE_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dPEc_dKd_Kd0 < 0.0)) + MKE_src = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) endif - MKE_src = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) - ! This block checks out different cases to determine Kd at the present interface. - if ((PE_chg_g0 < 0.0) .or. ((vstar == 0.0) .and. (dPEc_dKd_Kd0 < 0.0))) then + if (convectively_unstable) then ! This column is convectively unstable. if (PE_chg_max <= 0.0) then ! Does MKE_src need to be included in the calculation of vstar here? @@ -1280,14 +1542,14 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, mixvel(K) = vstar if (CS%orig_PE_calc) then - call find_PE_chg_orig(Kd(K)*dt_h, h(k), hp_a, dTe_term, dSe_term, & + call find_PE_chg_orig(Kd(K)*dt_h, h(k), hp_a(k-1), dTe_term, dSe_term, & dT_km1_t2, dS_km1_t2, dT_to_dPE(k), dS_to_dPE(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), & pres_Z(K), dT_to_dColHt(k), dS_to_dColHt(k), & dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & PE_chg=dPE_conv) else - call find_PE_chg(0.0, Kd(K)*dt_h, hp_a, h(k), & + call find_PE_chg(0.0, Kd(K)*dt_h, hp_a(k-1), h(k), & Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), & pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & @@ -1316,6 +1578,31 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif Kddt_h(K) = Kd(K) * dt_h + + elseif (no_MKE_conversion) then ! (PE_chg_max >= 0.0) and use the diffusivity from find_Kd_from_PE_chg. + ! Kd(K) and TKE_used were already set by find_Kd_from_PE_chg. + + ! frac_in_BL = min((TKE_used / PE_chg_g0), 1.0) + if (sfc_connected) MLD_output = MLD_output + frac_in_BL*dz(k) + if (frac_in_BL < 1.0) sfc_disconnect = .true. + + ! Reduce the mechanical and convective TKE proportionately. + TKE_reduc = 0.0 ! tot_TKE could be 0 if Convectively_stable is false. + if ((tot_TKE > 0.0) .and. (tot_TKE > TKE_used)) TKE_reduc = (tot_TKE - TKE_used) / tot_TKE + + ! All TKE should have been consumed. + if (CS%TKE_diagnostics) then + eCD%dTKE_mixing = eCD%dTKE_mixing - TKE_used * I_dtdiag + eCD%dTKE_conv_decay = eCD%dTKE_conv_decay + & + (1.0-TKE_reduc)*(CS%nstar-nstar_FC) * conv_PErel * I_dtdiag + endif + + tot_TKE = tot_TKE - TKE_used + mech_TKE = TKE_reduc*mech_TKE + conv_PErel = TKE_reduc*conv_PErel + + Kddt_h(K) = Kd(K) * dt_h + elseif (tot_TKE + (MKE_src - PE_chg_g0) >= 0.0) then ! This column is convectively stable and there is energy to support the suggested ! mixing. Keep that estimate. @@ -1366,19 +1653,19 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif do itt=1,max_itt if (CS%orig_PE_calc) then - call find_PE_chg_orig(Kddt_h_guess, h(k), hp_a, dTe_term, dSe_term, & + call find_PE_chg_orig(Kddt_h_guess, h(k), hp_a(k-1), dTe_term, dSe_term, & dT_km1_t2, dS_km1_t2, dT_to_dPE(k), dS_to_dPE(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), & pres_Z(K), dT_to_dColHt(k), dS_to_dColHt(k), & dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & PE_chg=PE_chg, dPEc_dKd=dPEc_dKd ) else - call find_PE_chg(0.0, Kddt_h_guess, hp_a, h(k), & + call find_PE_chg(0.0, Kddt_h_guess, hp_a(k-1), h(k), & Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE(k), dS_to_dPE(k), & pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & dT_to_dColHt(k), dS_to_dColHt(k), & - PE_chg=dPE_conv, dPEc_dKd=dPEc_dKd) + PE_chg=PE_chg, dPEc_dKd=dPEc_dKd) endif MKE_src = dMKE_max * (1.0 - exp(-MKE2_Hharm * Kddt_h_guess)) dMKE_src_dK = dMKE_max * MKE2_Hharm * exp(-MKE2_Hharm * Kddt_h_guess) @@ -1445,14 +1732,14 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, Kddt_h(K) = Kd(K) * dt_h ! At this point, the final value of Kddt_h(K) is known, so the ! estimated properties for layer k-1 can be calculated. - b1 = 1.0 / (hp_a + Kddt_h(K)) + b1 = 1.0 / (hp_a(k-1) + Kddt_h(K)) c1(K) = Kddt_h(K) * b1 if (CS%orig_PE_calc) then dTe(k-1) = b1 * ( Kddt_h(K)*(T0(k)-T0(k-1)) + dTe_t2 ) dSe(k-1) = b1 * ( Kddt_h(K)*(S0(k)-S0(k-1)) + dSe_t2 ) endif - hp_a = h(k) + (hp_a * b1) * Kddt_h(K) + hp_a(k) = h(k) + (hp_a(k-1) * b1) * Kddt_h(K) dT_to_dPE_a(k) = dT_to_dPE(k) + c1(K)*dT_to_dPE_a(k-1) dS_to_dPE_a(k) = dS_to_dPE(k) + c1(K)*dS_to_dPE_a(k-1) dT_to_dColHt_a(k) = dT_to_dColHt(k) + c1(K)*dT_to_dColHt_a(k-1) @@ -1476,7 +1763,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, endif if (calc_Te) then - if (k==2) then + if (K==2) then Te(1) = b1*(h(1)*T0(1)) Se(1) = b1*(h(1)*S0(1)) else @@ -1489,7 +1776,7 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, if (debug) then ! Complete the tridiagonal solve for Te. - b1 = 1.0 / hp_a + b1 = 1.0 / hp_a(nz) Te(nz) = b1 * (h(nz) * T0(nz) + Kddt_h(nz) * Te(nz-1)) Se(nz) = b1 * (h(nz) * S0(nz) + Kddt_h(nz) * Se(nz-1)) dT_expect(nz) = Te(nz) - T0(nz) ; dS_expect(nz) = Se(nz) - S0(nz) @@ -1508,25 +1795,39 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, enddo mixing_debug = dPE_debug * I_dtdiag endif - k = nz ! This is here to allow a breakpoint to be set. - !/BGR - ! The following lines are used for the iteration - ! note the iteration has been altered to use the value predicted by - ! the TKE threshold (ML_DEPTH). This is because the MSTAR - ! is now dependent on the ML, and therefore the ML needs to be estimated - ! more precisely than the grid spacing. - - !New method uses ML_DEPTH as computed in ePBL routine + + if (OBL_it >= CS%Max_MLD_Its) exit + + ! The following lines are used for the iteration. Note the iteration has been altered + ! to use the value predicted by the TKE threshold (ML_DEPTH). This is because the MSTAR + ! is now dependent on the ML, and therefore the ML needs to be estimated more precisely + ! than the grid spacing. + + ! New method uses ML_DEPTH as computed in ePBL routine MLD_found = MLD_output - if (MLD_found - MLD_guess > CS%MLD_tol) then - min_MLD = MLD_guess ; dMLD_min = MLD_found - MLD_guess - elseif (abs(MLD_found - MLD_guess) < CS%MLD_tol) then - OBL_converged = .true. ! Break convergence loop - else ! We know this guess was too deep - max_MLD = MLD_guess ; dMLD_max = MLD_found - MLD_guess ! < -CS%MLD_tol + + ! Find out whether to do another iteration and the new bounds on it. + if (CS%MLD_iter_bug) then + ! There is a bug in the logic here if (MLD_found - MLD_guess == CS%MLD_tol). + if (MLD_found - MLD_guess > CS%MLD_tol) then + min_MLD = MLD_guess ; dMLD_min = MLD_found - MLD_guess + elseif (abs(MLD_found - MLD_guess) < CS%MLD_tol) then + exit ! Break the MLD convergence loop + else ! We know this guess was too deep + max_MLD = MLD_guess ; dMLD_max = MLD_found - MLD_guess ! < -CS%MLD_tol + endif + else + if (abs(MLD_found - MLD_guess) < CS%MLD_tol) then + exit ! Break the MLD convergence loop + elseif (MLD_found > MLD_guess) then ! This guess was too shallow + min_MLD = MLD_guess ; dMLD_min = MLD_found - MLD_guess + else ! We know this guess was too deep + max_MLD = MLD_guess ; dMLD_max = MLD_found - MLD_guess ! < -CS%MLD_tol + endif endif - if (.not.OBL_converged) then ; if (CS%MLD_bisection) then + if (OBL_it < CS%Max_MLD_Its) then + if (CS%MLD_bisection) then ! For the next pass, guess the average of the minimum and maximum values. MLD_guess = 0.5*(min_MLD + max_MLD) else ! Try using the false position method or the returned value instead of simple bisection. @@ -1535,22 +1836,18 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, ! Both bounds have valid change estimates and are probably in the range of possible outputs. MLD_guess = (dMLD_min*max_MLD - dMLD_max*min_MLD) / (dMLD_min - dMLD_max) elseif ((MLD_found > min_MLD) .and. (MLD_found < max_MLD)) then - ! The output MLD_found is an interesting guess, as it likely to bracket the true solution + ! The output MLD_found is an interesting guess, as it is likely to bracket the true solution ! along with the previous value of MLD_guess and to be close to the solution. MLD_guess = MLD_found else ! Bisect if the other guesses would be out-of-bounds. This does not happen much. MLD_guess = 0.5*(min_MLD + max_MLD) endif - endif ; endif - endif - if ((OBL_converged) .or. (OBL_it==CS%Max_MLD_Its)) then - if (report_avg_its) then - CS%sum_its(1) = CS%sum_its(1) + real_to_EFP(real(OBL_it)) - CS%sum_its(2) = CS%sum_its(2) + real_to_EFP(1.0) endif - exit endif + enddo ! Iteration loop for converged boundary layer thickness. + + eCD%OBL_its = min(OBL_it, CS%Max_MLD_Its) if (CS%Use_LT) then eCD%LA = LA ; eCD%LAmod = LAmod ; eCD%mstar = mstar_total ; eCD%mstar_LT = mstar_LT else @@ -1561,150 +1858,1159 @@ subroutine ePBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, TKE_forcing, end subroutine ePBL_column -!> This subroutine calculates the change in potential energy and or derivatives -!! for several changes in an interface's diapycnal diffusivity times a timestep. -subroutine find_PE_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, & - dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, & - pres_Z, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, & - PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0, PE_ColHt_cor) - real, intent(in) :: Kddt_h0 !< The previously used diffusivity at an interface times - !! the time step and divided by the average of the - !! thicknesses around the interface [H ~> m or kg m-2]. - real, intent(in) :: dKddt_h !< The trial change in the diffusivity at an interface times - !! the time step and divided by the average of the - !! thicknesses around the interface [H ~> m or kg m-2]. - real, intent(in) :: hp_a !< The effective pivot thickness of the layer above the - !! interface, given by h_k plus a term that - !! is a fraction (determined from the tridiagonal solver) of - !! Kddt_h for the interface above [H ~> m or kg m-2]. - real, intent(in) :: hp_b !< The effective pivot thickness of the layer below the - !! interface, given by h_k plus a term that - !! is a fraction (determined from the tridiagonal solver) of - !! Kddt_h for the interface above [H ~> m or kg m-2]. - real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer - !! above, including implicit mixing effects with other - !! yet higher layers [C H ~> degC m or degC kg m-2]. - real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer - !! above, including implicit mixing effects with other - !! yet higher layers [S H ~> ppt m or ppt kg m-2]. - real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer - !! below, including implicit mixing effects with other - !! yet lower layers [C H ~> degC m or degC kg m-2]. - real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer - !! below, including implicit mixing effects with other - !! yet lower layers [S H ~> ppt m or ppt kg m-2]. - real, intent(in) :: dT_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating - !! a layer's temperature change to the change in column potential - !! energy, including all implicit diffusive changes in the - !! temperatures of all the layers above [R Z3 T-2 C-1 ~> J m-2 degC-1]. - real, intent(in) :: dS_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating - !! a layer's salinity change to the change in column potential - !! energy, including all implicit diffusive changes in the - !! salinities of all the layers above [R Z3 T-2 S-1 ~> J m-2 ppt-1]. - real, intent(in) :: dT_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating - !! a layer's temperature change to the change in column potential - !! energy, including all implicit diffusive changes in the - !! temperatures of all the layers below [R Z3 T-2 C-1 ~> J m-2 degC-1]. - real, intent(in) :: dS_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating - !! a layer's salinity change to the change in column potential - !! energy, including all implicit diffusive changes in the - !! salinities of all the layers below [R Z3 T-2 S-1 ~> J m-2 ppt-1]. - real, intent(in) :: pres_Z !< The rescaled hydrostatic interface pressure, which relates - !! the changes in column thickness to the energy that is radiated - !! as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3]. - real, intent(in) :: dT_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dT) relating - !! a layer's temperature change to the change in column - !! height, including all implicit diffusive changes - !! in the temperatures of all the layers above [Z C-1 ~> m degC-1]. - real, intent(in) :: dS_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dS) relating - !! a layer's salinity change to the change in column - !! height, including all implicit diffusive changes - !! in the salinities of all the layers above [Z S-1 ~> m ppt-1]. - real, intent(in) :: dT_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dT) relating - !! a layer's temperature change to the change in column - !! height, including all implicit diffusive changes - !! in the temperatures of all the layers below [Z C-1 ~> m degC-1]. - real, intent(in) :: dS_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dS) relating - !! a layer's salinity change to the change in column - !! height, including all implicit diffusive changes - !! in the salinities of all the layers below [Z S-1 ~> m ppt-1]. - - real, intent(out) :: PE_chg !< The change in column potential energy from applying - !! Kddt_h at the present interface [R Z3 T-2 ~> J m-2]. - real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h - !! [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. - real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could - !! be realized by applying a huge value of Kddt_h at the - !! present interface [R Z3 T-2 ~> J m-2]. - real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with Kddt_h in the - !! limit where Kddt_h = 0 [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. - real, optional, intent(out) :: PE_ColHt_cor !< The correction to PE_chg that is made due to a net - !! change in the column height [R Z3 T-2 ~> J m-2]. - - ! Local variables - real :: hps ! The sum of the two effective pivot thicknesses [H ~> m or kg m-2]. - real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term [H2 ~> m2 or kg2 m-4]. - real :: dT_c ! The core term in the expressions for the temperature changes [C H2 ~> degC m2 or degC kg2 m-4]. - real :: dS_c ! The core term in the expressions for the salinity changes [S H2 ~> ppt m2 or ppt kg2 m-4]. - real :: PEc_core ! The diffusivity-independent core term in the expressions - ! for the potential energy changes [R Z2 T-2 ~> J m-3]. - real :: ColHt_core ! The diffusivity-independent core term in the expressions - ! for the column height changes [H Z ~> m2 or kg m-1]. - real :: ColHt_chg ! The change in the column height [H ~> m or kg m-2]. - real :: y1_3 ! A local temporary term in [H-3 ~> m-3 or m6 kg-3]. - real :: y1_4 ! A local temporary term in [H-4 ~> m-4 or m8 kg-4]. - - ! The expression for the change in potential energy used here is derived - ! from the expression for the final estimates of the changes in temperature - ! and salinities, and then extensively manipulated to get it into its most - ! succinct form. The derivation is not necessarily obvious, but it demonstrably - ! works by comparison with separate calculations of the energy changes after - ! the tridiagonal solver for the final changes in temperature and salinity are - ! applied. - - hps = hp_a + hp_b - bdt1 = hp_a * hp_b + Kddt_h0 * hps - dT_c = hp_a * Th_b - hp_b * Th_a - dS_c = hp_a * Sh_b - hp_b * Sh_a - PEc_core = hp_b * (dT_to_dPE_a * dT_c + dS_to_dPE_a * dS_c) - & - hp_a * (dT_to_dPE_b * dT_c + dS_to_dPE_b * dS_c) - ColHt_core = hp_b * (dT_to_dColHt_a * dT_c + dS_to_dColHt_a * dS_c) - & - hp_a * (dT_to_dColHt_b * dT_c + dS_to_dColHt_b * dS_c) - - ! Find the change in column potential energy due to the change in the - ! diffusivity at this interface by dKddt_h. - y1_3 = dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) - PE_chg = PEc_core * y1_3 - ColHt_chg = ColHt_core * y1_3 - if (ColHt_chg < 0.0) PE_chg = PE_chg - pres_Z * ColHt_chg - if (present(PE_ColHt_cor)) PE_ColHt_cor = -pres_Z * min(ColHt_chg, 0.0) +!> This subroutine determines the diffusivities from a bottom boundary layer version of +!! the integrated energetics mixed layer model for a single column of water. +subroutine ePBL_BBL_column(h, dz, u, v, T0, S0, dSV_dT, dSV_dS, SpV_dt, absf, & + dt, Kd, BBL_TKE_in, u_star_BBL, Kd_BBL, BBLD_io, mixvel_BBL, mixlen_BBL, GV, US, CS, eCD) + type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure. + real, dimension(SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]. + real, dimension(SZK_(GV)), intent(in) :: dz !< The vertical distance across layers [Z ~> m]. + real, dimension(SZK_(GV)), intent(in) :: u !< Zonal velocities interpolated to h points + !! [L T-1 ~> m s-1]. + real, dimension(SZK_(GV)), intent(in) :: v !< Zonal velocities interpolated to h points + !! [L T-1 ~> m s-1]. + real, dimension(SZK_(GV)), intent(in) :: T0 !< The initial layer temperatures [C ~> degC]. + real, dimension(SZK_(GV)), intent(in) :: S0 !< The initial layer salinities [S ~> ppt]. - if (present(dPEc_dKd)) then - ! Find the derivative of the potential energy change with dKddt_h. - y1_4 = 1.0 / (bdt1 + dKddt_h * hps)**2 - dPEc_dKd = PEc_core * y1_4 - ColHt_chg = ColHt_core * y1_4 - if (ColHt_chg < 0.0) dPEc_dKd = dPEc_dKd - pres_Z * ColHt_chg - endif + real, dimension(SZK_(GV)), intent(in) :: dSV_dT !< The partial derivative of in-situ specific + !! volume with potential temperature + !! [R-1 C-1 ~> m3 kg-1 degC-1]. + real, dimension(SZK_(GV)), intent(in) :: dSV_dS !< The partial derivative of in-situ specific + !! volume with salinity [R-1 S-1 ~> m3 kg-1 ppt-1]. + real, dimension(SZK_(GV)+1), intent(in) :: SpV_dt !< Specific volume interpolated to interfaces + !! divided by dt (if non-Boussinesq) or + !! 1.0 / (dt * Rho0), in [R-1 T-1 ~> m3 kg-1 s-1], + !! used to convert local TKE into a turbulence + !! velocity cubed. + real, intent(in) :: absf !< The absolute value of the Coriolis parameter [T-1 ~> s-1]. + real, intent(in) :: dt !< Time increment [T ~> s]. + real, dimension(SZK_(GV)+1), & + intent(in) :: Kd !< The diffusivities at interfaces due to previously + !! applied mixing processes [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(in) :: BBL_TKE_in !< The mechanically generated turbulent + !! kinetic energy available for bottom boundary + !! layer mixing within a time step [R Z3 T-2 ~> J m-2]. + real, intent(in) :: u_star_BBL !< The bottom boundary layer friction velocity + !! in thickuness flux units [H T-1 ~> m s-1 or kg m-2 s-1] + real, dimension(SZK_(GV)+1), & + intent(out) :: Kd_BBL !< The bottom boundary layer contribution to diffusivities + !! at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(inout) :: BBLD_io !< A first guess at the bottom boundary layer depth on input, and + !! the calculated bottom boundary layer depth on output [Z ~> m] + real, dimension(SZK_(GV)+1), & + intent(out) :: mixvel_BBL !< The profile of boundary layer turbulent mixing + !! velocities [Z T-1 ~> m s-1] + real, dimension(SZK_(GV)+1), & + intent(out) :: mixlen_BBL !< The profile of bottom boundary layer turbulent + !! mixing lengths [Z ~> m] + type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type + type(energetic_PBL_CS), intent(in) :: CS !< Energetic PBL control structure + type(ePBL_column_diags), intent(inout) :: eCD !< A container for passing around diagnostics. - if (present(dPE_max)) then - ! This expression is the limit of PE_chg for infinite dKddt_h. - y1_3 = 1.0 / (bdt1 * hps) - dPE_max = PEc_core * y1_3 - ColHt_chg = ColHt_core * y1_3 - if (ColHt_chg < 0.0) dPE_max = dPE_max - pres_Z * ColHt_chg - endif +! This subroutine determines the contributions from diffusivities in a single column from a +! bottom-boundary layer adaptation of the integrated energetics planetary boundary layer (ePBL) +! model. It accounts for the possibility that the surface boundary diffusivities have already +! been determined. All calculations are done implicitly, and there is no stability limit on the +! time step. Only mechanical mixing in the bottom boundary layer is considered. (Geothermal heat +! fluxes are addressed elsewhere in the MOM6 code, and convection throughout the water column is +! handled by the surface version of ePBL.) There is no conversion of released mean kinetic energy +! into bottom boundary layer turbulent kinetic energy (at least for now), apart from the explicit +! energy that is supplied as an argument to this routine. - if (present(dPEc_dKd_0)) then - ! This expression is the limit of dPEc_dKd for dKddt_h = 0. - y1_4 = 1.0 / bdt1**2 - dPEc_dKd_0 = PEc_core * y1_4 - ColHt_chg = ColHt_core * y1_4 - if (ColHt_chg < 0.0) dPEc_dKd_0 = dPEc_dKd_0 - pres_Z * ColHt_chg - endif + ! Local variables + real, dimension(SZK_(GV)+1) :: & + pres_Z, & ! Interface pressures with a rescaling factor to convert interface height + ! movements into changes in column potential energy [R Z2 T-2 ~> kg m-1 s-2]. + dztop_dztot ! The distance from the surface divided by the thickness of the + ! water column [nondim]. + real :: mech_BBL_TKE ! The mechanically generated turbulent kinetic energy available for + ! bottom boundary layer mixing within a time step [R Z3 T-2 ~> J m-2]. + real :: TKE_eff_avail ! The turbulent kinetic energy that is effectively available to drive mixing + ! after any effects of exponentially decay have been taken into account + ! [R Z3 T-2 ~> J m-2] + real :: TKE_eff_used ! The amount of TKE_eff_avail that has been used to drive mixing [R Z3 T-2 ~> J m-2] + real :: htot ! The total thickness of the layers above an interface [H ~> m or kg m-2]. + real :: dztot ! The total depth of the layers above an interface [Z ~> m]. + real :: uhtot ! The depth integrated zonal velocities in the layers above [H L T-1 ~> m2 s-1 or kg m-1 s-1] + real :: vhtot ! The depth integrated meridional velocities in the layers above [H L T-1 ~> m2 s-1 or kg m-1 s-1] + real :: Idecay_len_TKE ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1]. + real :: dz_sum ! The total thickness of the water column [Z ~> m]. + + real, dimension(SZK_(GV)) :: & + dT_to_dColHt, & ! Partial derivative of the total column height with the temperature changes + ! within a layer [Z C-1 ~> m degC-1]. + dS_to_dColHt, & ! Partial derivative of the total column height with the salinity changes + ! within a layer [Z S-1 ~> m ppt-1]. + dT_to_dPE, & ! Partial derivatives of column potential energy with the temperature + ! changes within a layer, in [R Z3 T-2 C-1 ~> J m-2 degC-1]. + dS_to_dPE, & ! Partial derivatives of column potential energy with the salinity changes + ! within a layer, in [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + dT_to_dColHt_a, & ! Partial derivative of the total column height with the temperature changes + ! within a layer, including the implicit effects of mixing with layers higher + ! in the water column [Z C-1 ~> m degC-1]. + dS_to_dColHt_a, & ! Partial derivative of the total column height with the salinity changes + ! within a layer, including the implicit effects of mixing with layers higher + ! in the water column [Z S-1 ~> m ppt-1]. + dT_to_dPE_a, & ! Partial derivatives of column potential energy with the temperature changes + ! within a layer, including the implicit effects of mixing with layers higher + ! in the water column [R Z3 T-2 C-1 ~> J m-2 degC-1]. + dS_to_dPE_a, & ! Partial derivative of column potential energy with the salinity changes + ! within a layer, including the implicit effects of mixing with layers higher + ! in the water column [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + dT_to_dColHt_b, & ! Partial derivative of the total column height with the temperature changes + ! within a layer, including the implicit effects of mixing with layers deeper + ! in the water column [Z C-1 ~> m degC-1]. + dS_to_dColHt_b, & ! Partial derivative of the total column height with the salinity changes + ! within a layer, including the implicit effects of mixing with layers deeper + ! in the water column [Z S-1 ~> m ppt-1]. + dT_to_dPE_b, & ! Partial derivatives of column potential energy with the temperature changes + ! within a layer, including the implicit effects of mixing with layers deeper + ! in the water column [R Z3 T-2 C-1 ~> J m-2 degC-1]. + dS_to_dPE_b, & ! Partial derivative of column potential energy with the salinity changes + ! within a layer, including the implicit effects of mixing with layers deeper + ! in the water column [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + c1, & ! c1 is used by the tridiagonal solver [nondim]. + Te, & ! Estimated final values of T in the column [C ~> degC]. + Se, & ! Estimated final values of S in the column [S ~> ppt]. + dTe, & ! Running (1-way) estimates of temperature change [C ~> degC]. + dSe, & ! Running (1-way) estimates of salinity change [S ~> ppt]. + hp_a, & ! An effective pivot thickness of the layer including the effects + ! of coupling with layers above [H ~> m or kg m-2]. This is the first term + ! in the denominator of b1 in a downward-oriented tridiagonal solver. + hp_b, & ! An effective pivot thickness of the layer including the effects + ! of coupling with layers below [H ~> m or kg m-2]. This is the first term + ! in the denominator of b1 in an upward-oriented tridiagonal solver. + Th_a, & ! An effective temperature times a thickness in the layer above, including implicit + ! mixing effects with other yet higher layers [C H ~> degC m or degC kg m-2]. + Sh_a, & ! An effective salinity times a thickness in the layer above, including implicit + ! mixing effects with other yet higher layers [S H ~> ppt m or ppt kg m-2]. + Th_b, & ! An effective temperature times a thickness in the layer below, including implicit + ! mixing effects with other yet lower layers [C H ~> degC m or degC kg m-2]. + Sh_b ! An effective salinity times a thickness in the layer below, including implicit + ! mixing effects with other yet lower layers [S H ~> ppt m or ppt kg m-2]. + real, dimension(SZK_(GV)+1) :: & + MixLen_shape, & ! A nondimensional shape factor for the mixing length that + ! gives it an appropriate asymptotic value at the bottom of + ! the boundary layer [nondim]. + h_dz_int, & ! The ratio of the layer thicknesses over the vertical distances + ! across the layers surrounding an interface [H Z-1 ~> nondim or kg m-3] + Kddt_h ! The total diapycnal diffusivity at an interface times a timestep divided by the + ! average thicknesses around an interface [H ~> m or kg m-2]. + real :: b1 ! b1 is inverse of the pivot used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1]. + real :: h_neglect ! A thickness that is so small it is usually lost + ! in roundoff and can be neglected [H ~> m or kg m-2]. + real :: dz_neglect ! A vertical distance that is so small it is usually lost + ! in roundoff and can be neglected [Z ~> m]. + real :: dMass ! The mass per unit area within a layer [Z R ~> kg m-2]. + real :: dPres ! The hydrostatic pressure change across a layer [R Z2 T-2 ~> Pa = J m-3]. + + real :: dt_h ! The timestep divided by the averages of the vertical distances around + ! a layer [T Z-1 ~> s m-1]. + real :: dz_top ! The distance from the surface [Z ~> m]. + real :: dz_rsum ! The distance from the seafloor [Z ~> m]. + real :: I_dzsum ! The inverse of dz_sum [Z-1 ~> m-1]. + real :: I_BBLD ! The inverse of the current value of BBLD [Z-1 ~> m-1]. + real :: dz_tt ! The distance from the surface or up to the next interface + ! that did not exhibit turbulent mixing from this scheme plus + ! a surface mixing roughness length given by dz_tt_min [Z ~> m]. + real :: dz_tt_min ! A surface roughness length [Z ~> m]. + real :: C1_3 ! = 1/3 [nondim] + real :: vstar ! An in-situ turbulent velocity [Z T-1 ~> m s-1]. + real :: BBLD_output ! The bottom boundary layer depth output from this routine [Z ~> m] + real :: hbs_here ! The local minimum of dztop_dztot and MixLen_shape [nondim] + real :: TKE_used ! The TKE used to support mixing at an interface [R Z3 T-2 ~> J m-2]. + real :: Kd_guess0 ! A first guess of the diapycnal diffusivity [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real :: PE_chg_g0 ! The potential energy change when Kd is Kd_guess0 [R Z3 T-2 ~> J m-2] + real :: Kddt_h_g0 ! The first guess of the change in diapycnal diffusivity times a timestep + ! divided by the average thicknesses around an interface [H ~> m or kg m-2]. + real :: Kddt_h_prev ! The diapycnal diffusivity before modification times a timestep divided + ! by the average thicknesses around an interface [H ~> m or kg m-2]. + real :: dPEc_dKd_Kd0 ! The partial derivative of PE change with Kddt_h(K) + ! for very small values of Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. + real :: PE_chg ! The change in potential energy due to mixing at an + ! interface [R Z3 T-2 ~> J m-2], positive for the column increasing + ! in potential energy (i.e., consuming TKE). + real :: TKE_left ! The amount of turbulent kinetic energy left for the most + ! recent guess at Kddt_h(K) [R Z3 T-2 ~> J m-2]. + real :: dPEc_dKd ! The partial derivative of PE_chg with Kddt_h(K) [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. + real :: TKE_left_min, TKE_left_max ! Maximum and minimum values of TKE_left [R Z3 T-2 ~> J m-2]. + real :: Kddt_h_max, Kddt_h_min ! Maximum and minimum values of Kddt_h(K) [H ~> m or kg m-2]. + real :: Kddt_h_guess ! A guess at the value of Kddt_h(K) [H ~> m or kg m-2]. + real :: Kddt_h_next ! The next guess at the value of Kddt_h(K) [H ~> m or kg m-2]. + real :: dKddt_h ! The change between guesses at Kddt_h(K) [H ~> m or kg m-2]. + real :: dKddt_h_Newt ! The change between guesses at Kddt_h(K) with Newton's method [H ~> m or kg m-2]. + real :: Kddt_h_newt ! The Newton's method next guess for Kddt_h(K) [H ~> m or kg m-2]. + real :: exp_kh ! The nondimensional decay of TKE across a layer [nondim]. + real :: frac_in_BL ! The fraction of the energy required to support dKd_max that is suppiled by + ! max_PE_chg, used here to determine a fractional layer contribution to the + ! boundary layer thickness [nondim] + real :: TKE_rescale ! The effective fractional increase in energy available to + ! mixing at this interface once its exponential decay is accounted for [nondim] + logical :: use_Newt ! Use Newton's method for the next guess at Kddt_h(K). + logical :: convectively_stable ! If true the water column is convectively stable at this interface. + logical :: bot_connected ! If true the ocean is actively turbulent from the present + ! interface all the way down to the seafloor. + logical :: bot_disconnect ! If true, any turbulence has become disconnected + ! from the bottom. + + ! The following is only used for diagnostics. + real :: I_dtdiag ! = 1.0 / dt [T-1 ~> s-1]. + + real :: BBLD_guess, BBLD_found ! Bottom boundary layer depth guessed/found for iteration [Z ~> m] + real :: min_BBLD, max_BBLD ! Iteration bounds on BBLD [Z ~> m], which are adjusted at each step + real :: dBBLD_min ! The change in diagnosed mixed layer depth when the guess is min_BLD [Z ~> m] + real :: dBBLD_max ! The change in diagnosed mixed layer depth when the guess is max_BLD [Z ~> m] + logical :: BBL_converged ! Flag for convergence of BBLD + integer :: BBL_it ! Iteration counter + + real :: Surface_Scale ! Surface decay scale for vstar [nondim] + logical :: debug ! This is used as a hard-coded value for debugging. + logical :: no_MKE_conversion ! If true, there is conversion of MKE to TKE in this routine. + + ! The following arrays are used only for debugging purposes. + real :: dPE_debug ! An estimate of the potential energy change [R Z3 T-2 ~> J m-2] + real :: mixing_debug ! An estimate of the rate of change of potential energy due to mixing [R Z3 T-3 ~> W m-2] + real, dimension(20) :: TKE_left_itt ! The value of TKE_left after each iteration [R Z3 T-2 ~> J m-2] + real, dimension(20) :: PE_chg_itt ! The value of PE_chg after each iteration [R Z3 T-2 ~> J m-2] + real, dimension(20) :: Kddt_h_itt ! The value of Kddt_h_guess after each iteration [H ~> m or kg m-2] + real, dimension(20) :: dPEa_dKd_itt ! The value of dPEc_dKd after each iteration [R Z3 T-2 H-1 ~> J m-3 or J kg-1] +! real, dimension(20) :: MKE_src_itt ! The value of MKE_src after each iteration [R Z3 T-2 ~> J m-2] + real, dimension(SZK_(GV)) :: dT_expect !< Expected temperature changes [C ~> degC] + real, dimension(SZK_(GV)) :: dS_expect !< Expected salinity changes [S ~> ppt] + real, dimension(SZK_(GV)) :: mech_BBL_TKE_k ! The mechanically generated turbulent kinetic energy + ! available for bottom boundary mixing over a time step for each layer [R Z3 T-2 ~> J m-2]. + integer, dimension(SZK_(GV)) :: num_itts + + integer :: k, nz, itt, max_itt + + nz = GV%ke + + debug = .false. ! Change this hard-coded value for debugging. + no_MKE_conversion = ((CS%direct_calc) ) ! .and. (CS%MKE_to_TKE_effic == 0.0)) + + ! Add bottom boundary layer mixing if there is energy to support it. + if ((CS%ePBL_BBL_effic <= 0.0) .or. (BBL_TKE_in <= 0.0)) then + ! There is no added bottom boundary layer mixing. + BBLD_io = 0.0 + Kd_BBL(:) = 0.0 + mixvel_BBL(:) = 0.0 ; mixlen_BBL(:) = 0.0 + eCD%BBL_its = 0 + if (CS%TKE_diagnostics) then + eCD%dTKE_BBL_mixing = 0.0 ; eCD%dTKE_BBL_decay = 0.0 ; eCD%dTKE_BBL = 0.0 + ! eCD%dTKE_BBL_MKE = 0.0 + endif + return + else + ! There will be added bottom boundary layer mixing. + + h_neglect = GV%H_subroundoff + dz_neglect = GV%dZ_subroundoff + + C1_3 = 1.0 / 3.0 + I_dtdiag = 1.0 / dt + max_itt = 20 + dz_tt_min = 0.0 + + ! The next two blocks of code could be shared with ePBL_column. + + ! Set up fields relating a layer's temperature and salinity changes to potential energy changes. + pres_Z(1) = 0.0 + do k=1,nz + dMass = GV%H_to_RZ * h(k) + dPres = US%L_to_Z**2 * GV%g_Earth * dMass + dT_to_dPE(k) = (dMass * (pres_Z(K) + 0.5*dPres)) * dSV_dT(k) + dS_to_dPE(k) = (dMass * (pres_Z(K) + 0.5*dPres)) * dSV_dS(k) + dT_to_dColHt(k) = dMass * dSV_dT(k) + dS_to_dColHt(k) = dMass * dSV_dS(k) + + pres_Z(K+1) = pres_Z(K) + dPres + enddo + + if (GV%Boussinesq) then + do K=1,nz+1 ; h_dz_int(K) = GV%Z_to_H ; enddo + else + h_dz_int(1) = (h(1) + h_neglect) / (dz(1) + dz_neglect) + do K=2,nz + h_dz_int(K) = (h(k-1) + h(k) + h_neglect) / (dz(k-1) + dz(k) + dz_neglect) + enddo + h_dz_int(nz+1) = (h(nz) + h_neglect) / (dz(nz) + dz_neglect) + endif + ! The two previous blocks of code could be shared with ePBL_column. + + ! Determine the total thickness (dz_sum) and the fractional distance from the top (dztop_dztot). + dz_sum = 0.0 ; do k=1,nz ; dz_sum = dz_sum + dz(k) ; enddo + I_dzsum = 0.0 ; if (dz_sum > 0.0) I_dzsum = 1.0 / dz_sum + dz_top = 0.0 + dztop_dztot(nz+1) = 0.0 + do k=1,nz + dz_top = dz_top + dz(k) + dztop_dztot(K) = dz_top * I_dzsum + enddo + + ! Set terms from a tridiagonal solver based on the previously determined diffusivities. + Kddt_h(1) = 0.0 + hp_a(1) = h(1) + dT_to_dPE_a(1) = dT_to_dPE(1) ; dT_to_dColHt_a(1) = dT_to_dColHt(1) + dS_to_dPE_a(1) = dS_to_dPE(1) ; dS_to_dColHt_a(1) = dS_to_dColHt(1) + do K=2,nz + dt_h = dt / max(0.5*(dz(k-1)+dz(k)), 1e-15*dz_sum) + Kddt_h(K) = Kd(K) * dt_h + b1 = 1.0 / (hp_a(k-1) + Kddt_h(K)) + c1(K) = Kddt_h(K) * b1 + hp_a(k) = h(k) + (hp_a(k-1) * b1) * Kddt_h(K) + dT_to_dPE_a(k) = dT_to_dPE(k) + c1(K)*dT_to_dPE_a(k-1) + dS_to_dPE_a(k) = dS_to_dPE(k) + c1(K)*dS_to_dPE_a(k-1) + dT_to_dColHt_a(k) = dT_to_dColHt(k) + c1(K)*dT_to_dColHt_a(k-1) + dS_to_dColHt_a(k) = dS_to_dColHt(k) + c1(K)*dS_to_dColHt_a(k-1) + if (K<=2) then + Te(k-1) = b1*(h(k-1)*T0(k-1)) ; Se(k-1) = b1*(h(k-1)*S0(k-1)) + Th_a(k-1) = h(k-1) * T0(k-1) ; Sh_a(k-1) = h(k-1) * S0(k-1) + else + Te(k-1) = b1 * (h(k-1) * T0(k-1) + Kddt_h(K-1) * Te(k-2)) + Se(k-1) = b1 * (h(k-1) * S0(k-1) + Kddt_h(K-1) * Se(k-2)) + Th_a(k-1) = h(k-1) * T0(k-1) + Kddt_h(K-1) * Te(k-2) + Sh_a(k-1) = h(k-1) * S0(k-1) + Kddt_h(K-1) * Se(k-2) + endif + enddo + Kddt_h(nz+1) = 0.0 + if (debug) then + ! Complete the tridiagonal solve for Te and Se, which may be useful for debugging. + b1 = 1.0 / hp_a(nz) + Te(nz) = b1 * (h(nz) * T0(nz) + Kddt_h(nz) * Te(nz-1)) + Se(nz) = b1 * (h(nz) * S0(nz) + Kddt_h(nz) * Se(nz-1)) + do k=nz-1,1,-1 + Te(k) = Te(k) + c1(K+1)*Te(k+1) + Se(k) = Se(k) + c1(K+1)*Se(k+1) + enddo + endif + + BBLD_guess = BBLD_io + + !/The following lines are for the iteration over BBLD + ! max_BBLD will initialized as ocean bottom depth + max_BBLD = 0.0 ; do k=1,nz ; max_BBLD = max_BBLD + dz(k) ; enddo + ! min_BBLD will be initialized to 0. + min_BBLD = 0.0 + ! Set values of the wrong signs to indicate that these changes are not based on valid estimates + dBBLD_min = -1.0*US%m_to_Z ; dBBLD_max = 1.0*US%m_to_Z + + ! If no first guess is provided for BBLD, try the middle of the water column + if (BBLD_guess <= min_BBLD) BBLD_guess = 0.5 * (min_BBLD + max_BBLD) + + ! Iterate to determine a converged EPBL bottom boundary layer depth. + do BBL_it=1,CS%max_BBLD_its + + if (debug) then ; mech_BBL_TKE_k(:) = 0.0 ; endif + + ! Reset BBL_depth + BBLD_output = dz(nz) + bot_connected = .true. + + mech_BBL_TKE = BBL_TKE_in + + if (CS%TKE_diagnostics) then + ! eCD%dTKE_BBL_MKE = 0.0 + eCD%dTKE_BBL_mixing = 0.0 + eCD%dTKE_BBL_decay = 0.0 + eCD%dTKE_BBL = mech_BBL_TKE * I_dtdiag + endif + + ! Store in 1D arrays for output. + do K=1,nz+1 ; mixvel_BBL(K) = 0.0 ; mixlen_BBL(K) = 0.0 ; enddo + + ! Determine the mixing shape function MixLen_shape. + if ((.not.CS%Use_BBLD_iteration) .or. & + (CS%transLay_scale >= 1.0) .or. (CS%transLay_scale < 0.0) ) then + do K=1,nz+1 + MixLen_shape(K) = 1.0 + enddo + elseif (BBLD_guess <= 0.0) then + if (CS%transLay_scale > 0.0) then ; do K=1,nz+1 + MixLen_shape(K) = CS%transLay_scale + enddo ; else ; do K=1,nz+1 + MixLen_shape(K) = 1.0 + enddo ; endif + else + ! Reduce the mixing length based on BBLD, with a quadratic + ! expression that follows KPP. + I_BBLD = 1.0 / BBLD_guess + dz_rsum = 0.0 + MixLen_shape(nz+1) = 1.0 + if (CS%MixLenExponent_BBL==2.0) then + do K=nz,1,-1 + dz_rsum = dz_rsum + dz(k) + MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & + (max(0.0, (BBLD_guess - dz_rsum)*I_BBLD) )**2 + enddo + elseif (CS%MixLenExponent_BBL==1.0) then + do K=nz,1,-1 + dz_rsum = dz_rsum + dz(k) + MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & + (max(0.0, (BBLD_guess - dz_rsum)*I_BBLD) ) + enddo + else ! (CS%MixLenExponent_BBL /= 2.0 or 1.0) then + do K=nz,1,-1 + dz_rsum = dz_rsum + dz(k) + MixLen_shape(K) = CS%transLay_scale + (1.0 - CS%transLay_scale) * & + (max(0.0, (BBLD_guess - dz_rsum)*I_BBLD) )**CS%MixLenExponent_BBL + enddo + endif + endif + + Kd_BBL(nz+1) = 0.0 ; Kddt_h(nz+1) = 0.0 + hp_b(nz) = h(nz) + dT_to_dPE_b(nz) = dT_to_dPE(nz) ; dT_to_dColHt_b(nz) = dT_to_dColHt(nz) + dS_to_dPE_b(nz) = dS_to_dPE(nz) ; dS_to_dColHt_b(nz) = dS_to_dColHt(nz) + + htot = h(nz) ; dztot = dz(nz) ; uhtot = u(nz)*h(nz) ; vhtot = v(nz)*h(nz) + + if (debug) then + mech_BBL_TKE_k(nz) = mech_BBL_TKE + num_itts(:) = -1 + endif + + Idecay_len_TKE = (CS%TKE_decay_BBL * absf) / u_star_BBL + do K=nz,2,-1 + ! Apply dissipation to the TKE, here applied as an exponential decay + ! due to 3-d turbulent energy being lost to inefficient rotational modes. + ! The following form is often used for mechanical stirring from the surface. + ! There could be several different "flavors" of TKE that decay at different rates. + exp_kh = 1.0 + if (Idecay_len_TKE > 0.0) exp_kh = exp(-h(k)*Idecay_len_TKE) + if (CS%TKE_diagnostics) & + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay + (exp_kh-1.0) * mech_BBL_TKE * I_dtdiag + mech_BBL_TKE = mech_BBL_TKE * exp_kh + + if (debug) then + mech_BBL_TKE_k(K) = mech_BBL_TKE + endif + + ! Precalculate some temporary expressions that are independent of Kddt_h(K). + dt_h = dt / max(0.5*(dz(k-1)+dz(k)), 1e-15*dz_sum) + + ! This tests whether the layers above and below this interface are in + ! a convectively stable configuration, without considering any effects of + ! mixing at higher interfaces. It is an approximation to the more + ! complete test dPEc_dKd_Kd0 >= 0.0, that would include the effects of + ! mixing across interface K+1. The dT_to_dColHt here are effectively + ! mass-weighted estimates of dSV_dT. + Convectively_stable = ( 0.0 <= & + ( (dT_to_dColHt(k) + dT_to_dColHt(k-1) ) * (T0(k-1)-T0(k)) + & + (dS_to_dColHt(k) + dS_to_dColHt(k-1) ) * (S0(k-1)-S0(k)) ) ) + + if ((mech_BBL_TKE <= 0.0) .and. Convectively_stable) then + ! Energy is already exhausted, so set Kd_BBL = 0 and cycle or exit? + mech_BBL_TKE = 0.0 + Kd_BBL(K) = 0.0 ; Kddt_h(K) = Kd(K) * dt_h + bot_disconnect = .true. + ! if (.not.debug) exit + + else ! mech_BBL_TKE > 0.0 or this is a potentially convectively unstable profile. + bot_disconnect = .false. + + ! Precalculate some more temporary expressions that are independent of Kddt_h(K). + if (K>=nz) then + Th_b(k) = h(k) * T0(k) ; Sh_b(k) = h(k) * S0(k) + else + Th_b(k) = h(k) * T0(k) + Kddt_h(K+1) * Te(k+1) + Sh_b(k) = h(k) * S0(k) + Kddt_h(K+1) * Se(k+1) + endif + + ! Using Pr=1 and the diffusivity at the upper interface (once it is + ! known), determine how much resolved mean kinetic energy (MKE) will be + ! extracted within a timestep and add a fraction CS%MKE_to_TKE_effic of + ! this to the mTKE budget available for mixing in the next layer. + ! This is not enabled yet for BBL mixing. + ! if ((CS%MKE_to_TKE_effic > 0.0) .and. (htot*h(k-1) > 0.0)) then + ! ! This is the energy that would be available from homogenizing the + ! ! velocities between layer k-1 and the layers below. + ! dMKE_max = (US%L_to_Z**2*GV%H_to_RZ * CS%MKE_to_TKE_effic) * 0.5 * & + ! (h(k-1) / ((htot + h(k-1))*htot)) * & + ! ((uhtot-u(k-1)*htot)**2 + (vhtot-v(k-1)*htot)**2) + ! ! A fraction (1-exp(Kddt_h*MKE2_Hharm)) of this energy would be + ! ! extracted by mixing with a finite viscosity. + ! MKE2_Hharm = (htot + h(k-1) + 2.0*h_neglect) / & + ! ((htot+h_neglect) * (h(k-1)+h_neglect)) + ! else + ! dMKE_max = 0.0 + ! MKE2_Hharm = 0.0 + ! endif + + ! At this point, Kddt_h(K) will be unknown because its value may depend + ! on how much energy is available. + dz_tt = dztot + dz_tt_min + if (mech_BBL_TKE > 0.0) then + if (CS%wT_scheme_BBL==wT_from_cRoot_TKE) then + vstar = CS%vstar_scale_fac_BBL * cuberoot(SpV_dt(K)*mech_BBL_TKE) + elseif (CS%wT_scheme_BBL==wT_from_RH18) then + Surface_Scale = max(0.05, 1.0 - dztot / BBLD_guess) + vstar = (CS%vstar_scale_fac_BBL * Surface_Scale) * ( CS%vstar_surf_fac_BBL*u_star_BBL/h_dz_int(K) ) + endif + hbs_here = min(dztop_dztot(K), MixLen_shape(K)) + mixlen_BBL(K) = max(CS%min_BBL_mix_len, ((dz_tt*hbs_here)*vstar) / & + ((CS%Ekman_scale_coef_BBL * absf) * (dz_tt*hbs_here) + vstar)) + Kd_guess0 = (h_dz_int(K)*vstar) * CS%vonKar * mixlen_BBL(K) + else + vstar = 0.0 ; Kd_guess0 = 0.0 + endif + mixvel_BBL(K) = vstar ! Track vstar + + TKE_rescale = 1.0 + if (CS%decay_adjusted_BBL_TKE) then + ! Add a scaling factor that accounts for the exponential decay of TKE from a + ! near-bottom source and the assumption that an increase in the diffusivity at an + ! interface causes a linearly increasing buoyancy flux going from 0 at the bottom + ! to a peak at the interface, and then going back to 0 atop the layer above. + TKE_rescale = exp_decay_TKE_adjust(htot, h(k-1), Idecay_len_TKE) + endif + + TKE_eff_avail = TKE_rescale*mech_BBL_TKE + + if (no_MKE_conversion) then + ! Without conversion from MKE to TKE, the updated diffusivity can be determined directly. + call find_Kd_from_PE_chg(Kd(K), Kd_guess0, dt_h, TKE_eff_avail, hp_a(k-1), hp_b(k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & + dT_to_dColHt_b(k), dS_to_dColHt_b(k), Kd_add=Kd_BBL(K), PE_chg=TKE_eff_used, & + frac_dKd_max_PE=frac_in_BL) + + ! Do not add energy if the column is convectively unstable. This was handled previously + ! for mixing from the surface. + if (TKE_eff_used < 0.0) TKE_eff_used = 0.0 + + ! Convert back to the TKE that has actually been used. + if (CS%decay_adjusted_BBL_TKE) then + if (TKE_rescale == 0.0) then ! This probably never occurs, even at roundoff. + TKE_used = mech_BBL_TKE ! All the energy was dissipated before it could mix. + else + TKE_used = TKE_eff_used / TKE_rescale + endif + else + TKE_used = TKE_eff_used + endif + + if (bot_connected) BBLD_output = BBLD_output + frac_in_BL*dz(k-1) + if (frac_in_BL < 1.0) bot_disconnect = .true. + + if (CS%TKE_diagnostics) then + eCD%dTKE_BBL_mixing = eCD%dTKE_BBL_mixing - TKE_eff_used * I_dtdiag + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay - (TKE_used-TKE_eff_used) * I_dtdiag + endif + + mech_BBL_TKE = mech_BBL_TKE - TKE_used + + Kddt_h(K) = (Kd(K) + Kd_BBL(K)) * dt_h + + else + Kddt_h_prev = Kd(K) * dt_h + Kddt_h_g0 = Kd_guess0 * dt_h + ! Find the change in PE with the guess at the added bottom boundary layer mixing. + call find_PE_chg(Kddt_h_prev, Kddt_h_g0, hp_a(k-1), hp_b(k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & + dT_to_dColHt_b(k), dS_to_dColHt_b(k), & + PE_chg=PE_chg_g0, dPEc_dKd_0=dPEc_dKd_Kd0 ) + + ! MKE_src = 0.0 ! Enable later?: = dMKE_max*(1.0 - exp(-Kddt_h_g0 * MKE2_Hharm)) + + ! Do not add energy if the column is convectively unstable. This was handled previously + ! for mixing from the surface. + if (PE_chg_g0 < 0.0) PE_chg_g0 = 0.0 + + ! This block checks out different cases to determine Kd at the present interface. + ! if (mech_BBL_TKE*TKE_rescale + (MKE_src - PE_chg_g0) >= 0.0) then + if (TKE_eff_avail - PE_chg_g0 >= 0.0) then + ! This column is convectively stable and there is energy to support the suggested + ! mixing, or it is convectively unstable. Keep this first estimate of Kd. + Kd_BBL(K) = Kd_guess0 + Kddt_h(K) = Kddt_h_prev + Kddt_h_g0 + + TKE_used = PE_chg_g0 / TKE_rescale + + if (CS%TKE_diagnostics) then + eCD%dTKE_BBL_mixing = eCD%dTKE_BBL_mixing - PE_chg_g0 * I_dtdiag +! eCD%dTKE_MKE = eCD%dTKE_MKE + MKE_src * I_dtdiag + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay - (TKE_used - PE_chg_g0) * I_dtdiag + endif + + ! mech_BBL_TKE = mech_BBL_TKE + MKE_src - TKE_used + mech_BBL_TKE = mech_BBL_TKE - TKE_used + if (bot_connected) then + BBLD_output = BBLD_output + dz(k-1) + endif + + elseif (TKE_eff_avail == 0.0) then + ! This can arise if there is no energy input to drive mixing or if there + ! is such strong decay that the mech_BBL_TKE becomes 0 via an underflow. + Kd_BBL(K) = 0.0 ; Kddt_h(K) = Kddt_h_prev + if (CS%TKE_diagnostics) then + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay - mech_BBL_TKE * I_dtdiag + endif + mech_BBL_TKE = 0.0 + bot_disconnect = .true. + else + ! There is not enough energy to support the mixing, so reduce the + ! diffusivity to what can be supported. + Kddt_h_max = Kddt_h_g0 ; Kddt_h_min = 0.0 + ! TKE_left_max = TKE_eff_avail + (MKE_src - PE_chg_g0) + TKE_left_max = TKE_eff_avail - PE_chg_g0 + TKE_left_min = TKE_eff_avail + + ! As a starting guess, take the minimum of a false position estimate + ! and a Newton's method estimate starting from dKddt_h = 0.0 + ! Enable conversion from MKE to TKE in the bottom boundary layer later? + ! Kddt_h_guess = TKE_eff_avail * Kddt_h_max / max( PE_chg_g0 - MKE_src, & + ! Kddt_h_max * (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) ) + Kddt_h_guess = TKE_eff_avail * Kddt_h_max / max( PE_chg_g0, Kddt_h_max * dPEc_dKd_Kd0 ) + ! The above expression is mathematically the same as the following + ! except it is not susceptible to division by zero when + ! dPEc_dKd_Kd0 = dMKE_max = 0 . + ! Kddt_h_guess = TKE_eff_avail * min( Kddt_h_max / (PE_chg_g0 - MKE_src), & + ! 1.0 / (dPEc_dKd_Kd0 - dMKE_max * MKE2_Hharm) ) + if (debug) then + TKE_left_itt(:) = 0.0 ; dPEa_dKd_itt(:) = 0.0 ; PE_chg_itt(:) = 0.0 + Kddt_h_itt(:) = 0.0 ! ; MKE_src_itt(:) = 0.0 + endif + do itt=1,max_itt + call find_PE_chg(Kddt_h_prev, Kddt_h_guess, hp_a(k-1), hp_b(k), & + Th_a(k-1), Sh_a(k-1), Th_b(k), Sh_b(k), & + dT_to_dPE_a(k-1), dS_to_dPE_a(k-1), dT_to_dPE_b(k), dS_to_dPE_b(k), & + pres_Z(K), dT_to_dColHt_a(k-1), dS_to_dColHt_a(k-1), & + dT_to_dColHt_b(k), dS_to_dColHt_b(k), & + PE_chg=PE_chg, dPEc_dKd=dPEc_dKd) + ! Enable conversion from MKE to TKE in the bottom boundary layer later? + ! MKE_src = dMKE_max * (1.0 - exp(-MKE2_Hharm * Kddt_h_guess)) + ! dMKE_src_dK = dMKE_max * MKE2_Hharm * exp(-MKE2_Hharm * Kddt_h_guess) + + ! TKE_left = TKE_eff_avail + (MKE_src - PE_chg) + TKE_left = TKE_eff_avail - PE_chg + if (debug .and. itt<=20) then + Kddt_h_itt(itt) = Kddt_h_guess ! ; MKE_src_itt(itt) = MKE_src + PE_chg_itt(itt) = PE_chg ; dPEa_dKd_itt(itt) = dPEc_dKd + TKE_left_itt(itt) = TKE_left + endif + ! Store the new bounding values, bearing in mind that min and max + ! here refer to Kddt_h and dTKE_left/dKddt_h < 0: + if (TKE_left >= 0.0) then + Kddt_h_min = Kddt_h_guess ; TKE_left_min = TKE_left + else + Kddt_h_max = Kddt_h_guess ; TKE_left_max = TKE_left + endif + + ! Try to use Newton's method, but if it would go outside the bracketed + ! values use the false-position method instead. + use_Newt = .true. + ! if (dPEc_dKd - dMKE_src_dK <= 0.0) then + if (dPEc_dKd <= 0.0) then + use_Newt = .false. + else + ! dKddt_h_Newt = TKE_left / (dPEc_dKd - dMKE_src_dK) + dKddt_h_Newt = TKE_left / dPEc_dKd + Kddt_h_Newt = Kddt_h_guess + dKddt_h_Newt + if ((Kddt_h_Newt > Kddt_h_max) .or. (Kddt_h_Newt < Kddt_h_min)) & + use_Newt = .false. + endif + + if (use_Newt) then + Kddt_h_next = Kddt_h_guess + dKddt_h_Newt + dKddt_h = dKddt_h_Newt + else + Kddt_h_next = (TKE_left_max * Kddt_h_min - Kddt_h_max * TKE_left_min) / & + (TKE_left_max - TKE_left_min) + dKddt_h = Kddt_h_next - Kddt_h_guess + endif + + if ((abs(dKddt_h) < 1e-9*Kddt_h_guess) .or. (itt==max_itt)) then + ! Use the old value so that the energy calculation does not need to be repeated. + if (debug) num_itts(K) = itt + exit + else + Kddt_h_guess = Kddt_h_next + endif + enddo ! Inner iteration loop on itt. + Kd_BBL(K) = Kddt_h_guess / dt_h + Kddt_h(K) = (Kd(K) + Kd_BBL(K)) * dt_h + + ! All TKE should have been consumed. + if (CS%TKE_diagnostics) then + ! eCD%dTKE_BBL_mixing = eCD%dTKE_BBL_mixing - (TKE_eff_avail + MKE_src) * I_dtdiag + ! eCD%dTKE_BBL_MKE = eCD%dTKE_BBL_MKE + MKE_src * I_dtdiag + eCD%dTKE_BBL_mixing = eCD%dTKE_BBL_mixing - TKE_eff_avail * I_dtdiag + eCD%dTKE_BBL_decay = eCD%dTKE_BBL_decay - (mech_BBL_TKE-TKE_eff_avail) * I_dtdiag + endif + + if (bot_connected) BBLD_output = BBLD_output + (PE_chg / PE_chg_g0) * dz(k-1) + + mech_BBL_TKE = 0.0 + bot_disconnect = .true. + endif ! End of convective or forced mixing cases to determine Kd. + endif + + Kddt_h(K) = (Kd(K) + Kd_BBL(K)) * dt_h + endif ! tot_TKT > 0.0 branch. Kddt_h(K) has been set. + + ! At this point, the final value of Kddt_h(K) is known, so the + ! estimated properties for layer k can be calculated. + b1 = 1.0 / (hp_b(k) + Kddt_h(K)) + c1(K) = Kddt_h(K) * b1 + + hp_b(k-1) = h(k-1) + (hp_b(k) * b1) * Kddt_h(K) + dT_to_dPE_b(k-1) = dT_to_dPE(k-1) + c1(K)*dT_to_dPE_b(k) + dS_to_dPE_b(k-1) = dS_to_dPE(k-1) + c1(K)*dS_to_dPE_b(k) + dT_to_dColHt_b(k-1) = dT_to_dColHt(k-1) + c1(K)*dT_to_dColHt_b(k) + dS_to_dColHt_b(k-1) = dS_to_dColHt(k-1) + c1(K)*dS_to_dColHt_b(k) + + ! Store integrated velocities and thicknesses for MKE conversion calculations. + if (bot_disconnect) then + ! There is no turbulence at this interface, so restart the running sums. + uhtot = u(k-1)*h(k-1) + vhtot = v(k-1)*h(k-1) + htot = h(k-1) + dztot = dz(k-1) + bot_connected = .false. + else + uhtot = uhtot + u(k-1)*h(k-1) + vhtot = vhtot + v(k-1)*h(k-1) + htot = htot + h(k-1) + dztot = dztot + dz(k-1) + endif + + if (K==nz) then + Te(k) = b1*(h(k)*T0(k)) + Se(k) = b1*(h(k)*S0(k)) + else + Te(k) = b1 * (h(k) * T0(k) + Kddt_h(K+1) * Te(k+1)) + Se(k) = b1 * (h(k) * S0(k) + Kddt_h(K+1) * Se(k+1)) + endif + enddo + Kd_BBL(1) = 0.0 + + if (debug) then + ! Complete the tridiagonal solve for Te with a downward pass. + b1 = 1.0 / hp_b(1) + Te(1) = b1 * (h(1) * T0(1) + Kddt_h(2) * Te(2)) + Se(1) = b1 * (h(1) * S0(1) + Kddt_h(2) * Se(2)) + dT_expect(1) = Te(1) - T0(1) ; dS_expect(1) = Se(1) - S0(1) + do k=2,nz + Te(k) = Te(k) + c1(K)*Te(k-1) + Se(k) = Se(k) + c1(K)*Se(k-1) + dT_expect(k) = Te(k) - T0(k) ; dS_expect(k) = Se(k) - S0(k) + enddo + + dPE_debug = 0.0 + do k=1,nz + dPE_debug = dPE_debug + (dT_to_dPE(k) * (Te(k) - T0(k)) + & + dS_to_dPE(k) * (Se(k) - S0(k))) + enddo + mixing_debug = dPE_debug * I_dtdiag + endif + + ! Skip the rest of the contents of the do loop if there are no more BBL depth iterations. + if (BBL_it >= CS%max_BBLD_its) exit + + ! The following lines are used for the iteration to determine the boundary layer depth. + ! Note that the iteration uses the value predicted by the TKE threshold (BBL_DEPTH), + ! because the mixing length shape is dependent on the BBL depth, and therefore the BBL depth + ! should be estimated more precisely than the grid spacing. + + ! New method uses BBL_DEPTH as computed in ePBL routine + BBLD_found = BBLD_output + if (abs(BBLD_found - BBLD_guess) < CS%BBLD_tol) then + exit ! Break the BBL depth convergence loop + elseif (BBLD_found > BBLD_guess) then + min_BBLD = BBLD_guess ; dBBLD_min = BBLD_found - BBLD_guess + else ! We know this guess was too deep + max_BBLD = BBLD_guess ; dBBLD_max = BBLD_found - BBLD_guess ! <= -CS%BBLD_tol + endif + + ! Try using the false position method or the returned value instead of simple bisection. + ! Taking the occasional step with BBLD_output empirically helps to converge faster. + if ((dBBLD_min > 0.0) .and. (dBBLD_max < 0.0) .and. (BBL_it > 2) .and. (mod(BBL_it-1,4) > 0)) then + ! Both bounds have valid change estimates and are probably in the range of possible outputs. + BBLD_guess = (dBBLD_min*max_BBLD - dBBLD_max*min_BBLD) / (dBBLD_min - dBBLD_max) + elseif ((BBLD_found > min_BBLD) .and. (BBLD_found < max_BBLD)) then + ! The output BBLD_found is an interesting guess, as it is likely to bracket the true solution + ! along with the previous value of BBLD_guess and to be close to the solution. + BBLD_guess = BBLD_found + else ! Bisect if the other guesses would be out-of-bounds. This does not happen much. + BBLD_guess = 0.5*(min_BBLD + max_BBLD) + endif + + enddo ! Iteration loop for converged boundary layer thickness. + + eCD%BBL_its = min(BBL_it, CS%max_BBLD_its) + + BBLD_io = BBLD_output + endif + +end subroutine ePBL_BBL_column + +!> Determine a scaling factor that accounts for the exponential decay of turbulent kinetic energy +!! from a boundary source and the assumption that an increase in the diffusivity at an interface +!! causes a linearly increasing buoyancy flux going from 0 at the bottom to a peak at the interface, +!! and then going back to 0 atop the layer above. Where this factor increases the available mixing +!! TKE, it is only compensating for the fact that the TKE has already been reduced by the same +!! exponential decay rate. ha and hb must be non-negative, and this function generally increases +!! with hb and decreases with ha. +!! +!! Exp_decay_TKE_adjust is coded to have a lower bound of 1e-30 on the return value. For large +!! values of ha*Idecay, the return value is about 0.5*ka*(ha+hb)*Idecay**2 * exp(-ha*Idecay), but +!! return values of less than 1e-30 are deliberately reset to 1e-30. For relatively large values +!! of hb*Idecay, the return value increases linearly with hb. When Idecay ~= 0, the return value +!! is close to 1. +function exp_decay_TKE_adjust(hb, ha, Idecay) result(TKE_to_PE_scale) + real, intent(in) :: hb !< The thickness over which the buoyancy flux varies on the + !! near-boundary side of an interface (e.g., a well-mixed bottom + !! boundary layer thickness) [H ~> m or kg m-2] + real, intent(in) :: ha !< The thickness of the layer on the opposite side of an interface from + !! the boundary [H ~> m or kg m-2] + real, intent(in) :: Idecay !< The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1] + real :: TKE_to_PE_scale !< The effective fractional change in energy available to + !! drive mixing at this interface once the exponential decay of TKE + !! is accounted for [nondim]. TKE_to_PE_scale is always positive. + + real :: khb ! The thickness on the boundary side times the TKE decay rate [nondim] + real :: kha ! The thickness away from from the boundary times the TKE decay rate [nondim] + real, parameter :: C1_3 = 1.0/3.0 ! A rational constant [nondim] + + khb = abs(hb*Idecay) + kha = abs(ha*Idecay) + + ! For large enough kha that exp(kha) > 1.0e17*kha: + ! TKE_to_PE_scale = (0.5 * (khb + kha) * kha) * exp(-kha) > (0.5 * kha**2) * exp(-kha) + ! To keep TKE_to_PE_scale > -1e30 and avoid overflow in the exp(), keep kha < kha_max_30, where: + ! kha_max_30 = ln(0.5*1e30) + 2.0 * ln(kha_max_30) ~= 68.3844 + 2.0 * ln(68.3844+8.6895)) + ! If kha_max = 77.0739, (0.5 * kha_max**2) * exp(-kha_max) = 1.0e-30. + + if (kha > 77.0739) then + TKE_to_PE_scale = 1.0e-30 + elseif ((kha > 2.2e-4) .and. (khb > 2.2e-4)) then + ! This is the usual case, derived from integrals of z exp(z) over the layers above and below. + ! TKE_to_PE_scale = (0.5 * (khb + kha)) / & + ! ((exp(-khb) - (1.0 - khb)) / khb + (exp(kha) - (1.0 + kha)) / kha) + TKE_to_PE_scale = (0.5 * (khb + kha) * (kha * khb)) / & + (kha * (exp(-khb) - (1.0 - khb)) + khb * (exp(kha) - (1.0 + kha))) + elseif (khb > 2.2e-4) then + ! For small values of kha, approximate (exp(kha) - (1.0 + hha)) by the first two + ! terms of its Taylor series: 0.5*kha**2 + C1_6*kha**3 + ... + kha**n/n! + ... + ! which is more accurate when kha**4/24. < 1e-16 or kha < ~ 2.21e-4. + TKE_to_PE_scale = (0.5 * (khb + kha) * khb) / & + ((exp(-khb) - (1.0 - khb)) + 0.5*(khb * kha) * (1.0 + C1_3*kha)) + elseif (kha > 2.2e-4) then + ! Use a Taylor series expansion for small values of khb + TKE_to_PE_scale = (0.5 * (khb + kha) * kha) / & + (0.5 * (kha * khb) * (1.0 - C1_3*Khb) + (exp(kha) - (1.0 + kha))) + else ! (kha < 2.2e-4) .and. (khb < 2.2e-4) - use Taylor series approximations for both + TKE_to_PE_scale = 1.0 / (1.0 + C1_3*(kha - khb)) + endif + + if (TKE_to_PE_scale < 1.0e-30) TKE_to_PE_scale = 1.0e-30 + + ! For kha >> 1: + ! TKE_to_PE_scale = (0.5 * (khb + kha) * kha) * exp(-kha) + + ! For khb >> 1: + ! TKE_to_PE_scale = (0.5 * (khb + kha) * (kha * khb)) / & + ! (khb * exp(kha) - (kha + khb))) + ! For khb >> 1 and khb >> kha: + ! TKE_to_PE_scale = (0.5 * (kha * khb)) / (exp(kha) - 1.0)) + +end function exp_decay_TKE_adjust + +!> This subroutine calculates the change in potential energy and or derivatives +!! for several changes in an interface's diapycnal diffusivity times a timestep. +subroutine find_PE_chg(Kddt_h0, dKddt_h, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, & + dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, & + pres_Z, dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, & + PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0, PE_ColHt_cor) + real, intent(in) :: Kddt_h0 !< The previously used diffusivity at an interface times + !! the time step and divided by the average of the + !! thicknesses around the interface [H ~> m or kg m-2]. + real, intent(in) :: dKddt_h !< The trial change in the diffusivity at an interface times + !! the time step and divided by the average of the + !! thicknesses around the interface [H ~> m or kg m-2]. + real, intent(in) :: hp_a !< The effective pivot thickness of the layer above the + !! interface, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface above [H ~> m or kg m-2]. + real, intent(in) :: hp_b !< The effective pivot thickness of the layer below the + !! interface, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface below [H ~> m or kg m-2]. + real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers [C H ~> degC m or degC kg m-2]. + real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers [S H ~> ppt m or ppt kg m-2]. + real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers [C H ~> degC m or degC kg m-2]. + real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers [S H ~> ppt m or ppt kg m-2]. + real, intent(in) :: dT_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column potential + !! energy, including all implicit diffusive changes in the + !! temperatures of all the layers above [R Z3 T-2 C-1 ~> J m-2 degC-1]. + real, intent(in) :: dS_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column potential + !! energy, including all implicit diffusive changes in the + !! salinities of all the layers above [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + real, intent(in) :: dT_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column potential + !! energy, including all implicit diffusive changes in the + !! temperatures of all the layers below [R Z3 T-2 C-1 ~> J m-2 degC-1]. + real, intent(in) :: dS_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column potential + !! energy, including all implicit diffusive changes in the + !! salinities of all the layers below [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + real, intent(in) :: pres_Z !< The rescaled hydrostatic interface pressure, which relates + !! the changes in column thickness to the energy that is radiated + !! as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3]. + real, intent(in) :: dT_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers above [Z C-1 ~> m degC-1]. + real, intent(in) :: dS_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers above [Z S-1 ~> m ppt-1]. + real, intent(in) :: dT_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers below [Z C-1 ~> m degC-1]. + real, intent(in) :: dS_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers below [Z S-1 ~> m ppt-1]. + + real, intent(out) :: PE_chg !< The change in column potential energy from applying + !! dKddt_h at the present interface [R Z3 T-2 ~> J m-2]. + real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with dKddt_h + !! [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. + real, optional, intent(out) :: dPE_max !< The maximum change in column potential energy that could + !! be realized by applying a huge value of dKddt_h at the + !! present interface [R Z3 T-2 ~> J m-2]. + real, optional, intent(out) :: dPEc_dKd_0 !< The partial derivative of PE_chg with dKddt_h in the + !! limit where dKddt_h = 0 [R Z3 T-2 H-1 ~> J m-3 or J kg-1]. + real, optional, intent(out) :: PE_ColHt_cor !< The correction to PE_chg that is made due to a net + !! change in the column height [R Z3 T-2 ~> J m-2]. + + ! Local variables + real :: hps ! The sum of the two effective pivot thicknesses [H ~> m or kg m-2]. + real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term [H2 ~> m2 or kg2 m-4]. + real :: dT_c ! The core term in the expressions for the temperature changes [C H2 ~> degC m2 or degC kg2 m-4]. + real :: dS_c ! The core term in the expressions for the salinity changes [S H2 ~> ppt m2 or ppt kg2 m-4]. + real :: PEc_core ! The diffusivity-independent core term in the expressions + ! for the potential energy changes [R Z2 T-2 ~> J m-3]. + real :: ColHt_core ! The diffusivity-independent core term in the expressions + ! for the column height changes [H Z ~> m2 or kg m-1]. + real :: ColHt_chg ! The change in the column height [H ~> m or kg m-2]. + real :: y1_3 ! A local temporary term in [H-3 ~> m-3 or m6 kg-3]. + real :: y1_4 ! A local temporary term in [H-4 ~> m-4 or m8 kg-4]. + + ! The expression for the change in potential energy used here is derived + ! from the expression for the final estimates of the changes in temperature + ! and salinities, and then extensively manipulated to get it into its most + ! succinct form. The derivation is not necessarily obvious, but it demonstrably + ! works by comparison with separate calculations of the energy changes after + ! the tridiagonal solver for the final changes in temperature and salinity are + ! applied. + + hps = hp_a + hp_b + bdt1 = hp_a * hp_b + Kddt_h0 * hps + dT_c = hp_a * Th_b - hp_b * Th_a + dS_c = hp_a * Sh_b - hp_b * Sh_a + PEc_core = hp_b * (dT_to_dPE_a * dT_c + dS_to_dPE_a * dS_c) - & + hp_a * (dT_to_dPE_b * dT_c + dS_to_dPE_b * dS_c) + ColHt_core = hp_b * (dT_to_dColHt_a * dT_c + dS_to_dColHt_a * dS_c) - & + hp_a * (dT_to_dColHt_b * dT_c + dS_to_dColHt_b * dS_c) + + ! Find the change in column potential energy due to the change in the + ! diffusivity at this interface by dKddt_h. + y1_3 = dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) + PE_chg = PEc_core * y1_3 + ColHt_chg = ColHt_core * y1_3 + if (ColHt_chg < 0.0) PE_chg = PE_chg - pres_Z * ColHt_chg + + if (present(PE_ColHt_cor)) PE_ColHt_cor = -pres_Z * min(ColHt_chg, 0.0) + + if (present(dPEc_dKd)) then + ! Find the derivative of the potential energy change with dKddt_h. + y1_4 = 1.0 / (bdt1 + dKddt_h * hps)**2 + dPEc_dKd = PEc_core * y1_4 + ColHt_chg = ColHt_core * y1_4 + if (ColHt_chg < 0.0) dPEc_dKd = dPEc_dKd - pres_Z * ColHt_chg + endif + + if (present(dPE_max)) then + ! This expression is the limit of PE_chg for infinite dKddt_h. + y1_3 = 1.0 / (bdt1 * hps) + dPE_max = PEc_core * y1_3 + ColHt_chg = ColHt_core * y1_3 + if (ColHt_chg < 0.0) dPE_max = dPE_max - pres_Z * ColHt_chg + endif + + if (present(dPEc_dKd_0)) then + ! This expression is the limit of dPEc_dKd for dKddt_h = 0. + y1_4 = 1.0 / bdt1**2 + dPEc_dKd_0 = PEc_core * y1_4 + ColHt_chg = ColHt_core * y1_4 + if (ColHt_chg < 0.0) dPEc_dKd_0 = dPEc_dKd_0 - pres_Z * ColHt_chg + endif + +end subroutine find_PE_chg + + +!> This subroutine directly calculates the an increment in the diapycnal diffusivity based on the +!! change in potential energy within a timestep, subject to bounds on the possible change in +!! diffusivity, returning both the added diffusivity and the realized potential energy change, and +!! optionally also the maximum change in potential energy that would be realized for an infinitely +!! large diffusivity. +subroutine find_Kd_from_PE_chg(Kd_prev, dKd_max, dt_h, max_PE_chg, hp_a, hp_b, Th_a, Sh_a, Th_b, Sh_b, & + dT_to_dPE_a, dS_to_dPE_a, dT_to_dPE_b, dS_to_dPE_b, pres_Z, & + dT_to_dColHt_a, dS_to_dColHt_a, dT_to_dColHt_b, dS_to_dColHt_b, & + Kd_add, PE_chg, dPE_max, frac_dKd_max_PE) + real, intent(in) :: Kd_prev !< The previously used diffusivity at an interface + !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(in) :: dKd_max !< The maximum change in the diffusivity at an interface + !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(in) :: dt_h !< The time step and divided by the average of the + !! thicknesses around the interface [T Z-1 ~> s m-1]. + real, intent(in) :: max_PE_chg !< The maximum change in the column potential energy due to + !! additional mixing at an interface [R Z3 T-2 ~> J m-2]. + + real, intent(in) :: hp_a !< The effective pivot thickness of the layer above the + !! interface, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface above [H ~> m or kg m-2]. + real, intent(in) :: hp_b !< The effective pivot thickness of the layer below the + !! interface, given by h_k plus a term that + !! is a fraction (determined from the tridiagonal solver) of + !! Kddt_h for the interface below [H ~> m or kg m-2]. + real, intent(in) :: Th_a !< An effective temperature times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers [C H ~> degC m or degC kg m-2]. + real, intent(in) :: Sh_a !< An effective salinity times a thickness in the layer + !! above, including implicit mixing effects with other + !! yet higher layers [S H ~> ppt m or ppt kg m-2]. + real, intent(in) :: Th_b !< An effective temperature times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers [C H ~> degC m or degC kg m-2]. + real, intent(in) :: Sh_b !< An effective salinity times a thickness in the layer + !! below, including implicit mixing effects with other + !! yet lower layers [S H ~> ppt m or ppt kg m-2]. + real, intent(in) :: dT_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column potential + !! energy, including all implicit diffusive changes in the + !! temperatures of all the layers above [R Z3 T-2 C-1 ~> J m-2 degC-1]. + real, intent(in) :: dS_to_dPE_a !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column potential + !! energy, including all implicit diffusive changes in the + !! salinities of all the layers above [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + real, intent(in) :: dT_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dT) relating + !! a layer's temperature change to the change in column potential + !! energy, including all implicit diffusive changes in the + !! temperatures of all the layers below [R Z3 T-2 C-1 ~> J m-2 degC-1]. + real, intent(in) :: dS_to_dPE_b !< A factor (pres_lay*mass_lay*dSpec_vol/dS) relating + !! a layer's salinity change to the change in column potential + !! energy, including all implicit diffusive changes in the + !! salinities of all the layers below [R Z3 T-2 S-1 ~> J m-2 ppt-1]. + real, intent(in) :: pres_Z !< The rescaled hydrostatic interface pressure, which relates + !! the changes in column thickness to the energy that is radiated + !! as gravity waves and unavailable to drive mixing [R Z2 T-2 ~> J m-3]. + real, intent(in) :: dT_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers above [Z C-1 ~> m degC-1]. + real, intent(in) :: dS_to_dColHt_a !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers above [Z S-1 ~> m ppt-1]. + real, intent(in) :: dT_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dT) relating + !! a layer's temperature change to the change in column + !! height, including all implicit diffusive changes + !! in the temperatures of all the layers below [Z C-1 ~> m degC-1]. + real, intent(in) :: dS_to_dColHt_b !< A factor (mass_lay*dSColHtc_vol/dS) relating + !! a layer's salinity change to the change in column + !! height, including all implicit diffusive changes + !! in the salinities of all the layers below [Z S-1 ~> m ppt-1]. + real, intent(out) :: Kd_add !< The additional diffusivity at an interface + !! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]. + real, intent(out) :: PE_chg !< The realized change in the column potential energy due to + !! additional mixing at an interface [R Z3 T-2 ~> J m-2]. + real, optional, & + intent(out) :: dPE_max !< The maximum change in column potential energy that could + !! be realized by applying a huge value of dKddt_h at the + !! present interface [R Z3 T-2 ~> J m-2]. + real, optional, & + intent(out) :: frac_dKd_max_PE !< The fraction of the energy required to support dKd_max + !! that is supplied by max_PE_chg [nondim] + + ! Local variables + real :: Kddt_h0 ! The previously used diffusivity at an interface times the time step + ! and divided by the average of the thicknesses around the + ! interface [H ~> m or kg m-2]. + real :: dKddt_h ! The upper bound on the change in the diffusivity at an interface times + ! the time step and divided by the average of the thicknesses around + ! the interface [H ~> m or kg m-2]. + real :: hps ! The sum of the two effective pivot thicknesses [H ~> m or kg m-2]. + real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term [H2 ~> m2 or kg2 m-4]. + real :: dT_c ! The core term in the expressions for the temperature changes [C H2 ~> degC m2 or degC kg2 m-4]. + real :: dS_c ! The core term in the expressions for the salinity changes [S H2 ~> ppt m2 or ppt kg2 m-4]. + real :: PEc_core ! The diffusivity-independent core term in the expressions + ! for the potential energy changes [R Z2 T-2 ~> J m-3]. + real :: ColHt_core ! The diffusivity-independent core term in the expressions + ! for the column height changes [H Z ~> m2 or kg m-1]. + + ! The expression for the change in potential energy used here is derived from the expression + ! for the final estimates of the changes in temperature and salinities, which is then + ! extensively manipulated to get it into its most succinct form. It is the same as the + ! expression that appears in find_PE_chg. + + Kddt_h0 = Kd_prev * dt_h + hps = hp_a + hp_b + bdt1 = hp_a * hp_b + Kddt_h0 * hps + dT_c = hp_a * Th_b - hp_b * Th_a + dS_c = hp_a * Sh_b - hp_b * Sh_a + PEc_core = hp_b * (dT_to_dPE_a * dT_c + dS_to_dPE_a * dS_c) - & + hp_a * (dT_to_dPE_b * dT_c + dS_to_dPE_b * dS_c) + ColHt_core = hp_b * (dT_to_dColHt_a * dT_c + dS_to_dColHt_a * dS_c) - & + hp_a * (dT_to_dColHt_b * dT_c + dS_to_dColHt_b * dS_c) + if (ColHt_core < 0.0) PEc_core = PEc_core - pres_Z * ColHt_core + + ! Find the change in column potential energy due to the change in the + ! diffusivity at this interface by dKd_max, and use this to dermine which limit applies. + dKddt_h = dKd_max * dt_h + if ( (PEc_core * dKddt_h <= max_PE_chg * (bdt1 * (bdt1 + dKddt_h * hps))) .or. (PEc_core <= 0.0) ) then + ! There is more than enough energy available to support the maximum permitted diffusivity. + Kd_add = dKd_max + PE_chg = PEc_core * dKddt_h / (bdt1 * (bdt1 + dKddt_h * hps)) + if (present(frac_dKd_max_PE)) frac_dKd_max_PE = 1.0 + else + ! Mixing is constrained by the available energy, so solve the following for Kd_add: + ! max_PE_chg = PEc_core * Kd_add * dt_h / (bdt1 * (bdt1 + Kd_add * dt_h * hps)) + ! It has been verified that the two branches are continuous. + Kd_add = (bdt1**2 * max_PE_chg) / (dt_h * (PEc_core - bdt1 * hps * max_PE_chg)) + PE_chg = max_PE_chg + if (present(frac_dKd_max_PE)) & + frac_dKd_max_PE = (PE_chg * (bdt1 * (bdt1 + dKddt_h * hps))) / (PEc_core * dKddt_h) + endif + + ! Note that the derivative of PE_chg with dKddt_h is monotonic: + ! dPE_chg_dKd = PEc_core * ( (bdt1 * (bdt1 + dKddt_h * hps)) - bdtl * hps * dKddt_h ) / & + ! (bdt1 * (bdt1 + dKddt_h * hps))**2 + ! dPE_chg_dKd = PEc_core / (bdt1 + dKddt_h * hps)**2 + + ! This expression is the limit of PE_chg for infinite dKddt_h. + if (present(dPE_max)) dPE_max = PEc_core / (bdt1 * hps) + +end subroutine find_Kd_from_PE_chg -end subroutine find_PE_chg !> This subroutine calculates the change in potential energy and or derivatives !! for several changes in an interface's diapycnal diffusivity times a timestep @@ -2074,12 +3380,15 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) ! This include declares and sets the variable "version". # include "version_variable.h" character(len=40) :: mdl = "MOM_energetic_PBL" ! This module's name. - character(len=20) :: tmpstr + character(len=20) :: tmpstr ! A string that is parsed for parameter settings + character(len=20) :: vel_scale_str ! A string that is parsed for velocity scale parameter settings + character(len=120) :: diff_text ! A clause describing parameter setting that differ. real :: omega_frac_dflt ! The default for omega_frac [nondim] integer :: isd, ied, jsd, jed integer :: mstar_mode, LT_enhance, wT_mode integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. - logical :: use_temperature, use_omega + logical :: use_omega + logical :: no_BBL ! If true, EPBL_BBL_EFFIC < 0 and bottom boundary layer mixing is not enabled. logical :: use_la_windsea isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed @@ -2092,6 +3401,9 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) !/1. General ePBL settings + call get_param(param_file, mdl, "DEBUG", CS%debug, & + "If true, write out verbose debugging data.", & + default=.false., debuggingParam=.true.) call get_param(param_file, mdl, "OMEGA", CS%omega, & "The rotation rate of the earth.", & units="s-1", default=7.2921e-5, scale=US%T_to_S) @@ -2101,7 +3413,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "scale for turbulence.", default=.false., do_not_log=.true.) omega_frac_dflt = 0.0 if (use_omega) then - call MOM_error(WARNING, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.") + call MOM_error(WARNING, "ML_USE_OMEGA is deprecated; use ML_OMEGA_FRAC=1.0 instead.") omega_frac_dflt = 1.0 endif call get_param(param_file, mdl, "ML_OMEGA_FRAC", CS%omega_frac, & @@ -2130,10 +3442,10 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701) call get_param(param_file, mdl, "EPBL_ORIGINAL_PE_CALC", CS%orig_PE_calc, & - "If true, the ePBL code uses the original form of the "//& - "potential energy change code. Otherwise, the newer "//& - "version that can work with successive increments to the "//& - "diffusivity in upward or downward passes is used.", default=.true.) + "If true, the ePBL code uses the original form of the potential energy change "//& + "code. Otherwise, the newer version that can work with successive increments "//& + "to the diffusivity in upward or downward passes is used.", & + default=.true.) ! Change the default to .false.? call get_param(param_file, mdl, "MKE_TO_TKE_EFFIC", CS%MKE_to_TKE_effic, & "The efficiency with which mean kinetic energy released "//& @@ -2141,16 +3453,21 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "is converted to turbulent kinetic energy.", & units="nondim", default=0.0) call get_param(param_file, mdl, "TKE_DECAY", CS%TKE_decay, & - "TKE_DECAY relates the vertical rate of decay of the "//& - "TKE available for mechanical entrainment to the natural "//& - "Ekman depth.", units="nondim", default=2.5) + "TKE_DECAY relates the vertical rate of decay of the TKE available "//& + "for mechanical entrainment to the natural Ekman depth.", & + units="nondim", default=2.5) + call get_param(param_file, mdl, "DIRECT_EPBL_MIXING_CALC", CS%direct_calc, & + "If true and there is no conversion from mean kinetic energy to ePBL turbulent "//& + "kinetic energy, use a direct calculation of the diffusivity that is supported "//& + "by a given energy input instead of the more general but slower iterative solver.", & + default=.false., do_not_log=(CS%MKE_to_TKE_effic>0.0)) !/2. Options related to setting MSTAR call get_param(param_file, mdl, "EPBL_MSTAR_SCHEME", tmpstr, & "EPBL_MSTAR_SCHEME selects the method for setting mstar. Valid values are: \n"//& "\t CONSTANT - Use a fixed mstar given by MSTAR \n"//& - "\t OM4 - Use L_Ekman/L_Obukhov in the sabilizing limit, as in OM4 \n"//& + "\t OM4 - Use L_Ekman/L_Obukhov in the stabilizing limit, as in OM4 \n"//& "\t REICHL_H18 - Use the scheme documented in Reichl & Hallberg, 2018.", & default=CONSTANT_STRING, do_not_log=.true.) call get_param(param_file, mdl, "MSTAR_MODE", mstar_mode, default=-1) @@ -2209,7 +3526,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) units="nondim", default=0.085, do_not_log=(CS%mstar_scheme/=MStar_from_Ekman)) ! mstar_scheme==MStar_from_RH18 options call get_param(param_file, mdl, "RH18_MSTAR_CN1", CS%RH18_mstar_cn1,& - "MSTAR_N coefficient 1 (outter-most coefficient for fit). "//& + "MSTAR_N coefficient 1 (outer-most coefficient for fit). "//& "The value of 0.275 is given in RH18. Increasing this "//& "coefficient increases MSTAR for all values of Hf/ust, but more "//& "effectively at low values (weakly developed OSBLs).", & @@ -2276,10 +3593,14 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "mixed layer depth. Otherwise use the false position after a maximum and minimum "//& "bound have been evaluated and the returned value or bisection before this.", & default=.false., do_not_log=.not.CS%Use_MLD_iteration) + call get_param(param_file, mdl, "EPBL_MLD_ITER_BUG", CS%MLD_iter_bug, & + "If true, use buggy logic that gives the wrong bounds for the next iteration "//& + "when successive guesses increase by exactly EPBL_MLD_TOLERANCE.", & + default=.true., do_not_log=.not.CS%Use_MLD_iteration) ! The default should be changed to .false. call get_param(param_file, mdl, "EPBL_MLD_MAX_ITS", CS%max_MLD_its, & "The maximum number of iterations that can be used to find a self-consistent "//& "mixed layer depth. If EPBL_MLD_BISECTION is true, the maximum number "//& - "iteractions needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.", & + "of iterations needed is set by Depth/2^MAX_ITS < EPBL_MLD_TOLERANCE.", & default=20, do_not_log=.not.CS%Use_MLD_iteration) if (.not.CS%Use_MLD_iteration) CS%Max_MLD_Its = 1 call get_param(param_file, mdl, "EPBL_MIN_MIX_LEN", CS%min_mix_len, & @@ -2294,7 +3615,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) units="nondim", default=2.0) !/ Turbulent velocity scale in mixing coefficient - call get_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", tmpstr, & + call get_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", vel_scale_str, & "Selects the method for translating TKE into turbulent velocities. "//& "Valid values are: \n"//& "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//& @@ -2303,31 +3624,31 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) default=ROOT_TKE_STRING, do_not_log=.true.) call get_param(param_file, mdl, "EPBL_VEL_SCALE_MODE", wT_mode, default=-1) if (wT_mode == 0) then - tmpstr = ROOT_TKE_STRING + vel_scale_str = ROOT_TKE_STRING call MOM_error(WARNING, "Use EPBL_VEL_SCALE_SCHEME = CUBE_ROOT_TKE instead of the archaic EPBL_VEL_SCALE_MODE = 0.") elseif (wT_mode == 1) then - tmpstr = RH18_STRING + vel_scale_str = RH18_STRING call MOM_error(WARNING, "Use EPBL_VEL_SCALE_SCHEME = REICHL_H18 instead of the archaic EPBL_VEL_SCALE_MODE = 1.") elseif (wT_mode >= 2) then call MOM_error(FATAL, "An unrecognized value of the obsolete parameter EPBL_VEL_SCALE_MODE was specified.") endif - call log_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", tmpstr, & + call log_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", vel_scale_str, & "Selects the method for translating TKE into turbulent velocities. "//& "Valid values are: \n"//& "\t CUBE_ROOT_TKE - A constant times the cube root of remaining TKE. \n"//& "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//& "\t documented in Reichl & Hallberg, 2018.", & default=ROOT_TKE_STRING) - tmpstr = uppercase(tmpstr) - select case (tmpstr) + vel_scale_str = uppercase(vel_scale_str) + select case (vel_scale_str) case (ROOT_TKE_STRING) CS%wT_scheme = wT_from_cRoot_TKE case (RH18_STRING) CS%wT_scheme = wT_from_RH18 case default - call MOM_mesg('energetic_PBL_init: EPBL_VEL_SCALE_SCHEME ="'//trim(tmpstr)//'"', 0) + call MOM_mesg('energetic_PBL_init: EPBL_VEL_SCALE_SCHEME ="'//trim(vel_scale_str)//'"', 0) call MOM_error(FATAL, "energetic_PBL_init: Unrecognized setting "// & - "EPBL_VEL_SCALE_SCHEME = "//trim(tmpstr)//" found in input file.") + "EPBL_VEL_SCALE_SCHEME = "//trim(vel_scale_str)//" found in input file.") end select call get_param(param_file, mdl, "WSTAR_USTAR_COEF", CS%wstar_ustar_coef, & @@ -2343,6 +3664,77 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "The proportionality times ustar to set vstar at the surface.", & units="nondim", default=1.2) + !/ Bottom boundary layer mixing related options + call get_param(param_file, mdl, "EPBL_BBL_EFFIC", CS%ePBL_BBL_effic, & + "The efficiency of bottom boundary layer mixing via ePBL. Setting this to a "//& + "value that is greater than 0 to enable bottom boundary layer mixing from EPBL.", & + units="nondim", default=0.0) + no_BBL = (CS%ePBL_BBL_effic<=0.0) + call get_param(param_file, mdl, "USE_BBLD_ITERATION", CS%Use_BBLD_iteration, & + "A logical that specifies whether or not to use the distance to the top of the "//& + "actively turbulent bottom boundary layer to help set the EPBL length scale.", & + default=.true., do_not_log=no_BBL) + call get_param(param_file, mdl, "TKE_DECAY_BBL", CS%TKE_decay_BBL, & + "TKE_DECAY_BBL relates the vertical rate of decay of the TKE available for "//& + "mechanical entrainment in the bottom boundary layer to the natural Ekman depth.", & + units="nondim", default=CS%TKE_decay, do_not_log=no_BBL) + call get_param(param_file, mdl, "MIX_LEN_EXPONENT_BBL", CS%MixLenExponent_BBL, & + "The exponent applied to the ratio of the distance to the top of the BBL "//& + "and the total BBL depth which determines the shape of the mixing length. "//& + "This is only used if USE_MLD_ITERATION is True.", & + units="nondim", default=2.0, do_not_log=(no_BBL.or.(.not.CS%Use_BBLD_iteration))) + call get_param(param_file, mdl, "EPBL_MIN_BBL_MIX_LEN", CS%min_BBL_mix_len, & + "The minimum mixing length scale that will be used by ePBL for bottom boundary "//& + "layer mixing. Choosing (0) does not set a minimum.", & + units="meter", default=CS%min_mix_len, scale=US%m_to_Z, do_not_log=no_BBL) + call get_param(param_file, mdl, "EPBL_BBLD_TOLERANCE", CS%BBLD_tol, & + "The tolerance for the iteratively determined bottom boundary layer depth. "//& + "This is only used with USE_MLD_ITERATION.", & + units="meter", default=US%Z_to_m*CS%MLD_tol, scale=US%m_to_Z, & + do_not_log=(no_BBL.or.(.not.CS%Use_MLD_iteration))) + call get_param(param_file, mdl, "EPBL_BBLD_MAX_ITS", CS%max_BBLD_its, & + "The maximum number of iterations that can be used to find a self-consistent "//& + "bottom boundary layer depth.", & + default=CS%max_MLD_its, do_not_log=(no_BBL.or.(.not.CS%Use_MLD_iteration))) + if (.not.CS%Use_MLD_iteration) CS%max_BBLD_its = 1 + + call get_param(param_file, mdl, "EPBL_BBL_VEL_SCALE_SCHEME", tmpstr, & + "Selects the method for translating bottom boundary layer TKE into turbulent velocities. "//& + "Valid values are: \n"//& + "\t CUBE_ROOT_TKE - A constant times the cube root of remaining BBL TKE. \n"//& + "\t REICHL_H18 - Use the scheme based on a combination of w* and v* as \n"//& + "\t documented in Reichl & Hallberg, 2018.", & + default=vel_scale_str, do_not_log=no_BBL) + select case (tmpstr) + case (ROOT_TKE_STRING) + CS%wT_scheme_BBL = wT_from_cRoot_TKE + case (RH18_STRING) + CS%wT_scheme_BBL = wT_from_RH18 + case default + call MOM_mesg('energetic_PBL_init: EPBL_BBL_VEL_SCALE_SCHEME ="'//trim(tmpstr)//'"', 0) + call MOM_error(FATAL, "energetic_PBL_init: Unrecognized setting "// & + "EPBL_BBL_VEL_SCALE_SCHEME = "//trim(tmpstr)//" found in input file.") + end select + call get_param(param_file, mdl, "EPBL_BBL_VEL_SCALE_FACTOR", CS%vstar_scale_fac_BBL, & + "An overall nondimensional scaling factor for wT in the bottom boundary layer. "//& + "Making this larger increases the bottom boundary layer diffusivity.", & + units="nondim", default=CS%vstar_scale_fac, do_not_log=no_BBL) + call get_param(param_file, mdl, "VSTAR_BBL_SURF_FAC", CS%vstar_surf_fac_BBL,& + "The proportionality times ustar to set vstar in the bottom boundary layer.", & + units="nondim", default=CS%vstar_surf_fac, do_not_log=(no_BBL.or.(CS%wT_scheme_BBL/=wT_from_RH18))) + call get_param(param_file, mdl, "EKMAN_SCALE_COEF_BBL", CS%Ekman_scale_coef_BBL, & + "A nondimensional scaling factor controlling the inhibition of the diffusive "//& + "length scale by rotation in the bottom boundary layer. Making this larger "//& + "decreases the bottom boundary layer diffusivity.", & + units="nondim", default=CS%Ekman_scale_coef, do_not_log=no_BBL) + + call get_param(param_file, mdl, "DECAY_ADJUSTED_BBL_TKE", CS%decay_adjusted_BBL_TKE, & + "If true, include an adjustment factor in the bottom boundary layer energetics "//& + "that accounts for an exponential decay of TKE from a near-bottom source and "//& + "an assumed piecewise linear profile of the buoyancy flux response to a change "//& + "in a diffusivity.", & + default=.false., do_not_log=no_BBL) + !/ Options related to Langmuir turbulence call get_param(param_file, mdl, "USE_LA_LI2016", use_LA_Windsea, & "A logical to use the Li et al. 2016 (submitted) formula to "//& @@ -2360,7 +3752,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "Valid values are: \n"//& "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//& "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//& - "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", & + "\t ADDITIVE - Add a Langmuir turbulence contribution to mstar to other contributions", & default=NONE_STRING, do_not_log=.true.) call get_param(param_file, mdl, "LT_ENHANCE", LT_enhance, default=-1) if (LT_ENHANCE == 0) then @@ -2384,7 +3776,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "Valid values are: \n"//& "\t NONE - Do not do any extra mixing due to Langmuir turbulence \n"//& "\t RESCALE - Use a multiplicative rescaling of mstar to account for Langmuir turbulence \n"//& - "\t ADDITIVE - Add a Langmuir turblence contribution to mstar to other contributions", & + "\t ADDITIVE - Add a Langmuir turbulence contribution to mstar to other contributions", & default=NONE_STRING) tmpstr = uppercase(tmpstr) select case (tmpstr) @@ -2404,7 +3796,7 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) "Coefficient for Langmuir enhancement of mstar", & units="nondim", default=0.447, do_not_log=(CS%LT_enhance_form==No_Langmuir)) call get_param(param_file, mdl, "LT_ENHANCE_EXP", CS%LT_ENHANCE_EXP, & - "Exponent for Langmuir enhancementt of mstar", & + "Exponent for Langmuir enhancement of mstar", & units="nondim", default=-1.33, do_not_log=(CS%LT_enhance_form==No_Langmuir)) call get_param(param_file, mdl, "LT_MOD_LAC1", CS%LaC_MLDoEK, & "Coefficient for modification of Langmuir number due to "//& @@ -2428,6 +3820,27 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) units="nondim", default=0.95, do_not_log=(CS%LT_enhance_form==No_Langmuir)) endif + !/ Options for documenting differences from parameter choices + call get_param(param_file, mdl, "EPBL_OPTIONS_DIFF", CS%options_diff, & + "If positive, this is a coded integer indicating a pair of settings whose "//& + "differences are diagnosed in a passive diagnostic mode via extra calls to "//& + "ePBL_column. If this is 0 or negative no extra calls occur.", & + default=0) + if (CS%options_diff > 0) then + if (CS%options_diff == 1) then + diff_text = "EPBL_ORIGINAL_PE_CALC settings" + elseif (CS%options_diff == 2) then + diff_text = "EPBL_ANSWER_DATE settings" + elseif (CS%options_diff == 3) then + diff_text = "DIRECT_EPBL_MIXING_CALC settings" + elseif (CS%options_diff == 4) then + diff_text = "BBL DIRECT_EPBL_MIXING_CALC settings" + elseif (CS%options_diff == 5) then + diff_text = "BBL DECAY_ADJUSTED_BBL_TKE settings" + else + diff_text = "unchanged settings" + endif + endif !/ Logging parameters ! This gives a minimum decay scale that is typically much less than Angstrom. @@ -2441,32 +3854,54 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) !/ Checking output flags CS%id_ML_depth = register_diag_field('ocean_model', 'ePBL_h_ML', diag%axesT1, & - Time, 'Surface boundary layer depth', 'm', conversion=US%Z_to_m, & + Time, 'Surface boundary layer depth', units='m', conversion=US%Z_to_m, & cmor_long_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme') ! This is an alias for the same variable as ePBL_h_ML CS%id_hML_depth = register_diag_field('ocean_model', 'h_ML', diag%axesT1, & - Time, 'Surface mixed layer depth based on active turbulence', 'm', conversion=US%Z_to_m) + Time, 'Surface mixed layer depth based on active turbulence', units='m', conversion=US%Z_to_m) CS%id_TKE_wind = register_diag_field('ocean_model', 'ePBL_TKE_wind', diag%axesT1, & - Time, 'Wind-stirring source of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Wind-stirring source of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_MKE = register_diag_field('ocean_model', 'ePBL_TKE_MKE', diag%axesT1, & - Time, 'Mean kinetic energy source of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Mean kinetic energy source of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_conv = register_diag_field('ocean_model', 'ePBL_TKE_conv', diag%axesT1, & - Time, 'Convective source of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Convective source of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_forcing = register_diag_field('ocean_model', 'ePBL_TKE_forcing', diag%axesT1, & Time, 'TKE consumed by mixing surface forcing or penetrative shortwave radation '//& - 'through model layers', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + 'through model layers', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_mixing = register_diag_field('ocean_model', 'ePBL_TKE_mixing', diag%axesT1, & - Time, 'TKE consumed by mixing that deepens the mixed layer', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'TKE consumed by mixing that deepens the mixed layer', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_mech_decay = register_diag_field('ocean_model', 'ePBL_TKE_mech_decay', diag%axesT1, & - Time, 'Mechanical energy decay sink of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Mechanical energy decay sink of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_TKE_conv_decay = register_diag_field('ocean_model', 'ePBL_TKE_conv_decay', diag%axesT1, & - Time, 'Convective energy decay sink of mixed layer TKE', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + Time, 'Convective energy decay sink of mixed layer TKE', units='W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_Mixing_Length = register_diag_field('ocean_model', 'Mixing_Length', diag%axesTi, & - Time, 'Mixing Length that is used', 'm', conversion=US%Z_to_m) + Time, 'Mixing Length that is used', units='m', conversion=US%Z_to_m) CS%id_Velocity_Scale = register_diag_field('ocean_model', 'Velocity_Scale', diag%axesTi, & - Time, 'Velocity Scale that is used.', 'm s-1', conversion=US%Z_to_m*US%s_to_T) + Time, 'Velocity Scale that is used.', units='m s-1', conversion=US%Z_to_m*US%s_to_T) CS%id_MSTAR_mix = register_diag_field('ocean_model', 'MSTAR', diag%axesT1, & Time, 'Total mstar that is used.', 'nondim') + if (CS%ePBL_BBL_effic > 0.0) then + CS%id_Kd_BBL = register_diag_field('ocean_model', 'Kd_ePBL_BBL', diag%axesTi, & + Time, 'ePBL bottom boundary layer diffusivity', units='m2 s-1', conversion=GV%HZ_T_to_m2_s) + CS%id_BBL_Mix_Length = register_diag_field('ocean_model', 'BBL_Mixing_Length', diag%axesTi, & + Time, 'ePBL bottom boundary layer mixing length', units='m', conversion=US%Z_to_m) + CS%id_BBL_Vel_Scale = register_diag_field('ocean_model', 'BBL_Velocity_Scale', diag%axesTi, & + Time, 'ePBL bottom boundary layer velocity scale', units='m s-1', conversion=US%Z_to_m*US%s_to_T) + CS%id_BBL_depth = register_diag_field('ocean_model', 'h_BBL', diag%axesT1, & + Time, 'Bottom boundary layer depth based on active turbulence', units='m', conversion=US%Z_to_m) + CS%id_ustar_BBL = register_diag_field('ocean_model', 'ePBL_ustar_BBL', diag%axesT1, & + Time, 'The bottom boundary layer friction velocity', units='m s-1', conversion=GV%H_to_m*US%s_to_T) + CS%id_BBL_decay_scale = register_diag_field('ocean_model', 'BBL_decay_scale', diag%axesT1, & + Time, 'The bottom boundary layer TKE decay lengthscale', units='m', conversion=GV%H_to_m) + CS%id_TKE_BBL = register_diag_field('ocean_model', 'ePBL_BBL_TKE', diag%axesT1, & + Time, 'The source of TKE for the bottom boundary layer', units='W m-2', conversion=US%RZ3_T3_to_W_m2) + CS%id_TKE_BBL_mixing = register_diag_field('ocean_model', 'ePBL_BBL_TKE_mixing', diag%axesT1, & + Time, 'TKE consumed by mixing that thickens the bottom boundary layer', & + units='W m-2', conversion=US%RZ3_T3_to_W_m2) + CS%id_TKE_BBL_decay = register_diag_field('ocean_model', 'ePBL_BBL_TKE_decay', diag%axesT1, & + Time, 'Energy decay sink of mixed layer TKE in the bottom boundary layer', & + units='W m-2', conversion=US%RZ3_T3_to_W_m2) + endif if (CS%use_LT) then CS%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, & Time, 'Langmuir number.', 'nondim') @@ -2476,37 +3911,33 @@ subroutine energetic_PBL_init(Time, G, GV, US, param_file, diag, CS) Time, 'Increase in mstar due to Langmuir Turbulence.', 'nondim') endif - call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, & - "If true, temperature and salinity are used as state "//& - "variables.", default=.true.) + if (CS%options_diff > 0) then + CS%id_opt_diff_Kd_ePBL = register_diag_field('ocean_model', 'ePBL_opt_diff_Kd_ePBL', diag%axesTi, & + Time, 'Change in ePBL diapycnal diffusivity at interfaces due to '//trim(diff_text), & + units='m2 s-1', conversion=GV%HZ_T_to_m2_s) + CS%id_opt_maxdiff_Kd_ePBL = register_diag_field('ocean_model', 'ePBL_opt_maxdiff_Kd_ePBL', diag%axesT1, & + Time, 'Column maximum change in ePBL diapycnal diffusivity at interfaces due to '//trim(diff_text), & + units='m2 s-1', conversion=GV%HZ_T_to_m2_s) + CS%id_opt_diff_hML_depth = register_diag_field('ocean_model', 'ePBL_opt_diff_h_ML', diag%axesT1, Time, & + 'Change in surface or bottom boundary layer depth based on active turbulence due to '//trim(diff_text), & + units='m', conversion=US%Z_to_m) + endif if (report_avg_its) then CS%sum_its(1) = real_to_EFP(0.0) ; CS%sum_its(2) = real_to_EFP(0.0) + CS%sum_its_BBL(1) = real_to_EFP(0.0) ; CS%sum_its_BBL(2) = real_to_EFP(0.0) endif - if (max(CS%id_TKE_wind, CS%id_TKE_MKE, CS%id_TKE_conv, & - CS%id_TKE_mixing, CS%id_TKE_mech_decay, CS%id_TKE_forcing, & - CS%id_TKE_conv_decay) > 0) then - call safe_alloc_alloc(CS%diag_TKE_wind, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_MKE, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_conv, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_forcing, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_mixing, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_mech_decay, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%diag_TKE_conv_decay, isd, ied, jsd, jed) - - CS%TKE_diagnostics = .true. + CS%TKE_diagnostics = (max(CS%id_TKE_wind, CS%id_TKE_MKE, CS%id_TKE_conv, & + CS%id_TKE_mixing, CS%id_TKE_mech_decay, CS%id_TKE_forcing, & + CS%id_TKE_conv_decay) > 0) + if (CS%ePBL_BBL_effic > 0.0) then + CS%TKE_diagnostics = CS%TKE_diagnostics .or. & + (max(CS%id_TKE_BBL, CS%id_TKE_BBL_mixing, CS%id_TKE_BBL_decay) > 0) endif - if (CS%id_Velocity_Scale>0) call safe_alloc_alloc(CS%Velocity_Scale, isd, ied, jsd, jed, GV%ke+1) - if (CS%id_Mixing_Length>0) call safe_alloc_alloc(CS%Mixing_Length, isd, ied, jsd, jed, GV%ke+1) call safe_alloc_alloc(CS%ML_depth, isd, ied, jsd, jed) - if (max(CS%id_mstar_mix, CS%id_LA, CS%id_LA_mod, CS%id_MSTAR_LT ) >0) then - call safe_alloc_alloc(CS%Mstar_mix, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%LA, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%LA_MOD, isd, ied, jsd, jed) - call safe_alloc_alloc(CS%MSTAR_LT, isd, ied, jsd, jed) - endif + call safe_alloc_alloc(CS%BBL_depth, isd, ied, jsd, jed) end subroutine energetic_PBL_init @@ -2518,26 +3949,20 @@ subroutine energetic_PBL_end(CS) real :: avg_its ! The averaged number of iterations used by ePBL [nondim] if (allocated(CS%ML_depth)) deallocate(CS%ML_depth) - if (allocated(CS%LA)) deallocate(CS%LA) - if (allocated(CS%LA_MOD)) deallocate(CS%LA_MOD) - if (allocated(CS%MSTAR_MIX)) deallocate(CS%MSTAR_MIX) - if (allocated(CS%MSTAR_LT)) deallocate(CS%MSTAR_LT) - if (allocated(CS%diag_TKE_wind)) deallocate(CS%diag_TKE_wind) - if (allocated(CS%diag_TKE_MKE)) deallocate(CS%diag_TKE_MKE) - if (allocated(CS%diag_TKE_conv)) deallocate(CS%diag_TKE_conv) - if (allocated(CS%diag_TKE_forcing)) deallocate(CS%diag_TKE_forcing) - if (allocated(CS%diag_TKE_mixing)) deallocate(CS%diag_TKE_mixing) - if (allocated(CS%diag_TKE_mech_decay)) deallocate(CS%diag_TKE_mech_decay) - if (allocated(CS%diag_TKE_conv_decay)) deallocate(CS%diag_TKE_conv_decay) - if (allocated(CS%Mixing_Length)) deallocate(CS%Mixing_Length) - if (allocated(CS%Velocity_Scale)) deallocate(CS%Velocity_Scale) + if (allocated(CS%BBL_depth)) deallocate(CS%BBL_depth) if (report_avg_its) then call EFP_sum_across_PEs(CS%sum_its, 2) - avg_its = EFP_to_real(CS%sum_its(1)) / EFP_to_real(CS%sum_its(2)) write (mesg,*) "Average ePBL iterations = ", avg_its call MOM_mesg(mesg) + + if (CS%ePBL_BBL_effic > 0.0) then + call EFP_sum_across_PEs(CS%sum_its_BBL, 2) + avg_its = EFP_to_real(CS%sum_its_BBL(1)) / EFP_to_real(CS%sum_its_BBL(2)) + write (mesg,*) "Average ePBL BBL iterations = ", avg_its + call MOM_mesg(mesg) + endif endif end subroutine energetic_PBL_end diff --git a/src/parameterizations/vertical/MOM_set_diffusivity.F90 b/src/parameterizations/vertical/MOM_set_diffusivity.F90 index d15a364267..0e4213c0a6 100644 --- a/src/parameterizations/vertical/MOM_set_diffusivity.F90 +++ b/src/parameterizations/vertical/MOM_set_diffusivity.F90 @@ -182,6 +182,7 @@ module MOM_set_diffusivity integer :: id_Kd_quad = -1, id_Kd_itidal = -1, id_Kd_Froude = -1, id_Kd_slope = -1 integer :: id_prof_leak = -1, id_prof_quad = -1, id_prof_itidal= -1 integer :: id_prof_Froude= -1, id_prof_slope = -1, id_bbl_thick = -1, id_kbbl = -1 + integer :: id_Kd_Work_added = -1 !>@} end type set_diffusivity_CS @@ -192,7 +193,8 @@ module MOM_set_diffusivity N2_3d => NULL(), & !< squared buoyancy frequency at interfaces [T-2 ~> s-2] Kd_user => NULL(), & !< user-added diffusivity at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kd_BBL => NULL(), & !< BBL diffusivity at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] - Kd_work => NULL(), & !< layer integrated work by diapycnal mixing [R Z3 T-3 ~> W m-2] + Kd_Work => NULL(), & !< layer integrated work by diapycnal mixing [R Z3 T-3 ~> W m-2] + Kd_Work_added => NULL(), & !< layer integrated work by added mixing [R Z3 T-3 ~> W m-2] maxTKE => NULL(), & !< energy required to entrain to h_max [H Z2 T-3 ~> m3 s-3 or W m-2] Kd_bkgnd => NULL(), & !< Background diffusivity at interfaces [H Z T-1 ~> m2 s-1 or kg m-1 s-1] Kv_bkgnd => NULL(), & !< Viscosity from background diffusivity at interfaces [H Z T-1 ~> m2 s-1 or Pa s] @@ -359,7 +361,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i if (CS%id_N2 > 0) allocate(dd%N2_3d(isd:ied,jsd:jed,nz+1), source=0.0) if (CS%id_Kd_user > 0) allocate(dd%Kd_user(isd:ied,jsd:jed,nz+1), source=0.0) - if (CS%id_Kd_work > 0) allocate(dd%Kd_work(isd:ied,jsd:jed,nz), source=0.0) + if (CS%id_Kd_Work > 0) allocate(dd%Kd_Work(isd:ied,jsd:jed,nz), source=0.0) + if (CS%id_Kd_Work_added > 0) allocate(dd%Kd_Work_added(isd:ied,jsd:jed,nz), source=0.0) if (CS%id_maxTKE > 0) allocate(dd%maxTKE(isd:ied,jsd:jed,nz), source=0.0) if (CS%id_TKE_to_Kd > 0) allocate(dd%TKE_to_Kd(isd:ied,jsd:jed,nz), source=0.0) if ((CS%double_diffusion) .and. (CS%id_KT_extra > 0)) & @@ -549,7 +552,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i enddo ; enddo endif - if (CS%ML_radiation .or. CS%use_tidal_mixing .or. associated(dd%Kd_work)) then + if (CS%ML_radiation .or. CS%use_tidal_mixing .or. associated(dd%Kd_Work)) then call thickness_to_dz(h, tv, dz, j, G, GV) endif @@ -662,7 +665,7 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i enddo ; enddo endif - if (associated(dd%Kd_work)) then + if (associated(dd%Kd_Work)) then do k=1,nz ; do i=is,ie dd%Kd_Work(i,j,k) = GV%H_to_RZ * Kd_lay_2d(i,k) * N2_lay(i,k) * dz(i,k) ! Watt m-2 = kg s-3 enddo ; enddo @@ -675,6 +678,12 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i enddo ; enddo endif + if (associated(dd%Kd_Work_added)) then + do k=1,nz ; do i=is,ie + dd%Kd_Work_added(i,j,k) = GV%H_to_RZ * CS%Kd_add * N2_lay(i,k) * dz(i,k) ! Watt m-2 = kg s-3 + enddo ; enddo + endif + ! Copy the 2-d slices into the 3-d array that is exported; this was done above for Kd_int. if (present(Kd_lay)) then ; do k=1,nz ; do i=is,ie Kd_lay(i,j,k) = Kd_lay_2d(i,k) @@ -754,12 +763,13 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i if (CS%use_tidal_mixing) & call post_tidal_diagnostics(G, GV, h, CS%tidal_mixing) - if (CS%id_N2 > 0) call post_data(CS%id_N2, dd%N2_3d, CS%diag) - if (CS%id_Kd_Work > 0) call post_data(CS%id_Kd_Work, dd%Kd_Work, CS%diag) - if (CS%id_maxTKE > 0) call post_data(CS%id_maxTKE, dd%maxTKE, CS%diag) - if (CS%id_TKE_to_Kd > 0) call post_data(CS%id_TKE_to_Kd, dd%TKE_to_Kd, CS%diag) + if (CS%id_N2 > 0) call post_data(CS%id_N2, dd%N2_3d, CS%diag) + if (CS%id_Kd_Work > 0) call post_data(CS%id_Kd_Work, dd%Kd_Work, CS%diag) + if (CS%id_Kd_Work_added > 0) call post_data(CS%id_Kd_Work_added, dd%Kd_Work_added, CS%diag) + if (CS%id_maxTKE > 0) call post_data(CS%id_maxTKE, dd%maxTKE, CS%diag) + if (CS%id_TKE_to_Kd > 0) call post_data(CS%id_TKE_to_Kd, dd%TKE_to_Kd, CS%diag) - if (CS%id_Kd_user > 0) call post_data(CS%id_Kd_user, dd%Kd_user, CS%diag) + if (CS%id_Kd_user > 0) call post_data(CS%id_Kd_user, dd%Kd_user, CS%diag) ! double diffusive mixing if (CS%double_diffusion) then @@ -773,7 +783,8 @@ subroutine set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, optics, visc, dt, Kd_i if (CS%id_Kd_BBL > 0) call post_data(CS%id_Kd_BBL, dd%Kd_BBL, CS%diag) if (associated(dd%N2_3d)) deallocate(dd%N2_3d) - if (associated(dd%Kd_work)) deallocate(dd%Kd_work) + if (associated(dd%Kd_Work)) deallocate(dd%Kd_Work) + if (associated(dd%Kd_Work_added)) deallocate(dd%Kd_Work_added) if (associated(dd%Kd_user)) deallocate(dd%Kd_user) if (associated(dd%maxTKE)) deallocate(dd%maxTKE) if (associated(dd%TKE_to_Kd)) deallocate(dd%TKE_to_Kd) @@ -2517,6 +2528,8 @@ subroutine set_diffusivity_init(Time, G, GV, US, param_file, diag, CS, int_tide_ if (CS%use_tidal_mixing) then CS%id_Kd_Work = register_diag_field('ocean_model', 'Kd_Work', diag%axesTL, Time, & 'Work done by Diapycnal Mixing', 'W m-2', conversion=US%RZ3_T3_to_W_m2) + CS%id_Kd_Work_added = register_diag_field('ocean_model', 'Kd_Work_added', diag%axesTL, Time, & + 'Work done by additional mixing Kd_add', 'W m-2', conversion=US%RZ3_T3_to_W_m2) CS%id_maxTKE = register_diag_field('ocean_model', 'maxTKE', diag%axesTL, Time, & 'Maximum layer TKE', 'm3 s-3', conversion=(GV%H_to_m*US%Z_to_m**2*US%s_to_T**3)) CS%id_TKE_to_Kd = register_diag_field('ocean_model', 'TKE_to_Kd', diag%axesTL, Time, & diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index b96b363e3a..aca19e86d0 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -14,9 +14,7 @@ module Idealized_hurricane ! w/ T/S initializations in CVMix_tests (which should be moved ! into the main state_initialization to their utility ! for multiple example cases). -! To do -! 1. Remove the legacy SCM_idealized_hurricane_wind_forcing code -! +! December 2024: Removed the legacy subroutine SCM_idealized_hurricane_wind_forcing use MOM_error_handler, only : MOM_error, FATAL use MOM_file_parser, only : get_param, log_version, param_file_type @@ -37,8 +35,6 @@ module Idealized_hurricane ! hurricane wind profile. public idealized_hurricane_wind_forcing !Public interface to update the idealized ! hurricane wind profile. -public SCM_idealized_hurricane_wind_forcing !Public interface to the legacy idealized - ! hurricane wind profile for SCM. !> Container for parameters describing idealized wind structure type, public :: idealized_hurricane_CS ; private @@ -656,196 +652,6 @@ subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx TY = US%L_to_Z * CS%rho_a * Cd * du10 * dV end subroutine idealized_hurricane_wind_profile -!> This subroutine is primarily needed as a legacy for reproducing answers. -!! It is included as an additional subroutine rather than padded into the previous -!! routine with flags to ease its eventual removal. Its functionality is replaced -!! with the new routines and it can be deleted when answer changes are acceptable. -subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) - type(surface), intent(in) :: sfc_state !< Surface state structure - type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces - type(time_type), intent(in) :: day !< Time in days - type(ocean_grid_type), intent(inout) :: G !< Grid structure - type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type - type(idealized_hurricane_CS), pointer :: CS !< Container for SCM parameters - ! Local variables - integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq - integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB - real :: du10 ! The magnitude of the difference between the 10 m wind and the ocean flow [L T-1 ~> m s-1] - real :: U10 ! The 10 m wind speed [L T-1 ~> m s-1] - real :: A ! The radius of the maximum winds raised to the power given by B, used in the - ! wind profile expression, in [km^B] - real :: B ! A power used in the wind profile expression [nondim] - real :: C ! A temporary variable in units of the square root of a specific volume [sqrt(m3 kg-1)] - real :: rad ! The distance from the hurricane center [L ~> m] - real :: radius10 ! The distance from the hurricane center to its edge [L ~> m] - real :: rkm ! The distance from the hurricane center, sometimes scaled to km [L ~> m] or [1000 L ~> km] - real :: f_local ! The local Coriolis parameter [T-1 ~> s-1] - real :: xx ! x-position [L ~> m] - real :: t0 ! Time at which the eye crosses the origin [T ~> s] - real :: dP ! The pressure difference across the hurricane [R L2 T-2 ~> Pa] - real :: rB ! The distance from the center raised to the power given by B, in [m^B] - ! or [km^B] if BR_Bench is true. - real :: Cd ! Air-sea drag coefficient [nondim] - real :: Uocn, Vocn ! Surface ocean velocity components [L T-1 ~> m s-1] - real :: dU, dV ! Air-sea differential motion [L T-1 ~> m s-1] - ! Wind angle variables - real :: Alph ! The wind inflow angle (positive outward) [radians] - real :: Rstr ! A function of the position normalized by the radius of maximum winds [nondim] - real :: A0 ! The axisymmetric inflow angle [degrees] - real :: A1 ! The inflow angle asymmetry [degrees] - real :: P1 ! The angle difference between the translation direction and the inflow direction [radians] - real :: Adir ! The angle of the direction from the center to a point [radians] - real :: transdir ! Translation direction [radians] - real :: V_TS, U_TS ! Components of the translation speed [L T-1 ~> m s-1] - - ! Bounds for loops and memory allocation - is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec - Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB - isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed - IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB - - ! Allocate the forcing arrays, if necessary. - call allocate_mech_forcing(G, forces, stress=.true., ustar=.true., tau_mag=.true.) - - ! Implementing Holland (1980) parameteric wind profile - !------------------------------------------------------------| - t0 = 129600.*US%s_to_T ! TC 'eye' crosses (0,0) at 36 hours | - transdir = CS%pi ! translation direction (-x) | - !------------------------------------------------------------| - dP = CS%pressure_ambient - CS%pressure_central - if (CS%answer_date < 20190101) then - C = CS%max_windspeed / sqrt( US%R_to_kg_m3*dP ) - B = C**2 * US%R_to_kg_m3*CS%rho_a * exp(1.0) - if (CS%BR_Bench) then ! rho_a reset to value used in generated wind for benchmark test - B = C**2 * 1.2 * exp(1.0) - endif - else - B = (CS%max_windspeed**2 / dP ) * CS%rho_a * exp(1.0) - endif - - if (CS%BR_Bench) then - A = (US%L_to_m*CS%rad_max_wind / 1000.)**B - else - A = (US%L_to_m*CS%rad_max_wind)**B - endif - ! f_local = f(x,y), but in the SCM it is constant - if (CS%BR_Bench) then ! (CS%SCM_mode) then - f_local = CS%f_column - else - f_local = G%CoriolisBu(is,js) - endif - - ! Calculate x position relative to hurricane center as a function of time. - xx = (t0 - time_type_to_real(day)*US%s_to_T) * CS%hurr_translation_spd * cos(transdir) - rad = sqrt((xx**2) + (CS%dy_from_center**2)) - - ! rkm - rad converted to km for Holland prof. - ! used in km due to error, correct implementation should - ! not need rkm, but to match winds w/ experiment this must - ! be maintained. Causes winds far from storm center to be a - ! couple of m/s higher than the correct Holland prof. - if (CS%BR_Bench) then - rkm = rad/1000. - rB = (US%L_to_m*rkm)**B - else - ! if not comparing to benchmark, then use correct Holland prof. - rkm = rad - rB = (US%L_to_m*rad)**B - endif - - ! Calculate U10 in the interior (inside of the hurricane edge radius), - ! while adjusting U10 to 0 outside of the ambient wind radius. - if (rad > 0.001*CS%rad_max_wind .AND. rad < CS%rad_edge*CS%rad_max_wind) then - U10 = sqrt( A*B*dP*exp(-A/rB)/(CS%rho_a*rB) + 0.25*(rkm*f_local)**2 ) - 0.5*rkm*f_local - elseif (rad > CS%rad_edge*CS%rad_max_wind .AND. rad < CS%rad_ambient*CS%rad_max_wind) then - radius10 = CS%rad_max_wind*CS%rad_edge - if (CS%BR_Bench) then - rkm = radius10/1000. - rB = (US%L_to_m*rkm)**B - else - rkm = radius10 - rB = (US%L_to_m*radius10)**B - endif - if (CS%edge_taper_bug) rad = radius10 - U10 = ( sqrt( A*B*dP*exp(-A/rB)/(CS%rho_a*rB) + 0.25*(rkm*f_local)**2 ) - 0.5*rkm*f_local) & - * (CS%rad_ambient - rad/CS%rad_max_wind)/(CS%rad_ambient - CS%rad_edge) - else - U10 = 0. - endif - Adir = atan2(CS%dy_from_center,xx) - - ! Wind angle model following Zhang and Ulhorn (2012) - ! ALPH is inflow angle positive outward. - RSTR = min(CS%rad_edge, rad / CS%rad_max_wind) - A0 = CS%A0_Rnorm*RSTR + CS%A0_speed*CS%max_windspeed + CS%A0_0 - A1 = -A0*(CS%A1_Rnorm*RSTR + CS%A1_speed*CS%hurr_translation_spd + CS%A1_0) - P1 = (CS%P1_Rnorm*RSTR + CS%P1_speed*CS%hurr_translation_spd + CS%P1_0) * CS%pi/180. - ALPH = A0 - A1*cos( (TRANSDIR - ADIR ) - P1) - if (rad > CS%rad_edge*CS%rad_max_wind .AND. rad < CS%rad_ambient*CS%rad_max_wind) then - ALPH = ALPH* (CS%rad_ambient - rad/CS%rad_max_wind) / (CS%rad_ambient - CS%rad_edge) - elseif (rad > CS%rad_ambient*CS%rad_max_wind) then - ALPH = 0.0 - endif - ALPH = ALPH * CS%Deg2Rad - - ! Prepare for wind calculation - ! X_TS is component of translation speed added to wind vector - ! due to background steering wind. - U_TS = CS%hurr_translation_spd*0.5*cos(transdir) - V_TS = CS%hurr_translation_spd*0.5*sin(transdir) - - ! Set the surface wind stresses, in [R L Z T-2 ~> Pa]. A positive taux - ! accelerates the ocean to the (pseudo-)east. - ! The i-loop extends to is-1 so that taux can be used later in the - ! calculation of ustar - otherwise the lower bound would be Isq. - do j=js,je ; do I=is-1,Ieq - ! Turn off surface current for stress calculation to be - ! consistent with test case. - Uocn = 0. ! sfc_state%u(I,j) - Vocn = 0. ! 0.25*( (sfc_state%v(i,J) + sfc_state%v(i+1,J-1)) + & - ! (sfc_state%v(i+1,J) + sfc_state%v(i,J-1)) ) - ! Wind vector calculated from location/direction (sin/cos flipped b/c - ! cyclonic wind is 90 deg. phase shifted from position angle). - dU = U10*sin(Adir - CS%pi - Alph) - Uocn + U_TS - dV = U10*cos(Adir - Alph) - Vocn + V_TS - !/----------------------------------------------------| - ! Add a simple drag coefficient as a function of U10 | - !/----------------------------------------------------| - du10 = sqrt((du**2) + (dv**2)) - Cd = simple_wind_scaled_Cd(u10, du10, CS) - - forces%taux(I,j) = CS%rho_a * US%L_to_Z * G%mask2dCu(I,j) * Cd*du10*dU - enddo ; enddo - - ! See notes above - do J=js-1,Jeq ; do i=is,ie - Uocn = 0. ! 0.25*( (sfc_state%u(I,j) + sfc_state%u(I-1,j+1)) + & - ! (sfc_state%u(I-1,j) + sfc_state%u(I,j+1)) ) - Vocn = 0. ! sfc_state%v(i,J) - dU = U10*sin(Adir - CS%pi - Alph) - Uocn + U_TS - dV = U10*cos(Adir-Alph) - Vocn + V_TS - du10 = sqrt((du**2) + (dv**2)) - Cd = simple_wind_scaled_Cd(u10, du10, CS) - forces%tauy(I,j) = CS%rho_a * US%L_to_Z * G%mask2dCv(I,j) * Cd*dU10*dV - enddo ; enddo - - ! Set the surface friction velocity [Z T-1 ~> m s-1]. ustar is always positive. - if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie - ! This expression can be changed if desired, but need not be. - forces%ustar(i,j) = G%mask2dT(i,j) * sqrt(US%L_to_Z * (CS%gustiness/CS%Rho0 + & - sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + & - 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))/CS%Rho0)) - enddo ; enddo ; endif - - !> Set magnitude of the wind stress [R L Z T-2 ~> Pa] - if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie - forces%tau_mag(i,j) = G%mask2dT(i,j) * (CS%gustiness + & - sqrt(0.5*((forces%taux(I-1,j)**2) + (forces%taux(I,j)**2)) + & - 0.5*((forces%tauy(i,J-1)**2) + (forces%tauy(i,J)**2)))) - enddo ; enddo ; endif - -end subroutine SCM_idealized_hurricane_wind_forcing - !> This function returns the air-sea drag coefficient using a simple function of the air-sea velocity difference. function simple_wind_scaled_Cd(u10, du10, CS) result(Cd) real, intent(in) :: U10 !< The 10 m wind speed [L T-1 ~> m s-1]