MOM_energetic_PBL.F90

!> Energetically consistent planetary boundary layer parameterization
module mom_energetic_pbl

! This file is part of MOM6. See LICENSE.md for the license.

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_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
use mom_error_handler, only : mom_error, fatal, warning, mom_mesg
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_forcing_type, only : forcing
use mom_grid, only : ocean_grid_type
use mom_string_functions, only : uppercase
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_wave_interface, only: wave_parameters_cs, get_langmuir_number

implicit none ; private

#include <MOM_memory.h>

public energetic_pbl, energetic_pbl_init, energetic_pbl_end
public energetic_pbl_get_mld

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

!> This control structure holds parameters for the MOM_energetic_PBL module
type, public :: energetic_pbl_cs ; private

  !/ Constants
  real    :: vonkar = 0.41   !< The von Karman coefficient.  This should be a runtime parameter,
                             !! but because it is set to 0.4 at runtime in KPP it might change answers.
  real    :: omega           !< The Earth's rotation rate [T-1 ~> s-1].
  real    :: omega_frac      !< When setting the decay scale for turbulence, use this fraction of
                             !! the absolute rotation rate blended with the local value of f, as
                             !! sqrt((1-omega_frac)*f^2 + omega_frac*4*omega^2) [nondim].

  !/ Convection related terms
  real    :: nstar           !< The fraction of the TKE input to the mixed layer available to drive
                             !! entrainment [nondim]. This quantity is the vertically integrated
                             !! buoyancy production minus the vertically integrated dissipation of
                             !! TKE produced by buoyancy.

  !/ Mixing Length terms
  logical :: use_mld_iteration !< If true, use the proximity to the bottom of the actively turbulent
                             !! surface boundary layer to constrain the mixing lengths.
  logical :: mld_iteration_guess !< False to default to guessing half the
                             !! ocean depth for the first iteration.
  logical :: mld_bisection   !< If true, use bisection with the iterative determination of the
                             !! 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.
  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.
                             !! 1 is law-of-the-wall at top and bottom,
                             !! 2 is more KPP like.
  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].
  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
                             !! diffusive length scale by rotation.  Making this larger decreases
                             !! the diffusivity in the planetary boundary layer.
  real    :: translay_scale  !< A scale for the mixing length in the transition layer
                             !! at the edge of the boundary layer as a fraction of the
                             !! boundary layer thickness.  The default is 0, but a
                             !! value of 0.1 might be better justified by observations.
  real    :: mld_tol         !< A tolerance for determining the boundary layer thickness when
                             !! Use_MLD_iteration is true [Z ~> m].
  real    :: min_mix_len     !< The minimum mixing length scale that will be used by ePBL [Z ~> m].
                             !! The default (0) does not set a minimum.

  !/ Velocity scale terms
  integer :: wt_scheme       !< An enumerated value indicating the method for finding the 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    :: wstar_ustar_coef !< A ratio relating the efficiency with which convectively released
                             !! energy is converted to a turbulent velocity, relative to
                             !! mechanically forced turbulent kinetic energy [nondim].
                             !! 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.

  !mstar related options
  integer :: mstar_scheme    !< An encoded integer to determine which formula is used to set mstar
  logical :: mstar_flatcap=.true. !< Set false to use asymptotic mstar cap.
  real    :: mstar_cap       !< Since MSTAR is restoring undissipated energy to mixing,
                             !! there must be a cap on how large it can be.  This
                             !! is definitely a function of latitude (Ekman limit),
                             !! but will be taken as constant for now.

  !/ vertical decay related options
  real    :: tke_decay       !< The ratio of the natural Ekman depth to the TKE decay scale [nondim].

  !/ mstar_scheme == 0
  real    :: fixed_mstar     !< Mstar is the ratio of the friction velocity cubed to the TKE available to
                             !! drive entrainment, nondimensional. This quantity is the vertically
                             !! integrated shear production minus the vertically integrated
                             !! dissipation of TKE produced by shear.  This value is used if the option
                             !! for using a fixed mstar is used.

  !/ mstar_scheme == 2
  real :: c_ek = 0.17        !< MSTAR Coefficient in rotation limit for mstar_scheme=OM4
  real :: mstar_coef = 0.3   !< MSTAR coefficient in rotation/stabilizing balance for mstar_scheme=OM4

  !/ mstar_scheme == 3
  real    :: rh18_mstar_cn1  !< MSTAR_N coefficient 1 (outter-most coefficient for fit).
                             !! Value of 0.275 in RH18.  Increasing this
                             !! coefficient increases mechanical mixing for all values of Hf/ust,
                             !! but is most effective at low values (weakly developed OSBLs).
  real    :: rh18_mstar_cn2  !< MSTAR_N coefficient 2 (coefficient outside of exponential decay).
                             !! Value of 8.0 in RH18.  Increasing this coefficient increases MSTAR
                             !! for all values of HF/ust, with a consistent affect across
                             !! a wide range of Hf/ust.
  real    :: rh18_mstar_cn3  !< MSTAR_N coefficient 3 (exponential decay coefficient). Value of
                             !! -5.0 in RH18.  Increasing this increases how quickly the value
                             !! of MSTAR decreases as Hf/ust increases.
  real    :: rh18_mstar_cs1  !< MSTAR_S coefficient for RH18 in stabilizing limit.
                             !! Value of 0.2 in RH18.
  real    :: rh18_mstar_cs2  !< MSTAR_S exponent for RH18 in stabilizing limit.
                             !! Value of 0.4 in RH18.

  !/ Coefficient for shear/convective turbulence interaction
  real :: mstar_convect_coef !< Factor to reduce mstar when statically unstable.

  !/ Langmuir turbulence related parameters
  logical :: use_lt = .false. !< Flag for using LT in Energy calculation
  integer :: lt_enhance_form !< Integer for Enhancement functional form (various options)
  real    :: lt_enhance_coef !< Coefficient in fit for Langmuir Enhancment
  real    :: lt_enhance_exp  !< Exponent in fit for Langmuir Enhancement
  real :: lac_mldoek         !< Coefficient for Langmuir number modification based on the ratio of
                             !! the mixed layer depth over the Ekman depth.
  real :: lac_mldoob_stab    !< Coefficient for Langmuir number modification based on the ratio of
                             !! the mixed layer depth over the Obukov depth with stablizing forcing.
  real :: lac_ekoob_stab     !< Coefficient for Langmuir number modification based on the ratio of
                             !! the Ekman depth over the Obukov depth with stablizing forcing.
  real :: lac_mldoob_un      !< Coefficient for Langmuir number modification based on the ratio of
                             !! the mixed layer depth over the Obukov depth with destablizing forcing.
  real :: lac_ekoob_un       !< Coefficient for Langmuir number modification based on the ratio of
                             !! the Ekman depth over the Obukov depth with destablizing forcing.
  real :: max_enhance_m = 5. !< The maximum allowed LT enhancement to the mixing.

  !/ Others
  type(time_type), pointer :: time=>null() !< A pointer to the ocean model's clock.

  logical :: tke_diagnostics = .false. !< If true, diagnostics of the TKE budget are being calculated.
  logical :: answers_2018    !< If true, use the order of arithmetic and expressions that recover the
                             !! answers from the end of 2018.  Otherwise, use updated and more robust
                             !! forms of the same expressions.
  logical :: orig_pe_calc    !< If true, the ePBL code uses the original form of the
                             !! potential energy change code.  Otherwise, it uses a newer version
                             !! that can work with successive increments to the diffusivity in
                             !! upward or downward passes.
  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 [Z ~> m].

  ! These are terms in the mixed layer TKE budget, all in [R Z3 T-3 ~> W m-2 = kg s-3].
  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]

  type(efp_type), dimension(2) :: sum_its !< 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_mstar_mix = -1, id_la_mod = -1, id_la = -1, id_mstar_lt = -1
  !>@}
end type energetic_pbl_cs

!>@{ Enumeration values for mstar_Scheme
integer, parameter :: use_fixed_mstar = 0  !< The value of mstar_scheme to use a constant mstar
integer, parameter :: mstar_from_ekman = 2 !< The value of mstar_scheme to base mstar on the ratio
                                           !! of the Ekman layer depth to the Obukhov depth
integer, parameter :: mstar_from_rh18 = 3  !< The value of mstar_scheme to base mstar of of RH18
integer, parameter :: no_langmuir = 0      !< The value of LT_ENHANCE_FORM not use Langmuir turbolence.
integer, parameter :: langmuir_rescale = 2 !< The value of LT_ENHANCE_FORM to use a multiplicative
                                           !! rescaling of mstar to account for Langmuir turbulence.
integer, parameter :: langmuir_add = 3     !< The value of LT_ENHANCE_FORM to add a contribution to
                                           !! mstar from Langmuir turblence to other contributions.
integer, parameter :: wt_from_croot_tke = 0 !< Use a constant times the cube root of remaining TKE
                                           !! to calculate the turbulent velocity.
integer, parameter :: wt_from_rh18 = 1     !< Use a scheme based on a combination of w* and v* as
                                           !! documented in Reichl & Hallberg (2018) to calculate
                                           !! the turbulent velocity.
character*(20), parameter :: constant_string = "CONSTANT"
character*(20), parameter :: om4_string = "OM4"
character*(20), parameter :: rh18_string = "REICHL_H18"
character*(20), parameter :: root_tke_string = "CUBE_ROOT_TKE"
character*(20), parameter :: none_string = "NONE"
character*(20), parameter :: rescaled_string = "RESCALE"
character*(20), parameter :: additive_string = "ADDITIVE"
!>@}

logical :: report_avg_its = .false.  !< Report the average number of ePBL iterations for debugging.

!> A type for conveniently passing around ePBL diagnostics for a column.
type, public :: epbl_column_diags ; private
  !>@{ 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
  real :: dtke_mke, dtke_mech_decay, dtke_conv_decay
  !>@}
  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]
  real, allocatable, dimension(:) :: dt_expect !< Expected temperature changes [degC]
  real, allocatable, dimension(:) :: ds_expect !< Expected salinity changes [ppt]
end type epbl_column_diags

contains

