MOM_energetic_PBL.F90

Go to the documentation of this file.

module mom_energetic_pbl
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  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.                                                         *
!*                                                                     *
!*  Macros written all in capital letters are defined in MOM_memory.h. *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v                                        *
!*    j    x ^ x ^ x   At >:  u                                        *
!*    j    > o > o >   At o:  h, T, S, Kd, etc.                        *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
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
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_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
! use MOM_EOS, only : calculate_density, calculate_density_derivs

implicit none ; private

#include <MOM_memory.h>

public energetic_pbl, energetic_pbl_init, energetic_pbl_end
public energetic_pbl_get_mld

type, public :: energetic_pbl_cs ; private
  real    :: mstar           ! 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.
  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.
  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    :: tke_decay       ! The ratio of the natural Ekman depth to the TKE
                             ! decay scale, nondimensional.
  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    :: Hmix_min        ! The minimum mixed layer thickness in m.
  real    :: ustar_min       ! A minimum value of ustar to avoid numerical
                             ! problems, in m s-1.  If the value is small enough,
                             ! this should not affect the solution.
  real    :: omega           !   The Earth's rotation rate, in 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-of)*f^2 + of*4*omega^2).
  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_scale_fac ! An overall nondimensional scaling factor
                             ! for vstar.  Making this larger increases the
                             ! diffusivity.
  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, in m.
  real    :: min_mix_len     ! The minimum mixing length scale that will be
                             ! used by ePBL, in m.  The default (0) does not
                             ! set a minimum.
  real    :: n2_dissipation_scale_neg
  real    :: n2_dissipation_scale_pos
                             ! A nondimensional scaling factor controlling the
                             ! loss of TKE due to enhanced dissipation in the presence
                             ! of stratification.  This dissipation is applied to the
                             ! available TKE which includes both that generated at the
                             ! surface and that generated at depth.  It may be important
                             ! to distinguish which TKE flavor that this dissipation
                             ! applies to in subsequent revisions of this code.
                             ! "_Neg" and "_Pos" refer to which scale is applied as a
                             ! function of negative or positive local buoyancy.
  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.
  real    :: mstar_slope     ! Slope of the function which relates the shear production
                             ! to the mixing layer depth, Ekman depth, and Monin-Obukhov
                             ! depth.
  real    :: mstar_xint      ! Value where MSTAR function transitions from linear
                             ! to decay toward MSTAR->0 at fully developed Ekman depth.
  real    :: mstar_xint_up   ! Similar but for transition to asymptotic cap.
  real    :: mstar_at_xint   ! Intercept value of MSTAR at value where function
                             ! changes to linear transition.
  real    :: lt_enhance_coef ! Coefficient in fit for Langmuir Enhancment
  real    :: lt_enhance_exp  ! Exponent in fit for Langmuir Enhancement
  real :: mstar_n = -2.      ! Exponent in decay at negative and positive limits of MLD_over_STAB
  real :: mstar_a,mstar_a2   ! MSTAR_A and MSTAR_B are coefficients in asymptote toward limits.
  real :: mstar_b,mstar_b2   !  These are computed to match the function value and slope at both
                             !  ends of the linear fit within the well constrained region.
  real :: c_ek = 0.17        ! MSTAR Coefficient in rotation limit for mstar_mode=2
  real :: mstar_coef = 0.3   ! MSTAR coefficient in rotation/stabilizing balance for mstar_mode=2
  real :: lac_mldoek         ! Coefficients for Langmuir number modification based on
  real :: lac_mldoob_stab    !  length scale ratios, MLD is boundary, EK is Ekman,
  real :: lac_ekoob_stab     !  and OB is Obukhov, the "o" in the name is for division.
  real :: lac_mldoob_un      !  Stab/un are for stable (pos) and unstable (neg) Obukhov depths
  real :: lac_ekoob_un       !   ...
  real :: ladepthratio=0.04  ! The ratio of OBL depth to average Stokes drift over
  real :: max_enhance_m = 5. ! The maximum allowed LT enhancement to the mixing.
  real :: cnv_mst_fac        ! Factor to reduce mstar when statically unstable.
  type(time_type), pointer :: time ! A pointer to the ocean model's clock.
  integer :: lt_enhance_form = 0 ! Option for Langmuir enhancement form
  integer :: mstar_mode = 0  ! An integer to determine which formula is used to
                             !  set mstar
  integer :: const_mstar=0,mld_o_obukhov=1,ekman_o_obukhov=2
  logical :: mstar_flatcap=.true. !Set false to use asymptotic mstar cap.
  logical :: tke_diagnostics = .false.
  logical :: use_la_windsea = .false.
  logical :: orig_pe_calc = .true.
  logical :: use_mld_iteration=.false. ! False to use old ePBL method.
  logical :: orig_mld_iteration=.false. ! False to use old MLD value
  logical :: mld_iteration_guess=.false. ! False to default to guessing half the
                                         ! ocean depth for the iteration.
  logical :: mixing_diagnostics = .false. ! Will be true when outputing mixing
                                          !  length and velocity scale
  type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                             ! timing of diagnostic output.

! These are terms in the mixed layer TKE budget, all in J m-2 = kg s-2.
  real, allocatable, dimension(:,:) :: &
    diag_tke_wind, &   ! The wind source of TKE.
    diag_tke_mke, &    ! The resolved KE source of TKE.
    diag_tke_conv, &   ! The convective source of TKE.
    diag_tke_forcing, & ! The TKE sink required to mix surface
                       ! penetrating shortwave heating.
    diag_tke_mech_decay, & ! The decay of mechanical TKE.
    diag_tke_conv_decay, & ! The decay of convective TKE.
    diag_tke_mixing,&  ! The work done by TKE to deepen
                       ! the mixed layer.
    ! Additional output parameters also 2d
    ml_depth, &        ! The mixed layer depth in m.
    ml_depth2, &       ! The mixed layer depth in m.
    enhance_m, &       ! The enhancement to the turbulent velocity scale (non-dim)
    mstar_mix, &       ! Mstar used in EPBL
    mld_ekman, &       ! MLD over Ekman length
    mld_obukhov, &     ! MLD over Obukhov length
    ekman_obukhov, &   ! Ekman over Obukhov length
    la, &              ! Langmuir number
    la_mod             ! Modified Langmuir number

  real, allocatable, dimension(:,:,:) :: &
    velocity_scale, & ! The velocity scale used in getting Kd
    mixing_length     ! The length scale used in getting Kd
  integer :: id_ml_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_hsfc_used = -1
  integer :: id_mixing_length = -1, id_velocity_scale = -1
  integer :: id_osbl = -1, id_lt_enhancement = -1, id_mstar_mix = -1
  integer :: id_mld_ekman, id_mld_obukhov, id_ekman_obukhov
  integer :: id_la, id_la_mod
end type energetic_pbl_cs

integer :: num_msg = 0, max_msg = 2

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, CS, &
                         dSV_dT, dSV_dS, TKE_forced, Buoy_Flux, dt_diag, last_call, &
                         dT_expected, dS_expected )
  type(ocean_grid_type),   intent(inout) :: G      !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: GV     !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(inout) :: h_3d   !< Layer thickness, in m or kg m-2.
                                                   !! (Intent in/out) The units of h are referred
                                                   !! to as H below.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: u_3d   !< Zonal velocities interpolated to h points,
                                                   !! m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: v_3d   !< Zonal velocities interpolated to h points,
                                                   !! 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,
                                                   !! in m3 kg-1 K-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                           intent(in)    :: dSV_dS !< The partial derivative of in-situ specific
                                                   !! volume with salinity, in 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
                                                   !! through each layer, in 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, in s.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                           intent(out)   :: Kd_int !< The diagnosed diffusivities at interfaces,
                                                   !! in 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. in m2/s3.
  real,          optional, intent(in)    :: dt_diag   !< The diagnostic time step, which may be less
                                                   !! than dt if there are two callse to
                                                   !! mixedlayer, in 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, dS_expected

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

