MOM_entrain_diffusive.F90

!> Diapycnal mixing and advection in isopycnal mode
module mom_entrain_diffusive

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

use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : diag_ctrl, time_type
use mom_error_handler, only : mom_error, is_root_pe, fatal, warning, note
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_forcing_type, only : forcing
use mom_grid, only : ocean_grid_type
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 entrainment_diffusive, entrain_diffusive_init, entrain_diffusive_end

! 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 holding parametes for the MOM_entrain_diffusive module
type, public :: entrain_diffusive_cs ; private
  logical :: bulkmixedlayer  !< If true, a refined bulk mixed layer is used with
                             !! GV%nk_rho_varies variable density mixed & buffer layers.
  integer :: max_ent_it      !< The maximum number of iterations that may be used to
                             !! calculate the diapycnal entrainment.
  real    :: tolerance_ent   !< The tolerance with which to solve for entrainment values
                             !! [H ~> m or kg m-2].
  real    :: rho_sig_off     !< The offset between potential density and a sigma value [R ~> kg m-3]
  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
                             !! regulate the timing of diagnostic output.
  integer :: id_kd = -1      !< Diagnostic ID for diffusivity
  integer :: id_diff_work = -1 !< Diagnostic ID for mixing work
end type entrain_diffusive_cs

contains

!> This subroutine calculates ea and eb, the rates at which a layer entrains
!! from the layers above and below.  The entrainment rates are proportional to
!! the buoyancy flux in a layer and inversely proportional to the density
!! differences between layers.  The scheme that is used here is described in
!! detail in Hallberg, Mon. Wea. Rev. 2000.
subroutine entrainment_diffusive(h, tv, fluxes, dt, G, GV, US, CS, ea, eb, &
                                 kb_out, Kd_Lay, Kd_int)
  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),SZJ_(G),SZK_(G)),   &
                              intent(in)  :: h  !< Layer thicknesses [H ~> m or kg m-2].
  type(thermo_var_ptrs),      intent(in)  :: tv !< A structure containing pointers to any available
                                                !! thermodynamic fields. Absent fields have NULL
                                                !! ptrs.
  type(forcing),              intent(in)  :: fluxes !< A structure of surface fluxes that may
                                                !! be used.
  real,                       intent(in)  :: dt !< The time increment [T ~> s].
  type(entrain_diffusive_cs), pointer     :: CS !< The control structure returned by a previous
                                                !! call to entrain_diffusive_init.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   &
                              intent(out) :: ea !< The amount of fluid entrained from the layer
                                                !! above within this time step [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   &
                              intent(out) :: eb !< The amount of fluid entrained from the layer
                                                !! below within this time step [H ~> m or kg m-2].
  integer, dimension(SZI_(G),SZJ_(G)),        &
                  optional, intent(inout) :: kb_out !< The index of the lightest layer denser than
                                                !! the buffer layer.
  ! At least one of the two following arguments must be present.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   &
                  optional, intent(in)    :: Kd_Lay !< The diapycnal diffusivity of layers
                                                !! [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
                  optional, intent(in)    :: Kd_int !< The diapycnal diffusivity of interfaces
                                                !! [Z2 T-1 ~> m2 s-1].

!   This subroutine calculates ea and eb, the rates at which a layer entrains
! from the layers above and below.  The entrainment rates are proportional to
! the buoyancy flux in a layer and inversely proportional to the density
! differences between layers.  The scheme that is used here is described in
! detail in Hallberg, Mon. Wea. Rev. 2000.

  real, dimension(SZI_(G),SZK_(G)) :: &
    dtKd    ! The layer diapycnal diffusivity times the time step [H2 ~> m2 or kg2 m-4].
  real, dimension(SZI_(G),SZK_(G)+1) :: &
    dtKd_int ! The diapycnal diffusivity at the interfaces times the time step [H2 ~> m2 or kg2 m-4]
  real, dimension(SZI_(G),SZK_(G)) :: &
    F, &    ! The density flux through a layer within a time step divided by the
            ! density difference across the interface below the layer [H ~> m or kg m-2].
    maxf, & ! maxF is the maximum value of F that will not deplete all of the
            ! layers above or below a layer within a timestep [H ~> m or kg m-2].
    minf, & ! minF is the minimum flux that should be expected in the absence of
            ! interactions between layers [H ~> m or kg m-2].
    fprev, &! The previous estimate of F [H ~> m or kg m-2].
    dfdfm, &! The partial derivative of F with respect to changes in F of the
            ! neighboring layers.  [nondim]
    h_guess ! An estimate of the layer thicknesses after entrainment, but
            ! before the entrainments are adjusted to drive the layer
            ! densities toward their target values [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)+1) :: &
    Ent_bl  ! The average entrainment upward and downward across
            ! each interface around the buffer layers [H ~> m or kg m-2].
  real, allocatable, dimension(:,:,:) :: &
    Kd_eff, &     ! The effective diffusivity that actually applies to each
                  ! layer after the effects of boundary conditions are
                  ! considered [Z2 T-1 ~> m2 s-1].
    diff_work     ! The work actually done by diffusion across each
                  ! interface [R Z3 T-3 ~> W m-2].  Sum vertically for the total work.

  real :: hm, fm, fr, fk  ! Work variables with units of H, H, H, and H2.

  real :: b1(szi_(g))         ! b1 and c1 are variables used by the
  real :: c1(szi_(g),szk_(g)) ! tridiagonal solver.

  real, dimension(SZI_(G)) :: &
    htot, &       ! The total thickness above or below a layer [H ~> m or kg m-2].
    Rcv, &        ! Value of the coordinate variable (potential density)
                  ! based on the simulated T and S and P_Ref [R ~> kg m-3].
    pres, &       ! Reference pressure (P_Ref) [R L2 T-2 ~> Pa].
    eakb, &       ! The entrainment from above by the layer below the buffer
                  ! layer (i.e. layer kb) [H ~> m or kg m-2].
    ea_kbp1, &    ! The entrainment from above by layer kb+1 [H ~> m or kg m-2].
    eb_kmb, &     ! The entrainment from below by the deepest buffer layer [H ~> m or kg m-2].
    ds_kb, &      ! The reference potential density difference across the
                  ! interface between the buffer layers and layer kb [R ~> kg m-3].
    ds_anom_lim, &! The amount by which dS_kb is reduced when limits are
                  ! applied [R ~> kg m-3].
    i_dskbp1, &   ! The inverse of the potential density difference across the
                  ! interface below layer kb [R-1 ~> m3 kg-1].
    dtkd_kb, &    ! The diapycnal diffusivity in layer kb times the time step
                  ! [H2 ~> m2 or kg2 m-4].
    maxf_correct, & ! An amount by which to correct maxF due to excessive
                  ! surface heat loss [H ~> m or kg m-2].
    zeros, &      ! An array of all zeros. (Usually used with [H ~> m or kg m-2].)
    max_eakb, &   ! The maximum value of eakb that might be realized [H ~> m or kg m-2].
    min_eakb, &   ! The minimum value of eakb that might be realized [H ~> m or kg m-2].
    err_max_eakb0, & ! The value of error returned by determine_Ea_kb
    err_min_eakb0, & ! when eakb = min_eakb and max_eakb and ea_kbp1 = 0.
    err_eakb0, &  ! A value of error returned by determine_Ea_kb.
    f_kb, &       ! The value of F in layer kb, or equivalently the entrainment
                  ! from below by layer kb [H ~> m or kg m-2].
    dfdfm_kb, &   ! The partial derivative of F with fm [nondim]. See dFdfm.
    maxf_kb, &    ! The maximum value of F_kb that might be realized [H ~> m or kg m-2].
    eakb_maxf, &  ! The value of eakb that gives F_kb=maxF_kb [H ~> m or kg m-2].
    f_kb_maxent   ! The value of F_kb when eakb = max_eakb [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)) :: &
    Sref, &  ! The reference potential density of the mixed and buffer layers,
             ! and of the two lightest interior layers (kb and kb+1) copied
             ! into layers kmb+1 and kmb+2 [R ~> kg m-3].
    h_bl     ! The thicknesses of the mixed and buffer layers, and of the two
             ! lightest interior layers (kb and kb+1) copied into layers kmb+1
             ! and kmb+2 [H ~> m or kg m-2].

  real, dimension(SZI_(G),SZK_(G)) :: &
    ds_dsp1, &      ! The coordinate variable (sigma-2) difference across an
                    ! interface divided by the difference across the interface
                    ! below it. [nondim]
    dsp1_ds, &      ! The inverse coordinate variable (sigma-2) difference
                    ! across an interface times the difference across the
                    ! interface above it. [nondim]
    i2p2dsp1_ds, &  ! 1 / (2 + 2 * ds_k+1 / ds_k). [nondim]
    grats           ! 2*(2 + ds_k+1 / ds_k + ds_k / ds_k+1) =
                    !       4*ds_Lay*(1/ds_k + 1/ds_k+1). [nondim]

  real :: dRHo      ! The change in locally referenced potential density between
                    ! the layers above and below an interface [R ~> kg m-3].
  real :: g_2dt     ! 0.5 * G_Earth / dt, times unit conversion factors
                    ! [m3 H-2 s-2 T-1 ~> m s-3 or m7 kg-2 s-3].
  real, dimension(SZI_(G)) :: &
    pressure, &      ! The pressure at an interface [R L2 T-2 ~> Pa].
    T_eos, S_eos, &  ! The potential temperature and salinity at which to
                     ! evaluate dRho_dT and dRho_dS [degC] and [ppt].
    drho_dt, drho_ds ! The partial derivatives of potential density with temperature and
                     ! salinity, [R degC-1 ~> kg m-3 degC-1] and [R ppt-1 ~> kg m-3 ppt-1].

  real :: tolerance  ! The tolerance within which E must be converged [H ~> m or kg m-2].
  real :: Angstrom   ! The minimum layer thickness [H ~> m or kg m-2].
  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 :: F_cor      ! A correction to the amount of F that is used to
                     ! entrain from the layer above [H ~> m or kg m-2].
  real :: Kd_here    ! The effective diapycnal diffusivity [H2 s-1 ~> m2 s-1 or kg2 m-4 s-1].
  real :: h_avail    ! The thickness that is available for entrainment [H ~> m or kg m-2].
  real :: dS_kb_eff  ! The value of dS_kb after limiting is taken into account.
  real :: Rho_cor    ! The depth-integrated potential density anomaly that
                     ! needs to be corrected for [H R ~> kg m-2 or kg2 m-5].
  real :: ea_cor     ! The corrective adjustment to eakb [H ~> m or kg m-2].
  real :: h1         ! The layer thickness after entrainment through the
                     ! interface below is taken into account [H ~> m or kg m-2].
  real :: Idt        ! The inverse of the time step [T-1 ~> s-1].

  logical :: do_any
  logical :: do_entrain_eakb    ! True if buffer layer is entrained
  logical :: do_i(szi_(g)), did_i(szi_(g)), reiterate
  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
  integer :: it, i, j, k, is, ie, js, je, nz, K2, kmb
  integer :: kb(szi_(g))  ! The value of kb in row j.
  integer :: kb_min       ! The minimum value of kb in the current j-row.
  integer :: kb_min_act   ! The minimum active value of kb in the current j-row.
  integer :: is1, ie1     ! The minimum and maximum active values of i in the current j-row.
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  angstrom = gv%Angstrom_H
  h_neglect = gv%H_subroundoff

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

  if (.not.(present(kd_lay) .or. present(kd_int))) call mom_error(fatal, &
      "MOM_entrain_diffusive: Either Kd_Lay or Kd_int must be present in call.")

  if ((.not.cs%bulkmixedlayer .and. .not.associated(fluxes%buoy)) .and. &
      (associated(fluxes%lprec) .or. associated(fluxes%evap) .or. &
       associated(fluxes%sens) .or. associated(fluxes%sw))) then
    if (is_root_pe()) call mom_error(note, "Calculate_Entrainment: &
          &The code to handle evaporation and precipitation without &
          &a bulk mixed layer has not been implemented.")
    if (is_root_pe()) call mom_error(fatal, &
         "Either define BULKMIXEDLAYER in MOM_input or use fluxes%buoy &
         &and a linear equation of state to drive the model.")
  endif

  tolerance = cs%Tolerance_Ent
  kmb = gv%nk_rho_varies
  k2 = max(kmb+1,2) ; kb_min = k2
  if (.not. cs%bulkmixedlayer) then
    kb(:) = 1
    ! These lines fill in values that are arbitrary, but needed because
    ! they are used to normalize the buoyancy flux in layer nz.
    do i=is,ie ; ds_dsp1(i,nz) = 2.0 ; dsp1_ds(i,nz) = 0.5 ; enddo
  else
    kb(:) = 0
    do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; dsp1_ds(i,nz) = 0.0 ; enddo
  endif

  if (cs%id_diff_work > 0) allocate(diff_work(g%isd:g%ied,g%jsd:g%jed,nz+1))
  if (cs%id_Kd > 0)        allocate(kd_eff(g%isd:g%ied,g%jsd:g%jed,nz))

  if (associated(tv%eqn_of_state)) then
    pres(:) = tv%P_Ref
  else
    pres(:) = 0.0
  endif
  eosdom(:) = eos_domain(g%HI)

  !$OMP parallel do default(none) shared(is,ie,js,je,nz,Kd_Lay,G,GV,US,dt,CS,h,tv,   &
  !$OMP                                  kmb,Angstrom,fluxes,K2,h_neglect,tolerance, &
  !$OMP                                  ea,eb,Kd_int,Kd_eff,EOSdom,diff_work,g_2dt, kb_out) &
  !$OMP                     firstprivate(kb,ds_dsp1,dsp1_ds,pres,kb_min)             &
  !$OMP                          private(dtKd,dtKd_int,do_i,Ent_bl,dtKd_kb,h_bl,     &
  !$OMP                                  I2p2dsp1_ds,grats,htot,max_eakb,I_dSkbp1,   &
  !$OMP                                  zeros,maxF_kb,maxF,ea_kbp1,eakb,Sref,       &
  !$OMP                                  maxF_correct,do_any,do_entrain_eakb,        &
  !$OMP                                  err_min_eakb0,err_max_eakb0,eakb_maxF,      &
  !$OMP                                  min_eakb,err_eakb0,F,minF,hm,fk,F_kb_maxent,&
  !$OMP                                  F_kb,is1,ie1,kb_min_act,dFdfm_kb,b1,dFdfm,  &
  !$OMP                                  Fprev,fm,fr,c1,reiterate,eb_kmb,did_i,      &
  !$OMP                                  h_avail,h_guess,dS_kb,Rcv,F_cor,dS_kb_eff,  &
  !$OMP                                  Rho_cor,ea_cor,h1,Idt,Kd_here,pressure,     &
  !$OMP                                  T_eos,S_eos,dRho_dT,dRho_dS,dRho,dS_anom_lim)
  do j=js,je
    do i=is,ie ; kb(i) = 1 ; enddo

    if (present(kd_lay)) then
      do k=1,nz ; do i=is,ie
        dtkd(i,k) = gv%Z_to_H**2 * (dt * kd_lay(i,j,k))
      enddo ; enddo
      if (present(kd_int)) then
        do k=1,nz+1 ; do i=is,ie
          dtkd_int(i,k) = gv%Z_to_H**2 * (dt * kd_int(i,j,k))
        enddo ; enddo
      else
        do k=2,nz ; do i=is,ie
          dtkd_int(i,k) = gv%Z_to_H**2 * (0.5 * dt * (kd_lay(i,j,k-1) + kd_lay(i,j,k)))
        enddo ; enddo
      endif
    else ! Kd_int must be present, or there already would have been an error.
      do k=1,nz ; do i=is,ie
        dtkd(i,k) = gv%Z_to_H**2 * (0.5 * dt * (kd_int(i,j,k)+kd_int(i,j,k+1)))
      enddo ; enddo
      dO k=1,nz+1 ; do i=is,ie
        dtkd_int(i,k) = gv%Z_to_H**2 * (dt * kd_int(i,j,k))
      enddo ; enddo
    endif

    do i=is,ie ; do_i(i) = (g%mask2dT(i,j) > 0.5) ; enddo
    do i=is,ie ; ds_dsp1(i,nz) = 0.0 ; enddo
    do i=is,ie ; dsp1_ds(i,nz) = 0.0 ; enddo

    do k=2,nz-1 ; do i=is,ie
      ds_dsp1(i,k) = gv%g_prime(k) / gv%g_prime(k+1)
    enddo ; enddo

    if (cs%bulkmixedlayer) then
      !   This subroutine determines the averaged entrainment across each
      ! interface and causes thin and relatively light interior layers to be
      ! entrained by the deepest buffer layer.  This also determines kb.
      call set_ent_bl(h, dtkd_int, tv, kb, kmb, do_i, g, gv, us, cs, j, ent_bl, sref, h_bl)

      do i=is,ie
        dtkd_kb(i) = 0.0 ; if (kb(i) < nz) dtkd_kb(i) = dtkd(i,kb(i))
      enddo
    else
      do i=is,ie ; ent_bl(i,kmb+1) = 0.0 ; enddo
    endif

    do k=2,nz-1 ; do i=is,ie
      dsp1_ds(i,k) = 1.0 / ds_dsp1(i,k)
      i2p2dsp1_ds(i,k) = 0.5/(1.0+dsp1_ds(i,k))
      grats(i,k) = 2.0*(2.0+(dsp1_ds(i,k)+ds_dsp1(i,k)))
    enddo ; enddo

!   Determine the maximum flux, maxF, for each of the isopycnal layers.
!   Also determine when the fluxes start entraining
! from various buffer or mixed layers, where appropriate.
    if (cs%bulkmixedlayer) then
      kb_min = nz
      do i=is,ie
        htot(i) = h(i,j,1) - angstrom
      enddo
      do k=2,kmb ; do i=is,ie
        htot(i) = htot(i) + (h(i,j,k) - angstrom)
      enddo ; enddo
      do i=is,ie
        max_eakb(i) = max(ent_bl(i,kmb+1) + 0.5*htot(i), htot(i))
        i_dskbp1(i) = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))
        zeros(i) = 0.0
      enddo

      !   Find the maximum amount of entrainment from below that the top
      ! interior layer could exhibit (maxF(i,kb)), approximating that
      ! entrainment by (eakb*max(dS_kb/dSkbp1,0)).  eakb is in the range
      ! from 0 to max_eakb.
      call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, max_eakb, kmb, &
                        is, ie, g, gv, cs, maxf_kb, eakb_maxf, do_i, f_kb_maxent)
      do i=is,ie ; if (kb(i) <= nz) then
        maxf(i,kb(i)) = max(0.0, maxf_kb(i), f_kb_maxent(i))
        if ((maxf_kb(i) > f_kb_maxent(i)) .and. (eakb_maxf(i) >= htot(i))) &
          max_eakb(i) = eakb_maxf(i)
      endif ; enddo

      do i=is,ie ; ea_kbp1(i) = 0.0 ; eakb(i) = max_eakb(i) ; enddo
      call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, &
                           max_eakb, max_eakb, kmb, is, ie, do_i, g, gv, cs, eakb, &
                           error=err_max_eakb0, f_kb=f_kb)

      !   The maximum value of F(kb) between htot and max_eakb determines
      ! what maxF(kb+1) should be.
      do i=is,ie ; min_eakb(i) = min(htot(i), max_eakb(i)) ; enddo
      call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, min_eakb, max_eakb, &
                        kmb, is, ie, g, gv, cs, f_kb_maxent, do_i_in = do_i)

      do i=is,ie
        do_entrain_eakb = .false.
        ! If error_max_eakb0 < 0, then buffer layers are always all entrained
        if (do_i(i)) then ; if (err_max_eakb0(i) < 0.0) then
          do_entrain_eakb = .true.
        endif ; endif

        if (do_entrain_eakb) then
          eakb(i) = max_eakb(i) ; min_eakb(i) = max_eakb(i)
        else
          eakb(i) = 0.0 ; min_eakb(i) = 0.0
        endif
      enddo

      !   Find the amount of entrainment of the buffer layers that would occur
      ! if there were no entrainment by the deeper interior layers.  Also find
      ! how much entrainment of the deeper layers would occur.
      call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, &
                           zeros, max_eakb, kmb, is, ie, do_i, g, gv, cs, min_eakb, &
                           error=err_min_eakb0, f_kb=f_kb, err_max_eakb0=err_max_eakb0)
      ! Error_min_eakb0 should be ~0 unless error_max_eakb0 < 0.
      do i=is,ie ; if ((kb(i)<nz) .and. (kb_min>kb(i))) kb_min = kb(i) ; enddo
    else
      ! Without a bulk mixed layer, surface fluxes are applied in this
      ! subroutine.  (Otherwise, they are handled in mixedlayer.)
      !   Initially the maximum flux in layer zero is given by the surface
      ! buoyancy flux.  It will later be limited if the surface flux is
      ! too large.  Here buoy is the surface buoyancy flux.
      do i=is,ie
        maxf(i,1) = 0.0
        htot(i) = h(i,j,1) - angstrom
      enddo
      if (associated(fluxes%buoy)) then ; do i=is,ie
        maxf(i,1) = gv%Z_to_H * (dt*fluxes%buoy(i,j)) / gv%g_prime(2)
      enddo ; endif
    endif