!>    This subroutine determines the diffusivities from the integrated energetics
!!  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, &
                         dSV_dT, dSV_dS, TKE_forced, buoy_flux, dt_diag, last_call, &
                         dT_expected, dS_expected, 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.
  type(unit_scale_type),   intent(in)    :: us     !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: h_3d   !< Layer thicknesses [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: u_3d   !< Zonal velocities interpolated to h points
                                                   !! [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: v_3d   !< Zonal velocities interpolated to h points
                                                   !! [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: dsv_dt !< The partial derivative of in-situ specific
                                                   !! volume with potential temperature
                                                   !! [R-1 degC-1 ~> m3 kg-1 degC-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: dsv_ds !< The partial derivative of in-situ specific
                                                   !! volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: tke_forced !< The forcing requirements to homogenize the
                                                   !! forcing that has been applied to each layer
                                                   !! [R Z3 T-2 ~> J m-2].
  type(thermo_var_ptrs),   intent(inout) :: tv     !< A structure containing pointers to any
                                                   !! available thermodynamic fields. Absent fields
                                                   !! have NULL ptrs.
  type(forcing),           intent(inout) :: fluxes !< A structure containing pointers to any
                                                   !! possible forcing fields. Unused fields have
                                                   !! NULL ptrs.
  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
                                                   !! [Z2 s-1 ~> m2 s-1].
  type(energetic_pbl_cs),  pointer       :: cs     !< The control structure returned by a previous
                                                   !! call to mixedlayer_init.
  real, dimension(SZI_(G),SZJ_(G)), &
                           intent(in)    :: buoy_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3].
  real,          optional, intent(in)    :: dt_diag   !< The diagnostic time step, which may be less
                                                   !! than dt if there are two calls to mixedlayer [T ~> s].
  logical,       optional, intent(in)    :: last_call !< If true, this is the last call to
                                                   !! mixedlayer in the current time step, so
                                                   !! diagnostics will be written. The default
                                                   !! is .true.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                 optional, intent(out)   :: dt_expected !< The values of temperature change that
                                                   !! should be expected when the returned
                                                   !! diffusivities are applied [degC].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                 optional, intent(out)   :: ds_expected !< The values of salinity change that
                                                   !! should be expected when the returned
                                                   !! diffusivities are applied [ppt].
  type(wave_parameters_cs), &
                 optional, pointer       :: waves  !< Wave CS

!    This subroutine determines the diffusivities from the integrated energetics
!  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.
!
!    For each interior interface, first discard the TKE to account for mixing
! of shortwave radiation through the next denser cell.  Next drive mixing based
! on the local? values of ustar + wstar, subject to available energy.  This
! step sets the value of Kd(K).  Any remaining energy is then subject to decay
! before being handed off to the next interface.  mech_TKE and conv_PErel are treated
! separately for the purposes of decay, but are used proportionately to drive
! mixing.
!
!   The key parameters for the mixed layer are found in the control structure.
!   To use the classic constant mstar mixied layers choose MSTAR_SCHEME=CONSTANT.
! The key parameters then include mstar, nstar, TKE_decay, and conv_decay.
! For the Oberhuber (1993) mixed layer,the values of these are:
!      mstar = 1.25,  nstar = 1, TKE_decay = 2.5, conv_decay = 0.5
! TKE_decay is 1/kappa in eq. 28 of Oberhuber (1993), while conv_decay is 1/mu.
! For a traditional Kraus-Turner mixed layer, the values are:
!      mstar = 1.25, nstar = 0.4, TKE_decay = 0.0, conv_decay = 0.0

  ! Local variables
  real, dimension(SZI_(G),SZK_(GV)) :: &
    h_2d, &         ! A 2-d slice of the layer thickness [H ~> m or kg m-2].
    t_2d, &         ! A 2-d slice of the layer temperatures [degC].
    s_2d, &         ! A 2-d slice of the layer salinities [ppt].
    tke_forced_2d, & ! A 2-d slice of TKE_forced [R Z3 T-2 ~> J m-2].
    dsv_dt_2d, &    ! A 2-d slice of dSV_dT [R-1 degC-1 ~> m3 kg-1 degC-1].
    dsv_ds_2d, &    ! A 2-d slice of dSV_dS [R-1 ppt-1 ~> m3 kg-1 ppt-1].
    u_2d, &         ! A 2-d slice of the zonal velocity [L T-1 ~> m s-1].
    v_2d            ! A 2-d slice of the meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZK_(GV)+1) :: &
    kd_2d           ! A 2-d version of the diapycnal diffusivity [Z2 T-1 ~> m2 s-1].
  real, dimension(SZK_(GV)) :: &
    h, &            ! The layer thickness [H ~> m or kg m-2].
    t0, &           ! The initial layer temperatures [degC].
    s0, &           ! The initial layer salinities [ppt].
    dsv_dt_1d, &    ! The partial derivatives of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1].
    dsv_ds_1d, &    ! The partial derivatives of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
    tke_forcing, &  ! Forcing of the TKE in the layer coming from TKE_forced [R Z3 T-2 ~> J m-2].
    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 [Z2 T-1 ~> m2 s-1].
    mixvel, &       ! A turbulent mixing veloxity [Z T-1 ~> m s-1].
    mixlen          ! A turbulent mixing length [Z ~> m].
  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 :: absf      ! The absolute value of f [T-1 ~> s-1].
  real :: u_star    ! The surface friction velocity [Z T-1 ~> m s-1].
  real :: u_star_mean ! The surface friction without gustiness [Z T-1 ~> m 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].

! The following are only used for diagnostics.
  real :: dt__diag  ! A copy of dt_diag (if present) or dt [T ~> s].
  logical :: write_diags  ! If true, write out diagnostics with this step.
  logical :: reset_diags  ! If true, zero out the accumulated diagnostics.

  logical :: debug=.false.  ! Change this hard-coded value for debugging.
  type(epbl_column_diags) :: ecd ! A container for passing around diagnostics.

  integer :: i, j, k, is, ie, js, je, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

  if (.not. associated(cs)) call mom_error(fatal, "energetic_PBL: "//&
         "Module must be initialized before it is used.")

  if (.not. associated(tv%eqn_of_state)) call mom_error(fatal, &
      "energetic_PBL: Temperature, salinity and an equation of state "//&
      "must now be used.")
  if (.NOT. associated(fluxes%ustar)) call mom_error(fatal, &
      "energetic_PBL: No surface TKE fluxes (ustar) defined in mixedlayer!")
  debug = .false. ; if (present(dt_expected) .or. present(ds_expected)) debug = .true.

  if (debug) allocate(ecd%dT_expect(nz), ecd%dS_expect(nz))

  h_neglect = gv%H_subroundoff

  dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag
  write_diags = .true. ; if (present(last_call)) write_diags = last_call


  ! Determine whether to zero out diagnostics before accumulation.
  reset_diags = .true.
  if (present(dt_diag) .and. write_diags .and. (dt__diag > dt)) &
    reset_diags = .false.  ! This is the second call to mixedlayer.

  if (reset_diags) then
    if (cs%TKE_diagnostics) then
!!OMP parallel do default(none) shared(is,ie,js,je,CS)
      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
      enddo ; enddo
    endif
  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, &
!!OMP                                  CS,G,GV,US,fluxes,debug, &
!!OMP                                  TKE_forced,dSV_dT,dSV_dS,Kd_int)
  do j=js,je
    ! Copy the thicknesses and other fields to 2-d arrays.
    do k=1,nz ; do i=is,ie
      h_2d(i,k) = h_3d(i,j,k) ; u_2d(i,k) = u_3d(i,j,k) ; v_2d(i,k) = v_3d(i,j,k)
      t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
      tke_forced_2d(i,k) = tke_forced(i,j,k)
      dsv_dt_2d(i,k) = dsv_dt(i,j,k) ; dsv_ds_2d(i,k) = dsv_ds(i,j,k)
    enddo ; enddo

    !   Determine the initial mech_TKE and conv_PErel, including the energy required
    ! to mix surface heating through the topmost cell, the energy released by mixing
    ! surface cooling & brine rejection down through the topmost cell, and
    ! homogenizing the shortwave heating within that cell.  This sets the energy
    ! and ustar and wstar available to drive mixing at the first interior
    ! interface.
    do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then

      ! Copy the thicknesses and other fields to 1-d arrays.
      do k=1,nz
        h(k) = h_2d(i,k) + gv%H_subroundoff ; u(k) = u_2d(i,k) ; v(k) = v_2d(i,k)
        t0(k) = t_2d(i,k) ; s0(k) = s_2d(i,k) ; tke_forcing(k) =  tke_forced_2d(i,k)
        dsv_dt_1d(k) = dsv_dt_2d(i,k) ; dsv_ds_1d(k) = dsv_ds_2d(i,k)
      enddo
      do k=1,nz+1 ; kd(k) = 0.0 ; enddo

      ! Make local copies of surface forcing and process them.
      u_star = fluxes%ustar(i,j)
      u_star_mean = fluxes%ustar_gustless(i,j)
      b_flux = buoy_flux(i,j)
      if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
        if (fluxes%frac_shelf_h(i,j) > 0.0) &
          u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
                   fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
      endif
      if (u_star < cs%ustar_min) u_star = cs%ustar_min
      if (cs%omega_frac >= 1.0) then
        absf = 2.0*cs%omega
      else
        absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
                     (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
        if (cs%omega_frac > 0.0) &
          absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
      endif

      ! 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)

      call epbl_column(h, u, v, t0, s0, dsv_dt_1d, dsv_ds_1d, tke_forcing, b_flux, absf, &
                       u_star, u_star_mean, dt, mld_io, kd, mixvel, mixlen, gv, &
                       us, cs, ecd, dt_diag=dt_diag, waves=waves, g=g, i=i, j=j)


      ! 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

      if (present(dt_expected)) then
        do k=1,nz ; dt_expected(i,j,k) = ecd%dT_expect(k) ; enddo
      endif
      if (present(ds_expected)) then
        do k=1,nz ; ds_expected(i,j,k) = ecd%dS_expect(k) ; enddo
      endif

      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
      endif
      ! Write to 3-D for outputing Mixing length and velocity scale.
      if (cs%id_Mixing_Length>0) then ; do k=1,nz
        cs%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)
      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
    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

      if (present(dt_expected)) then
        do k=1,nz ; dt_expected(i,j,k) = 0.0 ; enddo
      endif
      if (present(ds_expected)) then
        do k=1,nz ; ds_expected(i,j,k) = 0.0 ; enddo
      endif
    endif ; enddo ! Close of i-loop - Note unusual loop order!

    do k=1,nz+1 ; do i=is,ie ; kd_int(i,j,k) = kd_2d(i,k) ; enddo ; enddo

  enddo ! j-loop

  if (write_diags) then
    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_mech_decay > 0) &
      call post_data(cs%id_TKE_mech_decay, cs%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)
  endif

  if (debug) deallocate(ecd%dT_expect, ecd%dS_expect)

end subroutine energetic_pbl



!> This subroutine determines the diffusivities from the integrated energetics
!!  mixed layer model for a single column of water.
subroutine epbl_column(h, u, v, T0, S0, dSV_dT, dSV_dS, TKE_forcing, B_flux, absf, &
                       u_star, u_star_mean, dt, MLD_io, Kd, mixvel, mixlen, GV, US, CS, eCD, &
                       dt_diag, Waves, G, i, j)
  type(verticalGrid_type), intent(in)    :: GV     !< The ocean's vertical grid structure.
  type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
  real, dimension(SZK_(GV)), intent(in)  :: h      !< Layer thicknesses [H ~> m or kg m-2].
  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 [degC].
  real, dimension(SZK_(GV)), intent(in)  :: S0     !< The initial layer salinities [ppt].

  real, dimension(SZK_(GV)), intent(in)  :: dSV_dT !< The partial derivative of in-situ specific
                                                   !! volume with potential temperature
                                                   !! [R-1 degC-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 ppt-1 ~> m3 kg-1 ppt-1].
  real, dimension(SZK_(GV)), intent(in)  :: TKE_forcing !< The forcing requirements to homogenize the
                                                   !! forcing that has been applied to each layer
                                                   !! [R Z3 T-2 ~> J m-2].
  real,                    intent(in)    :: B_flux !< The surface buoyancy flux [Z2 T-3 ~> m2 s-3]
  real,                    intent(in)    :: absf   !< The absolute value of the Coriolis parameter [T-1 ~> s-1].
  real,                    intent(in)    :: u_star !< The surface friction velocity [Z T-1 ~> m s-1].
  real,                    intent(in)    :: u_star_mean !< The surface friction velocity without any
                                                   !! contribution from unresolved gustiness  [Z T-1 ~> m s-1].
  real,                    intent(inout) :: MLD_io !< A first guess at the mixed layer depth on input, and
                                                   !! the calculated mixed layer depth on output [Z ~> m].
  real,                    intent(in)    :: dt     !< Time increment [T ~> s].
  real, dimension(SZK_(GV)+1), &
                           intent(out)   :: Kd     !< The diagnosed diffusivities at interfaces
                                                   !! [Z2 T-1 ~> m2 s-1].
  real, dimension(SZK_(GV)+1), &
                           intent(out)   :: mixvel !< The mixing velocity scale used in Kd
                                                   !! [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),  pointer       :: CS     !< The control structure returned by a previous
                                                   !! call to mixedlayer_init.
  type(ePBL_column_diags), intent(inout) :: eCD    !< A container for passing around diagnostics.
  real,          optional, intent(in)    :: dt_diag   !< The diagnostic time step, which may be less
                                                   !! than dt if there are two calls to mixedlayer [T ~> s].
  type(wave_parameters_CS), &
                 optional, pointer       :: Waves  !< Wave CS for Langmuir turbulence
  type(ocean_grid_type), &
                 optional, intent(inout) :: G      !< The ocean's grid structure.
  integer,       optional, intent(in)    :: i      !< The i-index to work on (used for Waves)
  integer,       optional, 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
