MOM_diabatic_aux.F90

!> Provides functions for some diabatic processes such as fraxil, brine rejection,
!! tendency due to surface flux divergence.
module mom_diabatic_aux

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

use mom_cpu_clock,     only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock,     only : clock_module_driver, clock_module, clock_routine
use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : diag_ctrl, time_type
use mom_eos,           only : calculate_density, calculate_tfreeze, eos_domain
use mom_eos,           only : calculate_specific_vol_derivs, calculate_density_derivs
use mom_error_handler, only : mom_error, fatal, warning, calltree_showquery
use mom_error_handler, only : calltree_enter, calltree_leave, calltree_waypoint
use mom_file_parser,   only : get_param, log_param, log_version, param_file_type
use mom_forcing_type,  only : forcing, extractfluxes1d, forcing_singlepointprint
use mom_grid,          only : ocean_grid_type
use mom_io,            only : slasher
use mom_opacity,       only : set_opacity, opacity_cs, extract_optics_slice, extract_optics_fields
use mom_opacity,       only : optics_type, optics_nbands, absorbremainingsw, sumswoverbands
use mom_tracer_flow_control, only : get_chl_from_model, tracer_flow_control_cs
use mom_unit_scaling,  only : unit_scale_type
use mom_variables,     only : thermo_var_ptrs ! , vertvisc_type, accel_diag_ptrs
use mom_verticalgrid,  only : verticalgrid_type
use time_interp_external_mod, only : init_external_field, time_interp_external
use time_interp_external_mod, only : time_interp_external_init

implicit none ; private

#include <MOM_memory.h>

public diabatic_aux_init, diabatic_aux_end
public make_frazil, adjust_salt, differential_diffuse_t_s, tridiagts
public find_uv_at_h, diagnosemldbydensitydifference, applyboundaryfluxesinout, set_pen_shortwave
public diagnosemldbyenergy

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

!> Control structure for diabatic_aux
type, public :: diabatic_aux_cs ; private
  logical :: do_rivermix = .false. !< Provide additional TKE to mix river runoff at the
                                   !! river mouths to a depth of "rivermix_depth"
  real    :: rivermix_depth = 0.0  !< The depth to which rivers are mixed if do_rivermix = T [Z ~> m].
  logical :: reclaim_frazil  !<   If true, try to use any frazil heat deficit to
                             !! to cool the topmost layer down to the freezing
                             !! point.  The default is false.
  logical :: pressure_dependent_frazil  !< If true, use a pressure dependent
                             !! freezing temperature when making frazil.  The
                             !! default is false, which will be faster but is
                             !! inappropriate with ice-shelf cavities.
  logical :: ignore_fluxes_over_land    !< If true, the model does not check
                             !! if fluxes are applied over land points. This
                             !! flag must be used when the ocean is coupled with
                             !! sea ice and ice shelves and use_ePBL = true.
  logical :: use_river_heat_content !< If true, assumes that ice-ocean boundary
                             !! has provided a river heat content. Otherwise, runoff
                             !! is added with a temperature of the local SST.
  logical :: use_calving_heat_content !< If true, assumes that ice-ocean boundary
                             !! has provided a calving heat content. Otherwise, calving
                             !! is added with a temperature of the local SST.
  logical :: var_pen_sw      !<   If true, use one of the CHL_A schemes to determine the
                             !! e-folding depth of incoming shortwave radiation.
  integer :: sbc_chl         !< An integer handle used in time interpolation of
                             !! chlorophyll read from a file.
  logical ::  chl_from_file  !< If true, chl_a is read from a file.

  type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
  type(diag_ctrl), pointer :: diag !< Structure used to regulate timing of diagnostic output

  ! Diagnostic handles
  integer :: id_createdh       = -1 !< Diagnostic ID of mass added to avoid grounding
  integer :: id_brine_lay      = -1 !< Diagnostic ID of which layer receives the brine
  integer :: id_pensw_diag     = -1 !< Diagnostic ID of Penetrative shortwave heating (flux convergence)
  integer :: id_penswflux_diag = -1 !< Diagnostic ID of Penetrative shortwave flux
  integer :: id_nonpensw_diag  = -1 !< Diagnostic ID of Non-penetrative shortwave heating
  integer :: id_chl            = -1 !< Diagnostic ID of chlorophyll-A handles for opacity

  ! Optional diagnostic arrays
  real, allocatable, dimension(:,:)   :: createdh       !< The amount of volume added in order to
                                                        !! avoid grounding [H T-1 ~> m s-1]
  real, allocatable, dimension(:,:,:) :: pensw_diag     !< Heating in a layer from convergence of
                                                        !! penetrative SW [Q R Z T-1 ~> W m-2]
  real, allocatable, dimension(:,:,:) :: penswflux_diag !< Penetrative SW flux at base of grid
                                                        !! layer [Q R Z T-1 ~> W m-2]
  real, allocatable, dimension(:,:)   :: nonpensw_diag  !< Non-downwelling SW radiation at ocean
                                                        !! surface [Q R Z T-1 ~> W m-2]

end type diabatic_aux_cs

!>@{ CPU time clock IDs
integer :: id_clock_uv_at_h, id_clock_frazil
!>@}

contains

!> Frazil formation keeps the temperature above the freezing point.
!! This subroutine warms any water that is colder than the (currently
!! surface) freezing point up to the freezing point and accumulates
!! the required heat (in [Q R Z ~> J m-2]) in tv%frazil.
subroutine make_frazil(h, tv, G, GV, US, CS, p_surf, halo)
  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]
  type(thermo_var_ptrs),   intent(inout) :: tv !< Structure containing pointers to any available
                                               !! thermodynamic fields.
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  type(diabatic_aux_cs),   intent(in)    :: cs !< The control structure returned by a previous
                                               !! call to diabatic_aux_init.
  real, dimension(SZI_(G),SZJ_(G)), &
                 optional, intent(in)    :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa].
  integer,       optional, intent(in)    :: halo !< Halo width over which to calculate frazil

  ! Local variables
  real, dimension(SZI_(G)) :: &
    fraz_col, & ! The accumulated heat requirement due to frazil [Q R Z ~> J m-2].
    t_freeze, & ! The freezing potential temperature at the current salinity [degC].
    ps          ! Surface pressure [R L2 T-2 ~> Pa]
  real, dimension(SZI_(G),SZK_(G)) :: &
    pressure    ! The pressure at the middle of each layer [R L2 T-2 ~> Pa].
  real :: h_to_rl2_t2  ! A conversion factor from thicknesses in H to pressure [R L2 T-2 H-1 ~> Pa m-1 or Pa m2 kg-1]
  real :: hc    ! A layer's heat capacity [Q R Z degC-1 ~> J m-2 degC-1].
  logical :: t_fr_set  ! True if the freezing point has been calculated for a
                       ! row of points.
  integer :: i, j, k, is, ie, js, je, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  if (present(halo)) then
    is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
  endif

  call cpu_clock_begin(id_clock_frazil)

  if (.not.cs%pressure_dependent_frazil) then
    do k=1,nz ; do i=is,ie ; pressure(i,k) = 0.0 ; enddo ; enddo
  else
    h_to_rl2_t2 = gv%H_to_RZ * gv%g_Earth
  endif
  !$OMP parallel do default(shared) private(fraz_col,T_fr_set,T_freeze,hc,ps)  &
  !$OMP                             firstprivate(pressure) ! pressure might be set above, so should be firstprivate
  do j=js,je
    ps(:) = 0.0
    if (PRESENT(p_surf)) then ; do i=is,ie
      ps(i) = p_surf(i,j)
    enddo ; endif

    do i=is,ie ; fraz_col(i) = 0.0 ; enddo

    if (cs%pressure_dependent_frazil) then
      do i=is,ie
        pressure(i,1) = ps(i) + (0.5*h_to_rl2_t2)*h(i,j,1)
      enddo
      do k=2,nz ; do i=is,ie
        pressure(i,k) = pressure(i,k-1) + &
          (0.5*h_to_rl2_t2) * (h(i,j,k) + h(i,j,k-1))
      enddo ; enddo
    endif

    if (cs%reclaim_frazil) then
      t_fr_set = .false.
      do i=is,ie ; if (tv%frazil(i,j) > 0.0) then
        if (.not.t_fr_set) then
          call calculate_tfreeze(tv%S(i:,j,1), pressure(i:,1), t_freeze(i:), &
                                 1, ie-i+1, tv%eqn_of_state, pres_scale=us%RL2_T2_to_Pa)
          t_fr_set = .true.
        endif

        if (tv%T(i,j,1) > t_freeze(i)) then
    ! If frazil had previously been formed, but the surface temperature is now
    ! above freezing, cool the surface layer with the frazil heat deficit.
          hc = (tv%C_p*gv%H_to_RZ) * h(i,j,1)
          if (tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i)) <= 0.0) then
            tv%T(i,j,1) = tv%T(i,j,1) - tv%frazil(i,j) / hc
            tv%frazil(i,j) = 0.0
          else
            tv%frazil(i,j) = tv%frazil(i,j) - hc * (tv%T(i,j,1) - t_freeze(i))
            tv%T(i,j,1) = t_freeze(i)
          endif
        endif
      endif ; enddo
    endif

    do k=nz,1,-1
      t_fr_set = .false.
      do i=is,ie
        if ((g%mask2dT(i,j) > 0.0) .and. &
            ((tv%T(i,j,k) < 0.0) .or. (fraz_col(i) > 0.0))) then
          if (.not.t_fr_set) then
            call calculate_tfreeze(tv%S(i:,j,k), pressure(i:,k), t_freeze(i:), &
                                   1, ie-i+1, tv%eqn_of_state, pres_scale=us%RL2_T2_to_Pa)
            t_fr_set = .true.
          endif

          hc = (tv%C_p*gv%H_to_RZ) * h(i,j,k)
          if (h(i,j,k) <= 10.0*(gv%Angstrom_H + gv%H_subroundoff)) then
            ! Very thin layers should not be cooled by the frazil flux.
            if (tv%T(i,j,k) < t_freeze(i)) then
              fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
              tv%T(i,j,k) = t_freeze(i)
            endif
          elseif ((fraz_col(i) > 0.0) .or. (tv%T(i,j,k) < t_freeze(i))) then
            if (fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k)) < 0.0) then
              tv%T(i,j,k) = tv%T(i,j,k) - fraz_col(i) / hc
              fraz_col(i) = 0.0
            else
              fraz_col(i) = fraz_col(i) + hc * (t_freeze(i) - tv%T(i,j,k))
              tv%T(i,j,k) = t_freeze(i)
            endif
          endif
        endif
      enddo
    enddo
    do i=is,ie
      tv%frazil(i,j) = tv%frazil(i,j) + fraz_col(i)
    enddo
  enddo
  call cpu_clock_end(id_clock_frazil)

end subroutine make_frazil