! Arguments: h_3d - Layer thickness, in m or kg m-2. (Intent in/out)
!                   The units of h are referred to as H below.
!  (in)      u_3d - Zonal velocities interpolated to h points, m s-1.
!  (in)      v_3d - Zonal velocities interpolated to h points, m s-1.
!  (in/out)  tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      dt - Time increment, in s.
!  (out)     Kd_int - The diagnosed diffusivities at interfaces, in m2 s-1.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 mixedlayer_init.
!  (in)      dSV_dT - The partial derivative of in-situ specific volume with
!                     potential temperature, in m3 kg-1 K-1.
!  (in)      dSV_dS - The partial derivative of in-situ specific volume with
!                     salinity, in m3 kg-1 ppt-1.
!  (in)      TKE_forced - The forcing requirements to homogenize the forcing
!                 that has been applied to each layer through each layer, in J m-2.
!  (in)      Buoy_Flux - The surface buoyancy flux. in m2/s3.
!  (in,opt)  dt_diag - The diagnostic time step, which may be less than dt
!                      if there are two callse to mixedlayer, in s.
!  (in,opt)  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),SZK_(GV)) :: &
    h, &            !   The layer thickness, in H (usually m or kg m-2).
    T, &            !   The layer temperatures, in deg C.
    S, &            !   The layer salinities, in psu.
    u, &            !   The zonal velocity, in m s-1.
    v               !   The meridional velocity, in m s-1.
  real, dimension(SZI_(G),SZK_(GV)+1) :: &
    Kd, &           ! The diapycnal diffusivity, in m2 s-1.
    pres, &         ! Interface pressures in Pa.
    hb_hs           ! The distance from the bottom over the thickness of the
                    ! water column, nondim.
  real, dimension(SZI_(G)) :: &
    mech_TKE, &     !   The mechanically generated turbulent kinetic energy
                    ! available for mixing over a time step, in J m-2 = kg s-2.
    conv_perel, &   ! The potential energy that has been convectively released
                    ! during this timestep, in J m-2 = kg s-2. A portion nstar_FC
                    ! of conv_PErel is available to drive mixing.
    htot, &         !   The total depth of the layers above an interface, in H.
    uhtot, &        !   The depth integrated zonal and meridional velocities in the
    vhtot, &        ! layers above, in H m s-1.
    mech_tke_top, & !    The value of mech_TKE at the top of the column, in J m-2.
    conv_perel_top, & !  The value of conv_PErel at the top of the column, in J m-2.

    idecay_len_tke, &  ! The inverse of a turbulence decay length scale, in H-1.
    h_sum, &        ! The total thickness of the water column, in H.
    absf            ! The absolute value of f, in s-1.


  real, dimension(SZI_(G),SZK_(GV)) :: &
    dT_to_dColHt, & ! Partial derivatives of the total column height with the temperature
    dS_to_dColHt, & ! and salinity changes within a layer, in m K-1 and m ppt-1.
    dT_to_dPE, &    ! Partial derivatives of column potential energy with the temperature
    dS_to_dPE, &    ! and salinity changes within a layer, in J m-2 K-1 and J m-2 ppt-1.
    dT_to_dColHt_a, & ! Partial derivatives of the total column height with the temperature
    dS_to_dColHt_a, & ! and salinity changes within a layer, including the implicit effects
                    ! of mixing with layers higher in the water colun, in m K-1 and m ppt-1.
    dt_to_dpe_a, &  ! Partial derivatives of column potential energy with the temperature
    ds_to_dpe_a     ! and salinity changes within a layer, including the implicit effects
                    ! of mixing with layers higher in the water column, in
                    ! units of J m-2 K-1 and J m-2 ppt-1.
  real, dimension(SZK_(GV)) :: &
    T0, S0, &       ! Initial values of T and S in the column, in K and ppt.
    Te, Se, &       ! Estimated final values of T and S in the column, in K and ppt.
    c1, &           ! c1 is used by the tridiagonal solver, ND.
    dTe, dSe        ! Running (1-way) estimates of temperature and salinity change.
  real, dimension(SZK_(GV)) :: &
    Th_a, &         ! An effective temperature times a thickness in the layer above,
                    ! including implicit mixing effects with other yet higher layers, in K H.
    sh_a, &         ! An effective salinity times a thickness in the layer above,
                    ! including implicit mixing effects with other yet higher layers, in K H.
    th_b, &         ! An effective temperature times a thickness in the layer below,
                    ! including implicit mixing effects with other yet lower layers, in K H.
    sh_b            ! An effective salinity times a thickness in the layer below,
                    ! including implicit mixing effects with other yet lower layers, in K H.
  real, dimension(SZI_(G)) :: &
    hp_a            ! An effective pivot thickness of the layer including the effects
                    ! of coupling with layers above, in H.  This is the first term
                    ! in the denominator of b1 in a downward-oriented tridiagonal solver.
  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, in H (m or kg m-2).
  real :: b1        ! b1 is inverse of the pivot used by the tridiagonal solver, in H-1.
  real :: h_neglect ! A thickness that is so small it is usually lost
                    ! in roundoff and can be neglected, in H.
  real :: dMass     ! The mass per unit area within a layer, in kg m-2.
  real :: dPres     ! The hydrostatic pressure change across a layer, in Pa.
  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, in J m-2 = kg s-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, in H-1.

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

  real :: C1_3      ! = 1/3.
  real :: vonKar    ! The vonKarman constant.
  real :: I_dtrho   ! 1.0 / (dt * Rho0) in m3 kg-1 s-1.  This is
                    ! used convert TKE back into ustar^3.
  real :: U_star    ! The surface friction velocity, in m s-1.
  real :: U_Star_Mean ! The surface friction without gustiness in m s-1.
  real :: vstar     ! An in-situ turbulent velocity, in m s-1.
  real :: Enhance_M ! An enhancement factor for vstar, based here on Langmuir impact.
  real :: LA        ! The Langmuir number (non-dim)
  real :: LAmod     ! A modified Langmuir number accounting for other parameters.
  real :: hbs_here  ! The local minimum of hb_hs and MixLen_shape, times a
                    ! conversion factor from H to M, in m H-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, in J m-2.
  real :: TKE_here  ! The total TKE at this point in the algorithm, in J m-2.
  real :: dT_km1_t2 ! A diffusivity-independent term related to the temperature
                    ! change in the layer above the interface, in K.
  real :: dS_km1_t2 ! A diffusivity-independent term related to the salinity
                    ! change in the layer above the interface, in ppt.
  real :: dTe_term  ! A diffusivity-independent term related to the temperature
                    ! change in the layer below the interface, in K H.
  real :: dSe_term  ! A diffusivity-independent term related to the salinity
                    ! change in the layer above the interface, in ppt H.
  real :: dTe_t2    ! A part of dTe_term, in K H.
  real :: dSe_t2    ! A part of dSe_term, in ppt H.
  real :: dPE_conv  ! The convective change in column potential energy, in J m-2.
  real :: MKE_src   ! The mean kinetic energy source of TKE due to Kddt_h(K), in J m-2.
  real :: dMKE_src_dK  ! The partial derivative of MKE_src with Kddt_h(K), in J m-2 H-1.
  real :: Kd_guess0, PE_chg_g0, dPEa_dKd_g0, Kddt_h_g0
  real :: PE_chg_max   ! The maximum PE change for very large values of Kddt_h(K).
  real :: dPEc_dKd_Kd0 ! The partial derivative of PE change with Kddt_h(K)
                       ! for very small values of Kddt_h(K), in J m-2 H-1.
  real :: PE_chg    ! The change in potential energy due to mixing at an
                    ! interface, in 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), in J m-2.
  real :: dPEc_dKd  ! The partial derivative of PE_chg with Kddt_h(K), in J m-2 H-1.
  real :: TKE_left_min, TKE_left_max, Kddt_h_max, Kddt_h_min
  real :: Kddt_h_guess ! A guess at the value of Kddt_h(K), in H.
  real :: Kddt_h_next  ! The next guess at the value of Kddt_h(K), in H.
  real :: dKddt_h      ! The change between guesses at Kddt_h(K), in H.
  real :: dKddt_h_Newt ! The change between guesses at Kddt_h(K) with Newton's method, in H.
  real :: Kddt_h_newt  ! The Newton's method next guess for Kddt_h(K), in H.
  real :: exp_kh    ! The nondimensional decay of TKE across a layer, ND.
  logical :: use_Newt  ! Use Newton's method for the next guess at Kddt_h(K).
  logical :: convectively_stable
  logical, dimension(SZI_(G)) :: &
    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 is only used as a diagnostic.
  real :: dt__diag  ! A copy of dt_diag (if present) or dt, in s.
  real :: IdtdR0    !  = 1.0 / (dt__diag * Rho0), in m3 kg-1 s-1.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    Hsfc_used   ! The thickness of the surface region after buffer layer
  logical :: write_diags  ! If true, write out diagnostics with this step.
  logical :: reset_diags  ! If true, zero out the accumulated diagnostics.
                ! detrainment, in units of m.
  ! Local column copies of energy change diagnostics, all in J m-2.
  real :: dTKE_conv, dTKE_forcing, dTKE_mixing
  real :: dTKE_MKE, dTKE_mech_decay, dTKE_conv_decay
  !----------------------------------------------------------------------
  !/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, in m.
  real :: max_MLD, min_MLD ! Iteration bounds, in m, which are adjusted at each step
                           !  - 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 rarely fail and the added
                           !    cost is likely not significant.  Additionally,
                           !    when it fails it does so in a reasonable
                           !    manner giving a usable guess. When it
                           !    does fail, it is due to convection within
                           !    the boundary.  Likely, a new method e.g.
                           !    surface_disconnect, can improve this.
  logical :: FIRST_OBL     ! Flag for computing "found" Mixing layer depth
  logical :: OBL_CONVERGED ! Flag for convergence of MLD
  integer :: OBL_IT        ! Iteration counter
!### These need to be made into run-time parameters.
  integer :: MAX_OBL_IT=20 ! Set maximum number of iterations.  Probably
                           !  best as an input parameter, but then may want
                           !  to use allocatable arrays if storing
                           !  guess/found (as diagnostic); skipping for now.
                           !  In reality, the maximum number of guesses
                           !  needed is set by:
                           !    DEPTH/2^M < DZ
                           !    where M is the number of guesses
                           !    e.g. M=12 for DEPTH=4000m and DZ=1m
  real, dimension(SZK_(GV)+1) :: Vstar_Used, &      ! 1D arrays used to store
                               Mixing_Length_Used   ! Vstar and Mixing_Length
  !/BGR - remaining variables are related to tracking iteration statistics.
  logical :: OBL_IT_STATS=.false. ! Flag for computing OBL iteration statistics
  REAL :: ITguess(20), ITresult(20),ITmax(20),ITmin(20) ! Flag for storing guess/result
                                                        ! should have dim=MAX_OBL_IT
  integer, save :: MAXIT=0   ! Stores maximum number of iterations
  integer, save :: MINIT=1e8 ! Stores minimum number of iterations
  integer, save :: SUMIT=0   ! Stores total iterations (summed over all)
  integer, save :: NUMIT=0   ! Stores number of times iterated
                             !e.g. Average iterations = SUMIT/NUMIT
  integer, save :: CONVERGED!
  integer, save :: NOTCONVERGED!
  !-End BGR iteration parameters-----------------------------------------
  real :: N2_dissipation
  real :: Bf_STABLE ! Buoyancy flux, capped at 0 (negative only)
  real :: Bf_UNSTABLE ! Buoyancy flux, floored at 0 (positive only)
  real :: STAB_SCALE ! Composite of Stabilizing length scales:
                     !  Ekman scale and Monin-Obukhov scale.
  real :: iL_Ekman ! Inverse of Ekman length scale
  real :: iL_Obukhov ! Inverse of Obukhov length scale
  real :: MLD_o_Ekman          ! >
  real :: MLD_o_Obukhov_stab   ! Ratios of length scales where MLD is boundary layer depth
  real :: Ekman_o_Obukhov_stab ! >
  real :: MLD_o_Obukhov_un   ! Ratios of length scales where MLD is boundary layer depth
  real :: Ekman_o_Obukhov_un ! >

  real :: C_MO = 1. ! Constant in STAB_SCALE for Monin-Obukhov
  real :: C_EK = 2. ! Constant in STAB_SCALE for Ekman length
  real :: MLD_over_STAB ! Mixing layer depth divided by STAB_SCALE
  real :: MSTAR_MIX! The value of mstar (Proportionality of TKE to drive mixing to ustar
                    ! cubed) which is computed as a function of latitude, boundary layer depth,
                    ! and the Monin-Obukhov depth.
  real :: MSTAR_LT ! The added mstar contribution due to Langmuir turbulence
  real :: MSTAR_Conv_Adj ! Adjustment made to mstar due to convection reducing mechanical mixing.
  real :: MSTAR_STAB, MSTAR_ROT ! Mstar in each limit, max is used.
  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(SZI_(G),SZK_(GV)) :: &
    mech_TKE_k, conv_PErel_k
  real, dimension(SZK_(GV)) :: nstar_k
  integer, dimension(SZK_(GV)) :: num_itts

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

  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!")
  if (present(dt_expected) .or. present(ds_expected)) debug = .true.

  h_neglect = gv%H_subroundoff

  if(.not.cs%Use_MLD_Iteration) max_obl_it=1
  c1_3 = 1.0 / 3.0
  dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag
  idtdr0 = 1.0 / (dt__diag * gv%Rho0)
  write_diags = .true. ; if (present(last_call)) write_diags = last_call
  max_itt = 20

  h_tt_min = 0.0
  vonkar = 0.41
  mstar_mix=cs%MSTAR!Initialize to mstar
  i_dtrho = 0.0 ; if (dt*gv%Rho0 > 0.0) i_dtrho = 1.0 / (dt*gv%Rho0)

  ! 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
!!OMP parallel default(none) shared(is,ie,js,je,CS)
    if (cs%TKE_diagnostics) then
!!OMP do
      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_forcing(i,j) = 0.0
      enddo ; enddo
    endif
    if (cs%Mixing_Diagnostics) then
      cs%Mixing_Length(:,:,:) = 0.0
      cs%Velocity_Scale(:,:,:) = 0.0
    endif
!!OMP end parallel
  endif