!  have already been applied.  All calculations are done implicitly, and there
!  is no stability limit on the time step.
!
!    For each interior interface, first discard the TKE to account for mixing
! of shortwave radiation through the next denser cell.  Next drive mixing based
! on the local? values of ustar + wstar, subject to available energy.  This
! step sets the value of Kd(K).  Any remaining energy is then subject to decay
! before being handed off to the next interface.  mech_TKE and conv_PErel are treated
! separately for the purposes of decay, but are used proportionately to drive
! mixing.

  ! 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].
    hb_hs           ! The distance from the bottom over the thickness of the
                    ! water column [nondim].
  real :: mech_TKE  !   The mechanically generated turbulent kinetic energy
                    ! available for mixing over a time step [R Z3 T-2 ~> J m-2].
  real :: conv_PErel ! The potential energy that has been convectively released
                    ! during this timestep [R Z3 T-2 ~> J m-2]. A portion nstar_FC
                    ! of conv_PErel is available to drive mixing.
  real :: htot      !   The total depth of the layers above an interface [H ~> m or kg m-2].
  real :: uhtot     !   The depth integrated zonal and meridional velocities in the
  real :: vhtot     ! 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 :: h_sum     ! The total thickness of the water column [H ~> m or kg m-2].

  real, dimension(SZK_(GV)) :: &
    dT_to_dColHt, & ! Partial derivative of the total column height with the temperature changes
                    ! within a layer [Z degC-1 ~> m degC-1].
    ds_to_dcolht, & ! Partial derivative of the total column height with the salinity changes
                    ! within a layer  [Z ppt-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 degC-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 ppt-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 degC-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 ppt-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 degC-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 ppt-1 ~> J m-2 ppt-1].
    c1, &           ! c1 is used by the tridiagonal solver [nondim].
    te, &           ! Estimated final values of T in the column [degC].
    se, &           ! Estimated final values of S in the column [ppt].
    dte, &          ! Running (1-way) estimates of temperature change [degC].
    dse, &          ! Running (1-way) estimates of salinity change [ppt].
    th_a, &         ! An effective temperature times a thickness in the layer above, including implicit
                    ! mixing effects with other yet higher layers [degC 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 [ppt 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 [degC 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 [ppt 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 assymptotic value at the bottom of
                    ! the boundary layer.
    kddt_h          ! The diapycnal diffusivity times a timestep divided by the
                    ! average thicknesses around a layer [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 :: 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 :: dMKE_max  ! The maximum amount of mean kinetic energy that could be
                    ! converted to turbulent kinetic energy if the velocity in
                    ! the layer below an interface were homogenized with all of
                    ! the water above the interface [R Z3 T-2 ~> J m-2].
  real :: MKE2_Hharm ! Twice the inverse of the harmonic mean of the thickness
                    ! of a layer and the thickness of the water above, used in
                    ! the MKE conversion equation [H-1 ~> m-1 or m2 kg-1].

  real :: dt_h      ! The timestep divided by the averages of the thicknesses around
                    ! a layer, times a thickness conversion factor [H T m-2 ~> s m-1 or kg s m-4].
  real :: h_bot     ! The distance from the bottom [H ~> m or kg m-2].
  real :: h_rsum    ! The running sum of h from the top [Z ~> m].
  real :: I_hs      ! The inverse of h_sum [H-1 ~> m-1 or m2 kg-1].
  real :: I_MLD     ! The inverse of the current value of MLD [Z-1 ~> m-1].
  real :: h_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 h_tt_min [H ~> m or kg m-2].
  real :: h_tt_min  ! A surface roughness length [H ~> m or kg m-2].

  real :: C1_3      ! = 1/3.
  real :: I_dtrho   ! 1.0 / (dt * Rho0) times conversion factors in [m3 Z-3 R-1 T2 s-3 ~> m3 kg-1 s-1].
                    ! This is used convert TKE back into ustar^3.
  real :: vstar     ! An in-situ turbulent velocity [Z T-1 ~> m s-1].
  real :: mstar_total ! The value of mstar used in ePBL [nondim]
  real :: mstar_LT  ! An addition to mstar due to Langmuir turbulence [nondim] (output for diagnostic)
  real :: MLD_output ! The mixed layer depth output from this routine [Z ~> m].
  real :: LA        ! The value of the Langmuir number [nondim]
  real :: LAmod     ! The modified Langmuir number by convection [nondim]
  real :: hbs_here  ! The local minimum of hb_hs and MixLen_shape, times a
                    ! conversion factor from H to Z [Z H-1 ~> 1 or m3 kg-1].
  real :: nstar_FC  ! The fraction of conv_PErel that can be converted to mixing [nondim].
  real :: TKE_reduc ! The fraction by which TKE and other energy fields are
                    ! reduced to support mixing [nondim]. between 0 and 1.
  real :: tot_TKE   ! The total TKE available to support mixing at interface K [R Z3 T-2 ~> J m-2].
  real :: TKE_here  ! The total TKE at this point in the algorithm [R Z3 T-2 ~> J m-2].
  real :: dT_km1_t2 ! A diffusivity-independent term related to the temperature
                    ! change in the layer above the interface [degC].
  real :: dS_km1_t2 ! A diffusivity-independent term related to the salinity
                    ! change in the layer above the interface [ppt].
  real :: dTe_term  ! A diffusivity-independent term related to the temperature
                    ! change in the layer below the interface [degC H ~> degC m or degC kg m-2].
  real :: dSe_term  ! A diffusivity-independent term related to the salinity
                    ! change in the layer above the interface [ppt H ~> ppt m or ppt kg m-2].
  real :: dTe_t2    ! A part of dTe_term [degC H ~> degC m or degC kg m-2].
  real :: dSe_t2    ! A part of dSe_term [ppt H ~> ppt m or ppt kg m-2].
  real :: dPE_conv  ! The convective change in column potential energy [R Z3 T-2 ~> J m-2].
  real :: MKE_src   ! The mean kinetic energy source of TKE due to Kddt_h(K) [R Z3 T-2 ~> J m-2].
  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 [Z2 T-1 ~> m2 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 :: 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].
  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 :: vstar_unit_scale ! A unit converion factor for turbulent velocities [Z T-1 s m-1 ~> 1]
  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 :: sfc_connected   ! If true the ocean is actively turbulent from the present
                    ! interface all the way up to the surface.
  logical :: sfc_disconnect ! If true, any turbulence has become disconnected
                    ! from the surface.

! The following are only used for diagnostics.
  real :: dt__diag  ! A copy of dt_diag (if present) or dt [T ~> s].
  real :: I_dtdiag  !  = 1.0 / dt__diag [T-1 ~> s-1].

  !----------------------------------------------------------------------
  !/BGR added Aug24,2016 for adding iteration to get boundary layer depth
  !    - needed to compute new mixing length.
  real :: MLD_guess, MLD_found ! Mixing Layer depth guessed/found for iteration [Z ~> m].
  real :: min_MLD   ! Iteration bounds [Z ~> m], which are adjusted at each step
  real :: max_MLD   !  - These are initialized based on surface/bottom
                    !  1. The iteration guesses a value (possibly from prev step or neighbor).
                    !  2. The iteration checks if value is converged, too shallow, or too deep.
                    !  3. Based on result adjusts the Max/Min and searches through the water column.
                    !  - If using an accurate guess the iteration is very quick (e.g. if MLD doesn't
                    !    change over timestep).  Otherwise it takes 5-10 passes, but has a high
                    !    convergence rate.  Other iteration may be tried, but this method seems to
                    !    fail very rarely and the added cost is likely not significant.
                    !    Additionally, when it fails to converge it does so in a reasonable
                    !    manner giving a usable guess. When it does fail, it is due to convection
                    !    within the boundary layer.  Likely, a new method e.g. surface_disconnect,
                    !    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 :: FIRST_OBL     ! Flag for computing "found" Mixing layer depth
  logical :: OBL_converged ! Flag for convergence of MLD
  integer :: OBL_it        ! Iteration counter

  real :: Surface_Scale ! Surface decay scale for vstar
  logical :: calc_dT_expect ! If true calculate the expected changes in temperature and salinity.
  logical :: calc_Te        ! If true calculate the expected final temperature and salinity values.
  logical :: debug=.false.  ! Change this hard-coded value for debugging.

  !  The following arrays are used only for debugging purposes.
  real :: dPE_debug, mixing_debug, taux2, tauy2
  real, dimension(20) :: TKE_left_itt, PE_chg_itt, Kddt_h_itt, dPEa_dKd_itt, MKE_src_itt
  real, dimension(SZK_(GV)) :: mech_TKE_k, conv_PErel_k, nstar_k
  integer, dimension(SZK_(GV)) :: num_itts

  integer :: k, nz, itt, max_itt

  nz = gv%ke

  if (.not. associated(cs)) call mom_error(fatal, "energetic_PBL: "//&
         "Module must be initialized before it is used.")

  calc_dt_expect = debug ; if (allocated(ecd%dT_expect) .or. allocated(ecd%dS_expect)) calc_dt_expect = .true.
  calc_te = (calc_dt_expect .or. (.not.cs%orig_PE_calc))

  h_neglect = gv%H_subroundoff

  c1_3 = 1.0 / 3.0
  dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag
  i_dtdiag = 1.0 / dt__diag
  max_itt = 20

  h_tt_min = 0.0
  i_dtrho = 0.0 ; if (dt*gv%Rho0 > 0.0) i_dtrho = (us%Z_to_m**3*us%s_to_T**3) / (dt*gv%Rho0)
  vstar_unit_scale = us%m_to_Z * us%T_to_s

  mld_guess = mld_io

!   Determine the initial mech_TKE and conv_PErel, including the energy required
! to mix surface heating through the topmost cell, the energy released by mixing
! surface cooling & brine rejection down through the topmost cell, and
! homogenizing the shortwave heating within that cell.  This sets the energy
! and ustar and wstar available to drive mixing at the first interior
! interface.

  do k=1,nz+1 ; kd(k) = 0.0 ; enddo

  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

  ! Determine the total thickness (h_sum) and the fractional distance from the bottom (hb_hs).
  h_sum = h_neglect ; do k=1,nz ; h_sum = h_sum + h(k) ; enddo
  i_hs = 0.0 ; if (h_sum > 0.0) i_hs = 1.0 / h_sum
  h_bot = 0.0
  hb_hs(nz+1) = 0.0
  do k=nz,1,-1
    h_bot = h_bot + h(k)
    hb_hs(k) = h_bot * i_hs
  enddo

  mld_output = h(1)*gv%H_to_Z

  !/The following lines are for the iteration over MLD
  ! max_MLD will initialized as ocean bottom depth
  max_mld = 0.0 ; do k=1,nz ; max_mld = max_mld + h(k)*gv%H_to_Z ; enddo
  ! min_MLD will be initialized to 0.
  min_mld = 0.0
  ! Set values of the wrong signs to indicate that these changes are not based on valid estimates
  dmld_min = -1.0*us%m_to_Z ; dmld_max = 1.0*us%m_to_Z

  ! If no first guess is provided for MLD, try the middle of the water column
  if (mld_guess <= min_mld) mld_guess = 0.5 * (min_mld + max_mld)

  ! 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.

      if (debug) then ; mech_tke_k(:) = 0.0 ; conv_perel_k(:) = 0.0 ; endif

      ! Reset ML_depth
      mld_output = h(1)*gv%H_to_Z
      sfc_connected = .true.

      !/ Here we get MStar, which is the ratio of convective TKE driven mixing to UStar**3
      if (cs%Use_LT) then
        call get_langmuir_number(la, g, gv, us, abs(mld_guess), u_star_mean, i, j, &
                                 h=h, u_h=u, v_h=v, waves=waves)
        call find_mstar(cs, us, b_flux, u_star, u_star_mean, mld_guess, absf, &
                        mstar_total, langmuir_number=la, convect_langmuir_number=lamod,&
                        mstar_lt=mstar_lt)
      else
        call find_mstar(cs, us, b_flux, u_star, u_star_mean, mld_guess, absf, mstar_total)
      endif

      !/ Apply MStar to get mech_TKE
      if ((cs%answers_2018) .and. (cs%mstar_scheme==use_fixed_mstar)) then
        mech_tke = (dt*mstar_total*gv%Rho0) * u_star**3
      else
        mech_tke = mstar_total * (dt*gv%Rho0* u_star**3)
      endif

      if (cs%TKE_diagnostics) then
        ecd%dTKE_conv = 0.0 ; ecd%dTKE_mixing = 0.0
        ecd%dTKE_MKE = 0.0 ; ecd%dTKE_mech_decay = 0.0 ; ecd%dTKE_conv_decay = 0.0

        ecd%dTKE_wind = mech_tke * i_dtdiag
        if (tke_forcing(1) <= 0.0) then
          ecd%dTKE_forcing = max(-mech_tke, tke_forcing(1)) * i_dtdiag
          ! eCD%dTKE_unbalanced = min(0.0, TKE_forcing(1) + mech_TKE) * I_dtdiag
        else
          ecd%dTKE_forcing = cs%nstar*tke_forcing(1) * i_dtdiag
          ! eCD%dTKE_unbalanced = 0.0
        endif
      endif

      if (tke_forcing(1) <= 0.0) then
        mech_tke = mech_tke + tke_forcing(1)
        if (mech_tke < 0.0) mech_tke = 0.0
        conv_perel = 0.0
      else
        conv_perel = tke_forcing(1)
      endif


      ! Store in 1D arrays for output.
      do k=1,nz+1 ; mixvel(k) = 0.0 ; mixlen(k) = 0.0 ; enddo

      ! Determine the mixing shape function MixLen_shape.
      if ((.not.cs%Use_MLD_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 (mld_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 MLD, with a quadratic
        ! expression that follows KPP.
        i_mld = 1.0 / mld_guess
        h_rsum = 0.0
        mixlen_shape(1) = 1.0
        do k=2,nz+1
          h_rsum = h_rsum + h(k-1)*gv%H_to_Z
          if (cs%MixLenExponent==2.0) then
            mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
                 (max(0.0, (mld_guess - h_rsum)*i_mld) )**2 ! CS%MixLenExponent
          else
            mixlen_shape(k) = cs%transLay_scale + (1.0 - cs%transLay_scale) * &
                 (max(0.0, (mld_guess - h_rsum)*i_mld) )**cs%MixLenExponent
          endif
        enddo
      endif

      kd(1) = 0.0 ; kddt_h(1) = 0.0
      hp_a = 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)

      htot = h(1) ; uhtot = u(1)*h(1) ; vhtot = v(1)*h(1)

      if (debug) then
        mech_tke_k(1) = mech_tke ; conv_perel_k(1) = conv_perel
        nstar_k(:) = 0.0 ; nstar_k(1) = cs%nstar ; num_itts(:) = -1
      endif

      do k=2,nz
        ! Apply dissipation to the TKE, here applied as an exponential decay
        ! due to 3-d turbulent energy being lost to inefficient rotational modes.

        !   There should be several different "flavors" of TKE that decay at
        ! different rates.  The following form is often used for mechanical
        ! stirring from the surface, perhaps due to breaking surface gravity
        ! waves and wind-driven turbulence.
        idecay_len_tke = (cs%TKE_decay * absf / u_star) * gv%H_to_Z
        exp_kh = 1.0
        if (idecay_len_tke > 0.0) exp_kh = exp(-h(k-1)*idecay_len_tke)
        if (cs%TKE_diagnostics) &
          ecd%dTKE_mech_decay = ecd%dTKE_mech_decay + (exp_kh-1.0) * mech_tke * i_dtdiag
        mech_tke = mech_tke * exp_kh

        !   Accumulate any convectively released potential energy to contribute
        ! to wstar and to drive penetrating convection.
        if (tke_forcing(k) > 0.0) then
          conv_perel = conv_perel + tke_forcing(k)
          if (cs%TKE_diagnostics) &
            ecd%dTKE_forcing = ecd%dTKE_forcing + cs%nstar*tke_forcing(k) * i_dtdiag
        endif

        if (debug) then
          mech_tke_k(k) = mech_tke ; conv_perel_k(k) = conv_perel
        endif

        !  Determine the total energy
        nstar_fc = cs%nstar
        if (cs%nstar * conv_perel > 0.0) then
          ! Here nstar is a function of the natural Rossby number 0.2/(1+0.2/Ro), based
          ! on a curve fit from the data of Wang (GRL, 2003).
          ! Note:         Ro = 1.0 / sqrt(0.5 * dt * Rho0 * (absf*htot)**3 / conv_PErel)
          nstar_fc = cs%nstar * conv_perel / (conv_perel + 0.2 * &
                     sqrt(0.5 * dt * gv%Rho0 * (absf*(htot*gv%H_to_Z))**3 * conv_perel))
        endif

        if (debug) nstar_k(k) = nstar_fc

        tot_tke = mech_tke + nstar_fc * conv_perel

        !   For each interior interface, first discard the TKE to account for
        ! mixing of shortwave radiation through the next denser cell.
        if (tke_forcing(k) < 0.0) then
          if (tke_forcing(k) + tot_tke < 0.0) then
            ! The shortwave requirements deplete all the energy in this layer.
            if (cs%TKE_diagnostics) then
              ecd%dTKE_mixing = ecd%dTKE_mixing + tot_tke * i_dtdiag
              ecd%dTKE_forcing = ecd%dTKE_forcing - tot_tke * i_dtdiag
              ! eCD%dTKE_unbalanced = eCD%dTKE_unbalanced + (TKE_forcing(k) + tot_TKE) * I_dtdiag
              ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + (cs%nstar-nstar_fc) * conv_perel * i_dtdiag
            endif
            tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
          else
            ! Reduce the mechanical and convective TKE proportionately.
            tke_reduc = (tot_tke + tke_forcing(k)) / tot_tke
            if (cs%TKE_diagnostics) then
              ecd%dTKE_mixing = ecd%dTKE_mixing - tke_forcing(k) * i_dtdiag
              ecd%dTKE_forcing = ecd%dTKE_forcing + tke_forcing(k) * 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 = tke_reduc*tot_tke   ! = tot_TKE + TKE_forcing(k)
            mech_tke = tke_reduc*mech_tke
            conv_perel = tke_reduc*conv_perel
          endif
        endif

        ! Precalculate some temporary expressions that are independent of Kddt_h(K).
        if (cs%orig_PE_calc) then
          if (k==2) then
            dte_t2 = 0.0 ; dse_t2 = 0.0
          else
            dte_t2 = kddt_h(k-1) * ((t0(k-2) - t0(k-1)) + dte(k-2))
            dse_t2 = kddt_h(k-1) * ((s0(k-2) - s0(k-1)) + dse(k-2))
          endif
        endif
        dt_h = (gv%Z_to_H**2*dt) / max(0.5*(h(k-1)+h(k)), 1e-15*h_sum)

        !   This tests whether the layers above and below this interface are in
        ! a convetively 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-weigted 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_tke + conv_perel) <= 0.0 .and. convectively_stable) then
          ! Energy is already exhausted, so set Kd = 0 and cycle or exit?
          tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
          kd(k) = 0.0 ; kddt_h(k) = 0.0
          sfc_disconnect = .true.
          ! if (.not.debug) exit

         !   The estimated properties for layer k-1 can be calculated, using
         ! greatly simplified expressions when Kddt_h = 0.  This enables the
         ! 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
          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)
          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)
          ds_to_dcolht_a(k) = ds_to_dcolht(k)

        else ! tot_TKE > 0.0 or this is a potentially convectively unstable profile.
          sfc_disconnect = .false.

          ! Precalculate some more temporary expressions that are independent of
          ! Kddt_h(K).
          if (cs%orig_PE_calc) then
            if (k==2) then
              dt_km1_t2 = (t0(k)-t0(k-1))
              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))
              ds_km1_t2 = (s0(k)-s0(k-1)) - &
                    (kddt_h(k-1) / hp_a) * ((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))
          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)
            else
              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
            th_b(k) = h(k) * t0(k) ; sh_b(k) = h(k) * s0(k)
          endif

          !   Using Pr=1 and the diffusivity at the bottom 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.

          if ((cs%MKE_to_TKE_effic > 0.0) .and. (htot*h(k) > 0.0)) then
            ! This is the energy that would be available from homogenizing the
            ! velocities between layer k and the layers above.
            dmke_max = (us%L_to_Z**2*gv%H_to_RZ * cs%MKE_to_TKE_effic) * 0.5 * &
                (h(k) / ((htot + h(k))*htot)) * &
                ((uhtot-u(k)*htot)**2 + (vhtot-v(k)*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) + 2.0*h_neglect) / &
                         ((htot+h_neglect) * (h(k)+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.  mech_TKE might be negative due to
          ! contributions from TKE_forced.
          h_tt = htot + h_tt_min
          tke_here = mech_tke + cs%wstar_ustar_coef*conv_perel
          if (tke_here > 0.0) then
            if (cs%wT_scheme==wt_from_croot_tke) then
              vstar = cs%vstar_scale_fac * vstar_unit_scale * (i_dtrho*tke_here)**c1_3
            elseif (cs%wT_scheme==wt_from_rh18) then
              surface_scale = max(0.05, 1.0 - htot/mld_guess)
              vstar = cs%vstar_scale_fac * surface_scale * (cs%vstar_surf_fac*u_star + &
                        vstar_unit_scale * (cs%wstar_ustar_coef*conv_perel*i_dtrho)**c1_3)
            endif
            hbs_here = gv%H_to_Z * min(hb_hs(k), mixlen_shape(k))
            mixlen(k) = max(cs%min_mix_len, ((h_tt*hbs_here)*vstar) / &
                ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar))
            !Note setting Kd_guess0 to vstar * CS%vonKar * mixlen(K) here will
            ! change the answers.  Therefore, skipping that.
            if (.not.cs%Use_MLD_iteration) then
              kd_guess0 = vstar * cs%vonKar * ((h_tt*hbs_here)*vstar) / &
                ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar)
            else
              kd_guess0 = vstar * cs%vonKar * mixlen(k)
            endif
          else
            vstar = 0.0 ; kd_guess0 = 0.0
          endif
          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, &
                     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 )
          else
            call find_pe_chg(0.0, kddt_h_g0, hp_a, 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 )
          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
            ! 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?
              tke_here = mech_tke + cs%wstar_ustar_coef*(conv_perel-pe_chg_max)
              if (tke_here > 0.0) then
                if (cs%wT_scheme==wt_from_croot_tke) then
                  vstar = cs%vstar_scale_fac * vstar_unit_scale * (i_dtrho*tke_here)**c1_3
                elseif (cs%wT_scheme==wt_from_rh18) then
                  surface_scale = max(0.05, 1. - htot/mld_guess)
                  vstar = cs%vstar_scale_fac * surface_scale * (cs%vstar_surf_fac*u_star + &
                                  vstar_unit_scale * (cs%wstar_ustar_coef*conv_perel*i_dtrho)**c1_3)
                endif
                hbs_here = gv%H_to_Z * min(hb_hs(k), mixlen_shape(k))
                mixlen(k) = max(cs%min_mix_len, ((h_tt*hbs_here)*vstar) / &
                    ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar))
                if (.not.cs%Use_MLD_iteration) then
                ! Note again (as prev) that using mixlen here
                !  instead of redoing the computation will change answers...
                  kd(k) = vstar * cs%vonKar *  ((h_tt*hbs_here)*vstar) / &
                        ((cs%Ekman_scale_coef * absf) * (h_tt*hbs_here) + vstar)
                else
                  kd(k) = vstar * cs%vonKar * mixlen(k)
                endif
              else
                vstar = 0.0 ; kd(k) = 0.0
              endif
              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, &
                         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), &
                         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)
              endif
              ! Should this be iterated to convergence for Kd?
              if (dpe_conv > 0.0) then
                kd(k) = kd_guess0 ; dpe_conv = pe_chg_g0
              else
                mke_src = dmke_max*(1.0 - exp(-(kd(k)*dt_h) * mke2_hharm))
              endif
            else
              ! The energy change does not vary monotonically with Kddt_h.  Find the maximum?
              kd(k) = kd_guess0 ; dpe_conv = pe_chg_g0
            endif

            conv_perel = conv_perel - dpe_conv
            mech_tke = mech_tke + mke_src
            if (cs%TKE_diagnostics) then
              ecd%dTKE_conv = ecd%dTKE_conv - cs%nstar*dpe_conv * i_dtdiag
              ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
            endif
            if (sfc_connected) then
              mld_output = mld_output + gv%H_to_Z * h(k)
            endif

            kddt_h(k) = kd(k) * dt_h
          elseif (tot_tke + (mke_src - pe_chg_g0) >= 0.0) then
            ! This column is convctively stable and there is energy to support the suggested
            ! mixing.  Keep that estimate.
            kd(k) = kd_guess0
            kddt_h(k) = kddt_h_g0

            ! Reduce the mechanical and convective TKE proportionately.
            tot_tke = tot_tke + mke_src
            tke_reduc = 0.0   ! tot_TKE could be 0 if Convectively_stable is false.
            if (tot_tke > 0.0) tke_reduc = (tot_tke - pe_chg_g0) / tot_tke
            if (cs%TKE_diagnostics) then
              ecd%dTKE_mixing = ecd%dTKE_mixing - pe_chg_g0 * i_dtdiag
              ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * 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 = tke_reduc*tot_tke
            mech_tke = tke_reduc*(mech_tke + mke_src)
            conv_perel = tke_reduc*conv_perel
            if (sfc_connected) then
              mld_output = mld_output + gv%H_to_Z * h(k)
            endif

          elseif (tot_tke == 0.0) then
            ! This can arise if nstar_FC = 0, but it is not common.
            kd(k) = 0.0 ; kddt_h(k) = 0.0
            tot_tke = 0.0 ; conv_perel = 0.0 ; mech_tke = 0.0
            sfc_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 = tot_tke + (mke_src - pe_chg_g0)
            tke_left_min = tot_tke

            ! As a starting guess, take the minimum of a false position estimate
            ! and a Newton's method estimate starting from Kddt_h = 0.0.
            kddt_h_guess = tot_tke * kddt_h_max / max( pe_chg_g0  - mke_src, &
                             kddt_h_max * (dpec_dkd_kd0 - dmke_max * mke2_hharm) )
            ! 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 = tot_TKE * 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
              mke_src_itt(:) = 0.0 ; kddt_h_itt(:) = 0.0
            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, &
                         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), &
                         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)
              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)

              tke_left = tot_tke + (mke_src - 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
                use_newt = .false.
              else
                dkddt_h_newt = tke_left / (dpec_dkd - dmke_src_dk)
                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(k) = kddt_h_guess / dt_h ; kddt_h(k) = kd(k) * dt_h

            ! All TKE should have been consumed.
            if (cs%TKE_diagnostics) then
              ecd%dTKE_mixing = ecd%dTKE_mixing - (tot_tke + mke_src) * i_dtdiag
              ecd%dTKE_MKE = ecd%dTKE_MKE + mke_src * i_dtdiag
              ecd%dTKE_conv_decay = ecd%dTKE_conv_decay + &
                  (cs%nstar-nstar_fc) * conv_perel * i_dtdiag
            endif

            if (sfc_connected) mld_output = mld_output + &
                 (pe_chg / (pe_chg_g0)) * gv%H_to_Z * h(k)

            tot_tke = 0.0 ; mech_tke = 0.0 ; conv_perel = 0.0
            sfc_disconnect = .true.
          endif ! End of convective or forced mixing cases to determine Kd.

          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))
          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)
          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)

        endif  ! tot_TKT > 0.0 branch.  Kddt_h(K) has been set.

        ! Store integrated velocities and thicknesses for MKE conversion calculations.
        if (sfc_disconnect) then
          ! There is no turbulence at this interface, so zero out the running sums.
          uhtot = u(k)*h(k)
          vhtot = v(k)*h(k)
          htot  = h(k)
          sfc_connected = .false.
        else
          uhtot = uhtot + u(k)*h(k)
          vhtot = vhtot + v(k)*h(k)
          htot  = htot + h(k)
        endif

        if (calc_te) then
          if (k==2) then
            te(1) = b1*(h(1)*t0(1))
            se(1) = b1*(h(1)*s0(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))
          endif
        endif
      enddo
      kd(nz+1) = 0.0

      if (calc_dt_expect) then
        ! Complete the tridiagonal solve for Te.
        b1 = 1.0 / hp_a
        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))
        ecd%dT_expect(nz) = te(nz) - t0(nz) ; ecd%dS_expect(nz) = se(nz) - s0(nz)
        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)
          ecd%dT_expect(k) = te(k) - t0(k) ; ecd%dS_expect(k) = se(k) - s0(k)
        enddo
      endif

      if (debug) then
        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
      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
      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
      endif

      if (.not.obl_converged) 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.
        ! Taking the occasional step with MLD_output empirically helps to converge faster.
        if ((dmld_min > 0.0) .and. (dmld_max < 0.0) .and. (obl_it > 2) .and. (mod(obl_it-1,4)>0)) then
          ! 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
          ! 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.
  if (cs%Use_LT) then
    ecd%LA = la ; ecd%LAmod = lamod ; ecd%mstar = mstar_total ; ecd%mstar_LT = mstar_lt
  else
    ecd%LA = 0.0 ; ecd%LAmod = 0.0 ; ecd%mstar = mstar_total ; ecd%mstar_LT = 0.0
  endif

  mld_io = mld_output