!> This subroutine applies double diffusion to T & S, assuming no diapycal mass
!! fluxes, using a simple triadiagonal solver.
subroutine differential_diffuse_t_s(h, T, S, Kd_T, Kd_S, dt, G, GV)
  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),SZJ_(G),SZK_(G)), &
                           intent(inout) :: t    !< Potential temperature [degC].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(inout) :: s    !< Salinity [PSU] or [gSalt/kg], generically [ppt].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
                           intent(inout)    :: kd_t !< The extra diffusivity of temperature due to
                                                 !! double diffusion relative to the diffusivity of
                                                 !! diffusivity of density [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
                           intent(in)    :: kd_s !< The extra diffusivity of salinity due to
                                                 !! double diffusion relative to the diffusivity of
                                                 !! diffusivity of density [Z2 T-1 ~> m2 s-1].
  real,                    intent(in)    :: dt   !<  Time increment [T ~> s].

  ! local variables
  real, dimension(SZI_(G)) :: &
    b1_t, b1_s, &  !  Variables used by the tridiagonal solvers of T & S [H ~> m or kg m-2].
    d1_t, d1_s     !  Variables used by the tridiagonal solvers [nondim].
  real, dimension(SZI_(G),SZK_(G)) :: &
    c1_t, c1_s     !  Variables used by the tridiagonal solvers [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZK_(G)+1) :: &
    mix_t, mix_s   !  Mixing distances in both directions across each interface [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 :: 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 :: i_h_int      ! The inverse of the thickness associated with an
                       ! interface [H-1 ~> m-1 or m2 kg-1].
  real :: b_denom_t    ! The first term in the denominators for the expressions
  real :: b_denom_s    ! for b1_T and b1_S, both [H ~> m or kg m-2].
  integer :: i, j, k, is, ie, js, je, nz

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

  !$OMP parallel do default(private) shared(is,ie,js,je,h,h_neglect,dt,Kd_T,Kd_S,G,GV,T,S,nz)
  do j=js,je
    do i=is,ie
      i_h_int = 1.0 / (0.5 * (h(i,j,1) + h(i,j,2)) + h_neglect)
      mix_t(i,2) = ((dt * kd_t(i,j,2)) * gv%Z_to_H**2) * i_h_int
      mix_s(i,2) = ((dt * kd_s(i,j,2)) * gv%Z_to_H**2) * i_h_int

      h_tr = h(i,j,1) + h_neglect
      b1_t(i) = 1.0 / (h_tr + mix_t(i,2))
      b1_s(i) = 1.0 / (h_tr + mix_s(i,2))
      d1_t(i) = h_tr * b1_t(i)
      d1_s(i) = h_tr * b1_s(i)
      t(i,j,1) = (b1_t(i)*h_tr)*t(i,j,1)
      s(i,j,1) = (b1_s(i)*h_tr)*s(i,j,1)
    enddo
    do k=2,nz-1 ; do i=is,ie
      ! Calculate the mixing across the interface below this layer.
      i_h_int = 1.0 / (0.5 * (h(i,j,k) + h(i,j,k+1)) + h_neglect)
      mix_t(i,k+1) = ((dt * kd_t(i,j,k+1)) * gv%Z_to_H**2) * i_h_int
      mix_s(i,k+1) = ((dt * kd_s(i,j,k+1)) * gv%Z_to_H**2) * i_h_int

      c1_t(i,k) = mix_t(i,k) * b1_t(i)
      c1_s(i,k) = mix_s(i,k) * b1_s(i)

      h_tr = h(i,j,k) + h_neglect
      b_denom_t = h_tr + d1_t(i)*mix_t(i,k)
      b_denom_s = h_tr + d1_s(i)*mix_s(i,k)
      b1_t(i) = 1.0 / (b_denom_t + mix_t(i,k+1))
      b1_s(i) = 1.0 / (b_denom_s + mix_s(i,k+1))
      d1_t(i) = b_denom_t * b1_t(i)
      d1_s(i) = b_denom_s * b1_s(i)

      t(i,j,k) = b1_t(i) * (h_tr*t(i,j,k) + mix_t(i,k)*t(i,j,k-1))
      s(i,j,k) = b1_s(i) * (h_tr*s(i,j,k) + mix_s(i,k)*s(i,j,k-1))
    enddo ; enddo
    do i=is,ie
      c1_t(i,nz) = mix_t(i,nz) * b1_t(i)
      c1_s(i,nz) = mix_s(i,nz) * b1_s(i)

      h_tr = h(i,j,nz) + h_neglect
      b1_t(i) = 1.0 / (h_tr + d1_t(i)*mix_t(i,nz))
      b1_s(i) = 1.0 / (h_tr + d1_s(i)*mix_s(i,nz))

      t(i,j,nz) = b1_t(i) * (h_tr*t(i,j,nz) + mix_t(i,nz)*t(i,j,nz-1))
      s(i,j,nz) = b1_s(i) * (h_tr*s(i,j,nz) + mix_s(i,nz)*s(i,j,nz-1))
    enddo
    do k=nz-1,1,-1 ; do i=is,ie
      t(i,j,k) = t(i,j,k) + c1_t(i,k+1)*t(i,j,k+1)
      s(i,j,k) = s(i,j,k) + c1_s(i,k+1)*s(i,j,k+1)
    enddo ; enddo
  enddo
end subroutine differential_diffuse_t_s

!> This subroutine keeps salinity from falling below a small but positive threshold.
!! This usually occurs when the ice model attempts to extract more salt then
!! is actually available to it from the ocean.
subroutine adjust_salt(h, tv, G, GV, CS, halo)
  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]
  type(thermo_var_ptrs),   intent(inout) :: tv   !< Structure containing pointers to any
                                                 !! available thermodynamic fields.
  type(diabatic_aux_cs),   intent(in)    :: cs   !< The control structure returned by a previous
                                                 !! call to diabatic_aux_init.
  integer,       optional, intent(in)    :: halo !< Halo width over which to work

  ! local variables
  real :: salt_add_col(szi_(g),szj_(g)) !< The accumulated salt requirement [ppt R Z ~> gSalt m-2]
  real :: s_min      !< The minimum salinity [ppt].
  real :: mc         !< A layer's mass [R Z ~> kg m-2].
  integer :: i, j, k, is, ie, js, je, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  if (present(halo)) then
    is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
  endif

!  call cpu_clock_begin(id_clock_adjust_salt)

  s_min = tv%min_salinity

  salt_add_col(:,:) = 0.0

  !$OMP parallel do default(shared) private(mc)
  do j=js,je
    do k=nz,1,-1 ; do i=is,ie
      if ( (g%mask2dT(i,j) > 0.0) .and. &
           ((tv%S(i,j,k) < s_min) .or. (salt_add_col(i,j) > 0.0)) ) then
        mc = gv%H_to_RZ * h(i,j,k)
        if (h(i,j,k) <= 10.0*gv%Angstrom_H) then
          ! Very thin layers should not be adjusted by the salt flux
          if (tv%S(i,j,k) < s_min) then
            salt_add_col(i,j) = salt_add_col(i,j) +  mc * (s_min - tv%S(i,j,k))
            tv%S(i,j,k) = s_min
          endif
        elseif (salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k)) <= 0.0) then
          tv%S(i,j,k) = tv%S(i,j,k) - salt_add_col(i,j) / mc
          salt_add_col(i,j) = 0.0
        else
          salt_add_col(i,j) = salt_add_col(i,j) + mc * (s_min - tv%S(i,j,k))
          tv%S(i,j,k) = s_min
        endif
      endif
    enddo ; enddo
    do i=is,ie
      tv%salt_deficit(i,j) = tv%salt_deficit(i,j) + salt_add_col(i,j)
    enddo
  enddo
!  call cpu_clock_end(id_clock_adjust_salt)

end subroutine adjust_salt