! The following code calculates the maximum flux, maxF, for the interior
! layers.
    do k=kb_min,nz-1 ; do i=is,ie
      if ((k == kb(i)+1) .and. cs%bulkmixedlayer) then
        maxf(i,k) = ds_dsp1(i,k)*(f_kb_maxent(i) + htot(i))
        htot(i) = htot(i) + (h(i,j,k) - angstrom)
      elseif (k > kb(i)) then
        maxf(i,k) = ds_dsp1(i,k)*(maxf(i,k-1) + htot(i))
        htot(i) = htot(i) + (h(i,j,k) - angstrom)
      endif
    enddo ; enddo
    do i=is,ie
      maxf(i,nz) = 0.0
      if (.not.cs%bulkmixedlayer) then
        maxf_correct(i) = max(0.0, -(maxf(i,nz-1) + htot(i)))
      endif
      htot(i) = h(i,j,nz) - angstrom
    enddo
    if (.not.cs%bulkmixedlayer) then
      do_any = .false. ; do i=is,ie ; if (maxf_correct(i) > 0.0) do_any = .true. ; enddo
      if (do_any) then
        do k=nz-1,1,-1 ; do i=is,ie
          maxf(i,k) = maxf(i,k) + maxf_correct(i)
          maxf_correct(i) = maxf_correct(i) * dsp1_ds(i,k)
        enddo ; enddo
      endif
    endif
    do k=nz-1,kb_min,-1 ; do i=is,ie ; if (do_i(i)) then
      if (k >= kb(i)) then
        maxf(i,k) = min(maxf(i,k),dsp1_ds(i,k+1)*maxf(i,k+1) + htot(i))
        htot(i) = htot(i) + (h(i,j,k) - angstrom)
      endif
      if (k == kb(i)) then
        if ((maxf(i,k) < f_kb(i)) .or. (maxf(i,k) < maxf_kb(i)) &
            .and. (eakb_maxf(i) <= max_eakb(i))) then
          !   In this case, too much was being entrained by the topmost interior
          ! layer, even with the minimum initial estimate.  The buffer layer
          ! will always entrain the maximum amount.
          f_kb(i) = maxf(i,k)
          if ((f_kb(i) <= maxf_kb(i)) .and. (eakb_maxf(i) <= max_eakb(i))) then
            eakb(i) = eakb_maxf(i)
          else
            eakb(i) = max_eakb(i)
          endif
          call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &
                             g, gv, cs, eakb, angstrom)
          if ((eakb(i) < max_eakb(i)) .or. (eakb(i) < min_eakb(i))) then
            call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, zeros, &
                                 eakb, eakb, kmb, i, i, do_i, g, gv, cs, eakb, &
                                 error=err_eakb0)
            if (eakb(i) < max_eakb(i)) then
              max_eakb(i) = eakb(i) ; err_max_eakb0(i) = err_eakb0(i)
            endif
            if (eakb(i) < min_eakb(i)) then
              min_eakb(i) = eakb(i) ; err_min_eakb0(i) = err_eakb0(i)
            endif
          endif
        endif
      endif
    endif ; enddo ; enddo
    if (.not.cs%bulkmixedlayer) then
      do i=is,ie
        maxf(i,1) = min(maxf(i,1),dsp1_ds(i,2)*maxf(i,2) + htot(i))
      enddo
    endif

!   The following code provides an initial estimate of the flux in
! each layer, F.  The initial guess for the layer diffusive flux is
! the smaller of a forward discretization or the maximum diffusive
! flux starting from zero thickness in one time step without
! considering adjacent layers.
    do i=is,ie
      f(i,1) = maxf(i,1)
      f(i,nz) = maxf(i,nz) ; minf(i,nz) = 0.0
    enddo
    do k=nz-1,k2,-1
      do i=is,ie
        if ((k==kb(i)) .and. (do_i(i))) then
          eakb(i) = min_eakb(i)
          minf(i,k) = 0.0
        elseif ((k>kb(i)) .and. (do_i(i))) then
!   Here the layer flux is estimated, assuming no entrainment from
! the surrounding layers.  The estimate is a forward (steady) flux,
! limited by the maximum flux for a layer starting with zero
! thickness.  This is often a good guess and leads to few iterations.
          hm = h(i,j,k) + h_neglect
  !   Note: Tried sqrt((0.5*ds_dsp1(i,k))*dtKd(i,k)) for the second limit,
  ! but it usually doesn't work as well.
          f(i,k) = min(maxf(i,k), sqrt(ds_dsp1(i,k)*dtkd(i,k)), &
                       0.5*(ds_dsp1(i,k)+1.0) * (dtkd(i,k) / hm))

!   Calculate the minimum flux that can be expected if there is no entrainment
! from the neighboring layers.  The 0.9 is used to give used to give a 10%
! known error tolerance.
          fk = dtkd(i,k) * grats(i,k)
          minf(i,k) = min(maxf(i,k), &
                          0.9*(i2p2dsp1_ds(i,k) * fk / (hm + sqrt(hm*hm + fk))))
          if (k==kb(i)) minf(i,k) = 0.0 ! BACKWARD COMPATIBILITY - DELETE LATER?
        else
          f(i,k) = 0.0
          minf(i,k) = 0.0
        endif
      enddo ! end of i loop
    enddo ! end of k loop

    ! This is where the fluxes are actually calculated.

    is1 = ie+1 ; ie1 = is-1
    do i=is,ie ; if (do_i(i)) then ; is1 = i ; exit ; endif ; enddo
    do i=ie,is,-1 ; if (do_i(i)) then ; ie1 = i ; exit ; endif ; enddo

    if (cs%bulkmixedlayer) then
      kb_min_act = nz
      do i=is,ie
        if (do_i(i) .and. (kb(i) < kb_min_act)) kb_min_act = kb(i)
      enddo
      !   Solve for the entrainment rate from above in the topmost interior
      ! layer, eakb, such that
      !   eakb*dS_implicit = dt*Kd*dS_layer_implicit / h_implicit.
      do i=is1,ie1
        ea_kbp1(i) = 0.0
        if (do_i(i) .and. (kb(i) < nz)) &
          ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*f(i,kb(i)+1)
      enddo
      call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &
                           max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &
                           err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
                           dfdfm_kb=dfdfm_kb)
    else
      kb_min_act = kb_min
    endif

    do it=0,cs%max_ent_it-1
      do i=is1,ie1 ; if (do_i(i)) then
        if (.not.cs%bulkmixedlayer) f(i,1) = min(f(i,1),maxf(i,1))
        b1(i) = 1.0
      endif ; enddo ! end of i loop

     ! F_kb has already been found for this iteration, either above or at
     ! the end of the code for the previous iteration.
      do k=kb_min_act,nz-1 ; do i=is1,ie1 ; if (do_i(i) .and. (k>=kb(i))) then
        ! Calculate the flux in layer k.
        if (cs%bulkmixedlayer .and. (k==kb(i))) then
          f(i,k) = f_kb(i)
          dfdfm(i,k) = dfdfm_kb(i)
        else ! k > kb(i)
          fprev(i,k) = f(i,k)
          fm = (f(i,k-1) - h(i,j,k)) + dsp1_ds(i,k+1)*f(i,k+1)
          fk = grats(i,k)*dtkd(i,k)
          fr = sqrt(fm*fm + fk)

          if (fm>=0) then
            f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fm+fr))
          else
            f(i,k) = min(maxf(i,k), i2p2dsp1_ds(i,k) * (fk / (-1.0*fm+fr)))
          endif

          if ((f(i,k) >= maxf(i,k)) .or. (fr == 0.0)) then
            dfdfm(i,k) = 0.0
          else
            dfdfm(i,k) = i2p2dsp1_ds(i,k) * ((fr + fm) / fr)
          endif

          if (k > k2) then
            ! This is part of a tridiagonal solver for the actual flux.
            c1(i,k) = dfdfm(i,k-1)*(dsp1_ds(i,k)*b1(i))
            b1(i) = 1.0 / (1.0 - c1(i,k)*dfdfm(i,k))
            f(i,k) = min(b1(i)*(f(i,k)-fprev(i,k)) + fprev(i,k), maxf(i,k))
            if (f(i,k) >= maxf(i,k)) dfdfm(i,k) = 0.0
          endif
        endif
      endif ; enddo ; enddo

      do k=nz-2,kb_min_act,-1 ; do i=is1,ie1
        if (do_i(i) .and. (k > kb(i))) &
          f(i,k) = min((f(i,k)+c1(i,k+1)*(f(i,k+1)-fprev(i,k+1))),maxf(i,k))
      enddo ; enddo

      if (cs%bulkmixedlayer) then
        do i=is1,ie1
          if (do_i(i) .and. (kb(i) < nz)) then
            ! F will be increased to minF later.
            ea_kbp1(i) = dsp1_ds(i,kb(i)+1)*max(f(i,kb(i)+1), minf(i,kb(i)+1))
          else
            ea_kbp1(i) = 0.0
          endif
        enddo
        call determine_ea_kb(h_bl, dtkd_kb, sref, i_dskbp1, ent_bl, ea_kbp1, min_eakb, &
                             max_eakb, kmb, is1, ie1, do_i, g, gv, cs, eakb, f_kb=f_kb, &
                             err_max_eakb0=err_max_eakb0, err_min_eakb0=err_min_eakb0, &
                             dfdfm_kb=dfdfm_kb)
        do i=is1,ie1
          if (do_i(i) .and. (kb(i) < nz)) f(i,kb(i)) = f_kb(i)
        enddo
      endif

