MOM_bulk_mixed_layer.F90

Go to the documentation of this file.

module mom_bulk_mixed_layer
!***********************************************************************
!*                   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, 1997 - 2005.                                   *
!*                                                                     *
!*    This file contains the subroutine (bulkmixedlayer) that          *
!*  implements a Kraus-Turner-like bulk mixed layer, based on the work *
!*  of various people, as described in the review paper by Niiler and  *
!*  Kraus (1979), with particular attention to the form proposed by    *
!*  Oberhuber (JPO, 1993, 808-829), with an extension to a refied bulk *
!*  mixed layer as described in Hallberg (Aha Huliko'a, 2003).  The    *
!*  physical processes portrayed in this subroutine include convective *
!*  adjustment and mixed layer entrainment and detrainment.            *
!*  Penetrating shortwave radiation and an exponential decay of TKE    *
!*  fluxes are also supported by this subroutine.  Several constants   *
!*  can alternately be set to give a traditional Kraus-Turner mixed    *
!*  layer scheme, although that is not the preferred option.  The      *
!*  physical processes and arguments are described in detail below.    *
!*                                                                     *
!*  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, buoy, T, S, eaml, ebml, 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, diag_update_remap_grids
use mom_domains,       only : create_group_pass, do_group_pass, group_pass_type
use mom_error_handler, only : mom_error, fatal, warning
use mom_file_parser,   only : get_param, log_param, log_version, param_file_type
use mom_forcing_type,  only : extractfluxes1d, forcing
use mom_grid,          only : ocean_grid_type
use mom_shortwave_abs, only : absorbremainingsw, optics_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 bulkmixedlayer, bulkmixedlayer_init

type, public :: bulkmixedlayer_cs ; private
  integer :: nkml            ! The number of layers in the mixed layer.
  integer :: nkbl            ! The number of buffer layers.
  integer :: nsw             ! The number of bands of penetrating shortwave radiation.
  real    :: mstar           ! The ratio of the friction velocity cubed to the
                             ! TKE input to the mixed layer, nondimensional.
  real    :: nstar           ! The fraction of the TKE input to the mixed layer
                             ! available to drive entrainment, nondim.
  real    :: nstar2          ! The fraction of potential energy released by
                             ! convective adjustment that drives entrainment, ND.
  logical :: absorb_all_sw   ! If true, all shortwave radiation is absorbed by the
                             ! ocean, instead of passing through to the bottom mud.
  real    :: tke_decay       ! The ratio of the natural Ekman depth to the TKE
                             ! decay scale, nondimensional.
  real    :: bulk_ri_ml      ! The efficiency with which mean kinetic energy
                             ! released by mechanically forced entrainment of
                             ! the mixed layer is converted to TKE, nondim.
  real    :: bulk_ri_convective ! The efficiency with which convectively
                             ! released mean kinetic energy becomes TKE, nondim.
  real    :: hmix_min        ! The minimum mixed layer thickness in m.
  real    :: h_limit_fluxes  ! When the total ocean depth is less than this
                             ! value, in m, scale away all surface forcing to
                             ! avoid boiling the ocean.
  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    :: dt_ds_wt        !   When forced to extrapolate T & S to match the
                             ! layer densities, this factor (in deg C / PSU) is
                             ! combined with the derivatives of density with T & S
                             ! to determines what direction is orthogonal to
                             ! density contours.  It should be a typical value of
                             ! (dR/dS) / (dR/dT) in oceanic profiles.
                             ! 6 K psu-1 might be reasonable.
  real    :: bl_extrap_lim   ! A limit on the density range over which
                             ! extrapolation can occur when detraining from the
                             ! buffer layers, relative to the density range
                             ! within the mixed and buffer layers, when the
                             ! detrainment is going into the lightest interior
                             ! layer, nondimensional.
  logical :: ml_resort       !   If true, resort the layers by density, rather than
                             ! doing convective adjustment.
  integer :: ml_presort_nz_conv_adj ! If ML_resort is true, do convective
                             ! adjustment on this many layers (starting from the
                             ! top) before sorting the remaining layers.
  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).
  logical :: correct_absorption ! If true, the depth at which penetrating
                             ! shortwave radiation is absorbed is corrected by
                             ! moving some of the heating upward in the water
                             ! column.  The default is false.
  logical :: resolve_ekman   !   If true, the nkml layers in the mixed layer are
                             ! chosen to optimally represent the impact of the
                             ! Ekman transport on the mixed layer TKE budget.
  type(time_type), pointer :: time ! A pointer to the ocean model's clock.
  logical :: tke_diagnostics = .false.
  logical :: do_rivermix = .false. ! Provide additional TKE to mix river runoff
                                   ! at the river mouths to "rivermix_depth" meters
  real    :: rivermix_depth = 0.0  ! Used if "do_rivermix" = T
  logical :: limit_det       ! If true, limit the extent of buffer layer
                             ! detrainment to be consistent with neighbors.
  real    :: lim_det_dh_sfc  ! The fractional limit in the change between grid
                             ! points of the surface region (mixed & buffer
                             ! layer) thickness, nondim.  0.5 by default.
  real    :: lim_det_dh_bathy ! The fraction of the total depth by which the
                             ! thickness of the surface region (mixed & buffer
                             ! layer) is allowed to change between grid points.
                             ! Nondimensional, 0.2 by default.
  logical :: use_river_heat_content ! If true, use the fluxes%runoff_Hflx field
                             ! to set the heat carried by runoff, instead of
                             ! using SST for temperature of liq_runoff
  logical :: use_calving_heat_content ! Use SST for temperature of froz_runoff
  logical :: salt_reject_below_ml ! It true, add salt below mixed layer (layer mode only)

  type(diag_ctrl), pointer :: diag ! A structure that is used to regulate the
                             ! timing of diagnostic output.
  real    :: allowed_t_chg   ! The amount by which temperature is allowed
                             ! to exceed previous values during detrainment, K.
  real    :: allowed_s_chg   ! The amount by which salinity is allowed
                             ! to exceed previous values during detrainment, PSU.

! These are terms in the mixed layer TKE budget, all in m3 s-2.
  real, allocatable, dimension(:,:) :: &
    ml_depth, &        ! The mixed layer depth in m.
    diag_tke_wind, &   ! The wind source of TKE.
    diag_tke_ribulk, & ! The resolved KE source of TKE.
    diag_tke_conv, &   ! The convective source of TKE.
    diag_tke_pen_sw, & ! The TKE sink required to mix
                       ! 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.
    diag_tke_conv_s2, &! The convective source of TKE due to
                       ! to mixing in sigma2.
    diag_pe_detrain, & ! The spurious source of potential
                       ! energy due to mixed layer
                       ! detrainment, W m-2.
    diag_pe_detrain2   ! The spurious source of potential
                       ! energy due to mixed layer only
                       ! detrainment, W m-2.
  logical :: allow_clocks_in_omp_loops  ! If true, clocks can be called
                                        ! from inside loops that can be threaded.
                                        ! To run with multiple threads, set to False.
  type(group_pass_type) :: pass_h_sum_hmbl_prev ! For group halo pass
  integer :: id_ml_depth = -1, id_tke_wind = -1, id_tke_mixing = -1
  integer :: id_tke_ribulk = -1, id_tke_conv = -1, id_tke_pen_sw = -1
  integer :: id_tke_mech_decay = -1, id_tke_conv_decay = -1, id_tke_conv_s2 = -1
  integer :: id_pe_detrain = -1, id_pe_detrain2 = -1, id_h_mismatch = -1
  integer :: id_hsfc_used = -1, id_hsfc_max = -1, id_hsfc_min = -1
end type bulkmixedlayer_cs

integer :: id_clock_detrain=0, id_clock_mech=0, id_clock_conv=0, id_clock_adjustment=0
integer :: id_clock_eos=0, id_clock_resort=0, id_clock_pass=0

integer :: num_msg = 0, max_msg = 2

contains

!>    This subroutine partially steps the bulk mixed layer model.
!!  The following processes are executed, in the order listed.
!!    1. Undergo convective adjustment into mixed layer.
!!    2. Apply surface heating and cooling.
!!    3. Starting from the top, entrain whatever fluid the TKE budget
!!       permits.  Penetrating shortwave radiation is also applied at
!!      this point.
!!    4. If there is any unentrained fluid that was formerly in the
!!      mixed layer, detrain this fluid into the buffer layer.  This
!!     is equivalent to the mixed layer detraining to the Monin-
!!       Obukhov depth.
!!    5. Divide the fluid in the mixed layer evenly into CS%nkml pieces.
!!    6. Split the buffer layer if appropriate.
!! Layers 1 to nkml are the mixed layer, nkml+1 to nkml+nkbl are the
!! buffer layers. The results of this subroutine are mathematically
!! identical if there are multiple pieces of the mixed layer with
!! the same density or if there is just a single layer. There is no
!! stability limit on the time step.
!!
!!   The key parameters for the mixed layer are found in the control structure.
!! These include mstar, nstar, nstar2, pen_SW_frac, pen_SW_scale, and TKE_decay.
!!   For the Oberhuber (1993) mixed layer, the values of these are:
!!      pen_SW_frac = 0.42, pen_SW_scale = 15.0 m, 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.
!!  Conv_decay has been eliminated in favor of the well-calibrated form for the
!!  efficiency of penetrating convection from Wang (2003).
!!    For a traditional Kraus-Turner mixed layer, the values are:
!!      pen_SW_frac = 0.0, pen_SW_scale = 0.0 m, mstar = 1.25,
!!      nstar = 0.4, TKE_decay = 0.0, conv_decay = 0.0
subroutine bulkmixedlayer(h_3d, u_3d, v_3d, tv, fluxes, dt, ea, eb, G, GV, CS, &
                          optics, Hml, aggregate_FW_forcing, dt_diag, last_call)
  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.
  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)), &
                              intent(inout) :: ea     !< The amount of fluid moved downward into a
                                                      !! layer; this should be increased due to
                                                      !! mixed layer detrainment, in the same units
                                                      !! as h - usually m or kg m-2 (i.e., H).
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                              intent(inout) :: eb     !< The amount of fluid moved upward into a
                                                      !! layer; this should be increased due to
                                                      !! mixed layer entrainment, in the same units
                                                      !! as h - usually m or kg m-2 (i.e., H).
  type(bulkmixedlayer_cs),    pointer       :: CS     !< The control structure returned by a
                                                      !! previous call to mixedlayer_init.
  type(optics_type),          pointer       :: optics !< The structure containing the inverse of the
                                                      !! vertical absorption decay scale for
                                                      !! penetrating shortwave radiation, in m-1.
  real, dimension(:,:),       pointer       :: Hml    !< active mixed layer depth
  logical,                    intent(in)    :: aggregate_FW_forcing
  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.

!    This subroutine partially steps the bulk mixed layer model.
!  The following processes are executed, in the order listed.
!    1. Undergo convective adjustment into mixed layer.
!    2. Apply surface heating and cooling.
!    3. Starting from the top, entrain whatever fluid the TKE budget
!       permits.  Penetrating shortwave radiation is also applied at
!       this point.
!    4. If there is any unentrained fluid that was formerly in the
!       mixed layer, detrain this fluid into the buffer layer.  This
!       is equivalent to the mixed layer detraining to the Monin-
!       Obukhov depth.
!    5. Divide the fluid in the mixed layer evenly into CS%nkml pieces.
!    6. Split the buffer layer if appropriate.
! Layers 1 to nkml are the mixed layer, nkml+1 to nkml+nkbl are the
! buffer layers. The results of this subroutine are mathematically
! identical if there are multiple pieces of the mixed layer with
! the same density or if there is just a single layer. There is no
! stability limit on the time step.
!
!   The key parameters for the mixed layer are found in the control structure.
! These include mstar, nstar, nstar2, pen_SW_frac, pen_SW_scale, and TKE_decay.
!   For the Oberhuber (1993) mixed layer, the values of these are:
!      pen_SW_frac = 0.42, pen_SW_scale = 15.0 m, 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.
!  Conv_decay has been eliminated in favor of the well-calibrated form for the
!  efficiency of penetrating convection from Wang (2003).
!    For a traditional Kraus-Turner mixed layer, the values are:
!      pen_SW_frac = 0.0, pen_SW_scale = 0.0 m, 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.
!  (in/out)  ea - The amount of fluid moved downward into a layer; this should
!                 be increased due to mixed layer detrainment, in the same units
!                 as h - usually m or kg m-2 (i.e., H).
!  (in/out)  eb - The amount of fluid moved upward into a layer; this should
!                 be increased due to mixed layer entrainment, in the same units
!                 as h - usually m or kg m-2 (i.e., H).
!  (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)      optics - The structure containing the inverse of the vertical
!                     absorption decay scale for penetrating shortwave
!                     radiation, in m-1.
!  (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)) :: &
    eaml, &     !   The amount of fluid moved downward into a layer due to mixed
                ! mixed layer detrainment, in m. (I.e. entrainment from above.)
    ebml        !   The amount of fluid moved upward into a layer due to mixed
                ! mixed layer detrainment, in m. (I.e. entrainment from below.)

  ! If there is resorting, the vertical coordinate for these variables is the
  ! new, sorted index space.  Here layer 0 is an initially massless layer that
  ! will be used to hold the new mixed layer properties.
  real, dimension(SZI_(G),SZK0_(GV)) :: &
    h, &        !   The layer thickness, in m or kg m-2.
    T, &        !   The layer temperatures, in deg C.
    S, &        !   The layer salinities, in psu.
    R0, &       !   The potential density referenced to the surface, in kg m-3.
    Rcv         !   The coordinate variable potential density, in kg m-3.
  real, dimension(SZI_(G),SZK_(GV)) :: &
    u, &        !   The zonal velocity, in m s-1.
    v, &        !   The meridional velocity, in m s-1.
    h_orig, &   !   The original thickness in m or kg m-2.
    d_eb, &     !   The downward increase across a layer in the entrainment from
                ! below, in H.  The sign convention is that positive values of
                ! d_eb correspond to a gain in mass by a layer by upward motion.
    d_ea, &     !   The upward increase across a layer in the entrainment from
                ! above, in H.  The sign convention is that positive values of
                ! d_ea mean a net gain in mass by a layer from downward motion.
    eps         ! The (small) thickness that must remain in a layer, in H.
  integer, dimension(SZI_(G),SZK_(GV)) :: &
    ksort       !   The sorted k-index that each original layer goes to.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    h_miss      !   The summed absolute mismatch, in m.
  real, dimension(SZI_(G)) :: &
    TKE, &      !   The turbulent kinetic energy available for mixing over a
                ! time step, in m3 s-2.
    conv_en, &  !   The turbulent kinetic energy source due to mixing down to
                ! the depth of free convection, in m3 s-2.
    htot, &     !   The total depth of the layers being considered for
                ! entrainment, in H.
    r0_tot, &   !   The integrated potential density referenced to the surface
                ! of the layers which are fully entrained, in H kg m-3.
    rcv_tot, &  !   The integrated coordinate value potential density of the
                ! layers that are fully entrained, in H kg m-3.
    ttot, &     !   The integrated temperature of layers which are fully
                ! entrained, in H K.
    stot, &     !   The integrated salt of layers which are fully entrained,
                ! in H PSU.
    uhtot, &    !   The depth integrated zonal and meridional velocities in the
    vhtot, &    ! mixed layer, in H m s-1.

    netmassinout, &  ! The net mass flux (if non-Boussinsq) or volume flux (if
                     ! Boussinesq - i.e. the fresh water flux (P+R-E)) into the
                     ! ocean over a time step, in H.
    netmassout,   &  ! The mass flux (if non-Boussinsq) or volume flux (if
                     ! Boussinesq) over a time step from evaporating fresh water (H)
    net_heat, & !   The net heating at the surface over a time step in K H.  Any
                ! penetrating shortwave radiation is not included in Net_heat.
    net_salt, & ! The surface salt flux into the ocean over a time step, psu H.
    idecay_len_tke, &  ! The inverse of a turbulence decay length scale, in H-1.
    p_ref, &    !   Reference pressure for the potential density governing mixed
                ! layer dynamics, almost always 0 (or 1e5) Pa.
    p_ref_cv, & !   Reference pressure for the potential density which defines
                ! the coordinate variable, set to P_Ref, in Pa.
    dr0_dt, &   !   Partial derivative of the mixed layer potential density with
                ! temperature, in units of kg m-3 K-1.
    drcv_dt, &  !   Partial derivative of the coordinate variable potential
                ! density in the mixed layer with temperature, in kg m-3 K-1.
    dr0_ds, &   !   Partial derivative of the mixed layer potential density with
                ! salinity, in units of kg m-3 psu-1.
    drcv_ds, &  !   Partial derivative of the coordinate variable potential
                ! density in the mixed layer with salinity, in kg m-3 psu-1.
    tke_river   ! The turbulent kinetic energy available for mixing at rivermouths over a
                ! time step, in m3 s-2.

  real, dimension(max(CS%nsw,1),SZI_(G)) :: &
    Pen_SW_bnd  !   The penetrating fraction of the shortwave heating integrated
                ! over a time step in each band, in K H.
  real, dimension(max(CS%nsw,1),SZI_(G),SZK_(GV)) :: &
    opacity_band ! The opacity in each band, in H-1. The indicies are band, i, k.

  real :: cMKE(2,szi_(g)) ! Coefficients of HpE and HpE^2 in calculating the
                          ! denominator of MKE_rate, in m-1 and m-2.
  real :: Irho0         ! 1.0 / rho_0
  real :: Inkml, Inkmlm1!  1.0 / REAL(nkml) and  1.0 / REAL(nkml-1)
  real :: Ih            !   The inverse of a thickness, in H-1.
  real :: Idt           !   The inverse of the timestep in s-1.
  real :: Idt_diag      !   The inverse of the timestep used for diagnostics in s-1.
  real :: RmixConst

  real, dimension(SZI_(G)) :: &
    dKE_FC, &   !   The change in mean kinetic energy due to free convection,
                ! in m3 s-2.
    h_ca        !   The depth to which convective adjustment has gone in H.
  real, dimension(SZI_(G),SZK_(GV)) :: &
    dKE_CA, &   !   The change in mean kinetic energy due to convective
                ! adjustment, in m3 s-2.
    ctke        !   The turbulent kinetic energy source due to convective
                ! adjustment, m3 s-2.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    Hsfc_max, & ! The thickness of the surface region (mixed and buffer layers)
                ! after entrainment but before any buffer layer detrainment, m.
    hsfc_used, & ! The thickness of the surface region after buffer layer
                ! detrainment, in units of m.
    hsfc_min, & ! The minimum thickness of the surface region based on the
                ! new mixed layer depth and the previous thickness of the
                ! neighboring water columns, in m.
    h_sum, &    ! The total thickness of the water column, in H.
    hmbl_prev   ! The previous thickness of the mixed and buffer layers, in H.
  real, dimension(SZI_(G)) :: &
    Hsfc, &     !   The thickness of the surface region (mixed and buffer
                ! layers before detrainment in to the interior, in H.
    max_bl_det  !   If non-negative, the maximum amount of entrainment from
                ! the buffer layers that will be allowed this time step, in H.
  real :: dHsfc, dHD ! Local copies of nondimensional parameters.
  real :: H_nbr ! A minimum thickness based on neighboring thicknesses, in H.

  real :: absf_x_H  ! The absolute value of f times the mixed layer thickness,
                    ! in units of m s-1.
  real :: kU_star   ! Ustar times the Von Karmen constant, in m s-1.
  real :: dt__diag ! A copy of dt_diag (if present) or dt, in s.
  logical :: write_diags  ! If true, write out diagnostics with this step.
  logical :: reset_diags  ! If true, zero out the accumulated diagnostics.
  integer :: i, j, k, is, ie, js, je, nz, nkmb, n
  integer :: nsw    ! The number of bands of penetrating shortwave radiation.

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

  if (.not. associated(cs)) call mom_error(fatal, "MOM_mixed_layer: "//&
         "Module must be initialized before it is used.")
  if (gv%nkml < 1) return

  if (.not. ASSOCIATED(tv%eqn_of_state)) call mom_error(fatal, &
      "MOM_mixed_layer: Temperature, salinity and an equation of state "//&
      "must now be used.")
  if (.NOT. ASSOCIATED(fluxes%ustar)) call mom_error(fatal, &
      "MOM_mixed_layer: No surface TKE fluxes (ustar) defined in mixedlayer!")

  nkmb = cs%nkml+cs%nkbl
  inkml = 1.0 / REAL(cs%nkml)
  if (cs%nkml > 1) inkmlm1 = 1.0 / REAL(cs%nkml-1)

  irho0 = 1.0 / gv%Rho0
  dt__diag = dt ; if (present(dt_diag)) dt__diag = dt_diag
  idt = 1.0/dt
  idt_diag = 1.0 / dt__diag
  write_diags = .true. ; if (present(last_call)) write_diags = last_call

  p_ref(:) = 0.0 ; p_ref_cv(:) = tv%P_Ref


  nsw = cs%nsw

  if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
!$OMP parallel default(none) shared(is,ie,js,je,nkmb,h_sum,hmbl_prev,h_3d,nz)
!$OMP do
    do j=js-1,je+1 ; do i=is-1,ie+1
      h_sum(i,j) = 0.0 ; hmbl_prev(i,j) = 0.0
    enddo ; enddo
!$OMP do
    do j=js-1,je+1
      do k=1,nkmb ; do i=is-1,ie+1
        h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
        hmbl_prev(i,j) = hmbl_prev(i,j) + h_3d(i,j,k)
      enddo ; enddo
      do k=nkmb+1,nz ; do i=is-1,ie+1
        h_sum(i,j) = h_sum(i,j) + h_3d(i,j,k)
      enddo ; enddo
    enddo
!$OMP end parallel

    call cpu_clock_begin(id_clock_pass)
    call create_group_pass(cs%pass_h_sum_hmbl_prev, h_sum,g%Domain)
    call create_group_pass(cs%pass_h_sum_hmbl_prev, hmbl_prev,g%Domain)
    call do_group_pass(cs%pass_h_sum_hmbl_prev, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  ! 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_RiBulk(i,j) = 0.0
        cs%diag_TKE_conv(i,j) = 0.0 ; cs%diag_TKE_pen_SW(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_conv_s2(i,j) = 0.0
      enddo ; enddo
    endif
    if (ALLOCATED(cs%diag_PE_detrain)) then
!$OMP do
      do j=js,je ; do i=is,ie
        cs%diag_PE_detrain(i,j) = 0.0
      enddo ; enddo
    endif
    if (ALLOCATED(cs%diag_PE_detrain2)) then
!$OMP do
      do j=js,je ; do i=is,ie
        cs%diag_PE_detrain2(i,j) = 0.0
      enddo ; enddo
    endif
!$OMP end parallel
  endif

  if (cs%ML_resort) then
    do i=is,ie ; h_ca(i) = 0.0 ; enddo
    do k=1,nz ; do i=is,ie ; dke_ca(i,k) = 0.0 ; ctke(i,k) = 0.0 ; enddo ; enddo
  endif
  max_bl_det(:) = -1

!$OMP parallel default(none) shared(is,ie,js,je,nz,h_3d,u_3d,v_3d,nkmb,G,GV,nsw,optics,      &
!$OMP                               CS,tv,fluxes,Irho0,dt,Idt_diag,Ih,write_diags,           &
!$OMP                               hmbl_prev,h_sum,Hsfc_min,Hsfc_max,dt__diag,              &
!$OMP                               Hsfc_used,Inkmlm1,Inkml,ea,eb,h_miss,Hml,                &
!$OMP                               id_clock_EOS,id_clock_resort,id_clock_adjustment,        &
!$OMP                               id_clock_conv,id_clock_mech,id_clock_detrain,aggregate_FW_forcing ) &
!$OMP                  firstprivate(dKE_CA,cTKE,h_CA,max_BL_det,p_ref,p_ref_cv)              &
!$OMP                       private(h,h_orig,u,v,eps,T,S,opacity_band,d_ea,d_eb,             &
!$OMP                               dR0_dT,dR0_dS,dRcv_dT,dRcv_dS,R0,Rcv,ksort,              &
!$OMP                               RmixConst,TKE_river,netMassInOut, NetMassOut,            &
!$OMP                               Net_heat, Net_salt, htot,TKE,Pen_SW_bnd,Ttot,Stot, uhtot,&
!$OMP                               vhtot, R0_tot, Rcv_tot,Conv_en,dKE_FC,Idecay_len_TKE,    &
!$OMP                               cMKE,Hsfc,dHsfc,dHD,H_nbr,kU_Star,absf_x_H,              &
!$OMP                               ebml,eaml)
!$OMP do
  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) ; u(i,k) = u_3d(i,j,k) ; v(i,k) = v_3d(i,j,k)
      h_orig(i,k) = h_3d(i,j,k)
      eps(i,k) = 0.0 ; if (k > nkmb) eps(i,k) = gv%Angstrom
      t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
      do n=1,nsw
        opacity_band(n,i,k) = gv%H_to_m*optics%opacity_band(n,i,j,k)
      enddo
    enddo ; enddo

    do k=1,nz ; do i=is,ie
      d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0
    enddo ; enddo

    if(id_clock_eos>0) call cpu_clock_begin(id_clock_eos)
    ! Calculate an estimate of the mid-mixed layer pressure (in Pa)
    do i=is,ie ; p_ref(i) = 0.0 ; enddo
    do k=1,cs%nkml ; do i=is,ie
      p_ref(i) = p_ref(i) + 0.5*gv%H_to_Pa*h(i,k)
    enddo ; enddo
    call calculate_density_derivs(t(:,1), s(:,1), p_ref, dr0_dt, dr0_ds, &
                                  is, ie-is+1, tv%eqn_of_state)
    call calculate_density_derivs(t(:,1), s(:,1), p_ref_cv, drcv_dt, drcv_ds, &
                                  is, ie-is+1, tv%eqn_of_state)
    do k=1,nz
      call calculate_density(t(:,k), s(:,k), p_ref, r0(:,k), is, ie-is+1, &
                             tv%eqn_of_state)
      call calculate_density(t(:,k), s(:,k), p_ref_cv, rcv(:,k), is, &
                             ie-is+1, tv%eqn_of_state)
    enddo
    if(id_clock_eos>0) call cpu_clock_end(id_clock_eos)

    if (cs%ML_resort) then
      if(id_clock_resort>0) call cpu_clock_begin(id_clock_resort)
      if (cs%ML_presort_nz_conv_adj > 0) &
        call convective_adjustment(h(:,1:), u, v, r0(:,1:), rcv(:,1:), t(:,1:), &
                                   s(:,1:), eps, d_eb, dke_ca, ctke, j, g, gv, cs, &
                                   cs%ML_presort_nz_conv_adj)

      call sort_ml(h(:,1:), r0(:,1:), eps, g, gv, cs, ksort)
      if(id_clock_resort>0) call cpu_clock_end(id_clock_resort)
    else
      do k=1,nz ; do i=is,ie ; ksort(i,k) = k ; enddo ; enddo

      if(id_clock_adjustment>0) call cpu_clock_begin(id_clock_adjustment)
      !  Undergo instantaneous entrainment into the buffer layers and mixed layers
      ! to remove hydrostatic instabilities.  Any water that is lighter than
      ! currently in the mixed or buffer layer is entrained.
      call convective_adjustment(h(:,1:), u, v, r0(:,1:), rcv(:,1:), t(:,1:), &
                                 s(:,1:), eps, d_eb, dke_ca, ctke, j, g, gv, cs)
      do i=is,ie ; h_ca(i) = h(i,1) ; enddo

      if(id_clock_adjustment>0) call cpu_clock_end(id_clock_adjustment)
    endif

    if (associated(fluxes%lrunoff) .and. cs%do_rivermix) then

      ! Here we add an additional source of TKE to the mixed layer where river
      ! is present to simulate unresolved estuaries. The TKE input is diagnosed
      ! as follows:
      !   TKE_river[m3 s-3] = 0.5*rivermix_depth*g*Irho0*drho_ds*
      !                       River*(Samb - Sriver) = CS%mstar*U_star^3
      ! where River is in units of m s-1.
      ! Samb = Ambient salinity at the mouth of the estuary
      ! rivermix_depth =  The prescribed depth over which to mix river inflow
      ! drho_ds = The gradient of density wrt salt at the ambient surface salinity.
      ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
      rmixconst = 0.5*cs%rivermix_depth*gv%g_Earth*irho0**2
      do i=is,ie
        tke_river(i) = max(0.0, rmixconst*dr0_ds(i)* &
            (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * s(i,1))
      enddo
    else
      do i=is,ie ; tke_river(i) = 0.0 ; enddo
    endif


    if(id_clock_conv>0) call cpu_clock_begin(id_clock_conv)

    ! The surface forcing is contained in the fluxes type.
    ! We aggregate the thermodynamic forcing for a time step into the following:
    ! netMassInOut = water (H units) added/removed via surface fluxes
    ! netMassOut   = water (H units) removed via evaporating surface fluxes
    ! net_heat     = heat (degC * H) via surface fluxes
    ! net_salt     = salt ( g(salt)/m2 for non-Bouss and ppt*m/s for Bouss ) via surface fluxes
    ! Pen_SW_bnd   = components to penetrative shortwave radiation
    call extractfluxes1d(g, gv, fluxes, optics, nsw, j, dt, &
                  cs%H_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
                  h(:,1:), t(:,1:), netmassinout, netmassout, net_heat, net_salt, pen_sw_bnd,&
                  tv, aggregate_fw_forcing)

    ! This subroutine causes the mixed layer to entrain to depth of free convection.
    call mixedlayer_convection(h(:,1:), d_eb, htot, ttot, stot, uhtot, vhtot, &
                               r0_tot, rcv_tot, u, v, t(:,1:), s(:,1:),       &
                               r0(:,1:), rcv(:,1:), eps,                      &
                               dr0_dt, drcv_dt, dr0_ds, drcv_ds,              &
                               netmassinout, netmassout, net_heat, net_salt,  &
                               nsw, pen_sw_bnd, opacity_band, conv_en,        &
                               dke_fc, j, ksort, g, gv, cs, tv, fluxes, dt,       &
                               aggregate_fw_forcing)

    if(id_clock_conv>0) call cpu_clock_end(id_clock_conv)

    !   Now the mixed layer undergoes mechanically forced entrainment.
    ! The mixed layer may entrain down to the Monin-Obukhov depth if the
    ! surface is becoming lighter, and is effectively detraining.

    !    First the TKE at the depth of free convection that is available
    !  to drive mixing is calculated.
    if(id_clock_mech>0) call cpu_clock_begin(id_clock_mech)

    call find_starting_tke(htot, h_ca, fluxes, conv_en, ctke, dke_fc, dke_ca, &
                           tke, tke_river, idecay_len_tke, cmke, dt, idt_diag, &
                           j, ksort, g, gv, cs)

    ! Here the mechanically driven entrainment occurs.
    call mechanical_entrainment(h(:,1:), d_eb, htot, ttot, stot, uhtot, vhtot, &
                                r0_tot, rcv_tot, u, v, t(:,1:), s(:,1:), r0(:,1:), rcv(:,1:), eps, dr0_dt, drcv_dt, &
                                cmke, idt_diag, nsw, pen_sw_bnd, opacity_band, tke, &
                                idecay_len_tke, j, ksort, g, gv, cs)

    call absorbremainingsw(g, gv, h(:,1:), opacity_band, nsw, j, dt, cs%H_limit_fluxes, &
                           cs%correct_absorption, cs%absorb_all_SW, &
                           t(:,1:), pen_sw_bnd, eps, ksort, htot, ttot)

    if (cs%TKE_diagnostics) then ; do i=is,ie
      cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) - idt_diag*tke(i)
    enddo ; endif
    if(id_clock_mech>0) call cpu_clock_end(id_clock_mech)

    ! Calculate the homogeneous mixed layer properties and store them in layer 0.
    do i=is,ie ; if (htot(i) > 0.0) then
      ih = 1.0 / htot(i)
      r0(i,0) = r0_tot(i) * ih ; rcv(i,0) = rcv_tot(i) * ih
      t(i,0) = ttot(i) * ih ; s(i,0) = stot(i) * ih
      h(i,0) = htot(i)
    else ! This may not ever be needed?
      t(i,0) = t(i,1) ; s(i,0) = s(i,1) ; r0(i,0) = r0(i,1) ; rcv(i,0) = rcv(i,1)
      h(i,0) = htot(i)
    endif ; enddo
    if (write_diags .and. ALLOCATED(cs%ML_depth)) then ; do i=is,ie
      cs%ML_depth(i,j) = h(i,0) * gv%H_to_m
    enddo ; endif
    if (ASSOCIATED(hml)) then ; do i=is,ie
      hml(i,j) = g%mask2dT(i,j) * (h(i,0) * gv%H_to_m)
    enddo ; endif

! At this point, return water to the original layers, but constrained to
! still be sorted.  After this point, all the water that is in massive
! interior layers will be denser than water remaining in the mixed- and
! buffer-layers.  To achieve this, some of these variable density layers
! might be split between two isopycnal layers that are denser than new
! mixed layer or any remaining water from the old mixed- or buffer-layers.
! Alternately, if there are fewer than nkbl of the old buffer or mixed layers
! with any mass, relatively light interior layers might be transferred to
! these unused layers (but not currently in the code).

    if (cs%ML_resort) then
      if(id_clock_resort>0) call cpu_clock_begin(id_clock_resort)
      call resort_ml(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), gv%Rlay, eps, &
                     d_ea, d_eb, ksort, g, gv, cs, dr0_dt, dr0_ds, drcv_dt, drcv_ds)
      if(id_clock_resort>0) call cpu_clock_end(id_clock_resort)
    endif

    if (cs%limit_det .or. (cs%id_Hsfc_max > 0) .or. (cs%id_Hsfc_min > 0)) then
      do i=is,ie ; hsfc(i) = h(i,0) ; enddo
      do k=1,nkmb ; do i=is,ie ; hsfc(i) = hsfc(i) + h(i,k) ; enddo ; enddo

      if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
        dhsfc = cs%lim_det_dH_sfc ; dhd = cs%lim_det_dH_bathy
        do i=is,ie
          h_nbr = min(dhsfc*max(hmbl_prev(i-1,j), hmbl_prev(i+1,j), &
                                hmbl_prev(i,j-1), hmbl_prev(i,j+1)), &
                      max(hmbl_prev(i-1,j) - dhd*min(h_sum(i,j),h_sum(i-1,j)), &
                          hmbl_prev(i+1,j) - dhd*min(h_sum(i,j),h_sum(i+1,j)), &
                          hmbl_prev(i,j-1) - dhd*min(h_sum(i,j),h_sum(i,j-1)), &
                          hmbl_prev(i,j+1) - dhd*min(h_sum(i,j),h_sum(i,j+1))) )

          hsfc_min(i,j) = gv%H_to_m*max(h(i,0), min(hsfc(i), h_nbr))

          if (cs%limit_det) max_bl_det(i) = max(0.0, hsfc(i)-h_nbr)
        enddo
      endif

      if (cs%id_Hsfc_max > 0) then ; do i=is,ie
        hsfc_max(i,j) = hsfc(i)*gv%H_to_m
      enddo ; endif
    endif

! Move water left in the former mixed layer into the buffer layer and
! from the buffer layer into the interior.  These steps might best be
! treated in conjuction.
    if(id_clock_detrain>0) call cpu_clock_begin(id_clock_detrain)
    if (cs%nkbl == 1) then
      call mixedlayer_detrain_1(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), &
                                gv%Rlay, dt, dt__diag, d_ea, d_eb, j, g, gv, cs, &
                                drcv_dt, drcv_ds, max_bl_det)
    elseif (cs%nkbl == 2) then
      call mixedlayer_detrain_2(h(:,0:), t(:,0:), s(:,0:), r0(:,0:), rcv(:,0:), &
                                gv%Rlay, dt, dt__diag, d_ea, j, g, gv, cs, &
                                dr0_dt, dr0_ds, drcv_dt, drcv_ds, max_bl_det)
    else ! CS%nkbl not = 1 or 2
      ! This code only works with 1 or 2 buffer layers.
      call mom_error(fatal, "MOM_mixed_layer: CS%nkbl must be 1 or 2 for now.")
    endif
    if(id_clock_detrain>0) call cpu_clock_end(id_clock_detrain)


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