!!OMP parallel do default(none) shared(js,je,nz,is,ie,h_3d,u_3d,v_3d,tv,dt,      &
!!OMP                                  CS,G,GV,fluxes,IdtdR0,                    &
!!OMP                                  TKE_forced,debug,H_neglect,dSV_dT,        &
!!OMP                                  dSV_dS,I_dtrho,C1_3,h_tt_min,vonKar,     &
!!OMP                                  max_itt,Kd_int)                           &
!!OMP                           private(i,j,k,h,u,v,T,S,Kd,mech_TKE_k,conv_PErel_k,    &
!!OMP                                   U_Star,absf,mech_TKE,conv_PErel,nstar_k, &
!!OMP                                   h_sum,I_hs,h_bot,hb_hs,T0,S0,num_itts,   &
!!OMP                                   pres,dMass,dPres,dT_to_dPE,dS_to_dPE,    &
!!OMP                                   dT_to_dColHt,dS_to_dColHt,Kddt_h,hp_a,   &
!!OMP                                   Th_a,Sh_a,Th_b,Sh_b,dT_to_dPE_a,htot,    &
!!OMP                                   dT_to_dColHt_a,dS_to_dColHt_a,uhtot,vhtot, &
!!OMP                                   Idecay_len_TKE,exp_kh,nstar_FC,tot_TKE,  &
!!OMP                                   TKE_reduc,dTe_t2,dSe_t2,dTe,dSe,dt_h,    &
!!OMP                                   Convectively_stable,sfc_disconnect,b1,   &
!!OMP                                   c1,dT_km1_t2,dS_km1_t2,dTe_term,         &
!!OMP                                   dSe_term,MKE2_Hharm,vstar,h_tt,h_rsum,   &
!!OMP                                   Kd_guess0,Kddt_h_g0,dPEc_dKd_Kd0,        &
!!OMP                                   PE_chg_max,dPEa_dKd_g0,PE_chg_g0,        &
!!OMP                                   MKE_src,dPE_conv,Kddt_h_max,Kddt_h_min,  &
!!OMP                                   dTKE_conv, dTKE_forcing, dTKE_mixing,    &
!!OMP                                   dTKE_MKE,dTKE_mech_decay,dTKE_conv_decay,&
!!OMP                                   TKE_left_max,TKE_left_min,Kddt_h_guess,  &
!!OMP                                   TKE_left_itt,dPEa_dKd_itt,PE_chg_itt,    &
!!OMP                                   MKE_src_itt,Kddt_h_itt,dPEc_dKd,PE_chg,  &
!!OMP                                   dMKE_src_dK,TKE_left,use_Newt,           &
!!OMP                                   dKddt_h_Newt,Kddt_h_Newt,Kddt_h_next,    &
!!OMP                                   dKddt_h,Te,Se,Hsfc_used,dS_to_dPE_a,     &
!!OMP                                   dMKE_max,sfc_connected,TKE_here)
  do j=js,je
    ! Copy the thicknesses and other fields to 2-d arrays.
    do k=1,nz ; do i=is,ie
      h(i,k) = h_3d(i,j,k) + h_neglect ; u(i,k) = u_3d(i,j,k) ; v(i,k) = v_3d(i,j,k)
      t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
      kd(i,k) = 0.0
    enddo ; enddo
    do i=is,ie
       cs%ML_depth(i,j) = h(i,1)*gv%H_to_m
       !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m
       sfc_connected(i) = .true.
    enddo

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

    !   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


      u_star = fluxes%ustar(i,j)
      taux2 = 0.
      if ((g%mask2dCu(i-1,j) + g%mask2dCu(i,j)) > 0) &
        taux2 = (g%mask2dCu(i-1,j)*fluxes%taux(i-1,j)**2 + &
                 g%mask2dCu(i,j)*fluxes%taux(i,j)**2) / (g%mask2dCu(i-1,j) + g%mask2dCu(i,j))
      tauy2 = 0.0
      if ((g%mask2dCv(i,j-1) + g%mask2dCv(i,j)) > 0) &
        tauy2 = (g%mask2dCv(i,j-1)*fluxes%tauy(i,j-1)**2 + &
                 g%mask2dCv(i,j)*fluxes%tauy(i,j)**2) / (g%mask2dCv(i,j-1) + g%mask2dCv(i,j))
      u_star_mean = sqrt(sqrt(taux2 + tauy2)/gv%rho0)
      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
      ! Computing Bf w/ limiters.
      bf_stable = max(0.0,buoy_flux(i,j)) ! Positive for stable
      bf_unstable = min(0.0,buoy_flux(i,j)) ! Negative for unstable
      if (cs%omega_frac >= 1.0) then ; absf(i) = 2.0*cs%omega
      else
        absf(i) = 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(i) = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf(i)**2)
      endif
      ! Computing stability scale which correlates with TKE for mixing, where
      ! TKE for mixing = TKE production minus TKE dissipation
      stab_scale = u_star**2 / ( vonkar * ( c_mo * bf_stable/u_star -  c_ek * u_star * absf(i)))
      ! Inverse of Ekman and Obukhov
      il_ekman   = absf(i)/u_star
      il_obukhov = buoy_flux(i,j)*vonkar/u_star**3

      if (cs%Mstar_Mode.eq.cs%CONST_MSTAR) then
        mech_tke(i) = (dt*cs%mstar*gv%Rho0)*((u_star**3))
        conv_perel(i) = 0.0

        if (cs%TKE_diagnostics) then
          cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + mech_tke(i) * idtdr0
          if (tke_forced(i,j,1) <= 0.0) then
            cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + &
                 max(-mech_tke(i), tke_forced(i,j,1)) * idtdr0
            ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + &
            !     min(0.0, TKE_forced(i,j,1) + mech_TKE(i)) * IdtdR0
          else
            cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + cs%nstar*tke_forced(i,j,1) * idtdr0
          endif
        endif

        if (tke_forced(i,j,1) <= 0.0) then
          mech_tke(i) = mech_tke(i) + tke_forced(i,j,1)
          if (mech_tke(i) < 0.0) mech_tke(i) = 0.0
        else
          conv_perel(i) = conv_perel(i) + tke_forced(i,j,1)
        endif

      endif

!    endif ; enddo

!    do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then

      h_sum(i) = h_neglect ; do k=1,nz ; h_sum(i) = h_sum(i) + h(i,k) ; enddo
      i_hs = 0.0 ; if (h_sum(i) > 0.0) i_hs = 1.0 / h_sum(i)

      h_bot = 0.0 ; hb_hs(i,nz+1) = 0.0
      do k=nz,1,-1
        h_bot = h_bot + h(i,k)
        hb_hs(i,k) = h_bot * i_hs
      enddo

      pres(i,1) = 0.0
      do k=1,nz
        dmass = gv%H_to_kg_m2 * h(i,k)
        dpres = gv%g_Earth * dmass
        dt_to_dpe(i,k) = (dmass * (pres(i,k) + 0.5*dpres)) * dsv_dt(i,j,k)
        ds_to_dpe(i,k) = (dmass * (pres(i,k) + 0.5*dpres)) * dsv_ds(i,j,k)
        dt_to_dcolht(i,k) = dmass * dsv_dt(i,j,k)
        ds_to_dcolht(i,k) = dmass * dsv_ds(i,j,k)

        pres(i,k+1) = pres(i,k) + dpres
      enddo

!    endif ; enddo

    ! Note the outer i-loop and inner k-loop loop order!!!