! Determine whether to do another iteration.
      if (it < cs%max_ent_it-1) then

        reiterate = .false.
        if (cs%bulkmixedlayer) then ; do i=is1,ie1 ; if (do_i(i)) then
          eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
        endif ; enddo ; endif
        do i=is1,ie1
          did_i(i) = do_i(i) ; do_i(i) = .false.
        enddo
        do k=kb_min_act,nz-1 ; do i=is1,ie1
          if (did_i(i) .and. (k >= kb(i))) then
            if (f(i,k) < minf(i,k)) then
              f(i,k) = minf(i,k)
              do_i(i) = .true. ; reiterate = .true.
            elseif (k > kb(i)) then
              if ((abs(f(i,k) - fprev(i,k)) > tolerance) .or. &
                  ((h(i,j,k) + ((1.0+dsp1_ds(i,k))*f(i,k) - &
                   (f(i,k-1) + dsp1_ds(i,k+1)*f(i,k+1)))) < 0.9*angstrom)) then
                do_i(i) = .true. ; reiterate = .true.
              endif
            else ! (k == kb(i))
! A more complicated test is required for the layer beneath the buffer layer,
! since its flux may be partially used to entrain directly from the mixed layer.
! Negative fluxes should not occur with the bulk mixed layer.
              if (h(i,j,k) + ((f(i,k) + eakb(i)) - &
                  (eb_kmb(i) + dsp1_ds(i,k+1)*f(i,k+1))) < 0.9*angstrom) then
                do_i(i) = .true. ; reiterate = .true.
              endif
            endif
          endif
        enddo ; enddo
        if (.not.reiterate) exit
      endif ! end of if (it < CS%max_ent_it-1)
    enddo ! end of it loop