end subroutine epbl_column

!> This subroutine calculates the change in potential energy and or derivatives
!! for several changes in an interfaces'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 [degC 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 [degC H ~> degC m or degC kg m-2].
  real, intent(in)  :: Th_b     !< An effective temperature times a thickness in the layer
                                !! below, including implicit mixfing effects with other
                                !! yet lower layers [degC 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 [degC H ~> degC m or degC 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 degC-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 ppt-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 degC-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 ppt-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 degC-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 ppt-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 degC-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 ppt-1 ~> m ppt-1].

  real, optional, 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 realizedd 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].

  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 [degC H2 ~> degC m2 or degC kg2 m-4].
  real :: dS_c ! The core term in the expressions for the salinity changes [ppt 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
  ! succint 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)

  if (present(pe_chg)) then
    ! 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)
  elseif (present(pe_colht_cor)) then
    y1_3 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
    pe_colht_cor = -pres_z * min(colht_core * y1_3, 0.0)
  endif

  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 calculates the change in potential energy and or derivatives
!! for several changes in an interfaces's diapycnal diffusivity times a timestep
!! using the original form used in the first version of ePBL.
subroutine find_pe_chg_orig(Kddt_h, h_k, b_den_1, dTe_term, dSe_term, &
                       dT_km1_t2, dS_km1_t2, dT_to_dPE_k, dS_to_dPE_k, &
                       dT_to_dPEa, dS_to_dPEa, pres_Z, dT_to_dColHt_k, &
                       dS_to_dColHt_k, dT_to_dColHta, dS_to_dColHta, &
                       PE_chg, dPEc_dKd, dPE_max, dPEc_dKd_0)
  real, intent(in)  :: Kddt_h   !< 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)  :: h_k      !< The thickness of the layer below the interface [H ~> m or kg m-2].
  real, intent(in)  :: b_den_1  !< The first term in the denominator of the pivot
                                !! for the tridiagonal solver, 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)  :: dTe_term !< A diffusivity-independent term related to the temperature change
                                !! in the layer below the interface [degC H ~> degC m or degC kg m-2].
  real, intent(in)  :: dSe_term !< A diffusivity-independent term related to the salinity change
                                !! in the layer below the interface [ppt H ~> ppt m or ppt kg m-2].
  real, intent(in)  :: dT_km1_t2 !< A diffusivity-independent term related to the
                                 !! temperature change in the layer above the interface [degC].
  real, intent(in)  :: dS_km1_t2 !< A diffusivity-independent term related to the
                                 !! salinity change in the layer above the interface [ppt].
  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_dPE_k !< 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 degC-1 ~> J m-2 degC-1].
  real, intent(in)  :: dS_to_dPE_k !< 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
                                 !! in the salinities of all the layers below [R Z3 T-2 ppt-1 ~> J m-2 ppt-1].
  real, intent(in)  :: dT_to_dPEa !< 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 degC-1 ~> J m-2 degC-1].
  real, intent(in)  :: dS_to_dPEa !< 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 ppt-1 ~> J m-2 ppt-1].
  real, intent(in)  :: dT_to_dColHt_k !< 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 degC-1 ~> m degC-1].
  real, intent(in)  :: dS_to_dColHt_k !< 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 ppt-1 ~> m ppt-1].
  real, intent(in)  :: dT_to_dColHta !< 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 degC-1 ~> m degC-1].
  real, intent(in)  :: dS_to_dColHta !< 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 ppt-1 ~> m ppt-1].

  real, optional, 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 realizedd 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].