!> This is a simple tri-diagonal solver for T and S.
!! "Simple" means it only uses arrays hold, ea and eb.
subroutine tridiagts(G, GV, is, ie, js, je, hold, ea, eb, T, S)
  type(ocean_grid_type),                    intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),                  intent(in)    :: gv   !< The ocean's vertical grid structure
  integer,                                  intent(in)    :: is   !< The start i-index to work on.
  integer,                                  intent(in)    :: ie   !< The end i-index to work on.
  integer,                                  intent(in)    :: js   !< The start j-index to work on.
  integer,                                  intent(in)    :: je   !< The end j-index to work on.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: hold !< The layer thicknesses before entrainment,
                                                                  !! [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: 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(in)    :: eb !< The amount of fluid entrained from the layer
                                                                !! below within this time step [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: t  !< Layer potential temperatures [degC].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: s  !< Layer salinities [ppt].

  ! Local variables
  real :: b1(szib_(g)), d1(szib_(g)) ! b1, c1, and d1 are variables used by the
  real :: c1(szib_(g),szk_(g))       ! tridiagonal solver.
  real :: h_tr, b_denom_1
  integer :: i, j, k

  !$OMP parallel do default(shared) private(h_tr,b1,d1,c1,b_denom_1)
  do j=js,je
    do i=is,ie
      h_tr = hold(i,j,1) + gv%H_subroundoff
      b1(i) = 1.0 / (h_tr + eb(i,j,1))
      d1(i) = h_tr * b1(i)
      t(i,j,1) = (b1(i)*h_tr)*t(i,j,1)
      s(i,j,1) = (b1(i)*h_tr)*s(i,j,1)
    enddo
    do k=2,g%ke ; do i=is,ie
      c1(i,k) = eb(i,j,k-1) * b1(i)
      h_tr = hold(i,j,k) + gv%H_subroundoff
      b_denom_1 = h_tr + d1(i)*ea(i,j,k)
      b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
      d1(i) = b_denom_1 * b1(i)
      t(i,j,k) = b1(i) * (h_tr*t(i,j,k) + ea(i,j,k)*t(i,j,k-1))
      s(i,j,k) = b1(i) * (h_tr*s(i,j,k) + ea(i,j,k)*s(i,j,k-1))
    enddo ; enddo
    do k=g%ke-1,1,-1 ; do i=is,ie
      t(i,j,k) = t(i,j,k) + c1(i,k+1)*t(i,j,k+1)
      s(i,j,k) = s(i,j,k) + c1(i,k+1)*s(i,j,k+1)
    enddo ; enddo
  enddo
end subroutine tridiagts

!>   This subroutine calculates u_h and v_h (velocities at thickness
!! points), optionally using the entrainment amounts passed in as arguments.
subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, US, ea, eb)
  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(SZIB_(G),SZJ_(G),SZK_(G)), &
                             intent(in)  :: u    !< The zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                             intent(in)  :: v    !< The meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                             intent(in)  :: h    !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                             intent(out)   :: u_h !< Zonal velocity interpolated to h points [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                             intent(out)   :: v_h !< Meridional velocity interpolated to h points [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                     optional, intent(in)  :: ea !< The amount of fluid entrained from the layer
                                                 !! above within this time step [H ~> m or kg m-2].
                                                 !! Omitting ea is the same as setting it to 0.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                     optional, intent(in)  :: eb !< The amount of fluid entrained from the layer
                                                 !! below within this time step [H ~> m or kg m-2].
                                                 !! Omitting eb is the same as setting it to 0.

  ! local variables
  real :: b_denom_1    ! The first term in the denominator of b1 [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 :: b1(szi_(g)), d1(szi_(g)), c1(szi_(g),szk_(g))
  real :: a_n(szi_(g)), a_s(szi_(g))  ! Fractional weights of the neighboring
  real :: a_e(szi_(g)), a_w(szi_(g))  ! velocity points, ~1/2 in the open
                                      ! ocean, nondimensional.
  real :: sum_area, idenom
  logical :: mix_vertically
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
  call cpu_clock_begin(id_clock_uv_at_h)
  h_neglect = gv%H_subroundoff

  mix_vertically = present(ea)
  if (present(ea) .neqv. present(eb)) call mom_error(fatal, &
      "find_uv_at_h: Either both ea and eb or neither one must be present "// &
      "in call to find_uv_at_h.")
!$OMP parallel do default(none) shared(is,ie,js,je,G,GV,mix_vertically,h,h_neglect, &
!$OMP                                  eb,u_h,u,v_h,v,nz,ea)                     &
!$OMP                          private(sum_area,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1)
  do j=js,je
    do i=is,ie
      sum_area = g%areaCu(i-1,j) + g%areaCu(i,j)
      if (sum_area>0.0) then
        idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
        a_w(i) = g%areaCu(i-1,j) * idenom
        a_e(i) = g%areaCu(i,j) * idenom
      else
        a_w(i) = 0.0 ; a_e(i) = 0.0
      endif

      sum_area = g%areaCv(i,j-1) + g%areaCv(i,j)
      if (sum_area>0.0) then
        idenom = sqrt(0.5*g%IareaT(i,j) / sum_area)
        a_s(i) = g%areaCv(i,j-1) * idenom
        a_n(i) = g%areaCv(i,j) * idenom
      else
        a_s(i) = 0.0 ; a_n(i) = 0.0
      endif
    enddo

    if (mix_vertically) then
      do i=is,ie
        b_denom_1 = h(i,j,1) + h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
        d1(i) = b_denom_1 * b1(i)
        u_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_e(i)*u(i,j,1) + a_w(i)*u(i-1,j,1))
        v_h(i,j,1) = (h(i,j,1)*b1(i)) * (a_n(i)*v(i,j,1) + a_s(i)*v(i,j-1,1))
      enddo
      do k=2,nz ; do i=is,ie
        c1(i,k) = eb(i,j,k-1) * b1(i)
        b_denom_1 = h(i,j,k) + d1(i)*ea(i,j,k) + h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
        d1(i) = b_denom_1 * b1(i)
        u_h(i,j,k) = (h(i,j,k) * (a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)) + &
                      ea(i,j,k)*u_h(i,j,k-1))*b1(i)
        v_h(i,j,k) = (h(i,j,k) * (a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)) + &
                      ea(i,j,k)*v_h(i,j,k-1))*b1(i)
      enddo ; enddo
      do k=nz-1,1,-1 ; do i=is,ie
        u_h(i,j,k) = u_h(i,j,k) + c1(i,k+1)*u_h(i,j,k+1)
        v_h(i,j,k) = v_h(i,j,k) + c1(i,k+1)*v_h(i,j,k+1)
      enddo ; enddo
    else
      do k=1,nz ; do i=is,ie
        u_h(i,j,k) = a_e(i)*u(i,j,k) + a_w(i)*u(i-1,j,k)
        v_h(i,j,k) = a_n(i)*v(i,j,k) + a_s(i)*v(i,j-1,k)
      enddo ; enddo
    endif
  enddo

  call cpu_clock_end(id_clock_uv_at_h)
end subroutine find_uv_at_h


subroutine set_pen_shortwave(optics, fluxes, G, GV, US, CS, opacity_CSp, tracer_flow_CSp)
  type(optics_type),       pointer       :: optics !< An optics structure that has will contain
                                                   !! information about shortwave fluxes and absorption.
  type(forcing),           intent(inout) :: fluxes !< points to forcing fields
                                                   !! unused fields have NULL ptrs
  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(diabatic_aux_cs),   pointer       :: cs     !< Control structure for diabatic_aux
  type(opacity_cs),        pointer       :: opacity_csp !< The control structure for the opacity module.
  type(tracer_flow_control_cs), pointer  :: tracer_flow_csp !< A pointer to the control structure
                                                   !! organizing the tracer modules.

  ! Local variables
  real, dimension(SZI_(G),SZJ_(G))          :: chl_2d !< Vertically uniform chlorophyll-A concentractions [mg m-3]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: chl_3d !< The chlorophyll-A concentractions of each layer [mg m-3]
  character(len=128) :: mesg
  integer :: i, j, k, is, ie, js, je
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  if (.not.associated(optics)) return

  if (cs%var_pen_sw) then
    if (cs%chl_from_file) then
      ! Only the 2-d surface chlorophyll can be read in from a file.  The
      ! same value is assumed for all layers.
      call time_interp_external(cs%sbc_chl, cs%Time, chl_2d)
      do j=js,je ; do i=is,ie
        if ((g%mask2dT(i,j) > 0.5) .and. (chl_2d(i,j) < 0.0)) then
          write(mesg,'(" Time_interp negative chl of ",(1pe12.4)," at i,j = ",&
                    & 2(i3), "lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
                     chl_2d(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j)
          call mom_error(fatal, "MOM_diabatic_aux set_pen_shortwave: "//trim(mesg))
        endif
      enddo ; enddo

      if (cs%id_chl > 0) call post_data(cs%id_chl, chl_2d, cs%diag)

      call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
                       fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, us, opacity_csp, chl_2d=chl_2d)
    else
      if (.not.associated(tracer_flow_csp)) call mom_error(fatal, &
        "The tracer flow control structure must be associated when the model sets "//&
        "the chlorophyll internally in set_pen_shortwave.")
      call get_chl_from_model(chl_3d, g, tracer_flow_csp)

      if (cs%id_chl > 0) call post_data(cs%id_chl, chl_3d(:,:,1), cs%diag)

      call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
                       fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, us, opacity_csp, chl_3d=chl_3d)
    endif
  else
    call set_opacity(optics, fluxes%sw, fluxes%sw_vis_dir, fluxes%sw_vis_dif, &
                     fluxes%sw_nir_dir, fluxes%sw_nir_dif, g, gv, us, opacity_csp)
  endif

end subroutine set_pen_shortwave


!> Diagnose a mixed layer depth (MLD) determined by a given density difference with the surface.
!> This routine is appropriate in MOM_diabatic_driver due to its position within the time stepping.
subroutine diagnosemldbydensitydifference(id_MLD, h, tv, densityDiff, G, GV, US, diagPtr, &
                                          id_N2subML, id_MLDsq, dz_subML)
  type(ocean_grid_type),   intent(in) :: g           !< Grid type
  type(verticalgrid_type), intent(in) :: gv          !< ocean vertical grid structure
  type(unit_scale_type),   intent(in) :: us          !< A dimensional unit scaling type
  integer,                 intent(in) :: id_mld      !< Handle (ID) of MLD diagnostic
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: h           !< Layer thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),   intent(in) :: tv          !< Structure containing pointers to any
                                                     !! available thermodynamic fields.
  real,                    intent(in) :: densitydiff !< Density difference to determine MLD [R ~> kg m-3]
  type(diag_ctrl),         pointer    :: diagptr     !< Diagnostics structure
  integer,       optional, intent(in) :: id_n2subml  !< Optional handle (ID) of subML stratification
  integer,       optional, intent(in) :: id_mldsq    !< Optional handle (ID) of squared MLD
  real,          optional, intent(in) :: dz_subml    !< The distance over which to calculate N2subML
                                                     !! or 50 m if missing [Z ~> m]

  ! Local variables
  real, dimension(SZI_(G)) :: deltarhoatkm1, deltarhoatk ! Density differences [R ~> kg m-3].
  real, dimension(SZI_(G)) :: pref_mld, pref_n2 ! Reference pressures [R L2 T-2 ~> Pa].
  real, dimension(SZI_(G)) :: h_subml, dh_n2    ! Summed thicknesses used in N2 calculation [H ~> m].
  real, dimension(SZI_(G)) :: t_subml, t_deeper ! Temperatures used in the N2 calculation [degC].
  real, dimension(SZI_(G)) :: s_subml, s_deeper ! Salinities used in the N2 calculation [ppt].
  real, dimension(SZI_(G)) :: rho_subml, rho_deeper ! Densities used in the N2 calculation [R ~> kg m-3].
  real, dimension(SZI_(G)) :: dk, dkm1          ! Depths [Z ~> m].
  real, dimension(SZI_(G)) :: rhosurf          ! Density used in finding the mixedlayer depth [R ~> kg m-3].
  real, dimension(SZI_(G), SZJ_(G)) :: mld     ! Diagnosed mixed layer depth [Z ~> m].
  real, dimension(SZI_(G), SZJ_(G)) :: submln2 ! Diagnosed stratification below ML [T-2 ~> s-2].
  real, dimension(SZI_(G), SZJ_(G)) :: mld2    ! Diagnosed MLD^2 [Z2 ~> m2].
  logical, dimension(SZI_(G)) :: n2_region_set ! If true, all necessary values for calculating N2
                                               ! have been stored already.
  real :: ge_rho0          ! The gravitational acceleration divided by a mean density [Z T-2 R-1 ~> m4 s-2 kg-1].
  real :: dh_subml         ! Depth below ML over which to diagnose stratification [H ~> m].
  real :: afac             ! A nondimensional factor [nondim]
  real :: ddrho            ! A density difference [R ~> kg m-3]
  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
  integer :: i, j, is, ie, js, je, k, nz, id_n2, id_sq

  id_n2 = -1 ; if (PRESENT(id_n2subml)) id_n2 = id_n2subml

  id_sq = -1 ; if (PRESENT(id_mldsq)) id_sq = id_mldsq

  ge_rho0 = us%L_to_Z**2*gv%g_Earth / (gv%Rho0)
  dh_subml = 50.*gv%m_to_H  ; if (present(dz_subml)) dh_subml = gv%Z_to_H*dz_subml

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

  pref_mld(:) = 0.0
  eosdom(:) = eos_domain(g%HI)
  do j=js,je
    do i=is,ie ; dk(i) = 0.5 * h(i,j,1) * gv%H_to_Z ; enddo ! Depth of center of surface layer
    call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pref_mld, rhosurf, tv%eqn_of_state, eosdom)
    do i=is,ie
      deltarhoatk(i) = 0.
      mld(i,j) = 0.
      if (id_n2>0) then
        submln2(i,j) = 0.0
        h_subml(i) = h(i,j,1) ; dh_n2(i) = 0.0
        t_subml(i) = 0.0  ; s_subml(i) = 0.0 ; t_deeper(i) = 0.0 ; s_deeper(i) = 0.0
        n2_region_set(i) = (g%mask2dT(i,j)<0.5) ! Only need to work on ocean points.
      endif
    enddo
    do k=2,nz
      do i=is,ie
        dkm1(i) = dk(i) ! Depth of center of layer K-1
        dk(i) = dk(i) + 0.5 * ( h(i,j,k) + h(i,j,k-1) ) * gv%H_to_Z ! Depth of center of layer K
      enddo

      ! Prepare to calculate stratification, N2, immediately below the mixed layer by finding
      ! the cells that extend over at least dz_subML.
      if (id_n2>0) then
         do i=is,ie
          if (mld(i,j)==0.0) then  ! Still in the mixed layer.
            h_subml(i) = h_subml(i) + h(i,j,k)
          elseif (.not.n2_region_set(i)) then ! This block is below the mixed layer, but N2 has not been found yet.
            if (dh_n2(i)==0.0) then ! Record the temperature, salinity, pressure, immediately below the ML
              t_subml(i) = tv%T(i,j,k) ; s_subml(i) = tv%S(i,j,k)
              h_subml(i) = h_subml(i) + 0.5 * h(i,j,k) ! Start midway through this layer.
              dh_n2(i) = 0.5 * h(i,j,k)
            elseif (dh_n2(i) + h(i,j,k) < dh_subml) then
              dh_n2(i) = dh_n2(i) + h(i,j,k)
            else  ! This layer includes the base of the region where N2 is calculated.
              t_deeper(i) = tv%T(i,j,k) ; s_deeper(i) = tv%S(i,j,k)
              dh_n2(i) = dh_n2(i) + 0.5 * h(i,j,k)
              n2_region_set(i) = .true.
            endif
          endif
        enddo ! i-loop
      endif ! id_N2>0

      ! Mixed-layer depth, using sigma-0 (surface reference pressure)
      do i=is,ie ; deltarhoatkm1(i) = deltarhoatk(i) ; enddo ! Store value from previous iteration of K
      call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_mld, deltarhoatk, tv%eqn_of_state, eosdom)
      do i = is, ie
        deltarhoatk(i) = deltarhoatk(i) - rhosurf(i) ! Density difference between layer K and surface
        ddrho = deltarhoatk(i) - deltarhoatkm1(i)
        if ((mld(i,j)==0.) .and. (ddrho>0.) .and. &
            (deltarhoatkm1(i)<densitydiff) .and. (deltarhoatk(i)>=densitydiff)) then
          afac = ( densitydiff - deltarhoatkm1(i) ) / ddrho
          mld(i,j) = dk(i) * afac + dkm1(i) * (1. - afac)
        endif
        if (id_sq > 0) mld2(i,j) = mld(i,j)**2
      enddo ! i-loop
    enddo ! k-loop
    do i=is,ie
      if ((mld(i,j)==0.) .and. (deltarhoatk(i)<densitydiff)) mld(i,j) = dk(i) ! Assume mixing to the bottom
    enddo

    if (id_n2>0) then  ! Now actually calculate stratification, N2, below the mixed layer.
      do i=is,ie ; pref_n2(i) = (gv%g_Earth * gv%H_to_RZ) * (h_subml(i) + 0.5*dh_n2(i)) ; enddo
      ! if ((.not.N2_region_set(i)) .and. (dH_N2(i) > 0.5*dH_subML)) then
      !    ! Use whatever stratification we can, measured over whatever distance is available?
      !    T_deeper(i) = tv%T(i,j,nz) ; S_deeper(i) = tv%S(i,j,nz)
      !    N2_region_set(i) = .true.
      ! endif
      call calculate_density(t_subml, s_subml, pref_n2, rho_subml, tv%eqn_of_state, eosdom)
      call calculate_density(t_deeper, s_deeper, pref_n2, rho_deeper, tv%eqn_of_state, eosdom)
      do i=is,ie ; if ((g%mask2dT(i,j)>0.5) .and. n2_region_set(i)) then
        submln2(i,j) =  ge_rho0 * (rho_deeper(i) - rho_subml(i)) / (gv%H_to_z * dh_n2(i))
      endif ; enddo
    endif
  enddo ! j-loop

  if (id_mld > 0) call post_data(id_mld, mld, diagptr)
  if (id_n2 > 0)  call post_data(id_n2, submln2, diagptr)
  if (id_sq > 0)  call post_data(id_sq, mld2, diagptr)