! This is the end of the section that might be iterated.


    if (it == (cs%max_ent_it)) then
      !   Limit the flux so that the layer below is not depleted.
      ! This should only be applied to the last iteration.
      do i=is1,ie1 ; if (do_i(i)) then
        f(i,nz-1) = max(f(i,nz-1), min(minf(i,nz-1), 0.0))
        if (kb(i) >= nz-1) then ; ea_kbp1(i) = 0.0 ; endif
      endif ; enddo
      do k=nz-2,kb_min_act,-1 ; do i=is1,ie1 ; if (do_i(i)) then
        if (k>kb(i)) then
          f(i,k) = min(max(minf(i,k),f(i,k)), (dsp1_ds(i,k+1)*f(i,k+1) + &
               max(((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + &
                    (h(i,j,k+1) - angstrom)), 0.5*(h(i,j,k+1)-angstrom))))
        elseif (k==kb(i)) then
          ea_kbp1(i) = dsp1_ds(i,k+1)*f(i,k+1)
          h_avail = dsp1_ds(i,k+1)*f(i,k+1) + max(0.5*(h(i,j,k+1)-angstrom), &
               ((f(i,k+1)-dsp1_ds(i,k+2)*f(i,k+2)) + (h(i,j,k+1) - angstrom)))
          if ((f(i,k) > 0.0) .and. (f(i,k) > h_avail)) then
            f_kb(i) = max(0.0, h_avail)
            f(i,k) = f_kb(i)
            if ((f_kb(i) < maxf_kb(i)) .and. (eakb_maxf(i) <= eakb(i))) &
              eakb(i) = eakb_maxf(i)
            call f_kb_to_ea_kb(h_bl, sref, ent_bl, i_dskbp1, f_kb, kmb, i, &
                               g, gv, cs, eakb)
          endif
        endif
      endif ; enddo ; enddo


      if (cs%bulkmixedlayer) then ; do i=is1,ie1
        if (do_i(i) .and. (kb(i) < nz)) then
          h_avail = eakb(i) + max(0.5*(h_bl(i,kmb+1)-angstrom), &
                (f_kb(i)-ea_kbp1(i)) + (h_bl(i,kmb+1)-angstrom))
          ! Ensure that 0 < eb_kmb < h_avail.
          ent_bl(i,kmb+1) = min(ent_bl(i,kmb+1),0.5*(eakb(i) + h_avail))

          eb_kmb(i) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
        endif
      enddo ; endif

     !  Limit the flux so that the layer above is not depleted.
      do k=kb_min_act+1,nz-1 ; do i=is1,ie1 ; if (do_i(i)) then
        if ((.not.cs%bulkmixedlayer) .or. (k > kb(i)+1)) then
          f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + &
              dsp1_ds(i,k-1)*f(i,k-1)) - f(i,k-2)) + (h(i,j,k-1) - angstrom)))
          f(i,k) = max(f(i,k),min(minf(i,k),0.0))
        elseif (k == kb(i)+1) then
          f(i,k) = min(f(i,k), ds_dsp1(i,k)*( ((f(i,k-1) + eakb(i)) - &
              eb_kmb(i)) + (h(i,j,k-1) - angstrom)))
          f(i,k) = max(f(i,k),min(minf(i,k),0.0))
        endif
      endif ; enddo ; enddo
    endif ! (it == (CS%max_ent_it))

    call f_to_ent(f, h, kb, kmb, j, g, gv, cs, dsp1_ds, eakb, ent_bl, ea, eb)

    !   Calculate the layer thicknesses after the entrainment to constrain the
    ! corrective fluxes.
    if (associated(tv%eqn_of_state)) then
      do i=is,ie
        h_guess(i,1) = (h(i,j,1) - angstrom) + (eb(i,j,1) - ea(i,j,2))
        h_guess(i,nz) = (h(i,j,nz) - angstrom) + (ea(i,j,nz) - eb(i,j,nz-1))
        if (h_guess(i,1) < 0.0) h_guess(i,1) = 0.0
        if (h_guess(i,nz) < 0.0) h_guess(i,nz) = 0.0
      enddo
      do k=2,nz-1 ; do i=is,ie
        h_guess(i,k) = (h(i,j,k) - angstrom) + ((ea(i,j,k) - eb(i,j,k-1)) + &
                   (eb(i,j,k) - ea(i,j,k+1)))
        if (h_guess(i,k) < 0.0) h_guess(i,k) = 0.0
      enddo ; enddo
      if (cs%bulkmixedlayer) then
        call determine_dskb(h_bl, sref, ent_bl, eakb, is, ie, kmb, g, gv, &
                            .true., ds_kb, ds_anom_lim=ds_anom_lim)
        do k=nz-1,kb_min,-1
          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, rcv, tv%eqn_of_state, eosdom)
          do i=is,ie
            if ((k>kb(i)) .and. (f(i,k) > 0.0)) then
              ! Within a time step, a layer may entrain no more than its
              ! thickness for correction.  This limitation should apply
              ! extremely rarely, but precludes undesirable behavior.
              !  Note: Corrected a sign/logic error & factor of 2 error, and
              !    the layers tracked the target density better, mostly due to
              !    the factor of 2 error.
              f_cor = h(i,j,k) * min(1.0 , max(-ds_dsp1(i,k), &
                          (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )

              ! Ensure that (1) Entrainments are positive, (2) Corrections in
              ! a layer cannot deplete the layer itself (very generously), and
              ! (3) a layer can take no more than a quarter the mass of its
              ! neighbor.
              if (f_cor > 0.0) then
                f_cor = min(f_cor, 0.9*f(i,k), ds_dsp1(i,k)*0.5*h_guess(i,k), &
                            0.25*h_guess(i,k+1))
              else
                f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &
                             0.25*ds_dsp1(i,k)*h_guess(i,k-1) )
              endif

              ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor
              eb(i,j,k) = eb(i,j,k) + f_cor
            elseif ((k==kb(i)) .and. (f(i,k) > 0.0)) then
              !   Rho_cor is the density anomaly that needs to be corrected,
              ! taking into account that the true potential density of the
              ! deepest buffer layer is not exactly what is returned as dS_kb.
              ds_kb_eff = 2.0*ds_kb(i) - ds_anom_lim(i) ! Could be negative!!!
              rho_cor = h(i,j,k) * (gv%Rlay(k)-rcv(i)) + eakb(i)*ds_anom_lim(i)

              ! Ensure that  -.9*eakb < ea_cor < .9*eakb
              if (abs(rho_cor) < abs(0.9*eakb(i)*ds_kb_eff)) then
                ea_cor = -rho_cor / ds_kb_eff
              else
                ea_cor = sign(0.9*eakb(i),-rho_cor*ds_kb_eff)
              endif

              if (ea_cor > 0.0) then
                ! Ensure that -F_cor < 0.5*h_guess
                ea_cor = min(ea_cor, 0.5*(max_eakb(i) - eakb(i)), &
                             0.5*h_guess(i,k) / (ds_kb(i) * i_dskbp1(i)))
              else
                ! Ensure that -ea_cor < 0.5*h_guess & F_cor < 0.25*h_guess(k+1)
                ea_cor = -min(-ea_cor, 0.5*h_guess(i,k), &
                              0.25*h_guess(i,k+1) / (ds_kb(i) * i_dskbp1(i)))
              endif

              ea(i,j,k) = ea(i,j,k) + ea_cor
              eb(i,j,k) = eb(i,j,k) - (ds_kb(i) * i_dskbp1(i)) * ea_cor
            elseif (k < kb(i)) then
              ! Repetative, unless ea(kb) has been corrected.
              ea(i,j,k) = ea(i,j,k+1)
            endif
          enddo
        enddo
        do k=kb_min-1,k2,-1 ; do i=is,ie
          ea(i,j,k) = ea(i,j,k+1)
        enddo ; enddo

        ! Repetative, unless ea(kb) has been corrected.
        k=kmb
        do i=is,ie
          ! Do not adjust eb through the base of the buffer layers, but it
          ! may be necessary to change entrainment from above.
          h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))
          ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
        enddo
        do k=kmb-1,2,-1 ; do i=is,ie
          ! Determine the entrainment from below for each buffer layer.
          eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)

          ! Determine the entrainment from above for each buffer layer.
          h1 = (h(i,j,k) - angstrom) + (eb(i,j,k) - ea(i,j,k+1))
          ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
        enddo ; enddo
        do i=is,ie
          eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)
        enddo

      else ! not bulkmixedlayer
        do k=k2,nz-1;
          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, rcv, tv%eqn_of_state, eosdom)
          do i=is,ie ; if (f(i,k) > 0.0) then
            ! Within a time step, a layer may entrain no more than
            ! its thickness for correction.  This limitation should
            ! apply extremely rarely, but precludes undesirable
            ! behavior.
            f_cor = h(i,j,k) * min(dsp1_ds(i,k) , max(-1.0, &
                       (gv%Rlay(k) - rcv(i)) / (gv%Rlay(k+1)-gv%Rlay(k))) )

            ! Ensure that (1) Entrainments are positive, (2) Corrections in
            ! a layer cannot deplete the layer itself (very generously), and
            ! (3) a layer can take no more than a quarter the mass of its
            ! neighbor.
            if (f_cor >= 0.0) then
              f_cor = min(f_cor, 0.9*f(i,k), 0.5*dsp1_ds(i,k)*h_guess(i,k), &
                          0.25*h_guess(i,k+1))
            else
              f_cor = -min(-f_cor, 0.9*f(i,k), 0.5*h_guess(i,k), &
                           0.25*ds_dsp1(i,k)*h_guess(i,k-1) )
            endif

            ea(i,j,k) = ea(i,j,k) - dsp1_ds(i,k)*f_cor
            eb(i,j,k) = eb(i,j,k) + f_cor
          endif ; enddo
        enddo
      endif

    endif   !  associated(tv%eqn_of_state))

    if (cs%id_Kd > 0) then
      idt = gv%H_to_Z**2 / dt
      do k=2,nz-1 ; do i=is,ie
        if (k<kb(i)) then ; kd_here = 0.0 ; else
          kd_here = f(i,k) * ( h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
              (eb(i,j,k) - ea(i,j,k+1))) ) / (i2p2dsp1_ds(i,k) * grats(i,k))
        endif

        kd_eff(i,j,k) = max(dtkd(i,k), kd_here)*idt
      enddo ; enddo
      do i=is,ie
        kd_eff(i,j,1) = dtkd(i,1)*idt
        kd_eff(i,j,nz) = dtkd(i,nz)*idt
      enddo
    endif

    if (cs%id_diff_work > 0) then
      g_2dt = 0.5 * gv%H_to_Z**2*us%L_to_Z**2 * (gv%g_Earth / dt)
      do i=is,ie ; diff_work(i,j,1) = 0.0 ; diff_work(i,j,nz+1) = 0.0 ; enddo
      if (associated(tv%eqn_of_state)) then
        if (associated(fluxes%p_surf)) then
          do i=is,ie ; pressure(i) = fluxes%p_surf(i,j) ; enddo
        else
          do i=is,ie ; pressure(i) = 0.0 ; enddo
        endif
        do k=2,nz
          do i=is,ie ; pressure(i) = pressure(i) + (gv%g_Earth*gv%H_to_RZ)*h(i,j,k-1) ; enddo
          do i=is,ie
            if (k==kb(i)) then
              t_eos(i) = 0.5*(tv%T(i,j,kmb) + tv%T(i,j,k))
              s_eos(i) = 0.5*(tv%S(i,j,kmb) + tv%S(i,j,k))
            else
              t_eos(i) = 0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
              s_eos(i) = 0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
            endif
          enddo
          call calculate_density_derivs(t_eos, s_eos, pressure, drho_dt, drho_ds, &
                                        tv%eqn_of_state, eosdom)
          do i=is,ie
            if ((k>kmb) .and. (k<kb(i))) then ; diff_work(i,j,k) = 0.0
            else
              if (k==kb(i)) then
                drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,kmb)) + &
                       drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,kmb))
              else
                drho = drho_dt(i) * (tv%T(i,j,k)-tv%T(i,j,k-1)) + &
                       drho_ds(i) * (tv%S(i,j,k)-tv%S(i,j,k-1))
              endif
              diff_work(i,j,k) = g_2dt * drho * &
                   (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
                    eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
            endif
          enddo
        enddo
      else
        do k=2,nz ; do i=is,ie
          diff_work(i,j,k) = g_2dt * (gv%Rlay(k)-gv%Rlay(k-1)) * &
               (ea(i,j,k) * (h(i,j,k) + ea(i,j,k)) + &
                eb(i,j,k-1)*(h(i,j,k-1) + eb(i,j,k-1)))
        enddo ; enddo
      endif
    endif

    if (present(kb_out)) then
      do i=is,ie ; kb_out(i,j) = kb(i) ; enddo
    endif

  enddo ! end of j loop

! Offer diagnostic fields for averaging.
  if (cs%id_Kd > 0) call post_data(cs%id_Kd, kd_eff, cs%diag)
  if (cs%id_Kd > 0) deallocate(kd_eff)
  if (cs%id_diff_work > 0) call post_data(cs%id_diff_work, diff_work, cs%diag)
  if (cs%id_diff_work > 0) deallocate(diff_work)

end subroutine entrainment_diffusive

!>   This subroutine calculates the actual entrainments (ea and eb) and the
!! amount of surface forcing that is applied to each layer if there is no bulk
!! mixed layer.
subroutine f_to_ent(F, h, kb, kmb, j, G, GV, CS, dsp1_ds, eakb, Ent_bl, ea, eb, do_i_in)
  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_(G)), intent(in)    :: F    !< The density flux through a layer within
                                                          !! a time step divided by the density
                                                          !! difference across the interface below
                                                          !! the layer [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                    intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2]
  integer, dimension(SZI_(G)),      intent(in)    :: kb   !< The index of the lightest layer denser than
                                                          !! the deepest buffer layer.
  integer,                          intent(in)    :: kmb  !< The number of mixed and buffer layers.
  integer,                          intent(in)    :: j    !< The meridional index upon which to work.
  type(entrain_diffusive_cs),       intent(in)    :: CS   !< This module's control structure.
  real, dimension(SZI_(G),SZK_(G)), intent(in)    :: dsp1_ds !< The ratio of coordinate variable
                                                          !! differences across the interfaces below
                                                          !! a layer over the difference across the
                                                          !! interface above the layer.
  real, dimension(SZI_(G)),         intent(in)    :: eakb !< The entrainment from above by the layer
                                                          !! below the buffer layer [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)), intent(in)    :: Ent_bl !< The average entrainment upward and
                                                          !! downward across each interface around
                                                          !! the buffer layers [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                    intent(inout) :: ea   !< The amount of fluid entrained from the layer
                                                          !! above within this time step [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                                    intent(inout) :: eb   !< The amount of fluid entrained from the layer
                                                          !! below within this time step [H ~> m or kg m-2].
  logical, dimension(SZI_(G)), &
                          optional, intent(in)    :: do_i_in !< Indicates which i-points to work on.
!   This subroutine calculates the actual entrainments (ea and eb) and the
! amount of surface forcing that is applied to each layer if there is no bulk
! mixed layer.

  real :: h1        ! The thickness in excess of the minimum that will remain
                    ! after exchange with the layer below [H ~> m or kg m-2].
  logical :: do_i(szi_(g))
  integer :: i, k, is, ie, nz

  is = g%isc ; ie = g%iec ; nz = g%ke

  if (present(do_i_in)) then
    do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
    do i=g%isc,g%iec ; if (do_i(i)) then
      is = i ; exit
    endif ; enddo
    do i=g%iec,g%isc,-1 ; if (do_i(i)) then
      ie = i ; exit
    endif ; enddo
  else
    do i=is,ie ; do_i(i) = .true. ; enddo
  endif

  do i=is,ie
    ea(i,j,nz) = 0.0 ; eb(i,j,nz) = 0.0
  enddo
  if (cs%bulkmixedlayer) then
    do i=is,ie
      eb(i,j,kmb) = max(2.0*ent_bl(i,kmb+1) - eakb(i), 0.0)
    enddo
    do k=nz-1,kmb+1,-1 ; do i=is,ie ; if (do_i(i)) then
      if (k > kb(i)) then
        ! With a bulk mixed layer, surface buoyancy fluxes are applied
        ! elsewhere, so F should always be nonnegative.
        ea(i,j,k) = dsp1_ds(i,k)*f(i,k)
        eb(i,j,k) = f(i,k)
      elseif (k == kb(i)) then
        ea(i,j,k) = eakb(i)
        eb(i,j,k) = f(i,k)
      elseif (k == kb(i)-1) then
        ea(i,j,k) = ea(i,j,k+1)
        eb(i,j,k) = eb(i,j,kmb)
      else
        ea(i,j,k) = ea(i,j,k+1)
        !   Add the entrainment of the thin interior layers to eb going
        ! up into the buffer layer.
        eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom_H)
      endif
    endif ; enddo ; enddo
    k = kmb
    do i=is,ie ; if (do_i(i)) then
      ! Adjust the previously calculated entrainment from below by the deepest
      ! buffer layer to account for entrainment of thin interior layers .
      if (kb(i) > kmb+1) &
        eb(i,j,k) = eb(i,j,k+1) + max(0.0, h(i,j,k+1) - gv%Angstrom_H)

      ! Determine the entrainment from above for each buffer layer.
      h1 = (h(i,j,k) - gv%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
      ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
    endif ; enddo
    do k=kmb-1,2,-1 ; do i=is,ie ; if (do_i(i)) then
      ! Determine the entrainment from below for each buffer layer.
      eb(i,j,k) = max(2.0*ent_bl(i,k+1) - ea(i,j,k+1), 0.0)

      ! Determine the entrainment from above for each buffer layer.
      h1 = (h(i,j,k) - gv%Angstrom_H) + (eb(i,j,k) - ea(i,j,k+1))
      ea(i,j,k) = max(ent_bl(i,k), ent_bl(i,k)-0.5*h1, -h1)
!       if (h1 >= 0.0) then ;                     ea(i,j,k) = Ent_bl(i,K)
!       elseif (Ent_bl(i,K)+0.5*h1 >= 0.0) then ; ea(i,j,k) = Ent_bl(i,K)-0.5*h1
!       else ;                                    ea(i,j,k) = -h1 ; endif
    endif ; enddo ; enddo
    do i=is,ie ; if (do_i(i)) then
      eb(i,j,1) = max(2.0*ent_bl(i,2) - ea(i,j,2), 0.0)
      ea(i,j,1) = 0.0
    endif ; enddo
  else                                          ! not BULKMIXEDLAYER
    ! Calculate the entrainment by each layer from above and below.
    ! Entrainment is always positive, but F may be negative due to
    ! surface buoyancy fluxes.
    do i=is,ie
      ea(i,j,1) = 0.0 ; eb(i,j,1) = max(f(i,1),0.0)
      ea(i,j,2) = dsp1_ds(i,2)*f(i,2) - min(f(i,1),0.0)
    enddo

    do k=2,nz-1 ; do i=is,ie
      eb(i,j,k) = max(f(i,k),0.0)
      ea(i,j,k+1) = dsp1_ds(i,k+1)*f(i,k+1) - (f(i,k)-eb(i,j,k))
      if (ea(i,j,k+1) < 0.0) then
        eb(i,j,k) = eb(i,j,k) - ea(i,j,k+1)
        ea(i,j,k+1) = 0.0
      endif
    enddo ; enddo
  endif                                         ! end BULKMIXEDLAYER
end subroutine f_to_ent

!>   This subroutine sets the average entrainment across each of the interfaces
!! between buffer layers within a timestep. It also causes thin and relatively
!! light interior layers to be entrained by the deepest buffer layer.
!! Also find the initial coordinate potential densities (Sref) of each layer.
subroutine set_ent_bl(h, dtKd_int, tv, kb, kmb, do_i, G, GV, US, CS, j, Ent_bl, Sref, h_bl)
  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),SZJ_(G),SZK_(G)), &
                                    intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZK_(G)+1),       &
                                    intent(in)    :: dtKd_int !< The diapycnal diffusivity across
                                                          !! each interface times the time step
                                                          !! [H2 ~> m2 or kg2 m-4].
  type(thermo_var_ptrs),            intent(in)    :: tv   !< A structure containing pointers to any
                                                          !! available thermodynamic fields. Absent
                                                          !! fields have NULL ptrs.
  integer, dimension(SZI_(G)),      intent(inout) :: kb   !< The index of the lightest layer denser
                                                          !! than the buffer layer or 1 if there is
                                                          !! no buffer layer.
  integer,                          intent(in)    :: kmb  !< The number of mixed and buffer layers.
  logical, dimension(SZI_(G)),      intent(in)    :: do_i !< A logical variable indicating which
                                                          !! i-points to work on.
  type(unit_scale_type),            intent(in)    :: US   !< A dimensional unit scaling type
  type(entrain_diffusive_cs),       pointer       :: CS   !< This module's control structure.
  integer,                          intent(in)    :: j    !< The meridional index upon which to work.
  real, dimension(SZI_(G),SZK_(G)+1),       &
                                    intent(out)   :: Ent_bl !< The average entrainment upward and
                                                          !! downward across each interface around
                                                          !! the buffer layers [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)), intent(out)   :: Sref !< The coordinate potential density minus
                                                          !! 1000 for each layer [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(G)), intent(out)   :: h_bl !< The thickness of each layer [H ~> m or kg m-2].

!   This subroutine sets the average entrainment across each of the interfaces
! between buffer layers within a timestep. It also causes thin and relatively
! light interior layers to be entrained by the deepest buffer layer.
!   Also find the initial coordinate potential densities (Sref) of each layer.
! Does there need to be limiting when the layers below are all thin?

  ! Local variables
  real, dimension(SZI_(G)) :: &
    b1, d1, &   ! Variables used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1] and [nondim].
    Rcv, &      ! Value of the coordinate variable (potential density)
                ! based on the simulated T and S and P_Ref [R ~> kg m-3].
    pres, &     ! Reference pressure (P_Ref) [R L2 T-2 ~> Pa].
    frac_rem, & ! The fraction of the diffusion remaining [nondim].
    h_interior  ! The interior thickness available for entrainment [H ~> m or kg m-2].
  real, dimension(SZI_(G), SZK_(G)) :: &
    S_est ! An estimate of the coordinate potential density - 1000 after
          ! entrainment for each layer [R ~> kg m-3].
  real :: max_ent  ! The maximum possible entrainment [H ~> m or kg m-2].
  real :: dh       ! An available thickness [H ~> m or kg m-2].
  real :: Kd_x_dt  ! The diffusion that remains after thin layers are
                   ! entrained [H2 ~> m2 or kg2 m-4].
  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].
  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
  integer :: i, k, is, ie, nz
  is = g%isc ; ie = g%iec ; nz = g%ke

!  max_ent = 1.0e14*GV%Angstrom_H ! This is set to avoid roundoff problems.
  max_ent = 1.0e4*gv%m_to_H
  h_neglect = gv%H_subroundoff

  do i=is,ie ; pres(i) = tv%P_Ref ; enddo
  eosdom(:) = eos_domain(g%HI)
  do k=1,kmb
    call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pres, rcv, tv%eqn_of_state, eosdom)
    do i=is,ie
      h_bl(i,k) = h(i,j,k) + h_neglect
      sref(i,k) = rcv(i) - cs%Rho_sig_off
    enddo
  enddo

  do i=is,ie
    h_interior(i) = 0.0 ; ent_bl(i,1) = 0.0