!   This subroutine determines the total potential energy change due to mixing
! at an interface, including all of the implicit effects of the prescribed
! mixing at interfaces above.  Everything here is derived by careful manipulation
! of the robust tridiagonal solvers used for tracers by MOM6.  The results are
! positive for mixing in a stably stratified environment.
!   The comments describing these arguments are for a downward mixing pass, but
! this routine can also be used for an upward pass with the sense of direction
! reversed.

  real :: b1            ! b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
  real :: b1Kd          ! Temporary array [nondim]
  real :: ColHt_chg     ! The change in column thickness [Z ~> m].
  real :: dColHt_max    ! The change in column thickness for infinite diffusivity [Z ~> m].
  real :: dColHt_dKd    ! The partial derivative of column thickness with Kddt_h [Z H-1 ~> 1 or m3 kg-2].
  real :: dT_k, dT_km1  ! Temporary arrays [degC].
  real :: dS_k, dS_km1  ! Temporary arrays [ppt].
  real :: I_Kr_denom    ! Temporary array [H-2 ~> m-2 or m4 kg-2]
  real :: dKr_dKd       ! Nondimensional temporary array.
  real :: ddT_k_dKd, ddT_km1_dKd ! Temporary arrays [degC H-1 ~> m-1 or m2 kg-1].
  real :: ddS_k_dKd, ddS_km1_dKd ! Temporary arrays [ppt H-1 ~> ppt m-1 or ppt m2 kg-1].

  b1 = 1.0 / (b_den_1 + kddt_h)
  b1kd = kddt_h*b1

  ! Start with the temperature change in layer k-1 due to the diffusivity at
  ! interface K without considering the effects of changes in layer k.

  ! Calculate the change in PE due to the diffusion at interface K
  ! if Kddt_h(K+1) = 0.
  i_kr_denom = 1.0 / (h_k*b_den_1 + (b_den_1 + h_k)*kddt_h)

  dt_k = (kddt_h*i_kr_denom) * dte_term
  ds_k = (kddt_h*i_kr_denom) * dse_term

  if (present(pe_chg)) then
    ! Find the change in energy due to diffusion with strength Kddt_h at this interface.
    ! Increment the temperature changes in layer k-1 due the changes in layer k.
    dt_km1 = b1kd * ( dt_k + dt_km1_t2 )
    ds_km1 = b1kd * ( ds_k + ds_km1_t2 )
    pe_chg = (dt_to_dpe_k * dt_k + dt_to_dpea * dt_km1) + &
             (ds_to_dpe_k * ds_k + ds_to_dpea * ds_km1)
    colht_chg = (dt_to_dcolht_k * dt_k + dt_to_dcolhta * dt_km1) + &
                (ds_to_dcolht_k * ds_k + ds_to_dcolhta * ds_km1)
    if (colht_chg < 0.0) pe_chg = pe_chg - pres_z * colht_chg
  endif

  if (present(dpec_dkd)) then
    ! Find the derivatives of the temperature and salinity changes with Kddt_h.
    dkr_dkd = (h_k*b_den_1) * i_kr_denom**2

    ddt_k_dkd = dkr_dkd * dte_term
    dds_k_dkd = dkr_dkd * dse_term
    ddt_km1_dkd = (b1**2 * b_den_1) * ( dt_k + dt_km1_t2 ) + b1kd * ddt_k_dkd
    dds_km1_dkd = (b1**2 * b_den_1) * ( ds_k + ds_km1_t2 ) + b1kd * dds_k_dkd

    ! Calculate the partial derivative of Pe_chg with Kddt_h.
    dpec_dkd = (dt_to_dpe_k * ddt_k_dkd + dt_to_dpea * ddt_km1_dkd) + &
               (ds_to_dpe_k * dds_k_dkd + ds_to_dpea * dds_km1_dkd)
    dcolht_dkd = (dt_to_dcolht_k * ddt_k_dkd + dt_to_dcolhta * ddt_km1_dkd) + &
                 (ds_to_dcolht_k * dds_k_dkd + ds_to_dcolhta * dds_km1_dkd)
    if (dcolht_dkd < 0.0) dpec_dkd = dpec_dkd - pres_z * dcolht_dkd
  endif

  if (present(dpe_max)) then
    ! This expression is the limit of PE_chg for infinite Kddt_h.
    dpe_max = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) + &
              ((dt_to_dpe_k + dt_to_dpea) * dte_term + &
               (ds_to_dpe_k + ds_to_dpea) * dse_term) / (b_den_1 + h_k)
    dcolht_max = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) + &
              ((dt_to_dcolht_k + dt_to_dcolhta) * dte_term + &
               (ds_to_dcolht_k + ds_to_dcolhta) * dse_term) / (b_den_1 + h_k)
    if (dcolht_max < 0.0) dpe_max = dpe_max - pres_z*dcolht_max
  endif

  if (present(dpec_dkd_0)) then
    ! This expression is the limit of dPEc_dKd for Kddt_h = 0.
    dpec_dkd_0 = (dt_to_dpea * dt_km1_t2 + ds_to_dpea * ds_km1_t2) / (b_den_1) + &
                 (dt_to_dpe_k * dte_term + ds_to_dpe_k * dse_term) / (h_k*b_den_1)
    dcolht_dkd = (dt_to_dcolhta * dt_km1_t2 + ds_to_dcolhta * ds_km1_t2) / (b_den_1) + &
                 (dt_to_dcolht_k * dte_term + ds_to_dcolht_k * dse_term) / (h_k*b_den_1)
    if (dcolht_dkd < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres_z*dcolht_dkd
  endif

end subroutine find_pe_chg_orig

!> This subroutine finds the Mstar value for ePBL
subroutine find_mstar(CS, US, Buoyancy_Flux, UStar, UStar_Mean,&
                      BLD, Abs_Coriolis, MStar, Langmuir_Number,&
                      MStar_LT, Convect_Langmuir_Number)
  type(energetic_PBL_CS), pointer    :: CS    !< Energetic_PBL control structure.
  type(unit_scale_type), intent(in)  :: US    !< A dimensional unit scaling type
  real,                  intent(in)  :: UStar !< ustar w/ gustiness [Z T-1 ~> m s-1]
  real,                  intent(in)  :: UStar_Mean !< ustar w/o gustiness [Z T-1 ~> m s-1]
  real,                  intent(in)  :: Abs_Coriolis !< abolute value of the Coriolis parameter [T-1 ~> s-1]
  real,                  intent(in)  :: Buoyancy_Flux !< Buoyancy flux [Z2 T-3 ~> m2 s-3]
  real,                  intent(in)  :: BLD   !< boundary layer depth [Z ~> m]
  real,                  intent(out) :: Mstar !< Ouput mstar (Mixing/ustar**3) [nondim]
  real,        optional, intent(in)  :: Langmuir_Number !< Langmuir number [nondim]
  real,        optional, intent(out) :: MStar_LT !< Mstar increase due to Langmuir turbulence [nondim]
  real,        optional, intent(out) :: Convect_Langmuir_number !< Langmuir number including buoyancy flux [nondim]

  !/ Variables used in computing mstar
  real :: MSN_term       ! Temporary terms [nondim]
  real :: MSCR_term1, MSCR_term2 ! Temporary terms [Z3 T-3 ~> m3 s-3]
  real :: MStar_Conv_Red ! Adjustment made to mstar due to convection reducing mechanical mixing [nondim]
  real :: MStar_S, MStar_N ! Mstar in (S)tabilizing/(N)ot-stabilizing buoyancy flux [nondim]

  !/  Integer options for how to find mstar

  !/

  if (cs%mstar_scheme == use_fixed_mstar) then
    mstar = cs%Fixed_MStar
  !/ 1. Get mstar
  elseif (cs%mstar_scheme == mstar_from_ekman) then

    if (cs%answers_2018) then
      ! The limit for the balance of rotation and stabilizing is f(L_Ekman,L_Obukhov)
      mstar_s = cs%MStar_coef*sqrt(max(0.0,buoyancy_flux) / ustar**2 / &
                    (abs_coriolis + 1.e-10*us%T_to_s) )
      ! The limit for rotation (Ekman length) limited mixing
      mstar_n =  cs%C_Ek * log( max( 1., ustar / (abs_coriolis + 1.e-10*us%T_to_s) / bld ) )
    else
      ! The limit for the balance of rotation and stabilizing is f(L_Ekman,L_Obukhov)
      mstar_s = cs%MSTAR_COEF*sqrt(max(0.0, buoyancy_flux) / (ustar**2 * max(abs_coriolis, 1.e-20*us%T_to_s)))
      ! The limit for rotation (Ekman length) limited mixing
      mstar_n = 0.0
      if (ustar > abs_coriolis * bld) mstar_n = cs%C_EK * log(ustar / (abs_coriolis * bld))
    endif

    ! Here 1.25 is about .5/von Karman, which gives the Obukhov limit.
    mstar = max(mstar_s, min(1.25, mstar_n))
    if (cs%MStar_Cap > 0.0) mstar = min( cs%MStar_Cap,mstar )
  elseif ( cs%mstar_scheme == mstar_from_rh18 ) then
    if (cs%answers_2018) then
      mstar_n = cs%RH18_MStar_cn1 * ( 1.0 - 1.0 / ( 1. + cs%RH18_MStar_cn2 * &
                exp( cs%RH18_mstar_CN3 * bld * abs_coriolis / ustar) ) )
    else
      msn_term = cs%RH18_MStar_cn2 * exp( cs%RH18_mstar_CN3 * bld * abs_coriolis / ustar)
      mstar_n = (cs%RH18_MStar_cn1 *  msn_term) / ( 1. + msn_term)
    endif
    mstar_s = cs%RH18_MStar_CS1 * ( max(0.0, buoyancy_flux)**2 * bld / &
             ( ustar**5 * max(abs_coriolis,1.e-20*us%T_to_s) ) )**cs%RH18_mstar_cs2
    mstar = mstar_n + mstar_s
  endif

  !/ 2. Adjust mstar to account for convective turbulence
  if (cs%answers_2018) then
    mstar_conv_red = 1. - cs%MStar_Convect_coef * (-min(0.0,buoyancy_flux) + 1.e-10*us%T_to_s**3*us%m_to_Z**2) / &
                         ( (-min(0.0,buoyancy_flux) + 1.e-10*us%T_to_s**3*us%m_to_Z**2) + &
                         2.0 *mstar * ustar**3 / bld )
  else
    mscr_term1 = -bld * min(0.0, buoyancy_flux)
    mscr_term2 = 2.0*mstar * ustar**3
    if ( abs(mscr_term2) > 0.0) then
      mstar_conv_red = ((1.-cs%mstar_convect_coef) * mscr_term1 + mscr_term2) / (mscr_term1 + mscr_term2)
    else
      mstar_conv_red = 1.-cs%mstar_convect_coef
    endif
  endif

  !/3. Combine various mstar terms to get final value
  mstar = mstar * mstar_conv_red

  if (present(langmuir_number)) then
    !### In this call, ustar was previously ustar_mean.  Is this change deliberate, Brandon?  -RWH
    call mstar_langmuir(cs, us, abs_coriolis, buoyancy_flux, ustar, bld, langmuir_number, mstar, &
                        mstar_lt, convect_langmuir_number)
  endif

end subroutine find_mstar

!> This subroutine modifies the Mstar value if the Langmuir number is present
subroutine mstar_langmuir(CS, US, Abs_Coriolis, Buoyancy_Flux, UStar, BLD, Langmuir_Number, &
                          Mstar, MStar_LT, Convect_Langmuir_Number)
  type(energetic_PBL_CS), pointer    :: CS    !< Energetic_PBL control structure.
  type(unit_scale_type), intent(in)  :: US    !< A dimensional unit scaling type
  real,                  intent(in)  :: Abs_Coriolis !< Absolute value of the Coriolis parameter [T-1 ~> s-1]
  real,                  intent(in)  :: Buoyancy_Flux !< Buoyancy flux [Z2 T-3 ~> m2 s-3]
  real,                  intent(in)  :: UStar !< Surface friction velocity with? gustiness [Z T-1 ~> m s-1]
  real,                  intent(in)  :: BLD   !< boundary layer depth [Z ~> m]
  real,                  intent(inout) :: Mstar !< Input/output mstar (Mixing/ustar**3) [nondim]
  real,                  intent(in)  :: Langmuir_Number !< Langmuir number [nondim]
  real,                  intent(out) :: MStar_LT !< Mstar increase due to Langmuir turbulence [nondim]
  real,                  intent(out) :: Convect_Langmuir_number !< Langmuir number including buoyancy flux [nondim]

  !/
  real, parameter :: Max_ratio = 1.0e16  ! The maximum value of a nondimensional ratio.
  real :: enhance_mstar ! A multiplicative scaling of mstar due to Langmuir turbulence.
  real :: mstar_LT_add ! A value that is added to mstar due to Langmuir turbulence.
  real :: iL_Ekman    ! Inverse of Ekman length scale [Z-1 ~> m-1].
  real :: iL_Obukhov  ! Inverse of Obukhov length scale [Z-1 ~> m-1].
  real :: I_ustar     ! The Adcroft reciprocal of ustar [T Z-1 ~> s m-1]
  real :: I_f         ! The Adcroft reciprocal of the Coriolis parameter [T ~> s]
  real :: MLD_Ekman          ! The ratio of the mixed layer depth to the Ekman layer depth [nondim].
  real :: Ekman_Obukhov      ! The Ekman layer thickness divided by the Obukhov depth [nondim].
  real :: MLD_Obukhov        ! The mixed layer depth divided by the Obukhov depth [nondim].
  real :: MLD_Obukhov_stab   ! Ratios of length scales where MLD is boundary layer depth [nondim].
  real :: Ekman_Obukhov_stab ! >
  real :: MLD_Obukhov_un     ! Ratios of length scales where MLD is boundary layer depth
  real :: Ekman_Obukhov_un   ! >

  ! Set default values for no Langmuir effects.
  enhance_mstar = 1.0 ; mstar_lt_add = 0.0

  if (cs%LT_Enhance_Form /= no_langmuir) then
    ! a. Get parameters for modified LA
    if (cs%answers_2018) then
      il_ekman   = abs_coriolis / ustar
      il_obukhov = buoyancy_flux*cs%vonkar / ustar**3
      ekman_obukhov_stab = abs(max(0., il_obukhov / (il_ekman + 1.e-10*us%Z_to_m)))
      ekman_obukhov_un = abs(min(0., il_obukhov / (il_ekman + 1.e-10*us%Z_to_m)))
      mld_obukhov_stab = abs(max(0., bld*il_obukhov))
      mld_obukhov_un = abs(min(0., bld*il_obukhov))
      mld_ekman = abs( bld*il_ekman )
    else
      ekman_obukhov = max_ratio ; mld_obukhov = max_ratio ; mld_ekman = max_ratio
      i_f = 0.0 ; if (abs(abs_coriolis) > 0.0) i_f = 1.0 / abs_coriolis
      i_ustar = 0.0 ; if (abs(ustar) > 0.0) i_ustar = 1.0 / ustar
      if (abs(buoyancy_flux*cs%vonkar) < max_ratio*(abs_coriolis * ustar**2)) &
        ekman_obukhov = abs(buoyancy_flux*cs%vonkar) * (i_f * i_ustar**2)
      if (abs(bld*buoyancy_flux*cs%vonkar) < max_ratio*ustar**3) &
        mld_obukhov = abs(bld*buoyancy_flux*cs%vonkar) * i_ustar**3
      if (bld*abs_coriolis < max_ratio*ustar) &
        mld_ekman = bld*abs_coriolis * i_ustar

      if (buoyancy_flux > 0.0) then
        ekman_obukhov_stab = ekman_obukhov ; ekman_obukhov_un = 0.0
        mld_obukhov_stab = mld_obukhov ; mld_obukhov_un = 0.0
      else
        ekman_obukhov_un = ekman_obukhov ; ekman_obukhov_stab = 0.0
        mld_obukhov_un = mld_obukhov ; mld_obukhov_stab = 0.0
      endif
    endif

    ! b. Adjust LA based on various parameters.
    !    Assumes linear factors based on length scale ratios to adjust LA
    !    Note when these coefficients are set to 0 recovers simple LA.
    convect_langmuir_number = langmuir_number * &
                    ( (1.0 + max(-0.5, cs%LaC_MLDoEK * mld_ekman)) + &
                   ((cs%LaC_EKoOB_stab * ekman_obukhov_stab + cs%LaC_EKoOB_un * ekman_obukhov_un) + &
                    (cs%LaC_MLDoOB_stab * mld_obukhov_stab  + cs%LaC_MLDoOB_un * mld_obukhov_un)) )

    if (cs%LT_Enhance_Form == langmuir_rescale) then
      ! Enhancement is multiplied (added mst_lt set to 0)
      enhance_mstar = min(cs%Max_Enhance_M, &
                          (1. + cs%LT_ENHANCE_COEF * convect_langmuir_number**cs%LT_ENHANCE_EXP) )
    elseif (cs%LT_ENHANCE_Form == langmuir_add) then
      ! or Enhancement is additive (multiplied enhance_m set to 1)
      mstar_lt_add = cs%LT_ENHANCE_COEF * convect_langmuir_number**cs%LT_ENHANCE_EXP
    endif
  endif

  mstar_lt = (enhance_mstar - 1.0)*mstar + mstar_lt_add  ! Diagnose the full increase in mstar.
  mstar = mstar*enhance_mstar + mstar_lt_add

end subroutine mstar_langmuir


!> Copies the ePBL active mixed layer depth into MLD, in units of [Z ~> m] unless other units are specified.
subroutine energetic_pbl_get_mld(CS, MLD, G, US, m_to_MLD_units)
  type(energetic_pbl_cs),           pointer     :: cs  !< Control structure for ePBL
  type(ocean_grid_type),            intent(in)  :: g   !< Grid structure
  type(unit_scale_type),            intent(in)  :: us  !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: mld !< Depth of ePBL active mixing layer [Z ~> m] or other units
  real,                   optional, intent(in)  :: m_to_mld_units !< A conversion factor from meters
                                                       !! to the desired units for MLD
  ! Local variables
  real :: scale  ! A dimensional rescaling factor
  integer :: i,j

  scale = 1.0 ; if (present(m_to_mld_units)) scale = us%Z_to_m * m_to_mld_units

  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    mld(i,j) = scale*cs%ML_Depth(i,j)
  enddo ; enddo

end subroutine energetic_pbl_get_mld


!> This subroutine initializes the energetic_PBL module
subroutine energetic_pbl_init(Time, G, GV, US, param_file, diag, CS)
  type(time_type), target, intent(in)    :: time !< The current model time
  type(ocean_grid_type),   intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type), intent(in)    :: gv   !< The ocean's vertical grid structure
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic output
  type(energetic_pbl_cs),  pointer       :: cs   !< A pointer that is set to point to the control
                                                 !! structure for this module
  ! Local variables
  ! 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
  real :: omega_frac_dflt
  integer :: isd, ied, jsd, jed
  integer :: mstar_mode, lt_enhance, wt_mode
  logical :: default_2018_answers
  logical :: use_temperature, use_omega
  logical :: use_la_windsea
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  if (associated(cs)) then
    call mom_error(warning, "mixedlayer_init called with an associated"//&
                            "associated control structure.")
    return
  else ; allocate(cs) ; endif

  cs%diag => diag
  cs%Time => time