end subroutine diagnosemldbydensitydifference

!> Diagnose a mixed layer depth (MLD) determined by the depth a given energy value would mix.
!> This routine is appropriate in MOM_diabatic_driver due to its position within the time stepping.
subroutine diagnosemldbyenergy(id_MLD, h, tv, G, GV, US, Mixing_Energy, diagPtr)
  ! Author: Brandon Reichl
  ! Date: October 2, 2020
  ! //
  ! *Note that gravity is assumed constant everywhere and divided out of all calculations.
  !
  ! This code has been written to step through the columns layer by layer, summing the PE
  ! change inferred by mixing the layer with all layers above.  When the change exceeds a
  ! threshold (determined by input array Mixing_Energy), the code needs to solve for how far
  ! into this layer the threshold PE change occurs (assuming constant density layers).
  ! This is expressed here via solving the function F(X) = 0 where:
  ! F(X) = 0.5 * ( Ca*X^3/(D1+X) + Cb*X^2/(D1+X) + Cc*X/(D1+X) + Dc/(D1+X)
  !                + Ca2*X^2 + Cb2*X + Cc2)
  ! where all coefficients are determined by the previous mixed layer depth, the
  ! density of the previous mixed layer, the present layer thickness, and the present
  ! layer density.  This equation is worked out by computing the total PE assuming constant
  ! density in the mixed layer as well as in the remaining part of the present layer that is
  ! not mixed.
  ! To solve for X in this equation a Newton's method iteration is employed, which
  ! converges extremely quickly (usually 1 guess) since this equation turns out being rather
  ! lienar for PE change with increasing X.
  ! Input parameters:
  integer, dimension(3),   intent(in) :: id_mld      !< Energy output diag IDs
  type(ocean_grid_type),   intent(in) :: g           !< Grid type
  type(verticalgrid_type), intent(in) :: gv          !< ocean vertical grid structure
  type(unit_scale_type),   intent(in) :: us          !< A dimensional unit scaling type
  real, dimension(3),      intent(in) :: mixing_energy !< Energy values for up to 3 MLDs
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: h           !< Layer thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),   intent(in) :: tv          !< Structure containing pointers to any
                                                     !! available thermodynamic fields.
  type(diag_ctrl),         pointer    :: diagptr     !< Diagnostics structure

  ! Local variables
  real, dimension(SZI_(G), SZJ_(G),3) :: mld     ! Diagnosed mixed layer depth [Z ~> m].
  real, dimension(SZK_(G)) :: z_l, z_u, dz, rho_c, pref_mld
  real, dimension(3) :: pe_threshold

  real :: ig, e_g
  real :: pe_threshold_fraction, pe, pe_mixed, pe_mixed_tst
  real :: rhodz_ml, h_ml, rhodz_ml_tst, h_ml_tst
  real :: rho_ml

  real :: r1, d1, r2, d2
  real :: ca, cb,d ,cc, cd, ca2, cb2, c, cc2
  real :: gx, gpx, hx, ihx, hpx, ix, ipx, fgx, fpx, x, x2

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

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

  pref_mld(:) = 0.0
  mld(:,:,:) = 0.0
  pe_threshold_fraction = 1.e-4 !Fixed threshold of 0.01%, could be runtime.

  do im=1,3
    pe_threshold(im) = mixing_energy(im)/gv%g_earth
  enddo

  do j=js,je; do i=is,ie
    if (g%mask2dT(i,j) > 0.0) then

      call calculate_density(tv%T(i,j,:), tv%S(i,j,:), pref_mld, rho_c, 1, nz, &
                             tv%eqn_of_state, scale=us%kg_m3_to_R)

      do k=1,nz
        dz(k) = h(i,j,k) * gv%H_to_Z
      enddo
      z_u(1) = 0.0
      z_l(1) = -dz(1)
      do k=2,nz
        z_u(k) = z_l(k-1)
        z_l(k) = z_l(k-1)-dz(k)
      enddo

      do im=1,3

        ! Initialize these for each column-wise calculation
        pe = 0.0
        rhodz_ml = 0.0
        h_ml = 0.0
        rhodz_ml_tst = 0.0
        h_ml_tst = 0.0
        pe_mixed = 0.0

        do k=1,nz

          ! This is the unmixed PE cumulative sum from top down
          pe = pe + 0.5 * rho_c(k) * (z_u(k)**2 - z_l(k)**2)

          ! This is the depth and integral of density
          h_ml_tst = h_ml + dz(k)
          rhodz_ml_tst = rhodz_ml + rho_c(k) * dz(k)

          ! The average density assuming all layers including this were mixed
          rho_ml = rhodz_ml_tst/h_ml_tst

          ! The PE assuming all layers including this were mixed
          ! Note that 0. could be replaced with "Surface", which doesn't have to be 0
          ! but 0 is a good reference value.
          pe_mixed_tst = 0.5 * rho_ml * (0.**2 - (0. - h_ml_tst)**2)

          ! Check if we supplied enough energy to mix to this layer
          if (pe_mixed_tst - pe <= pe_threshold(im)) then
            h_ml = h_ml_tst
            rhodz_ml = rhodz_ml_tst

          else ! If not, we need to solve where the energy ran out
            ! This will be done with a Newton's method iteration:

            r1 = rhodz_ml / h_ml ! The density of the mixed layer (not including this layer)
            d1 = h_ml ! The thickness of the mixed layer (not including this layer)
            r2 = rho_c(k) ! The density of this layer
            d2 = dz(k) ! The thickness of this layer

            ! This block could be used to calculate the function coefficients if
            ! we don't reference all values to a surface designated as z=0
            ! S = Surface
            ! Ca  = -(R2)
            ! Cb  = -( (R1*D1) + R2*(2.*D1-2.*S) )
            ! D   = D1**2. - 2.*D1*S
            ! Cc  = -( R1*D1*(2.*D1-2.*S) + R2*D )
            ! Cd  = -(R1*D1*D)
            ! Ca2 = R2
            ! Cb2 = R2*(2*D1-2*S)
            ! C   = S**2 + D2**2 + D1**2 - 2*D1*S - 2.*D2*S +2.*D1*D2
            ! Cc2 = R2*(D+S**2-C)
            !
            ! If the surface is S = 0, it simplifies to:
            ca  = -r2
            cb  = -(r1 * d1 + r2 * (2. * d1))
            d   = d1**2
            cc  = -(r1 * d1 * (2. * d1) + (r2 * d))
            cd  = -r1 * (d1 * d)
            ca2 = r2
            cb2 = r2 * (2. * d1)
            c   = d2**2 + d1**2 + 2. * (d1 * d2)
            cc2 = r2 * (d - c)

            ! First guess for an iteration using Newton's method
            x = dz(k) * 0.5

            it=0
            do while(it<10)!We can iterate up to 10 times
              ! We are trying to solve the function:
              ! F(x) = G(x)/H(x)+I(x)
              ! for where F(x) = PE+PE_threshold, or equivalently for where
              ! F(x) = G(x)/H(x)+I(x) - (PE+PE_threshold) = 0
              ! We also need the derivative of this function for the Newton's method iteration
              ! F'(x) = (G'(x)H(x)-G(x)H'(x))/H(x)^2 + I'(x)
              ! G and its derivative
              gx = 0.5 * (ca * (x*x*x) + cb * x**2 + cc * x + cd)
              gpx = 0.5 * (3. * (ca * x**2) + 2. * (cb * x) + cc)
              ! H, its inverse, and its derivative
              hx = d1 + x
              ihx = 1. / hx
              hpx = 1.
              ! I and its derivative
              ix = 0.5 * (ca2 * x**2 + cb2 * x + cc2)
              ipx = 0.5 * (2. * ca2 * x + cb2)

              ! The Function and its derivative:
              pe_mixed = gx * ihx + ix
              fgx = pe_mixed - (pe + pe_threshold(im))
              fpx = (gpx * hx - hpx * gx) * ihx**2 + ipx

              ! Check if our solution is within the threshold bounds, if not update
              ! using Newton's method.  This appears to converge almost always in
              ! one step because the function is very close to linear in most applications.
              if (abs(fgx) > pe_threshold(im) * pe_threshold_fraction) then
                x2 = x - fgx / fpx
                it = it + 1
                if (x2 < 0. .or. x2 > dz(k)) then
                  ! The iteration seems to be robust, but we need to do something *if*
                  ! things go wrong... How should we treat failed iteration?
                  ! Present solution: Stop trying to compute and just say we can't mix this layer.
                  x=0
                  exit
                else
                  x = x2
                endif
              else
                exit! Quit the iteration
              endif
            enddo
            h_ml = h_ml + x
            exit! Quit looping through the column
          endif
        enddo
        mld(i,j,im) = h_ml
      enddo
    else
      mld(i,j,:) = 0.0
    endif
  enddo ; enddo

  if (id_mld(1) > 0) call post_data(id_mld(1), mld(:,:,1), diagptr)
  if (id_mld(2) > 0) call post_data(id_mld(2), mld(:,:,2), diagptr)
  if (id_mld(3) > 0) call post_data(id_mld(3), mld(:,:,3), diagptr)