!     if (kb(i) > nz) Ent_bl(i,Kmb+1) = 0.0
  enddo

  do k=2,kmb ; do i=is,ie
    if (do_i(i)) then
      ent_bl(i,k) = min(2.0 * dtkd_int(i,k) / (h(i,j,k-1) + h(i,j,k) + h_neglect), &
                        max_ent)
    else ; ent_bl(i,k) = 0.0 ; endif
  enddo ; enddo

  !   Determine the coordinate density of the bottommost buffer layer if there
  ! is no entrainment from the layers below.  This is a partial solver, based
  ! on the first pass of a tridiagonal solver, as the values in the upper buffer
  ! layers are not needed.

  do i=is,ie
    b1(i) = 1.0 / (h_bl(i,1) + ent_bl(i,2))
    d1(i) = h_bl(i,1) * b1(i)  ! = 1.0 - Ent_bl(i,2)*b1(i)
    s_est(i,1) = (h_bl(i,1)*sref(i,1)) * b1(i)
  enddo
  do k=2,kmb-1 ; do i=is,ie
    b1(i) = 1.0 / ((h_bl(i,k) + ent_bl(i,k+1)) + d1(i)*ent_bl(i,k))
    d1(i) = (h_bl(i,k) + d1(i)*ent_bl(i,k)) * b1(i)  ! = 1.0 - Ent_bl(i,K+1)*b1(i)
    s_est(i,k) = (h_bl(i,k)*sref(i,k) + ent_bl(i,k)*s_est(i,k-1)) * b1(i)
  enddo ; enddo
  do i=is,ie
    s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &
                   (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))
    frac_rem(i) = 1.0
  enddo

  !   Entrain any thin interior layers that are lighter (in the coordinate
  ! potential density) than the deepest buffer layer will be, and adjust kb.
  do i=is,ie ; kb(i) = nz+1 ; if (do_i(i)) kb(i) = kmb+1 ; enddo

  do k=kmb+1,nz ; do i=is,ie ; if (do_i(i)) then
    if ((k == kb(i)) .and. (s_est(i,kmb) > (gv%Rlay(k) - cs%Rho_sig_off))) then
      if (4.0*dtkd_int(i,kmb+1)*frac_rem(i) > &
          (h_bl(i,kmb) + h(i,j,k)) * (h(i,j,k) - gv%Angstrom_H)) then
        ! Entrain this layer into the buffer layer and move kb down.
        dh = max((h(i,j,k) - gv%Angstrom_H), 0.0)
        if (dh > 0.0) then
          frac_rem(i) = frac_rem(i) - ((h_bl(i,kmb) + h(i,j,k)) * dh) / &
                                       (4.0*dtkd_int(i,kmb+1))
          sref(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + dh*(gv%Rlay(k)-cs%Rho_sig_off)) / &
                        (h_bl(i,kmb) + dh)
          h_bl(i,kmb) = h_bl(i,kmb) + dh
          s_est(i,kmb) = (h_bl(i,kmb)*sref(i,kmb) + ent_bl(i,kmb)*s_est(i,kmb-1)) / &
                         (h_bl(i,kmb) + d1(i)*ent_bl(i,kmb))
        endif
        kb(i) = kb(i) + 1
      endif
    endif
  endif ; enddo ; enddo

 !    This is where variables are be set up with a different vertical grid
 !  in which the (newly?) massless layers are taken out.
  do k=nz,kmb+1,-1 ; do i=is,ie
    if (k >= kb(i)) h_interior(i) = h_interior(i) + (h(i,j,k)-gv%Angstrom_H)
    if (k==kb(i)) then
      h_bl(i,kmb+1) = h(i,j,k) ; sref(i,kmb+1) = gv%Rlay(k) - cs%Rho_sig_off
    elseif (k==kb(i)+1) then
      h_bl(i,kmb+2) = h(i,j,k) ; sref(i,kmb+2) = gv%Rlay(k) - cs%Rho_sig_off
    endif
  enddo ; enddo
  do i=is,ie ; if (kb(i) >= nz) then
    h_bl(i,kmb+1) = h(i,j,nz)
    sref(i,kmb+1) = gv%Rlay(nz) - cs%Rho_sig_off
    h_bl(i,kmb+2) = gv%Angstrom_H
    sref(i,kmb+2) = sref(i,kmb+1) + (gv%Rlay(nz) - gv%Rlay(nz-1))
  endif ; enddo

  !   Perhaps we should revisit the way that the average entrainment between the
  ! buffer layer and the interior is calculated so that it is not unduly
  ! limited when the layers are less than sqrt(Kd * dt) thick?
  do i=is,ie ; if (do_i(i)) then
    kd_x_dt = frac_rem(i) * dtkd_int(i,kmb+1)
    if ((kd_x_dt <= 0.0) .or. (h_interior(i) <= 0.0)) then
      ent_bl(i,kmb+1) = 0.0
    else
      !   If the combined layers are exceptionally thin, use sqrt(Kd*dt) as the
      ! estimate of the thickness in the denominator of the thickness diffusion.
      ent_bl(i,kmb+1) = min(0.5*h_interior(i), sqrt(kd_x_dt), &
                            kd_x_dt / (0.5*(h_bl(i,kmb) + h_bl(i,kmb+1))))
    endif
  else
    ent_bl(i,kmb+1) = 0.0
  endif ; enddo

end subroutine set_ent_bl

!>   This subroutine determines the reference density difference between the
!! bottommost buffer layer and the first interior after the mixing between mixed
!! and buffer layers and mixing with the layer below. Within the mixed and buffer
!! layers, entrainment from the layer above is increased when it is necessary to
!! keep the layers from developing a negative thickness; otherwise it equals
!! Ent_bl.  At each interface, the upward and downward fluxes average out to
!! Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl.
!!   The density difference across the first interior layer may also be returned.
!! It could also be limited to avoid negative values or values that greatly
!! exceed the density differences across an interface.
!!   Additionally, the partial derivatives of dSkb and dSlay with E_kb could
!! also be returned.
subroutine determine_dskb(h_bl, Sref, Ent_bl, E_kb, is, ie, kmb, G, GV, limit, &
                          dSkb, ddSkb_dE, dSlay, ddSlay_dE, dS_anom_lim, do_i_in)
  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_(G)),   intent(in)    :: h_bl   !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZK_(G)),   intent(in)    :: Sref   !< Reference potential density [R ~> kg m-3]
  real, dimension(SZI_(G),SZK_(G)),   intent(in)    :: Ent_bl !< The average entrainment upward and
                                                              !! downward across each interface
                                                              !! around the buffer layers [H ~> m or kg m-2].
  real, dimension(SZI_(G)),           intent(in)    :: E_kb   !< The entrainment by the top interior
                                                              !! layer [H ~> m or kg m-2].
  integer,                            intent(in)    :: is     !< The start of the i-index range to work on.
  integer,                            intent(in)    :: ie     !< The end of the i-index range to work on.
  integer,                            intent(in)    :: kmb    !< The number of mixed and buffer layers.
  logical,                            intent(in)    :: limit  !< If true, limit dSkb and dSlay to
                                                              !! avoid negative values.
  real, dimension(SZI_(G)),           intent(inout) :: dSkb   !< The limited potential density
                                                              !! difference across the interface
                                                              !! between the bottommost buffer layer
                                                              !! and the topmost interior layer. [R ~> kg m-3]
                                                              !! dSkb > 0.
  real, dimension(SZI_(G)), optional, intent(inout) :: ddSkb_dE !< The partial derivative of dSkb
                                                              !! with E [R H-1 ~> kg m-4 or m-1].
  real, dimension(SZI_(G)), optional, intent(inout) :: dSlay  !< The limited potential density
                                                              !! difference across the topmost
                                                              !! interior layer. 0 < dSkb [R ~> kg m-3]
  real, dimension(SZI_(G)), optional, intent(inout) :: ddSlay_dE !< The partial derivative of dSlay
                                                              !! with E [R H-1 ~> kg m-4 or m-1].
  real, dimension(SZI_(G)), optional, intent(inout) :: dS_anom_lim !< A limiting value to use for
                                                              !! the density anomalies below the
                                                              !! buffer layer [R ~> kg m-3].
  logical, dimension(SZI_(G)), optional, intent(in) :: do_i_in !< If present, determines which
                                                              !! columns are worked on.

! Note that dSkb, ddSkb_dE, dSlay, ddSlay_dE, and dS_anom_lim are declared
! intent inout  because they should not change where do_i_in is false.

!   This subroutine determines the reference density difference between the
! bottommost buffer layer and the first interior after the mixing between mixed
! and buffer layers and mixing with the layer below. Within the mixed and buffer
! layers, entrainment from the layer above is increased when it is necessary to
! keep the layers from developing a negative thickness; otherwise it equals
! Ent_bl.  At each interface, the upward and downward fluxes average out to
! Ent_bl, unless entrainment by the layer below is larger than twice Ent_bl.
!   The density difference across the first interior layer may also be returned.
! It could also be limited to avoid negative values or values that greatly
! exceed the density differences across an interface.
!   Additionally, the partial derivatives of dSkb and dSlay with E_kb could
! also be returned.

  ! Local variables
  real, dimension(SZI_(G),SZK_(G)) :: &
    b1, c1, &       ! b1 and c1 are variables used by the tridiagonal solver.
    S, dS_dE, &     ! The coordinate density [R ~> kg m-3] and its derivative with E.
    ea, dea_dE, &   ! The entrainment from above and its derivative with E.
    eb, deb_dE      ! The entrainment from below and its derivative with E.
  real :: deriv_dSkb(szi_(g))
  real :: d1(szi_(g))  ! d1 = 1.0-c1 is also used by the tridiagonal solver.
  real :: src       ! A source term for dS_dR.
  real :: h1        ! The thickness in excess of the minimum that will remain
                    ! after exchange with the layer below [H ~> m or kg m-2].
  logical, dimension(SZI_(G)) :: do_i
  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 :: h_tr      ! h_tr is h at tracer points with a tiny thickness
                    ! added to ensure positive definiteness [H ~> m or kg m-2].
  real :: b_denom_1 ! The first term in the denominator of b1 [H ~> m or kg m-2].
  real :: rat
  real :: dS_kbp1, IdS_kbp1
  real :: deriv_dSLay
  real :: Inv_term     ! [nondim]
  real :: f1, df1_drat ! Temporary variables [nondim].
  real :: z, dz_drat, f2, df2_dz, expz ! Temporary variables [nondim].
  real :: eps_dSLay, eps_dSkb ! Small nondimensional constants.
  integer :: i, k

  if (present(ddslay_de) .and. .not.present(dslay)) call mom_error(fatal, &
      "In deterimine_dSkb, ddSLay_dE may only be present if dSlay is.")

  h_neglect = gv%H_subroundoff

  do i=is,ie
    ea(i,kmb+1) = e_kb(i) ; dea_de(i,kmb+1) = 1.0
    s(i,kmb+1) = sref(i,kmb+1) ; ds_de(i,kmb+1) = 0.0
    b1(i,kmb+1) = 0.0
    d1(i) = 1.0
    do_i(i) = .true.
  enddo
  if (present(do_i_in)) then
    do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
  endif
  do k=kmb,1,-1 ; do i=is,ie
    if (do_i(i)) then
      ! The do_i test here is only for efficiency.
    ! Determine the entrainment from below for each buffer layer.
      if (2.0*ent_bl(i,k+1) > ea(i,k+1)) then
        eb(i,k) = 2.0*ent_bl(i,k+1) - ea(i,k+1) ; deb_de(i,k) = -dea_de(i,k+1)
      else
        eb(i,k) = 0.0 ; deb_de(i,k) = 0.0
      endif

      ! Determine the entrainment from above for each buffer layer.
      h1 = (h_bl(i,k) - gv%Angstrom_H) + (eb(i,k) - ea(i,k+1))
      if (h1 >= 0.0) then
        ea(i,k) = ent_bl(i,k) ; dea_de(i,k) = 0.0
      elseif (ent_bl(i,k) + 0.5*h1 >= 0.0) then
        ea(i,k) = ent_bl(i,k) - 0.5*h1
        dea_de(i,k) = 0.5*(dea_de(i,k+1) - deb_de(i,k))
      else
        ea(i,k) = -h1
        dea_de(i,k) = dea_de(i,k+1) - deb_de(i,k)
      endif
    else
      ea(i,k) = 0.0 ; dea_de(i,k) = 0.0 ; eb(i,k) = 0.0 ; deb_de(i,k) = 0.0
    endif

    ! This is the first-pass of a tridiagonal solver for S.
    h_tr = h_bl(i,k) + h_neglect
    c1(i,k) = ea(i,k+1) * b1(i,k+1)
    b_denom_1 = (h_tr + d1(i)*eb(i,k))
    b1(i,k) = 1.0 / (b_denom_1 + ea(i,k))
    d1(i) = b_denom_1 * b1(i,k)

    s(i,k) = (h_tr*sref(i,k) + eb(i,k)*s(i,k+1)) * b1(i,k)
  enddo ; enddo
  do k=2,kmb ; do i=is,ie
    s(i,k) = s(i,k) + c1(i,k-1)*s(i,k-1)
  enddo ; enddo

  if (present(ddskb_de) .or. present(ddslay_de)) then
    ! These two tridiagonal solvers cannot be combined because the solutions for
    ! S are required as a source for dS_dE.
    do k=kmb,2,-1 ; do i=is,ie
      if (do_i(i) .and. (dea_de(i,k) - deb_de(i,k) > 0.0)) then
        src = (((s(i,k+1) - sref(i,k)) * (h_bl(i,k) + h_neglect) + &
                (s(i,k+1) - s(i,k-1)) * ea(i,k)) * deb_de(i,k) - &
               ((sref(i,k) - s(i,k-1)) * h_bl(i,k) + &
                (s(i,k+1) - s(i,k-1)) * eb(i,k)) * dea_de(i,k)) / &
              ((h_bl(i,k) + h_neglect + ea(i,k)) + eb(i,k))
      else ; src = 0.0 ; endif
      ds_de(i,k) = (src + eb(i,k)*ds_de(i,k+1)) * b1(i,k)
    enddo ; enddo
    do i=is,ie
      if (do_i(i) .and. (deb_de(i,1) < 0.0)) then
        src = (((s(i,2) - sref(i,1)) * (h_bl(i,1) + h_neglect)) * deb_de(i,1)) / &
              (h_bl(i,1) + h_neglect + eb(i,1))
      else ; src = 0.0 ; endif
      ds_de(i,1) = (src + eb(i,1)*ds_de(i,2)) * b1(i,1)
    enddo
    do k=2,kmb ; do i=is,ie
      ds_de(i,k) = ds_de(i,k) + c1(i,k-1)*ds_de(i,k-1)
    enddo ; enddo
  endif

  ! Now, apply any limiting and return the requested variables.

  eps_dskb = 1.0e-6   ! Should be a small, nondimensional, positive number.
  if (.not.limit) then
    do i=is,ie ; if (do_i(i)) then
      dskb(i) = sref(i,kmb+1) - s(i,kmb)
    endif ; enddo
    if (present(ddskb_de)) then ; do i=is,ie ; if (do_i(i)) then
      ddskb_de(i) = -1.0*ds_de(i,kmb)
    endif ; enddo ; endif

    if (present(dslay)) then ; do i=is,ie ; if (do_i(i)) then
      dslay(i) = 0.5 * (sref(i,kmb+2) - s(i,kmb))
    endif ; enddo ; endif
    if (present(ddslay_de)) then ; do i=is,ie ; if (do_i(i)) then
      ddslay_de(i) = -0.5*ds_de(i,kmb)
    endif ; enddo ; endif
  else
    do i=is,ie ; if (do_i(i)) then
      ! Need to ensure that 0 < dSkb <= S_kb - Sbl
      if (sref(i,kmb+1) - s(i,kmb) < eps_dskb*(sref(i,kmb+2) - sref(i,kmb+1))) then
        dskb(i) = eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) ; deriv_dskb(i) = 0.0
      else
        dskb(i) = sref(i,kmb+1) - s(i,kmb) ; deriv_dskb(i) = -1.0
      endif
      if (present(ddskb_de)) ddskb_de(i) = deriv_dskb(i)*ds_de(i,kmb)
    endif ; enddo

    if (present(dslay)) then
      dz_drat = 1000.0    ! The limit of large dz_drat the same as choosing a
                          ! Heaviside function.
      eps_dslay = 1.0e-10 ! Should be ~= GV%Angstrom_H / sqrt(Kd*dt)
      do i=is,ie ; if (do_i(i)) then
        ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)
        ids_kbp1 = 1.0 / (sref(i,kmb+2) - sref(i,kmb+1))
        rat = (sref(i,kmb+1) - s(i,kmb)) * ids_kbp1
        ! Need to ensure that 0 < dSLay <= 2*dSkb
        if (rat < 0.5) then
          ! The coefficients here are chosen so that at rat = 0.5, the value (1.5)
          ! and first derivative (-0.5) match with the "typical" case (next).
          ! The functional form here is arbitrary.
          !   f1 provides a reasonable profile that matches the value and derivative
          ! of the "typical" case at rat = 0.5, and has a maximum of less than 2.
          inv_term = 1.0 / (1.0-rat)
          f1 = 2.0 - 0.125*(inv_term**2)
          df1_drat = - 0.25*(inv_term**3)

          !   f2 ensures that dSLay goes to 0 rapidly if rat is significantly
          ! negative.
          z = dz_drat * rat + 4.0 ! The 4 here gives f2(0) = 0.982.
          if (z >= 18.0) then ; f2 = 1.0 ; df2_dz = 0.0
          elseif (z <= -58.0) then ; f2 = eps_dslay ; df2_dz = 0.0
          else
            expz = exp(z) ; inv_term = 1.0 / (1.0 + expz)
            f2 = (eps_dslay + expz) * inv_term
            df2_dz = (1.0 - eps_dslay) * expz * inv_term**2
          endif

          dslay(i) = dskb(i) * f1 * f2
          deriv_dslay = deriv_dskb(i) * (f1 * f2) - (dskb(i)*ids_kbp1) * &
                            (df1_drat*f2 + f1 * dz_drat * df2_dz)
        elseif (dskb(i) <= 3.0*ds_kbp1) then
          ! This is the "typical" case.
          dslay(i) = 0.5 * (dskb(i) + ds_kbp1)
          deriv_dslay = 0.5 * deriv_dskb(i) ! = -0.5
        else
          dslay(i) = 2.0*ds_kbp1
          deriv_dslay = 0.0
        endif
        if (present(ddslay_de)) ddslay_de(i) = deriv_dslay*ds_de(i,kmb)
      endif ; enddo
    endif ! present(dSlay)
  endif ! Not limited.

  if (present(ds_anom_lim)) then ; do i=is,ie ; if (do_i(i)) then
    ds_anom_lim(i) = max(0.0, eps_dskb * (sref(i,kmb+2) - sref(i,kmb+1)) - &
                              (sref(i,kmb+1) - s(i,kmb)) )
  endif ; enddo ; endif