! Now set the properties of the layers in the mixed layer in the original
! 3-d variables.
    if (cs%Resolve_Ekman .and. (cs%nkml>1)) then
      ! The thickness of the topmost piece of the mixed layer is given by
      ! h_1 = H / (3 + sqrt(|f|*H^2/2*nu_max)), which asymptotes to the Ekman
      ! layer depth and 1/3 of the mixed layer depth.  This curve has been
      ! determined to maximize the impact of the Ekman transport in the mixed
      ! layer TKE budget with nkml=2.  With nkml=3, this should also be used,
      ! as the third piece will then optimally describe mixed layer
      ! restratification.  For nkml>=4 the whole strategy should be revisited.
      do i=is,ie
        ku_star = 0.41*fluxes%ustar(i,j) ! Maybe could be replaced with u*+w*?
        if (associated(fluxes%ustar_shelf) .and. &
            associated(fluxes%frac_shelf_h)) then
          if (fluxes%frac_shelf_h(i,j) > 0.0) &
            ku_star = (1.0 - fluxes%frac_shelf_h(i,j)) * ku_star + &
                      fluxes%frac_shelf_h(i,j) * (0.41*fluxes%ustar_shelf(i,j))
        endif
        absf_x_h = 0.25 * h(i,0) * &
            ((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 the mixed layer vertical viscosity specification is changed in
        ! MOM_vert_friction.F90, this line will have to be modified accordingly.
        h_3d(i,j,1) = h(i,0) / (3.0 + sqrt(absf_x_h*(absf_x_h + 2.0*ku_star) / &
                                         (ku_star**2)) )
        do k=2,cs%nkml
          ! The other layers are evenly distributed through the mixed layer.
          h_3d(i,j,k) = (h(i,0)-h_3d(i,j,1)) * inkmlm1
          d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
          d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
        enddo
      enddo
    else
      do i=is,ie
        h_3d(i,j,1) = h(i,0) * inkml
      enddo
      do k=2,cs%nkml ; do i=is,ie
        h_3d(i,j,k) = h(i,0) * inkml
        d_ea(i,k) = d_ea(i,k) + h_3d(i,j,k)
        d_ea(i,1) = d_ea(i,1) - h_3d(i,j,k)
      enddo ; enddo
    endif
    do i=is,ie ; h(i,0) = 0.0 ; enddo
    do k=1,cs%nkml ; do i=is,ie
      tv%T(i,j,k) = t(i,0) ; tv%S(i,j,k) = s(i,0)
    enddo ; enddo

    ! These sum needs to be done in the original layer space.

    ! The treatment of layer 1 is atypical because evaporation shows up as
    ! negative ea(i,1), and because all precipitation goes straight into layer 1.
    ! The code is ordered so that any roundoff errors in ea are lost the surface.
!    do i=is,ie ; eaml(i,1) = 0.0 ; enddo
!    do k=2,nz ; do i=is,ie ; eaml(i,k) = eaml(i,k-1) - d_ea(i,k-1) ; enddo ; enddo
!    do i=is,ie ; eaml(i,1) = netMassInOut(i) ; enddo


    do i=is,ie
! eaml(i,nz) is derived from  h(i,nz) - h_orig(i,nz) = eaml(i,nz) - ebml(i,nz-1)
      ebml(i,nz) = 0.0
      eaml(i,nz) = (h(i,nz) - h_orig(i,nz)) - d_eb(i,nz)
    enddo
    do k=nz-1,1,-1 ; do i=is,ie
      ebml(i,k) = ebml(i,k+1) - d_eb(i,k+1)
      eaml(i,k) = eaml(i,k+1) + d_ea(i,k)
    enddo ; enddo
    do i=is,ie ; eaml(i,1) = netmassinout(i) ; enddo

    ! Copy the interior thicknesses and other fields back to the 3-d arrays.
    do k=cs%nkml+1,nz ; do i=is,ie
      h_3d(i,j,k) = h(i,k); tv%T(i,j,k) = t(i,k) ; tv%S(i,j,k) = s(i,k)
    enddo ; enddo

    do k=1,nz ; do i=is,ie
      ea(i,j,k) = ea(i,j,k) + eaml(i,k)
      eb(i,j,k) = eb(i,j,k) + ebml(i,k)
    enddo ; enddo

    if (cs%id_h_mismatch > 0) then
      do i=is,ie
        h_miss(i,j) = abs(h_3d(i,j,1) - (h_orig(i,1) + &
          (eaml(i,1) + (ebml(i,1) - eaml(i,1+1)))))
      enddo
      do k=2,nz-1 ; do i=is,ie
        h_miss(i,j) = h_miss(i,j) + abs(h_3d(i,j,k) - (h_orig(i,k) + &
          ((eaml(i,k) - ebml(i,k-1)) + (ebml(i,k) - eaml(i,k+1)))))
      enddo ; enddo
      do i=is,ie
        h_miss(i,j) = h_miss(i,j) + abs(h_3d(i,j,nz) - (h_orig(i,nz) + &
          ((eaml(i,nz) - ebml(i,nz-1)) + ebml(i,nz))))
        h_miss(i,j) = gv%H_to_m * h_miss(i,j)
      enddo
    endif

  enddo ! j loop

  ! Whenever thickness changes let the diag manager know, target grids
  ! for vertical remapping may need to be regenerated.
  ! This needs to happen after the H update and before the next post_data.
  call diag_update_remap_grids(cs%diag)

!$OMP end parallel

  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_RiBulk > 0) &
      call post_data(cs%id_TKE_RiBulk, cs%diag_TKE_RiBulk, 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_pen_SW > 0) &
      call post_data(cs%id_TKE_pen_SW, cs%diag_TKE_pen_SW, 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_TKE_conv_s2 > 0) &
      call post_data(cs%id_TKE_conv_s2, cs%diag_TKE_conv_s2, cs%diag)
    if (cs%id_PE_detrain > 0) &
      call post_data(cs%id_PE_detrain, cs%diag_PE_detrain, cs%diag)
    if (cs%id_PE_detrain2 > 0) &
      call post_data(cs%id_PE_detrain2, cs%diag_PE_detrain2, cs%diag)
    if (cs%id_h_mismatch > 0) &
      call post_data(cs%id_h_mismatch, h_miss, cs%diag)
    if (cs%id_Hsfc_used > 0) &
      call post_data(cs%id_Hsfc_used, hsfc_used, cs%diag)
    if (cs%id_Hsfc_max > 0) &
      call post_data(cs%id_Hsfc_max, hsfc_max, cs%diag)
    if (cs%id_Hsfc_min > 0) &
      call post_data(cs%id_Hsfc_min, hsfc_min, cs%diag)
  endif

end subroutine bulkmixedlayer

!>   This subroutine does instantaneous convective entrainment into the buffer
!! layers and mixed layers to remove hydrostatic instabilities.  Any water that
!! is lighter than currently in the mixed- or buffer- layer is entrained.
subroutine convective_adjustment(h, u, v, R0, Rcv, T, S, eps, d_eb, &
                                 dKE_CA, cTKE, j, G, GV, CS, nz_conv)
  type(ocean_grid_type),             intent(in)    :: G    !< The ocean's grid structure.
  type(verticalgrid_type),           intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: h    !< 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),SZK_(GV)), intent(inout) :: u    !< Zonal velocities interpolated to h
                                                           !! points, m s-1.
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: v    !< Zonal velocities interpolated to h
                                                           !! points, m s-1.
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: T    !< Layer temperatures, in deg C.
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: S
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: R0   !< Potential density referenced to
                                                           !! surface pressure, in kg m-3.
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: Rcv  !< The coordinate defining potential
                                                           !! density, in kg m-3.
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer
                                                           !! in the entrainment from below, in H.
                                                           !! Positive values go with mass gain by
                                                           !! a layer.
  real, dimension(SZI_(G),SZK_(GV)), intent(in)    :: eps
  real, dimension(SZI_(G),SZK_(GV)), intent(out)   :: dKE_CA !< The vertically integrated change in
                                                           !! kinetic energy due to convective
                                                           !! adjustment, in m3 s-2.
  real, dimension(SZI_(G),SZK_(GV)), intent(out)   :: cTKE !< The buoyant turbulent kinetic energy
                                                           !! source due to convective adjustment,
                                                           !! in m3 s-2.
  integer,                           intent(in)    :: j    !< The j-index to work on.
  type(bulkmixedlayer_cs),           pointer       :: CS   !< The control structure for this module.
  integer,                 optional, intent(in)    :: nz_conv !< If present, the number of layers
                                                           !! over which to do convective adjustment
                                                           !! (perhaps CS%nkml).

!   This subroutine does instantaneous convective entrainment into the buffer
! layers and mixed layers to remove hydrostatic instabilities.  Any water that
! is lighter than currently in the mixed- or buffer- layer is entrained.

! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)  The units
!                of h are referred to as H below.
!  (in/out)  u - Zonal velocities interpolated to h points, m s-1.
!  (in/out)  v - Zonal velocities interpolated to h points, m s-1.
!  (in/out)  R0 - Potential density referenced to surface pressure, in kg m-3.
!  (in/out)  Rcv - The coordinate defining potential density, in kg m-3.
!  (in/out)  T - Layer temperatures, in deg C.
!  (in/out)  S - Layer salinities, in psu.
!  (in/out)  d_eb - The downward increase across a layer in the entrainment from
!                   below, in H. Positive values go with mass gain by a layer.
!  (out)     dKE_CA - The vertically integrated change in kinetic energy due
!                     to convective adjustment, in m3 s-2.
!  (out)     cTKE - The buoyant turbulent kinetic energy source due to
!                   convective adjustment, in m3 s-2.
!  (in)      j - The j-index to work on.
!  (in)      ksort - The density-sorted k-indicies.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure for this module.
!  (in,opt)  nz_conv - If present, the number of layers over which to do
!                      convective adjustment (perhaps CS%nkml).
  real, dimension(SZI_(G)) :: &
    htot, &     !   The total depth of the layers being considered for
                ! entrainment, in H.
    r0_tot, &   !   The integrated potential density referenced to the surface
                ! of the layers which are fully entrained, in H kg m-3.
    rcv_tot, &  !   The integrated coordinate value potential density of the
                ! layers that are fully entrained, in H kg m-3.
    ttot, &     !   The integrated temperature of layers which are fully
                ! entrained, in H K.
    stot, &     !   The integrated salt of layers which are fully entrained,
                ! in H PSU.
    uhtot, &    !   The depth integrated zonal and meridional velocities in
    vhtot, &    ! the mixed layer, in H m s-1.
    ke_orig, &  !   The total mean kinetic energy in the mixed layer before
                ! convection, H m2 s-2.
    h_orig_k1   !   The depth of layer k1 before convective adjustment, in H.
  real :: h_ent !   The thickness from a layer that is entrained, in H.
  real :: Ih    !   The inverse of a thickness, in H-1.
  real :: g_H2_2Rho0  ! Half the gravitational acceleration times the
                      ! square of the conversion from H to m divided
                      ! by the mean density, in m6 s-2 H-2 kg-1.
  integer :: is, ie, nz, i, k, k1, nzc, nkmb

  is = g%isc ; ie = g%iec ; nz = gv%ke
  g_h2_2rho0 = (gv%g_Earth * gv%H_to_m**2) / (2.0 * gv%Rho0)
  nzc = nz ; if (present(nz_conv)) nzc = nz_conv
  nkmb = cs%nkml+cs%nkbl

!   Undergo instantaneous entrainment into the buffer layers and mixed layers
! to remove hydrostatic instabilities.  Any water that is lighter than currently
! in the layer is entrained.
  do k1=min(nzc-1,nkmb),1,-1
    do i=is,ie
      h_orig_k1(i) = h(i,k1)
      ke_orig(i) = 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)
      uhtot(i) = h(i,k1)*u(i,k1) ; vhtot(i) = h(i,k1)*v(i,k1)
      r0_tot(i) = r0(i,k1) * h(i,k1)
      ctke(i,k1) = 0.0 ; dke_ca(i,k1) = 0.0

      rcv_tot(i) = rcv(i,k1) * h(i,k1)
      ttot(i) = t(i,k1) * h(i,k1) ; stot(i) = s(i,k1) * h(i,k1)
    enddo
    do k=k1+1,nzc
      do i=is,ie
        if ((h(i,k) > eps(i,k)) .and. (r0_tot(i) > h(i,k1)*r0(i,k))) then
          h_ent = h(i,k)-eps(i,k)
          ctke(i,k1) = ctke(i,k1) + (h_ent * g_h2_2rho0 * &
                   (r0_tot(i) - h(i,k1)*r0(i,k)) * cs%nstar2)
          if (k < nkmb) then
            ctke(i,k1) = ctke(i,k1) + ctke(i,k)
            dke_ca(i,k1) = dke_ca(i,k1) + dke_ca(i,k)
          endif
          r0_tot(i) = r0_tot(i) + h_ent * r0(i,k)
          ke_orig(i) = ke_orig(i) + 0.5*h_ent* &
              (u(i,k)*u(i,k) + v(i,k)*v(i,k))
          uhtot(i) = uhtot(i) + h_ent*u(i,k)
          vhtot(i) = vhtot(i) + h_ent*v(i,k)

          rcv_tot(i) = rcv_tot(i) + h_ent * rcv(i,k)
          ttot(i) = ttot(i) + h_ent * t(i,k)
          stot(i) = stot(i) + h_ent * s(i,k)
          h(i,k1) = h(i,k1) + h_ent ; h(i,k) = eps(i,k)

          d_eb(i,k) = d_eb(i,k) - h_ent
          d_eb(i,k1) = d_eb(i,k1) + h_ent
        endif
      enddo
    enddo
! Determine the temperature, salinity, and velocities of the mixed or buffer
! layer in question, if it has entrained.
    do i=is,ie ; if (h(i,k1) > h_orig_k1(i)) then
      ih = 1.0 / h(i,k1)
      r0(i,k1) = r0_tot(i) * ih
      u(i,k1) = uhtot(i) * ih ; v(i,k1) = vhtot(i) * ih
      dke_ca(i,k1) = dke_ca(i,k1) + gv%H_to_m * (cs%bulk_Ri_convective * &
           (ke_orig(i) - 0.5*h(i,k1)*(u(i,k1)**2 + v(i,k1)**2)))
      rcv(i,k1) = rcv_tot(i) * ih
      t(i,k1) = ttot(i) * ih ; s(i,k1) = stot(i) * ih
    endif ; enddo
  enddo
! If lower mixed or buffer layers are massless, give them the properties of the
! layer above.
  do k=2,min(nzc,nkmb) ; do i=is,ie ; if (h(i,k) == 0.0) then
    r0(i,k) = r0(i,k-1)
    rcv(i,k) = rcv(i,k-1) ; t(i,k) = t(i,k-1) ; s(i,k) = s(i,k-1)
  endif ; enddo ; enddo

end subroutine convective_adjustment