!    do i=is,ie ; if (G%mask2dT(i,j) > 0.5) then
      do k=1,nz ; t0(k) = t(i,k) ; s0(k) = s(i,k) ; enddo

      ! Store the initial mechanical TKE and convectively released PE to
      ! enable multiple iterations.
      mech_tke_top(i) = mech_tke(i) ; conv_perel_top(i) = conv_perel(i)

      !/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(i,k)*gv%H_to_m ; enddo
      min_mld = 0.0 !min_MLD will initialize as 0.
      !/BGR: May add user-input bounds for max/min MLD

      !/BGR: Add MLD_guess based on stored previous value.
      !      note that this is different from ML_Depth already
      !      computed by EPBL, need to figure out why.
      if (cs%MLD_iteration_guess .and. cs%ML_Depth2(i,j) > 1.) then
        !If prev value is present use for guess.
        mld_guess=cs%ML_Depth2(i,j)
      else
        !Otherwise guess middle of water column (or stab_scale if smaller).
        mld_guess = 0.5 * (min_mld+max_mld)
      endif

      ! Iterate up to MAX_OBL_IT times to determine a converged EPBL depth.
      obl_converged = .false.
      ! Initialize ENHANCE_M to 1 and mstar_lt to 0
      enhance_m=1.e0
      mstar_lt = 0.0
      do obl_it=1,max_obl_it ; if (.not. obl_converged) then

        ! Reset ML_depth
        cs%ML_depth(i,j) = h(i,1)*gv%H_to_m
        !CS%ML_depth2(i,j) = h(i,1)*GV%H_to_m

        sfc_connected(i) = .true.

        if (cs%Mstar_Mode.gt.0) then
        ! Note the value of mech_TKE(i) now must be iterated over, so it is moved here
        ! First solve for the TKE to PE length scale
          if (cs%MSTAR_MODE.eq.cs%MLD_o_OBUKHOV) then
            mld_over_stab = mld_guess / stab_scale - cs%MSTAR_XINT
            if ((mld_over_stab) .le. 0.0) then
              !Asymptote to 0 as MLD_over_Stab -> -infinity (always)
              mstar_mix = (cs%MSTAR_B*(mld_over_stab)+cs%MSTAR_A)**(cs%MSTAR_N)
            else
              if (cs%MSTAR_CAP>=0.) then
                if (cs%MSTAR_FLATCAP .OR. (mld_over_stab .le.cs%MSTAR_XINT_UP)) then
                !If using flat cap (or if using asymptotic cap
                !   but within linear regime we can make use of same code)
                  mstar_mix = min(cs%MSTAR_CAP, &
                     cs%MSTAR_SLOPE*(mld_over_stab)+cs%MSTAR_AT_XINT)
                else
                  !Asymptote to MSTAR_CAP as MLD_over_Stab -> infinity
                  mstar_mix = cs%MSTAR_CAP - &
                       (cs%MSTAR_B2*(mld_over_stab-cs%MSTAR_XINT_UP)&
                       +cs%MSTAR_A2)**(cs%MSTAR_N)
                endif
              else
                !No cap if negative cap value given.
                mstar_mix = cs%MSTAR_SLOPE*(mld_over_stab)+cs%MSTAR_AT_XINT
              endif
            endif
          elseif (cs%MSTAR_MODE.eq.cs%EKMAN_o_OBUKHOV) then
            mstar_stab = cs%MSTAR_COEF*sqrt(bf_stable/u_star**2/(absf(i)+1.e-10))
            mstar_rot =  cs%C_EK*log(max(1.,u_star/(absf(i)+1.e-10)/mld_guess))
            if ( cs%MSTAR_CAP.le.0.0) then !No cap.
              mstar_mix = max(mstar_stab,& ! 1st term if balance of rotation and stabilizing
                                    ! the balance is f(L_Ekman,L_Obukhov)
                              min(& ! 2nd term for forced stratification limited
                                 1.25,& !.5/von Karman (Obukhov limit)
                                 ! 3rd term for rotation (Ekman length) limited
                                 mstar_rot))
            else
              mstar_mix = min( & ! Sets a cap. The cap should be large and just
                               !  meant to be a safety net.
                           cs%MSTAR_CAP, &
                           max(mstar_stab,& ! 1st term if balance of rotation and stabilizing
                                 ! the balance is f(L_Ekman,L_Obukhov)
                              min(& ! 2nd term for forced stratification limited
                                 1.25,& !.5/von Karman (Obukhov limit)
                                 ! 3rd term for rotation (Ekman length) limited
                                 mstar_rot)))
            endif!cap for mstar_mode==2
          endif!mstar_mode==1 or ==2
          ! Adjustment for unstable buoyancy flux.
          !  Convection reduces mechanical mixing because there
          !  is less density gradient to mix. (Statically unstable near surface)
          mstar_conv_adj = 1. - cs%CNV_MST_FAC * (-bf_unstable+1.e-10) / &
                                 ( (-bf_unstable+1.e-10)+                 &
                                   2. *mstar_mix *u_star**3 / mld_guess )
          if (cs%Use_LA_windsea) then
            ! 1. Get LA
            call get_la_windsea( u_star_mean, mld_guess*cs%LaDepthRatio, gv, la)
            ! 2. Get parameters for modified LA
            mld_o_ekman = abs(mld_guess*il_ekman)
            mld_o_obukhov_stab = abs(max(0.,mld_guess*il_obukhov))
            mld_o_obukhov_un = abs(min(0.,mld_guess*il_obukhov))
            ekman_o_obukhov_stab = abs(max(0.,il_obukhov/(il_ekman+1.e-10)))
            ekman_o_obukhov_un = abs(min(0.,il_obukhov/(il_ekman+1.e-10)))
            ! 3. 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.
            lamod = la * (1.0 + max(-0.5,cs%LaC_MLDoEK * mld_o_ekman) +   &
                 cs%LaC_EKoOB_stab * ekman_o_obukhov_stab + &
                 cs%LaC_EKoOB_un * ekman_o_obukhov_un + &
                 cs%LaC_MLDoOB_stab * mld_o_obukhov_stab  + &
                 cs%LaC_MLDoOB_un * mld_o_obukhov_un )
            if (cs%LT_Enhance_Form==1) then
              !Original w'/ust scaling w/ Van Roekel's scaling
              enhance_m = (1+(1.4*la)**(-2)+(5.4*la)**(-4))**(1.5)
            elseif (cs%LT_Enhance_Form==2) then
              ! Enhancement is multiplied (added mst_lt set to 0)
              enhance_m = min(cs%Max_Enhance_M,(1.+cs%LT_ENHANCE_COEF*lamod**cs%LT_ENHANCE_EXP))
              mstar_lt = 0.0
            elseif (cs%LT_ENHANCE_Form == 3) then
              ! or Enhancement is additive (multiplied enhance_m set to 1)
              mstar_lt = cs%LT_ENHANCE_COEF * lamod**cs%LT_ENHANCE_EXP
              enhance_m = 1.0
            endif
          endif
          !Reset mech_tke and conv_perel values (based on new mstar)
          mech_tke(i) = ( mstar_mix * mstar_conv_adj * enhance_m + mstar_lt) * (dt*gv%Rho0*u_star**3)
          conv_perel(i) = 0.0
          if (cs%TKE_diagnostics) then
            cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + mech_tke(i) * idtdr0
            if (tke_forced(i,j,1) <= 0.0) then
              cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + &
                   max(-mech_tke(i), tke_forced(i,j,1)) * idtdr0
              ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + &
              !     min(0.0, TKE_forced(i,j,1) + mech_TKE(i)) * IdtdR0
            else
              cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + cs%nstar*tke_forced(i,j,1) * idtdr0
            endif
          endif

          if (tke_forced(i,j,1) <= 0.0) then
            mech_tke(i) = mech_tke(i) + tke_forced(i,j,1)
            if (mech_tke(i) < 0.0) mech_tke(i) = 0.0
          else
            conv_perel(i) = conv_perel(i) + tke_forced(i,j,1)
          endif
        else
          mech_tke(i) = mech_tke_top(i)*enhance_m ; conv_perel(i) = conv_perel_top(i)
        endif

        if (cs%TKE_diagnostics) then
          dtke_conv = 0.0 ; dtke_forcing = 0.0 ; dtke_mixing = 0.0
          dtke_mke = 0.0 ; dtke_mech_decay = 0.0 ; dtke_conv_decay = 0.0
        endif

        ! Store in 1D arrays cleared out each iteration.  Only write in
        !  3D arrays after convergence.
        do k=1,nz
          vstar_used(k) = 0.0 ; mixing_length_used(k) = 0.0
        enddo
        if (.not.cs%Use_MLD_Iteration) obl_converged = .true.

        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(i,k-1)
            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(i,1) = 0.0 ; kddt_h(1) = 0.0
        hp_a(i) = h(i,1)
        dt_to_dpe_a(i,1) = dt_to_dpe(i,1) ; dt_to_dcolht_a(i,1) = dt_to_dcolht(i,1)
        ds_to_dpe_a(i,1) = ds_to_dpe(i,1) ; ds_to_dcolht_a(i,1) = ds_to_dcolht(i,1)

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

        if (debug) then
          mech_tke_k(i,1) = mech_tke(i) ; conv_perel_k(i,1) = conv_perel(i)
          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(i) = (cs%TKE_decay * absf(i) / u_star) * gv%H_to_m
          exp_kh = 1.0
          if (idecay_len_tke(i) > 0.0) exp_kh = exp(-h(i,k-1)*idecay_len_tke(i))
          if (cs%TKE_diagnostics) &
            dtke_mech_decay = dtke_mech_decay + (exp_kh-1.0) * mech_tke(i) * idtdr0
          mech_tke(i) = mech_tke(i) * exp_kh

          !   Accumulate any convectively released potential energy to contribute
          ! to wstar and to drive penetrating convection.
          if (tke_forced(i,j,k) > 0.0) then
            conv_perel(i) = conv_perel(i) + tke_forced(i,j,k)
            if (cs%TKE_diagnostics) &
              dtke_forcing = dtke_forcing + cs%nstar*tke_forced(i,j,k) * idtdr0
          endif

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

          !  Determine the total energy
          nstar_fc = cs%nstar
          if (cs%nstar * conv_perel(i) > 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(i))**3 / conv_PErel(i))
            nstar_fc = cs%nstar * conv_perel(i) / (conv_perel(i) + 0.2 * &
                            sqrt(0.5 * dt * gv%Rho0 * (absf(i)*(htot(i)*gv%H_to_m))**3 * conv_perel(i)))
          endif
          if (debug) nstar_k(k) = nstar_fc

          tot_tke = mech_tke(i) + nstar_fc * conv_perel(i)

          !   For each interior interface, first discard the TKE to account for
          ! mixing of shortwave radiation through the next denser cell.
          if (tke_forced(i,j,k) < 0.0) then
            if (tke_forced(i,j,k) + tot_tke < 0.0) then
              ! The shortwave requirements deplete all the energy in this layer.
              if (cs%TKE_diagnostics) then
                dtke_mixing = dtke_mixing + tot_tke * idtdr0
                dtke_forcing = dtke_forcing - tot_tke * idtdr0
                ! dTKE_unbalanced_forcing = dTKE_unbalanced_forcing + &
                !     (TKE_forced(i,j,k) + tot_TKE) * IdtdR0
                dtke_conv_decay = dtke_conv_decay + &
                        (cs%nstar-nstar_fc) * conv_perel(i) * idtdr0
              endif
              tot_tke = 0.0 ; mech_tke(i) = 0.0 ; conv_perel(i) = 0.0
            else
              ! Reduce the mechanical and convective TKE proportionately.
              tke_reduc = (tot_tke + tke_forced(i,j,k)) / tot_tke
              if (cs%TKE_diagnostics) then
                dtke_mixing = dtke_mixing - tke_forced(i,j,k) * idtdr0
                dtke_forcing = dtke_forcing + tke_forced(i,j,k) * idtdr0
                dtke_conv_decay = dtke_conv_decay + &
                    (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel(i) * idtdr0
              endif
              tot_tke = tke_reduc*tot_tke   ! = tot_TKE + TKE_forced(i,j,k)
              mech_tke(i) = tke_reduc*mech_tke(i)
              conv_perel(i) = tke_reduc*conv_perel(i)
            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%m_to_H**2*dt) / max(0.5*(h(i,k-1)+h(i,k)), 1e-15*h_sum(i))

          !   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(i,k) + dt_to_dcolht(i,k-1) ) * (t0(k-1)-t0(k)) + &
              (ds_to_dcolht(i,k) + ds_to_dcolht(i,k-1) ) * (s0(k-1)-s0(k)) ) )

          if ((mech_tke(i) + conv_perel(i)) <= 0.0 .and. convectively_stable) then
            ! Energy is already exhausted, so set Kd = 0 and cycle or exit?
            tot_tke = 0.0 ; mech_tke(i) = 0.0 ; conv_perel(i) = 0.0
            kd(i,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(i)
            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(i) = h(i,k)
            dt_to_dpe_a(i,k) = dt_to_dpe(i,k)
            ds_to_dpe_a(i,k) = ds_to_dpe(i,k)
            dt_to_dcolht_a(i,k) = dt_to_dcolht(i,k)
            ds_to_dcolht_a(i,k) = ds_to_dcolht(i,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(i)) * ((t0(k-2) - t0(k-1)) + dte(k-2))
                ds_km1_t2 = (s0(k)-s0(k-1)) - &
                      (kddt_h(k-1) / hp_a(i)) * ((s0(k-2) - s0(k-1)) + dse(k-2))
              endif
              dte_term = dte_t2 + hp_a(i) * (t0(k-1)-t0(k))
              dse_term = dse_t2 + hp_a(i) * (s0(k-1)-s0(k))
            else
              if (k<=2) then
                th_a(k-1) = h(i,k-1) * t0(k-1) ; sh_a(k-1) = h(i,k-1) * s0(k-1)
              else
                th_a(k-1) = h(i,k-1) * t0(k-1) + kddt_h(k-1) * te(k-2)
                sh_a(k-1) = h(i,k-1) * s0(k-1) + kddt_h(k-1) * se(k-2)
              endif
              th_b(k) = h(i,k) * t0(k) ; sh_b(k) = h(i,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(i)*h(i,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 = (gv%H_to_kg_m2 * cs%MKE_to_TKE_effic) * 0.5 * &
                  (h(i,k) / ((htot(i) + h(i,k))*htot(i))) * &
                  ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)
              ! A fraction (1-exp(Kddt_h*MKE2_Hharm)) of this energy would be
              ! extracted by mixing with a finite viscosity.
              mke2_hharm = (htot(i) + h(i,k) + 2.0*h_neglect) / &
                           ((htot(i)+h_neglect) * (h(i,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(i) + h_tt_min
            tke_here = mech_tke(i) + cs%wstar_ustar_coef*conv_perel(i)
            if (tke_here > 0.0) then
              vstar = cs%vstar_scale_fac * (i_dtrho*tke_here)**c1_3
              hbs_here = gv%H_to_m * min(hb_hs(i,k), mixlen_shape(k))
              mixing_length_used(k) = max(cs%min_mix_len,((h_tt*hbs_here)*vstar) / &
                  ((cs%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar))
              !Note setting Kd_guess0 to Mixing_Length_Used(K) here will
              ! change the answers.  Therefore, skipping that.
              if (.not.cs%Use_MLD_Iteration) then
                 kd_guess0 = vstar * vonkar *  ((h_tt*hbs_here)*vstar) / &
                  ((cs%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar)
              else
                 kd_guess0 = vstar * vonkar * mixing_length_used(k)
              endif
            else
              vstar = 0.0 ; kd_guess0 = 0.0
            endif
            vstar_used(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(i,k), hp_a(i), dte_term, dse_term, &
                       dt_km1_t2, ds_km1_t2, dt_to_dpe(i,k), ds_to_dpe(i,k), &
                       dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), &
                       pres(i,k), dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                       dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                       pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_g0, dpe_max=pe_chg_max, &
                       dpec_dkd_0=dpec_dkd_kd0 )
            else
              call find_pe_chg(0.0, kddt_h_g0, hp_a(i), h(i,k), &
                         th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
                         dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), dt_to_dpe(i,k), ds_to_dpe(i,k), &
                         pres(i,k), dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                         dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                         pe_chg=pe_chg_g0, dpec_dkd=dpea_dkd_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))

            if (pe_chg_g0 .gt. 0.0) then
              !Negative buoyancy (increases PE)
              n2_dissipation = 1.+cs%N2_DISSIPATION_SCALE_NEG
            else
              !Positive buoyancy (decreases PE)
              n2_dissipation = 1.+cs%N2_DISSIPATION_SCALE_POS
            endif

            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(i) + cs%wstar_ustar_coef*(conv_perel(i)-pe_chg_max)
                if (tke_here > 0.0) then
                  vstar = cs%vstar_scale_fac * (i_dtrho*tke_here)**c1_3
                  hbs_here = gv%H_to_m * min(hb_hs(i,k), mixlen_shape(k))
                  mixing_length_used(k) = max(cs%min_mix_len,((h_tt*hbs_here)*vstar) / &
                      ((cs%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar))
                  if (.not.cs%Use_MLD_Iteration) then
                  ! Note again (as prev) that using Mixing_Length_Used here
                  !  instead of redoing the computation will change answers...
                     kd(i,k) = vstar * vonkar *  ((h_tt*hbs_here)*vstar) / &
                          ((cs%Ekman_scale_coef * absf(i)) * (h_tt*hbs_here) + vstar)
                  else
                     kd(i,k) = vstar * vonkar * mixing_length_used(k)
                  endif
                else
                  vstar = 0.0 ; kd(i,k) = 0.0
                endif
                vstar_used(k) = vstar

                if (cs%orig_PE_calc) then
                  call find_pe_chg_orig(kd(i,k)*dt_h, h(i,k), hp_a(i), dte_term, dse_term, &
                           dt_km1_t2, ds_km1_t2, dt_to_dpe(i,k), ds_to_dpe(i,k), &
                           dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), &
                           pres(i,k), dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                           dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                           pe_chg=dpe_conv)
                else
                  call find_pe_chg(0.0, kd(i,k)*dt_h, hp_a(i), h(i,k), &
                           th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
                           dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), dt_to_dpe(i,k), ds_to_dpe(i,k), &
                           pres(i,k), dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                           dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                           pe_chg=dpe_conv)
                endif
                ! Should this be iterated to convergence for Kd?
                if (dpe_conv > 0.0) then
                  kd(i,k) = kd_guess0 ; dpe_conv = pe_chg_g0
                else
                  mke_src = dmke_max*(1.0 - exp(-(kd(i,k)*dt_h) * mke2_hharm))
                endif
              else
                ! The energy change does not vary monotonically with Kddt_h.  Find the maximum?
                kd(i,k) = kd_guess0 ; dpe_conv = pe_chg_g0
              endif
              conv_perel(i) = conv_perel(i) - dpe_conv
              mech_tke(i) = mech_tke(i) + mke_src
              if (cs%TKE_diagnostics) then
                dtke_conv = dtke_conv - cs%nstar*dpe_conv * idtdr0
                dtke_mke = dtke_mke + mke_src * idtdr0
              endif
              if (sfc_connected(i)) then
                cs%ML_depth(i,j) = cs%ML_depth(i,j) + gv%H_to_m * h(i,k)
                !CS%ML_depth2(i,j) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k)
              endif

              kddt_h(k) = kd(i,k)*dt_h
            elseif (tot_tke + (mke_src - n2_dissipation*pe_chg_g0) >= 0.0) then
              ! There is energy to support the suggested mixing.  Keep that estimate.
              kd(i,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 - n2_dissipation*pe_chg_g0) &
                                             / tot_tke
              if (cs%TKE_diagnostics) then
                dtke_mixing = dtke_mixing - pe_chg_g0 * idtdr0
                dtke_mke = dtke_mke + mke_src * idtdr0
                dtke_conv_decay = dtke_conv_decay + &
                    (1.0-tke_reduc)*(cs%nstar-nstar_fc) * conv_perel(i) * idtdr0
              endif
              tot_tke = tke_reduc*tot_tke
              mech_tke(i) = tke_reduc*(mech_tke(i) + mke_src)
              conv_perel(i) = tke_reduc*conv_perel(i)
              if (sfc_connected(i)) then
                cs%ML_depth(i,j) = cs%ML_depth(i,j) + gv%H_to_m * h(i,k)
                !CS%ML_depth2(i,J) = CS%ML_depth2(i,J) + GV%H_to_m * h(i,k)
              endif
            elseif (tot_tke == 0.0) then
              ! This can arise if nstar_FC = 0.
              kd(i,k) = 0.0 ; kddt_h(k) = 0.0
              tot_tke = 0.0 ; conv_perel(i) = 0.0 ; mech_tke(i) = 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 - n2_dissipation*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( n2_dissipation*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(i,k), hp_a(i), dte_term, dse_term, &
                           dt_km1_t2, ds_km1_t2, dt_to_dpe(i,k), ds_to_dpe(i,k), &
                           dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), &
                           pres(i,k), dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                           dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                           pe_chg=pe_chg, dpec_dkd=dpec_dkd )
                else
                  call find_pe_chg(0.0, kddt_h_guess, hp_a(i), h(i,k), &
                           th_a(k-1), sh_a(k-1), th_b(k), sh_b(k), &
                           dt_to_dpe_a(i,k-1), ds_to_dpe_a(i,k-1), dt_to_dpe(i,k), ds_to_dpe(i,k), &
                           pres(i,k), dt_to_dcolht_a(i,k-1), ds_to_dcolht_a(i,k-1), &
                           dt_to_dcolht(i,k), ds_to_dcolht(i,k), &
                           pe_chg=dpe_conv)
                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 - n2_dissipation*pe_chg)
                if (debug) then
                  kddt_h_itt(itt) = kddt_h_guess ; mke_src_itt(itt) = mke_src
                  pe_chg_itt(itt) = n2_dissipation*pe_chg
                  tke_left_itt(itt) = tke_left
                  dpea_dkd_itt(itt) = dpec_dkd
                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*n2_dissipation - dmke_src_dk <= 0.0) then
                  use_newt = .false.
                else
                  dkddt_h_newt = tke_left / (dpec_dkd*n2_dissipation - 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
              kd(i,k) = kddt_h_guess / dt_h ; kddt_h(k) = kd(i,k)*dt_h

              ! All TKE should have been consumed.
              if (cs%TKE_diagnostics) then
                dtke_mixing = dtke_mixing - (tot_tke + mke_src) * idtdr0
                dtke_mke = dtke_mke + mke_src * idtdr0
                dtke_conv_decay = dtke_conv_decay + &
                    (cs%nstar-nstar_fc) * conv_perel(i) * idtdr0
              endif

              if (sfc_connected(i)) cs%ML_depth(i,j) = cs%ML_depth(i,j) + &
                   (pe_chg / pe_chg_g0) * gv%H_to_m * h(i,k)
              tot_tke = 0.0 ; mech_tke(i) = 0.0 ; conv_perel(i) = 0.0
              sfc_disconnect = .true.
            endif

            kddt_h(k) = kd(i,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(i) + 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(i) = h(i,k) + (hp_a(i) * b1) * kddt_h(k)
            dt_to_dpe_a(i,k) = dt_to_dpe(i,k) + c1(k)*dt_to_dpe_a(i,k-1)
            ds_to_dpe_a(i,k) = ds_to_dpe(i,k) + c1(k)*ds_to_dpe_a(i,k-1)
            dt_to_dcolht_a(i,k) = dt_to_dcolht(i,k) + c1(k)*dt_to_dcolht_a(i,k-1)
            ds_to_dcolht_a(i,k) = ds_to_dcolht(i,k) + c1(k)*ds_to_dcolht_a(i,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(i) = u(i,k)*h(i,k)
            vhtot(i) = v(i,k)*h(i,k)
            htot(i)  = h(i,k)
            sfc_connected(i) = .false.
          else
            uhtot(i) = uhtot(i) + u(i,k)*h(i,k)
            vhtot(i) = vhtot(i) + v(i,k)*h(i,k)
            htot(i)  = htot(i) + h(i,k)
          endif

          if (debug) then
            if (k==2) then
              te(1) = b1*(h(i,1)*t0(1))
              se(1) = b1*(h(i,1)*s0(1))
            else
              te(k-1) = b1 * (h(i,k-1) * t0(k-1) + kddt_h(k-1) * te(k-2))
              se(k-1) = b1 * (h(i,k-1) * s0(k-1) + kddt_h(k-1) * se(k-2))
            endif
          endif
        enddo
        kd(i,nz+1) = 0.0

        if (debug) then
          ! Complete the tridiagonal solve for Te.
          b1 = 1.0 / hp_a(i)
          te(nz) = b1 * (h(i,nz) * t0(nz) + kddt_h(nz) * te(nz-1))
          se(nz) = b1 * (h(i,nz) * s0(nz) + kddt_h(nz) * se(nz-1))
          do k=nz-1,1,-1
            te(k) = te(k) + c1(k+1)*te(k+1)
            se(k) = se(k) + c1(k+1)*se(k+1)
          enddo
        endif
        if (present(dt_expected)) then
          do k=1,nz ; dt_expected(i,j,k) = te(k) - t0(k) ; enddo
        endif
        if (present(ds_expected)) then
          do k=1,nz ; ds_expected(i,j,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(i,k) * (te(k) - t0(k)) + &
                                     ds_to_dpe(i,k) * (se(k) - s0(k)))
          enddo
          mixing_debug = dpe_debug * idtdr0
        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.
        !/
        itmax(obl_it) = max_mld       ! Track max    }
        itmin(obl_it) = min_mld       ! Track min    } For debug purpose
        itguess(obl_it) = mld_guess   ! Track guess  }
        !/
        mld_found=0.0 ; first_obl=.true.
        if (cs%Orig_MLD_iteration) then
          !This is how the iteration was original conducted
          do k=2,nz
            if (first_obl) then !Breaks when OBL found
              if (vstar_used(k) > 1.e-10 .and. k < nz) then
                mld_found = mld_found+h(i,k-1)*gv%H_to_m
              else
                first_obl = .false.
                if (mld_found-cs%MLD_tol > mld_guess) then
                  min_mld = mld_guess
                elseif ((mld_guess-mld_found) < max(cs%MLD_tol,h(i,k-1)*gv%H_to_m)) then
                  obl_converged = .true.!Break convergence loop
                  if (obl_it_stats) then !Compute iteration statistics
                    maxit = max(maxit,obl_it)
                    minit = min(minit,obl_it)
                    sumit = sumit+obl_it
                    numit = numit+1
                    print*,maxit,minit,sumit/numit
                  endif
                  cs%ML_Depth2(i,j) = mld_guess
                else
                  max_mld = mld_guess !We know this guess was too deep
                endif
              endif
            endif
          enddo
        else
          !New method uses ML_DEPTH as computed in ePBL routine
          mld_found=cs%ML_DEPTH(i,j)
          if (mld_found-cs%MLD_tol > mld_guess) then
            min_mld = mld_guess
          elseif (abs(mld_guess-mld_found) < (cs%MLD_tol)) then
            obl_converged = .true.!Break convergence loop
            if (obl_it_stats) then !Compute iteration statistics
              maxit = max(maxit,obl_it)
              minit = min(minit,obl_it)
              sumit = sumit+obl_it
              numit = numit+1
              print*,maxit,minit,sumit/numit
            endif
            cs%ML_Depth2(i,j) = mld_guess
          else
            max_mld = mld_guess !We know this guess was too deep
          endif
        endif
        ! For next pass, guess average of minimum and maximum values.
        mld_guess = min_mld*0.5 + max_mld*0.5
        itresult(obl_it) = mld_found
      endif ; enddo ! Iteration loop for converged boundary layer thickness.
      if (.not.obl_converged) then
        !/Temp output, warn that EPBL didn't converge
        !/Print guess/found for every iteration step
        !print*,'EPBL MLD DID NOT CONVERGE'
        notconverged=notconverged+1
        !do obl_it=1,max_obl_it
        !   print*,ITmin(obl_it),ITmax(obl_it)
        !   print*,obl_it,ITguess(obl_it),ITresult(obl_it)
        !enddo
        !Activate to print out some output when not converged
        !{
        !print*,'Min/Max: ',ITmin(50),ITmax(50)
        !print*,'Guess/result: ',ITguess(50),ITresult(50)
        !print*,'Stats on CPU: ',CONVERGED,NOTCONVERGED,&
        !     real(NOTCONVERGED)/real(CONVERGED)
        !}
        !stop !Kill if not converged during testing.
      else
        converged=converged+1
      endif

      if (cs%TKE_diagnostics) then
        cs%diag_TKE_MKE(i,j) = cs%diag_TKE_MKE(i,j) + dtke_mke
        cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + dtke_conv
        cs%diag_TKE_forcing(i,j) = cs%diag_TKE_forcing(i,j) + dtke_forcing
        cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) + dtke_mixing
        cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + dtke_mech_decay
        cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + dtke_conv_decay
       ! CS%diag_TKE_unbalanced_forcing(i,j) = CS%diag_TKE_unbalanced_forcing(i,j) + dTKE_unbalanced
      endif
      if (cs%Mixing_Diagnostics) then
        !Write to 3-D for outputing Mixing length and
        !  velocity scale.
        do k=1,nz
          cs%Mixing_Length(i,j,k) = mixing_length_used(k)
          cs%Velocity_Scale(i,j,k) = vstar_used(k)
        enddo
      endif
      if (allocated(cs%Enhance_M)) cs%Enhance_M(i,j) = enhance_m
      if (allocated(cs%mstar_mix)) cs%mstar_mix(i,j) = mstar_mix
      if (allocated(cs%MLD_Obukhov)) cs%MLD_Obukhov(i,j) = (mld_guess*il_obukhov)
      if (allocated(cs%MLD_Ekman)) cs%MLD_Ekman(i,j) = (mld_guess*il_ekman)
      if (allocated(cs%Ekman_Obukhov)) cs%Ekman_Obukhov(i,j) = (il_obukhov/(il_ekman+1.e-10))
      if (allocated(cs%La)) cs%La(i,j) = la
      if (allocated(cs%La_mod)) cs%La_mod(i,j) = lamod
    else
      ! For masked points, Kd_int must still be set (to 0) because it has intent(out).
      do k=1,nz+1
        kd(i,k) = 0.
      enddo
      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!

    if (cs%id_Hsfc_used > 0) then
      do i=is,ie ; hsfc_used(i,j) = h(i,1)*gv%H_to_m ; enddo
      do k=2,nz ; do i=is,ie
        if (kd(i,k) > 0.0) hsfc_used(i,j) = hsfc_used(i,j) + h(i,k)*gv%H_to_m
      enddo ; enddo
    endif

    do k=1,nz+1 ; do i=is,ie
      kd_int(i,j,k) = kd(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_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_Hsfc_used > 0) &
      call post_data(cs%id_Hsfc_used, hsfc_used, 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_OSBL >0) &
      call post_data(cs%id_OSBL, cs%ML_Depth2, cs%diag)
    if (cs%id_LT_Enhancement >0) &
      call post_data(cs%id_LT_Enhancement, cs%Enhance_M, cs%diag)
    if (cs%id_MSTAR_MIX >0) &
      call post_data(cs%id_MSTAR_MIX, cs%MSTAR_MIX, cs%diag)
    if (cs%id_MLD_OBUKHOV >0) &
      call post_data(cs%id_MLD_Obukhov, cs%MLD_OBUKHOV, cs%diag)
    if (cs%id_MLD_EKMAN >0) &
      call post_data(cs%id_MLD_Ekman, cs%MLD_EKMAN, cs%diag)
    if (cs%id_Ekman_Obukhov >0) &
      call post_data(cs%id_Ekman_Obukhov, cs%Ekman_Obukhov, 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)

  endif