end subroutine determine_dskb

!>   Given an entrainment from below for layer kb, determine a consistent
!! entrainment from above, such that dSkb * ea_kb = dSkbp1 * F_kb.  The input
!! value of ea_kb is both the maximum value that can be obtained and the first
!! guess of the iterations.  Ideally ea_kb should be an under-estimate
subroutine f_kb_to_ea_kb(h_bl, Sref, Ent_bl, I_dSkbp1, F_kb, kmb, i, &
                         G, GV, CS, ea_kb, tol_in)
  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_(G)), &
                            intent(in)    :: h_bl !< Layer thickness, with the top interior
                                                  !! layer at k-index kmb+1 [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)), &
                            intent(in)    :: Sref !< The coordinate reference potential density,
                                                  !! with the value of the topmost interior layer
                                                  !! at index kmb+1 [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(G)), &
                            intent(in)    :: Ent_bl !< The average entrainment upward and downward
                                                  !! across each interface around the buffer layers,
                                                  !! [H ~> m or kg m-2].
  real, dimension(SZI_(G)), intent(in)    :: I_dSkbp1 !< The inverse of the difference in reference
                                                  !! potential density across the base of the
                                                  !! uppermost interior layer [R-1 ~> m3 kg-1].
  real, dimension(SZI_(G)), intent(in)    :: F_kb !< The entrainment from below by the
                                                  !! uppermost interior layer [H ~> m or kg m-2]
  integer,                  intent(in)    :: kmb  !< The number of mixed and buffer layers.
  integer,                  intent(in)    :: i    !< The i-index to work on
  type(entrain_diffusive_cs), pointer     :: CS   !< This module's control structure.
  real, dimension(SZI_(G)), intent(inout) :: ea_kb !< The entrainment from above by the layer below
                                                  !! the buffer layer (i.e. layer kb) [H ~> m or kg m-2].
  real,           optional, intent(in)    :: tol_in !< A tolerance for the iterative determination
                                                  !! of the entrainment [H ~> m or kg m-2].

  real :: max_ea, min_ea
  real :: err, err_min, err_max
  real :: derr_dea
  real :: val, tolerance, tol1
  real :: ea_prev
  real :: dS_kbp1
  logical :: bisect_next, Newton
  real, dimension(SZI_(G)) :: dS_kb
  real, dimension(SZI_(G)) :: maxF, ent_maxF, zeros
  real, dimension(SZI_(G)) :: ddSkb_dE
  integer :: it
  integer, parameter :: MAXIT = 30

  ds_kbp1 = sref(i,kmb+2) - sref(i,kmb+1)
  max_ea = ea_kb(i) ; min_ea = 0.0
  val = ds_kbp1 * f_kb(i)
  err_min = -val

  tolerance = cs%Tolerance_Ent
  if (present(tol_in)) tolerance = tol_in
  bisect_next = .true.

  call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &
                      ds_kb, ddskb_de)

  err = ds_kb(i) * ea_kb(i) - val
  derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)
  ! Return if Newton's method on the first guess would give a tolerably small
  ! change in the value of ea_kb.
  if ((err <= 0.0) .and. (abs(err) <= tolerance*abs(derr_dea))) return

  if (err == 0.0) then ; return ! The exact solution on the first guess...
  elseif (err > 0.0) then ! The root is properly bracketed.
    max_ea = ea_kb(i) ; err_max = err
    !   Use Newton's method (if it stays bounded) or the false position method
    ! to find the next value.
    if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i) - min_ea) > err) .and. &
        (derr_dea*(max_ea - ea_kb(i)) > -1.0*err)) then
      ea_kb(i) = ea_kb(i) - err / derr_dea
    else ! Use the bisection for the next guess.
      ea_kb(i) = 0.5*(max_ea+min_ea)
    endif
  else
    !   Try to bracket the root first.  If unable to bracket the root, return
    ! the maximum.
    zeros(i) = 0.0
    call find_maxf_kb(h_bl, sref, ent_bl, i_dskbp1, zeros, ea_kb, &
                      kmb, i, i, g, gv, cs, maxf, ent_maxf, f_thresh = f_kb)
    err_max = ds_kbp1 * maxf(i) - val
    ! If err_max is negative, there is no good solution, so use the maximum
    ! value of F in the valid range.
    if (err_max <= 0.0) then
      ea_kb(i) = ent_maxf(i) ; return
    else
      max_ea = ent_maxf(i)
      ea_kb(i) = 0.5*(max_ea+min_ea) ! Use bisection for the next guess.
    endif
  endif

  ! Exit if the range between max_ea and min_ea already acceptable.
  ! if (abs(max_ea - min_ea) < 0.1*tolerance) return

  do it = 1, maxit
    call determine_dskb(h_bl, sref, ent_bl, ea_kb, i, i, kmb, g, gv, .true., &
                        ds_kb, ddskb_de)

    err = ds_kb(i) * ea_kb(i) - val
    derr_dea = ds_kb(i) + ddskb_de(i) * ea_kb(i)

    ea_prev = ea_kb(i)
    ! Use Newton's method or the false position method to find the next value.
    newton = .false.
    if (err > 0.0) then
      max_ea = ea_kb(i) ; err_max = err
      if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-min_ea) > err)) newton = .true.
    else
      min_ea = ea_kb(i) ; err_min = err
      if ((derr_dea > 0.0) .and. (derr_dea*(ea_kb(i)-max_ea) < err)) newton = .true.
    endif

    if (newton) then
      ea_kb(i) = ea_kb(i) - err / derr_dea
    elseif (bisect_next) then ! Use bisection to reduce the range.
      ea_kb(i) = 0.5*(max_ea+min_ea)
      bisect_next = .false.
    else  ! Use the false-position method for the next guess.
      ea_kb(i) = min_ea + (max_ea-min_ea) * (err_min/(err_min - err_max))
      bisect_next = .true.
    endif

    tol1 = tolerance ; if (err > 0.0) tol1 = 0.099*tolerance
    if (ds_kb(i) <= ds_kbp1) then
      if (abs(ea_kb(i) - ea_prev) <= tol1) return
    else
      if (ds_kbp1*abs(ea_kb(i) - ea_prev) <= ds_kb(i)*tol1) return
    endif
  enddo

end subroutine f_kb_to_ea_kb


!>  This subroutine determines the entrainment from above by the top interior
!! layer (labeled kb elsewhere) given an entrainment by the layer below it,
!! constrained to be within the provided bounds.
subroutine determine_ea_kb(h_bl, dtKd_kb, Sref, I_dSkbp1, Ent_bl, ea_kbp1, &
                           min_eakb, max_eakb, kmb, is, ie, do_i, G, GV, CS, Ent, &
                           error, err_min_eakb0, err_max_eakb0, F_kb, dFdfm_kb)
  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_(G)), intent(in)  :: h_bl     !< Layer thickness, with the top interior
                                                            !! layer at k-index kmb+1 [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)), intent(in)  :: Sref     !< The coordinate reference potential
                                                            !! density, with the value of the
                                                            !! topmost interior layer at layer
                                                            !! kmb+1 [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(G)), intent(in)  :: Ent_bl   !< The average entrainment upward and
                                                            !! downward across each interface around
                                                            !! the buffer layers [H ~> m or kg m-2].
  real, dimension(SZI_(G)),         intent(in)  :: I_dSkbp1 !< The inverse of the difference in
                                                            !! reference potential density across
                                                            !! the base of the uppermost interior
                                                            !! layer [R-1 ~> m3 kg-1].
  real, dimension(SZI_(G)),         intent(in)  :: dtKd_kb  !< The diapycnal diffusivity in the top
                                                            !! interior layer times the time step
                                                            !! [H2 ~> m2 or kg2 m-4].
  real, dimension(SZI_(G)),         intent(in)  :: ea_kbp1  !< The entrainment from above by layer
                                                            !! kb+1 [H ~> m or kg m-2].
  real, dimension(SZI_(G)),         intent(in)  :: min_eakb !< The minimum permissible rate of
                                                            !! entrainment [H ~> m or kg m-2].
  real, dimension(SZI_(G)),         intent(in)  :: max_eakb !< The maximum permissible rate of
                                                            !! entrainment [H ~> m or kg m-2].
  integer,                          intent(in)  :: kmb      !< The number of mixed and buffer layers.
  integer,                          intent(in)  :: is       !< The start of the i-index range to work on.
  integer,                          intent(in)  :: ie       !< The end of the i-index range to work on.
  logical, dimension(SZI_(G)),      intent(in)  :: do_i     !< A logical variable indicating which
                                                            !! i-points to work on.
  type(entrain_diffusive_cs),       pointer     :: CS       !< This module's control structure.
  real, dimension(SZI_(G)),         intent(inout) :: Ent    !< The entrainment rate of the uppermost
                                                            !! interior layer [H ~> m or kg m-2].
                                                            !! The input value is the first guess.
  real, dimension(SZI_(G)), optional, intent(out) :: error  !< The error (locally defined in this
                                                            !! routine) associated with the returned
                                                            !! solution.
  real, dimension(SZI_(G)), optional, intent(in)  :: err_min_eakb0 !< The errors (locally defined)
                                                            !! associated with min_eakb when ea_kbp1 = 0,
                                                            !! returned from a previous call to this fn.
  real, dimension(SZI_(G)), optional, intent(in)  :: err_max_eakb0 !< The errors (locally defined)
                                                            !! associated with min_eakb when ea_kbp1 = 0,
                                                            !! returned from a previous call to this fn.
  real, dimension(SZI_(G)), optional, intent(out) :: F_kb   !< The entrainment from below by the
                                                            !! uppermost interior layer
                                                            !! corresponding to the returned
                                                            !! value of Ent [H ~> m or kg m-2].
  real, dimension(SZI_(G)), optional, intent(out) :: dFdfm_kb !< The partial derivative of F_kb with
                                                            !! ea_kbp1 [nondim].