!>   This subroutine causes the mixed layer to entrain to the depth of free
!! convection.  The depth of free convection is the shallowest depth at which the
!! fluid is denser than the average of the fluid above.
subroutine mixedlayer_convection(h, d_eb, htot, Ttot, Stot, uhtot, vhtot,      &
                                 R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps,    &
                                 dR0_dT, dRcv_dT, dR0_dS, dRcv_dS,             &
                                 netMassInOut, netMassOut, Net_heat, Net_salt, &
                                 nsw, Pen_SW_bnd, opacity_band, Conv_en,       &
                                 dKE_FC, j, ksort, G, GV, CS, tv, fluxes, dt,      &
                                 aggregate_FW_forcing)
  type(ocean_grid_type),             intent(in)    :: G       !< The ocean's grid structure.
  type(verticalgrid_type),           intent(in)    :: GV      !< The ocean's vertical grid
                                                              !! structure.
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: h       !< 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),SZK_(GV)), intent(inout) :: d_eb    !< The downward increase across a
                                                              !! layer in the entrainment from below
                                                              !! , in H. Positive values go with
                                                              !! mass gain by a layer.
  real, dimension(SZI_(G)),          intent(out)   :: htot    !< The accumulated mixed layer
                                                              !! thickness, in H.
  real, dimension(SZI_(G)),          intent(out)   :: Ttot    !< The depth integrated mixed layer
                                                              !! temperature, in deg C H.
  real, dimension(SZI_(G)),          intent(out)   :: Stot    !< The depth integrated mixed layer
                                                              !! salinity, in psu H.
  real, dimension(SZI_(G)),          intent(out)   :: uhtot   !< The depth integrated mixed layer
                                                              !! zonal velocity, H m s-1.
  real, dimension(SZI_(G)),          intent(out)   :: vhtot   !< The integrated mixed layer
                                                              !! meridional velocity, H m s-1.
  real, dimension(SZI_(G)),          intent(out)   :: R0_tot  !< The integrated mixed layer
                                                              !! potential density referenced to 0
                                                              !! pressure, in kg m-2.
  real, dimension(SZI_(G)),          intent(out)   :: Rcv_tot !< The integrated mixed layer
                                                              !! coordinate variable potential
                                                              !! density, in kg m-2.
  real, dimension(SZI_(G),SZK_(GV)), intent(in)    :: u, v, T, S, R0, Rcv, eps
  real, dimension(SZI_(G)),          intent(in)    :: dR0_dT, dRcv_dT, dR0_dS, dRcv_dS
  real, dimension(SZI_(G)),          intent(in)    :: netMassInOut, netMassOut
  real, dimension(SZI_(G)),          intent(in)    :: Net_heat, Net_salt
  integer,                           intent(in)    :: nsw     !< The number of bands of penetrating
                                                              !! shortwave radiation.
  real, dimension(:,:),              intent(inout) :: Pen_SW_bnd !< The penetrating shortwave
                                                              !! heating at the sea surface in each
                                                              !! penetrating band, in K H,
                                                              !! size nsw x SZI_(G).
  real, dimension(:,:,:),            intent(in)    :: opacity_band
  real, dimension(SZI_(G)),          intent(out)   :: Conv_en !< The buoyant turbulent kinetic
                                                              !! energy source due to free
                                                              !! convection, in m3 s-2.
  real, dimension(SZI_(G)),          intent(out)   :: dKE_FC  !< The vertically integrated change
                                                              !! in kinetic energy due to free
                                                              !! convection, in m3 s-2.
  integer,                           intent(in)    :: j       !< The j-index to work on.
  integer, dimension(SZI_(G),SZK_(GV)), &
                                     intent(in)    :: ksort   !< The density-sorted k-indices.
  type(bulkmixedlayer_cs),           pointer       :: CS      !< The control structure for this
                                                              !! module.
  type(thermo_var_ptrs),             intent(inout) :: tv
  type(forcing),                     intent(inout) :: fluxes
  real,                              intent(in)    :: dt
  logical,                           intent(in)    :: aggregate_FW_forcing

!   This subroutine causes the mixed layer to entrain to the depth of free
! convection.  The depth of free convection is the shallowest depth at which the
! fluid is denser than the average of the fluid above.
! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)  The units
!                of h are referred to as H below.
!  (in/out)  d_eb - The downward increase across a layer in the entrainment from
!                   below, in H. Positive values go with mass gain by a layer.
!  (out)     htot - The accumulated mixed layer thickness, in H.
!  (out)     Ttot - The depth integrated mixed layer temperature, in deg C H.
!  (out)     Stot - The depth integrated mixed layer salinity, in psu H.
!  (out)     uhtot - The depth integrated mixed layer zonal velocity, H m s-1.
!  (out)     vhtot - The integrated mixed layer meridional velocity, H m s-1.
!  (out)     R0_tot - The integrated mixed layer potential density referenced
!                     to 0 pressure, in kg m-2.
!  (out)     Rcv_tot - The integrated mixed layer coordinate variable
!                      potential density, in kg m-2.
!  (in)      nsw - The number of bands of penetrating shortwave radiation.
!  (out)     Pen_SW_bnd - The penetrating shortwave heating at the sea surface
!                         in each penetrating band, in K H, size nsw x SZI_(G).
!  (out)     Conv_en - The buoyant turbulent kinetic energy source due to
!                      free convection, in m3 s-2.
!  (out)     dKE_FC - The vertically integrated change in kinetic energy due
!                     to free convection, in m3 s-2.
!  (in)      j - The j-index to work on.
!  (in)      ksort - The density-sorted k-indices.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure for this module.

  real, dimension(SZI_(G)) :: &
    massOutRem, &      !   Evaporation that remains to be supplied, in H.
    netMassIn          ! mass entering through ocean surface (H)
  real :: SW_trans     !   The fraction of shortwave radiation
                       ! that is not absorbed in a layer, ND.
  real :: Pen_absorbed !   The amount of penetrative shortwave radiation
                       ! that is absorbed in a layer, in units of K H.
  real :: h_avail      !   The thickness in a layer available for
                       ! entrainment, in H.
  real :: h_ent        !   The thickness from a layer that is entrained, in H.
  real :: T_precip     !   The temperature of the precipitation, in deg C.
  real :: C1_3, C1_6   !  1/3 and 1/6.
  real :: En_fn, Frac, x1 !  Nondimensional temporary variables.
  real :: dr, dr0      ! Temporary variables with units of kg m-3 H.
  real :: dr_ent, dr_comp ! Temporary variables with units of kg m-3 H.
  real :: dr_dh        ! The partial derivative of dr_ent with h_ent, in kg m-3.
  real :: h_min, h_max !   The minimum, maximum, and previous estimates for
  real :: h_prev       ! h_ent, in H.
  real :: h_evap       !   The thickness that is evaporated, in H.
  real :: dh_Newt      !   The Newton's method estimate of the change in
                       ! h_ent between iterations, in H.
  real :: g_H2_2Rho0   !   Half the gravitational acceleration times the
                       ! square of the conversion from H to m divided
                       ! by the mean density, in m6 s-2 H-2 kg-1.
  real :: Angstrom     !   The minimum layer thickness, in H.
  real :: opacity      !   The opacity converted to units of H-1.
  real :: sum_Pen_En   !   The potential energy change due to penetrating
                       ! shortwave radiation, integrated over a layer, in
                       ! H kg m-3.
  real :: Idt          ! 1.0/dt
  real :: netHeatOut   ! accumulated heat content of mass leaving ocean
  integer :: is, ie, nz, i, k, ks, itt, n
  real, dimension(max(nsw,1)) :: &
    C2, &              ! Temporary variable with units of kg m-3 H-1.
    r_SW_top           ! Temporary variables with units of H kg m-3.

  angstrom = gv%Angstrom
  c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0
  g_h2_2rho0 = (gv%g_Earth * gv%H_to_m**2) / (2.0 * gv%Rho0)
  idt        = 1.0/dt
  is = g%isc ; ie = g%iec ; nz = gv%ke

  do i=is,ie ; if (ksort(i,1) > 0) then
    k = ksort(i,1)

    if (aggregate_fw_forcing) then
      massoutrem(i) = 0.0
      if (netmassinout(i) < 0.0) massoutrem(i) = -netmassinout(i)
      netmassin(i) = netmassinout(i) + massoutrem(i)
    else
      massoutrem(i) = -netmassout(i)
      netmassin(i)  = netmassinout(i) - netmassout(i)
    endif

    ! htot is an Angstrom (taken from layer 1) plus any net precipitation.
    h_ent     = max(min(angstrom,h(i,k)-eps(i,k)),0.0)
    htot(i)   = h_ent + netmassin(i)
    h(i,k)    = h(i,k) - h_ent
    d_eb(i,k) = d_eb(i,k) - h_ent

    pen_absorbed = 0.0
    do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
      sw_trans        = exp(-htot(i)*opacity_band(n,i,k))
      pen_absorbed    = pen_absorbed + pen_sw_bnd(n,i) * (1.0-sw_trans)
      pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
    endif ; enddo

    ! Precipitation is assumed to have the same temperature and velocity
    ! as layer 1.  Because layer 1 might not be the topmost layer, this
    ! involves multiple terms.
    t_precip = t(i,1)
    ttot(i) = (net_heat(i) + (netmassin(i) * t_precip + h_ent * t(i,k))) + &
              pen_absorbed
    ! Net_heat contains both heat fluxes and the heat content of mass fluxes.
 !! Ttot(i) = netMassIn(i) * T_precip + h_ent * T(i,k)
 !! Ttot(i) = Net_heat(i) + Ttot(i)
 !! Ttot(i) = Ttot(i) + Pen_absorbed
  ! smg:
  ! Ttot(i)   = (Net_heat(i) + (h_ent * T(i,k))) + Pen_absorbed
    stot(i)   = h_ent*s(i,k) + net_salt(i)
    uhtot(i)  = u(i,1)*netmassin(i) + u(i,k)*h_ent
    vhtot(i)  = v(i,1)*netmassin(i) + v(i,k)*h_ent
    r0_tot(i) = (h_ent*r0(i,k) + netmassin(i)*r0(i,1)) + &
!                   dR0_dT(i)*netMassIn(i)*(T_precip - T(i,1)) + &
                (dr0_dt(i)*(net_heat(i) + pen_absorbed) - &
                 dr0_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
    rcv_tot(i) = (h_ent*rcv(i,k) + netmassin(i)*rcv(i,1)) + &
!                    dRcv_dT(i)*netMassIn(i)*(T_precip - T(i,1)) + &
                 (drcv_dt(i)*(net_heat(i) + pen_absorbed) - &
                  drcv_ds(i) * (netmassin(i) * s(i,1) - net_salt(i)))
    conv_en(i) = 0.0 ; dke_fc(i) = 0.0
    if(ASSOCIATED(fluxes%heat_content_massin))                            &
           fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) &
                       + t_precip * netmassin(i) * gv%H_to_kg_m2 * fluxes%C_p * idt
    if (ASSOCIATED(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
                         t_precip * netmassin(i) * gv%H_to_kg_m2
  endif ; enddo

  ! Now do netMassOut case in this block.
  ! At this point htot contains an Angstrom of fluid from layer 0 plus netMassIn.
  do ks=1,nz
    do i=is,ie ; if (ksort(i,ks) > 0) then
      k = ksort(i,ks)

      if ((htot(i) < angstrom) .and. (h(i,k) > eps(i,k))) then
        ! If less than an Angstrom was available from the layers above plus
        ! any precipitation, add more fluid from this layer.
        h_ent = min(angstrom-htot(i), h(i,k)-eps(i,k))
        htot(i) = htot(i) + h_ent
        h(i,k) = h(i,k) - h_ent
        d_eb(i,k) = d_eb(i,k) - h_ent

        r0_tot(i) = r0_tot(i) + h_ent*r0(i,k)
        uhtot(i) = uhtot(i) + h_ent*u(i,k)
        vhtot(i) = vhtot(i) + h_ent*v(i,k)

        rcv_tot(i) = rcv_tot(i) + h_ent*rcv(i,k)
        ttot(i) = ttot(i) + h_ent*t(i,k)
        stot(i) = stot(i) + h_ent*s(i,k)
      endif

      ! Water is removed from the topmost layers with any mass.
      ! We may lose layers if they are thin enough.
      ! The salt that is left behind goes into Stot.
      if ((massoutrem(i) > 0.0) .and. (h(i,k) > eps(i,k))) then
        if (massoutrem(i) > (h(i,k) - eps(i,k))) then
          h_evap = h(i,k) - eps(i,k)
          h(i,k) = eps(i,k)
          massoutrem(i) = massoutrem(i) - h_evap
        else
          h_evap = massoutrem(i)
          h(i,k) = h(i,k) - h_evap
          massoutrem(i) = 0.0
        endif

        stot(i) = stot(i) + h_evap*s(i,k)
        r0_tot(i) = r0_tot(i) + dr0_ds(i)*h_evap*s(i,k)
        rcv_tot(i) = rcv_tot(i) + drcv_ds(i)*h_evap*s(i,k)
        d_eb(i,k) = d_eb(i,k) - h_evap

        ! smg: when resolve the A=B code, we will set
        ! heat_content_massout = heat_content_massout - T(i,k)*h_evap*GV%H_to_kg_m2*fluxes%C_p*Idt
        ! by uncommenting the lines here.
        ! we will also then completely remove TempXpme from the model.
        if(ASSOCIATED(fluxes%heat_content_massout))                            &
           fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) &
                                    - t(i,k)*h_evap*gv%H_to_kg_m2 * fluxes%C_p * idt
        if (ASSOCIATED(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) - &
                                      t(i,k)*h_evap*gv%H_to_kg_m2

      endif

      ! The following section calculates how much fluid will be entrained.
      h_avail = h(i,k) - eps(i,k)
      if (h_avail > 0.0) then
        dr = r0_tot(i) - htot(i)*r0(i,k)
        h_ent = 0.0

        dr0 = dr
        do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
          dr0 = dr0 - (dr0_dt(i)*pen_sw_bnd(n,i)) * &
                      opacity_band(n,i,k)*htot(i)
        endif ; enddo

        ! Some entrainment will occur from this layer.
        if (dr0 > 0.0) then
          dr_comp = dr
          do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
            !   Compare the density at the bottom of a layer with the
            ! density averaged over the mixed layer and that layer.
            opacity = opacity_band(n,i,k)
            sw_trans = exp(-h_avail*opacity)
            dr_comp = dr_comp + (dr0_dt(i)*pen_sw_bnd(n,i)) * &
                ((1.0 - sw_trans) - opacity*(htot(i)+h_avail)*sw_trans)
          endif ; enddo
          if (dr_comp >= 0.0) then
            ! The entire layer is entrained.
            h_ent = h_avail
          else
            !  The layer is partially entrained.   Iterate to determine how much
            !  entrainment occurs.  Solve for the h_ent at which dr_ent = 0.

            ! Instead of assuming that the curve is linear between the two end
            ! points, assume that the change is concentrated near small values
            ! of entrainment.  On average, this saves about 1 iteration.
            frac = dr0 / (dr0 - dr_comp)
            h_ent = h_avail * frac*frac
            h_min = 0.0 ; h_max = h_avail

            do n=1,nsw
              r_sw_top(n) = dr0_dt(i) * pen_sw_bnd(n,i)
              c2(n) = r_sw_top(n) * opacity_band(n,i,k)**2
            enddo
            do itt=1,10
              dr_ent = dr ; dr_dh = 0.0
              do n=1,nsw
                opacity = opacity_band(n,i,k)
                sw_trans = exp(-h_ent*opacity)
                dr_ent = dr_ent + r_sw_top(n) * ((1.0 - sw_trans) - &
                           opacity*(htot(i)+h_ent)*sw_trans)
                dr_dh = dr_dh + c2(n) * (htot(i)+h_ent) * sw_trans
              enddo

              if (dr_ent > 0.0) then
                h_min = h_ent
              else
                h_max = h_ent
              endif

              dh_newt = -dr_ent / dr_dh
              h_prev = h_ent ; h_ent = h_prev+dh_newt
              if (h_ent > h_max) then
                h_ent = 0.5*(h_prev+h_max)
              else if (h_ent < h_min) then
                h_ent = 0.5*(h_prev+h_min)
              endif

              if (abs(dh_newt) < 0.2*angstrom) exit
            enddo

          endif

          !  Now that the amount of entrainment (h_ent) has been determined,
          !  calculate changes in various terms.
          sum_pen_en = 0.0 ; pen_absorbed = 0.0
          do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
            opacity = opacity_band(n,i,k)
            sw_trans = exp(-h_ent*opacity)

            x1 = h_ent*opacity
            if (x1 < 2.0e-5) then
              en_fn = (opacity*htot(i)*(1.0 - 0.5*(x1 - c1_3*x1)) + &
                       x1*x1*c1_6)
            else
              en_fn = ((opacity*htot(i) + 2.0) * &
                       ((1.0-sw_trans) / x1) - 1.0 + sw_trans)
            endif
            sum_pen_en = sum_pen_en - (dr0_dt(i)*pen_sw_bnd(n,i)) * en_fn

            pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
            pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
          endif ; enddo

          conv_en(i) = conv_en(i) + g_h2_2rho0 * h_ent * &
                       ( (r0_tot(i) - r0(i,k)*htot(i)) + sum_pen_en )

          r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
          stot(i) = stot(i) + h_ent * s(i,k)
          ttot(i) = ttot(i) + (h_ent * t(i,k) + pen_absorbed)
          rcv_tot(i) = rcv_tot(i) + (h_ent * rcv(i,k) + pen_absorbed*drcv_dt(i))
        endif ! dr0 > 0.0

        if (h_ent > 0.0) then
          if (htot(i) > 0.0) &
            dke_fc(i) = dke_fc(i) + cs%bulk_Ri_convective * 0.5 * &
              ((gv%H_to_m*h_ent) / (htot(i)*(h_ent+htot(i)))) * &
              ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)

          htot(i)  = htot(i)  + h_ent
          h(i,k) = h(i,k) - h_ent
          d_eb(i,k) = d_eb(i,k) - h_ent
          uhtot(i) = u(i,k)*h_ent ; vhtot(i) = v(i,k)*h_ent
        endif


      endif ! h_avail>0
    endif ; enddo ! i loop
  enddo ! k loop

end subroutine mixedlayer_convection

!>   This subroutine determines the TKE available at the depth of free
!! convection to drive mechanical entrainment.
subroutine find_starting_tke(htot, h_CA, fluxes, Conv_En, cTKE, dKE_FC, dKE_CA, &
                             TKE, TKE_river, Idecay_len_TKE, cMKE, dt, Idt_diag, &
                             j, ksort, G, GV, CS)
  type(ocean_grid_type),      intent(in)    :: G       !< The ocean's grid structure.
  type(verticalgrid_type),    intent(in)    :: GV      !< The ocean's vertical grid structure.
  real, dimension(SZI_(G)),   intent(in)    :: htot    !< The accumlated mixed layer thickness, in m
                                                       !! or kg m-2. (Intent in).
  real, dimension(SZI_(G)),   intent(in)    :: h_CA    !< The mixed layer depth after convective
                                                       !! adjustment, in H.
  type(forcing),              intent(in)    :: fluxes  !< A structure containing pointers to any
                                                       !! possible forcing fields.  Unused fields
                                                       !! have NULL ptrs.
  real, dimension(SZI_(G)),   intent(inout) :: Conv_En !< The buoyant turbulent kinetic energy
                                                       !! source due to free convection,
                                                       !! in m3 s-2.
  real, dimension(SZI_(G)),   intent(in)    :: dKE_FC  !< The vertically integrated change in
                                                       !! kinetic energy due to free convection,
                                                       !! in m3 s-2.
  real, dimension(SZI_(G),SZK_(GV)), &
                              intent(in)    :: cTKE    !< The buoyant turbulent kinetic energy
                                                       !! source due to convective adjustment,
                                                       !! in m3 s-2.
  real, dimension(SZI_(G),SZK_(GV)), &
                              intent(in)    :: dKE_CA  !< The vertically integrated change in
                                                       !! kinetic energy due to convective
                                                       !! adjustment, in m3 s-2.
  real, dimension(SZI_(G)),   intent(out)   :: TKE     !< The turbulent kinetic energy available for
                                                       !! mixing over a time step, in m3 s-2.
  real, dimension(SZI_(G)),   intent(out)   :: Idecay_len_TKE !< The inverse of the vertical decay
                                                       !! scale for TKE, in H-1.
  real, dimension(SZI_(G)),   intent(in)    :: TKE_river
  real, dimension(2,SZI_(G)), intent(out)   :: cMKE    !< Coefficients of HpE and HpE^2 in
                                                       !! calculating the denominator of MKE_rate,
                                                       !! in H-1 and H-2.
  real,                       intent(in)    :: dt      !< The time step in s.
  real,                       intent(in)    :: Idt_diag
  integer,                    intent(in)    :: j       !< The j-index to work on.
  integer, dimension(SZI_(G),SZK_(GV)), &
                              intent(in)    :: ksort   !< The density-sorted k-indicies.
  type(bulkmixedlayer_cs),    pointer       :: CS      !< The control structure for this module.

!   This subroutine determines the TKE available at the depth of free
! convection to drive mechanical entrainment.

! Arguments: htot - The accumlated mixed layer thickness, in m or kg m-2. (Intent in)
!                   The units of htot are referred to as H below.
!  (in)      h_CA - The mixed layer depth after convective adjustment, in H.
!  (in)      fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields have NULL ptrs.
!  (in)      Conv_en - The buoyant turbulent kinetic energy source due to
!                      free convection, in m3 s-2.
!  (in)      cTKE - The buoyant turbulent kinetic energy source due to
!                   convective adjustment, in m3 s-2.
!  (in)      dKE_FC - The vertically integrated change in kinetic energy due
!                     to free convection, in m3 s-2.
!  (in)      dKE_CA - The vertically integrated change in kinetic energy due
!                     to convective adjustment, in m3 s-2.
!  (out)     TKE - The turbulent kinetic energy available for mixing over a
!                  time step, in m3 s-2.
!  (out)     Idecay_len_TKE - The inverse of the vertical decay scale for
!                             TKE, in H-1.
!  (out)     cMKE - Coefficients of HpE and HpE^2 in calculating the
!                   denominator of MKE_rate, in H-1 and H-2.
!  (in)      dt - The time step in s.
!  (in)      j - The j-index to work on.
!  (in)      ksort - The density-sorted k-indicies.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure for this module.
  real :: dKE_conv  ! The change in mean kinetic energy due
                    ! to all convection, in m3 s-2.
  real :: nstar_FC  ! The effective efficiency with which the energy released by
                    ! free convection is converted to TKE, often ~0.2, ND.
  real :: nstar_CA  ! The effective efficiency with which the energy released by
                    ! convective adjustment is converted to TKE, often ~0.2, ND.
  real :: TKE_CA    ! The potential energy released by convective adjustment if
                    ! that release is positive, in m3 s2.
  real :: MKE_rate_CA ! MKE_rate for convective adjustment, ND, 0 to 1.
  real :: MKE_rate_FC ! MKE_rate for free convection, ND, 0 to 1.
  real :: totEn     ! The total potential energy released by convection, m3 s-2.
  real :: Ih        ! The inverse of a thickness, in H-1.
  real :: exp_kh    ! The nondimensional decay of TKE across a layer, ND.
  real :: absf      ! The absolute value of f averaged to thickness points, s-1.
  real :: U_star    ! The friction velocity in m s-1.
  real :: absf_Ustar  ! The absolute value of f divided by U_star, in m-1.
  real :: wind_TKE_src ! The surface wind source of TKE, in m3 s-3.
  real :: diag_wt   ! The ratio of the current timestep to the diagnostic
                    ! timestep (which may include 2 calls), ND.
  integer :: is, ie, nz, i

  is = g%isc ; ie = g%iec ; nz = gv%ke
  diag_wt = dt * idt_diag

  if (cs%omega_frac >= 1.0) absf = 2.0*cs%omega
  do i=is,ie
    u_star = fluxes%ustar(i,j)
    if (associated(fluxes%ustar_shelf) .and. associated(fluxes%frac_shelf_h)) then
      if (fluxes%frac_shelf_h(i,j) > 0.0) &
        u_star = (1.0 - fluxes%frac_shelf_h(i,j)) * u_star + &
                  fluxes%frac_shelf_h(i,j) * fluxes%ustar_shelf(i,j)
    endif

    if (u_star < cs%ustar_min) u_star = cs%ustar_min
    if (cs%omega_frac < 1.0) then
      absf = 0.25*((abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i-1,j-1))) + &
                   (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i-1,j))))
      if (cs%omega_frac > 0.0) &
        absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
    endif
    absf_ustar = absf / u_star
    idecay_len_tke(i) = (absf_ustar * cs%TKE_decay) * gv%H_to_m

!    The first number in the denominator could be anywhere up to 16.  The
!  value of 3 was chosen to minimize the time-step dependence of the amount
!  of shear-driven mixing in 10 days of a 1-degree global model, emphasizing
!  the equatorial areas.  Although it is not cast as a parameter, it should
!  be considered an empirical parameter, and it might depend strongly on the
!  number of sublayers in the mixed layer and their locations.
!  The 0.41 is VonKarman's constant.  This equation assumes that small & large
!  scales contribute to mixed layer deepening at similar rates, even though
!  small scales are dissipated more rapidly (implying they are less efficient).
!     Ih = 1.0/(16.0*0.41*U_star*dt)
    ih = gv%H_to_m/(3.0*0.41*u_star*dt)
    cmke(1,i) = 4.0 * ih ; cmke(2,i) = (absf_ustar*gv%H_to_m) * ih

    if (idecay_len_tke(i) > 0.0) then
      exp_kh = exp(-htot(i)*idecay_len_tke(i))
    else
      exp_kh = 1.0
    endif