end subroutine energetic_pbl

!> 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, 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, 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, in units of H (m or kg-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, in units of H (m or kg-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, in H.
  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, in H.
  real, intent(in)  :: Th_a     !< An effective temperature times a thickness in the layer
                                !! above, including implicit mixing effects with other
                                !! yet higher layers, in K H.
  real, intent(in)  :: Sh_a     !< An effective salinity times a thickness in the layer
                                !! above, including implicit mixing effects with other
                                !! yet higher layers, in K H.
  real, intent(in)  :: Th_b     !< An effective temperature times a thickness in the layer
                                !! below, including implicit mixing effects with other
                                !! yet lower layers, in K H.
  real, intent(in)  :: Sh_b     !< An effective salinity times a thickness in the layer
                                !! below, including implicit mixing effects with other
                                !! yet lower layers, in K H.
  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, in J m-2 K-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, in 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, in J m-2 K-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, in J m-2 ppt-1.
  real, intent(in)  :: pres     !< The hydrostatic interface pressure, which is used to relate
                                !! the changes in column thickness to the energy that is radiated
                                !! as gravity waves and unavailable to drive mixing, in Pa.
  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, in m K-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, in 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, in m K-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, in m ppt-1.

  real, optional, intent(out) :: PE_chg   !< The change in column potential energy from applying
                                          !! Kddt_h at the present interface, in J m-2.
  real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h,
                                          !! in units of J m-2 H-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, in 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, in J m-2 H-1.
  real, optional, intent(out) :: ColHt_cor  !< The correction to PE_chg that is made due to a net
                                            !! change in the column height, in J m-2.

  real :: hps ! The sum of the two effective pivot thicknesses, in H.
  real :: bdt1 ! A product of the two pivot thicknesses plus a diffusive term, in H2.
  real :: dT_c ! The core term in the expressions for the temperature changes, in K H2.
  real :: dS_c ! The core term in the expressions for the salinity changes, in psu H2.
  real :: PEc_core ! The diffusivity-independent core term in the expressions
                   ! for the potential energy changes, J m-3.
  real :: ColHt_core ! The diffusivity-independent core term in the expressions
                     ! for the column height changes, J m-3.
  real :: ColHt_chg  ! The change in the column height, in m.
  real :: y1   ! A local temporary term, in units of H-3 or H-4 in various contexts.

  !   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 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
    pe_chg = pec_core * y1
    colht_chg = colht_core * y1
    if (colht_chg < 0.0) pe_chg = pe_chg - pres * colht_chg
    if (present(colht_cor)) colht_cor = -pres * min(colht_chg, 0.0)
  else if (present(colht_cor)) then
    y1 = dkddt_h / (bdt1 * (bdt1 + dkddt_h * hps))
    colht_cor = -pres * min(colht_core * y1, 0.0)
  endif

  if (present(dpec_dkd)) then
    ! Find the derivative of the potential energy change with dKddt_h.
    y1 = 1.0 / (bdt1 + dkddt_h * hps)**2
    dpec_dkd = pec_core * y1
    colht_chg = colht_core * y1
    if (colht_chg < 0.0) dpec_dkd = dpec_dkd - pres * colht_chg
  endif

  if (present(dpe_max)) then
    ! This expression is the limit of PE_chg for infinite dKddt_h.
    y1 = 1.0 / (bdt1 * hps)
    dpe_max = pec_core * y1
    colht_chg = colht_core * y1
    if (colht_chg < 0.0) dpe_max = dpe_max - pres * colht_chg
  endif

  if (present(dpec_dkd_0)) then
    ! This expression is the limit of dPEc_dKd for dKddt_h = 0.
    y1 = 1.0 / bdt1**2
    dpec_dkd_0 = pec_core * y1
    colht_chg = colht_core * y1
    if (colht_chg < 0.0) dpec_dkd_0 = dpec_dkd_0 - pres * 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, 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, in units of H (m or kg-2).
  real, intent(in)  :: h_k      !< The thickness of the layer below the interface, in H.
  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, in H.
  real, intent(in)  :: dTe_term !< A diffusivity-independent term related to the
                                !! temperature change in the layer below the interface, in K H.
  real, intent(in)  :: dSe_term !< A diffusivity-independent term related to the
                                !! salinity change in the layer below the interface, in ppt H.
  real, intent(in)  :: dT_km1_t2 !< A diffusivity-independent term related to the
                                 !! temperature change in the layer above the interface, in K.
  real, intent(in)  :: dS_km1_t2 !< A diffusivity-independent term related to the
                                 !! salinity change in the layer above the interface, in ppt.
  real, intent(in)  :: pres      !< The hydrostatic interface pressure, which is used to relate
                                 !! the changes in column thickness to the energy that is radiated
                                 !! as gravity waves and unavailable to drive mixing, in Pa.
  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, in J m-2 K-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 salinities of all the layers below, in 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, in J m-2 K-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, in 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, in m K-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, in 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, in m K-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, in m ppt-1.

  real, optional, intent(out) :: PE_chg   !< The change in column potential energy from applying
                                          !! Kddt_h at the present interface, in J m-2.
  real, optional, intent(out) :: dPEc_dKd !< The partial derivative of PE_chg with Kddt_h,
                                          !! in units of J m-2 H-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, in 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, in J m-2 H-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, in H-1.
  real :: b1Kd          ! Temporary array (nondim.)
  real :: ColHt_chg     ! The change in column thickness in m.
  real :: dColHt_max    ! The change in column thickess for infinite diffusivity, in m.
  real :: dColHt_dKd    ! The partial derivative of column thickess with diffusivity, in s m-1.
  real :: dT_k, dT_km1  ! Temporary arrays in K.
  real :: dS_k, dS_km1  ! Temporary arrays in ppt.
  real :: I_Kr_denom, dKr_dKd   ! Temporary arrays in H-2 and nondim.
  real :: ddT_k_dKd, ddT_km1_dKd ! Temporary arrays in K H-1.
  real :: ddS_k_dKd, ddS_km1_dKd ! Temporary arrays in ppt H-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 * 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 * 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*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*dcolht_dkd
  endif

end subroutine find_pe_chg_orig

!> Copies the ePBL active mixed layer depth into MLD
subroutine energetic_pbl_get_mld(CS, MLD, G)
  type(energetic_pbl_cs),           pointer     :: CS  !< Control structure for ePBL
  type(ocean_grid_type),            intent(in)  :: G   !< Grid structure
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: MLD !< Depth of ePBL active mixing layer
  ! Local variables
  integer :: i,j
  do j = g%jsc, g%jec ; do i = g%isc, g%iec
    mld(i,j) = cs%ML_depth(i,j)
  enddo ; enddo
end subroutine energetic_pbl_get_mld

!> Computes wind speed from ustar_air based on COARE 3.5 Cd relationship
subroutine ust_2_u10_coare3p5(USTair,U10,GV)
  real, intent(in)  :: USTair
  type(verticalgrid_type), intent(in) :: GV   !< The ocean's vertical grid structure
  real, intent(out) :: U10
  real, parameter :: vonkar = 0.4
  real, parameter :: nu=1e-6
  real :: z0sm, z0, z0rough, u10a, alpha, CD
  integer :: CT

  ! Uses empirical formula for z0 to convert ustar_air to u10 based on the
  !  COARE 3.5 paper (Edson et al., 2013)
  !alpha=m*U10+b
  !Note in Edson et al. 2013, eq. 13 m is given as 0.017.  However,
  ! m=0.0017 reproduces the curve in their figure 6.

  z0sm = 0.11 * nu / ustair; !Compute z0smooth from ustar guess
  u10 = ustair/sqrt(0.001);  !Guess for u10
  u10a = 1000;

  ct=0
  do while (abs(u10a/u10-1.)>0.001)
    ct=ct+1
    u10a = u10
    alpha = min(0.028,0.0017 * u10 - 0.005)
    z0rough = alpha * ustair**2/gv%g_Earth ! Compute z0rough from ustar guess
    z0=z0sm+z0rough
    cd = ( vonkar / log(10/z0) )**2 ! Compute CD from derived roughness
    u10 = ustair/sqrt(cd);!Compute new u10 from derived CD, while loop
                       ! ends and checks for convergence...CT counter
                       ! makes sure loop doesn't run away if function
                       ! doesn't converge.  This code was produced offline
                       ! and converged rapidly (e.g. 2 cycles)
                       ! for ustar=0.0001:0.0001:10.
    if (ct>20) then
      u10 = ustair/sqrt(0.0015) ! I don't expect to get here, but just
                              !  in case it will output a reasonable value.
      exit
    endif
  enddo
  return
end subroutine ust_2_u10_coare3p5

subroutine get_la_windsea(ustar, hbl, GV, LA)
! Original description:
! This function returns the enhancement factor, given the 10-meter
! wind (m/s), friction velocity (m/s) and the boundary layer depth (m).
! Update (Jan/25):
! Converted from function to subroutine, now returns Langmuir number.
! Computs 10m wind internally, so only ustar and hbl need passed to
! subroutine.
!
! Qing Li, 160606
! BGR port from CVMix to MOM6 Jan/25/2017
! BGR change output to LA from Efactor
! BGR remove u10 input

! Input
  real, intent(in) :: &
       ! water-side surface friction velocity (m/s)
       ustar, &
       ! boundary layer depth (m)
       hbl
  type(verticalgrid_type), intent(in) :: GV   !< The ocean's vertical grid structure
  real, intent(out) :: LA
! Local variables
  ! parameters
  real, parameter :: &
       ! ratio of U19.5 to U10 (Holthuijsen, 2007)
       u19p5_to_u10 = 1.075, &
       ! ratio of mean frequency to peak frequency for
       ! Pierson-Moskowitz spectrum (Webb, 2011)
       fm_to_fp = 1.296, &
       ! ratio of surface Stokes drift to U10
       us_to_u10 = 0.0162, &
       ! loss ratio of Stokes transport
       r_loss = 0.667
  real :: us, hm0, fm, fp, vstokes, kphil, kstar
  real :: z0, z0i, r1, r2, r3, r4, tmp, us_sl, lasl_sqr_i
  real :: pi, u10
  pi = 4.0*atan(1.0)
  ! Computing u10 based on u_star and COARE 3.5 relationships
  call ust_2_u10_coare3p5(ustar*sqrt(gv%Rho0/1.225),u10,gv)
  if (u10 .gt. 0.0 .and. ustar .gt. 0.0) then
    ! surface Stokes drift
    us = us_to_u10*u10
    !
    ! significant wave height from Pierson-Moskowitz
    ! spectrum (Bouws, 1998)
    hm0 = 0.0246 *u10**2
    !
    ! peak frequency (PM, Bouws, 1998)
    tmp = 2.0 * pi * u19p5_to_u10 * u10
    fp = 0.877 * gv%g_Earth / tmp
    !
    ! mean frequency
    fm = fm_to_fp * fp
    !
    ! total Stokes transport (a factor r_loss is applied to account
    !  for the effect of directional spreading, multidirectional waves
    !  and the use of PM peak frequency and PM significant wave height
    !  on estimating the Stokes transport)
    vstokes = 0.125 * pi * r_loss * fm * hm0**2
    !
    ! the general peak wavenumber for Phillips' spectrum
    ! (Breivik et al., 2016) with correction of directional spreading
    kphil = 0.176 * us / vstokes
    !
    ! surface layer averaged Stokes dirft with Stokes drift profile
    ! estimated from Phillips' spectrum (Breivik et al., 2016)
    ! the directional spreading effect from Webb and Fox-Kemper, 2015
    ! is also included
    kstar = kphil * 2.56
    ! surface layer
    !z0 = 0.2 * abs(hbl)
    !BGR hbl now adjusted by averaging ratio before function call.
    z0 = abs(hbl)
    z0i = 1.0 / z0
    ! term 1 to 4
    r1 = ( 0.151 / kphil * z0i -0.84 ) &
         * ( 1.0 - exp(-2.0 * kphil * z0) )
    r2 = -( 0.84 + 0.0591 / kphil * z0i ) &
         *sqrt( 2.0 * pi * kphil * z0 ) &
         *erfc( sqrt( 2.0 * kphil * z0 ) )
    r3 = ( 0.0632 / kstar * z0i + 0.125 ) &
         * (1.0 - exp(-2.0 * kstar * z0) )
    r4 = ( 0.125 + 0.0946 / kstar * z0i ) &
         *sqrt( 2.0 * pi *kstar * z0) &
         *erfc( sqrt( 2.0 * kstar * z0 ) )
    us_sl = us * (0.715 + r1 + r2 + r3 + r4)
    !
    la = sqrt(ustar / us_sl)
  else
    la=1.e8
  endif
endsubroutine get_la_windsea

subroutine energetic_pbl_init(Time, G, GV, param_file, diag, CS)
  type(time_type), target, intent(in)    :: 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(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(diag_ctrl), target, intent(inout) :: diag
  type(energetic_pbl_cs), pointer        :: CS
! Arguments: Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      diag - A structure that is used to regulate diagnostic output.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                  for this module
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_energetic_PBL"  ! This module's name.
  real :: omega_frac_dflt
  integer :: isd, ied, jsd, jed
  logical :: use_temperature, use_omega
  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, "")

  call get_param(param_file, mdl, "MSTAR_MODE", cs%mstar_mode, &
                 "An integer switch for how to compute MSTAR. \n"//&
                 "    0 for constant MSTAR\n"//&
                 "    1 for MSTAR w/ MLD in stabilizing limit\n"//&
                 "    2 for MSTAR w/ L_E/L_O in stabilizing limit.",&
                 "units=nondim",default=0)
  call get_param(param_file, mdl, "MSTAR", cs%mstar, &
                 "The ratio of the friction velocity cubed to the TKE \n"//&
                 "input to the mixed layer.", "units=nondim", default=1.2)
  call get_param(param_file, mdl, "MIX_LEN_EXPONENT", cs%MixLenExponent, &
                 "The exponent applied to the ratio of the distance to the MLD \n"//&
                 "and the MLD depth which determines the shape of the mixing length.",&
                 "units=nondim", default=2.0)
  call get_param(param_file, mdl, "MSTAR_CAP", cs%mstar_cap, &
                 "Maximum value of mstar allowed in model if non-negative\n"//&
                 "(used if MSTAR_MODE>0).",&
                 "units=nondim", default=-1.0)
  call get_param(param_file, mdl, "MSTAR_CONV_ADJ", cs%cnv_mst_fac, &
                 "Factor used for reducing mstar during convection \n"//&
                 " due to reduction of stable density gradient.",&
                 "units=nondim", default=0.0)
  call get_param(param_file, mdl, "MSTAR_SLOPE", cs%mstar_slope, &
                 "The slope of the linear relationship between mstar \n"//&
                 "and the length scale ratio (used if MSTAR_MODE=1).",&
                 "units=nondim", default=0.85)
  call get_param(param_file, mdl, "MSTAR_XINT", cs%mstar_xint, &
                 "The value of the length scale ratio where the mstar \n"//&
                 "is linear above (used if MSTAR_MODE=1).",&
                 "units=nondim", default=-0.3)
  call get_param(param_file, mdl, "MSTAR_AT_XINT", cs%mstar_at_xint, &
                 "The value of mstar at MSTAR_XINT \n"//&
                 "(used if MSTAR_MODE=1).",&
                 "units=nondim", default=0.095)
  call get_param(param_file, mdl, "MSTAR_FLATCAP", cs%MSTAR_FLATCAP, &
                 "Set false to use asymptotic cap, defaults to true.\n"//&
                 "(used only if MSTAR_MODE=1)"&
                 ,"units=nondim",default=.true.)
  call get_param(param_file, mdl, "MSTAR2_COEF1", cs%MSTAR_COEF, &
                 "Coefficient in computing mstar when rotation and \n"//&
                 " stabilizing effects are both important (used if MSTAR_MODE=2)"&
                  ,"units=nondim",default=0.3)
  call get_param(param_file, mdl, "MSTAR2_COEF2", cs%C_EK, &
                 "Coefficient in computing mstar when only rotation limits \n"//&
                 " the total mixing. (used only if MSTAR_MODE=2)"&
                  ,"units=nondim",default=0.085)
  call get_param(param_file, mdl, "NSTAR", cs%nstar, &
                 "The portion of the buoyant potential energy imparted by \n"//&
                 "surface fluxes that is available to drive entrainment \n"//&
                 "at the base of mixed layer when that energy is positive.", &
                 units="nondim", default=0.2)
  call get_param(param_file, mdl, "MKE_TO_TKE_EFFIC", cs%MKE_to_TKE_effic, &
                 "The efficiency with which mean kinetic energy released \n"//&
                 "by mechanically forced entrainment of the mixed layer \n"//&
                 "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 \n"//&
                 "TKE available for mechanical entrainment to the natural \n"//&
                 "Ekman depth.", units="nondim", default=2.5)
!  call get_param(param_file, mdl, "HMIX_MIN", CS%Hmix_min, &
!                 "The minimum mixed layer depth if the mixed layer depth \n"//&
!                 "is determined dynamically.", units="m", default=0.0)

  call get_param(param_file, mdl, "OMEGA",cs%omega,              &
                 "The rotation rate of the earth.", units="s-1", &
                 default=7.2921e-5)
  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega,                  &
                 "If true, use the absolute rotation rate instead of the \n"//&
                 "vertical component of rotation when setting the decay \n"// &
                 "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 \n"//  &
                 "fraction of the absolute rotation rate blended with the \n"//&
                 "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, "WSTAR_USTAR_COEF", cs%wstar_ustar_coef,     &
                 "A ratio relating the efficiency with which convectively \n"//&
                 "released energy is converted to a turbulent velocity, \n"//  &
                 "relative to mechanically forced TKE. Making this larger \n"//&
                 "increases the BL diffusivity", units="nondim", default=1.0)
  call get_param(param_file, mdl, "VSTAR_SCALE_FACTOR", cs%vstar_scale_fac, &
                 "An overall nondimensional scaling factor for v*. \n"//    &
                 "Making this larger decreases the PBL diffusivity.",       &
                 units="nondim", default=1.0)
  call get_param(param_file, mdl, "EKMAN_SCALE_COEF", cs%Ekman_scale_coef,           &
                 "A nondimensional scaling factor controlling the inhibition \n"//   &
                 "of the diffusive length scale by rotation. Making this larger \n"//&
                 "decreases the PBL diffusivity.", units="nondim", default=1.0)
  call get_param(param_file, mdl, "USE_MLD_ITERATION", cs%USE_MLD_ITERATION,    &
                 "A logical that specifies whether or not to use the \n"//      &
                 "distance to the bottom of the actively turblent boundary \n"//&
                 "layer to help set the EPBL length scale.", default=.false.)
  call get_param(param_file, mdl, "ORIG_MLD_ITERATION", cs%ORIG_MLD_ITERATION,  &
                 "A logical that specifies whether or not to use the \n"//      &
                 "old method for determining MLD depth in iteration, which \n"//&
                 "is limited to resolution.", default=.true.)
  call get_param(param_file, mdl, "MLD_ITERATION_GUESS", cs%MLD_ITERATION_GUESS,       &
                 "A logical that specifies whether or not to use the \n"//             &
                 "previous timestep MLD as a first guess in the MLD iteration.\n"//    &
                 "The default is false to facilitate reproducibility.", default=.false.)
  call get_param(param_file, mdl, "EPBL_MLD_TOLERANCE", cs%MLD_tol,         &
                 "The tolerance for the iteratively determined mixed \n"//  &
                 "layer depth.  This is only used with USE_MLD_ITERATION.", &
                 units="meter", default=1.0)
  call get_param(param_file, mdl, "EPBL_MIN_MIX_LEN", cs%min_mix_len,    &
                 "The minimum mixing length scale that will be used \n"//&
                 "by ePBL.  The default (0) does not set a minimum.",    &
                 units="meter", default=0.0)
  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 \n"//        &
                 "potential energy change code.  Otherwise, the newer \n"//         &
                 "version that can work with successive increments to the \n"//     &
                 "diffusivity in upward or downward passes is used.", default=.true.)
  call get_param(param_file, mdl, "EPBL_TRANSITION_SCALE", cs%transLay_scale, &
                 "A scale for the mixing length in the transition layer \n"// &
                 "at the edge of the boundary layer as a fraction of the \n"//&
                 "boundary layer thickness.  The default is 0.1.", &
                 units="nondim", default=0.1)
  if ( cs%USE_MLD_ITERATION .and. abs(cs%transLay_scale-0.5).ge.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, "N2_DISSIPATION_POS", cs%N2_Dissipation_Scale_Pos, &
                 "A scale for the dissipation of TKE due to stratification \n"//     &
                 "in the boundary layer, applied when local stratification \n"//     &
                 "is positive.  The default is 0, but should probably be ~0.4.",     &
                 units="nondim", default=0.0)
  call get_param(param_file, mdl, "N2_DISSIPATION_NEG", cs%N2_Dissipation_Scale_Neg,&
                 "A scale for the dissipation of TKE due to stratification \n"//    &
                 "in the boundary layer, applied when local stratification \n"//    &
                 "is negative.  The default is 0, but should probably be ~1.",      &
                 units="nondim", default=0.0)
   call get_param(param_file, mdl, "USE_LA_LI2016", cs%USE_LA_Windsea,            &
                 "A logical to use the Li et al. 2016 (submitted) formula to \n"//&
                 " determine the Langmuir number.",                               &
                 units="nondim", default=.false.)
   call get_param(param_file, mdl, "LA_DEPTH_RATIO", cs%LaDepthRatio,                &
                 "The depth (normalized by BLD) to average Stokes drift over in \n"//&
                 " Lanmguir number calculation, where La = sqrt(ust/Stokes).",       &
                 units="nondim",default=0.04)
   call get_param(param_file, mdl, "LT_ENHANCE", cs%LT_ENHANCE_FORM,        &
                 "Integer for Langmuir number mode. \n"//                   &
                 " *Requires USE_LA_LI2016 to be set to True. \n"//         &
                 "Options: 0 - No Langmuir \n"//                            &
                 "         1 - Van Roekel et al. 2014/Li et al., 2016  \n"//&
                 "         2 - Multiplied w/ adjusted La. \n"//             &
                 "         3 - Added w/ adjusted La.",                      &
                 units="nondim", default=0)
   call get_param(param_file, mdl, "LT_ENHANCE_COEF", cs%LT_ENHANCE_COEF, &
                 "Coefficient for Langmuir enhancement if LT_ENHANCE > 1",&
                 units="nondim", default=0.447)
   call get_param(param_file, mdl, "LT_ENHANCE_EXP", cs%LT_ENHANCE_EXP, &
                 "Exponent for Langmuir enhancement if LT_ENHANCE > 1", &
                 units="nondim", default=-1.33)
   call get_param(param_file, mdl, "LT_MOD_LAC1", cs%LaC_MLDoEK,             &
                 "Coefficient for modification of Langmuir number due to\n"//&
                 " MLD approaching Ekman depth if LT_ENHANCE=2.",            &
                 units="nondim", default=-0.87)
   call get_param(param_file, mdl, "LT_MOD_LAC2", cs%LaC_MLDoOB_stab,        &
                 "Coefficient for modification of Langmuir number due to\n"//&
                 " MLD approaching stable Obukhov depth if LT_ENHANCE=2.",   &
                  units="nondim", default=0.0)
   call get_param(param_file, mdl, "LT_MOD_LAC3", cs%LaC_MLDoOB_un,          &
                 "Coefficient for modification of Langmuir number due to\n"//&
                 " MLD approaching unstable Obukhov depth if LT_ENHANCE=2.", &
                  units="nondim", default=0.0)
   call get_param(param_file, mdl, "LT_MOD_LAC4", cs%Lac_EKoOB_stab,         &
                 "Coefficient for modification of Langmuir number due to\n"//&
                 " ratio of Ekman to stable Obukhov depth if LT_ENHANCE=2.", &
                  units="nondim", default=0.95)
   call get_param(param_file, mdl, "LT_MOD_LAC5", cs%Lac_EKoOB_un,            &
                 "Coefficient for modification of Langmuir number due to\n"// &
                 " ratio of Ekman to unstable Obukhov depth if LT_ENHANCE=2.",&
                  units="nondim", default=0.95)
   if (cs%LT_ENHANCE_FORM.gt.0 .and. (.not.cs%USE_LA_Windsea)) then
      call mom_error(fatal, "If flag USE_LA_LI2016 is false, then "//&
                 " LT_ENHANCE must be 0.")
  endif
  ! 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_m*gv%H_subroundoff)
  call log_param(param_file, mdl, "EPBL_USTAR_MIN", cs%ustar_min, &
                 "The (tiny) minimum friction velocity used within the \n"//&
                 "ePBL code, derived from OMEGA and ANGSTROM.", units="meter second-1")

  cs%id_ML_depth = register_diag_field('ocean_model', 'ePBL_h_ML', diag%axesT1, &
      time, 'Surface mixed layer depth', 'meter')
  cs%id_TKE_wind = register_diag_field('ocean_model', 'ePBL_TKE_wind', diag%axesT1, &
      time, 'Wind-stirring source of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_MKE = register_diag_field('ocean_model', 'ePBL_TKE_MKE', diag%axesT1, &
      time, 'Mean kinetic energy source of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_conv = register_diag_field('ocean_model', 'ePBL_TKE_conv', diag%axesT1, &
      time, 'Convective source of mixed layer TKE', 'meter3 second-3')
  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', 'meter3 second-3')
  cs%id_TKE_mixing = register_diag_field('ocean_model', 'ePBL_TKE_mixing', diag%axesT1, &
      time, 'TKE consumed by mixing that deepens the mixed layer', 'meter3 second-3')
  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', 'meter3 second-3')
  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', 'meter3 second-3')
  cs%id_Hsfc_used = register_diag_field('ocean_model', 'ePBL_Hs_used', diag%axesT1, &
      time, 'Surface region thickness that is used', 'meter')
  cs%id_Mixing_Length = register_diag_field('ocean_model', 'Mixing_Length', diag%axesTi, &
      time, 'Mixing Length that is used', 'meter')
  cs%id_Velocity_Scale = register_diag_field('ocean_model', 'Velocity_Scale', diag%axesTi, &
      time, 'Velocity Scale that is used.', 'meter second-1')
  cs%id_LT_enhancement = register_diag_field('ocean_model', 'LT_Enhancement', diag%axesT1, &
      time, 'LT enhancement that is used.', 'non-dim')
  cs%id_MSTAR_mix = register_diag_field('ocean_model', 'MSTAR', diag%axesT1, &
      time, 'MSTAR that is used.', 'non-dim')
  cs%id_OSBL = register_diag_field('ocean_model', 'ePBL_OSBL', diag%axesT1, &
      time, 'Boundary layer depth from the iteration.', 'meter')
  cs%id_mld_ekman = register_diag_field('ocean_model', 'MLD_EKMAN', diag%axesT1, &
      time, 'Boundary layer depth over Ekman length.', 'meter')
  cs%id_mld_obukhov = register_diag_field('ocean_model', 'MLD_OBUKHOV', diag%axesT1, &
      time, 'Boundary layer depth over Obukhov length.', 'meter')
  cs%id_ekman_obukhov = register_diag_field('ocean_model', 'EKMAN_OBUKHOV', diag%axesT1, &
      time, 'Ekman length over Obukhov length.', 'meter')
  cs%id_LA = register_diag_field('ocean_model', 'LA', diag%axesT1, &
      time, 'Langmuir number.', 'non-dim')
  cs%id_LA_mod = register_diag_field('ocean_model', 'LA_MOD', diag%axesT1, &
      time, 'Modified Langmuir number.', 'non-dim')


  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
                 "If true, temperature and salinity are used as state \n"//&
                 "variables.", default=.true.)

  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_Mixing_Length>0) .or. (cs%id_Velocity_Scale>0)) then
    call safe_alloc_alloc(cs%Velocity_Scale,isd,ied,jsd,jed,gv%ke+1)
    call safe_alloc_alloc(cs%Mixing_Length,isd,ied,jsd,jed,gv%ke+1)
    cs%Velocity_Scale(:,:,:) = 0.0
    cs%Mixing_Length(:,:,:) = 0.0
    cs%mixing_diagnostics = .true.
  endif
  call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)
  call safe_alloc_alloc(cs%ML_depth2, isd, ied, jsd, jed)
  if (max(cs%id_LT_Enhancement, cs%id_mstar_mix,cs%id_mld_ekman, &
       cs%id_ekman_obukhov, cs%id_mld_obukhov, cs%id_LA, cs%id_LA_mod)>0) then
    call safe_alloc_alloc(cs%Mstar_mix, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%Enhance_M, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%MLD_EKMAN, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%MLD_OBUKHOV, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%EKMAN_OBUKHOV, 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)
  endif

  !Fitting coefficients to asymptote twoard 0 as MLD -> Ekman depth
  cs%MSTAR_A = cs%MSTAR_AT_XINT**(1./cs%MSTAR_N)
  cs%MSTAR_B = cs%MSTAR_SLOPE / (cs%MSTAR_N*cs%MSTAR_A**(cs%MSTAR_N-1.))
  !Fitting coefficients to asymptote toward MSTAR_CAP
  !*Fixed to begin asymptote at MSTAR_CAP-0.5 toward MSTAR_CAP
  cs%MSTAR_A2 = 0.5**(1./cs%MSTAR_N)
  cs%MSTAR_B2 = -cs%MSTAR_SLOPE / (cs%MSTAR_N*cs%MSTAR_A2**(cs%MSTAR_N-1))
  !Compute value of X (referenced to MSTAR_XINT) where transition
  ! to asymptotic regime based on value of X where MSTAR=MSTAR_CAP-0.5
  cs%MSTAR_XINT_UP = (cs%MSTAR_CAP-0.5-cs%MSTAR_AT_XINT)/cs%MSTAR_SLOPE

end subroutine energetic_pbl_init

subroutine energetic_pbl_end(CS)
  type(energetic_pbl_cs), pointer :: CS

  if (.not.associated(cs)) return

  if (allocated(cs%ML_depth))            deallocate(cs%ML_depth)
  if (allocated(cs%ML_depth2))           deallocate(cs%ML_depth2)
  if (allocated(cs%Enhance_M))           deallocate(cs%Enhance_M)
  if (allocated(cs%MLD_EKMAN))           deallocate(cs%MLD_EKMAN)
  if (allocated(cs%MLD_OBUKHOV))         deallocate(cs%MLD_OBUKHOV)
  if (allocated(cs%EKMAN_OBUKHOV))       deallocate(cs%EKMAN_OBUKHOV)
  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%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)

  deallocate(cs)

end subroutine energetic_pbl_end

end module mom_energetic_pbl