end subroutine diagnosemldbyenergy

!> Update the thickness, temperature, and salinity due to thermodynamic
!! boundary forcing (contained in fluxes type) applied to h, tv%T and tv%S,
!! and calculate the TKE implications of this heating.
subroutine applyboundaryfluxesinout(CS, G, GV, US, dt, fluxes, optics, nsw, h, tv, &
                                    aggregate_FW_forcing, evap_CFL_limit, &
                                    minimum_forcing_depth, cTKE, dSV_dT, dSV_dS, &
                                    SkinBuoyFlux )
  type(diabatic_aux_cs),   pointer       :: cs !< Control structure for diabatic_aux
  type(ocean_grid_type),   intent(in)    :: g  !< Grid structure
  type(verticalgrid_type), intent(in)    :: gv !< ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: us !< A dimensional unit scaling type
  real,                    intent(in)    :: dt !< Time-step over which forcing is applied [T ~> s]
  type(forcing),           intent(inout) :: fluxes !< Surface fluxes container
  type(optics_type),       pointer       :: optics !< Optical properties container
  integer,                 intent(in)    :: nsw !< The number of frequency bands of penetrating
                                                !! shortwave radiation
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(inout) :: h  !< Layer thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),   intent(inout) :: tv !< Structure containing pointers to any
                                               !! available thermodynamic fields.
  logical,                 intent(in)    :: aggregate_fw_forcing !< If False, treat in/out fluxes separately.
  real,                    intent(in)    :: evap_cfl_limit !< The largest fraction of a layer that
                                               !! can be evaporated in one time-step [nondim].
  real,                    intent(in)    :: minimum_forcing_depth !< The smallest depth over which
                                               !! heat and freshwater fluxes is applied [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: ctke !< Turbulent kinetic energy requirement to mix
                                               !! forcing through each layer [R Z3 T-2 ~> J m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: dsv_dt !< Partial derivative of specific volume with
                                               !! potential temperature [R-1 degC-1 ~> m3 kg-1 degC-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: dsv_ds !< Partial derivative of specific volume with
                                               !! salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
  real, dimension(SZI_(G),SZJ_(G)), &
                   optional, intent(out) :: skinbuoyflux !< Buoyancy flux at surface [Z2 T-3 ~> m2 s-3].

  ! Local variables
  integer, parameter :: maxgroundings = 5
  integer :: numberofgroundings, iground(maxgroundings), jground(maxgroundings)
  real :: h_limit_fluxes
  real :: iforcingdepthscale
  real :: idt        ! The inverse of the timestep [T-1 ~> s-1]
  real :: dthickness, dtemp, dsalt
  real :: fractionofforcing, hold, ithickness
  real :: rivermixconst  ! A constant used in implementing river mixing [R Z2 T-1 ~> Pa s].

  real, dimension(SZI_(G)) :: &
    d_pres,       &  ! pressure change across a layer [R L2 T-2 ~> Pa]
    p_lay,        &  ! average pressure in a layer [R L2 T-2 ~> Pa]
    pres,         &  ! pressure at an interface [R L2 T-2 ~> Pa]
    netmassinout, &  ! surface water fluxes [H ~> m or kg m-2] over time step
    netmassin,    &  ! mass entering ocean surface [H ~> m or kg m-2] over a time step
    netmassout,   &  ! mass leaving ocean surface [H ~> m or kg m-2] over a time step
    netheat,      &  ! heat via surface fluxes excluding Pen_SW_bnd and netMassOut
                     ! [degC H ~> degC m or degC kg m-2]
    netsalt,      &  ! surface salt flux ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
                     ! [ppt H ~> ppt m or ppt kg m-2]
    nonpensw,     &  ! non-downwelling SW, which is absorbed at ocean surface
                     ! [degC H ~> degC m or degC kg m-2]
    surfpressure, &  ! Surface pressure (approximated as 0.0) [R L2 T-2 ~> Pa]
    drhodt,       &  ! change in density per change in temperature [R degC-1 ~> kg m-3 degC-1]
    drhods,       &  ! change in density per change in salinity [R ppt-1 ~> kg m-3 ppt-1]
    netheat_rate, &  ! netheat but for dt=1 [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
    netsalt_rate, &  ! netsalt but for dt=1 (e.g. returns a rate)
                     ! [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
    netmassinout_rate! netmassinout but for dt=1 [H T-1 ~> m s-1 or kg m-2 s-1]
  real, dimension(SZI_(G), SZK_(G)) :: &
    h2d, &           ! A 2-d copy of the thicknesses [H ~> m or kg m-2]
    t2d, &           ! A 2-d copy of the layer temperatures [degC]
    pen_tke_2d, &    ! The TKE required to homogenize the heating by shortwave radiation within
                     ! a layer [R Z3 T-2 ~> J m-2]
    dsv_dt_2d        ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
  real, dimension(SZI_(G)) :: &
    netpen_rate      ! The surface penetrative shortwave heating rate summed over all bands
                     ! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
  real, dimension(max(nsw,1),SZI_(G)) :: &
    pen_sw_bnd, &    ! The penetrative shortwave heating integrated over a timestep by band
                     ! [degC H ~> degC m or degC kg m-2]
    pen_sw_bnd_rate  ! The penetrative shortwave heating rate by band
                     ! [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
  real, dimension(max(nsw,1),SZI_(G),SZK_(G)) :: &
    opacityband      ! The opacity (inverse of the exponential absorption length) of each frequency
                     ! band of shortwave radation in each layer [H-1 ~> m-1 or m2 kg-1]
  real, dimension(maxGroundings) :: hgrounding ! Thickness added by each grounding event [H ~> m or kg m-2]
  real    :: temp_in, salin_in
  real    :: g_hconv2 ! A conversion factor for use in the TKE calculation
                      ! in units of [Z3 R2 T-2 H-2 ~> kg2 m-5 s-2 or m s-2].
  real    :: gorho    ! g_Earth times a unit conversion factor divided by density
                      ! [Z T-2 R-1 ~> m4 s-2 kg-1]
  logical :: calculate_energetics
  logical :: calculate_buoyancy
  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
  integer :: i, j, is, ie, js, je, k, nz, n, nb
  character(len=45) :: mesg

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

  idt = 1.0 / dt

  calculate_energetics = (present(ctke) .and. present(dsv_dt) .and. present(dsv_ds))
  calculate_buoyancy = present(skinbuoyflux)
  if (calculate_buoyancy) skinbuoyflux(:,:) = 0.0
  if (present(ctke)) ctke(:,:,:) = 0.0
  g_hconv2 = (us%L_to_Z**2*gv%g_Earth * gv%H_to_RZ) * gv%H_to_RZ
  eosdom(:) = eos_domain(g%HI)

  ! Only apply forcing if fluxes%sw is associated.
  if (.not.associated(fluxes%sw) .and. .not.calculate_energetics) return

  if (calculate_buoyancy) then
    surfpressure(:) = 0.0
    gorho = us%L_to_Z**2*gv%g_Earth / gv%Rho0
  endif

  ! H_limit_fluxes is used by extractFluxes1d to scale down fluxes if the total
  ! depth of the ocean is vanishing. It does not (yet) handle a value of zero.
  ! To accommodate vanishing upper layers, we need to allow for an instantaneous
  ! distribution of forcing over some finite vertical extent. The bulk mixed layer
  ! code handles this issue properly.
  h_limit_fluxes = max(gv%Angstrom_H, 1.e-30*gv%m_to_H)

  ! diagnostic to see if need to create mass to avoid grounding
  if (cs%id_createdH>0) cs%createdH(:,:) = 0.
  numberofgroundings = 0

  !$OMP parallel do default(none) shared(is,ie,js,je,nz,h,tv,nsw,G,GV,US,optics,fluxes,    &
  !$OMP                                  H_limit_fluxes,numberOfGroundings,iGround,jGround,&
  !$OMP                                  nonPenSW,hGrounding,CS,Idt,aggregate_FW_forcing,  &
  !$OMP                                  minimum_forcing_depth,evap_CFL_limit,dt,EOSdom,   &
  !$OMP                                  calculate_buoyancy,netPen_rate,SkinBuoyFlux,GoRho, &
  !$OMP                                  calculate_energetics,dSV_dT,dSV_dS,cTKE,g_Hconv2) &
  !$OMP                          private(opacityBand,h2d,T2d,netMassInOut,netMassOut,      &
  !$OMP                                  netHeat,netSalt,Pen_SW_bnd,fractionOfForcing,     &
  !$OMP                                  IforcingDepthScale,                               &
  !$OMP                                  dThickness,dTemp,dSalt,hOld,Ithickness,           &
  !$OMP                                  netMassIn,pres,d_pres,p_lay,dSV_dT_2d,            &
  !$OMP                                  netmassinout_rate,netheat_rate,netsalt_rate,      &
  !$OMP                                  drhodt,drhods,pen_sw_bnd_rate,                    &
  !$OMP                                  pen_TKE_2d,Temp_in,Salin_in,RivermixConst)        &
  !$OMP                     firstprivate(SurfPressure)
  do j=js,je
  ! Work in vertical slices for efficiency

    ! Copy state into 2D-slice arrays
    do k=1,nz ; do i=is,ie
      h2d(i,k) = h(i,j,k)
      t2d(i,k) = tv%T(i,j,k)
    enddo ; enddo

    if (calculate_energetics) then
      ! The partial derivatives of specific volume with temperature and
      ! salinity need to be precalculated to avoid having heating of
      ! tiny layers give nonsensical values.
      if (associated(tv%p_surf)) then
        do i=is,ie ; pres(i) = tv%p_surf(i,j) ; enddo
      else
        do i=is,ie ; pres(i) = 0.0 ; enddo
      endif
      do k=1,nz
        do i=is,ie
          d_pres(i) = (gv%g_Earth * gv%H_to_RZ) * h2d(i,k)
          p_lay(i) = pres(i) + 0.5*d_pres(i)
          pres(i) = pres(i) + d_pres(i)
        enddo
        call calculate_specific_vol_derivs(t2d(:,k), tv%S(:,j,k), p_lay(:), &
                 dsv_dt(:,j,k), dsv_ds(:,j,k), tv%eqn_of_state, eosdom)
        do i=is,ie ; dsv_dt_2d(i,k) = dsv_dt(i,j,k) ; enddo
      enddo
      pen_tke_2d(:,:) = 0.0
    endif

    ! Nothing more is done on this j-slice if there is no buoyancy forcing.
    if (.not.associated(fluxes%sw)) cycle

    if (nsw>0) call extract_optics_slice(optics, j, g, gv, opacity=opacityband, opacity_scale=(1.0/gv%m_to_H))

    ! The surface forcing is contained in the fluxes type.
    ! We aggregate the thermodynamic forcing for a time step into the following:
    ! netMassInOut = surface water fluxes [H ~> m or kg m-2] over time step
    !              = lprec + fprec + vprec + evap + lrunoff + frunoff
    !                note that lprec generally has sea ice melt/form included.
    ! netMassOut   = net mass leaving ocean surface [H ~> m or kg m-2] over a time step.
    !                netMassOut < 0 means mass leaves ocean.
    ! netHeat      = heat via surface fluxes [degC H ~> degC m or degC kg m-2], excluding the part
    !                contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0.
    ! netSalt      = surface salt fluxes [ppt H ~> dppt m or gSalt m-2]
    ! Pen_SW_bnd   = components to penetrative shortwave radiation split according to bands.
    !                This field provides that portion of SW from atmosphere that in fact
    !                enters to the ocean and participates in pentrative SW heating.
    ! nonpenSW     = non-downwelling SW flux, which is absorbed in ocean surface
    !                (in tandem w/ LW,SENS,LAT); saved only for diagnostic purposes.

    !----------------------------------------------------------------------------------------
    !BGR-June 26, 2017{
    !Temporary action to preserve answers while fixing a bug.
    ! To fix a bug in a diagnostic calculation, applyboundaryfluxesinout now returns
    !  the surface buoyancy flux. Previously, extractbuoyancyflux2d was called, meaning
    !  a second call to extractfluxes1d (causing the diagnostic net_heat to be incorrect).
    !  Note that this call to extract buoyancyflux2d was AFTER applyboundaryfluxesinout,
    !  which means it used the T/S fields after this routine.  Therefore, the surface
    !  buoyancy flux is computed here at the very end of this routine for legacy reasons.
    !  A few specific notes follow:
    !     1) The old method did not included river/calving contributions to heat flux.  This
    !        is kept consistent here via commenting code in the present extractFluxes1d <_rate>
    !        outputs, but we may reconsider this approach.
    !     2) The old method computed the buoyancy flux rate directly (by setting dt=1), instead
    !        of computing the integrated value (and dividing by dt). Hence the required
    !        additional outputs from extractFluxes1d.
    !          *** This is because: A*dt/dt =/=  A due to round off.
    !     3) The old method computed buoyancy flux after this routine, meaning the returned
    !        surface fluxes (from extractfluxes1d) must be recorded for use later in the code.
    !        We could (and maybe should) move that loop up to before the surface fluxes are
    !        applied, but this will change answers.
    !     For all these reasons we compute additional values of <_rate> which are preserved
    !     for the buoyancy flux calculation and reproduce the old answers.
    !   In the future this needs more detailed investigation to make sure everything is
    !   consistent and correct. These details shouldnt significantly effect climate,
    !   but do change answers.
    !-----------------------------------------------------------------------------------------
    if (calculate_buoyancy) then
      call extractfluxes1d(g, gv, us, fluxes, optics, nsw, j, dt,          &
                  h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
                  h2d, t2d, netmassinout, netmassout, netheat, netsalt,                   &
                  pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw,                &
                  net_heat_rate=netheat_rate, net_salt_rate=netsalt_rate,                 &
                  netmassinout_rate=netmassinout_rate, pen_sw_bnd_rate=pen_sw_bnd_rate)
    else
      call extractfluxes1d(g, gv, us, fluxes, optics, nsw, j, dt,          &
                  h_limit_fluxes, cs%use_river_heat_content, cs%use_calving_heat_content, &
                  h2d, t2d, netmassinout, netmassout, netheat, netsalt,                   &
                  pen_sw_bnd, tv, aggregate_fw_forcing, nonpensw=nonpensw)
   endif
    ! ea is for passive tracers
    do i=is,ie
    !  ea(i,j,1) = netMassInOut(i)
      if (aggregate_fw_forcing) then
        netmassout(i) = netmassinout(i)
        netmassin(i) = 0.
      else
        netmassin(i) = netmassinout(i) - netmassout(i)
      endif
      if (g%mask2dT(i,j)>0.0) then
        fluxes%netMassOut(i,j) = netmassout(i)
        fluxes%netMassIn(i,j) = netmassin(i)
      else
        fluxes%netMassOut(i,j) = 0.0
        fluxes%netMassIn(i,j) = 0.0
      endif
    enddo

    ! Apply the surface boundary fluxes in three steps:
    ! A/ update mass, temp, and salinity due to all terms except mass leaving
    !    ocean (and corresponding outward heat content), and ignoring penetrative SW.
    ! B/ update mass, salt, temp from mass leaving ocean.
    ! C/ update temp due to penetrative SW
    do i=is,ie
      if (g%mask2dT(i,j)>0.) then

        ! A/ Update mass, temp, and salinity due to incoming mass flux.
        do k=1,1

          ! Change in state due to forcing
          dthickness = netmassin(i) ! Since we are adding mass, we can use all of it
          dtemp = 0.
          dsalt = 0.

          ! Update the forcing by the part to be consumed within the present k-layer.
          ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
          netmassin(i) = netmassin(i) - dthickness
          ! This line accounts for the temperature of the mass exchange
          temp_in = t2d(i,k)
          salin_in = 0.0
          dtemp = dtemp + dthickness*temp_in

          ! Diagnostics of heat content associated with mass fluxes
          if (associated(fluxes%heat_content_massin))                             &
            fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) +   &
                         t2d(i,k) * max(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
          if (associated(fluxes%heat_content_massout))                            &
            fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
                         t2d(i,k) * min(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
          if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
                         t2d(i,k) * dthickness * gv%H_to_RZ

          ! Determine the energetics of river mixing before updating the state.
          if (calculate_energetics .and. associated(fluxes%lrunoff) .and. cs%do_rivermix) then
            ! Here we add an additional source of TKE to the mixed layer where river
            ! is present to simulate unresolved estuaries. The TKE input, TKE_river in
            ! [Z3 T-3 ~> m3 s-3], is diagnosed as follows:
            !   TKE_river = 0.5*rivermix_depth*g*(1/rho)*drho_ds*
            !               River*(Samb - Sriver) = CS%mstar*U_star^3
            ! where River is in units of [Z T-1 ~> m s-1].
            ! Samb = Ambient salinity at the mouth of the estuary
            ! rivermix_depth =  The prescribed depth over which to mix river inflow
            ! drho_ds = The gradient of density wrt salt at the ambient surface salinity.
            ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
            if (gv%Boussinesq) then
              rivermixconst = -0.5*(cs%rivermix_depth*dt) * ( us%L_to_Z**2*gv%g_Earth ) * gv%Rho0
            else
              rivermixconst = -0.5*(cs%rivermix_depth*dt) * gv%Rho0 * ( us%L_to_Z**2*gv%g_Earth )
            endif
            ctke(i,j,k) = ctke(i,j,k) + max(0.0, rivermixconst*dsv_ds(i,j,1) * &
                            (fluxes%lrunoff(i,j) + fluxes%frunoff(i,j)) * tv%S(i,j,1))
          endif

          ! Update state
          hold     = h2d(i,k)               ! Keep original thickness in hand
          h2d(i,k) = h2d(i,k) + dthickness  ! New thickness
          if (h2d(i,k) > 0.0) then
            if (calculate_energetics .and. (dthickness > 0.)) then
              ! Calculate the energy required to mix the newly added water over
              ! the topmost grid cell.
              ctke(i,j,k) = ctke(i,j,k) + 0.5*g_hconv2*(hold*dthickness) * &
                 ((t2d(i,k) - temp_in) * dsv_dt(i,j,k) + (tv%S(i,j,k) - salin_in) * dsv_ds(i,j,k))
            endif
            ithickness  = 1.0/h2d(i,k)      ! Inverse new thickness
            ! The "if"s below avoid changing T/S by roundoff unnecessarily
            if (dthickness /= 0. .or. dtemp /= 0.) t2d(i,k)    = (hold*t2d(i,k)    + dtemp)*ithickness
            if (dthickness /= 0. .or. dsalt /= 0.) tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness

          endif

        enddo ! k=1,1

        ! B/ Update mass, salt, temp from mass leaving ocean and other fluxes of heat and salt.
        do k=1,nz

          ! Place forcing into this layer if this layer has nontrivial thickness.
          ! For layers thin relative to 1/IforcingDepthScale, then distribute
          ! forcing into deeper layers.
          iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth - netmassout(i) )
          ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
          fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale)

          ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we
          ! limit the forcing applied to this cell, leaving the remaining forcing to
          ! be distributed downwards.
          if (-fractionofforcing*netmassout(i) > evap_cfl_limit*h2d(i,k)) then
            fractionofforcing = -evap_cfl_limit*h2d(i,k)/netmassout(i)
          endif

          ! Change in state due to forcing

          dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) )
          dtemp      = fractionofforcing*netheat(i)
          !   ### The 0.9999 here should become a run-time parameter?
          dsalt = max( fractionofforcing*netsalt(i), -0.9999*h2d(i,k)*tv%S(i,j,k))

          ! Update the forcing by the part to be consumed within the present k-layer.
          ! If fractionOfForcing = 1, then new netMassOut vanishes.
          netmassout(i) = netmassout(i) - dthickness
          netheat(i) = netheat(i) - dtemp
          netsalt(i) = netsalt(i) - dsalt

          ! This line accounts for the temperature of the mass exchange
          dtemp = dtemp + dthickness*t2d(i,k)

          ! Diagnostics of heat content associated with mass fluxes
          if (associated(fluxes%heat_content_massin)) &
            fluxes%heat_content_massin(i,j) = fluxes%heat_content_massin(i,j) + &
                         t2d(i,k) * max(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
          if (associated(fluxes%heat_content_massout)) &
            fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
                         t2d(i,k) * min(0.,dthickness) * gv%H_to_RZ * fluxes%C_p * idt
          if (associated(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
                         t2d(i,k) * dthickness * gv%H_to_RZ

          ! Update state by the appropriate increment.
          hold     = h2d(i,k)               ! Keep original thickness in hand
          h2d(i,k) = h2d(i,k) + dthickness  ! New thickness

          if (h2d(i,k) > 0.) then
            if (calculate_energetics) then
              ! Calculate the energy required to mix the newly added water over the topmost grid
              ! cell, assuming that the fluxes of heat and salt and rejected brine are initially
              ! applied in vanishingly thin layers at the top of the layer before being mixed
              ! throughout the layer.  Note that dThickness is always <= 0 here, and that
              ! negative cTKE is a deficit that will need to be filled later.
              ctke(i,j,k) = ctke(i,j,k) - (0.5*h2d(i,k)*g_hconv2) * &
                            ((dtemp - dthickness*t2d(i,k)) * dsv_dt(i,j,k) + &
                             (dsalt - dthickness*tv%S(i,j,k)) * dsv_ds(i,j,k))
            endif
            ithickness  = 1.0/h2d(i,k) ! Inverse of new thickness
            t2d(i,k)    = (hold*t2d(i,k) + dtemp)*ithickness
            tv%S(i,j,k) = (hold*tv%S(i,j,k) + dsalt)*ithickness
          elseif (h2d(i,k) < 0.0) then ! h2d==0 is a special limit that needs no extra handling
            call forcing_singlepointprint(fluxes,g,i,j,'applyBoundaryFluxesInOut (h<0)')
            write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
            write(0,*) 'applyBoundaryFluxesInOut(): netT,netS,netH=',netheat(i),netsalt(i),netmassinout(i)
            write(0,*) 'applyBoundaryFluxesInOut(): dT,dS,dH=',dtemp,dsalt,dthickness
            write(0,*) 'applyBoundaryFluxesInOut(): h(n),h(n+1),k=',hold,h2d(i,k),k
            call mom_error(fatal, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
                           "Complete mass loss in column!")
          endif

        enddo ! k

      ! Check if trying to apply fluxes over land points
      elseif ((abs(netheat(i))+abs(netsalt(i))+abs(netmassin(i))+abs(netmassout(i)))>0.) then

        if (.not. cs%ignore_fluxes_over_land) then
           call forcing_singlepointprint(fluxes,g,i,j,'applyBoundaryFluxesInOut (land)')
           write(0,*) 'applyBoundaryFluxesInOut(): lon,lat=',g%geoLonT(i,j),g%geoLatT(i,j)
           write(0,*) 'applyBoundaryFluxesInOut(): netHeat,netSalt,netMassIn,netMassOut=',&
                   netheat(i),netsalt(i),netmassin(i),netmassout(i)

           call mom_error(fatal, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
                                 "Mass loss over land?")
        endif

      endif

      ! If anything remains after the k-loop, then we have grounded out, which is a problem.
      if (netmassin(i)+netmassout(i) /= 0.0) then
!$OMP critical
        numberofgroundings = numberofgroundings +1
        if (numberofgroundings<=maxgroundings) then
          iground(numberofgroundings) = i ! Record i,j location of event for
          jground(numberofgroundings) = j ! warning message
          hgrounding(numberofgroundings) = netmassin(i)+netmassout(i)
        endif
!$OMP end critical
        if (cs%id_createdH>0) cs%createdH(i,j) = cs%createdH(i,j) - (netmassin(i)+netmassout(i))/dt
      endif

    enddo ! i

    ! Step C/ in the application of fluxes
    ! Heat by the convergence of penetrating SW.
    ! SW penetrative heating uses the updated thickness from above.

    ! Save temperature before increment with SW heating
    ! and initialize CS%penSWflux_diag to zero.
    if (cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0) then
      do k=1,nz ; do i=is,ie
        cs%penSW_diag(i,j,k)     = t2d(i,k)
        cs%penSWflux_diag(i,j,k) = 0.0
      enddo ; enddo
      k=nz+1 ; do i=is,ie
        cs%penSWflux_diag(i,j,k) = 0.0
      enddo
    endif

    if (calculate_energetics) then
      call absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h_limit_fluxes, &
                             .false., .true., t2d, pen_sw_bnd, tke=pen_tke_2d, dsv_dt=dsv_dt_2d)
      k = 1 ! For setting break-points.
      do k=1,nz ; do i=is,ie
        ctke(i,j,k) = ctke(i,j,k) + pen_tke_2d(i,k)
      enddo ; enddo
    else
      call absorbremainingsw(g, gv, us, h2d, opacityband, nsw, optics, j, dt, h_limit_fluxes, &
                             .false., .true., t2d, pen_sw_bnd)
    endif


    ! Step D/ copy updated thickness and temperature
    ! 2d slice now back into model state.
    do k=1,nz ; do i=is,ie
      h(i,j,k)    = h2d(i,k)
      tv%T(i,j,k) = t2d(i,k)
    enddo ; enddo

    ! Diagnose heating [Q R Z T-1 ~> W m-2] applied to a grid cell from SW penetration
    ! Also diagnose the penetrative SW heat flux at base of layer.
    if (cs%id_penSW_diag > 0 .or. cs%id_penSWflux_diag > 0) then

      ! convergence of SW into a layer
      do k=1,nz ; do i=is,ie
        cs%penSW_diag(i,j,k) = (t2d(i,k)-cs%penSW_diag(i,j,k))*h(i,j,k) * idt * tv%C_p * gv%H_to_RZ
      enddo ; enddo

      ! Perform a cumulative sum upwards from bottom to
      ! diagnose penetrative SW flux at base of tracer cell.
      ! CS%penSWflux_diag(i,j,k=1)    is penetrative shortwave at top of ocean.
      ! CS%penSWflux_diag(i,j,k=kbot+1) is zero, since assume no SW penetrates rock.
      ! CS%penSWflux_diag = rsdo  and CS%penSW_diag = rsdoabsorb
      ! rsdoabsorb(k) = rsdo(k) - rsdo(k+1), so that rsdo(k) = rsdo(k+1) + rsdoabsorb(k)
      if (cs%id_penSWflux_diag > 0) then
        do k=nz,1,-1 ; do i=is,ie
          cs%penSWflux_diag(i,j,k) = cs%penSW_diag(i,j,k) + cs%penSWflux_diag(i,j,k+1)
        enddo ; enddo
      endif

    endif

    ! Fill CS%nonpenSW_diag
    if (cs%id_nonpenSW_diag > 0) then
      do i=is,ie
        cs%nonpenSW_diag(i,j) = nonpensw(i) * idt * tv%C_p * gv%H_to_RZ
      enddo
    endif

    ! BGR: Get buoyancy flux to return for ePBL
    !  We want the rate, so we use the rate values returned from extractfluxes1d.
    !  Note that the *dt values could be divided by dt here, but
    !  1) Answers will change due to round-off
    !  2) Be sure to save their values BEFORE fluxes are used.
    if (calculate_buoyancy) then
      drhodt(:) = 0.0
      drhods(:) = 0.0
      netpen_rate(:) = 0.0
      ! Sum over bands and attenuate as a function of depth.
      ! netPen_rate is the netSW as a function of depth, but only the surface value is used here,
      ! in which case the values of dt, h, optics and H_limit_fluxes are irrelevant.  Consider
      ! writing a shorter and simpler variant to handle this very limited case.
      ! call sumSWoverBands(G, GV, US, h2d(:,:), optics_nbands(optics), optics, j, dt, &
      !                     H_limit_fluxes, .true., pen_SW_bnd_rate, netPen)
      do i=is,ie ; do nb=1,nsw ; netpen_rate(i) = netpen_rate(i) + pen_sw_bnd_rate(nb,i) ; enddo ; enddo

      ! Density derivatives
      if (associated(tv%p_surf)) then ; do i=is,ie ; surfpressure(i) = tv%p_surf(i,j) ; enddo ; endif
      call calculate_density_derivs(t2d(:,1), tv%S(:,j,1), surfpressure, drhodt, drhods, &
                                    tv%eqn_of_state, eosdom)
      ! 1. Adjust netSalt to reflect dilution effect of FW flux
      ! 2. Add in the SW heating for purposes of calculating the net
      ! surface buoyancy flux affecting the top layer.
      ! 3. Convert to a buoyancy flux, excluding penetrating SW heating
      !    BGR-Jul 5, 2017: The contribution of SW heating here needs investigated for ePBL.
      do i=is,ie
        skinbuoyflux(i,j) = - gorho * gv%H_to_Z * &
            (drhods(i) * (netsalt_rate(i) - tv%S(i,j,1)*netmassinout_rate(i)) + &
             drhodt(i) * ( netheat_rate(i) + netpen_rate(i)) ) ! [Z2 T-3 ~> m2 s-3]
      enddo
    endif

  enddo ! j-loop finish

  ! Post the diagnostics
  if (cs%id_createdH       > 0) call post_data(cs%id_createdH      , cs%createdH      , cs%diag)
  if (cs%id_penSW_diag     > 0) call post_data(cs%id_penSW_diag    , cs%penSW_diag    , cs%diag)
  if (cs%id_penSWflux_diag > 0) call post_data(cs%id_penSWflux_diag, cs%penSWflux_diag, cs%diag)
  if (cs%id_nonpenSW_diag  > 0) call post_data(cs%id_nonpenSW_diag , cs%nonpenSW_diag , cs%diag)

! The following check will be ignored if ignore_fluxes_over_land = true
  if (numberofgroundings>0 .and. .not. cs%ignore_fluxes_over_land) then
    do i = 1, min(numberofgroundings, maxgroundings)
      call forcing_singlepointprint(fluxes,g,iground(i),jground(i),'applyBoundaryFluxesInOut (grounding)')
      write(mesg(1:45),'(3es15.3)') g%geoLonT( iground(i), jground(i) ), &
                             g%geoLatT( iground(i), jground(i)), hgrounding(i)*gv%H_to_m
      call mom_error(warning, "MOM_diabatic_driver.F90, applyBoundaryFluxesInOut(): "//&
                              "Mass created. x,y,dh= "//trim(mesg), all_print=.true.)
    enddo

    if (numberofgroundings - maxgroundings > 0) then
      write(mesg, '(i4)') numberofgroundings - maxgroundings
      call mom_error(warning, "MOM_diabatic_driver:F90, applyBoundaryFluxesInOut(): "//&
                              trim(mesg) // " groundings remaining")
    endif
  endif

end subroutine applyboundaryfluxesinout

!> This subroutine initializes the parameters and control structure of the diabatic_aux module.
subroutine diabatic_aux_init(Time, G, GV, US, param_file, diag, CS, useALEalgorithm, use_ePBL)
  type(time_type), target, intent(in)    :: time !< The current model time.
  type(ocean_grid_type),   intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type), intent(in)    :: gv   !< The ocean's vertical grid structure
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(diag_ctrl), target, intent(inout) :: diag !< A structure used to regulate diagnostic output
  type(diabatic_aux_cs),   pointer       :: cs   !< A pointer to the control structure for the
                                                 !! diabatic_aux module, which is initialized here.
  logical,                 intent(in)    :: usealealgorithm !< If true, use the ALE algorithm rather
                                                 !! than layered mode.
  logical,                 intent(in)    :: use_epbl !< If true, use the implicit energetics planetary
                                                 !! boundary layer scheme to determine the diffusivity
                                                 !! in the surface boundary layer.

! This "include" declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl  = "MOM_diabatic_aux" ! This module's name.
  character(len=48)  :: thickness_units
  character(len=200) :: inputdir   ! The directory where NetCDF input files
  character(len=240) :: chl_filename ! A file from which chl_a concentrations are to be read.
  character(len=128) :: chl_file ! Data containing chl_a concentrations. Used
                                 ! when var_pen_sw is defined and reading from file.
  character(len=32)  :: chl_varname ! Name of chl_a variable in chl_file.
  logical :: use_temperature     ! True if thermodynamics are enabled.
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz, nbands
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

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

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

! Set default, read and log parameters
  call log_version(param_file, mdl, version, &
                   "The following parameters are used for auxiliary diabatic processes.")

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

  call get_param(param_file, mdl, "RECLAIM_FRAZIL", cs%reclaim_frazil, &
                 "If true, try to use any frazil heat deficit to cool any "//&
                 "overlying layers down to the freezing point, thereby "//&
                 "avoiding the creation of thin ice when the SST is above "//&
                 "the freezing point.", default=.true.)
  call get_param(param_file, mdl, "PRESSURE_DEPENDENT_FRAZIL", &
                                cs%pressure_dependent_frazil, &
                 "If true, use a pressure dependent freezing temperature "//&
                 "when making frazil. The default is false, which will be "//&
                 "faster but is inappropriate with ice-shelf cavities.", &
                 default=.false.)

  if (use_epbl) then
    call get_param(param_file, mdl, "IGNORE_FLUXES_OVER_LAND", cs%ignore_fluxes_over_land,&
         "If true, the model does not check if fluxes are being applied "//&
         "over land points. This is needed when the ocean is coupled "//&
         "with ice shelves and sea ice, since the sea ice mask needs to "//&
         "be different than the ocean mask to avoid sea ice formation "//&
         "under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
    call get_param(param_file, mdl, "DO_RIVERMIX", cs%do_rivermix, &
                 "If true, apply additional mixing wherever there is "//&
                 "runoff, so that it is mixed down to RIVERMIX_DEPTH "//&
                 "if the ocean is that deep.", default=.false.)
    if (cs%do_rivermix) &
      call get_param(param_file, mdl, "RIVERMIX_DEPTH", cs%rivermix_depth, &
                 "The depth to which rivers are mixed if DO_RIVERMIX is "//&
                 "defined.", units="m", default=0.0, scale=us%m_to_Z)
  else
    cs%do_rivermix = .false. ; cs%rivermix_depth = 0.0 ; cs%ignore_fluxes_over_land = .false.
  endif

  if (gv%nkml == 0) then
    call get_param(param_file, mdl, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
                   "If true, use the fluxes%runoff_Hflx field to set the "//&
                   "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
                   default=.false.)
    call get_param(param_file, mdl, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
                   "If true, use the fluxes%calving_Hflx field to set the "//&
                   "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
                   default=.false.)
  else
    cs%use_river_heat_content = .false.
    cs%use_calving_heat_content = .false.
  endif

  if (usealealgorithm) then
    cs%id_createdH = register_diag_field('ocean_model',"created_H",diag%axesT1, &
        time, "The volume flux added to stop the ocean from drying out and becoming negative in depth", &
        "m s-1", conversion=gv%H_to_m*us%s_to_T)
    if (cs%id_createdH>0) allocate(cs%createdH(isd:ied,jsd:jed))

    ! diagnostic for heating of a grid cell from convergence of SW heat into the cell
    cs%id_penSW_diag = register_diag_field('ocean_model', 'rsdoabsorb',                     &
          diag%axesTL, time, 'Convergence of Penetrative Shortwave Flux in Sea Water Layer',&
          'W m-2', conversion=us%QRZ_T_to_W_m2, &
          standard_name='net_rate_of_absorption_of_shortwave_energy_in_ocean_layer', v_extensive=.true.)

    ! diagnostic for penetrative SW heat flux at top interface of tracer cell (nz+1 interfaces)
    ! k=1 gives penetrative SW at surface; SW(k=nz+1)=0 (no penetration through rock).
    cs%id_penSWflux_diag = register_diag_field('ocean_model', 'rsdo',                               &
          diag%axesTi, time, 'Downwelling Shortwave Flux in Sea Water at Grid Cell Upper Interface',&
          'W m-2', conversion=us%QRZ_T_to_W_m2, standard_name='downwelling_shortwave_flux_in_sea_water')

    ! need both arrays for the SW diagnostics (one for flux, one for convergence)
    if (cs%id_penSW_diag>0 .or. cs%id_penSWflux_diag>0) then
      allocate(cs%penSW_diag(isd:ied,jsd:jed,nz)) ; cs%penSW_diag(:,:,:) = 0.0
      allocate(cs%penSWflux_diag(isd:ied,jsd:jed,nz+1)) ; cs%penSWflux_diag(:,:,:) = 0.0
    endif

    ! diagnostic for non-downwelling SW radiation (i.e., SW absorbed at ocean surface)
    cs%id_nonpenSW_diag = register_diag_field('ocean_model', 'nonpenSW',                       &
          diag%axesT1, time,                                                                   &
          'Non-downwelling SW radiation (i.e., SW absorbed in ocean surface with LW,SENS,LAT)',&
          'W m-2', conversion=us%QRZ_T_to_W_m2, &
          standard_name='nondownwelling_shortwave_flux_in_sea_water')
    if (cs%id_nonpenSW_diag > 0) then
      allocate(cs%nonpenSW_diag(isd:ied,jsd:jed)) ; cs%nonpenSW_diag(:,:) = 0.0
    endif
  endif

  if (use_temperature) then
    call get_param(param_file, mdl, "VAR_PEN_SW", cs%var_pen_sw, &
                   "If true, use one of the CHL_A schemes specified by "//&
                   "OPACITY_SCHEME to determine the e-folding depth of "//&
                   "incoming short wave radiation.", default=.false.)
    if (cs%var_pen_sw) then

      call get_param(param_file, mdl, "CHL_FROM_FILE", cs%chl_from_file, &
                   "If true, chl_a is read from a file.", default=.true.)
      if (cs%chl_from_file) then
        call time_interp_external_init()

        call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
        call get_param(param_file, mdl, "CHL_FILE", chl_file, &
                   "CHL_FILE is the file containing chl_a concentrations in "//&
                   "the variable CHL_A. It is used when VAR_PEN_SW and "//&
                   "CHL_FROM_FILE are true.", fail_if_missing=.true.)
        chl_filename = trim(slasher(inputdir))//trim(chl_file)
        call log_param(param_file, mdl, "INPUTDIR/CHL_FILE", chl_filename)
        call get_param(param_file, mdl, "CHL_VARNAME", chl_varname, &
                   "Name of CHL_A variable in CHL_FILE.", default='CHL_A')
        cs%sbc_chl = init_external_field(chl_filename, trim(chl_varname), domain=g%Domain%mpp_domain)
      endif

      cs%id_chl = register_diag_field('ocean_model', 'Chl_opac', diag%axesT1, time, &
          'Surface chlorophyll A concentration used to find opacity', 'mg m-3')
    endif
  endif

  id_clock_uv_at_h = cpu_clock_id('(Ocean find_uv_at_h)', grain=clock_routine)
  id_clock_frazil  = cpu_clock_id('(Ocean frazil)', grain=clock_routine)

end subroutine diabatic_aux_init

!> This subroutine initializes the control structure and any related memory
!! for the diabatic_aux module.
subroutine diabatic_aux_end(CS)
  type(diabatic_aux_cs), pointer :: cs !< The control structure returned by a previous
                                       !! call to diabatic_aux_init; it is deallocated here.

  if (.not.associated(cs)) return

  if (cs%id_createdH       >0) deallocate(cs%createdH)
  if (cs%id_penSW_diag     >0) deallocate(cs%penSW_diag)
  if (cs%id_penSWflux_diag >0) deallocate(cs%penSWflux_diag)
  if (cs%id_nonpenSW_diag  >0) deallocate(cs%nonpenSW_diag)

  if (associated(cs)) deallocate(cs)

end subroutine diabatic_aux_end

!> \namespace mom_diabatic_aux
!!
!!    This module contains the subroutines that, along with the
!!  subroutines that it calls, implements diapycnal mass and momentum
!!  fluxes and a bulk mixed layer.  The diapycnal diffusion can be
!!  used without the bulk mixed layer.
!!
!!    diabatic first determines the (diffusive) diapycnal mass fluxes
!!  based on the convergence of the buoyancy fluxes within each layer.
!!  The dual-stream entrainment scheme of MacDougall and Dewar (JPO,
!!  1997) is used for combined diapycnal advection and diffusion,
!!  calculated implicitly and potentially with the Richardson number
!!  dependent mixing, as described by Hallberg (MWR, 2000). 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 estimated flux in the
!!  layer below.  This flux is subject to the following 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 thick-
!!  ness.  If there is a bulk mixed layer, the buffer layer is treat-
!!  ed as a fixed density layer with vanishingly small diffusivity.
!!
!!    diabatic takes 5 arguments:  the two velocities (u and v), the
!!  thicknesses (h), a structure containing the forcing fields, and
!!  the length of time over which to act (dt).  The velocities and
!!  thickness are taken as inputs and modified within the subroutine.
!!  There is no limit on the time step.

end module mom_diabatic_aux