! 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 * (absf*htot(i))**3 / totEn)
    if (conv_en(i) < 0.0) conv_en(i) = 0.0
    if (ctke(i,1) > 0.0) then ; tke_ca = ctke(i,1) ; else ; tke_ca = 0.0 ; endif
    if ((htot(i) >= h_ca(i)) .or. (tke_ca == 0.0)) then
      toten = conv_en(i) + tke_ca

      if (toten > 0.0) then
        nstar_fc = cs%nstar * toten / (toten + 0.2 * &
                        sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_m))**3 * toten))
      else
        nstar_fc = cs%nstar
      endif
      nstar_ca = nstar_fc
    else
      ! This reconstructs the Buoyancy flux within the topmost htot of water.
      if (conv_en(i) > 0.0) then
        toten = conv_en(i) + tke_ca * (htot(i) / h_ca(i))
       nstar_fc = cs%nstar * toten / (toten + 0.2 * &
                        sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_m))**3 * toten))
      else
        nstar_fc = cs%nstar
      endif

      toten = conv_en(i) + tke_ca
      if (tke_ca > 0.0) then
        nstar_ca = cs%nstar * toten / (toten + 0.2 * &
                        sqrt(0.5 * dt * (absf*(h_ca(i)*gv%H_to_m))**3 * toten))
      else
        nstar_ca = cs%nstar
      endif
    endif

    if (dke_fc(i) + dke_ca(i,1) > 0.0) then
      if (htot(i) >= h_ca(i)) then
        mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
        mke_rate_ca = mke_rate_fc
      else
        mke_rate_fc = 1.0 / (1.0 + htot(i)*(cmke(1,i) + cmke(2,i)*htot(i)) )
        mke_rate_ca = 1.0 / (1.0 + h_ca(i)*(cmke(1,i) + cmke(2,i)*h_ca(i)) )
      endif
    else
      ! This branch just saves unnecessary calculations.
      mke_rate_fc = 1.0 ; mke_rate_ca = 1.0
    endif

    dke_conv = dke_ca(i,1) * mke_rate_ca + dke_fc(i) * mke_rate_fc
! At this point, it is assumed that cTKE is positive and stored in TKE_CA!
! Note: Removed factor of 2 in u*^3 terms.
    tke(i) = (dt*cs%mstar)*((u_star*u_star*u_star)*exp_kh) + &
             (exp_kh * dke_conv + nstar_fc*conv_en(i) + nstar_ca*tke_ca)
! Add additional TKE at river mouths

    if (cs%do_rivermix) then
      tke(i) = tke(i) + tke_river(i)*dt*exp_kh
    endif

    if (cs%TKE_diagnostics) then
      wind_tke_src = cs%mstar*(u_star*u_star*u_star) * diag_wt
      cs%diag_TKE_wind(i,j) = cs%diag_TKE_wind(i,j) + &
          wind_tke_src + tke_river(i) * diag_wt
      cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + dke_conv*idt_diag
      cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
          (exp_kh-1.0)*(wind_tke_src + dke_conv*idt_diag)
      cs%diag_TKE_conv(i,j) = cs%diag_TKE_conv(i,j) + &
          idt_diag*(nstar_fc*conv_en(i) + nstar_ca*tke_ca)
      cs%diag_TKE_conv_decay(i,j) = cs%diag_TKE_conv_decay(i,j) + &
          idt_diag*((cs%nstar-nstar_fc)*conv_en(i) + (cs%nstar-nstar_ca)*tke_ca)
      cs%diag_TKE_conv_s2(i,j) = cs%diag_TKE_conv_s2(i,j) + &
          idt_diag*(ctke(i,1)-tke_ca)
    endif
  enddo

end subroutine find_starting_tke

!> This subroutine calculates mechanically driven entrainment.
subroutine mechanical_entrainment(h, d_eb, htot, Ttot, Stot, uhtot, vhtot, &
                                  R0_tot, Rcv_tot, u, v, T, S, R0, Rcv, eps, &
                                  dR0_dT, dRcv_dT, cMKE, Idt_diag, nsw, &
                                  Pen_SW_bnd, opacity_band, TKE, &
                                  Idecay_len_TKE, j, ksort, G, GV, CS)
  type(ocean_grid_type),             intent(in)    :: G       !< The ocean's grid structure.
  type(verticalgrid_type),           intent(in)    :: GV      !< The ocean's vertical grid
                                                              !! structure.
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: h       !< 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),SZK_(GV)), intent(inout) :: d_eb    !< The downward increase across a
                                                              !! layer in the entrainment from
                                                              !! below, in H. Positive values go
                                                              !! with mass gain by a layer.
  real, dimension(SZI_(G)),          intent(inout) :: htot    !< The accumlated mixed layer
                                                              !! thickness, in H.
  real, dimension(SZI_(G)),          intent(inout) :: Ttot    !< The depth integrated mixed layer
                                                              !! temperature, in deg C H.
  real, dimension(SZI_(G)),          intent(inout) :: Stot    !< The depth integrated mixed layer
                                                              !! salinity, in psu H.
  real, dimension(SZI_(G)),          intent(inout) :: uhtot   !< The depth integrated mixed layer
                                                              !! zonal velocity, H m s-1.
  real, dimension(SZI_(G)),          intent(inout) :: vhtot   !< The integrated mixed layer
                                                              !! meridional velocity, H m s-1.
  real, dimension(SZI_(G)),          intent(inout) :: R0_tot  !< The integrated mixed layer
                                                              !! potential density referenced to 0
                                                              !! pressure, in H kg m-3.
  real, dimension(SZI_(G)),          intent(inout) :: Rcv_tot !< The integrated mixed layer
                                                              !! coordinate variable potential
                                                              !! density, in H kg m-3.
  real, dimension(SZI_(G),SZK_(GV)), intent(in)    :: u, v, T, S, R0, Rcv, eps
  real, dimension(SZI_(G)),          intent(in)    :: dR0_dT, dRcv_dT
  real, dimension(2,SZI_(G)),        intent(in)    :: cMKE
  real,                              intent(in)    :: Idt_diag
  integer,                           intent(in)    :: nsw     !< The number of bands of penetrating
                                                              !! shortwave radiation.
  real, dimension(:,:),              intent(inout) :: Pen_SW_bnd !< The penetrating shortwave
                                                              !! heating at the sea surface in each
                                                              !! penetrating band, in K H,
                                                              !! size nsw x SZI_(G).
  real, dimension(:,:,:),            intent(in)    :: opacity_band
  real, dimension(SZI_(G)),          intent(inout) :: TKE     !< The turbulent kinetic energy
                                                              !! available for mixing over a time
                                                              !! step, in m3 s-2.
  real, dimension(SZI_(G)),          intent(inout) :: Idecay_len_TKE
  integer,                           intent(in)    :: j       !< The j-index to work on.
  integer, dimension(SZI_(G),SZK_(GV)), &
                                     intent(in)    :: ksort !< The density-sorted k-indicies.
  type(bulkmixedlayer_cs),           pointer       :: CS    !< The control structure for this
                                                            !! module.

! This subroutine calculates mechanically driven entrainment.
! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)  The units
!                of h are referred to as H below.
!  (in/out)  d_eb - The downward increase across a layer in the entrainment from
!                   below, in H. Positive values go with mass gain by a layer.
!  (in/out)  htot - The accumlated mixed layer thickness, in H.
!  (in/out)  Ttot - The depth integrated mixed layer temperature, in deg C H.
!  (in/out)  Stot - The depth integrated mixed layer salinity, in psu H.
!  (in/out)  uhtot - The depth integrated mixed layer zonal velocity, H m s-1.
!  (in/out)  vhtot - The integrated mixed layer meridional velocity, H m s-1.
!  (in/out)  R0_tot - The integrated mixed layer potential density referenced
!                  to 0 pressure, in H kg m-3.
!  (in/out)  Rcv_tot - The integrated mixed layer coordinate variable
!                   potential density, in H kg m-3.
!  (in)      nsw - The number of bands of penetrating shortwave radiation.
!  (in/out)  Pen_SW_bnd - The penetrating shortwave heating at the sea surface
!                         in each penetrating band, in K H, size nsw x SZI_(G).
!  (in/out)  TKE - The turbulent kinetic energy available for mixing over a
!               time step, in m3 s-2.
!  (in)      j - The j-index to work on.
!  (in)      ksort - The density-sorted k-indicies.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure for this module.

  real :: SW_trans  !   The fraction of shortwave radiation that is not
                    ! absorbed in a layer, nondimensional.
  real :: Pen_absorbed  !   The amount of penetrative shortwave radiation
                        ! that is absorbed in a layer, in units of K m.
  real :: h_avail   ! The thickness in a layer available for entrainment in H.
  real :: h_ent     ! The thickness from a layer that is entrained, in H.
  real :: h_min, h_max ! Limits on the solution for h_ent, in H.
  real :: dh_Newt      !   The Newton's method estimate of the change in
                       ! h_ent between iterations, in H.
  real :: MKE_rate  !   The fraction of the energy in resolved shears
                    ! within the mixed layer that will be eliminated
                    ! within a timestep, nondim, 0 to 1.
  real :: HpE       !   The current thickness plus entrainment, H.
  real :: g_H_2Rho0   !   Half the gravitational acceleration times the
                      ! conversion from H to m divided by the mean density,
                      ! in m5 s-2 H-1 kg-1.
  real :: TKE_full_ent  ! The TKE remaining if a layer is fully entrained,
                        ! in units of m3 s-2.
  real :: dRL       ! Work required to mix water from the next layer
                    ! across the mixed layer, in m2 s-2.
  real :: Pen_En_Contrib  ! Penetrating SW contributions to the changes in
                          ! TKE, divided by layer thickness in m, in m2 s-2.
  real :: C1        ! A temporary variable in units of m2 s-2.
  real :: dMKE      ! A temporary variable related to the release of mean
                    ! kinetic energy, with units of H m3 s-2.
  real :: TKE_ent   ! The TKE that remains if h_ent were entrained, in m3 s-2.
  real :: TKE_ent1  ! The TKE that would remain, without considering the
                    ! release of mean kinetic energy, in m3 s2.
  real :: dTKE_dh   ! The partial derivative of TKE with h_ent, in m3 s-2 H-1.
  real :: Pen_dTKE_dh_Contrib ! The penetrating shortwave contribution to
                    ! dTKE_dh, in m2 s-2.
  real :: EF4_val   ! The result of EF4() (see later), 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 :: dEF4_dh   ! The partial derivative of EF4 with h, in H-2.
  real :: Pen_En1   ! A nondimensional temporary variable.
  real :: kh, exp_kh  ! Nondimensional temporary variables related to the.
  real :: f1_kh       ! fractional decay of TKE across a layer.
  real :: x1, e_x1      !   Nondimensional temporary variables related to
  real :: f1_x1, f2_x1  ! the relative decay of TKE and SW radiation across
  real :: f3_x1         ! a layer, and exponential-related functions of x1.
  real :: E_HxHpE   ! Entrainment divided by the product of the new and old
                    ! thicknesses, in H-1.
  real :: Hmix_min  ! The minimum mixed layer depth in H.
  real :: H_to_m    ! Local copies of unit conversion factors.
  real :: opacity
  real :: C1_3, C1_6, C1_24   !  1/3, 1/6, and 1/24.
  integer :: is, ie, nz, i, k, ks, itt, n

  c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0 ; c1_24 = 1.0/24.0
  h_to_m = gv%H_to_m
  g_h_2rho0 = (gv%g_Earth * h_to_m) / (2.0 * gv%Rho0)
  hmix_min = cs%Hmix_min * gv%m_to_H
  h_neglect = gv%H_subroundoff
  is = g%isc ; ie = g%iec ; nz = gv%ke

  do ks=1,nz

    do i=is,ie ; if (ksort(i,ks) > 0) then
      k = ksort(i,ks)

      h_avail = h(i,k) - eps(i,k)
      if ((h_avail > 0.) .and. ((tke(i) > 0.) .or. (htot(i) < hmix_min))) then
        drl = g_h_2rho0 * (r0(i,k)*htot(i) - r0_tot(i) )
        dmke = (h_to_m * cs%bulk_Ri_ML) * 0.5 * &
            ((uhtot(i)-u(i,k)*htot(i))**2 + (vhtot(i)-v(i,k)*htot(i))**2)

! Find the TKE that would remain if the entire layer were entrained.
        kh = idecay_len_tke(i)*h_avail ; exp_kh = exp(-kh)
        if (kh >= 2.0e-5) then ; f1_kh = (1.0-exp_kh) / kh
        else ; f1_kh = (1.0 - kh*(0.5 - c1_6*kh)) ; endif

        pen_en_contrib = 0.0
        do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
          opacity = opacity_band(n,i,k)
! Two different forms are used here to make sure that only negative
! values are taken into exponentials to avoid excessively large
! numbers.  They are, of course, mathematically identical.
          if (idecay_len_tke(i) > opacity) then
            x1 = (idecay_len_tke(i) - opacity) * h_avail
            if (x1 >= 2.0e-5) then
              e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
              f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
            else
              f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
              f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
            endif

            pen_en1 = exp(-opacity*h_avail) * &
               ((1.0+opacity*htot(i))*f1_x1 + opacity*h_avail*f3_x1)
          else
            x1 = (opacity - idecay_len_tke(i)) * h_avail
            if (x1 >= 2.0e-5) then
              e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
              f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
            else
              f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
              f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
            endif

            pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
                                 opacity*h_avail*f2_x1)
          endif
          pen_en_contrib = pen_en_contrib + &
            (g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)) * (pen_en1 - f1_kh)
        endif ; enddo

        hpe = htot(i)+h_avail
        mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
        ef4_val = ef4(htot(i)+h_neglect,h_avail,idecay_len_tke(i))
        tke_full_ent = (exp_kh*tke(i) - (h_avail*h_to_m)*(drl*f1_kh + pen_en_contrib)) + &
            mke_rate*dmke*ef4_val
        if ((tke_full_ent >= 0.0) .or. (h_avail+htot(i) <= hmix_min)) then
          ! The layer will be fully entrained.
          h_ent = h_avail

          if (cs%TKE_diagnostics) then
            e_hxhpe = h_ent / ((htot(i)+h_neglect)*(htot(i)+h_ent+h_neglect))
            cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
                idt_diag * ((exp_kh-1.0)*tke(i) + &
                            (h_ent*h_to_m)*drl*(1.0-f1_kh) + &
                            mke_rate*dmke*(ef4_val-e_hxhpe))
            cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
                idt_diag*(h_to_m*h_ent)*drl
            cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
                idt_diag*(h_to_m*h_ent)*pen_en_contrib
            cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
                idt_diag*mke_rate*dmke*e_hxhpe
          endif

          tke(i) = tke_full_ent
          if (tke(i) <= 0.0) tke(i) = 1e-300
        else
! The layer is only partially entrained.  The amount that will be
! entrained is determined iteratively.  No further layers will be
! entrained.
          h_min = 0.0 ; h_max = h_avail
          if (tke(i) <= 0.0) then
            h_ent = 0.0
          else
            h_ent = h_avail * tke(i) / (tke(i) - tke_full_ent)

            do itt=1,15
              ! Evaluate the TKE that would remain if h_ent were entrained.

              kh = idecay_len_tke(i)*h_ent ; exp_kh = exp(-kh)
              if (kh >= 2.0e-5) then
                f1_kh = (1.0-exp_kh) / kh
              else
                f1_kh = (1.0 - kh*(0.5 - c1_6*kh))
              endif


              pen_en_contrib = 0.0 ; pen_dtke_dh_contrib = 0.0
              do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
                ! Two different forms are used here to make sure that only negative
                ! values are taken into exponentials to avoid excessively large
                ! numbers.  They are, of course, mathematically identical.
                opacity = opacity_band(n,i,k)
                sw_trans = exp(-h_ent*opacity)
                if (idecay_len_tke(i) > opacity) then
                  x1 = (idecay_len_tke(i) - opacity) * h_ent
                  if (x1 >= 2.0e-5) then
                    e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
                    f3_x1 = ((e_x1-(1.0-x1))/(x1*x1))
                  else
                    f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
                    f3_x1 = (0.5 - x1*(c1_6 - c1_24*x1))
                  endif
                  pen_en1 = sw_trans * ((1.0+opacity*htot(i))*f1_x1 + &
                                          opacity*h_ent*f3_x1)
                else
                  x1 = (opacity - idecay_len_tke(i)) * h_ent
                  if (x1 >= 2.0e-5) then
                    e_x1 = exp(-x1) ; f1_x1 = ((1.0-e_x1)/(x1))
                    f2_x1 = ((1.0-(1.0+x1)*e_x1)/(x1*x1))
                  else
                    f1_x1 = (1.0 - x1*(0.5 - c1_6*x1))
                    f2_x1 = (0.5 - x1*(c1_3 - 0.125*x1))
                  endif

                  pen_en1 = exp_kh * ((1.0+opacity*htot(i))*f1_x1 + &
                                        opacity*h_ent*f2_x1)
                endif
                c1 = g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)
                pen_en_contrib = pen_en_contrib + c1*(pen_en1 - f1_kh)
                pen_dtke_dh_contrib = pen_dtke_dh_contrib + &
                           c1*((1.0-sw_trans) - opacity*(htot(i) + h_ent)*sw_trans)
              endif ; enddo ! (Pen_SW_bnd(n,i) > 0.0)

              tke_ent1 = exp_kh*tke(i) - (h_ent*h_to_m)*(drl*f1_kh + pen_en_contrib)
              ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i),def4_dh)
              hpe = htot(i)+h_ent
              mke_rate = 1.0/(1.0 + (cmke(1,i)*hpe + cmke(2,i)*hpe**2))
              tke_ent = tke_ent1 + dmke*ef4_val*mke_rate
              ! TKE_ent is the TKE that would remain if h_ent were entrained.

              dtke_dh = ((-idecay_len_tke(i)*tke_ent1 - drl*h_to_m) + &
                         pen_dtke_dh_contrib*h_to_m) + dmke * mke_rate* &
                       (def4_dh - ef4_val*mke_rate*(cmke(1,i)+2.0*cmke(2,i)*hpe))
              !  dh_Newt = -TKE_ent / dTKE_dh
              ! Bisect if the Newton's method prediction is outside of the bounded range.
              if (tke_ent > 0.0) then
                if ((h_max-h_ent)*(-dtke_dh) > tke_ent) then
                  dh_newt = -tke_ent / dtke_dh
                else
                  dh_newt = 0.5*(h_max-h_ent)
                endif
                h_min = h_ent
              else
                if ((h_min-h_ent)*(-dtke_dh) < tke_ent) then
                  dh_newt = -tke_ent / dtke_dh
                else
                  dh_newt = 0.5*(h_min-h_ent)
                endif
                h_max = h_ent
              endif
              h_ent = h_ent + dh_newt

              if (abs(dh_newt) < 0.2*gv%Angstrom) exit
            enddo
          endif

          if (h_ent < hmix_min-htot(i)) h_ent = hmix_min - htot(i)

          if (cs%TKE_diagnostics) then
            hpe = htot(i)+h_ent
            mke_rate = 1.0/(1.0 + cmke(1,i)*hpe + cmke(2,i)*hpe**2)
            ef4_val = ef4(htot(i)+h_neglect,h_ent,idecay_len_tke(i))

            e_hxhpe = h_ent / ((htot(i)+h_neglect)*(hpe+h_neglect))
            cs%diag_TKE_mech_decay(i,j) = cs%diag_TKE_mech_decay(i,j) + &
                idt_diag * ((exp_kh-1.0)*tke(i) + &
                            (h_ent*h_to_m)*drl*(1.0-f1_kh) + &
                             dmke*mke_rate*(ef4_val-e_hxhpe))
            cs%diag_TKE_mixing(i,j) = cs%diag_TKE_mixing(i,j) - &
                idt_diag*(h_ent*h_to_m)*drl
            cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
                idt_diag*(h_ent*h_to_m)*pen_en_contrib
            cs%diag_TKE_RiBulk(i,j) = cs%diag_TKE_RiBulk(i,j) + &
                idt_diag*dmke*mke_rate*e_hxhpe
          endif

          tke(i) = 0.0
        endif ! TKE_full_ent > 0.0

        pen_absorbed = 0.0
        do n=1,nsw ; if (pen_sw_bnd(n,i) > 0.0) then
          sw_trans = exp(-h_ent*opacity_band(n,i,k))
          pen_absorbed = pen_absorbed + pen_sw_bnd(n,i) * (1.0 - sw_trans)
          pen_sw_bnd(n,i) = pen_sw_bnd(n,i) * sw_trans
        endif ; enddo

        htot(i)   = htot(i)   + h_ent
        r0_tot(i) = r0_tot(i) + (h_ent * r0(i,k) + pen_absorbed*dr0_dt(i))
        h(i,k)    = h(i,k)    - h_ent
        d_eb(i,k) = d_eb(i,k) - h_ent

        stot(i)    = stot(i)    + h_ent * s(i,k)
        ttot(i)    = ttot(i)    + (h_ent * t(i,k) + pen_absorbed)
        rcv_tot(i) = rcv_tot(i) + (h_ent*rcv(i,k) + pen_absorbed*drcv_dt(i))

        uhtot(i) = uhtot(i) + u(i,k)*h_ent
        vhtot(i) = vhtot(i) + v(i,k)*h_ent
      endif ! h_avail > 0.0 .AND TKE(i) > 0.0

    endif ; enddo ! i loop
  enddo ! k loop

end subroutine mechanical_entrainment

!> This subroutine generates an array of indices that are sorted by layer
!! density.
subroutine sort_ml(h, R0, eps, G, GV, CS, ksort)
  type(ocean_grid_type),                intent(in)  :: G     !< The ocean's grid structure.
  type(verticalgrid_type),              intent(in)  :: GV    !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZK_(GV)),    intent(in)  :: h     !< 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),SZK_(GV)),    intent(in)  :: R0    !< The potential density used to sort
                                                             !! the layers, in kg m-3.
  real, dimension(SZI_(G),SZK_(GV)),    intent(in)  :: eps   !< The (small) thickness that must
                                                             !! remain in each layer, in H.
  type(bulkmixedlayer_cs),              pointer     :: CS    !< The control structure returned by a
                                                             !! previous call to mixedlayer_init.
  integer, dimension(SZI_(G),SZK_(GV)), intent(out) :: ksort !< The k-index to use in the sort.

!   This subroutine generates an array of indices that are sorted by layer
! density.

! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)  The units
!                of h are referred to as H below.
!  (in)      R0 - The potential density used to sort the layers, in kg m-3.
!  (in)      eps - The (small) thickness that must remain in each layer, in H.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      j - The meridional row to work on.
!  (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.
!  (out)     ksort - The k-index to use in the sort.
  real :: R0sort(szi_(g),szk_(gv))
  integer :: nsort(szi_(g))
  logical :: done_sorting(szi_(g))
  integer :: i, k, ks, is, ie, nz, nkmb

  is = g%isc ; ie = g%iec ; nz = gv%ke
  nkmb = cs%nkml+cs%nkbl

  !   Come up with a sorted index list of layers with increasing R0.
  ! Assume that the layers below nkmb are already stably stratified.
  ! Only layers that are thicker than eps are in the list.  Extra elements
  ! have an index of -1.

  !   This is coded using straight insertion, on the assumption that the
  ! layers are usually in the right order (or massless) anyway.

  do k=1,nz ; do i=is,ie ; ksort(i,k) = -1 ; enddo ; enddo

  do i=is,ie ; nsort(i) = 0 ; done_sorting(i) = .false. ; enddo
  do k=1,nz ; do i=is,ie ; if (h(i,k) > eps(i,k)) then
    if (done_sorting(i)) then ; ks = nsort(i) ; else
      do ks=nsort(i),1,-1
        if (r0(i,k) >= r0sort(i,ks)) exit
        r0sort(i,ks+1) = r0sort(i,ks) ; ksort(i,ks+1) = ksort(i,ks)
      enddo
      if ((k > nkmb) .and. (ks == nsort(i))) done_sorting(i) = .true.
    endif

    ksort(i,ks+1) = k
    r0sort(i,ks+1) = r0(i,k)
    nsort(i) = nsort(i) + 1
  endif ; enddo ; enddo

end subroutine sort_ml

!>   This subroutine actually moves properties between layers to achieve a
!! resorted state, with all of the resorted water either moved into the correct
!! interior layers or in the top nkmb layers.
subroutine resort_ml(h, T, S, R0, Rcv, RcvTgt, eps, d_ea, d_eb, ksort, G, GV, CS, &
                     dR0_dT, dR0_dS, dRcv_dT, dRcv_dS)
  type(ocean_grid_type),                intent(in)    :: G       !< The ocean's grid structure.
  type(verticalgrid_type),              intent(in)    :: GV      !< The ocean's vertical grid
                                                                 !! structure.
  real, dimension(SZI_(G),SZK0_(GV)),   intent(inout) :: h       !< Layer thickness, in m or kg m-2.
                                                                 !! (Intent in/out)  The units of h
                                                                 !! are referred to as H below.
                                                                 !! Layer 0 is the new mixed layer.
  real, dimension(SZI_(G),SZK0_(GV)),   intent(inout) :: T       !< Layer temperatures, in deg C.
  real, dimension(SZI_(G),SZK0_(GV)),   intent(inout) :: S       !< Layer salinities, in psu.
  real, dimension(SZI_(G),SZK0_(GV)),   intent(inout) :: R0      !< Potential density referenced to
                                                                 !! surface pressure, in kg m-3.
  real, dimension(SZI_(G),SZK0_(GV)),   intent(inout) :: Rcv     !< The coordinate defining
                                                                 !! potential density, in kg m-3.
  real, dimension(SZK_(GV)),            intent(in)    :: RcvTgt  !< The target value of Rcv for each
                                                                 !! layer, in kg m-3.
  real, dimension(SZI_(G),SZK_(GV)),    intent(inout) :: eps     !< The (small) thickness that must
                                                                 !! remain in each layer, in H.
  real, dimension(SZI_(G),SZK_(GV)),    intent(inout) :: d_ea    !< The upward increase across a
                                                                 !! layer in the entrainment from
                                                                 !! above, in m or kg m-2 (H).
                                                                 !!  Positive d_ea goes with layer
                                                                 !! thickness increases.
  real, dimension(SZI_(G),SZK_(GV)),    intent(inout) :: d_eb    !< The downward increase across a
                                                                 !! layer in the entrainment from
                                                                 !! below, in H. Positive values go
                                                                 !! with mass gain by a layer.
  integer, dimension(SZI_(G),SZK_(GV)), intent(in)    :: ksort   !< The density-sorted k-indicies.
  type(bulkmixedlayer_cs),              pointer       :: CS      !< The control structure for this
                                                                 !! module.
  real, dimension(SZI_(G)),             intent(in)    :: dR0_dT  !< The partial derivative of
                                                                 !! potential density referenced
                                                                 !! to the surface with potential
                                                                 !! temperature, in kg m-3 K-1.
  real, dimension(SZI_(G)),             intent(in)    :: dR0_dS  !< The partial derivative of
                                                                 !! cpotential density referenced
                                                                 !! to the surface with salinity,
                                                                 !! in kg m-3 psu-1.
  real, dimension(SZI_(G)),             intent(in)    :: dRcv_dT !< The partial derivative of
                                                                 !! coordinate defining potential
                                                                 !! density with potential
                                                                 !! temperature, in kg m-3 K-1.
  real, dimension(SZI_(G)),             intent(in)    :: dRcv_dS !< The partial derivative of
                                                                 !! coordinate defining potential
                                                                 !! density with salinity,
                                                                 !! in kg m-3 psu-1.

!   This subroutine actually moves properties between layers to achieve a
! resorted state, with all of the resorted water either moved into the correct
! interior layers or in the top nkmb layers.
!
! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)  The units of
!                h are referred to as H below.  Layer 0 is the new mixed layer.
!  (in/out)  T - Layer temperatures, in deg C.
!  (in/out)  S - Layer salinities, in psu.
!  (in/out)  R0 - Potential density referenced to surface pressure, in kg m-3.
!  (in/out)  Rcv - The coordinate defining potential density, in kg m-3.
!  (in)      RcvTgt - The target value of Rcv for each layer, in kg m-3.
!  (in)      eps - The (small) thickness that must remain in each layer, in H.
!  (in/out)  d_ea - The upward increase across a layer in the entrainment from
!                   above, in m or kg m-2 (H).  Positive d_ea goes with layer
!                   thickness increases.
!  (in/out)  d_eb - The downward increase across a layer in the entrainment from
!                   below, in H. Positive values go with mass gain by a layer.
!  (in)      ksort - The density-sorted k-indicies.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure for this module.
!  (in/out)  dR0_dT - The partial derivative of potential density referenced
!                     to the surface with potential temperature, in kg m-3 K-1.
!  (in/out)  dR0_dS - The partial derivative of cpotential density referenced
!                     to the surface with salinity, in kg m-3 psu-1.
!  (in/out)  dRcv_dT - The partial derivative of coordinate defining potential
!                      density with potential temperature, in kg m-3 K-1.
!  (in/out)  dRcv_dS - The partial derivative of coordinate defining potential
!                      density with salinity, in kg m-3 psu-1.