!  This subroutine determines the entrainment from above by the top interior
! layer (labeled kb elsewhere) given an entrainment by the layer below it,
! constrained to be within the provided bounds.

  ! Local variables
  real, dimension(SZI_(G)) :: &
    dS_kb, &                !   The coordinate-density difference between the
                            ! layer kb and deepest buffer layer, limited to
                            ! ensure that it is positive [R ~> kg m-3].
    ds_lay, &               !   The coordinate-density difference across layer
                            ! kb, limited to ensure that it is positive and not
                            ! too much bigger than dS_kb or dS_kbp1 [R ~> kg m-3].
    ddskb_de, ddslay_de, &  ! The derivatives of dS_kb and dS_Lay with E
                            ! [R H-1 ~> kg m-4 or m-1].
    derror_de, &            ! The derivative of err with E [H ~> m or kg m-2].
    err, &                  ! The "error" whose zero is being sought [H2 ~> m2 or kg2 m-4].
    e_min, e_max, &         ! The minimum and maximum values of E [H ~> m or kg m-2].
    error_mine, error_maxe  ! err when E = E_min or E = E_max [H2 ~> m2 or kg2 m-4].
  real :: err_est           ! An estimate of what err will be [H2 ~> m2 or kg2 m-4].
  real :: eL                ! 1 or 0, depending on whether increases in E lead
                            ! to decreases in the entrainment from below by the
                            ! deepest buffer layer [nondim].
  real :: fa                ! Temporary variable used to calculate err [nondim].
  real :: fk                ! Temporary variable used to calculate err [H2 ~> m2 or kg2 m-4].
  real :: fm, fr            ! Temporary variables used to calculate err [H ~> m or kg m-2].
  real :: tolerance         ! The tolerance within which E must be converged [H ~> m or kg m-2].
  real :: E_prev            ! The previous value of E [H ~> m or kg m-2].
  logical, dimension(SZI_(G)) :: false_position ! If true, the false position
                            ! method might be used for the next iteration.
  logical, dimension(SZI_(G)) :: redo_i ! If true, more work is needed on this column.
  logical :: do_any
  real :: large_err         ! A large error measure [H2 ~> m2 or kg2 m-4].
  integer :: i, it
  integer, parameter :: MAXIT = 30

  if (.not.cs%bulkmixedlayer) then
    call mom_error(fatal, "determine_Ea_kb should not be called "//&
                           "unless BULKMIXEDLAYER is defined.")
  endif
  tolerance = cs%Tolerance_Ent
  large_err = gv%m_to_H**2 * 1.0e30

  do i=is,ie ; redo_i(i) = do_i(i) ; enddo

  do i=is,ie ; if (do_i(i)) then
    ! The first guess of Ent was the value from the previous iteration.

    !   These were previously calculated and provide good limits and estimates
    ! of the errors there. By construction the errors increase with R*ea_kbp1.
    e_min(i) = min_eakb(i) ; e_max(i) = max_eakb(i)
    error_mine(i) = -large_err ; error_maxe(i) = large_err
    false_position(i) = .true. ! Used to alternate between false_position and
                               ! bisection when Newton's method isn't working.
    if (present(err_min_eakb0)) error_mine(i) = err_min_eakb0(i) - e_min(i) * ea_kbp1(i)
    if (present(err_max_eakb0)) error_maxe(i) = err_max_eakb0(i) - e_max(i) * ea_kbp1(i)

    if ((error_maxe(i) <= 0.0) .or. (error_mine(i) >= 0.0)) then
      ! The root is not bracketed and one of the limiting values should be used.
      if (error_maxe(i) <= 0.0) then
        ! The errors decrease with E*ea_kbp1, so E_max is the best solution.
        ent(i) = e_max(i) ; err(i) = error_maxe(i)
      else  ! error_minE >= 0 is equivalent to ea_kbp1 = 0.0.
        ent(i) = e_min(i) ; err(i) = error_mine(i)
      endif
      derror_de(i) = 0.0
      redo_i(i) = .false.
    endif
  endif ; enddo   ! End of i-loop

  do it = 1,maxit
    do_any = .false. ; do i=is,ie ; if (redo_i(i)) do_any = .true. ; enddo
    if (.not.do_any) exit
    call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., ds_kb, &
                        ddskb_de, ds_lay, ddslay_de, do_i_in = redo_i)
    do i=is,ie ; if (redo_i(i)) then
      !  The correct root is bracketed between E_min and E_max.
      ! Note the following limits:  Ent >= 0 ; fa > 1 ; fk > 0
      el = 0.0 ; if (2.0*ent_bl(i,kmb+1) >= ent(i)) el = 1.0
      fa = (1.0 + el) + ds_kb(i)*i_dskbp1(i)
      fk = dtkd_kb(i) * (ds_lay(i)/ds_kb(i))
      fm = (ea_kbp1(i) - h_bl(i,kmb+1)) + el*2.0*ent_bl(i,kmb+1)
      if (fm > -gv%Angstrom_H) fm = fm + gv%Angstrom_H  ! This could be smooth if need be.
      err(i) = (fa * ent(i)**2 - fm * ent(i)) - fk
      derror_de(i) = ((2.0*fa + (ddskb_de(i)*i_dskbp1(i))*ent(i))*ent(i) - fm) - &
          dtkd_kb(i) * (ddslay_de(i)*ds_kb(i) - ddskb_de(i)*ds_lay(i))/(ds_kb(i)**2)

      if (err(i) == 0.0) then
        redo_i(i) = .false. ; cycle
      elseif (err(i) > 0.0) then
        e_max(i) = ent(i) ; error_maxe(i) = err(i)
      else
        e_min(i) = ent(i) ; error_mine(i) = err(i)
      endif

      e_prev = ent(i)
      if ((it == 1) .or. (derror_de(i) <= 0.0)) then
        !   Assuming that the coefficients of the quadratic equation are correct
        ! will usually give a very good first guess.  Also, if derror_dE < 0.0,
        ! R is on the wrong side of the approximate parabola.  In either case,
        ! try assuming that the error is approximately a parabola and solve.
        fr = sqrt(fm**2 + 4.0*fa*fk)
        if (fm >= 0.0) then
          ent(i) = (fm + fr) / (2.0 * fa)
        else
          ent(i) = (2.0 * fk) / (fr - fm)
        endif
        ! But make sure that the root stays bracketed, bisecting if needed.
        if ((ent(i) > e_max(i)) .or. (ent(i) < e_min(i))) &
          ent(i) = 0.5*(e_max(i) + e_min(i))
      elseif (((e_max(i)-ent(i))*derror_de(i) > -err(i)) .and. &
              ((ent(i)-e_min(i))*derror_de(i) > err(i)) ) then
        ! Use Newton's method for the next estimate, provided it will
        ! remain bracketed between Rmin and Rmax.
        ent(i) = ent(i) - err(i) / derror_de(i)
      elseif (false_position(i) .and. &
              (error_maxe(i) - error_mine(i) < 0.9*large_err)) then
        ! Use the false postion method if there are decent error estimates.
        ent(i) = e_min(i) + (e_max(i)-e_min(i)) * &
                (-error_mine(i)/(error_maxe(i) - error_mine(i)))
        false_position(i) = .false.
      else ! Bisect as a last resort or if the false position method was used last.
        ent(i) = 0.5*(e_max(i) + e_min(i))
        false_position(i) = .true.
      endif

      if (abs(e_prev - ent(i)) < tolerance) then
        err_est = err(i) + (ent(i) - e_prev) * derror_de(i)
        if ((it > 1) .or. (err_est*err(i) <= 0.0) .or. &
            (abs(err_est) < abs(tolerance*derror_de(i)))) redo_i(i) = .false.
      endif

    endif ; enddo   ! End of i-loop
  enddo ! End of iterations to determine Ent(i).

  ! Update the value of dS_kb for consistency with Ent.
  if (present(f_kb) .or. present(dfdfm_kb)) &
    call determine_dskb(h_bl, sref, ent_bl, ent, is, ie, kmb, g, gv, .true., &
                        ds_kb, do_i_in = do_i)

  if (present(f_kb)) then ; do i=is,ie ; if (do_i(i)) then
    f_kb(i) = ent(i) * (ds_kb(i) * i_dskbp1(i))
  endif ; enddo ; endif
  if (present(error)) then ; do i=is,ie ; if (do_i(i)) then
    error(i) = err(i)
  endif ; enddo ; endif
  if (present(dfdfm_kb)) then ; do i=is,ie ; if (do_i(i)) then
    !   derror_dE and ddSkb_dE are _not_ recalculated here, since dFdfm_kb is
    ! only used in Newton's method, and slightly increasing the accuracy of the
    ! estimate is unlikely to speed convergence.
    if (derror_de(i) > 0.0) then
      dfdfm_kb(i) = ((ds_kb(i) + ent(i) * ddskb_de(i)) * i_dskbp1(i)) * &
                    (ent(i) / derror_de(i))
    else ! Use Adcroft's division by 0 convention.
      dfdfm_kb(i) = 0.0
    endif
  endif ; enddo ; endif

end subroutine determine_ea_kb

!> Maximize F = ent*ds_kb*I_dSkbp1 in the range min_ent < ent < max_ent.
subroutine find_maxf_kb(h_bl, Sref, Ent_bl, I_dSkbp1, min_ent_in, max_ent_in, &
                        kmb, is, ie, G, GV, CS, maxF, ent_maxF, do_i_in, &
                        F_lim_maxent, F_thresh)
  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_(G)), &
                              intent(in)  :: h_bl     !< Layer thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZK_(G)), &
                              intent(in)  :: Sref     !< Reference potential density [R ~> kg m-3].
  real, dimension(SZI_(G),SZK_(G)), &
                              intent(in)  :: Ent_bl   !< The average entrainment upward and
                                                      !! downward across each interface around
                                                      !! the buffer layers [H ~> m or kg m-2].
  real, dimension(SZI_(G)),   intent(in)  :: I_dSkbp1 !< The inverse of the difference in
                                                      !! reference potential density across the
                                                      !! base of the uppermost interior layer
                                                      !! [R-1 ~> m3 kg-1].
  real, dimension(SZI_(G)),   intent(in)  :: min_ent_in !< The minimum value of ent to search,
                                                      !! [H ~> m or kg m-2].
  real, dimension(SZI_(G)),   intent(in)  :: max_ent_in !< The maximum value of ent to search,
                                                      !! [H ~> m or kg m-2].
  integer,                    intent(in)  :: kmb      !< The number of mixed and buffer layers.
  integer,                    intent(in)  :: is       !< The start of the i-index range to work on.
  integer,                    intent(in)  :: ie       !< The end of the i-index range to work on.
  type(entrain_diffusive_cs), pointer     :: CS       !< This module's control structure.
  real, dimension(SZI_(G)),   intent(out) :: maxF     !< The maximum value of F
                                                      !! = ent*ds_kb*I_dSkbp1 found in the range
                                                      !! min_ent < ent < max_ent [H ~> m or kg m-2].
  real, dimension(SZI_(G)), &
                    optional, intent(out) :: ent_maxF !< The value of ent at that maximum [H ~> m or kg m-2].
  logical, dimension(SZI_(G)), &
                    optional, intent(in)  :: do_i_in  !< A logical array indicating which columns
                                                      !! to work on.
  real, dimension(SZI_(G)), &
                    optional, intent(out) :: F_lim_maxent !< If present, do not apply the limit in
                                                      !! finding the maximum value, but return the
                                                      !! limited value at ent=max_ent_in in this
                                                      !! array [H ~> m or kg m-2].
  real, dimension(SZI_(G)), &
                    optional, intent(in)  :: F_thresh !< If F_thresh is present, return the first
                                                      !! value found that has F > F_thresh, or
                                                      !! the maximum.

