MOM_bulk_mixed_layer.F90

!> Build mixed layer parameterization
module mom_bulk_mixed_layer

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

use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
use mom_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_opacity,       only : absorbremainingsw, optics_type, extract_optics_slice
use mom_unit_scaling,  only : unit_scale_type
use mom_variables,     only : thermo_var_ptrs
use mom_verticalgrid,  only : verticalgrid_type
use mom_eos,           only : calculate_density, calculate_density_derivs, eos_domain

implicit none ; private

#include <MOM_memory.h>

public bulkmixedlayer, bulkmixedlayer_init

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

!> The control structure with parameters for the MOM_bulk_mixed_layer module
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 [nondim].
  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 [H ~> m or kg m-2].
  real    :: h_limit_fluxes  !< When the total ocean depth is less than this
                             !! value [H ~> m or kg m-2], scale away all surface forcing to
                             !! avoid boiling the ocean.
  real    :: ustar_min       !< A minimum value of ustar to avoid numerical problems [Z T-1 ~> m s-1].
                             !! If the value is small enough, this should not affect the solution.
  real    :: omega           !<   The Earth's rotation rate [T-1 ~> s-1].
  real    :: dt_ds_wt        !<   When forced to extrapolate T & S to match the
                             !! layer densities, this factor (in degC / ppt) 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 degC ppt-1 might be reasonable.
  real    :: hbuffer_min     !< The minimum buffer layer thickness when the mixed layer
                             !! is very large [H ~> m or kg m-2].
  real    :: hbuffer_rel_min !< The minimum buffer layer thickness relative to the combined
                             !! mixed and buffer layer thicknesses when they are thin [nondim]
  real    :: bl_detrain_time !< A timescale that characterizes buffer layer detrainment
                             !! events [T ~> s].
  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  [nondim].
  real :: bl_split_rho_tol   !< The fractional tolerance for matching layer target densities
                             !! when splitting layers to deal with massive interior layers
                             !! that are lighter than one of the mixed or buffer layers [nondim].
  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 => null() !< A pointer to the ocean model's clock.
  logical :: tke_diagnostics = .false. !< If true, calculate extensive diagnostics of the TKE budget
  logical :: do_rivermix = .false. !< Provide additional TKE to mix river runoff
                             !! at the river mouths to rivermix_depth
  real    :: rivermix_depth = 0.0  !< The depth of mixing if do_rivermix is true [Z ~> m].
  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 :: convect_mom_bug !< If true, use code with a bug that causes a loss of momentum
                             !! conservation during mixedlayer convection.

  type(diag_ctrl), pointer :: diag => null() !< 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, ppt.

  ! These are terms in the mixed layer TKE budget, all in [Z L2 T-3 ~> m3 s-3] except as noted.
  real, allocatable, dimension(:,:) :: &
    ml_depth, &        !< The mixed layer depth [H ~> m or kg m-2].
    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 [R Z L2 T-3 ~> W m-2].
    diag_pe_detrain2   !< The spurious source of potential energy due to mixed layer only
                       !! detrainment [R Z L2 T-3 ~> 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

  !>@{ Diagnostic IDs
  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

!>@{ CPU clock IDs
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
!>@}

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, US, 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.
  type(unit_scale_type),      intent(in)    :: us     !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                              intent(inout) :: h_3d   !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                              intent(in)    :: u_3d   !< Zonal velocities interpolated to h points
                                                      !! [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                              intent(in)    :: v_3d   !< Zonal velocities interpolated to h points
                                                      !! [L T-1 ~> m s-1].
  type(thermo_var_ptrs),      intent(inout) :: tv     !< A structure containing pointers to any
                                                      !! available thermodynamic fields. Absent
                                                      !! fields have NULL ptrs.
  type(forcing),              intent(inout) :: fluxes !< A structure containing pointers to any
                                                      !! possible forcing fields.  Unused fields
                                                      !! have NULL ptrs.
  real,                       intent(in)    :: dt     !< Time increment [T ~> s].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                              intent(inout) :: ea     !< The amount of fluid moved downward into a
                                                      !! layer; this should be increased due to
                                                      !! mixed layer detrainment [H ~> m or kg m-2].
  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 [H ~> m or kg m-2].
  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 [m-1].
  real, dimension(:,:),       pointer       :: hml    !< Active mixed layer depth [Z ~> m].
  logical,                    intent(in)    :: aggregate_fw_forcing !< If true, the net incoming and
                                                     !! outgoing surface freshwater fluxes are
                                                     !! combined before being applied, instead of
                                                     !! being applied separately.
  real,             optional, intent(in)    :: dt_diag  !< The diagnostic time step,
                                                      !! which may be less than dt if there are
                                                      !! two callse to mixedlayer [T ~> s].
  logical,          optional, intent(in)    :: last_call !< if true, this is the last call
                                                      !! to mixedlayer in the current time step, so
                                                      !! diagnostics will be written. The default is
                                                      !! .true.

  ! Local variables
  real, dimension(SZI_(G),SZK_(GV)) :: &
    eaml, &     !   The amount of fluid moved downward into a layer due to mixed
                ! layer detrainment [H ~> m or kg m-2]. (I.e. entrainment from above.)
    ebml        !   The amount of fluid moved upward into a layer due to mixed
                ! layer detrainment [H ~> m or kg m-2]. (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 [H ~> m or kg m-2].
    t, &        !   The layer temperatures [degC].
    s, &        !   The layer salinities [ppt].
    r0, &       !   The potential density referenced to the surface [R ~> kg m-3].
    rcv         !   The coordinate variable potential density [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(GV)) :: &
    u, &        !   The zonal velocity [L T-1 ~> m s-1].
    v, &        !   The meridional velocity [L T-1 ~> m s-1].
    h_orig, &   !   The original thickness [H ~> m or kg m-2].
    d_eb, &     !   The downward increase across a layer in the entrainment from
                ! below [H ~> m or kg m-2].  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 [H ~> m or kg m-2].  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 [H ~> m or kg m-2].
  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 [Z ~> m].
  real, dimension(SZI_(G)) :: &
    tke, &      !   The turbulent kinetic energy available for mixing over a
                ! time step [Z L2 T-2 ~> m3 s-2].
    conv_en, &  !   The turbulent kinetic energy source due to mixing down to
                ! the depth of free convection [Z L2 T-2 ~> m3 s-2].
    htot, &     !   The total depth of the layers being considered for
                ! entrainment [H ~> m or kg m-2].
    r0_tot, &   !   The integrated potential density referenced to the surface
                ! of the layers which are fully entrained [H R ~> kg m-2 or kg2 m-5].
    rcv_tot, &  !   The integrated coordinate value potential density of the
                ! layers that are fully entrained [H R ~> kg m-2 or kg2 m-5].
    ttot, &     !   The integrated temperature of layers which are fully
                ! entrained [degC H ~> degC m or degC kg m-2].
    stot, &     !   The integrated salt of layers which are fully entrained
                ! [H ppt ~> m ppt or ppt kg m-2].
    uhtot, &    !   The depth integrated zonal and meridional velocities in the
    vhtot, &    ! mixed layer [H L T-1 ~> m2 s-1 or kg m-1 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 [H ~> m or kg m-2].
    netmassout,   &  ! The mass flux (if non-Boussinesq) or volume flux (if Boussinesq)
                     ! over a time step from evaporating fresh water [H ~> m or kg m-2]
    net_heat, & !   The net heating at the surface over a time step [degC H ~> degC m or degC kg m-2].
                ! Any penetrating shortwave radiation is not included in Net_heat.
    net_salt, & ! The surface salt flux into the ocean over a time step, ppt H.
    idecay_len_tke, &  ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
    p_ref, &    !   Reference pressure for the potential density governing mixed
                ! layer dynamics, almost always 0 (or 1e5) [R L2 T-2 ~> Pa].
    p_ref_cv, & !   Reference pressure for the potential density which defines
                ! the coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].
    dr0_dt, &   !   Partial derivative of the mixed layer potential density with
                ! temperature [R degC-1 ~> kg m-3 degC-1].
    drcv_dt, &  !   Partial derivative of the coordinate variable potential
                ! density in the mixed layer with temperature [R degC-1 ~> kg m-3 degC-1].
    dr0_ds, &   !   Partial derivative of the mixed layer potential density with
                ! salinity [R ppt-1 ~> kg m-3 ppt-1].
    drcv_ds, &  !   Partial derivative of the coordinate variable potential
                ! density in the mixed layer with salinity [R ppt-1 ~> kg m-3 ppt-1].
    tke_river   ! The source of turbulent kinetic energy available for mixing
                ! at rivermouths [Z L2 T-3 ~> m3 s-3].

  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 [degC H ~> degC m or degC kg m-2].
  real, dimension(max(CS%nsw,1),SZI_(G),SZK_(GV)) :: &
    opacity_band ! The opacity in each band [H-1 ~> m-1 or m2 kg-1]. The indicies are band, i, k.

  real :: cmke(2,szi_(g)) ! Coefficients of HpE and HpE^2 used in calculating the
                          ! denominator of MKE_rate; the two elements have differing
                          ! units of [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
  real :: irho0         ! 1.0 / rho_0 [R-1 ~> m3 kg-1]
  real :: inkml, inkmlm1!  1.0 / REAL(nkml) and  1.0 / REAL(nkml-1)
  real :: ih            !   The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
  real :: idt_diag      !   The inverse of the timestep used for diagnostics [T-1 ~> s-1].
  real :: rmixconst

  real, dimension(SZI_(G)) :: &
    dke_fc, &   !   The change in mean kinetic energy due to free convection
                ! [Z L2 T-2 ~> m3 s-2].
    h_ca        !   The depth to which convective adjustment has gone [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(GV)) :: &
    dke_ca, &   !   The change in mean kinetic energy due to convective
                ! adjustment [Z L2 T-2 ~> m3 s-2].
    ctke        !   The turbulent kinetic energy source due to convective
                ! adjustment [Z L2 T-2 ~> 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 [Z ~> m].
    hsfc_used, & ! The thickness of the surface region after buffer layer
                ! detrainment [Z ~> 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 [Z ~> m].
    h_sum, &    ! The total thickness of the water column [H ~> m or kg m-2].
    hmbl_prev   ! The previous thickness of the mixed and buffer layers [H ~> m or kg m-2].
  real, dimension(SZI_(G)) :: &
    hsfc, &     !   The thickness of the surface region (mixed and buffer
                ! layers before detrainment in to the interior [H ~> m or kg m-2].
    max_bl_det  !   If non-negative, the maximum amount of entrainment from
                ! the buffer layers that will be allowed this time step [H ~> m or kg m-2].
  real :: dhsfc, dhd ! Local copies of nondimensional parameters.
  real :: h_nbr ! A minimum thickness based on neighboring thicknesses [H ~> m or kg m-2].

  real :: absf_x_h  ! The absolute value of f times the mixed layer thickness [Z T-1 ~> m s-1].
  real :: ku_star   ! Ustar times the Von Karmen constant [Z T-1 ~> m s-1].
  real :: dt__diag  ! A recaled copy of dt_diag (if present) or dt [T ~> s].
  logical :: write_diags  ! If true, write out diagnostics with this step.
  logical :: reset_diags  ! If true, zero out the accumulated diagnostics.
  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
  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_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 do default(shared)
    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 parallel do default(shared)
    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

    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
    if (cs%TKE_diagnostics) then
      !$OMP parallel do default(shared)
      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 parallel do default(shared)
      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 parallel do default(shared)
      do j=js,je ; do i=is,ie
        cs%diag_PE_detrain2(i,j) = 0.0
      enddo ; enddo
    endif
  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
  eosdom(:) = eos_domain(g%HI)

  !$OMP parallel default(shared) firstprivate(dKE_CA,cTKE,h_CA,max_BL_det,p_ref,p_ref_cv) &
  !$OMP                 private(h,u,v,h_orig,eps,T,S,opacity_band,d_ea,d_eb,R0,Rcv,ksort, &
  !$OMP                         dR0_dT,dR0_dS,dRcv_dT,dRcv_dS,htot,Ttot,Stot,TKE,Conv_en, &
  !$OMP                         RmixConst,TKE_river,Pen_SW_bnd,netMassInOut,NetMassOut,   &
  !$OMP                         Net_heat,Net_salt,uhtot,vhtot,R0_tot,Rcv_tot,dKE_FC,      &
  !$OMP                         Idecay_len_TKE,cMKE,Hsfc,dHsfc,dHD,H_nbr,kU_Star,         &
  !$OMP                         absf_x_H,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_H
      t(i,k) = tv%T(i,j,k) ; s(i,k) = tv%S(i,j,k)
    enddo ; enddo
    if (nsw>0) call extract_optics_slice(optics, j, g, gv, opacity=opacity_band, opacity_scale=gv%H_to_m)

    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 [R L2 T-2 ~> Pa]
    if (associated(tv%p_surf)) then
      do i=is,ie ; p_ref(i) = tv%p_surf(i,j) ; enddo
    else
      do i=is,ie ; p_ref(i) = 0.0 ; enddo
    endif
    do k=1,cs%nkml ; do i=is,ie
      p_ref(i) = p_ref(i) + 0.5*(gv%H_to_RZ*gv%g_Earth)*h(i,k)
    enddo ; enddo
    call calculate_density_derivs(t(:,1), s(:,1), p_ref, dr0_dt, dr0_ds, tv%eqn_of_state, eosdom)
    call calculate_density_derivs(t(:,1), s(:,1), p_ref_cv, drcv_dt, drcv_ds, tv%eqn_of_state, eosdom)
    do k=1,nz
      call calculate_density(t(:,k), s(:,k), p_ref, r0(:,k), tv%eqn_of_state, eosdom)
      call calculate_density(t(:,k), s(:,k), p_ref_cv, rcv(:,k), tv%eqn_of_state, eosdom)
    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, us, 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, us, 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 ~> m or kg m-2] added/removed via surface fluxes
    ! netMassOut   = water [H ~> m or kg m-2] removed via evaporating surface fluxes
    ! net_heat     = heat via surface fluxes [degC H ~> degC m or degC kg m-2]
    ! net_salt     = salt via surface fluxes [ppt H ~> dppt m or gSalt m-2]
    ! Pen_SW_bnd   = components to penetrative shortwave radiation
    call extractfluxes1d(g, gv, us, 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, us, 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 effecti1336vely 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, us, 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, us, cs)

    call absorbremainingsw(g, gv, us, h(:,1:), opacity_band, nsw, optics, 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  ! Rescale the diagnostic.
    enddo ; endif
    if (associated(hml)) then ; do i=is,ie
      hml(i,j) = g%mask2dT(i,j) * (h(i,0) * gv%H_to_Z) ! Rescale the diagnostic for output.
    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_Z * 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) = gv%H_to_Z * hsfc(i)
      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, us, 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, us, 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) = gv%H_to_Z * h(i,0) ; enddo
      do k=cs%nkml+1,nkmb ; do i=is,ie
        hsfc_used(i,j) = hsfc_used(i,j) + gv%H_to_Z * h(i,k)
      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 * gv%H_to_Z * 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) = gv%H_to_Z * 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) + gv%H_to_Z * 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) + gv%H_to_Z * abs(h_3d(i,j,nz) - (h_orig(i,nz) + &
          ((eaml(i,nz) - ebml(i,nz-1)) + ebml(i,nz))))
      enddo
    endif

  enddo ! j loop
  !$OMP end parallel

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


  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, US, 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 [H ~> m or kg m-2].
                                                           !! 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 [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: v    !< Zonal velocities interpolated to h
                                                           !! points [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: T    !< Layer temperatures [degC].
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: S    !< Layer salinities [ppt].
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: R0   !< Potential density referenced to
                                                           !! surface pressure [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: Rcv  !< The coordinate defining potential
                                                           !! density [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(GV)), intent(inout) :: d_eb !< The downward increase across a layer
                                                           !! in the entrainment from below [H ~> m or kg m-2].
                                                           !! Positive values go with mass gain by
                                                           !! a layer.
  real, dimension(SZI_(G),SZK_(GV)), intent(in)    :: eps  !< The negligibly small amount of water
                                                           !! that will be left in each layer [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(GV)), intent(out)   :: dKE_CA !< The vertically integrated change in
                                                           !! kinetic energy due to convective
                                                           !! adjustment [Z L2 T-2 ~> m3 s-2].
  real, dimension(SZI_(G),SZK_(GV)), intent(out)   :: cTKE !< The buoyant turbulent kinetic energy
                                                           !! source due to convective adjustment
                                                           !! [Z L2 T-2 ~> m3 s-2].
  integer,                           intent(in)    :: j    !< The j-index to work on.
  type(unit_scale_type),             intent(in)    :: US   !< A dimensional unit scaling type
  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.

  ! Local variables
  real, dimension(SZI_(G)) :: &
    htot, &     !   The total depth of the layers being considered for
                ! entrainment [H ~> m or kg m-2].
    r0_tot, &   !   The integrated potential density referenced to the surface
                ! of the layers which are fully entrained [H R ~> kg m-2 or kg2 m-5].
    rcv_tot, &  !   The integrated coordinate value potential density of the
                ! layers that are fully entrained [H R ~> kg m-2 or kg2 m-5].
    ttot, &     !   The integrated temperature of layers which are fully
                ! entrained [degC H ~> degC m or degC kg m-2].
    stot, &     !   The integrated salt of layers which are fully entrained
                ! [H ppt ~> m ppt or ppt kg m-2].
    uhtot, &    !   The depth integrated zonal and meridional velocities in
    vhtot, &    ! the mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
    ke_orig, &  !   The total mean kinetic energy in the mixed layer before
                ! convection, [H L2 T-2 ~> H m2 s-2].
    h_orig_k1   !   The depth of layer k1 before convective adjustment [H ~> m or kg m-2].
  real :: h_ent !   The thickness from a layer that is entrained [H ~> m or kg m-2].
  real :: Ih    !   The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
  real :: g_H2_2Rho0  !   Half the gravitational acceleration times the square of
                      ! the conversion from H to Z divided by the mean density,
                      ! in [L2 Z T-3 H-2 R-1 ~> m4 s-3 kg-1 or m10 s-3 kg-3].
  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_Z**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_Z * (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, US, 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 [H ~> m or kg m-2].
                                                   !! 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
                                                   !! layer in the entrainment from below [H ~> m or kg m-2].
                                                   !! Positive values go with mass gain by a layer.
  real, dimension(SZI_(G)), intent(out)   :: htot  !< The accumulated mixed layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G)), intent(out)   :: Ttot  !< The depth integrated mixed layer temperature
                                                   !! [degC H ~> degC m or degC kg m-2].
  real, dimension(SZI_(G)), intent(out)   :: Stot  !< The depth integrated mixed layer salinity
                                                   !! [ppt H ~> ppt m or ppt kg m-2].
  real, dimension(SZI_(G)), intent(out)   :: uhtot !< The depth integrated mixed layer zonal
                                                   !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
  real, dimension(SZI_(G)), intent(out)   :: vhtot !< The integrated mixed layer meridional
                                                   !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
  real, dimension(SZI_(G)), intent(out)   :: R0_tot !< The integrated mixed layer potential density referenced
                                                   !! to 0 pressure [H R ~> kg m-2 or kg2 m-5].
  real, dimension(SZI_(G)), intent(out)   :: Rcv_tot !< The integrated mixed layer coordinate
                                                   !! variable potential density [H R ~> kg m-2 or kg2 m-5].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: u     !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: v     !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: T     !< Layer temperatures [degC].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: S     !< Layer salinities [ppt].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: R0    !< Potential density referenced to
                                                   !! surface pressure [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: Rcv   !< The coordinate defining potential
                                                   !! density [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: eps   !< The negligibly small amount of water
                                                   !! that will be left in each layer [H ~> m or kg m-2].
  real, dimension(SZI_(G)), intent(in)    :: dR0_dT  !< The partial derivative of R0 with respect to
                                                   !! temperature [R degC-1 ~> kg m-3 degC-1].
  real, dimension(SZI_(G)), intent(in)    :: dRcv_dT !< The partial derivative of Rcv with respect to
                                                   !! temperature [R degC-1 ~> kg m-3 degC-1].
  real, dimension(SZI_(G)), intent(in)    :: dR0_dS  !< The partial derivative of R0 with respect to
                                                   !! salinity [R ppt-1 ~> kg m-3 ppt-1].
  real, dimension(SZI_(G)), intent(in)    :: dRcv_dS !< The partial derivative of Rcv with respect to
                                                   !! salinity [R ppt-1 ~> kg m-3 ppt-1].
  real, dimension(SZI_(G)), intent(in)    :: netMassInOut !< The net mass flux (if non-Boussinesq)
                                                   !! or volume flux (if Boussinesq) into the ocean
                                                   !! within a time step [H ~> m or kg m-2]. (I.e. P+R-E.)
  real, dimension(SZI_(G)), intent(in)    :: netMassOut !< The mass or volume flux out of the ocean
                                                   !! within a time step [H ~> m or kg m-2].
  real, dimension(SZI_(G)), intent(in)    :: Net_heat !< The net heating at the surface over a time
                                                   !! step [degC H ~> degC m or degC kg m-2].  Any penetrating
                                                   !! shortwave radiation is not included in Net_heat.
  real, dimension(SZI_(G)), intent(in)    :: Net_salt !< The net surface salt flux into the ocean
                                                   !! over a time step [ppt H ~> ppt m or ppt kg m-2].
  integer,                  intent(in)    :: nsw   !< The number of bands of penetrating
                                                   !! shortwave radiation.
  real, dimension(max(nsw,1),SZI_(G)), intent(inout) :: Pen_SW_bnd !< The penetrating shortwave
                                                   !! heating at the sea surface in each penetrating
                                                   !! band [degC H ~> degC m or degC kg m-2].
  real, dimension(max(nsw,1),SZI_(G),SZK_(GV)), intent(in) :: opacity_band !< The opacity in each band of
                                                   !! penetrating shortwave radiation [H-1 ~> m-1 or m2 kg-1].
  real, dimension(SZI_(G)), intent(out)   :: Conv_En !< The buoyant turbulent kinetic energy source
                                                   !! due to free convection [Z L2 T-2 ~> m3 s-2].
  real, dimension(SZI_(G)), intent(out)   :: dKE_FC !< The vertically integrated change in kinetic
                                                   !! energy due to free convection [Z L2 T-2 ~> 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(unit_scale_type),    intent(in)    :: US    !< A dimensional unit scaling type
  type(bulkmixedlayer_CS),  pointer       :: CS    !< The control structure for this module.
  type(thermo_var_ptrs),    intent(inout) :: tv    !< A structure containing pointers to any
                                                   !! available thermodynamic fields. Absent
                                                   !! fields have NULL ptrs.
  type(forcing),            intent(inout) :: fluxes  !< A structure containing pointers to any
                                                   !! possible forcing fields.  Unused fields
                                                   !! have NULL ptrs.
  real,                     intent(in)    :: dt    !< Time increment [T ~> s].
  logical,                  intent(in)    :: aggregate_FW_forcing !< If true, the net incoming and
                                                   !! outgoing surface freshwater fluxes are
                                                   !! combined before being applied, instead of
                                                   !! being applied separately.

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

  ! Local variables
  real, dimension(SZI_(G)) :: &
    massOutRem, &      !   Evaporation that remains to be supplied [H ~> m or kg m-2].
    netMassIn          ! mass entering through ocean surface [H ~> m or kg m-2]
  real :: SW_trans     !   The fraction of shortwave radiation
                       ! that is not absorbed in a layer [nondim].
  real :: Pen_absorbed !   The amount of penetrative shortwave radiation
                       ! that is absorbed in a layer [degC H ~> degC m or degC kg m-2].
  real :: h_avail      !   The thickness in a layer available for
                       ! entrainment [H ~> m or kg m-2].
  real :: h_ent        !   The thickness from a layer that is entrained [H ~> m or kg m-2].
  real :: T_precip     !   The temperature of the precipitation [degC].
  real :: C1_3, C1_6   !  1/3 and 1/6.
  real :: En_fn, Frac, x1 !  Nondimensional temporary variables.
  real :: dr, dr0      ! Temporary variables [R H ~> kg m-2 or kg2 m-5].
  real :: dr_ent, dr_comp ! Temporary variables [R H ~> kg m-2 or kg2 m-5].
  real :: dr_dh        ! The partial derivative of dr_ent with h_ent [R ~> kg m-3].
  real :: h_min, h_max !   The minimum, maximum, and previous estimates for
  real :: h_prev       ! h_ent [H ~> m or kg m-2].
  real :: h_evap       !   The thickness that is evaporated [H ~> m or kg m-2].
  real :: dh_Newt      !   The Newton's method estimate of the change in
                       ! h_ent between iterations [H ~> m or kg m-2].
  real :: g_H2_2Rho0   !   Half the gravitational acceleration times the square of
                       ! the conversion from H to Z divided by the mean density,
                       ! [L2 Z T-3 H-2 R-1 ~> m4 s-3 kg-1 or m10 s-3 kg-3].
  real :: Angstrom     !   The minimum layer thickness [H ~> m or kg m-2].
  real :: opacity      !   The opacity converted to inverse thickness units [H-1 ~> m-1 or m2 kg-1]
  real :: sum_Pen_En   !   The potential energy change due to penetrating
                       ! shortwave radiation, integrated over a layer
                       ! [H R ~> kg m-2 or kg2 m-5].
  real :: Idt          ! 1.0/dt [T-1 ~> s-1]
  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 R H-1 ~> kg m-4 or m-1].
    r_SW_top           ! Temporary variables [H R ~> kg m-2 or kg2 m-5].

  angstrom = gv%Angstrom_H
  c1_3 = 1.0/3.0 ; c1_6 = 1.0/6.0
  g_h2_2rho0 = (gv%g_Earth * gv%H_to_Z**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_RZ * fluxes%C_p * idt
    if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
                         t_precip * netmassin(i) * gv%H_to_RZ
  else  ! This is a massless column, but zero out the summed variables anyway for safety.
    htot(i) = 0.0 ; ttot(i) = 0.0 ; stot(i) = 0.0 ; r0_tot(i) = 0.0 ; rcv_tot = 0.0
    uhtot(i) = 0.0 ; vhtot(i) = 0.0 ; conv_en(i) = 0.0 ; dke_fc(i) = 0.0
  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_RZ*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_RZ * fluxes%C_p * idt
        if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) - &
                                      t(i,k)*h_evap*gv%H_to_RZ

      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)
              elseif (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_Z*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
          if (cs%convect_mom_bug) then
            uhtot(i) = u(i,k)*h_ent ; vhtot(i) = v(i,k)*h_ent
          else
            uhtot(i) = uhtot(i) + h_ent*u(i,k) ; vhtot(i) = vhtot(i) + h_ent*v(i,k)
          endif
        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, US, 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.
  type(unit_scale_type),      intent(in)    :: US      !< A dimensional unit scaling type
  real, dimension(SZI_(G)),   intent(in)    :: htot    !< The accumulated mixed layer thickness
                                                       !! [H ~> m or kg m-2]
  real, dimension(SZI_(G)),   intent(in)    :: h_CA    !< The mixed layer depth after convective
                                                       !! adjustment [H ~> m or kg m-2].
  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 [Z L2 T-2 ~> m3 s-2].
  real, dimension(SZI_(G)),   intent(in)    :: dKE_FC  !< The vertically integrated change in
                                                       !! kinetic energy due to free convection
                                                       !! [Z L2 T-2 ~> m3 s-2].
  real, dimension(SZI_(G),SZK_(GV)), &
                              intent(in)    :: cTKE    !< The buoyant turbulent kinetic energy
                                                       !! source due to convective adjustment
                                                       !! [Z L2 T-2 ~> m3 s-2].
  real, dimension(SZI_(G),SZK_(GV)), &
                              intent(in)    :: dKE_CA  !< The vertically integrated change in
                                                       !! kinetic energy due to convective
                                                       !! adjustment [Z L2 T-2 ~> m3 s-2].
  real, dimension(SZI_(G)),   intent(out)   :: TKE     !< The turbulent kinetic energy available for
                                                       !! mixing over a time step [Z m2 T-2 ~> m3 s-2].
  real, dimension(SZI_(G)),   intent(out)   :: Idecay_len_TKE !< The inverse of the vertical decay
                                                       !! scale for TKE [H-1 ~> m-1 or m2 kg-1].
  real, dimension(SZI_(G)),   intent(in)    :: TKE_river !< The source of turbulent kinetic energy
                                                       !! available for driving mixing at river mouths
                                                       !! [Z L2 T-3 ~> m3 s-3].
  real, dimension(2,SZI_(G)), intent(out)   :: cMKE    !< Coefficients of HpE and HpE^2 in
                                                       !! calculating the denominator of MKE_rate,
                                                       !! [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
  real,                       intent(in)    :: dt      !< The time step [T ~> s].
  real,                       intent(in)    :: Idt_diag !< The inverse of the accumulated diagnostic
                                                       !! time interval [T-1 ~> s-1].
  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.

  ! Local variables
  real :: dKE_conv  ! The change in mean kinetic energy due to all convection [Z L2 T-2 ~> m3 s-2].
  real :: nstar_FC  ! The effective efficiency with which the energy released by
                    ! free convection is converted to TKE, often ~0.2 [nondim].
  real :: nstar_CA  ! The effective efficiency with which the energy released by
                    ! convective adjustment is converted to TKE, often ~0.2 [nondim].
  real :: TKE_CA    ! The potential energy released by convective adjustment if
                    ! that release is positive [Z L2 T-2 ~> m3 s-2].
  real :: MKE_rate_CA ! MKE_rate for convective adjustment [nondim], 0 to 1.
  real :: MKE_rate_FC ! MKE_rate for free convection [nondim], 0 to 1.
  real :: totEn_Z   ! The total potential energy released by convection, [Z3 T-2 ~> m3 s-2].
  real :: Ih        ! The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
  real :: exp_kh    ! The nondimensional decay of TKE across a layer [nondim].
  real :: absf      ! The absolute value of f averaged to thickness points [T-1 ~> s-1].
  real :: U_star    ! The friction velocity [Z T-1 ~> m s-1].
  real :: absf_Ustar  ! The absolute value of f divided by U_star [Z-1 ~> m-1].
  real :: wind_TKE_src ! The surface wind source of TKE [Z L2 T-3 ~> m3 s-3].
  real :: diag_wt   ! The ratio of the current timestep to the diagnostic
                    ! timestep (which may include 2 calls) [nondim].
  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_Z

!    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_Z/(3.0*0.41*u_star*dt)
    cmke(1,i) = 4.0 * ih ; cmke(2,i) = (absf_ustar*gv%H_to_Z) * 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_z = us%L_to_Z**2 * (conv_en(i) + tke_ca)

      if (toten_z > 0.0) then
        nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
                        sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_Z))**3 * toten_z))
      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_z = us%L_to_Z**2 * (conv_en(i) + tke_ca * (htot(i) / h_ca(i)) )
        nstar_fc = cs%nstar * toten_z / (toten_z + 0.2 * &
                        sqrt(0.5 * dt * (absf*(htot(i)*gv%H_to_Z))**3 * toten_z))
      else
        nstar_fc = cs%nstar
      endif

      toten_z = us%L_to_Z**2 * (conv_en(i) + tke_ca)
      if (tke_ca > 0.0) then
        nstar_ca = cs%nstar * toten_z / (toten_z + 0.2 * &
                        sqrt(0.5 * dt * (absf*(h_ca(i)*gv%H_to_Z))**3 * toten_z))
      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)*((us%Z_to_L**2*(u_star*u_star*u_star))*exp_kh) + &
             (exp_kh * dke_conv + nstar_fc*conv_en(i) + nstar_ca * tke_ca)

    if (cs%do_rivermix) then ! Add additional TKE at river mouths
      tke(i) = tke(i) + tke_river(i)*dt*exp_kh
    endif

    if (cs%TKE_diagnostics) then
      wind_tke_src = cs%mstar*(us%Z_to_L**2*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, US, 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.
  type(unit_scale_type),    intent(in)    :: US    !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(inout) :: h     !< Layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(inout) :: d_eb  !< The downward increase across a layer in the
                                                   !! layer in the entrainment from below [H ~> m or kg m-2].
                                                   !! Positive values go with mass gain by a layer.
  real, dimension(SZI_(G)), intent(inout) :: htot  !< The accumlated mixed layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G)), intent(inout) :: Ttot  !< The depth integrated mixed layer temperature
                                                   !! [degC H ~> degC m or degC kg m-2].
  real, dimension(SZI_(G)), intent(inout) :: Stot  !< The depth integrated mixed layer salinity
                                                   !! [ppt H ~> ppt m or ppt kg m-2].
  real, dimension(SZI_(G)), intent(inout) :: uhtot !< The depth integrated mixed layer zonal
                                                   !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
  real, dimension(SZI_(G)), intent(inout) :: vhtot !< The integrated mixed layer meridional
                                                   !! velocity [H L T-1 ~> m2 s-1 or kg m-1 s-1].
  real, dimension(SZI_(G)), intent(inout) :: R0_tot !< The integrated mixed layer potential density
                                                   !! referenced to 0 pressure [H R ~> kg m-2 or kg2 m-5].
  real, dimension(SZI_(G)), intent(inout) :: Rcv_tot !< The integrated mixed layer coordinate variable
                                                   !! potential density [H R ~> kg m-2 or kg2 m-5].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: u     !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: v     !< Zonal velocities interpolated to h points [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: T     !< Layer temperatures [degC].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: S     !< Layer salinities [ppt].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: R0    !< Potential density referenced to
                                                   !! surface pressure [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: Rcv   !< The coordinate defining potential
                                                   !! density [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(GV)), &
                            intent(in)    :: eps   !< The negligibly small amount of water
                                                   !! that will be left in each layer [H ~> m or kg m-2].
  real, dimension(SZI_(G)), intent(in)    :: dR0_dT  !< The partial derivative of R0 with respect to
                                                   !! temperature [R degC-1 ~> kg m-3 degC-1].
  real, dimension(SZI_(G)), intent(in)    :: dRcv_dT !< The partial derivative of Rcv with respect to
                                                   !! temperature [R degC-1 ~> kg m-3 degC-1].
  real, dimension(2,SZI_(G)), intent(in)  :: cMKE  !< Coefficients of HpE and HpE^2 used in calculating the
                                                   !! denominator of MKE_rate; the two elements have differing
                                                   !! units of [H-1 ~> m-1 or m2 kg-1] and [H-2 ~> m-2 or m4 kg-2].
  real,                     intent(in)    :: Idt_diag !< The inverse of the accumulated diagnostic
                                                   !! time interval [T-1 ~> s-1].
  integer,                  intent(in)    :: nsw   !< The number of bands of penetrating
                                                   !! shortwave radiation.
  real, dimension(max(nsw,1),SZI_(G)), intent(inout) :: Pen_SW_bnd !< The penetrating shortwave
                                                   !! heating at the sea surface in each penetrating
                                                   !! band [degC H ~> degC m or degC kg m-2].
  real, dimension(max(nsw,1),SZI_(G),SZK_(GV)), intent(in) :: opacity_band !< The opacity in each band of
                                                   !! penetrating shortwave radiation [H-1 ~> m-1 or m2 kg-1].
  real, dimension(SZI_(G)), intent(inout) :: TKE   !< The turbulent kinetic energy
                                                   !! available for mixing over a time
                                                   !! step [Z m2 T-2 ~> m3 s-2].
  real, dimension(SZI_(G)), intent(inout) :: Idecay_len_TKE !< The vertical TKE decay rate [H-1 ~> m-1 or m2 kg-1].
  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.

  ! Local variables
  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 [degC H ~> degC m or degC kg m-2].
  real :: h_avail   ! The thickness in a layer available for entrainment [H ~> m or kg m-2].
  real :: h_ent     ! The thickness from a layer that is entrained [H ~> m or kg m-2].
  real :: h_min, h_max ! Limits on the solution for h_ent [H ~> m or kg m-2].
  real :: dh_Newt      !   The Newton's method estimate of the change in
                       ! h_ent between iterations [H ~> m or kg m-2].
  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 ~> m or kg m-2].
  real :: g_H_2Rho0   !   Half the gravitational acceleration times the
                      ! conversion from H to m divided by the mean density,
                      ! in [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
  real :: TKE_full_ent  ! The TKE remaining if a layer is fully entrained
                        ! [Z L2 T-2 ~> m3 s-2].
  real :: dRL       ! Work required to mix water from the next layer
                    ! across the mixed layer [L2 T-2 ~> L2 s-2].
  real :: Pen_En_Contrib  ! Penetrating SW contributions to the changes in
                          ! TKE, divided by layer thickness in m [L2 T2 ~> m2 s-2].
  real :: Cpen1     ! A temporary variable [L2 T-2 ~> m2 s-2].
  real :: dMKE      ! A temporary variable related to the release of mean
                    ! kinetic energy [H Z L2 T-2 ~> m4 s-2 or kg m s-2]
  real :: TKE_ent   ! The TKE that remains if h_ent were entrained [Z L2 T-2 ~> m3 s-2].
  real :: TKE_ent1  ! The TKE that would remain, without considering the
                    ! release of mean kinetic energy [Z L2 T-2 ~> m3 s-2].
  real :: dTKE_dh   ! The partial derivative of TKE with h_ent [Z L2 T-2 H-1 ~> m2 s-2 or m5 s-2 kg-1].
  real :: Pen_dTKE_dh_Contrib ! The penetrating shortwave contribution to
                    ! dTKE_dh [L2 T-2 ~> m2 s-2].
  real :: EF4_val   ! The result of EF4() (see later) [H-1 ~> m-1 or m2 kg-1].
  real :: h_neglect ! A thickness that is so small it is usually lost
                    ! in roundoff and can be neglected [H ~> m or kg m-2].
  real :: dEF4_dh   ! The partial derivative of EF4 with h [H-2 ~> m-2 or m4 kg-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 [H-1 ~> m-1 or m2 kg-1].
  real :: Hmix_min  ! The minimum mixed layer depth [H ~> m or kg m-2].
  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
  g_h_2rho0 = (gv%g_Earth * gv%H_to_Z) / (2.0 * gv%Rho0)
  hmix_min = cs%Hmix_min
  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 = (gv%H_to_Z * 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*gv%H_to_Z)*(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*gv%H_to_Z)*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*(gv%H_to_Z*h_ent)*drl
            cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
                idt_diag*(gv%H_to_Z*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
          !### The minimum TKE value in this line may be problematically small.
          if (tke(i) <= 0.0) tke(i) = 1.0e-150*us%m_to_Z*us%m_s_to_L_T**2
        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
                cpen1 = g_h_2rho0*dr0_dt(i)*pen_sw_bnd(n,i)
                pen_en_contrib = pen_en_contrib + cpen1*(pen_en1 - f1_kh)
                pen_dtke_dh_contrib = pen_dtke_dh_contrib + &
                           cpen1*((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*gv%H_to_Z)*(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*gv%H_to_Z) + &
                         pen_dtke_dh_contrib*gv%H_to_Z) + 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_H) 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*gv%H_to_Z)*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*gv%H_to_Z)*drl
            cs%diag_TKE_pen_SW(i,j) = cs%diag_TKE_pen_SW(i,j) - &
                idt_diag*(h_ent*gv%H_to_Z)*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 [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(GV)),    intent(in)  :: R0    !< The potential density used to sort
                                                             !! the layers [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(GV)),    intent(in)  :: eps   !< The (small) thickness that must
                                                             !! remain in each layer [H ~> m or kg m-2].
  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.

  ! Local variables
  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 [H ~> m or kg m-2].
                                                                 !! Layer 0 is the new mixed layer.
  real, dimension(SZI_(G),SZK0_(GV)),   intent(inout) :: T       !< Layer temperatures [degC].
  real, dimension(SZI_(G),SZK0_(GV)),   intent(inout) :: S       !< Layer salinities [ppt].
  real, dimension(SZI_(G),SZK0_(GV)),   intent(inout) :: R0      !< Potential density referenced to
                                                                 !! surface pressure [R ~> kg m-3].
  real, dimension(SZI_(G),SZK0_(GV)),   intent(inout) :: Rcv     !< The coordinate defining
                                                                 !! potential density [R ~> kg m-3].
  real, dimension(SZK_(GV)),            intent(in)    :: RcvTgt  !< The target value of Rcv for each
                                                                 !! layer [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(GV)),    intent(inout) :: eps     !< The (small) thickness that must
                                                                 !! remain in each layer [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(GV)),    intent(inout) :: d_ea    !< The upward increase across a
                                                                 !! layer in the entrainment from
                                                                 !! above [H ~> m or kg m-2].
                                                                 !!  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 [H ~> m or kg m-2]. 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 [R degC-1 ~> kg m-3 degC-1].
  real, dimension(SZI_(G)),             intent(in)    :: dR0_dS  !< The partial derivative of
                                                                 !! cpotential density referenced
                                                                 !! to the surface with salinity,
                                                                 !! [R ppt-1 ~> kg m-3 ppt-1].
  real, dimension(SZI_(G)),             intent(in)    :: dRcv_dT !< The partial derivative of
                                                                 !! coordinate defining potential
                                                                 !! density with potential
                                                                 !! temperature [R degC-1 ~> kg m-3 degC-1].
  real, dimension(SZI_(G)),             intent(in)    :: dRcv_dS !< The partial derivative of
                                                                 !! coordinate defining potential
                                                                 !! density with salinity,
                                                                 !! [R ppt-1 ~> kg m-3 ppt-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

  ! Local variables
  real    :: h_move, h_tgt_old, I_hnew
  real    :: dT_dS_wt2, dT_dR, dS_dR, I_denom
  real    :: Rcv_int
  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

  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) < cs%BL_split_rho_tol*(rcvtgt(k) - rcvtgt(k-1))) &
            leave_in_layer = .true.
        elseif (k > nkmb) then
          if (rcv(i,k)-rcvtgt(k) < cs%BL_split_rho_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, US, 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 [H ~> m or kg m-2].
                                                            !!  Layer 0 is the new mixed layer.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T    !< Potential temperature [degC].
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S    !< Salinity [ppt].
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0   !< Potential density referenced to
                                                            !! surface pressure [R ~> kg m-3].
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv  !< The coordinate defining potential
                                                            !! density [R ~> kg m-3].
  real, dimension(SZK_(GV)),          intent(in)    :: RcvTgt  !< The target value of Rcv for each
                                                            !! layer [R ~> kg m-3].
  real,                               intent(in)    :: dt   !< Time increment [T ~> s].
  real,                               intent(in)    :: dt_diag !< The diagnostic time step [T ~> s].
  real, dimension(SZI_(G),SZK_(GV)),  intent(inout) :: d_ea !< The upward increase across a layer in
                                                            !! the entrainment from above
                                                            !! [H ~> m or kg m-2]. Positive d_ea
                                                            !! goes with layer thickness increases.
  integer,                            intent(in)    :: j    !< The meridional row to work on.
  type(unit_scale_type),              intent(in)    :: US   !< A dimensional unit scaling type
  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,
                                                            !! [R degC-1 ~> kg m-3 degC-1].
  real, dimension(SZI_(G)),           intent(in)    :: dR0_dS  !< The partial derivative of
                                                            !! cpotential density referenced to the
                                                            !! surface with salinity
                                                            !! [R ppt-1 ~> kg m-3 ppt-1].
  real, dimension(SZI_(G)),           intent(in)    :: dRcv_dT !< The partial derivative of
                                                            !! coordinate defining potential density
                                                            !! with potential temperature,
                                                            !! [R degC-1 ~> kg m-3 degC-1].
  real, dimension(SZI_(G)),           intent(in)    :: dRcv_dS !< The partial derivative of
                                                            !! coordinate defining potential density
                                                            !! with salinity [R ppt-1 ~> kg m-3 ppt-1].
  real, dimension(SZI_(G)),           intent(in)    :: max_BL_det !< If non-negative, the maximum
                                                            !! detrainment permitted from the buffer
                                                            !! layers [H ~> m or kg m-2].

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

  ! Local variables
  real :: h_to_bl                 ! The total thickness detrained to the buffer
                                  ! layers [H ~> m or kg m-2].
  real :: R0_to_bl                ! The depth integrated amount of R0 that is detrained to the
                                  ! buffer layer [H R ~> kg m-2 or kg2 m-5]
  real :: Rcv_to_bl               ! The depth integrated amount of Rcv that is detrained to the
                                  ! buffer layer [H R ~> kg m-2 or kg2 m-5]
  real :: T_to_bl                 ! The depth integrated amount of T that is detrained to the
                                  ! buffer layer [degC H ~> degC m or degC kg m-2]
  real :: S_to_bl                 ! The depth integrated amount of S that is detrained to the
                                  ! buffer layer [ppt H ~> ppt m or ppt kg m-2]
  real :: h_min_bl                ! The minimum buffer layer thickness [H ~> m or kg m-2].

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

  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 [H ~> m or kg m-2].
  real :: stays_min_merge         ! The minimum allowed value of stays_merge [H ~> m or kg m-2].

  real :: dR0_2dz, dRcv_2dz       ! Half the vertical gradients of R0 and Rcv [R H-1 ~> kg m-4 or m-1]
!  real :: dT_2dz, dS_2dz         ! Half the vertical gradients of T and S, in degC H-1, and ppt H-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 [H2 ~> m2 or kg2 m-4].
  real :: dPE_det, dPE_merge      ! The energy required to mix the detrained water
                                  ! into the buffer layer or the merge the two
                                  ! buffer layers [R H2 L2 Z-1 T-2 ~> J m-2 or J kg2 m-8].

  real :: h_from_ml               ! The amount of additional water that must be
                                  ! drawn from the mixed layer [H ~> m or kg m-2].
  real :: h_det_h2                ! The amount of detrained water and mixed layer
                                  ! water that will go directly into the lower
                                  ! buffer layer [H ~> m or kg m-2].
  real :: h_det_to_h2, h_ml_to_h2 ! All of the variables hA_to_hB are the thickness fluxes
  real :: h_det_to_h1, h_ml_to_h1 ! from one layer to another [H ~> m or kg m-2],
  real :: h1_to_h2, h1_to_k0      ! 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 [R ~> kg m-3], T [degC], and S [ppt].
  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 [R ~> 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 [R ~> kg m-3].
  real :: dSpice_2dz              ! The vertical gradient of spiciness used for
                                  ! advection [R H-1 ~> kg m-4 or m-1].

  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 :: dPE_time_ratio          ! Larger of 1 and the detrainment timescale over dt [nondim].
  real :: dT_dS_gauge, dS_dT_gauge ! The relative scales of temperature and
                                  ! salinity changes in defining spiciness, in
                                  ! [degC ppt-1] and [ppt degC-1].
  real :: I_denom                 ! A work variable with units of [ppt2 R-2 ~> ppt2 m6 kg-2].

  real :: g_2                     ! 1/2 g_Earth [L2 Z-1 T-2 ~> m s-2].
  real :: Rho0xG                  ! Rho0 times G_Earth [R L2 Z-1 T-2 ~> kg m-2 s-2].
  real :: I2Rho0                  ! 1 / (2 Rho0) [R-1 ~> m3 kg-1].
  real :: Idt_H2                  ! The square of the conversion from thickness to Z
                                  ! divided by the time step [Z2 H-2 T-1 ~> s-1 or m6 kg-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 [H ~> m or kg m-2].

  real :: s1en                    ! A work variable [H2 L2 kg m-1 T-3 ~> kg m3 s-3 or kg3 m-3 s-3].
  real :: s1, s2, bh0             ! Work variables [H ~> m or kg m-2].
  real :: s3sq                    ! A work variable [H2 ~> m2 or kg2 m-4].
  real :: I_ya, b1                ! Nondimensional work variables.
  real :: Ih, Ihdet, Ih1f, Ih2f   ! Assorted inverse thickness work variables,
  real :: Ihk0, Ihk1, Ih12        ! all in [H-1 ~> m-1 or m2 kg-1].
  real :: dR1, dR2, dR2b, dRk1    ! Assorted density difference work variables,
  real :: dR0, dR21, dRcv         ! all in [R ~> kg m-3].
  real :: dRcv_stays, dRcv_det, dRcv_lim
  real :: Angstrom                ! The minumum layer thickness [H ~> m or kg m-2].

  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_Z**2 / dt_diag
  i2rho0 = 0.5 / (gv%Rho0)
  angstrom = gv%Angstrom_H

  ! This is hard coding of arbitrary and dimensional numbers.
  dt_ds_gauge = cs%dT_dS_wt ; ds_dt_gauge = 1.0 / dt_ds_gauge
  num_events = 10.0

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

  if (dt < cs%BL_detrain_time) then ; dpe_time_ratio = cs%BL_detrain_time / (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(cs%Hbuffer_min, cs%Hbuffer_rel_min*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, US, 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 [H ~> m or kg m-2].
                                                            !! Layer 0 is the new mixed layer.
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: T    !< Potential temperature [degC].
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: S    !< Salinity [ppt].
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: R0   !< Potential density referenced to
                                                            !! surface pressure [R ~> kg m-3].
  real, dimension(SZI_(G),SZK0_(GV)), intent(inout) :: Rcv  !< The coordinate defining potential
                                                            !! density [R ~> kg m-3].
  real, dimension(SZK_(GV)),          intent(in)    :: RcvTgt !< The target value of Rcv for each
                                                            !! layer [R ~> kg m-3].
  real,                               intent(in)    :: dt   !< Time increment [T ~> s].
  real,                               intent(in)    :: dt_diag !< The accumulated time interval for
                                                            !! diagnostics [T ~> s].
  real, dimension(SZI_(G),SZK_(GV)),  intent(inout) :: d_ea !< The upward increase across a layer in
                                                            !! the entrainment from above
                                                            !! [H ~> m or kg m-2]. 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 [H ~> m or kg m-2].
                                                            !! Positive values go with mass gain by
                                                            !! a layer.
  integer,                            intent(in)    :: j    !< The meridional row to work on.
  type(unit_scale_type),              intent(in)    :: US   !< A dimensional unit scaling type
  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
                                                            !! [R degC-1 ~> kg m-3 degC-1].
  real, dimension(SZI_(G)),           intent(in)    :: dRcv_dS    !< The partial derivative of
                                                            !! coordinate defining potential density
                                                            !! with salinity [R ppt-1 ~> kg m-3 ppt-1].
  real, dimension(SZI_(G)),           intent(in)    :: max_BL_det !< If non-negative, the maximum
                                                            !! detrainment permitted from the buffer
                                                            !! layers [H ~> m or kg m-2].

  ! Local variables
  real :: Ih                  ! The inverse of a thickness [H-1 ~> m-1 or m2 kg-1].
  real :: h_ent               ! The thickness from a layer that is
                              ! entrained [H ~> m or kg m-2].
  real :: max_det_rem(SZI_(G)) ! Remaining permitted detrainment [H ~> m or kg m-2].
  real :: detrain(SZI_(G))    ! The thickness of fluid to detrain
                              ! from the mixed layer [H ~> m or kg m-2].
  real :: dT_dR, dS_dR, dRml, dR0_dRcv, dT_dS_wt2
  real :: I_denom             ! A work variable [ppt2 R-2 ~> ppt2 m6 kg-2].
  real :: Sdown, Tdown
  real :: dt_Time             ! The timestep divided by the detrainment timescale [nondim].
  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 [L2 Z T-3 H-2 R-1 ~> m4 s-3 kg-1 or m10 s-3 kg-3].
  real :: g_H2_2dt            ! Half the gravitational acceleration times the square of the
                              ! conversion from H to Z divided by the diagnostic time step
                              ! [L2 Z H-2 T-3 ~> m s-3 or m7 kg-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.")

  dt_time = dt / cs%BL_detrain_time
  g_h2_2rho0dt = (gv%g_Earth * gv%H_to_Z**2) / (2.0 * gv%Rho0 * dt_diag)
  g_h2_2dt = (gv%g_Earth * gv%H_to_Z**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

  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_H) .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_H) .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 initializes the MOM bulk mixed layer module.
subroutine bulkmixedlayer_init(Time, G, GV, US, param_file, diag, CS)
  type(time_type), target, intent(in)    :: time !< The model's clock with the current time.
  type(ocean_grid_type),   intent(in)    :: g    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: gv   !< The ocean's vertical grid structure.
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time
                                                 !! parameters.
  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
                                                 !! output.
  type(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 :: bl_detrain_time_dflt ! The default value for BUFFER_LAY_DETRAIN_TIME [s]
  real :: omega_frac_dflt, ustar_min_dflt, hmix_min_m
  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 "//&
                 "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 "//&
                 "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 "//&
                 "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 "//&
                 "surface fluxes that is available to drive entrainment "//&
                 "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 "//&
                 "by mechanically forced entrainment of the mixed layer "//&
                 "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 "//&
                 "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 "//&
                 "TKE available for mechanical entrainment to the natural "//&
                 "Ekman depth.", units="nondim", default=2.5)
  call get_param(param_file, mdl, "NSTAR2", cs%nstar2, &
                 "The portion of any potential energy released by "//&
                 "convective adjustment that is available to drive "//&
                 "entrainment at the base of mixed layer. By default "//&
                 "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 "//&
                 "kinetic energy is converted to turbulent kinetic "//&
                 "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 "//&
                 "is determined dynamically.", units="m", default=0.0, scale=gv%m_to_H, &
                 unscaled=hmix_min_m)

  call get_param(param_file, mdl, "LIMIT_BUFFER_DETRAIN", cs%limit_det, &
                 "If true, limit the detrainment from the buffer layers "//&
                 "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 "//&
                 "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 "//&
                 "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 "//&
                 "densities, this factor (in deg C / PSU) is combined "//&
                 "with the derivatives of density with T & S to determine "//&
                 "what direction is orthogonal to density contours. It "//&
                 "should be a typical value of (dR/dS) / (dR/dT) in "//&
                 "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 "//&
                 "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, or a negative "//&
                 "value not to apply this limit.", units="nondim", default = -1.0)
  call get_param(param_file, mdl, "BUFFER_LAYER_HMIN_THICK", cs%Hbuffer_min, &
                 "The minimum buffer layer thickness when the mixed layer is very thick.", &
                 units="m", default=5.0, scale=gv%m_to_H)
  call get_param(param_file, mdl, "BUFFER_LAYER_HMIN_REL", cs%Hbuffer_rel_min, &
                 "The minimum buffer layer thickness relative to the combined mixed "//&
                 "land buffer ayer thicknesses when they are thin.", &
                 units="nondim", default=0.1/cs%nkbl)
  bl_detrain_time_dflt = 4.0*3600.0 ; if (cs%nkbl==1) bl_detrain_time_dflt = 86400.0*30.0
  call get_param(param_file, mdl, "BUFFER_LAY_DETRAIN_TIME", cs%BL_detrain_time, &
                 "A timescale that characterizes buffer layer detrainment events.", &
                 units="s", default=bl_detrain_time_dflt, scale=us%s_to_T)
  call get_param(param_file, mdl, "BUFFER_SPLIT_RHO_TOL", cs%BL_split_rho_tol, &
                 "The fractional tolerance for matching layer target densities when splitting "//&
                 "layers to deal with massive interior layers that are lighter than one of the "//&
                 "mixed or buffer layers.", units="nondim", default=0.1)

  call get_param(param_file, mdl, "DEPTH_LIMIT_FLUXES", cs%H_limit_fluxes, &
                 "The surface fluxes are scaled away when the total ocean "//&
                 "depth is less than DEPTH_LIMIT_FLUXES.", &
                 units="m", default=0.1*hmix_min_m, scale=gv%m_to_H)
  call get_param(param_file, mdl, "OMEGA", cs%omega, &
                 "The rotation rate of the earth.", &
                 default=7.2921e-5, units="s-1", scale=us%T_to_s)
  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
                 "If true, use the absolute rotation rate instead of the "//&
                 "vertical component of rotation when setting the decay "//&
                 "scale for turbulence.", default=.false., do_not_log=.true.)
  omega_frac_dflt = 0.0
  if (use_omega) then
    call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
    omega_frac_dflt = 1.0
  endif
  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
                 "When setting the decay scale for turbulence, use this "//&
                 "fraction of the absolute rotation rate blended with the "//&
                 "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
                 units="nondim", default=omega_frac_dflt)
  call get_param(param_file, mdl, "ML_RESORT", cs%ML_resort, &
                 "If true, resort the topmost layers by potential density "//&
                 "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 "//&
                 "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*us%s_to_T*cs%omega*(gv%Angstrom_m + 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 "//&
                 "bulk mixed layer model in setting vertical TKE decay "//&
                 "scales. This must be greater than 0.", units="m s-1", &
                 default=ustar_min_dflt, scale=us%m_to_Z*us%T_to_s)
  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 "//&
                 "chosen to optimally represent the impact of the Ekman "//&
                 "transport on the mixed layer TKE budget.  Otherwise, "//&
                 "the sublayers are distributed uniformly through the "//&
                 "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 "//&
                 "radiation is absorbed is adjusted to match the average "//&
                 "heating depth of an exponential profile by moving some "//&
                 "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 wherever there is "//&
                 "runoff, so that it is mixed down to RIVERMIX_DEPTH, "//&
                 "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 "//&
                 "defined.", units="m", default=0.0, scale=us%m_to_Z)
  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 "//&
                 "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 "//&
                 "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
                 default=.false.)
  call get_param(param_file, mdl, "BULKML_CONV_MOMENTUM_BUG", cs%convect_mom_bug, &
                 "If true, use code with a bug that causes a loss of momentum conservation "//&
                 "during mixedlayer convection.", 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 "//&
                 "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', 'm')
  cs%id_TKE_wind = register_diag_field('ocean_model', 'TKE_wind', diag%axesT1, &
      time, 'Wind-stirring source of mixed layer TKE', &
      'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
  cs%id_TKE_RiBulk = register_diag_field('ocean_model', 'TKE_RiBulk', diag%axesT1, &
      time, 'Mean kinetic energy source of mixed layer TKE', &
      'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
  cs%id_TKE_conv = register_diag_field('ocean_model', 'TKE_conv', diag%axesT1, &
      time, 'Convective source of mixed layer TKE', &
      'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**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', &
      'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
  cs%id_TKE_mixing = register_diag_field('ocean_model', 'TKE_mixing', diag%axesT1, &
      time, 'TKE consumed by mixing that deepens the mixed layer', &
      'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**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', &
      'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**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', &
      'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**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', &
      'm3 s-3', conversion=us%Z_to_m*(us%L_to_m**2)*(us%s_to_T**3))
  cs%id_PE_detrain = register_diag_field('ocean_model', 'PE_detrain', diag%axesT1, &
      time, 'Spurious source of potential energy from mixed layer detrainment', &
      'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**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', &
      'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)
  cs%id_h_mismatch = register_diag_field('ocean_model', 'h_miss_ML', diag%axesT1, &
      time, 'Summed absolute mismatch in entrainment terms', 'm', conversion=us%Z_to_m)
  cs%id_Hsfc_used = register_diag_field('ocean_model', 'Hs_used', diag%axesT1, &
      time, 'Surface region thickness that is used', 'm', conversion=us%Z_to_m)
  cs%id_Hsfc_max = register_diag_field('ocean_model', 'Hs_max', diag%axesT1, &
      time, 'Maximum surface region thickness', 'm', conversion=us%Z_to_m)
  cs%id_Hsfc_min = register_diag_field('ocean_model', 'Hs_min', diag%axesT1, &
      time, 'Minimum surface region thickness', 'm', conversion=us%Z_to_m)
 !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 "//&
                 "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 "//&
                 "of the surface region (mixed & buffer layer) is allowed "//&
                 "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 "//&
                 "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/((Ht+E)*Ht).
function ef4(Ht, En, I_L, dR_de)
  real,           intent(in)    :: ht  !< Total thickness [H ~> m or kg m-2].
  real,           intent(in)    :: en  !< Entrainment [H ~> m or kg m-2].
  real,           intent(in)    :: i_l !< The e-folding scale [H-1 ~> m-1 or m2 kg-1]
  real, optional, intent(inout) :: dr_de !< The partial derivative of the result R with E [H-2 ~> m-2 or m4 kg-2].
  real :: ef4 !< The integral [H-1 ~> m-1 or m2 kg-1].

  ! Local variables
  real :: exp_lhpe ! A nondimensional exponential decay.
  real :: i_hpe    ! An inverse thickness plus entrainment [H-1 ~> m-1 or m2 kg-1].
  real :: res      ! The result of the integral above [H-1 ~> m-1 or m2 kg-1].

  exp_lhpe = exp(-i_l*(en+ht))
  i_hpe = 1.0/(ht+en)
  res = exp_lhpe * (en*i_hpe/ht - 0.5*i_l*log(ht*i_hpe) + 0.5*i_l*i_l*en)
  if (PRESENT(dr_de)) &
    dr_de = -i_l*res + exp_lhpe*(i_hpe*i_hpe + 0.5*i_l*i_hpe + 0.5*i_l*i_l)
  ef4 = res

end function ef4

!> \namespace mom_bulk_mixed_layer
!!
!! 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.

end module mom_bulk_mixed_layer