!   If there are no massive light layers above the deepest of the mixed- and
! buffer layers, do nothing (except perhaps to reshuffle these layers).
!   If there are nkbl or fewer layers above the deepest mixed- or buffer-
! layers, move them (in sorted order) into the buffer layers, even if they
! were previously interior layers.
!   If there are interior layers that are intermediate in density (both in-situ
! and the coordinate density (sigma-2)) between the newly forming mixed layer
! and a residual buffer- or mixed layer, and the number of massive layers above
! the deepest massive buffer or mixed layer is greater than nkbl, then split
! those buffer layers into peices that match the target density of the two
! nearest interior layers.
!   Otherwise, if there are more than nkbl+1 remaining massive layers
  real    :: h_move, h_tgt_old, I_hnew
  real    :: dT_dS_wt2, dT_dR, dS_dR, I_denom
  real    :: Rcv_int
  real    :: target_match_tol
  real    :: T_up, S_up, R0_up, I_hup, h_to_up
  real    :: T_dn, S_dn, R0_dn, I_hdn, h_to_dn
  real    :: wt_dn
  real    :: dR1, dR2
  real    :: dPE, hmin, min_dPE, min_hmin
  real, dimension(SZK_(GV)) :: &
    h_tmp, R0_tmp, T_tmp, S_tmp, Rcv_tmp
  integer :: ks_min
  logical :: sorted, leave_in_layer
  integer :: ks_deep(szi_(g)), k_count(szi_(g)), ks2_reverse(szi_(g), szk_(gv))
  integer :: ks2(szk_(gv))
  integer :: i, k, ks, is, ie, nz, k1, k2, k_tgt, k_src, k_int_top
  integer :: nks, nkmb, num_interior, top_interior_ks

  is = g%isc ; ie = g%iec ; nz = gv%ke
  nkmb = cs%nkml+cs%nkbl
  target_match_tol = 0.1 ! ### MAKE THIS A PARAMETER.

  dt_ds_wt2 = cs%dT_dS_wt**2

! Find out how many massive layers are above the deepest buffer or mixed layer.
  do i=is,ie ; ks_deep(i) = -1 ; k_count(i) = 0 ; enddo
  do ks=nz,1,-1 ; do i=is,ie ; if (ksort(i,ks) > 0) then
    k = ksort(i,ks)

    if (h(i,k) > eps(i,k)) then
      if (ks_deep(i) == -1) then
        if (k <= nkmb) then
          ks_deep(i) = ks ; k_count(i) = k_count(i) + 1
          ks2_reverse(i,k_count(i)) = k
        endif
      else
        k_count(i) = k_count(i) + 1
        ks2_reverse(i,k_count(i)) = k
      endif
    endif
  endif ; enddo ; enddo

  do i=is,ie ; if (k_count(i) > 1) then
    ! This column might need to be reshuffled.
    nks = k_count(i)

    ! Put ks2 in the right order and check whether reshuffling is needed.
    sorted = .true.
    ks2(nks) = ks2_reverse(i,1)
    do ks=nks-1,1,-1
      ks2(ks) = ks2_reverse(i,1+nks-ks)
      if (ks2(ks) > ks2(ks+1)) sorted = .false.
    enddo

    ! Go to the next column of no reshuffling is needed.
    if (sorted) cycle

    ! Find out how many interior layers are being reshuffled.  If none,
    ! then this is a simple swapping procedure.
    num_interior = 0 ; top_interior_ks = 0
    do ks=nks,1,-1 ; if (ks2(ks) > nkmb) then
      num_interior = num_interior+1 ; top_interior_ks = ks
    endif ; enddo

    if (num_interior >= 1) then
      ! Find the lightest interior layer with a target coordinate density
      ! greater than the newly forming mixed layer.
      do k=nkmb+1,nz ; if (rcv(i,0) < rcvtgt(k)) exit ; enddo
      k_int_top = k ; rcv_int = rcvtgt(k)

      i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
      dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
      ds_dr = drcv_ds(i) * i_denom


      ! Examine whether layers can be split out of existence.  Stop when there
      ! is a layer that cannot be handled this way, or when the topmost formerly
      ! interior layer has been dealt with.
      do ks = nks,top_interior_ks,-1
        k = ks2(ks)
        leave_in_layer = .false.
        if ((k > nkmb) .and. (rcv(i,k) <= rcvtgt(k))) then
          if (rcvtgt(k)-rcv(i,k) < target_match_tol*(rcvtgt(k) - rcvtgt(k-1))) &
            leave_in_layer = .true.
        elseif (k > nkmb) then
          if (rcv(i,k)-rcvtgt(k) < target_match_tol*(rcvtgt(k+1) - rcvtgt(k))) &
            leave_in_layer = .true.
        endif

        if (leave_in_layer) then
          ! Just drop this layer from the sorted list.
          nks = nks-1
        elseif (rcv(i,k) < rcv_int) then
          ! There is no interior layer with a target density that is intermediate
          ! between this layer and the mixed layer.
          exit
        else
          ! Try splitting the layer between two interior isopycnal layers.
          ! Find the target densities that bracket this layer.
          do k2=k_int_top+1,nz ; if (rcv(i,k) < rcvtgt(k2)) exit ; enddo
          if (k2>nz) exit

          ! This layer is bracketed in density between layers k2-1 and k2.

          dr1 = (rcvtgt(k2-1) - rcv(i,k)) ; dr2 = (rcvtgt(k2) - rcv(i,k))
          t_up = t(i,k) + dt_dr * dr1
          s_up = s(i,k) + ds_dr * dr1
          t_dn = t(i,k) + dt_dr * dr2
          s_dn = s(i,k) + ds_dr * dr2

          r0_up = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr1
          r0_dn = r0(i,k) + (dt_dr*dr0_dt(i) + ds_dr*dr0_ds(i)) * dr2

          ! Make sure the new properties are acceptable.
          if ((r0_up > r0(i,0)) .or. (r0_dn > r0(i,0))) &
            ! Avoid creating obviously unstable profiles.
            exit

          wt_dn = (rcv(i,k) - rcvtgt(k2-1)) / (rcvtgt(k2) - rcvtgt(k2-1))
          h_to_up = (h(i,k)-eps(i,k)) * (1.0 - wt_dn)
          h_to_dn = (h(i,k)-eps(i,k)) * wt_dn

          i_hup = 1.0 / (h(i,k2-1) + h_to_up)
          i_hdn = 1.0 / (h(i,k2) + h_to_dn)
          r0(i,k2-1) = (r0(i,k2)*h(i,k2-1) + r0_up*h_to_up) * i_hup
          r0(i,k2) = (r0(i,k2)*h(i,k2) + r0_dn*h_to_dn) * i_hdn

          t(i,k2-1) = (t(i,k2)*h(i,k2-1) + t_up*h_to_up) * i_hup
          t(i,k2) = (t(i,k2)*h(i,k2) + t_dn*h_to_dn) * i_hdn
          s(i,k2-1) = (s(i,k2)*h(i,k2-1) + s_up*h_to_up) * i_hup
          s(i,k2) = (s(i,k2)*h(i,k2) + s_dn*h_to_dn) * i_hdn
          rcv(i,k2-1) = (rcv(i,k2)*h(i,k2-1) + rcvtgt(k2-1)*h_to_up) * i_hup
          rcv(i,k2) = (rcv(i,k2)*h(i,k2) + rcvtgt(k2)*h_to_dn) * i_hdn

          h(i,k) = eps(i,k)
          h(i,k2) = h(i,k2) + h_to_dn
          h(i,k2-1) = h(i,k2-1) + h_to_up

          if (k > k2-1) then
            d_eb(i,k) = d_eb(i,k) - h_to_up
            d_eb(i,k2-1) = d_eb(i,k2-1) + h_to_up
          elseif (k < k2-1) then
            d_ea(i,k) = d_ea(i,k) - h_to_up
            d_ea(i,k2-1) = d_ea(i,k2-1) + h_to_up
          endif
          if (k > k2) then
            d_eb(i,k) = d_eb(i,k) - h_to_dn
            d_eb(i,k2) = d_eb(i,k2) + h_to_dn
          elseif (k < k2) then
            d_ea(i,k) = d_ea(i,k) - h_to_dn
            d_ea(i,k2) = d_ea(i,k2) + h_to_dn
          endif
          nks = nks-1
        endif
      enddo
    endif

    do while (nks > nkmb)
      ! Having already tried to move surface layers into the interior, there
      ! are still too many layers, and layers must be merged until nks=nkmb.
      ! Examine every merger of a layer with its neighbors, and merge the ones
      ! that increase the potential energy the least.  If there are layers
      ! with (apparently?) negative potential energy change, choose the one
      ! with the smallest total thickness.  Repeat until nkmb layers remain.
      ! Choose the smaller value for the remaining index for convenience.

      ks_min = -1 ; min_dpe = 1.0 ; min_hmin = 0.0
      do ks=1,nks-1
        k1 = ks2(ks) ; k2 = ks2(ks+1)
        dpe = max(0.0, (r0(i,k2)-r0(i,k1)) * h(i,k1) * h(i,k2))
        hmin = min(h(i,k1)-eps(i,k1), h(i,k2)-eps(i,k2))
        if ((ks_min < 0) .or. (dpe < min_dpe) .or. &
            ((dpe <= 0.0) .and. (hmin < min_hmin))) then
           ks_min = ks ; min_dpe = dpe ; min_hmin = hmin
        endif
      enddo

      ! Now merge the two layers that do the least damage.
      k1 = ks2(ks_min) ; k2 = ks2(ks_min+1)
      if (k1 < k2) then ; k_tgt = k1 ; k_src = k2
      else ; k_tgt = k2 ; k_src = k1 ; ks2(ks_min) = k2 ; endif

      h_tgt_old = h(i,k_tgt)
      h_move = h(i,k_src)-eps(i,k_src)
      h(i,k_src) = eps(i,k_src)
      h(i,k_tgt) = h(i,k_tgt) + h_move
      i_hnew = 1.0 / (h(i,k_tgt))
      r0(i,k_tgt) = (r0(i,k_tgt)*h_tgt_old + r0(i,k_src)*h_move) * i_hnew

      t(i,k_tgt) = (t(i,k_tgt)*h_tgt_old + t(i,k_src)*h_move) * i_hnew
      s(i,k_tgt) = (s(i,k_tgt)*h_tgt_old + s(i,k_src)*h_move) * i_hnew
      rcv(i,k_tgt) = (rcv(i,k_tgt)*h_tgt_old + rcv(i,k_src)*h_move) * i_hnew

      d_eb(i,k_src) = d_eb(i,k_src) - h_move
      d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move

      ! Remove the newly missing layer from the sorted list.
      do ks=ks_min+1,nks ; ks2(ks) = ks2(ks+1) ; enddo
      nks = nks-1
    enddo

    !   Check again whether the layers are sorted, and go on to the next column
    ! if they are.
    sorted = .true.
    do ks=1,nks-1 ; if (ks2(ks) > ks2(ks+1)) sorted = .false. ; enddo
    if (sorted) cycle

    if (nks > 1) then
      ! At this point, actually swap the properties of the layers, and place
      ! the remaining layers in order starting with nkmb.

      ! Save all the properties of the nkmb layers that might be replaced.
      do k=1,nkmb
        h_tmp(k) = h(i,k) ; r0_tmp(k) = r0(i,k)
        t_tmp(k) = t(i,k) ; s_tmp(k) = s(i,k) ; rcv_tmp(k) = rcv(i,k)

        h(i,k) = 0.0
      enddo

      do ks=nks,1,-1
        k_tgt = nkmb - nks + ks ; k_src = ks2(ks)
        if (k_tgt == k_src) then
          h(i,k_tgt) = h_tmp(k_tgt)  ! This layer doesn't move, so put the water back.
          cycle
        endif

        ! Note below that eps=0 for k<=nkmb.
        if (k_src > nkmb) then
          h_move = h(i,k_src)-eps(i,k_src)
          h(i,k_src) = eps(i,k_src)
          h(i,k_tgt) = h_move
          r0(i,k_tgt) = r0(i,k_src)

          t(i,k_tgt) = t(i,k_src) ; s(i,k_tgt) = s(i,k_src)
          rcv(i,k_tgt) = rcv(i,k_src)

          d_eb(i,k_src) = d_eb(i,k_src) - h_move
          d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_move
        else
          h(i,k_tgt) = h_tmp(k_src)
          r0(i,k_tgt) = r0_tmp(k_src)

          t(i,k_tgt) = t_tmp(k_src) ; s(i,k_tgt) = s_tmp(k_src)
          rcv(i,k_tgt) = rcv_tmp(k_src)

          if (k_src > k_tgt) then
            d_eb(i,k_src) = d_eb(i,k_src) - h_tmp(k_src)
            d_eb(i,k_tgt) = d_eb(i,k_tgt) + h_tmp(k_src)
          else
            d_ea(i,k_src) = d_ea(i,k_src) - h_tmp(k_src)
            d_ea(i,k_tgt) = d_ea(i,k_tgt) + h_tmp(k_src)
          endif
        endif
      enddo
    endif

  endif ; enddo

end subroutine resort_ml

!> This subroutine moves any water left in the former mixed layers into the
!! two buffer layers and may also move buffer layer water into the interior
!! isopycnal layers.
subroutine mixedlayer_detrain_2(h, T, S, R0, Rcv, RcvTgt, dt, dt_diag, d_ea, j, G, GV, CS, &
                                dR0_dT, dR0_dS, dRcv_dT, dRcv_dS, max_BL_det)
  type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure.
  type(verticalgrid_type),            intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h    !< Layer thickness, in m or kg m-2.
                                                            !! (Intent in/out)  The units of h are
                                                            !! referred to as H below.
                                                            !!  Layer 0 is the new mixed layer.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T    !< Potential temperature, in C.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S    !< Salinity, in psu.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0   !< Potential density referenced to
                                                            !! surface pressure, in kg m-3.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv  !< The coordinate defining potential
                                                            !! density, in kg m-3.
  real, dimension(SZK_(GV)),          intent(in)    :: RcvTgt  !< The target value of Rcv for each
                                                            !! layer, in kg m-3.
  real,                               intent(in)    :: dt   !< Time increment, in s.
  real,                               intent(in)    :: dt_diag !< The diagnostic time step, in s.
  real, dimension(SZI_(G),SZK_(GV)),  intent(inout) :: d_ea !< The upward increase across a layer in
                                                            !! the entrainment from above, in m or
                                                            !! kg m-2 (H).  Positive d_ea goes with
                                                            !! layer thickness increases.
  integer,                            intent(in)    :: j    !< The meridional row to work on.
  type(bulkmixedlayer_cs),            pointer       :: CS   !< The control structure returned by a
                                                            !! previous call to mixedlayer_init.
  real, dimension(SZI_(G)),           intent(in)    :: dR0_dT  !< The partial derivative of
                                                            !! potential density referenced to the
                                                            !! surface with potential temperature,
                                                            !! in kg m-3 K-1.
  real, dimension(SZI_(G)),           intent(in)    :: dR0_dS  !< The partial derivative of
                                                            !! cpotential density referenced to the
                                                            !! surface with salinity,
                                                            !! in kg m-3 psu-1.
  real, dimension(SZI_(G)),           intent(in)    :: dRcv_dT !< The partial derivative of
                                                            !! coordinate defining potential density
                                                            !! with potential temperature,
                                                            !! in kg m-3 K-1.
  real, dimension(SZI_(G)),           intent(in)    :: dRcv_dS !< The partial derivative of
                                                            !! coordinate defining potential density
                                                            !! with salinity, in kg m-3 psu-1.
  real, dimension(SZI_(G)),           intent(in)    :: max_BL_det !< If non-negative, the maximum
                                                            !! detrainment permitted from the buffer
                                                            !! layers, in H.

! This subroutine moves any water left in the former mixed layers into the
! two buffer layers and may also move buffer layer water into the interior
! isopycnal layers.

! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)  The units of
!                h are referred to as H below.  Layer 0 is the new mixed layer.
!  (in/out)  T - Potential temperature, in C.
!  (in/out)  S - Salinity, in psu.
!  (in/out)  R0 - Potential density referenced to surface pressure, in kg m-3.
!  (in/out)  Rcv - The coordinate defining potential density, in kg m-3.
!  (in)      RcvTgt - The target value of Rcv for each layer, in kg m-3.
!  (in)      dt - Time increment, in s.
!  (in)      dt_diag - The diagnostic time step, in s.
!  (in/out)  d_ea - The upward increase across a layer in the entrainment from
!                   above, in m or kg m-2 (H).  Positive d_ea goes with layer
!                   thickness increases.
!  (in)      j - The meridional row to work on.
!  (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)      max_BL_det - If non-negative, the maximum detrainment permitted
!                         from the buffer layers, in H.
!  (in)  dR0_dT - The partial derivative of potential density referenced
!                 to the surface with potential temperature, in kg m-3 K-1.
!  (in)  dR0_dS - The partial derivative of cpotential density referenced
!                 to the surface with salinity, in kg m-3 psu-1.
!  (in)  dRcv_dT - The partial derivative of coordinate defining potential
!                  density with potential temperature, in kg m-3 K-1.
!  (in)  dRcv_dS - The partial derivative of coordinate defining potential
!                  density with salinity, in kg m-3 psu-1.
  real :: h_to_bl                 ! The total thickness detrained to the buffer
                                  ! layers, in H (the units of h).
  real :: R0_to_bl, Rcv_to_bl     ! The depth integrated amount of R0, Rcv, T
  real :: T_to_bl, S_to_bl        ! and S that is detrained to the buffer layer,
                                  ! in H kg m-3, H kg m-3, K H, and psu H.

  real :: h_min_bl                ! The minimum buffer layer thickness, in H.
  real :: h_min_bl_thick          ! The minimum buffer layer thickness when the
                                  ! mixed layer is very large, in H.
  real :: h_min_bl_frac_ml = 0.05 ! The minimum buffer layer thickness relative
                                  ! to the total mixed layer thickness for thin
                                  ! mixed layers, nondim., maybe 0.1/CS%nkbl.

  real :: h1, h2                  ! Scalar variables holding the values of
                                  ! h(i,CS%nkml+1) and h(i,CS%nkml+2), in H.
  real :: h1_avail                ! The thickess of the upper buffer layer
                                  ! available to move into the lower buffer
                                  ! layer, in H.
  real :: stays                   ! stays is the thickness of the upper buffer
                                  ! layer that remains there, in units of H.
  real :: stays_min, stays_max    ! The minimum and maximum permitted values of
                                  ! stays, in units of H.

  logical :: mergeable_bl         ! If true, it is an option to combine the two
                                  ! buffer layers and create water that matches
                                  ! the target density of an interior layer.
  real :: stays_merge             ! If the two buffer layers can be combined
                                  ! stays_merge is the thickness of the upper
                                  ! layer that remains, in units of H.
  real :: stays_min_merge         ! The minimum allowed value of stays_merge in H.

  real :: dR0_2dz, dRcv_2dz       ! Half the vertical gradients of R0, Rcv, T, and