! Set default, read and log parameters
  call log_version(param_file, mdl, version, "")


!/1. General ePBL settings
  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)
  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
                 "If true, use the absolute rotation rate instead of the "//&
                 "vertical component of rotation when setting the decay "//&
                 "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.")
    omega_frac_dflt = 1.0
  endif
  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
                 "When setting the decay scale for turbulence, use this "//&
                 "fraction of the absolute rotation rate blended with the "//&
                 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
                 units="nondim", default=omega_frac_dflt)
  call get_param(param_file, mdl, "EKMAN_SCALE_COEF", cs%Ekman_scale_coef, &
                 "A nondimensional scaling factor controlling the inhibition "//&
                 "of the diffusive length scale by rotation. Making this larger "//&
                 "decreases the PBL diffusivity.", units="nondim", default=1.0)
  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
  call get_param(param_file, mdl, "EPBL_2018_ANSWERS", cs%answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)


  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.)

  call get_param(param_file, mdl, "MKE_TO_TKE_EFFIC", cs%MKE_to_TKE_effic, &
                 "The efficiency with which mean kinetic energy released "//&
                 "by mechanically forced entrainment of the mixed layer "//&
                 "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)


!/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 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)
  if (mstar_mode == 0) then
    tmpstr = constant_string
    call mom_error(warning, "Use EPBL_MSTAR_SCHEME = CONSTANT instead of the archaic MSTAR_MODE = 0.")
  elseif (mstar_mode == 1) then
    call mom_error(fatal, "You are using a legacy mstar mode in ePBL that has been phased out. "//&
                          "If you need to use this setting please report this error.  Also use "//&
                          "EPBL_MSTAR_SCHEME to specify the scheme for mstar.")
  elseif (mstar_mode == 2) then
    tmpstr = om4_string
    call mom_error(warning, "Use EPBL_MSTAR_SCHEME = OM4 instead of the archaic MSTAR_MODE = 2.")
  elseif (mstar_mode == 3) then
    tmpstr = rh18_string
    call mom_error(warning, "Use EPBL_MSTAR_SCHEME = REICHL_H18 instead of the archaic MSTAR_MODE = 3.")
  elseif (mstar_mode > 3) then
    call mom_error(fatal, "An unrecognized value of the obsolete parameter MSTAR_MODE was specified.")
  endif
  call log_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 REICHL_H18 - Use the scheme documented in Reichl & Hallberg, 2018.", &
                 default=constant_string)
  tmpstr = uppercase(tmpstr)
  select case (tmpstr)
    case (constant_string)
      cs%mstar_Scheme = use_fixed_mstar
    case (om4_string)
      cs%mstar_Scheme = mstar_from_ekman
    case (rh18_string)
      cs%mstar_Scheme = mstar_from_rh18
    case default
      call mom_mesg('energetic_PBL_init: EPBL_MSTAR_SCHEME ="'//trim(tmpstr)//'"', 0)
      call mom_error(fatal, "energetic_PBL_init: Unrecognized setting "// &
            "EPBL_MSTAR_SCHEME = "//trim(tmpstr)//" found in input file.")
  end select

  call get_param(param_file, mdl, "MSTAR", cs%fixed_mstar, &
                 "The ratio of the friction velocity cubed to the TKE input to the "//&
                 "mixed layer.  This option is used if EPBL_MSTAR_SCHEME = CONSTANT.", &
                 units="nondim", default=1.2, do_not_log=(cs%mstar_scheme/=use_fixed_mstar))
  call get_param(param_file, mdl, "MSTAR_CAP", cs%mstar_cap, &
                 "If this value is positive, it sets the maximum value of mstar "//&
                 "allowed in ePBL.  (This is not used if EPBL_MSTAR_SCHEME = CONSTANT).", &
                 units="nondim", default=-1.0, do_not_log=(cs%mstar_scheme==use_fixed_mstar))
  ! mstar_scheme==MStar_from_Ekman options
  call get_param(param_file, mdl, "MSTAR2_COEF1", cs%MSTAR_COEF, &
                 "Coefficient in computing mstar when rotation and stabilizing "//&
                 "effects are both important (used if EPBL_MSTAR_SCHEME = OM4).", &
                 units="nondim", default=0.3, do_not_log=(cs%mstar_scheme/=mstar_from_ekman))
  call get_param(param_file, mdl, "MSTAR2_COEF2", cs%C_EK, &
                 "Coefficient in computing mstar when only rotation limits "// &
                 "the total mixing (used if EPBL_MSTAR_SCHEME = OM4)", &
                 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). "//&
                 "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).", &
                 units="nondim", default=0.275, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
  call get_param(param_file, mdl, "RH18_MSTAR_CN2", cs%RH18_mstar_cn2,&
                 "MSTAR_N coefficient 2 (coefficient outside of exponential decay). "//&
                 "The value of 8.0 is given in RH18.  Increasing this coefficient "//&
                 "increases MSTAR for all values of HF/ust, with a much more even "//&
                 "effect across a wide range of Hf/ust than CN1.", &
                 units="nondim", default=8.0, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
  call get_param(param_file, mdl, "RH18_MSTAR_CN3", cs%RH18_mstar_CN3,&
                 "MSTAR_N coefficient 3 (exponential decay coefficient). "//&
                 "The value of -5.0 is given in RH18.  Increasing this increases how "//&
                 "quickly the value of MSTAR decreases as Hf/ust increases.", &
                  units="nondim", default=-5.0, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
  call get_param(param_file, mdl, "RH18_MSTAR_CS1", cs%RH18_mstar_cs1,&
                 "MSTAR_S coefficient for RH18 in stabilizing limit. "//&
                 "The value of 0.2 is given in RH18 and increasing it increases "//&
                 "MSTAR in the presence of a stabilizing surface buoyancy flux.", &
                 units="nondim", default=0.2, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))
  call get_param(param_file, mdl, "RH18_MSTAR_CS2", cs%RH18_mstar_cs2,&
                 "MSTAR_S exponent for RH18 in stabilizing limit. "//&
                 "The value of 0.4 is given in RH18 and increasing it increases MSTAR "//&
                 "exponentially in the presence of a stabilizing surface buoyancy flux.", &
                 units="nondim", default=0.4, do_not_log=(cs%mstar_scheme/=mstar_from_rh18))


!/ Convective turbulence related options
  call get_param(param_file, mdl, "NSTAR", cs%nstar, &
                 "The portion of the buoyant potential energy imparted by "//&
                 "surface fluxes that is available to drive entrainment "//&
                 "at the base of mixed layer when that energy is positive.", &
                 units="nondim", default=0.2)
  call get_param(param_file, mdl, "MSTAR_CONV_ADJ", cs%mstar_convect_coef, &
                 "Coefficient used for reducing mstar during convection "//&
                 "due to reduction of stable density gradient.", &
                 units="nondim", default=0.0)

!/ Mixing Length Options
  call get_param(param_file, mdl, "USE_MLD_ITERATION", cs%Use_MLD_iteration, &
                 "A logical that specifies whether or not to use the "//&
                 "distance to the bottom of the actively turbulent boundary "//&
                 "layer to help set the EPBL length scale.", default=.true.)
  call get_param(param_file, mdl, "EPBL_TRANSITION_SCALE", cs%transLay_scale, &
                 "A scale for the mixing length in the transition layer "//&
                 "at the edge of the boundary layer as a fraction of the "//&
                 "boundary layer thickness.", units="nondim", default=0.1)
  if ( cs%Use_MLD_iteration .and. abs(cs%transLay_scale-0.5) >= 0.5) then
    call mom_error(fatal, "If flag USE_MLD_ITERATION is true, then "//&
                 "EPBL_TRANSITION should be greater than 0 and less than 1.")
  endif

  call get_param(param_file, mdl, "MLD_ITERATION_GUESS", cs%MLD_ITERATION_GUESS, &
                 "If true, use the previous timestep MLD as a first guess in the MLD iteration, "//&
                 "otherwise use half the ocean depth as the first guess of the boundary layer "//&
                 "depth.  The default is false to facilitate reproducibility.", &
                 default=.false., do_not_log=.not.cs%Use_MLD_iteration)
  call get_param(param_file, mdl, "EPBL_MLD_TOLERANCE", cs%MLD_tol, &
                 "The tolerance for the iteratively determined mixed "//&
                 "layer depth.  This is only used with USE_MLD_ITERATION.", &
                 units="meter", default=1.0, scale=us%m_to_Z, do_not_log=.not.cs%Use_MLD_iteration)
  call get_param(param_file, mdl, "EPBL_MLD_BISECTION", cs%MLD_bisection, &
                 "If true, use bisection with the iterative determination of the self-consistent "//&
                 "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=.true., do_not_log=.not.cs%Use_MLD_iteration) !### The default should become 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.", &
                 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, &
                 "The minimum mixing length scale that will be used "//&
                 "by ePBL.  The default (0) does not set a minimum.", &
                 units="meter", default=0.0, scale=us%m_to_Z)

  call get_param(param_file, mdl, "MIX_LEN_EXPONENT", cs%MixLenExponent, &
                 "The exponent applied to the ratio of the distance to the MLD "//&
                 "and the MLD depth which determines the shape of the mixing length. "//&
                 "This is only used if USE_MLD_ITERATION is True.", &
                 units="nondim", default=2.0)

!/ Turbulent velocity scale in mixing coefficient
  call get_param(param_file, mdl, "EPBL_VEL_SCALE_SCHEME", tmpstr, &
                 "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, 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
    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
    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, &
                 "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)
    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_error(fatal, "energetic_PBL_init: Unrecognized setting "// &
            "EPBL_VEL_SCALE_SCHEME = "//trim(tmpstr)//" found in input file.")
  end select

  call get_param(param_file, mdl, "WSTAR_USTAR_COEF", cs%wstar_ustar_coef, &
                 "A ratio relating the efficiency with which convectively "//&
                 "released energy is converted to a turbulent velocity, "//&
                 "relative to mechanically forced TKE. Making this larger "//&
                 "increases the BL diffusivity", units="nondim", default=1.0)
  call get_param(param_file, mdl, "EPBL_VEL_SCALE_FACTOR", cs%vstar_scale_fac, &
                 "An overall nondimensional scaling factor for wT. "//&
                 "Making this larger increases the PBL diffusivity.", &
                 units="nondim", default=1.0)
  call get_param(param_file, mdl, "VSTAR_SURF_FAC", cs%vstar_surf_fac,&
                 "The proportionality times ustar to set vstar at the surface.", &
                 units="nondim", default=1.2)

  !/ 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 "//&
       "determine the Langmuir number.", units="nondim", default=.false.)
  ! Note this can be activated in other ways, but this preserves the old method.
  if (use_la_windsea) then
    cs%USE_LT = .true.
  else
    call get_param(param_file, mdl, "EPBL_LT", cs%USE_LT, &
                 "A logical to use a LT parameterization.", &
                 units="nondim", default=.false.)
  endif
  if (cs%USE_LT) then
    call get_param(param_file, mdl, "EPBL_LANGMUIR_SCHEME", tmpstr, &
                 "EPBL_LANGMUIR_SCHEME selects the method for including Langmuir turbulence. "//&
                 "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", &
                 default=none_string, do_not_log=.true.)
    call get_param(param_file, mdl, "LT_ENHANCE", lt_enhance, default=-1)
    if (lt_enhance == 0) then
      tmpstr = none_string
      call mom_error(warning, "Use EPBL_LANGMUIR_SCHEME = NONE instead of the archaic LT_ENHANCE = 0.")
    elseif (lt_enhance == 1) then
      call mom_error(fatal, "You are using a legacy LT_ENHANCE mode in ePBL that has been phased out. "//&
                            "If you need to use this setting please report this error.  Also use "//&
                            "EPBL_LANGMUIR_SCHEME to specify the scheme for mstar.")
    elseif (lt_enhance == 2) then
      tmpstr = rescaled_string
      call mom_error(warning, "Use EPBL_LANGMUIR_SCHEME = RESCALE instead of the archaic LT_ENHANCE = 2.")
    elseif (lt_enhance == 3) then
      tmpstr = additive_string
      call mom_error(warning, "Use EPBL_LANGMUIR_SCHEME = ADDITIVE instead of the archaic LT_ENHANCE = 3.")
    elseif (lt_enhance > 3) then
      call mom_error(fatal, "An unrecognized value of the obsolete parameter LT_ENHANCE was specified.")
    endif
    call log_param(param_file, mdl, "EPBL_LANGMUIR_SCHEME", tmpstr, &
                 "EPBL_LANGMUIR_SCHEME selects the method for including Langmuir turbulence. "//&
                 "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", &
                 default=none_string)
    tmpstr = uppercase(tmpstr)
    select case (tmpstr)
      case (none_string)
        cs%LT_enhance_form = no_langmuir
      case (rescaled_string)
        cs%LT_enhance_form = langmuir_rescale
      case (additive_string)
        cs%LT_enhance_form = langmuir_add
      case default
        call mom_mesg('energetic_PBL_init: EPBL_LANGMUIR_SCHEME ="'//trim(tmpstr)//'"', 0)
        call mom_error(fatal, "energetic_PBL_init: Unrecognized setting "// &
              "EPBL_LANGMUIR_SCHEME = "//trim(tmpstr)//" found in input file.")
    end select

    call get_param(param_file, mdl, "LT_ENHANCE_COEF", cs%LT_ENHANCE_COEF, &
                 "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", &
                 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 "//&
                 "MLD approaching Ekman depth.", &
                 units="nondim", default=-0.87,  do_not_log=(cs%LT_enhance_form==no_langmuir))
    call get_param(param_file, mdl, "LT_MOD_LAC2", cs%LaC_MLDoOB_stab, &
                 "Coefficient for modification of Langmuir number due to "//&
                 "MLD approaching stable Obukhov depth.", &
                 units="nondim", default=0.0,  do_not_log=(cs%LT_enhance_form==no_langmuir))
    call get_param(param_file, mdl, "LT_MOD_LAC3", cs%LaC_MLDoOB_un, &
                 "Coefficient for modification of Langmuir number due to "//&
                 "MLD approaching unstable Obukhov depth.", &
                 units="nondim", default=0.0,  do_not_log=(cs%LT_enhance_form==no_langmuir))
    call get_param(param_file, mdl, "LT_MOD_LAC4", cs%Lac_EKoOB_stab, &
                 "Coefficient for modification of Langmuir number due to "//&
                 "ratio of Ekman to stable Obukhov depth.", &
                 units="nondim", default=0.95,  do_not_log=(cs%LT_enhance_form==no_langmuir))
    call get_param(param_file, mdl, "LT_MOD_LAC5", cs%Lac_EKoOB_un, &
                 "Coefficient for modification of Langmuir number due to "//&
                 "ratio of Ekman to unstable Obukhov depth.", &
                 units="nondim", default=0.95,  do_not_log=(cs%LT_enhance_form==no_langmuir))
  endif


!/ Logging parameters
  ! This gives a minimum decay scale that is typically much less than Angstrom.
  cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_Z + gv%H_to_Z*gv%H_subroundoff)
  call log_param(param_file, mdl, "!EPBL_USTAR_MIN", cs%ustar_min*us%Z_to_m*us%s_to_T, &
                 "The (tiny) minimum friction velocity used within the "//&
                 "ePBL code, derived from OMEGA and ANGSTROM.", units="m s-1", &
                 like_default=.true.)


!/ 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, &
      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)
  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)
  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)
  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)
  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)
  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)
  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)
  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)
  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)
  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)
  cs%id_MSTAR_mix = register_diag_field('ocean_model', 'MSTAR', diag%axesT1, &
      time, 'Total mstar that is used.', 'nondim')

  if (cs%use_LT) then
    cs%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, &
        time, 'Langmuir number.', 'nondim')
    cs%id_LA_mod = register_diag_field('ocean_model', 'LA_MOD', diag%axesT1, &
        time, 'Modified Langmuir number.', 'nondim')
    cs%id_MSTAR_LT = register_diag_field('ocean_model', 'MSTAR_LT', diag%axesT1, &
        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 (report_avg_its) then
    cs%sum_its(1) = real_to_efp(0.0) ; cs%sum_its(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.
  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

end subroutine energetic_pbl_init

!> Clean up and deallocate memory associated with the energetic_PBL module.
subroutine energetic_pbl_end(CS)
  type(energetic_pbl_cs), pointer :: cs !< Energetic_PBL control structure that
                                        !! will be deallocated in this subroutine.

  character(len=256) :: mesg
  real :: avg_its

  if (.not.associated(cs)) return

  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 (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)
  endif

  deallocate(cs)

end subroutine energetic_pbl_end

!> \namespace MOM_energetic_PBL
!!
!! By Robert Hallberg, 2015.
!!
!!   This file contains the subroutine (energetic_PBL) that uses an
!! integrated boundary layer energy budget (like a bulk- or refined-
!! bulk mixed layer scheme), but instead of homogenizing this model
!! calculates a finite diffusivity and viscosity, which in this
!! regard is conceptually similar to what is done with KPP or various
!! two-equation closures.  However, the scheme that is implemented
!! here has the big advantage that is entirely implicit, but is
!! simple enough that it requires only a single vertical pass to
!! determine the diffusivity. The development of bulk mixed layer
!! models stems from the work of various people, as described in the
!! review paper by Niiler and Kraus (1979). The work here draws in
!! with particular on the form for TKE decay proposed by Oberhuber
!! (JPO, 1993, 808-829), with an extension to a refined bulk mixed
!! layer as described in Hallberg (Aha Huliko'a, 2003).  The physical
!! processes portrayed in this subroutine include convectively driven
!! mixing and mechanically driven mixing.  Unlike boundary-layer
!! mixing, stratified shear mixing is not a one-directional turbulent
!! process, and it is dealt with elsewhere in the MOM6 code within
!! the module MOM_kappa_shear.F90.  It is assumed that the heat,
!! mass, and salt fluxes have been applied elsewhere, but that their
!! implications for the integrated TKE budget have been captured in
!! an array that is provided as an argument to this subroutine.  This
!! is a full 3-d array due to the effects of penetrating shortwave
!! radiation.

end module mom_energetic_pbl