! Maximize F = ent*ds_kb*I_dSkbp1 in the range min_ent < ent < max_ent.
! ds_kb may itself be limited to positive values in determine_dSkb, which gives
! the prospect of two local maxima in the range - one at max_ent_in with that
! minimum value of ds_kb, and the other due to the unlimited (potentially
! negative) value.  It is faster to find the true maximum by first finding the
! unlimited maximum and comparing it to the limited value at max_ent_in.
  real, dimension(SZI_(G)) :: &
    ent, &
    minent, maxent, ent_best, &
    F_max_ent_in, &
    F_maxent, F_minent, F, F_best, &
    dF_dent, dF_dE_max, dF_dE_min, dF_dE_best, &
    dS_kb, dS_kb_lim, ddSkb_dE, dS_anom_lim, &
    chg_prev, chg_pre_prev
  real :: dF_dE_mean, maxslope, minslope
  real :: tolerance
  real :: ratio_select_end
  real :: rat, max_chg, min_chg, chg1, chg2, chg
  logical, dimension(SZI_(G)) :: do_i, last_it, need_bracket, may_use_best
  logical :: doany, OK1, OK2, bisect, new_min_bound
  integer :: i, it, is1, ie1
  integer, parameter :: MAXIT = 20

  tolerance = cs%Tolerance_Ent

  if (present(do_i_in)) then
    do i=is,ie ; do_i(i) = do_i_in(i) ; enddo
  else
    do i=is,ie ; do_i(i) = .true. ; enddo
  endif

  ! The most likely value is at max_ent.
  call determine_dskb(h_bl, sref, ent_bl, max_ent_in, is, ie, kmb, g, gv, .false., &
                      ds_kb, ddskb_de, ds_anom_lim=ds_anom_lim)
  ie1 = is-1 ; doany = .false.
  do i=is,ie
    ds_kb_lim(i) = ds_kb(i) + ds_anom_lim(i)
    f_max_ent_in(i) = max_ent_in(i)*ds_kb_lim(i)*i_dskbp1(i)
    maxent(i) = max_ent_in(i) ; minent(i) = min_ent_in(i)
    if ((abs(maxent(i) - minent(i)) < tolerance) .or. (.not.do_i(i))) then
      f_best(i) = max_ent_in(i)*ds_kb(i)*i_dskbp1(i)
      ent_best(i) = max_ent_in(i) ; ent(i) = max_ent_in(i)
      do_i(i) = .false.
    else
      f_maxent(i) = maxent(i) * ds_kb(i) * i_dskbp1(i)
      df_de_max(i) = (ds_kb(i) + maxent(i)*ddskb_de(i)) * i_dskbp1(i)
      doany = .true. ; last_it(i) = .false. ; need_bracket(i) = .true.
    endif
  enddo

  if (doany) then
    ie1 = is-1 ; do i=is,ie ; if (do_i(i)) ie1 = i ; enddo
    do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo
    ! Find the value of F and its derivative at min_ent.
    call determine_dskb(h_bl, sref, ent_bl, minent, is1, ie1, kmb, g, gv, .false., &
                        ds_kb, ddskb_de, do_i_in = do_i)
    do i=is1,ie1 ; if (do_i(i)) then
      f_minent(i) = minent(i) * ds_kb(i) * i_dskbp1(i)
      df_de_min(i) = (ds_kb(i) + minent(i)*ddskb_de(i)) * i_dskbp1(i)
    endif ; enddo

    ratio_select_end = 0.9
    do it=1,maxit
      ratio_select_end = 0.5*ratio_select_end
      do i=is1,ie1 ; if (do_i(i)) then
        if (need_bracket(i)) then
          df_de_mean = (f_maxent(i) - f_minent(i)) / (maxent(i) - minent(i))
          maxslope = max(df_de_mean, df_de_min(i), df_de_max(i))
          minslope = min(df_de_mean, df_de_min(i), df_de_max(i))
          if (f_minent(i) >= f_maxent(i)) then
            if (df_de_min(i) > 0.0) then ; rat = 0.02 ! A small step should bracket the soln.
            elseif (maxslope < ratio_select_end*minslope) then
              ! The maximum of F is at minent.
              f_best(i) = f_minent(i) ; ent_best(i) = minent(i) ; rat = 0.0
              do_i(i) = .false.
            else ; rat = 0.382 ; endif ! Use the golden ratio
          else
            if (df_de_max(i) < 0.0) then ; rat = 0.98 ! A small step should bracket the soln.
            elseif (minslope > ratio_select_end*maxslope) then
              ! The maximum of F is at maxent.
              f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i) ; rat = 1.0
              do_i(i) = .false.
            else ; rat = 0.618 ; endif ! Use the golden ratio
          endif

          if (rat >= 0.0) ent(i) = rat*maxent(i) + (1.0-rat)*minent(i)
          if (((maxent(i) - minent(i)) < tolerance) .or. (it==maxit)) &
            last_it(i) = .true.
        else ! The maximum is bracketed by minent, ent_best, and maxent.
          chg1 = 2.0*(maxent(i) - minent(i)) ; chg2 = chg1
          if (df_de_best(i) > 0) then
            max_chg = maxent(i) - ent_best(i) ; min_chg = 0.0
          else
            max_chg = 0.0 ; min_chg = minent(i) - ent_best(i) ! < 0
          endif
          if (max_chg - min_chg < 2.0*tolerance) last_it(i) = .true.
          if (df_de_max(i) /= df_de_best(i)) &
            chg1 = (maxent(i) - ent_best(i))*df_de_best(i) / &
                   (df_de_best(i) - df_de_max(i))
          if (df_de_min(i) /= df_de_best(i)) &
            chg2 = (minent(i) - ent_best(i))*df_de_best(i) / &
                   (df_de_best(i) - df_de_min(i))
          ok1 = ((chg1 < max_chg) .and. (chg1 > min_chg))
          ok2 = ((chg2 < max_chg) .and. (chg2 > min_chg))
          if (.not.(ok1 .or. ok2)) then ; bisect = .true. ; else
            if (ok1 .and. ok2) then ! Take the acceptable smaller change.
              chg = chg1 ; if (abs(chg2) < abs(chg1)) chg = chg2
            elseif (ok1) then ; chg = chg1
            else ; chg = chg2 ; endif
            if (abs(chg) > 0.5*abs(chg_pre_prev(i))) then ; bisect = .true.
            else ; bisect = .false. ; endif
          endif
          chg_pre_prev(i) = chg_prev(i)
          if (bisect) then
            if (df_de_best(i) > 0.0) then
              ent(i) = 0.5*(maxent(i) + ent_best(i))
              chg_prev(i) = 0.5*(maxent(i) - ent_best(i))
            else
              ent(i) = 0.5*(minent(i) + ent_best(i))
              chg_prev(i) = 0.5*(minent(i) - ent_best(i))
            endif
          else
            if (abs(chg) < tolerance) chg = sign(tolerance,chg)
            ent(i) = ent_best(i) + chg
            chg_prev(i) = chg
          endif
        endif
      endif ; enddo

      if (mod(it,3) == 0) then  ! Re-determine the loop bounds.
        ie1 = is-1 ; do i=is1,ie ; if (do_i(i)) ie1 = i ; enddo
        do i=ie1,is,-1 ; if (do_i(i)) is1 = i ; enddo
      endif

      call determine_dskb(h_bl, sref, ent_bl, ent, is1, ie1, kmb, g, gv, .false., &
                          ds_kb, ddskb_de, do_i_in = do_i)
      do i=is1,ie1 ; if (do_i(i)) then
        f(i) = ent(i)*ds_kb(i)*i_dskbp1(i)
        df_dent(i) = (ds_kb(i) + ent(i)*ddskb_de(i)) * i_dskbp1(i)
      endif ; enddo

      if (present(f_thresh)) then ; do i=is1,ie1 ; if (do_i(i)) then
        if (f(i) >= f_thresh(i)) then
          f_best(i) = f(i) ; ent_best(i) = ent(i) ; do_i(i) = .false.
        endif
      endif ; enddo ; endif

      doany = .false.
      do i=is1,ie1 ; if (do_i(i)) then
        if (.not.last_it(i)) doany = .true.
        if (last_it(i)) then
          if (need_bracket(i)) then
            if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i))) then
              f_best(i) = f(i) ; ent_best(i) = ent(i)
            elseif (f_maxent(i) > f_minent(i)) then
              f_best(i) = f_maxent(i) ; ent_best(i) = maxent(i)
            else
              f_best(i) = f_minent(i) ; ent_best(i) = minent(i)
            endif
          elseif (f(i) > f_best(i)) then
            f_best(i) = f(i) ; ent_best(i) = ent(i)
          endif
          do_i(i) = .false.
        elseif (need_bracket(i)) then
          if ((f(i) > f_maxent(i)) .and. (f(i) > f_minent(i))) then
            need_bracket(i) = .false. ! The maximum is now bracketed.
            chg_prev(i) = (maxent(i) - minent(i))
            chg_pre_prev(i) = 2.0*chg_prev(i)
            ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)
          elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i))) then
            new_min_bound = .true.  ! We have a new minimum bound.
          elseif ((f(i) <= f_maxent(i)) .and. (f(i) > f_minent(i))) then
            new_min_bound = .false. ! We have a new maximum bound.
          else ! This case would bracket a minimum.  Wierd.
             ! Unless the derivative indicates that there is a maximum near the
             ! lower bound, try keeping the end with the larger value of F
             ! in a tie keep the minimum as the answer here will be compared
             ! with the maximum input value later.
             new_min_bound = .true.
             if (df_de_min(i) > 0.0 .or. (f_minent(i) >= f_maxent(i))) &
               new_min_bound = .false.
          endif
          if (need_bracket(i)) then ! Still not bracketed.
            if (new_min_bound) then
              minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)
            else
              maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)
            endif
          endif
        else  ! The root was previously bracketed.
          if (f(i) >= f_best(i)) then ! There is a new maximum.
            if (ent(i) > ent_best(i)) then ! Replace minent with ent_prev.
              minent(i) = ent_best(i) ; f_minent(i) = f_best(i) ; df_de_min(i) = df_de_best(i)
            else ! Replace maxent with ent_best.
              maxent(i) = ent_best(i) ; f_maxent(i) = f_best(i) ; df_de_max(i) = df_de_best(i)
            endif
            ent_best(i) = ent(i) ; f_best(i) = f(i) ; df_de_best(i) = df_dent(i)
          else
            if (ent(i) < ent_best(i)) then ! Replace the minent with ent.
              minent(i) = ent(i) ; f_minent(i) = f(i) ; df_de_min(i) = df_dent(i)
            else ! Replace maxent with ent_prev.
              maxent(i) = ent(i) ; f_maxent(i) = f(i) ; df_de_max(i) = df_dent(i)
            endif
          endif
          if ((maxent(i) - minent(i)) <= tolerance) do_i(i) = .false. ! Done.
        endif ! need_bracket.
      endif ; enddo
      if (.not.doany) exit
    enddo
  endif

  if (present(f_lim_maxent)) then
    ! Return the unlimited maximum in maxF, and the limited value of F at maxent.
    do i=is,ie
      maxf(i) = f_best(i)
      f_lim_maxent(i) = f_max_ent_in(i)
      if (present(ent_maxf)) ent_maxf(i) = ent_best(i)
    enddo
  else
    ! Now compare the two? potential maxima using the limited value of dF_kb.
    doany = .false.
    do i=is,ie
      may_use_best(i) = (ent_best(i) /= max_ent_in(i))
      if (may_use_best(i)) doany = .true.
    enddo
    if (doany) then
      ! For efficiency, could save previous value of dS_anom_lim_best?
      call determine_dskb(h_bl, sref, ent_bl, ent_best, is, ie, kmb, g, gv, .true., &
                          ds_kb_lim)
      do i=is,ie
        f_best(i) = ent_best(i)*ds_kb_lim(i)*i_dskbp1(i)
        ! The second test seems necessary because of roundoff differences that
        ! can arise during compilation.
        if ((f_best(i) > f_max_ent_in(i)) .and. (may_use_best(i))) then
          maxf(i) = f_best(i)
          if (present(ent_maxf)) ent_maxf(i) = ent_best(i)
        else
          maxf(i) = f_max_ent_in(i)
          if (present(ent_maxf)) ent_maxf(i) = max_ent_in(i)
        endif
      enddo
    else
      ! All of the maxima are at the maximum entrainment.
      do i=is,ie ; maxf(i) = f_max_ent_in(i) ; enddo
      if (present(ent_maxf)) then
        do i=is,ie ; ent_maxf(i) = max_ent_in(i) ; enddo
      endif
    endif
  endif

end subroutine find_maxf_kb

!> This subroutine initializes the parameters and memory associated with the
!! entrain_diffusive module.
subroutine entrain_diffusive_init(Time, G, GV, US, param_file, diag, CS, just_read_params)
  type(time_type),         intent(in)    :: Time !< The current model time.
  type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time
                                                 !! parameters.
  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
                                                 !! output.
  type(entrain_diffusive_cs), pointer    :: CS   !< A pointer that is set to point to the control
                                                 !! structure.
  logical,       optional, intent(in)    :: just_read_params !< If present and true, this call will
                                                 !! only read parameters logging them or registering
                                                 !! any diagnostics

  ! Local variables
  real :: dt  ! The dynamics timestep, used here in the default for TOLERANCE_ENT, in MKS units [s]
  real :: Kd  ! A diffusivity used in the default for TOLERANCE_ENT, in MKS units [m2 s-1]
  logical :: just_read    ! If true, just read parameters but do nothing else.
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "MOM_entrain_diffusive" ! This module's name.

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

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  cs%diag => diag

  cs%bulkmixedlayer = (gv%nkml > 0)

! Set default, read and log parameters
  if (.not.just_read) call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "MAX_ENT_IT", cs%max_ent_it, &
                 "The maximum number of iterations that may be used to "//&
                 "calculate the interior diapycnal entrainment.", default=5, do_not_log=just_read)
  ! In this module, KD is only used to set the default for TOLERANCE_ENT. [m2 s-1]
  call get_param(param_file, mdl, "KD", kd, default=0.0)
  call get_param(param_file, mdl, "DT", dt, &
                 "The (baroclinic) dynamics time step.", units = "s", &
                 fail_if_missing=.true., do_not_log=just_read)
  call get_param(param_file, mdl, "TOLERANCE_ENT", cs%Tolerance_Ent, &
                 "The tolerance with which to solve for entrainment values.", &
                 units="m", default=max(100.0*gv%Angstrom_m,1.0e-4*sqrt(dt*kd)), scale=gv%m_to_H, &
                 do_not_log=just_read)

  cs%Rho_sig_off = 1000.0*us%kg_m3_to_R

  if (.not.just_read) then
    cs%id_Kd = register_diag_field('ocean_model', 'Kd_effective', diag%axesTL, time, &
        'Diapycnal diffusivity as applied', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
    cs%id_diff_work = register_diag_field('ocean_model', 'diff_work', diag%axesTi, time, &
        'Work actually done by diapycnal diffusion across each interface', &
        'W m-2', conversion=us%RZ3_T3_to_W_m2)
  endif

  if (just_read) deallocate(cs)

end subroutine entrain_diffusive_init

!> This subroutine cleans up and deallocates any memory associated with the
!! entrain_diffusive module.
subroutine entrain_diffusive_end(CS)
  type(entrain_diffusive_cs), pointer :: CS !< A pointer to the control structure for this
                                            !! module that will be deallocated.
  if (associated(cs)) deallocate(cs)

end subroutine entrain_diffusive_end

!> \namespace mom_entrain_diffusive
!!
!! By Robert Hallberg, September 1997 - July 2000
!!
!!   This file contains the subroutines that implement diapycnal
!! mixing and advection in isopycnal layers.  The main subroutine,
!! calculate_entrainment, returns the entrainment by each layer
!! across the interfaces above and below it.  These are calculated
!! subject to the constraints that no layers can be driven to neg-
!! ative thickness and that the each layer maintains its target
!! density, using the scheme described in Hallberg (MWR 2000). There
!! may or may not be a bulk mixed layer above the isopycnal layers.
!! The solution is iterated until the change in the entrainment
!! between successive iterations is less than some small tolerance.
!!
!!   The dual-stream entrainment scheme of MacDougall and Dewar
!! (JPO 1997) is used for combined diapycnal advection and diffusion,
!! modified as described in Hallberg (MWR 2000) to be solved
!! implicitly in time.  Any profile of diffusivities may be used.
!! Diapycnal advection is fundamentally the residual of diapycnal
!! diffusion, so the fully implicit upwind differencing scheme that
!! is used is entirely appropriate.  The downward buoyancy flux in
!! each layer is determined from an implicit calculation based on
!! the previously calculated flux of the layer above and an estim-
!! ated flux in the layer below.  This flux is subject to the foll-
!! owing conditions:  (1) the flux in the top and bottom layers are
!! set by the boundary conditions, and (2) no layer may be driven
!! below an Angstrom thickness.  If there is a bulk mixed layer, the
!! mixed and buffer layers are treated as Eulerian layers, whose
!! thicknesses only change due to entrainment by the interior layers.

end module mom_entrain_diffusive