!  real :: dT_2dz, dS_2dz          ! S, in kg m-4, kg m-4, K m-1, and psu m-1.
  real :: scale_slope             ! A nondimensional number < 1 used to scale down
                                  ! the slope within the upper buffer layer when
                                  ! water MUST be detrained to the lower layer.

  real :: dPE_extrap              ! The potential energy change due to dispersive
                                  ! advection or mixing layers, divided by
                                  ! rho_0*g, in units of H2.
  real :: dPE_det, dPE_merge      ! The energy required to mix the detrained water
                                  ! into the buffer layer or the merge the two
                                  ! buffer layers, both in units of J H2 m-4.

  real :: h_from_ml               ! The amount of additional water that must be
                                  ! drawn from the mixed layer, in H.
  real :: h_det_h2                ! The amount of detrained water and mixed layer
                                  ! water that will go directly into the lower
                                  ! buffer layer, in H.
  real :: h_det_to_h2, h_ml_to_h2 ! All of the variables hA_to_hB are the
  real :: h_det_to_h1, h_ml_to_h1 ! thickess fluxes from one layer to another,
  real :: h1_to_h2, h1_to_k0      ! in H, with h_det the detrained water, h_ml
  real :: h2_to_k1, h2_to_k1_rem  ! the actively mixed layer, h1 and h2 the upper
                                  ! and lower buffer layers, and k0 and k1 the
                                  ! interior layers that are just lighter and
                                  ! just denser than the lower buffer layer.

  real :: R0_det, T_det, S_det    ! Detrained values of R0, T, and S.
  real :: Rcv_stays, R0_stays     ! Values of Rcv and R0 that stay in a layer.
  real :: T_stays, S_stays        ! Values of T and S that stay in a layer.
  real :: dSpice_det, dSpice_stays! The spiciness difference between an original
                                  ! buffer layer and the water that moves into
                                  ! an interior layer or that stays in that
                                  ! layer, in kg m-3.
  real :: dSpice_lim, dSpice_lim2 ! Limits to the spiciness difference between
                                  ! the lower buffer layer and the water that
                                  ! moves into an interior layer, in kg m-3.
  real :: dSpice_2dz              ! The vertical gradient of spiciness used for
                                  ! advection, in kg m-4.

  real :: dPE_ratio               ! Multiplier of dPE_det at which merging is
                                  ! permitted - here (detrainment_per_day/dt)*30
                                  ! days?
  real :: num_events              ! The number of detrainment events over which
                                  ! to prefer merging the buffer layers.
  real :: detrainment_timescale   ! The typical timescale for a detrainment
                                  ! event, in s.
  real :: dPE_time_ratio          ! Larger of 1 and the detrainment_timescale
                                  ! over dt, nondimensional.
  real :: dT_dS_gauge, dS_dT_gauge ! The relative scales of temperature and
                                  ! salinity changes in defining spiciness, in
                                  ! K psu-1 and psu K-1.
  real :: I_denom                 ! A work variable with units of psu2 m6 kg-2.

  real :: G_2                     ! 1/2 G_Earth, in m s-2.
  real :: Rho0xG                  ! Rho0 times G_Earth, in kg m-2 s-2.
  real :: I2Rho0                  ! 1 / (2 Rho0), in m3 kg-1.
  real :: Idt_H2                  ! The square of the conversion from thickness
                                  ! to m divided by the time step in m2 H-2 s-1.
  logical :: stable_Rcv           ! If true, the buffer layers are stable with
                                  ! respect to the coordinate potential density.
  real :: h_neglect ! A thickness that is so small it is usually lost
                    ! in roundoff and can be neglected, in H.

  real :: s1en                    ! A work variable with units of H2 kg m s-3.
  real :: s1, s2, bh0             ! Work variables with units of H.
  real :: s3sq                    ! A work variable with units of H2.
  real :: I_ya, b1                ! Nondimensional work variables.
  real :: Ih, Ihdet, Ih1f, Ih2f   ! Assorted inverse thickness work variables,
  real :: Ihk0, Ihk1, Ih12        ! all with units of H-1.
  real :: dR1, dR2, dR2b, dRk1    ! Assorted density difference work variables,
  real :: dR0, dR21, dRcv         ! all with units of kg m-3.
  real :: dRcv_stays, dRcv_det, dRcv_lim
  real :: Angstrom                ! The minumum layer thickness, in H.

  real :: h2_to_k1_lim, T_new, S_new, T_max, T_min, S_max, S_min
  character(len=200) :: mesg

  integer :: i, k, k0, k1, is, ie, nz, kb1, kb2, nkmb
  is = g%isc ; ie = g%iec ; nz = gv%ke
  kb1 = cs%nkml+1; kb2 = cs%nkml+2
  nkmb = cs%nkml+cs%nkbl
  h_neglect = gv%H_subroundoff
  g_2 = 0.5*gv%g_Earth
  rho0xg = gv%Rho0 * gv%g_Earth
  idt_h2 = gv%H_to_m**2 / dt_diag
  i2rho0 = 0.5 / gv%Rho0
  angstrom = gv%Angstrom

  ! This is hard coding of arbitrary and dimensional numbers.
  h_min_bl_thick = 5.0 * gv%m_to_H
  dt_ds_gauge = cs%dT_dS_wt ; ds_dt_gauge = 1.0 /dt_ds_gauge
  num_events = 10.0
  detrainment_timescale = 4.0*3600.0

  if (cs%nkbl /= 2) call mom_error(fatal, "MOM_mixed_layer"// &
                        "CS%nkbl must be 2 in mixedlayer_detrain_2.")

  if (dt < detrainment_timescale) then ; dpe_time_ratio = detrainment_timescale/dt
  else ; dpe_time_ratio = 1.0 ; endif

  do i=is,ie

  ! Determine all of the properties being detrained from the mixed layer.

  ! As coded this has the k and i loop orders switched, but k is CS%nkml is
  ! often just 1 or 2, so this seems like it should not be a problem, especially
  ! since it means that a number of variables can now be scalars, not arrays.
    h_to_bl = 0.0 ; r0_to_bl = 0.0
    rcv_to_bl = 0.0 ; t_to_bl = 0.0 ; s_to_bl = 0.0

    do k=1,cs%nkml ; if (h(i,k) > 0.0) then
      h_to_bl = h_to_bl + h(i,k)
      r0_to_bl = r0_to_bl + r0(i,k)*h(i,k)

      rcv_to_bl = rcv_to_bl + rcv(i,k)*h(i,k)
      t_to_bl = t_to_bl + t(i,k)*h(i,k)
      s_to_bl = s_to_bl + s(i,k)*h(i,k)

      d_ea(i,k) = d_ea(i,k) - h(i,k)
      h(i,k) = 0.0
    endif ; enddo
    if (h_to_bl > 0.0) then ; r0_det = r0_to_bl / h_to_bl
    else ; r0_det = r0(i,0) ; endif

    ! This code does both downward detrainment from both the mixed layer and the
    ! buffer layers.
    !   Several considerations apply in detraining water into the interior:
    ! (1) Water only moves into the interior from the deeper buffer layer,
    !     so the deeper buffer layer must have some mass.
    ! (2) The upper buffer layer must have some mass so the extrapolation of
    !     density is meaningful (i.e. there is not detrainment from the buffer
    !     layers when there is strong mixed layer entrainment).
    ! (3) The lower buffer layer density extrapolated to its base with a
    !     linear fit between the two layers must exceed the density of the
    !     next denser interior layer.
    ! (4) The average extroplated coordinate density that is moved into the
    !     isopycnal interior matches the target value for that layer.
    ! (5) The potential energy change is calculated and might be used later
    !     to allow the upper buffer layer to mix more into the lower buffer
    !     layer.

    ! Determine whether more must be detrained from the mixed layer to keep a
    ! minimal amount of mass in the buffer layers.  In this case the 5% of the
    ! mixed layer thickness is hard-coded, but probably shouldn't be!
    h_min_bl = min(h_min_bl_thick,h_min_bl_frac_ml*h(i,0))

    stable_rcv = .true.
    if (((r0(i,kb2)-r0(i,kb1)) * (rcv(i,kb2)-rcv(i,kb1)) <= 0.0)) &
      stable_rcv = .false.

    h1 = h(i,kb1) ; h2 = h(i,kb2)

    h2_to_k1_rem = (h1 + h2) + h_to_bl
    if ((max_bl_det(i) >= 0.0) .and. (h2_to_k1_rem > max_bl_det(i))) &
      h2_to_k1_rem = max_bl_det(i)


    if ((h2 == 0.0) .and. (h1 > 0.0)) then
      ! The lower buffer layer has been eliminated either by convective
      ! adjustment or entrainment from the interior, and its current properties
      ! are not meaningful, but may later be used to determine the properties of
      ! waters moving into the lower buffer layer.  So the properties of the
      ! lower buffer layer are set to be between those of the upper buffer layer
      ! and the next denser interior layer, measured by R0.  This probably does
      ! not happen very often, so I am not too worried about the inefficiency of
      ! the following loop.
      do k1=kb2+1,nz ; if (h(i,k1) > 2.0*angstrom) exit ; enddo

      r0(i,kb2) = r0(i,kb1)

      rcv(i,kb2)=rcv(i,kb1) ; t(i,kb2)=t(i,kb1) ; s(i,kb2)=s(i,kb1)


      if (k1 <= nz) then ; if (r0(i,k1) >= r0(i,kb1)) then
        r0(i,kb2) = 0.5*(r0(i,kb1)+r0(i,k1))

        rcv(i,kb2) = 0.5*(rcv(i,kb1)+rcv(i,k1))
        t(i,kb2) = 0.5*(t(i,kb1)+t(i,k1))
        s(i,kb2) = 0.5*(s(i,kb1)+s(i,k1))
      endif ; endif
    endif ! (h2 = 0 && h1 > 0)

    dpe_extrap = 0.0 ; dpe_merge = 0.0
    mergeable_bl = .false.
    if ((h1 > 0.0) .and. (h2 > 0.0) .and. (h_to_bl > 0.0) .and. &
        (stable_rcv)) then
      ! Check whether it is permissible for the buffer layers to detrain
      ! into the interior isopycnal layers.

      ! Determine the layer that has the lightest target density that is
      ! denser than the lowermost buffer layer.
      do k1=kb2+1,nz ; if (rcvtgt(k1) >= rcv(i,kb2)) exit ; enddo ; k0 = k1-1
      dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)

      ! Use an energy-balanced combination of downwind advection into the next
      ! denser interior layer and upwind advection from the upper buffer layer
      ! into the lower one, each with an energy change that equals that required
      ! to mix the detrained water with the upper buffer layer.
      h1_avail = h1 - max(0.0,h_min_bl-h_to_bl)
      if ((k1<=nz) .and. (h2 > h_min_bl) .and. (h1_avail > 0.0) .and. &
          (r0(i,kb1) < r0(i,kb2)) .and. (h_to_bl*r0(i,kb1) > r0_to_bl)) then
        drk1 = (rcvtgt(k1) - rcv(i,kb2)) * (r0(i,kb2) - r0(i,kb1)) / &
                                           (rcv(i,kb2) - rcv(i,kb1))
        b1 = drk1 / (r0(i,kb2) - r0(i,kb1))
        ! b1 = RcvTgt(k1) - Rcv(i,kb2)) / (Rcv(i,kb2) - Rcv(i,kb1))

        ! Apply several limits to the detrainment.
        ! Entrain less than the mass in h2, and keep the base of the buffer
        ! layers from becoming shallower than any neighbors.
        h2_to_k1 = min(h2 - h_min_bl, h2_to_k1_rem)
        ! Balance downwind advection of density into the layer below the
        ! buffer layers with upwind advection from the layer above.
        if (h2_to_k1*(h1_avail + b1*(h1_avail + h2)) > h2*h1_avail) &
          h2_to_k1 = (h2*h1_avail) / (h1_avail + b1*(h1_avail + h2))
        if (h2_to_k1*(drk1 * h2) > (h_to_bl*r0(i,kb1) - r0_to_bl) * h1) &
          h2_to_k1 = (h_to_bl*r0(i,kb1) - r0_to_bl) * h1 / (drk1 * h2)

        if ((k1==kb2+1) .and. (cs%BL_extrap_lim > 0.)) then
          ! Simply do not detrain very light water into the lightest isopycnal
          ! coordinate layers if the density jump is too large.
          drcv_lim = rcv(i,kb2)-rcv(i,0)
          do k=1,kb2 ; drcv_lim = max(drcv_lim, rcv(i,kb2)-rcv(i,k)) ; enddo
          drcv_lim = cs%BL_extrap_lim*drcv_lim
          if ((rcvtgt(k1) - rcv(i,kb2)) >= drcv_lim) then
            h2_to_k1 = 0.0
          elseif ((rcvtgt(k1) - rcv(i,kb2)) > 0.5*drcv_lim) then
            h2_to_k1 = h2_to_k1 * (2.0 - 2.0*((rcvtgt(k1) - rcv(i,kb2)) / drcv_lim))
          endif
        endif

        drcv = (rcvtgt(k1) - rcv(i,kb2))

        ! Use 2nd order upwind advection of spiciness, limited by the values
        ! in deeper thick layers to determine the detrained temperature and
        ! salinity.
        dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
                      dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
                      (h2 - h2_to_k1) / (h1 + h2)
        dspice_lim = 0.0
        if (h(i,k1) > 10.0*angstrom) then
          dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
                       dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
          if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
        endif
        if (k1<nz) then ; if (h(i,k1+1) > 10.0*angstrom) then
          dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
                        dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
          if ((dspice_det*dspice_lim2 > 0.0) .and. &
              (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
        endif ; endif
        if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim

        i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
        t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
            (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
        s_det = s(i,kb2) + i_denom * &
            (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
        ! The detrained values of R0 are based on changes in T and S.
        r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
                             (s_det-s(i,kb2)) * dr0_ds(i)

        if (cs%BL_extrap_lim >= 0.) then
          ! Only do this detrainment if the new layer's temperature and salinity
          ! are not too far outside of the range of previous values.
          if (h(i,k1) > 10.0*angstrom) then
            t_min = min(t(i,kb1), t(i,kb2), t(i,k1)) - cs%Allowed_T_chg
            t_max = max(t(i,kb1), t(i,kb2), t(i,k1)) + cs%Allowed_T_chg
            s_min = min(s(i,kb1), s(i,kb2), s(i,k1)) - cs%Allowed_S_chg
            s_max = max(s(i,kb1), s(i,kb2), s(i,k1)) + cs%Allowed_S_chg
          else
            t_min = min(t(i,kb1), t(i,kb2)) - cs%Allowed_T_chg
            t_max = max(t(i,kb1), t(i,kb2)) + cs%Allowed_T_chg
            s_min = min(s(i,kb1), s(i,kb2)) - cs%Allowed_S_chg
            s_max = max(s(i,kb1), s(i,kb2)) + cs%Allowed_S_chg
          endif
          ihk1 = 1.0 / (h(i,k1) + h2_to_k1)
          t_new = (h(i,k1)*t(i,k1) + h2_to_k1*t_det) * ihk1
          s_new = (h(i,k1)*s(i,k1) + h2_to_k1*s_det) * ihk1
          ! A less restrictive limit might be used here.
          if ((t_new < t_min) .or. (t_new > t_max) .or. &
              (s_new < s_min) .or. (s_new > s_max)) &
            h2_to_k1 = 0.0
        endif

        h1_to_h2 = b1*h2*h2_to_k1 / (h2 - (1.0+b1)*h2_to_k1)

        ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
        ih2f = 1.0 / ((h(i,kb2) - h2_to_k1) + h1_to_h2)

        rcv(i,kb2) = ((h(i,kb2)*rcv(i,kb2) - h2_to_k1*rcvtgt(k1)) + &
                      h1_to_h2*rcv(i,kb1))*ih2f
        rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1

        t(i,kb2) = ((h(i,kb2)*t(i,kb2) - h2_to_k1*t_det) + &
                    h1_to_h2*t(i,kb1)) * ih2f
        t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1

        s(i,kb2) = ((h(i,kb2)*s(i,kb2) - h2_to_k1*s_det) + &
                    h1_to_h2*s(i,kb1)) * ih2f
        s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1

        ! Changes in R0 are based on changes in T and S.
        r0(i,kb2) = ((h(i,kb2)*r0(i,kb2) - h2_to_k1*r0_det) + &
                     h1_to_h2*r0(i,kb1)) * ih2f
        r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1

        h(i,kb1) = h(i,kb1) - h1_to_h2 ; h1 = h(i,kb1)
        h(i,kb2) = (h(i,kb2) - h2_to_k1) + h1_to_h2 ; h2 = h(i,kb2)
        h(i,k1) = h(i,k1) + h2_to_k1

        d_ea(i,kb1) = d_ea(i,kb1) - h1_to_h2
        d_ea(i,kb2) = (d_ea(i,kb2) - h2_to_k1) + h1_to_h2
        d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
        h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)

        !   The lower buffer layer has become lighter - it may be necessary to
        ! adjust k1 lighter.
        if ((k1>kb2+1) .and. (rcvtgt(k1-1) >= rcv(i,kb2))) then
          do k1=k1,kb2+1,-1 ; if (rcvtgt(k1-1) < rcv(i,kb2)) exit ; enddo
        endif
      endif

      k0 = k1-1
      dr1 = rcvtgt(k0)-rcv(i,kb1) ; dr2 = rcv(i,kb2)-rcvtgt(k0)

      if ((k0>kb2) .and. (dr1 > 0.0) .and. (h1 > h_min_bl) .and. &
          (h2*dr2 < h1*dr1) .and. (r0(i,kb2) > r0(i,kb1))) then
        ! An interior isopycnal layer (k0) is intermediate in density between
        ! the two buffer layers, and there can be detrainment. The entire
        ! lower buffer layer is combined with a portion of the upper buffer
        ! layer to match the target density of layer k0.
        stays_merge = 2.0*(h1+h2)*(h1*dr1 - h2*dr2) / &
                     ((dr1+dr2)*h1 + dr1*(h1+h2) + &
                      sqrt((dr2*h1-dr1*h2)**2 + 4*(h1+h2)*h2*(dr1+dr2)*dr2))

        stays_min_merge = max(h_min_bl, 2.0*h_min_bl - h_to_bl, &
                  h1 - (h1+h2)*(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1)))
        if ((stays_merge > stays_min_merge) .and. &
            (stays_merge + h2_to_k1_rem >= h1 + h2)) then
          mergeable_bl = .true.
          dpe_merge = g_2*(r0(i,kb2)-r0(i,kb1))*(h1-stays_merge)*(h2-stays_merge)
        endif
      endif

      if ((k1<=nz).and.(.not.mergeable_bl)) then
        ! Check whether linear extrapolation of density (i.e. 2nd order upwind
        ! advection) will allow some of the lower buffer layer to detrain into
        ! the next denser interior layer (k1).
        dr2b = rcvtgt(k1)-rcv(i,kb2) ; dr21 = rcv(i,kb2) - rcv(i,kb1)
        if (dr2b*(h1+h2) < h2*dr21) then
          ! Some of layer kb2 is denser than k1.
          h2_to_k1 = min(h2 - (h1+h2) * dr2b / dr21, h2_to_k1_rem)

          if (h2 > h2_to_k1) then
            drcv = (rcvtgt(k1) - rcv(i,kb2))

            ! Use 2nd order upwind advection of spiciness, limited by the values
            ! in deeper thick layers to determine the detrained temperature and
            ! salinity.
            dspice_det = (ds_dt_gauge*drcv_ds(i)*(t(i,kb2)-t(i,kb1)) - &
                          dt_ds_gauge*drcv_dt(i)*(s(i,kb2)-s(i,kb1))) * &
                          (h2 - h2_to_k1) / (h1 + h2)
            dspice_lim = 0.0
            if (h(i,k1) > 10.0*angstrom) then
              dspice_lim = ds_dt_gauge*drcv_ds(i)*(t(i,k1)-t(i,kb2)) - &
                           dt_ds_gauge*drcv_dt(i)*(s(i,k1)-s(i,kb2))
              if (dspice_det*dspice_lim <= 0.0) dspice_lim = 0.0
            endif
            if (k1<nz) then; if (h(i,k1+1) > 10.0*angstrom) then
              dspice_lim2 = ds_dt_gauge*drcv_ds(i)*(t(i,k1+1)-t(i,kb2)) - &
                            dt_ds_gauge*drcv_dt(i)*(s(i,k1+1)-s(i,kb2))
              if ((dspice_det*dspice_lim2 > 0.0) .and. &
                  (abs(dspice_lim2) > abs(dspice_lim))) dspice_lim = dspice_lim2
            endif; endif
            if (abs(dspice_det) > abs(dspice_lim)) dspice_det = dspice_lim

            i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
            t_det = t(i,kb2) + dt_ds_gauge * i_denom * &
                (dt_ds_gauge * drcv_dt(i) * drcv + drcv_ds(i) * dspice_det)
            s_det = s(i,kb2) + i_denom * &
                (drcv_ds(i) * drcv - dt_ds_gauge * drcv_dt(i) * dspice_det)
            ! The detrained values of R0 are based on changes in T and S.
            r0_det = r0(i,kb2) + (t_det-t(i,kb2)) * dr0_dt(i) + &
                                 (s_det-s(i,kb2)) * dr0_ds(i)

            ! Now that the properties of the detrained water are known,
            ! potentially limit the amount of water that is detrained to
            ! avoid creating unphysical properties in the remaining water.
            ih2f = 1.0 / (h2 - h2_to_k1)

            t_min = min(t(i,kb2), t(i,kb1)) - cs%Allowed_T_chg
            t_max = max(t(i,kb2), t(i,kb1)) + cs%Allowed_T_chg
            t_new = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
            if (t_new < t_min) then
              h2_to_k1_lim = h2 * (t(i,kb2) - t_min) / (t_det - t_min)
!               write(mesg,'("Low temperature limits det to ", &
!                    & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. T=", &
!                    & 5(1pe12.5))') &
!                    h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
!                    T_new, T(i,kb2), T(i,kb1), T_det, T_new-T_min
!               call MOM_error(WARNING, mesg)
              h2_to_k1 = h2_to_k1_lim
              ih2f = 1.0 / (h2 - h2_to_k1)
            elseif (t_new > t_max) then
              h2_to_k1_lim = h2 * (t(i,kb2) - t_max) / (t_det - t_max)
!               write(mesg,'("High temperature limits det to ", &
!                    & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. T=", &
!                    & 5(1pe12.5))') &
!                    h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
!                    T_new, T(i,kb2), T(i,kb1), T_det, T_new-T_max
!               call MOM_error(WARNING, mesg)
              h2_to_k1 = h2_to_k1_lim
              ih2f = 1.0 / (h2 - h2_to_k1)
            endif
            s_min = max(min(s(i,kb2), s(i,kb1)) - cs%Allowed_S_chg, 0.0)
            s_max = max(s(i,kb2), s(i,kb1)) + cs%Allowed_S_chg
            s_new = (h2*s(i,kb2) - h2_to_k1*s_det)*ih2f
            if (s_new < s_min) then
              h2_to_k1_lim = h2 * (s(i,kb2) - s_min) / (s_det - s_min)
!               write(mesg,'("Low salinity limits det to ", &
!                    & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. S=", &
!                    & 5(1pe12.5))') &
!                    h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
!                    S_new, S(i,kb2), S(i,kb1), S_det, S_new-S_min
!               call MOM_error(WARNING, mesg)
              h2_to_k1 = h2_to_k1_lim
              ih2f = 1.0 / (h2 - h2_to_k1)
            elseif (s_new > s_max) then
              h2_to_k1_lim = h2 * (s(i,kb2) - s_max) / (s_det - s_max)
!               write(mesg,'("High salinity limits det to ", &
!                    & 1pe12.5, " from ", 1pe12.5, " at ", 1pg11.4,"E, ",1pg11.4,"N. S=", &
!                    & 5(1pe12.5))') &
!                    h2_to_k1_lim, h2_to_k1, G%geoLonT(i,j), G%geoLatT(i,j), &
!                    S_new, S(i,kb2), S(i,kb1), S_det, S_new-S_max
!               call MOM_error(WARNING, mesg)
              h2_to_k1 = h2_to_k1_lim
              ih2f = 1.0 / (h2 - h2_to_k1)
            endif

            ihk1 = 1.0 / (h(i,k1) + h_neglect + h2_to_k1)
            rcv(i,k1) = ((h(i,k1)+h_neglect)*rcv(i,k1) + h2_to_k1*rcvtgt(k1)) * ihk1
            rcv(i,kb2) = rcv(i,kb2) - h2_to_k1*drcv*ih2f

            t(i,kb2) = (h2*t(i,kb2) - h2_to_k1*t_det)*ih2f
            t(i,k1) = ((h(i,k1)+h_neglect)*t(i,k1) + h2_to_k1*t_det) * ihk1

            s(i,kb2) = (h2*s(i,kb2) - h2_to_k1*s_det) * ih2f
            s(i,k1) = ((h(i,k1)+h_neglect)*s(i,k1) + h2_to_k1*s_det) * ihk1

            ! Changes in R0 are based on changes in T and S.
            r0(i,kb2) = (h2*r0(i,kb2) - h2_to_k1*r0_det) * ih2f
            r0(i,k1) = ((h(i,k1)+h_neglect)*r0(i,k1) + h2_to_k1*r0_det) * ihk1
          else
            ! h2==h2_to_k1 can happen if dR2b = 0 exactly, but this is very
            ! unlikely.  In this case the entirety of layer kb2 is detrained.
            h2_to_k1 = h2  ! These 2 lines are probably unnecessary.
            ihk1 = 1.0 / (h(i,k1) + h2)

            rcv(i,k1) = (h(i,k1)*rcv(i,k1) + h2*rcv(i,kb2)) * ihk1
            t(i,k1) = (h(i,k1)*t(i,k1) + h2*t(i,kb2)) * ihk1
            s(i,k1) = (h(i,k1)*s(i,k1) + h2*s(i,kb2)) * ihk1
            r0(i,k1) = (h(i,k1)*r0(i,k1) + h2*r0(i,kb2)) * ihk1
          endif

          h(i,k1) = h(i,k1) + h2_to_k1
          h(i,kb2) = h(i,kb2) - h2_to_k1 ; h2 = h(i,kb2)
          ! dPE_extrap should be positive here.
          dpe_extrap = i2rho0*(r0_det-r0(i,kb2))*h2_to_k1*h2

          d_ea(i,kb2) = d_ea(i,kb2) - h2_to_k1
          d_ea(i,k1) = d_ea(i,k1) + h2_to_k1
          h2_to_k1_rem = max(h2_to_k1_rem - h2_to_k1, 0.0)
        endif
      endif ! Detrainment by extrapolation.

    endif ! Detrainment to the interior at all.

    ! Does some of the detrained water go into the lower buffer layer?
    h_det_h2 = max(h_min_bl-(h1+h2), 0.0)
    if (h_det_h2 > 0.0) then
      ! Detrained water will go into both upper and lower buffer layers.
      ! h(kb2) will be h_min_bl, but h(kb1) may be larger if there was already
      ! ample detrainment; all water in layer kb1 moves into layer kb2.

      ! Determine the fluxes between the various layers.
      h_det_to_h2 = min(h_to_bl, h_det_h2)
      h_ml_to_h2 = h_det_h2 - h_det_to_h2
      h_det_to_h1 = h_to_bl - h_det_to_h2
      h_ml_to_h1 = max(h_min_bl-h_det_to_h1,0.0)

      ih = 1.0/h_min_bl;
      ihdet = 0.0 ; if (h_to_bl > 0.0) ihdet = 1.0 / h_to_bl
      ih1f = 1.0 / (h_det_to_h1 + h_ml_to_h1)

      r0(i,kb2) = ((h2*r0(i,kb2) + h1*r0(i,kb1)) + &
                   (h_det_to_h2*r0_to_bl*ihdet + h_ml_to_h2*r0(i,0))) * ih
      r0(i,kb1) = (h_det_to_h1*r0_to_bl*ihdet + h_ml_to_h1*r0(i,0)) * ih1f

      rcv(i,kb2) = ((h2*rcv(i,kb2) + h1*rcv(i,kb1)) + &
                    (h_det_to_h2*rcv_to_bl*ihdet + h_ml_to_h2*rcv(i,0))) * ih
      rcv(i,kb1) = (h_det_to_h1*rcv_to_bl*ihdet + h_ml_to_h1*rcv(i,0)) * ih1f

      t(i,kb2) = ((h2*t(i,kb2) + h1*t(i,kb1)) + &
                  (h_det_to_h2*t_to_bl*ihdet + h_ml_to_h2*t(i,0))) * ih
      t(i,kb1) = (h_det_to_h1*t_to_bl*ihdet + h_ml_to_h1*t(i,0)) * ih1f

      s(i,kb2) = ((h2*s(i,kb2) + h1*s(i,kb1)) + &
                  (h_det_to_h2*s_to_bl*ihdet + h_ml_to_h2*s(i,0))) * ih
      s(i,kb1) = (h_det_to_h1*s_to_bl*ihdet + h_ml_to_h1*s(i,0)) * ih1f

      ! Recall that h1 = h(i,kb1) & h2 = h(i,kb2).
      d_ea(i,1) = d_ea(i,1) - (h_ml_to_h1 + h_ml_to_h2)
      d_ea(i,kb1) = d_ea(i,kb1) + ((h_det_to_h1 + h_ml_to_h1) - h1)
      d_ea(i,kb2) = d_ea(i,kb2) + (h_min_bl - h2)

      h(i,kb1) = h_det_to_h1 + h_ml_to_h1 ; h(i,kb2) = h_min_bl
      h(i,0) = h(i,0) - (h_ml_to_h1 + h_ml_to_h2)


      if (ALLOCATED(cs%diag_PE_detrain) .or. ALLOCATED(cs%diag_PE_detrain2)) then
        r0_det = r0_to_bl*ihdet
        s1en = g_2 * idt_h2 * ( ((r0(i,kb2)-r0(i,kb1))*h1*h2 + &
            h_det_to_h2*( (r0(i,kb1)-r0_det)*h1 + (r0(i,kb2)-r0_det)*h2 ) + &
            h_ml_to_h2*( (r0(i,kb2)-r0(i,0))*h2 + (r0(i,kb1)-r0(i,0))*h1 + &
                         (r0_det-r0(i,0))*h_det_to_h2 ) + &
            h_det_to_h1*h_ml_to_h1*(r0_det-r0(i,0))) - 2.0*gv%Rho0*dpe_extrap )

        if (ALLOCATED(cs%diag_PE_detrain)) &
          cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + s1en

        if (ALLOCATED(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
            cs%diag_PE_detrain2(i,j) + s1en + idt_h2*rho0xg*dpe_extrap
      endif

    elseif ((h_to_bl > 0.0) .or. (h1 < h_min_bl) .or. (h2 < h_min_bl)) then
    ! Determine how much of the upper buffer layer will be moved into
    ! the lower buffer layer and the properties with which it is moving.
    ! This implementation assumes a 2nd-order upwind advection of density
    ! from the uppermost buffer layer into the next one down.
      h_from_ml = h_min_bl + max(h_min_bl-h2,0.0) - h1 - h_to_bl
      if (h_from_ml > 0.0) then
        ! Some water needs to be moved from the mixed layer so that the upper
        ! (and perhaps lower) buffer layers exceed their minimum thicknesses.
        dpe_extrap = dpe_extrap - i2rho0*h_from_ml*(r0_to_bl - r0(i,0)*h_to_bl)
        r0_to_bl = r0_to_bl + h_from_ml*r0(i,0)
        rcv_to_bl = rcv_to_bl + h_from_ml*rcv(i,0)
        t_to_bl = t_to_bl + h_from_ml*t(i,0)
        s_to_bl = s_to_bl + h_from_ml*s(i,0)

        h_to_bl = h_to_bl + h_from_ml
        h(i,0) = h(i,0) - h_from_ml
        d_ea(i,1) = d_ea(i,1) - h_from_ml
      endif

      ! The absolute value should be unnecessary and 1e9 is just a large number.
      b1 = 1.0e9
      if (r0(i,kb2) - r0(i,kb1) > 1.0e-9*abs(r0(i,kb1) - r0_det)) &
        b1 = abs(r0(i,kb1) - r0_det) / (r0(i,kb2) - r0(i,kb1))
      stays_min = max((1.0-b1)*h1 - b1*h2, 0.0, h_min_bl - h_to_bl)
      stays_max = h1 - max(h_min_bl-h2,0.0)

      scale_slope = 1.0
      if (stays_max <= stays_min) then
        stays = stays_max
        mergeable_bl = .false.
        if (stays_max < h1) scale_slope = (h1 - stays_min) / (h1 - stays_max)
      else
        ! There are numerous temporary variables used here that should not be
        ! used outside of this "else" branch: s1, s2, s3sq, I_ya, bh0
        bh0 = b1*h_to_bl
        i_ya =  (h1 + h2) / ((h1 + h2) + h_to_bl)
        ! s1 is the amount staying that minimizes the PE increase.
        s1 = 0.5*(h1 + (h2 - bh0) * i_ya) ; s2 = h1 - s1

        if (s2 < 0.0) then
          ! The energy released by detrainment from the lower buffer layer can be
          ! used to mix water from the upper buffer layer into the lower one.
          s3sq = i_ya*max(bh0*h1-dpe_extrap, 0.0)
        else
          s3sq = i_ya*(bh0*h1-min(dpe_extrap,0.0))
        endif

        if (s3sq == 0.0) then
          ! There is a simple, exact solution to the quadratic equation, namely:
          stays = h1 ! This will revert to stays_max later.
        elseif (s2*s2 <= s3sq) then
          ! There is no solution with 0 PE change - use the minimum energy input.
          stays = s1
        else
          ! The following choose the solutions that are continuous with all water
          ! staying in the upper buffer layer when there is no detrainment,
          ! namely the + root when s2>0 and the - root otherwise. They also
          ! carefully avoid differencing large numbers, using s2 = (h1-s).
          if (bh0 <= 0.0) then ; stays = h1
          elseif (s2 > 0.0) then
  !         stays = s + sqrt(s2*s2 - s3sq) ! Note that s2 = h1-s
            if (s1 >= stays_max) then ; stays = stays_max
            elseif (s1 >= 0.0) then ; stays = s1 + sqrt(s2*s2 - s3sq)
            else ; stays = (h1*(s2-s1) - s3sq) / (-s1 + sqrt(s2*s2 - s3sq))
            endif
          else
  !         stays = s - sqrt(s2*s2 - s3sq) ! Note that s2 = h1-s & stays_min >= 0
            if (s1 <= stays_min) then ; stays = stays_min
            else ; stays = (h1*(s1-s2) + s3sq) / (s1 + sqrt(s2*s2 - s3sq))
            endif
          endif
        endif

        ! Limit the amount that stays so that the motion of water is from the
        ! upper buffer layer into the lower, but no more than is in the upper
        ! layer, and the water left in the upper layer is no lighter than the
        ! detrained water.
        if (stays >= stays_max) then ; stays = stays_max
        elseif (stays < stays_min) then ; stays = stays_min
        endif
      endif

      dpe_det = g_2*((r0(i,kb1)*h_to_bl - r0_to_bl)*stays + &
                     (r0(i,kb2)-r0(i,kb1)) * (h1-stays) * &
                     (h2 - scale_slope*stays*((h1+h2)+h_to_bl)/(h1+h2)) ) - &
                rho0xg*dpe_extrap

      if (dpe_time_ratio*h_to_bl > h_to_bl+h(i,0)) then
        dpe_ratio = (h_to_bl+h(i,0)) / h_to_bl
      else
        dpe_ratio = dpe_time_ratio
      endif

      if ((mergeable_bl) .and. (num_events*dpe_ratio*dpe_det > dpe_merge)) then
        ! It is energetically preferable to merge the two buffer layers, detrain
        ! them into interior layer (k0), move the remaining upper buffer layer
        ! water into the lower buffer layer, and detrain undiluted into the
        ! upper buffer layer.
        h1_to_k0 = (h1-stays_merge)
        stays = max(h_min_bl-h_to_bl,0.0)
        h1_to_h2 = stays_merge - stays

        ihk0 = 1.0 / ((h1_to_k0 + h2) + h(i,k0))
        ih1f = 1.0 / (h_to_bl + stays); ih2f = 1.0 / h1_to_h2
        ih12 = 1.0 / (h1 + h2)

        drcv_2dz = (rcv(i,kb1) - rcv(i,kb2)) * ih12
        drcv_stays = drcv_2dz*(h1_to_k0 + h1_to_h2)
        drcv_det = - drcv_2dz*(stays + h1_to_h2)
        rcv(i,k0) = ((h1_to_k0*(rcv(i,kb1) + drcv_det) + &
                      h2*rcv(i,kb2)) + h(i,k0)*rcv(i,k0)) * ihk0
        rcv(i,kb2) = rcv(i,kb1) + drcv_2dz*(h1_to_k0-stays)
        rcv(i,kb1) = (rcv_to_bl + stays*(rcv(i,kb1) + drcv_stays)) * ih1f

        ! Use 2nd order upwind advection of spiciness, limited by the value in
        ! the water from the mixed layer to determine the temperature and
        ! salinity of the water that stays in the buffer layers.
        i_denom = 1.0 / (drcv_ds(i)**2 + (dt_ds_gauge*drcv_dt(i))**2)
        dspice_2dz = (ds_dt_gauge*drcv_ds(i)*(t(i,kb1)-t(i,kb2)) - &
                      dt_ds_gauge*drcv_dt(i)*(s(i,kb1)-s(i,kb2))) * ih12
        dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
                      dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) / h_to_bl
        if (dspice_lim * dspice_2dz <= 0.0) dspice_2dz = 0.0

        if (stays > 0.0) then
        ! Limit the spiciness of the water that stays in the upper buffer layer.
          if (abs(dspice_lim) < abs(dspice_2dz*(h1_to_k0 + h1_to_h2))) &
            dspice_2dz = dspice_lim/(h1_to_k0 + h1_to_h2)

          dspice_stays = dspice_2dz*(h1_to_k0 + h1_to_h2)
          t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
              (dt_ds_gauge * drcv_dt(i) * drcv_stays + drcv_ds(i) * dspice_stays)
          s_stays = s(i,kb1) + i_denom * &
              (drcv_ds(i) * drcv_stays - dt_ds_gauge * drcv_dt(i) * dspice_stays)
          ! The values of R0 are based on changes in T and S.
          r0_stays = r0(i,kb1) + (t_stays-t(i,kb1)) * dr0_dt(i) + &
                                 (s_stays-s(i,kb1)) * dr0_ds(i)
        else
          ! Limit the spiciness of the water that moves into the lower buffer layer.
          if (abs(dspice_lim) < abs(dspice_2dz*h1_to_k0)) &
            dspice_2dz = dspice_lim/h1_to_k0
          ! These will be multiplied by 0 later.
          t_stays = 0.0 ; s_stays = 0.0 ; r0_stays = 0.0
        endif

        dspice_det = - dspice_2dz*(stays + h1_to_h2)
        t_det = t(i,kb1) + dt_ds_gauge * i_denom * &
            (dt_ds_gauge * drcv_dt(i) * drcv_det + drcv_ds(i) * dspice_det)
        s_det = s(i,kb1) + i_denom * &
            (drcv_ds(i) * drcv_det - dt_ds_gauge * drcv_dt(i) * dspice_det)
        ! The values of R0 are based on changes in T and S.
        r0_det = r0(i,kb1) + (t_det-t(i,kb1)) * dr0_dt(i) + &
                             (s_det-s(i,kb1)) * dr0_ds(i)

        t(i,k0) = ((h1_to_k0*t_det + h2*t(i,kb2)) + h(i,k0)*t(i,k0)) * ihk0
        t(i,kb2) = (h1*t(i,kb1) - stays*t_stays - h1_to_k0*t_det) * ih2f
        t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f

        s(i,k0) = ((h1_to_k0*s_det + h2*s(i,kb2)) + h(i,k0)*s(i,k0)) * ihk0
        s(i,kb2) = (h1*s(i,kb1) - stays*s_stays - h1_to_k0*s_det) * ih2f
        s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f

        r0(i,k0) = ((h1_to_k0*r0_det + h2*r0(i,kb2)) + h(i,k0)*r0(i,k0)) * ihk0
        r0(i,kb2) = (h1*r0(i,kb1) - stays*r0_stays - h1_to_k0*r0_det) * ih2f
        r0(i,kb1) = (r0_to_bl + stays*r0_stays) * ih1f

!        ! The following is 2nd-order upwind advection without limiters.
!        dT_2dz = (T(i,kb1) - T(i,kb2)) * Ih12
!        T(i,k0) = (h1_to_k0*(T(i,kb1) - dT_2dz*(stays+h1_to_h2)) + &
!                     h2*T(i,kb2) + h(i,k0)*T(i,k0)) * Ihk0
!        T(i,kb2) = T(i,kb1) + dT_2dz*(h1_to_k0-stays)
!        T(i,kb1) = (T_to_bl + stays*(T(i,kb1) + &
!                      dT_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
!        dS_2dz = (S(i,kb1) - S(i,kb2)) * Ih12
!        S(i,k0) = (h1_to_k0*(S(i,kb1) - dS_2dz*(stays+h1_to_h2)) + &
!                     h2*S(i,kb2) + h(i,k0)*S(i,k0)) * Ihk0
!        S(i,kb2) = S(i,kb1) + dS_2dz*(h1_to_k0-stays)
!        S(i,kb1) = (S_to_bl + stays*(S(i,kb1) + &
!                      dS_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f
!        dR0_2dz = (R0(i,kb1) - R0(i,kb2)) * Ih12
!        R0(i,k0) = (h1_to_k0*(R0(i,kb1) - dR0_2dz*(stays+h1_to_h2)) + &
!                    h2*R0(i,kb2) + h(i,k0)*R0(i,k0)) * Ihk0
!        R0(i,kb2) = R0(i,kb1) + dR0_2dz*(h1_to_k0-stays)
!        R0(i,kb1) = (R0_to_bl + stays*(R0(i,kb1) + &
!                     dR0_2dz*(h1_to_k0 + h1_to_h2))) * Ih1f

        d_ea(i,kb1) = (d_ea(i,kb1) + h_to_bl) + (stays - h1)
        d_ea(i,kb2) = d_ea(i,kb2) + (h1_to_h2 - h2)
        d_ea(i,k0) = d_ea(i,k0) + (h1_to_k0 + h2)

        h(i,kb1) = stays + h_to_bl
        h(i,kb2) = h1_to_h2
        h(i,k0) = h(i,k0) + (h1_to_k0 + h2)
        if (ALLOCATED(cs%diag_PE_detrain)) &
          cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_merge
        if (ALLOCATED(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
             cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det+rho0xg*dpe_extrap)
      else ! Not mergeable_bl.
        ! There is no further detrainment from the buffer layers, and the
        ! upper buffer layer water is distributed optimally between the
        ! upper and lower buffer layer.
        h1_to_h2 = h1 - stays
        ih1f = 1.0 / (h_to_bl + stays) ; ih2f = 1.0 / (h2 + h1_to_h2)
        ih = 1.0 / (h1 + h2)
        dr0_2dz = (r0(i,kb1) - r0(i,kb2)) * ih
        r0(i,kb2) = (h2*r0(i,kb2) + h1_to_h2*(r0(i,kb1) - &
                     scale_slope*dr0_2dz*stays)) * ih2f
        r0(i,kb1) = (r0_to_bl + stays*(r0(i,kb1) + &
                        scale_slope*dr0_2dz*h1_to_h2)) * ih1f

        ! Use 2nd order upwind advection of spiciness, limited by the value
        ! in the detrained water to determine the detrained temperature and
        ! salinity.
        dr0 = scale_slope*dr0_2dz*h1_to_h2
        dspice_stays = (ds_dt_gauge*dr0_ds(i)*(t(i,kb1)-t(i,kb2)) - &
                        dt_ds_gauge*dr0_dt(i)*(s(i,kb1)-s(i,kb2))) * &
                        scale_slope*h1_to_h2 * ih
        if (h_to_bl > 0.0) then
          dspice_lim = (ds_dt_gauge*dr0_ds(i)*(t_to_bl-t(i,kb1)*h_to_bl) - &
                        dt_ds_gauge*dr0_dt(i)*(s_to_bl-s(i,kb1)*h_to_bl)) /&
                        h_to_bl
        else
          dspice_lim = ds_dt_gauge*dr0_ds(i)*(t(i,0)-t(i,kb1)) - &
                       dt_ds_gauge*dr0_dt(i)*(s(i,0)-s(i,kb1))
        endif
        if (dspice_stays*dspice_lim <= 0.0) then
          dspice_stays = 0.0
        elseif (abs(dspice_stays) > abs(dspice_lim)) then
          dspice_stays = dspice_lim
        endif
        i_denom = 1.0 / (dr0_ds(i)**2 + (dt_ds_gauge*dr0_dt(i))**2)
        t_stays = t(i,kb1) + dt_ds_gauge * i_denom * &
            (dt_ds_gauge * dr0_dt(i) * dr0 + dr0_ds(i) * dspice_stays)
        s_stays = s(i,kb1) + i_denom * &
            (dr0_ds(i) * dr0 - dt_ds_gauge * dr0_dt(i) * dspice_stays)
        ! The detrained values of Rcv are based on changes in T and S.
        rcv_stays = rcv(i,kb1) + (t_stays-t(i,kb1)) * drcv_dt(i) + &
                                 (s_stays-s(i,kb1)) * drcv_ds(i)

        t(i,kb2) = (h2*t(i,kb2) + h1*t(i,kb1) - t_stays*stays) * ih2f
        t(i,kb1) = (t_to_bl + stays*t_stays) * ih1f
        s(i,kb2) = (h2*s(i,kb2) + h1*s(i,kb1) - s_stays*stays) * ih2f
        s(i,kb1) = (s_to_bl + stays*s_stays) * ih1f
        rcv(i,kb2) = (h2*rcv(i,kb2) + h1*rcv(i,kb1) - rcv_stays*stays) * ih2f
        rcv(i,kb1) = (rcv_to_bl + stays*rcv_stays) * ih1f

!        ! The following is 2nd-order upwind advection without limiters.
!        dRcv_2dz = (Rcv(i,kb1) - Rcv(i,kb2)) * Ih
!        dRcv = scale_slope*dRcv_2dz*h1_to_h2
!        Rcv(i,kb2) = (h2*Rcv(i,kb2) + h1_to_h2*(Rcv(i,kb1) - &
!                      scale_slope*dRcv_2dz*stays)) * Ih2f
!        Rcv(i,kb1) = (Rcv_to_bl + stays*(Rcv(i,kb1) + dRcv)) * Ih1f
!        dT_2dz = (T(i,kb1) - T(i,kb2)) * Ih
!        T(i,kb2) = (h2*T(i,kb2) + h1_to_h2*(T(i,kb1) - &
!                      scale_slope*dT_2dz*stays)) * Ih2f
!        T(i,kb1) = (T_to_bl + stays*(T(i,kb1) + &
!                      scale_slope*dT_2dz*h1_to_h2)) * Ih1f
!        dS_2dz = (S(i,kb1) - S(i,kb2)) * Ih
!        S(i,kb2) = (h2*S(i,kb2) + h1_to_h2*(S(i,kb1) - &
!                      scale_slope*dS_2dz*stays)) * Ih2f
!        S(i,kb1) = (S_to_bl + stays*(S(i,kb1) + &
!                      scale_slope*dS_2dz*h1_to_h2)) * Ih1f

        d_ea(i,kb1) = d_ea(i,kb1) + ((stays - h1) + h_to_bl)
        d_ea(i,kb2) = d_ea(i,kb2) + h1_to_h2

        h(i,kb1) = stays + h_to_bl
        h(i,kb2) = h(i,kb2) + h1_to_h2

        if (ALLOCATED(cs%diag_PE_detrain)) &
          cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + idt_h2*dpe_det
        if (ALLOCATED(cs%diag_PE_detrain2)) cs%diag_PE_detrain2(i,j) = &
          cs%diag_PE_detrain2(i,j) + idt_h2*(dpe_det+rho0xg*dpe_extrap)
      endif
    endif ! End of detrainment...

  enddo ! i loop

end subroutine mixedlayer_detrain_2

!> This subroutine moves any water left in the former mixed layers into the
!! single buffer layers and may also move buffer layer water into the interior
!! isopycnal layers.
subroutine mixedlayer_detrain_1(h, T, S, R0, Rcv, RcvTgt, dt, dt_diag, d_ea, d_eb, &
                                j, G, GV, CS, dRcv_dT, dRcv_dS, max_BL_det)
  type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure.
  type(verticalgrid_type),            intent(in)    :: GV   !< The ocean's vertical grid structure.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: h    !< Layer thickness, in m or kg m-2.
                                                            !! (Intent in/out) The units of h are
                                                            !! referred to as H below. Layer 0 is
                                                            !! the new mixed layer.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T    !< Potential temperature, in C.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S    !< Salinity, in psu.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0   !< Potential density referenced to
                                                            !! surface pressure, in kg m-3.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv  !< The coordinate defining potential
                                                            !! density, in kg m-3.
  real, dimension(SZK_(GV)),          intent(in)    :: RcvTgt !< The target value of Rcv for each
                                                            !! layer, in kg m-3.
  real,                               intent(in)    :: dt   !< Time increment, in s.
  real,                               intent(in)    :: dt_diag
  real, dimension(SZI_(G),SZK_(GV)),  intent(inout) :: d_ea !< The upward increase across a layer in
                                                            !! the entrainment from above, in m or
                                                            !! kg m-2 (H). Positive d_ea goes with
                                                            !! layer thickness increases.
  real, dimension(SZI_(G),SZK_(GV)),  intent(inout) :: d_eb !< The downward increase across a layer
                                                            !! in the entrainment from below, in H.
                                                            !! Positive values go with mass gain by
                                                            !! a layer.
  integer,                            intent(in)    :: j    !< The meridional row to work on.
  type(bulkmixedlayer_cs),            pointer       :: CS   !< The control structure returned by a
                                                            !! previous call to mixedlayer_init.
  real, dimension(SZI_(G)),           intent(in)    :: dRcv_dT !< The partial derivative of
                                                            !! coordinate defining potential density
                                                            !! with potential temperature,
                                                            !! in kg m-3 K-1.
  real, dimension(SZI_(G)),           intent(in)    :: dRcv_dS    !< The partial derivative of
                                                            !! coordinate defining potential density
                                                            !! with salinity, in kg m-3 psu-1.
  real, dimension(SZI_(G)),           intent(in)    :: max_BL_det !< If non-negative, the maximum
                                                            !! detrainment permitted from the buffer
                                                            !! layers, in H.

! This subroutine moves any water left in the former mixed layers into the
! single buffer layers and may also move buffer layer water into the interior
! isopycnal layers.

! Arguments: h - Layer thickness, in m or kg m-2. (Intent in/out)  The units of
!                h are referred to as H below.  Layer 0 is the new mixed layer.
!  (in/out)  T - Potential temperature, in C.
!  (in/out)  S - Salinity, in psu.
!  (in/out)  R0 - Potential density referenced to surface pressure, in kg m-3.
!  (in/out)  Rcv - The coordinate defining potential density, in kg m-3.
!  (in)      RcvTgt - The target value of Rcv for each layer, in kg m-3.
!  (in)      dt - Time increment, in s.
!  (in/out)  d_ea - The upward increase across a layer in the entrainment from
!                   above, in m or kg m-2 (H).  Positive d_ea goes with layer
!                   thickness increases.
!  (in/out)  d_eb - The downward increase across a layer in the entrainment from
!                   below, in H. Positive values go with mass gain by a layer.
!  (in)      j - The meridional row to work on.
!  (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)      max_BL_det - If non-negative, the maximum detrainment permitted
!                         from the buffer layers, in H.
!  (in/out)  dRcv_dT - The partial derivative of coordinate defining potential
!                      density with potential temperature, in kg m-3 K-1.
!  (in/out)  dRcv_dS - The partial derivative of coordinate defining potential
!                      density with salinity, in kg m-3 psu-1.
  real :: Ih                  ! The inverse of a thickness, in H-1.
  real :: h_ent               ! The thickness from a layer that is
                              ! entrained, in H.
  real :: max_det_rem(szi_(g)) ! Remaining permitted detrainment, in H.
  real :: detrain(szi_(g))    ! The thickness of fluid to detrain
                              ! from the mixed layer, in H.
  real :: Idt                 ! The inverse of the timestep in s-1.
  real :: dT_dR, dS_dR, dRml, dR0_dRcv, dT_dS_wt2
  real :: I_denom             ! A work variable with units of psu2 m6 kg-2.
  real :: Sdown, Tdown
  real :: dt_Time, Timescale = 86400.0*30.0! *365.0/12.0
  real :: g_H2_2Rho0dt       ! Half the gravitational acceleration times the
                              ! square of the conversion from H to m divided
                              ! by the mean density times the time step, in m6 s-3 H-2 kg-1.
  real :: g_H2_2dt            ! Half the gravitational acceleration times the
                              ! square of the conversion from H to m divided
                              ! by the diagnostic time step, in m3 H-2 s-3.

  logical :: splittable_BL(szi_(g)), orthogonal_extrap
  real :: x1

  integer :: i, is, ie, k, k1, nkmb, nz
  is = g%isc ; ie = g%iec ; nz = gv%ke
  nkmb = cs%nkml+cs%nkbl
  if (cs%nkbl /= 1) call mom_error(fatal,"MOM_mixed_layer: "// &
                        "CS%nkbl must be 1 in mixedlayer_detrain_1.")
  idt = 1.0/dt
  dt_time = dt/timescale
  g_h2_2rho0dt = (gv%g_Earth * gv%H_to_m**2) / (2.0 * gv%Rho0 * dt_diag)
  g_h2_2dt = (gv%g_Earth * gv%H_to_m**2) / (2.0 * dt_diag)

  ! Move detrained water into the buffer layer.
  do k=1,cs%nkml
    do i=is,ie ; if (h(i,k) > 0.0) then
      ih = 1.0 / (h(i,nkmb) + h(i,k))
      if (cs%TKE_diagnostics) &
        cs%diag_TKE_conv_s2(i,j) =  cs%diag_TKE_conv_s2(i,j) + &
          g_h2_2rho0dt * h(i,k) * h(i,nkmb) * &
          (r0(i,nkmb) - r0(i,k))
      if (ALLOCATED(cs%diag_PE_detrain)) &
        cs%diag_PE_detrain(i,j) = cs%diag_PE_detrain(i,j) + &
            g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))
      if (ALLOCATED(cs%diag_PE_detrain2)) &
        cs%diag_PE_detrain2(i,j) = cs%diag_PE_detrain2(i,j) + &
            g_h2_2dt * h(i,k) * h(i,nkmb) * (r0(i,nkmb) - r0(i,k))

      r0(i,nkmb) = (r0(i,nkmb)*h(i,nkmb) + r0(i,k)*h(i,k)) * ih
      rcv(i,nkmb) = (rcv(i,nkmb)*h(i,nkmb) + rcv(i,k)*h(i,k)) * ih
      t(i,nkmb) = (t(i,nkmb)*h(i,nkmb) + t(i,k)*h(i,k)) * ih
      s(i,nkmb) = (s(i,nkmb)*h(i,nkmb) + s(i,k)*h(i,k)) * ih

      d_ea(i,k) = d_ea(i,k) - h(i,k)
      d_ea(i,nkmb) = d_ea(i,nkmb) + h(i,k)
      h(i,nkmb) = h(i,nkmb) + h(i,k)
      h(i,k) = 0.0
    endif ; enddo
  enddo

  do i=is,ie
    max_det_rem(i) = 10.0 * h(i,nkmb)
    if (max_bl_det(i) >= 0.0) max_det_rem(i) = max_bl_det(i)
  enddo

!   If the mixed layer was denser than the densest interior layer,
! but is now lighter than this layer, leaving a buffer layer that
! is denser than this layer, there are problems.  This should prob-
! ably be considered a case of an inadequate choice of resolution in
! density space and should be avoided.  To make the model run sens-
! ibly in this case, it will make the mixed layer denser while making
! the buffer layer the density of the densest interior layer (pro-
! vided that the this will not make the mixed layer denser than the
! interior layer).  Otherwise, make the mixed layer the same density
! as the densest interior layer and lighten the buffer layer with
! the released buoyancy.  With multiple buffer layers, much more
! graceful options are available.
  do i=is,ie ; if (h(i,nkmb) > 0.0) then
    if ((r0(i,0)<r0(i,nz)) .and. (r0(i,nz)<r0(i,nkmb))) then
      if ((r0(i,nz)-r0(i,0))*h(i,0) > &
          (r0(i,nkmb)-r0(i,nz))*h(i,nkmb)) then
        detrain(i) = (r0(i,nkmb)-r0(i,nz))*h(i,nkmb) / (r0(i,nkmb)-r0(i,0))
      else
        detrain(i) = (r0(i,nz)-r0(i,0))*h(i,0) / (r0(i,nkmb)-r0(i,0))
      endif

      d_eb(i,cs%nkml) = d_eb(i,cs%nkml) + detrain(i)
      d_ea(i,cs%nkml) = d_ea(i,cs%nkml) - detrain(i)
      d_eb(i,nkmb) = d_eb(i,nkmb) - detrain(i)
      d_ea(i,nkmb) = d_ea(i,nkmb) + detrain(i)

      if (ALLOCATED(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
        cs%diag_PE_detrain(i,j) + g_h2_2dt * detrain(i)* &
                     (h(i,0) + h(i,nkmb)) * (r0(i,nkmb) - r0(i,0))
      x1 = r0(i,0)
      r0(i,0) = r0(i,0) - detrain(i)*(r0(i,0)-r0(i,nkmb)) / h(i,0)
      r0(i,nkmb) = r0(i,nkmb) - detrain(i)*(r0(i,nkmb)-x1) / h(i,nkmb)
      x1 = rcv(i,0)
      rcv(i,0) = rcv(i,0) - detrain(i)*(rcv(i,0)-rcv(i,nkmb)) / h(i,0)
      rcv(i,nkmb) = rcv(i,nkmb) - detrain(i)*(rcv(i,nkmb)-x1) / h(i,nkmb)
      x1 = t(i,0)
      t(i,0) = t(i,0) - detrain(i)*(t(i,0)-t(i,nkmb)) / h(i,0)
      t(i,nkmb) = t(i,nkmb) - detrain(i)*(t(i,nkmb)-x1) / h(i,nkmb)
      x1 = s(i,0)
      s(i,0) = s(i,0) - detrain(i)*(s(i,0)-s(i,nkmb)) / h(i,0)
      s(i,nkmb) = s(i,nkmb) - detrain(i)*(s(i,nkmb)-x1) / h(i,nkmb)

    endif
  endif ; enddo

  ! Move water out of the buffer layer, if convenient.
!   Split the buffer layer if possible, and replace the buffer layer
! with a small amount of fluid from the mixed layer.
! This is the exponential-in-time splitting, circa 2005.
  do i=is,ie
    if (h(i,nkmb) > 0.0) then ; splittable_bl(i) = .true.
    else ; splittable_bl(i) = .false. ; endif
  enddo

  dt_ds_wt2 = cs%dT_dS_wt**2

!    dt_Time = dt/Timescale
  do k=nz-1,nkmb+1,-1 ; do i=is,ie
    if (splittable_bl(i)) then
      if (rcvtgt(k)<=rcv(i,nkmb)) then
! Estimate dR/drho, dTheta/dR, and dS/dR, where R is the coordinate variable
! and rho is in-situ (or surface) potential density.
! There is no "right" way to do this, so this keeps things reasonable, if
! slightly arbitrary.
        splittable_bl(i) = .false.

        k1 = k+1 ; orthogonal_extrap = .false.
        ! Here we try to find a massive layer to use for interpolating the
        ! temperature and salinity.  If none is available a pseudo-orthogonal
        ! extrapolation is used.  The 10.0 and 0.9 in the following are
        ! arbitrary but probably about right.
        if ((h(i,k+1) < 10.0*gv%Angstrom) .or. &
            ((rcvtgt(k+1)-rcv(i,nkmb)) >= 0.9*(rcv(i,k1) - rcv(i,0)))) then
          if (k>=nz-1) then ; orthogonal_extrap = .true.
          elseif ((h(i,k+2) <= 10.0*gv%Angstrom) .and. &
              ((rcvtgt(k+1)-rcv(i,nkmb)) < 0.9*(rcv(i,k+2)-rcv(i,0)))) then
            k1 = k+2
          else ; orthogonal_extrap = .true. ; endif
        endif

        if ((r0(i,0) >= r0(i,k1)) .or. (rcv(i,0) >= rcv(i,nkmb))) cycle
          ! In this case there is an inversion of in-situ density relative to
          ! the coordinate variable.  Do not detrain from the buffer layer.

        if (orthogonal_extrap) then
          ! 36 here is a typical oceanic value of (dR/dS) / (dR/dT) - it says
          ! that the relative weights of T & S changes is a plausible 6:1.
          ! Also, this was coded on Athena's 6th birthday!
          i_denom = 1.0 / (drcv_ds(i)**2 + dt_ds_wt2*drcv_dt(i)**2)
          dt_dr = dt_ds_wt2*drcv_dt(i) * i_denom
          ds_dr = drcv_ds(i) * i_denom
        else
          dt_dr = (t(i,0) - t(i,k1)) / (rcv(i,0) - rcv(i,k1))
          ds_dr = (s(i,0) - s(i,k1)) / (rcv(i,0) - rcv(i,k1))
        endif
        drml = dt_time * (r0(i,nkmb) - r0(i,0)) * &
               (rcv(i,0) - rcv(i,k1)) / (r0(i,0) - r0(i,k1))
        ! Once again, there is an apparent density inversion in Rcv.
        if (drml < 0.0) cycle
        dr0_drcv = (r0(i,0) - r0(i,k1)) / (rcv(i,0) - rcv(i,k1))

        if ((rcv(i,nkmb) - drml < rcvtgt(k)) .and. (max_det_rem(i) > h(i,nkmb))) then
          ! In this case, the buffer layer is split into two isopycnal layers.
          detrain(i) = h(i,nkmb)*(rcv(i,nkmb) - rcvtgt(k)) / &
                                  (rcvtgt(k+1) - rcvtgt(k))

          if (ALLOCATED(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
            cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * &
                 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcvtgt(k)) * dr0_drcv

          tdown = detrain(i) * (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
          t(i,k) = (h(i,k) * t(i,k) + &
                        (h(i,nkmb) * t(i,nkmb) - tdown)) / &
                       (h(i,k) + (h(i,nkmb) - detrain(i)))
          t(i,k+1) = (h(i,k+1) * t(i,k+1) + tdown)/ &
                          (h(i,k+1) + detrain(i))
          t(i,nkmb) = t(i,0)
          sdown = detrain(i) * (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
          s(i,k) = (h(i,k) * s(i,k) + &
                      (h(i,nkmb) * s(i,nkmb) - sdown)) / &
                      (h(i,k) + (h(i,nkmb) - detrain(i)))
          s(i,k+1) = (h(i,k+1) * s(i,k+1) + sdown)/ &
                         (h(i,k+1) + detrain(i))
          s(i,nkmb) = s(i,0)
          rcv(i,nkmb) = rcv(i,0)

          d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
          d_ea(i,k) = d_ea(i,k) + (h(i,nkmb) - detrain(i))
          d_ea(i,nkmb) = d_ea(i,nkmb) - h(i,nkmb)

          h(i,k+1) = h(i,k+1) + detrain(i)
          h(i,k) = h(i,k) + h(i,nkmb) - detrain(i)
          h(i,nkmb) = 0.0
        else
          ! Here only part of the buffer layer is moved into the interior.
          detrain(i) = h(i,nkmb) * drml / (rcvtgt(k+1) - rcv(i,nkmb) + drml)
          if (detrain(i) > max_det_rem(i)) detrain(i) = max_det_rem(i)
          ih = 1.0 / (h(i,k+1) + detrain(i))

          tdown = (t(i,nkmb) + dt_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
          t(i,nkmb) = t(i,nkmb) - dt_dr * drml
          t(i,k+1) = (h(i,k+1) * t(i,k+1) + detrain(i) * tdown) * ih
          sdown = (s(i,nkmb) + ds_dr*(rcvtgt(k+1)-rcv(i,nkmb)))
!  The following two expressions updating S(nkmb) are mathematically identical.
!            S(i,nkmb) = (h(i,nkmb) * S(i,nkmb) - detrain(i) * Sdown) / &
!                           (h(i,nkmb) - detrain(i))
          s(i,nkmb) = s(i,nkmb) - ds_dr * drml
          s(i,k+1) = (h(i,k+1) * s(i,k+1) + detrain(i) * sdown) * ih

          d_ea(i,k+1) = d_ea(i,k+1) + detrain(i)
          d_ea(i,nkmb) = d_ea(i,nkmb) - detrain(i)

          h(i,k+1) = h(i,k+1) + detrain(i)
          h(i,nkmb) = h(i,nkmb) - detrain(i)

          if (ALLOCATED(cs%diag_PE_detrain)) cs%diag_PE_detrain(i,j) = &
            cs%diag_PE_detrain(i,j) - g_h2_2dt * detrain(i) * dr0_drcv * &
                 (h(i,nkmb)-detrain(i)) * (rcvtgt(k+1) - rcv(i,nkmb) + drml)
        endif
      endif ! RcvTgt(k)<=Rcv(i,nkmb)
    endif ! splittable_BL
  enddo ; enddo ! i & k loops

!   The numerical behavior of the buffer layer is dramatically improved
! if it is always at least a small fraction (say 10%) of the thickness
! of the mixed layer.  As the physical distinction between the mixed
! and buffer layers is vague anyway, this seems hard to argue against.
  do i=is,ie
    if (h(i,nkmb) < 0.1*h(i,0)) then
      h_ent =  0.1*h(i,0) - h(i,nkmb)
      ih = 10.0/h(i,0)
      t(i,nkmb) = (h(i,nkmb)*t(i,nkmb) + h_ent*t(i,0)) * ih
      s(i,nkmb) = (h(i,nkmb)*s(i,nkmb) + h_ent*s(i,0)) * ih

      d_ea(i,1) = d_ea(i,1) - h_ent
      d_ea(i,nkmb) = d_ea(i,nkmb) + h_ent

      h(i,0) = h(i,0) - h_ent
      h(i,nkmb) = h(i,nkmb) + h_ent
    endif
  enddo

end subroutine mixedlayer_detrain_1

! #@# This subroutine needs a doxygen description.
subroutine bulkmixedlayer_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 !< A structure that is used to regulate diagnostic
                                                 !! output.
  type(bulkmixedlayer_cs), pointer       :: CS   !< A pointer that is set to point to the control
                                                 !! structure for this module.
! 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_mixed_layer"  ! This module's name.
  real :: omega_frac_dflt, ustar_min_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

  if (gv%nkml < 1) return

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

  cs%nkml = gv%nkml
  call log_param(param_file, mdl, "NKML", cs%nkml, &
                 "The number of sublayers within the mixed layer if \n"//&
                 "BULKMIXEDLAYER is true.", units="nondim", default=2)
  cs%nkbl = gv%nk_rho_varies - gv%nkml
  call log_param(param_file, mdl, "NKBL", cs%nkbl, &
                 "The number of variable density buffer layers if \n"//&
                 "BULKMIXEDLAYER is true.", units="nondim", default=2)
  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, "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.15)
  call get_param(param_file, mdl, "BULK_RI_ML", cs%bulk_Ri_ML, &
                 "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",&
                 fail_if_missing=.true.)
  call get_param(param_file, mdl, "ABSORB_ALL_SW", cs%absorb_all_sw, &
                 "If true,  all shortwave radiation is absorbed by the \n"//&
                 "ocean, instead of passing through to the bottom mud.", &
                 default=.false.)
  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, "NSTAR2", cs%nstar2, &
                 "The portion of any potential energy released by \n"//&
                 "convective adjustment that is available to drive \n"//&
                 "entrainment at the base of mixed layer. By default \n"//&
                 "NSTAR2=NSTAR.", units="nondim", default=cs%nstar)
  call get_param(param_file, mdl, "BULK_RI_CONVECTIVE", cs%bulk_Ri_convective, &
                 "The efficiency with which convectively released mean \n"//&
                 "kinetic energy is converted to turbulent kinetic \n"//&
                 "energy.  By default BULK_RI_CONVECTIVE=BULK_RI_ML.", &
                 units="nondim", default=cs%bulk_Ri_ML)
  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, "LIMIT_BUFFER_DETRAIN", cs%limit_det, &
                 "If true, limit the detrainment from the buffer layers \n"//&
                 "to not be too different from the neighbors.", default=.false.)
  call get_param(param_file, mdl, "ALLOWED_DETRAIN_TEMP_CHG", cs%Allowed_T_chg, &
                 "The amount by which temperature is allowed to exceed \n"//&
                 "previous values during detrainment.", units="K", default=0.5)
  call get_param(param_file, mdl, "ALLOWED_DETRAIN_SALT_CHG", cs%Allowed_S_chg, &
                 "The amount by which salinity is allowed to exceed \n"//&
                 "previous values during detrainment.", units="PSU", default=0.1)
  call get_param(param_file, mdl, "ML_DT_DS_WEIGHT", cs%dT_dS_wt, &
                 "When forced to extrapolate T & S to match the layer \n"//&
                 "densities, this factor (in deg C / PSU) is combined \n"//&
                 "with the derivatives of density with T & S to determine \n"//&
                 "what direction is orthogonal to density contours. It \n"//&
                 "should be a typical value of (dR/dS) / (dR/dT) in \n"//&
                 "oceanic profiles.", units="degC PSU-1", default=6.0)
  call get_param(param_file, mdl, "BUFFER_LAYER_EXTRAP_LIMIT", cs%BL_extrap_lim, &
                 "A limit on the density range over which extrapolation \n"//&
                 "can occur when detraining from the buffer layers, \n"//&
                 "relative to the density range within the mixed and \n"//&
                 "buffer layers, when the detrainment is going into the \n"//&
                 "lightest interior layer, nondimensional, or a negative \n"//&
                 "value not to apply this limit.", units="nondim", default = -1.0)
  call get_param(param_file, mdl, "DEPTH_LIMIT_FLUXES", cs%H_limit_fluxes, &
                 "The surface fluxes are scaled away when the total ocean \n"//&
                 "depth is less than DEPTH_LIMIT_FLUXES.", &
                 units="m", default=0.1*cs%Hmix_min)
  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, "ML_RESORT", cs%ML_resort, &
                 "If true, resort the topmost layers by potential density \n"//&
                 "before the mixed layer calculations.", default=.false.)
  if (cs%ML_resort) &
    call get_param(param_file, mdl, "ML_PRESORT_NK_CONV_ADJ", cs%ML_presort_nz_conv_adj, &
                 "Convectively mix the first ML_PRESORT_NK_CONV_ADJ \n"//&
                 "layers before sorting when ML_RESORT is true.", &
                 units="nondim", default=0, fail_if_missing=.true.) ! Fail added by AJA.
  ! This gives a minimum decay scale that is typically much less than Angstrom.
  ustar_min_dflt = 2e-4*cs%omega*(gv%Angstrom_z + gv%H_to_m*gv%H_subroundoff)
  call get_param(param_file, mdl, "BML_USTAR_MIN", cs%ustar_min, &
                 "The minimum value of ustar that should be used by the \n"//&
                 "bulk mixed layer model in setting vertical TKE decay \n"//&
                 "scales. This must be greater than 0.", units="m s-1", &
                 default=ustar_min_dflt)
  if (cs%ustar_min<=0.0) call mom_error(fatal, "BML_USTAR_MIN must be positive.")

  call get_param(param_file, mdl, "RESOLVE_EKMAN", cs%Resolve_Ekman, &
                 "If true, the NKML>1 layers in the mixed layer are \n"//&
                 "chosen to optimally represent the impact of the Ekman \n"//&
                 "transport on the mixed layer TKE budget.  Otherwise, \n"//&
                 "the sublayers are distributed uniformly through the \n"//&
                 "mixed layer.", default=.false.)
  call get_param(param_file, mdl, "CORRECT_ABSORPTION_DEPTH", cs%correct_absorption, &
                 "If true, the average depth at which penetrating shortwave \n"//&
                 "radiation is absorbed is adjusted to match the average \n"//&
                 "heating depth of an exponential profile by moving some \n"//&
                 "of the heating upward in the water column.", default=.false.)
  call get_param(param_file, mdl, "DO_RIVERMIX", cs%do_rivermix, &
                 "If true, apply additional mixing whereever there is \n"//&
                 "runoff, so that it is mixed down to RIVERMIX_DEPTH, \n"//&
                 "if the ocean is that deep.", default=.false.)
  if (cs%do_rivermix) &
    call get_param(param_file, mdl, "RIVERMIX_DEPTH", cs%rivermix_depth, &
                 "The depth to which rivers are mixed if DO_RIVERMIX is \n"//&
                 "defined.", units="m", default=0.0)
  call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
                 "If true, use the fluxes%runoff_Hflx field to set the \n"//&
                 "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
                 default=.false.)
  call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
                 "If true, use the fluxes%calving_Hflx field to set the \n"//&
                 "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
                 default=.false.)

  call get_param(param_file, mdl, "ALLOW_CLOCKS_IN_OMP_LOOPS", &
                 cs%allow_clocks_in_omp_loops, &
                 "If true, clocks can be called from inside loops that can \n"//&
                 "be threaded. To run with multiple threads, set to False.", &
                 default=.true.)

  cs%id_ML_depth = register_diag_field('ocean_model', 'h_ML', diag%axesT1, &
      time, 'Surface mixed layer depth', 'meter')
  cs%id_TKE_wind = register_diag_field('ocean_model', 'TKE_wind', diag%axesT1, &
      time, 'Wind-stirring source of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_RiBulk = register_diag_field('ocean_model', 'TKE_RiBulk', diag%axesT1, &
      time, 'Mean kinetic energy source of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_conv = register_diag_field('ocean_model', 'TKE_conv', diag%axesT1, &
      time, 'Convective source of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_pen_SW = register_diag_field('ocean_model', 'TKE_pen_SW', diag%axesT1, &
      time, 'TKE consumed by mixing penetrative shortwave radation through the mixed layer', 'meter3 second-3')
  cs%id_TKE_mixing = register_diag_field('ocean_model', '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', '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', 'TKE_conv_decay', diag%axesT1, &
      time, 'Convective energy decay sink of mixed layer TKE', 'meter3 second-3')
  cs%id_TKE_conv_s2 = register_diag_field('ocean_model', 'TKE_conv_s2', diag%axesT1, &
      time, 'Spurious source of mixed layer TKE from sigma2', 'meter3 second-3')
  cs%id_PE_detrain = register_diag_field('ocean_model', 'PE_detrain', diag%axesT1, &
      time, 'Spurious source of potential energy from mixed layer detrainment', 'Watt meter-2')
  cs%id_PE_detrain2 = register_diag_field('ocean_model', 'PE_detrain2', diag%axesT1, &
      time, 'Spurious source of potential energy from mixed layer only detrainment', 'Watt meter-2')
  cs%id_h_mismatch = register_diag_field('ocean_model', 'h_miss_ML', diag%axesT1, &
      time, 'Summed absolute mismatch in entrainment terms', 'meter')
  cs%id_Hsfc_used = register_diag_field('ocean_model', 'Hs_used', diag%axesT1, &
      time, 'Surface region thickness that is used', 'meter')
  cs%id_Hsfc_max = register_diag_field('ocean_model', 'Hs_max', diag%axesT1, &
      time, 'Maximum surface region thickness', 'meter')
  cs%id_Hsfc_min = register_diag_field('ocean_model', 'Hs_min', diag%axesT1, &
      time, 'Minimum surface region thickness', 'meter')
 !CS%lim_det_dH_sfc = 0.5 ; CS%lim_det_dH_bathy = 0.2 ! Technically these should not get used if limit_det is false?
  if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) then
    call get_param(param_file, mdl, "LIMIT_BUFFER_DET_DH_SFC", cs%lim_det_dH_sfc, &
                 "The fractional limit in the change between grid points \n"//&
                 "of the surface region (mixed & buffer layer) thickness.", &
                 units="nondim", default=0.5)
    call get_param(param_file, mdl, "LIMIT_BUFFER_DET_DH_BATHY", cs%lim_det_dH_bathy, &
                 "The fraction of the total depth by which the thickness \n"//&
                 "of the surface region (mixed & buffer layer) is allowed \n"//&
                 "to change between grid points.", units="nondim", default=0.2)
  endif

  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
                 "If true, temperature and salinity are used as state \n"//&
                 "variables.", default=.true.)
  cs%nsw = 0
  if (use_temperature) then
    call get_param(param_file, mdl, "PEN_SW_NBANDS", cs%nsw, default=1)
  endif


  if (max(cs%id_TKE_wind, cs%id_TKE_RiBulk, cs%id_TKE_conv, &
          cs%id_TKE_mixing, cs%id_TKE_pen_SW, cs%id_TKE_mech_decay, &
          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_RiBulk, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%diag_TKE_conv, isd, ied, jsd, jed)
    call safe_alloc_alloc(cs%diag_TKE_pen_SW, 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)
    call safe_alloc_alloc(cs%diag_TKE_conv_s2, isd, ied, jsd, jed)

    cs%TKE_diagnostics = .true.
  endif
  if (cs%id_PE_detrain > 0) call safe_alloc_alloc(cs%diag_PE_detrain, isd, ied, jsd, jed)
  if (cs%id_PE_detrain2 > 0) call safe_alloc_alloc(cs%diag_PE_detrain2, isd, ied, jsd, jed)
  if (cs%id_ML_depth > 0) call safe_alloc_alloc(cs%ML_depth, isd, ied, jsd, jed)

  if(cs%allow_clocks_in_omp_loops) then
    id_clock_detrain = cpu_clock_id('(Ocean mixed layer detrain)', grain=clock_routine)
    id_clock_mech = cpu_clock_id('(Ocean mixed layer mechanical entrainment)', grain=clock_routine)
    id_clock_conv = cpu_clock_id('(Ocean mixed layer convection)', grain=clock_routine)
    if (cs%ML_resort) then
      id_clock_resort = cpu_clock_id('(Ocean mixed layer resorting)', grain=clock_routine)
    else
      id_clock_adjustment = cpu_clock_id('(Ocean mixed layer convective adjustment)', grain=clock_routine)
    endif
    id_clock_eos = cpu_clock_id('(Ocean mixed layer EOS)', grain=clock_routine)
  endif

  if (cs%limit_det .or. (cs%id_Hsfc_min > 0)) &
      id_clock_pass = cpu_clock_id('(Ocean mixed layer halo updates)', grain=clock_routine)


!  if (CS%limit_det) then
!  endif

end subroutine bulkmixedlayer_init

!> This subroutine returns an approximation to the integral
!!   R = exp(-L*(H+E)) integral(LH to L(H+E)) L/(1-(1+x)exp(-x)) dx.
!! The approximation to the integrand is good to within -2% at x~.3
!! and +25% at x~3.5, but the exponential deemphasizes the importance of
!! large x.  When L=0, EF4 returns E/((H+E)*H).
function ef4(H, E, L, dR_de)
real, intent(in)              :: H !< Total thickness, in m or kg m-2. (Intent in) The units of h
                                   !! are referred to as H below.
real, intent(in)              :: E !< Entrainment, in units of H.
real, intent(in)              :: L !< The e-folding scale in H-1.
real, intent(inout), optional :: dR_de !< The partial derivative of the result R with E, in H-2.
real :: EF4
! This subroutine returns an approximation to the integral
!   R = exp(-L*(H+E)) integral(LH to L(H+E)) L/(1-(1+x)exp(-x)) dx.
! The approximation to the integrand is good to within -2% at x~.3
! and +25% at x~3.5, but the exponential deemphasizes the importance of
! large x.  When L=0, EF4 returns E/((H+E)*H).
!
! Arguments: h - Total thickness, in m or kg m-2. (Intent in)  The units
!                of h are referred to as H below.
!  (in)      E - Entrainment, in units of H.
!  (in)      L - The e-folding scale in H-1.
!  (out)     dR_de - the partial derivative of the result R with E, in H-2.
! (return value) R - The integral, in units of H-1.
  real :: exp_LHpE ! A nondimensional exponential decay.
  real :: I_HpE    ! An inverse thickness plus entrainment, in H-1.
  real :: R        ! The result of the integral above, in H-1.

  exp_lhpe = exp(-l*(e+h))
  i_hpe = 1.0/(h+e)
  r = exp_lhpe * (e*i_hpe/h - 0.5*l*log(h*i_hpe) + 0.5*l*l*e)
  if (PRESENT(dr_de)) &
    dr_de = -l*r + exp_lhpe*(i_hpe*i_hpe + 0.5*l*i_hpe + 0.5*l*l)
  ef4 = r

end function ef4

end module mom_bulk_mixed_layer