MOM_diabatic_aux.F90

Go to the documentation of this file.

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

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, April 1994 - July 2000                         *
!*     Alistair Adcroft, and Stephen Griffies                          *
!*                                                                     *
!*    This program contains the subroutine 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.                                *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v                                        *
!*    j    x ^ x ^ x   At >:  u                                        *
!*    j    > o > o >   At o:  h, T, S, buoy, ustar, ea, eb, etc.       *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_cpu_clock,     only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
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
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_version, param_file_type
use mom_forcing_type,  only : forcing, extractfluxes1d, forcing_singlepointprint
use mom_grid,          only : ocean_grid_type
use mom_io,            only : vardesc
use mom_shortwave_abs, only : absorbremainingsw, optics_type, sumswoverbands
use mom_variables,     only : thermo_var_ptrs, vertvisc_type! , accel_diag_ptrs
use mom_verticalgrid,  only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

public diabatic_aux_init, diabatic_aux_end
public make_frazil, adjust_salt, insert_brine, differential_diffuse_t_s, tridiagts
public find_uv_at_h, diagnosemldbydensitydifference, applyboundaryfluxesinout
!> 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 "rivermix_depth" meters
  real    :: rivermix_depth = 0.0  !< The depth to which rivers are mixed if
                                   !! do_rivermix = T, in 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.

  type(diag_ctrl), pointer :: diag !< Structure used to regulate timing of diagnostic output

  ! Diagnostic handles
  integer :: id_createdh       = -1
  integer :: id_brine_lay      = -1
  integer :: id_pensw_diag     = -1 !< Penetrative shortwave heating (flux convergence) diagnostic
  integer :: id_penswflux_diag = -1 !< Penetrative shortwave flux diagnostic
  integer :: id_nonpensw_diag  = -1 !< Non-penetrative shortwave heating diagnostic

  ! Optional diagnostic arrays
  real, allocatable, dimension(:,:)   :: createdh       !< The amount of volume added in order to avoid grounding (m/s)
  real, allocatable, dimension(:,:,:) :: pensw_diag     !< Heating in a layer from convergence of penetrative SW (W/m2)
  real, allocatable, dimension(:,:,:) :: penswflux_diag !< Penetrative SW flux at base of grid layer (W/m2)
  real, allocatable, dimension(:,:)   :: nonpensw_diag  !< Non-downwelling SW radiation (W/m2) at ocean surface

end type diabatic_aux_cs

integer :: id_clock_uv_at_h, id_clock_frazil

contains

subroutine make_frazil(h, tv, G, GV, CS, p_surf)
  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, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                 intent(inout) :: tv
  type(diabatic_aux_cs),                 intent(in)    :: CS
  real, dimension(SZI_(G),SZJ_(G)), intent(in), optional :: p_surf

!   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 J m-2) in tv%frazil.
!   The expression, below, for the freezing point of sea water comes
! from Millero (1978) via Appendix A of Gill, 1982.

! Arguments: h - Layer thickness, in m or kg m-2.
!  (in/out)  tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 diabatic_driver_init.
  real, dimension(SZI_(G)) :: &
    fraz_col, & ! The accumulated heat requirement due to frazil, in J.
    T_freeze, & ! The freezing potential temperature at the current salinity, C.
    ps          ! pressure
  real, dimension(SZI_(G),SZK_(G)) :: &
    pressure    ! The pressure at the middle of each layer in Pa.
  real :: hc    ! A layer's heat capacity in J m-2 K-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

  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
  endif
!$OMP parallel do default(none) shared(is,ie,js,je,CS,G,GV,h,nz,tv,p_surf) &
!$OMP                           private(fraz_col,T_fr_set,T_freeze,hc,ps,pressure)
  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*gv%H_to_Pa)*h(i,j,1)
      enddo
      do k=2,nz ; do i=is,ie
        pressure(i,k) = pressure(i,k-1) + &
          (0.5*gv%H_to_Pa) * (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)
          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_kg_m2) * 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)
            t_fr_set = .true.
          endif

          hc = (tv%C_p*gv%H_to_kg_m2) * h(i,j,k)
          if (h(i,j,k) <= 10.0*gv%Angstrom) 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
          else
            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

subroutine differential_diffuse_t_s(h, tv, visc, 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, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                 intent(inout) :: tv
  type(vertvisc_type),                   intent(in)    :: visc
  real,                                  intent(in)    :: dt

! This subroutine applies double diffusion to T & S, assuming no diapycal mass
! fluxes, using a simple triadiagonal solver.

! Arguments: h - Layer thickness, in m or kg m-2.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      visc - A structure containing vertical viscosities, bottom boundary
!                   layer properies, and related fields.
!  (in)      dt - Time increment, in s.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.

  real, dimension(SZI_(G)) :: &
    b1_T, b1_S, &  !  Variables used by the tridiagonal solvers of T & S, in H.
    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, in m or kg m-2.
  real, dimension(SZI_(G),SZK_(G)+1) :: &
    mix_T, mix_S   !  Mixing distances in both directions across each
                   !  interface, in m or kg m-2.
  real :: h_tr         ! h_tr is h at tracer points with a tiny thickness
                       ! added to ensure positive definiteness, in m or kg m-2.
  real :: h_neglect    ! A thickness that is so small it is usually lost
                       ! in roundoff and can be neglected, in m or kg m-2.
  real :: I_h_int      ! The inverse of the thickness associated with an
                       ! interface, in 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 in m or kg m-2.

  integer :: i, j, k, is, ie, js, je, nz
  real, pointer :: T(:,:,:), S(:,:,:), Kd_T(:,:,:), Kd_S(:,:,:)
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  h_neglect = gv%H_subroundoff

  if (.not.associated(tv%T)) call mom_error(fatal, &
      "differential_diffuse_T_S: Called with an unassociated tv%T")
  if (.not.associated(tv%S)) call mom_error(fatal, &
      "differential_diffuse_T_S: Called with an unassociated tv%S")
  if (.not.associated(visc%Kd_extra_T)) call mom_error(fatal, &
      "differential_diffuse_T_S: Called with an unassociated visc%Kd_extra_T")
  if (.not.associated(visc%Kd_extra_S)) call mom_error(fatal, &
      "differential_diffuse_T_S: Called with an unassociated visc%Kd_extra_S")

  t => tv%T ; s => tv%S
  kd_t => visc%Kd_extra_T ; kd_s => visc%Kd_extra_S
!$OMP parallel do default(none) shared(is,ie,js,je,h,h_neglect,dt,Kd_T,Kd_S,G,GV,T,S,nz) &
!$OMP                          private(I_h_int,mix_T,mix_S,h_tr,b1_T,b1_S, &
!$OMP                                  d1_T,d1_S,c1_T,c1_S,b_denom_T,b_denom_S)
  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%m_to_H**2) * i_h_int
      mix_s(i,2) = ((dt * kd_s(i,j,2)) * gv%m_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%m_to_H**2) * i_h_int
      mix_s(i,k+1) = ((dt * kd_s(i,j,k+1)) * gv%m_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

subroutine adjust_salt(h, tv, G, GV, CS)
  type(ocean_grid_type),                 intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type),               intent(in)    :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                 intent(inout) :: tv
  type(diabatic_aux_cs),                 intent(in)    :: CS

!  Keep salinity from falling below a small but positive threshold
!  This occurs when the ice model attempts to extract more salt then
!  is actually available to it from the ocean.

! Arguments: h - Layer thickness, in m.
!  (in/out)  tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 diabatic_driver_init.
  real :: salt_add_col(szi_(g),szj_(g)) ! The accumulated salt requirement
  real :: S_min      ! The minimum salinity
  real :: mc         ! A layer's mass 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

!  call cpu_clock_begin(id_clock_adjust_salt)

  s_min = 0.01

  salt_add_col(:,:) = 0.0

!$OMP parallel do default(none) shared(is,ie,js,je,nz,G,GV,tv,h,salt_add_col, S_min) &
!$OMP                          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_kg_m2 * h(i,j,k)
        if (h(i,j,k) <= 10.0*gv%Angstrom) 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
        else
          if (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
      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

subroutine insert_brine(h, tv, G, GV, fluxes, nkmb, CS, dt, id_brine_lay)
  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, in H (usually m or kg m-2)
  type(thermo_var_ptrs),                 intent(inout) :: tv
  type(forcing),                         intent(in)    :: fluxes
  integer,                               intent(in)    :: nkmb
  type(diabatic_aux_cs),                 intent(in)    :: CS
  real,                                  intent(in)    :: dt
  integer,                               intent(in)    :: id_brine_lay

! Insert salt from brine rejection into the first layer below
! the mixed layer which both contains mass and in which the
! change in layer density remains stable after the addition
! of salt via brine rejection.

! Arguments: h - Layer thickness, in m.
!  (in/out)  tv - A structure containing pointers to any available
!                 thermodynamic fields. Absent fields have NULL ptrs.
!  (in)      fluxes = A structure containing pointers to any possible
!                     forcing fields; unused fields have NULL ptrs.
!  (in)      nkmb - The number of layers in the mixed and buffer layers.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 diabatic_driver_init.

  real :: salt(szi_(g)) ! The amount of salt rejected from
                        !  sea ice. [grams]
  real :: dzbr(szi_(g)) ! cumulative depth over which brine is distributed
  real :: inject_layer(szi_(g),szj_(g)) ! diagnostic

  real :: p_ref_cv(szi_(g))
  real :: T(szi_(g),szk_(g))
  real :: S(szi_(g),szk_(g))
  real :: h_2d(szi_(g),szk_(g))
  real :: Rcv(szi_(g),szk_(g))
  real :: mc  ! A layer's mass in kg m-2 .
  real :: s_new,R_new,t0,scale, cdz
  integer :: i, j, k, is, ie, js, je, nz, ks

  real, parameter :: brine_dz = 1.0  ! minumum thickness over which to distribute brine
  real, parameter :: s_max = 45.0    ! salinity bound

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

  if (.not.ASSOCIATED(fluxes%salt_flux)) return

  p_ref_cv(:)  = tv%P_ref

  inject_layer = nz

  do j=js,je

    salt(:)=0.0 ; dzbr(:)=0.0

    do i=is,ie ; if (g%mask2dT(i,j) > 0.) then
      salt(i) = dt * (1000. * fluxes%salt_flux(i,j))
    endif ; enddo

    do k=1,nz
      do i=is,ie
        t(i,k)=tv%T(i,j,k); s(i,k)=tv%S(i,j,k)
        ! avoid very small thickness
        h_2d(i,k)=max(h(i,j,k), gv%Angstrom)
      enddo

      call calculate_density(t(:,k), s(:,k), p_ref_cv, rcv(:,k), is, &
                             ie-is+1, tv%eqn_of_state)
    enddo

    ! First, try to find an interior layer where inserting all the salt
    ! will not cause the layer to become statically unstable.
    ! Bias towards deeper layers.

    do k=nkmb+1,nz-1 ; do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
        mc = gv%H_to_kg_m2 * h_2d(i,k)
        s_new = s(i,k) + salt(i)/mc
        t0 = t(i,k)
        call calculate_density(t0,s_new,tv%P_Ref,r_new,tv%eqn_of_state)
        if (r_new < 0.5*(rcv(i,k)+rcv(i,k+1)) .and. s_new<s_max) then
          dzbr(i)=dzbr(i)+h_2d(i,k)
          inject_layer(i,j) = min(inject_layer(i,j),real(k))
        endif
      endif
    enddo ; enddo

    ! Then try to insert into buffer layers if they exist
    do k=nkmb,gv%nkml+1,-1 ; do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
        mc = gv%H_to_kg_m2 * h_2d(i,k)
        dzbr(i)=dzbr(i)+h_2d(i,k)
        inject_layer(i,j) = min(inject_layer(i,j),real(k))
      endif
    enddo ; enddo

    ! finally if unable to find a layer to insert, then place in mixed layer

    do k=1,gv%nkml ; do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .and. dzbr(i) < brine_dz .and. salt(i) > 0.) then
        mc = gv%H_to_kg_m2 * h_2d(i,k)
        dzbr(i)=dzbr(i)+h_2d(i,k)
        inject_layer(i,j) = min(inject_layer(i,j),real(k))
      endif
    enddo ; enddo


    do i=is,ie
      if ((g%mask2dT(i,j) > 0.0) .and. salt(i) > 0.) then
    !   if (dzbr(i)< brine_dz) call MOM_error(FATAL,"insert_brine: failed")
        ks=inject_layer(i,j)
        cdz=0.0
        do k=ks,nz
          mc = gv%H_to_kg_m2 * h_2d(i,k)
          scale = h_2d(i,k)/dzbr(i)
          cdz=cdz+h_2d(i,k)
          if (cdz > 1.0) exit
          tv%S(i,j,k) = tv%S(i,j,k) + scale*salt(i)/mc
        enddo
      endif
    enddo

  enddo

  if (cs%id_brine_lay > 0) call post_data(cs%id_brine_lay,inject_layer,cs%diag)

end subroutine insert_brine

subroutine tridiagts(G, GV, is, ie, js, je, hold, ea, eb, T, S)
! Simple tri-diagnonal solver for T and S
! "Simple" means it only uses arrays hold, ea and eb
  ! Arguments
  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, ie, js, je
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)    :: hold, ea, eb
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: T, S
  ! 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(none) shared(is,ie,js,je,G,GV,hold,eb,T,S,ea) &
!$OMP                          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


subroutine find_uv_at_h(u, v, h, u_h, v_h, G, GV, 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
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)  :: u    !< The zonal velocity, in m s-1
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)  :: v    !< The meridional velocity, in m s-1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)  :: h    !< Layer thicknesses, in H (usually m or kg m-2)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(out) :: u_h, v_h
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in), optional  :: ea, eb
!   This subroutine calculates u_h and v_h (velocities at thickness
! points), optionally using the entrainments (in m) passed in as arguments.

! Arguments: u - Zonal velocity, in m s-1.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m or kg m-2.
!  (out)     u_h - The zonal velocity at thickness points after
!                  entrainment, in m s-1.
!  (out)     v_h - The meridional velocity at thickness points after
!                  entrainment, in m s-1.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in, opt) ea - The amount of fluid entrained from the layer above within
!                 this time step, in units of m or kg m-2.  Omitting ea is the
!                 same as setting it to 0.
!  (in, opt) eb - The amount of fluid entrained from the layer below within
!                 this time step, in units of m or kg m-2.  Omitting eb is the
!                 same as setting it to 0.  ea and eb must either be both
!                 present or both absent.

  real :: b_denom_1    ! The first term in the denominator of b1 in m or kg m-2.
  real :: h_neglect    ! A thickness that is so small it is usually lost
                       ! in roundoff and can be neglected, in 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 :: s, 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 = g%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(s,Idenom,a_w,a_e,a_s,a_n,b_denom_1,b1,d1,c1)
  do j=js,je
    do i=is,ie
      s = g%areaCu(i-1,j)+g%areaCu(i,j)
      if (s>0.0) then
        idenom = sqrt(0.5*g%IareaT(i,j)/s)
        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

      s = g%areaCv(i,j-1)+g%areaCv(i,j)
      if (s>0.0) then
        idenom = sqrt(0.5*g%IareaT(i,j)/s)
        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


!> 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, diagPtr, id_N2subML, id_MLDsq)
  type(ocean_grid_type),                 intent(in) :: G           !< Grid type
  type(verticalgrid_type),               intent(in) :: GV          !< ocean vertical grid structure
  integer,                               intent(in) :: id_MLD      !< Handle (ID) of MLD diagnostic
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h        !< Layer thickness
  type(thermo_var_ptrs),                 intent(in) :: tv          !< Thermodynamics type
  real,                                  intent(in) :: densityDiff !< Density difference to determine MLD (kg/m3)
  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

  ! Local variables
  real, dimension(SZI_(G))          :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK, dK, dKm1, pRef_MLD
  real, dimension(SZI_(G))          :: rhoAtK, rho1, d1, pRef_N2 ! Used for N2
  real, dimension(SZI_(G), SZJ_(G)) :: MLD ! Diagnosed mixed layer depth
  real, dimension(SZI_(G), SZJ_(G)) :: subMLN2 ! Diagnosed stratification below ML
  real, dimension(SZI_(G), SZJ_(G)) :: MLD2 ! Diagnosed MLD^2
  real, parameter                   :: dz_subML = 50. ! Depth below ML over which to diagnose stratification (m)
  integer :: i, j, is, ie, js, je, k, nz, id_N2, id_SQ
  real :: aFac, ddRho

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

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

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  pref_mld(:) = 0. ; pref_n2(:) = 0.
  do j=js,je
    do i=is,ie ; dk(i) = 0.5 * h(i,j,1) * gv%H_to_m ; enddo ! Depth of center of surface layer
    call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pref_mld, rhosurf, is, ie-is+1, tv%eqn_of_state)
    do i=is,ie
      deltarhoatk(i) = 0.
      mld(i,j) = 0.
      if (id_n2>0) then
        submln2(i,j) = 0.
        rho1(i) = 0.
        d1(i) = 0.
        pref_n2(i) = gv%g_Earth * gv%Rho0 * h(i,j,1) * gv%H_to_m ! Boussinesq approximation!!!! ?????
        !### This should be: pRef_N2(i) = GV%g_Earth * GV%H_to_kg_m2 * h(i,j,1) ! This might change answers at roundoff.
      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_m ! Depth of center of layer K
      enddo

      ! Stratification, N2, immediately below the mixed layer, averaged over at least 50 m.
      if (id_n2>0) then
        do i=is,ie
          pref_n2(i) = pref_n2(i) + gv%g_Earth * gv%Rho0 * h(i,j,k) * gv%H_to_m ! Boussinesq approximation!!!! ?????
          !### This should be: pRef_N2(i) = pRev_N2(i) + GV%g_Earth * GV%H_to_kg_m2 * h(i,j,k) ! This might change answers at roundoff.
        enddo
        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_n2, rhoatk, is, ie-is+1, tv%eqn_of_state)
        do i=is,ie
          if (mld(i,j)>0. .and. submln2(i,j)==0.) then ! This block is below the mixed layer
            if (d1(i)==0.) then ! Record the density, depth and pressure, immediately below the ML
              rho1(i) = rhoatk(i)
              d1(i) = dk(i)
              !### It looks to me like there is bad logic here. - RWH
              ! Use pressure at the bottom of the upper layer used in calculating d/dz rho
              pref_n2(i) = pref_n2(i) + gv%g_Earth * gv%Rho0 * h(i,j,k) * gv%H_to_m ! Boussinesq approximation!!!! ?????
              !### This line should be: pRef_N2(i) = pRev_N2(i) + GV%g_Earth * GV%H_to_kg_m2 * h(i,j,k) ! This might change answers at roundoff.
            endif
            if (d1(i)>0. .and. dk(i)-d1(i)>=dz_subml) then
              submln2(i,j) = gv%g_Earth/ gv%Rho0 * (rho1(i)-rhoatk(i)) / (d1(i) - dk(i))
            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, is, ie-is+1, tv%eqn_of_state)
      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
   !  if (id_N2>0 .and. subMLN2(i,j)==0. .and. d1(i)>0. .and. dK(i)-d1(i)>0.) then
   !    ! Use what ever stratification we can, measured over what ever distance is available
   !    subMLN2(i,j) = GV%g_Earth/ GV%Rho0 * (rho1(i)-rhoAtK(i)) / (d1(i) - dK(i))
   !  endif
    enddo
  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

!> 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, dt, fluxes, optics, 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
  real,                                  intent(in)    :: dt !< Time-step over which forcing is applied (s)
  type(forcing),                         intent(inout) :: fluxes !< Surface fluxes container
  type(optics_type),                     pointer       :: optics !< Optical properties container
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h  !< Layer thickness in H units
  type(thermo_var_ptrs),                 intent(inout) :: tv !< Thermodynamics container
  !> If False, treat in/out fluxes separately.
  logical,                               intent(in)    :: aggregate_FW_forcing
  !> The largest fraction of a layer that can be evaporated in one time-step (non-dim).
  real,                                  intent(in)   :: evap_CFL_limit
  !> The smallest depth over which heat and freshwater fluxes is applied, in m.
  real,                                  intent(in)   :: minimum_forcing_depth
  !> Turbulent kinetic energy requirement to mix forcing through each layer, in W m-2
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: cTKE
  !> Partial derivative of specific volume with potential temperature, in m3 kg-1 K-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: dSV_dT
  !> Partial derivative of specific a volume with potential salinity, in m3 kg-1 / (g kg-1).
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), optional, intent(out) :: dSV_dS
  !> Buoyancy flux at surface in m2 s-3
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(out) :: SkinBuoyFlux

  ! Local variables
  integer, parameter :: maxGroundings = 5
  integer :: numberOfGroundings, iGround(maxgroundings), jGround(maxgroundings)
  real :: H_limit_fluxes, IforcingDepthScale, Idt
  real :: dThickness, dTemp, dSalt
  real :: fractionOfForcing, hOld, Ithickness
  real :: RivermixConst  ! A constant used in implementing river mixing, in Pa s.
  real, dimension(SZI_(G)) :: &
    d_pres,       &  ! pressure change across a layer (Pa)
    p_lay,        &  ! average pressure in a layer (Pa)
    pres,         &  ! pressure at an interface (Pa)
    netMassInOut, &  ! surface water fluxes (H units) over time step
    netMassIn,    &  ! mass entering ocean surface (H units) over a time step
    netMassOut,   &  ! mass leaving ocean surface (H units) over a time step
    netHeat,      &  ! heat (degC * H) via surface fluxes, excluding
                     ! Pen_SW_bnd and netMassOut
    netsalt,      &  ! surface salt flux ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
    nonpensw,     &  ! non-downwelling SW, which is absorbed at ocean surface
    surfpressure, &  ! Surface pressure (approximated as 0.0)
    drhodt,       &  ! change in density per change in temperature
    drhods,       &  ! change in density per change in salinity
    netheat_rate, &  ! netheat but for dt=1 (e.g. returns a rate)
    netsalt_rate, &  ! netsalt but for dt=1 (e.g. returns a rate)
    netmassinout_rate! netmassinout but for dt=1 (e.g. returns a rate)
  real, dimension(SZI_(G), SZK_(G))                     :: h2d, T2d
  real, dimension(SZI_(G), SZK_(G))                     :: pen_TKE_2d, dSV_dT_2d
  real, dimension(SZI_(G),SZK_(G)+1)                    :: netPen
  real, dimension(max(optics%nbands,1),SZI_(G))         :: Pen_SW_bnd, Pen_SW_bnd_rate
                                                           !^ _rate is w/ dt=1
  real, dimension(max(optics%nbands,1),SZI_(G),SZK_(G)) :: opacityBand
  real                                                  :: hGrounding(maxgroundings)
  real    :: Temp_in, Salin_in
  real    :: I_G_Earth, g_Hconv2
  real    :: GoRho
  logical :: calculate_energetics
  logical :: calculate_buoyancy
  integer :: i, j, is, ie, js, je, k, nz, n, nsw
  integer :: start, npts
  character(len=45) :: mesg

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

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

#define _OLD_ALG_
  nsw = optics%nbands
  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
  i_g_earth = 1.0 / gv%g_Earth
  g_hconv2 = gv%g_Earth * gv%H_to_kg_m2**2

  if (present(ctke)) ctke(:,:,:) = 0.0
  if (calculate_buoyancy) then
    surfpressure(:) = 0.0
    gorho       = gv%g_Earth / gv%Rho0
    start       = 1 + g%isc - g%isd
    npts        = 1 + g%iec - g%isc
  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, 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,optics,fluxes,dt,    &
!$OMP                                  H_limit_fluxes,IforcingDepthScale,                &
!$OMP                                  numberOfGroundings,iGround,jGround,nonPenSW,      &
!$OMP                                  hGrounding,CS,Idt,aggregate_FW_forcing,           &
!$OMP                                  minimum_forcing_depth,evap_CFL_limit,             &
!$OMP                                  calculate_buoyancy,netPen,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                                  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,SurfPressure,       &
!$OMP                                  start,npts,                                       &
!$OMP                                  pen_TKE_2d,Temp_in,Salin_in,RivermixConst)

  ! Work in vertical slices for efficiency
  do j=js,je

    ! 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)
      do n=1,nsw
        opacityband(n,i,k) = (1.0 / gv%m_to_H)*optics%opacity_band(n,i,j,k)
      enddo
    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.
      do i=is,ie ; pres(i) = 0.0 ; enddo ! Add surface pressure?
      do k=1,nz
        do i=is,ie
          d_pres(i) = gv%g_Earth * gv%H_to_kg_m2 * 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), is, ie-is+1, tv%eqn_of_state)
        do i=is,ie ; dsv_dt_2d(i,k) = dsv_dt(i,j,k) ; enddo
!        do i=is,ie
!          dT_to_dPE(i,k) = I_G_Earth * d_pres(i) * p_lay(i) * dSV_dT(i,j,k)
!          dS_to_dPE(i,k) = I_G_Earth * d_pres(i) * p_lay(i) * dSV_dS(i,j,k)
!        enddo
      enddo
      pen_tke_2d(:,:) = 0.0
    endif

    ! 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 units) 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 units) over a time step.
    !                netMassOut < 0 means mass leaves ocean.
    ! netHeat      = heat (degC * H) via surface fluxes, excluding the part
    !                contained in Pen_SW_bnd; and excluding heat_content of netMassOut < 0.
    ! netSalt      = surface salt fluxes ( g(salt)/m2 for non-Bouss and ppt*H for Bouss )
    ! 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, 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, 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_kg_m2 * 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_kg_m2 * fluxes%C_p * idt
          if (ASSOCIATED(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
                         t2d(i,k) * dthickness * gv%H_to_kg_m2

          ! 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 is diagnosed
            ! as follows:
            !   TKE_river[m3 s-3] = 0.5*rivermix_depth*g*(1/rho)*drho_ds*
            !                       River*(Samb - Sriver) = CS%mstar*U_star^3
            ! where River is in units of m s-1.
            ! Samb = Ambient salinity at the mouth of the estuary
            ! rivermix_depth =  The prescribed depth over which to mix river inflow
            ! drho_ds = The gradient of density wrt salt at the ambient surface salinity.
            ! Sriver = 0 (i.e. rivers are assumed to be pure freshwater)
            rivermixconst = -0.5*(cs%rivermix_depth*dt)*gv%m_to_H*gv%H_to_Pa

            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.  ###CHECK THE SIGNS!!!
              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*gv%m_to_H - 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) +   &
                         tv%T(i,j,k) * max(0.,dthickness) * gv%H_to_kg_m2 * fluxes%C_p * idt
          if (ASSOCIATED(fluxes%heat_content_massout))                            &
            fluxes%heat_content_massout(i,j) = fluxes%heat_content_massout(i,j) + &
                         tv%T(i,j,k) * min(0.,dthickness) * gv%H_to_kg_m2 * fluxes%C_p * idt
          if (ASSOCIATED(tv%TempxPmE)) tv%TempxPmE(i,j) = tv%TempxPmE(i,j) + &
                         tv%T(i,j,k) * dthickness * gv%H_to_kg_m2
!NOTE tv%T should be T2d

          ! 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. ###CHECK THE SIGNS!!!
              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, h2d, opacityband, nsw, 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, h2d, opacityband, nsw, 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 (W/m2) 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_kg_m2
      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)
      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(:,:) = 0.0
       ! Sum over bands and attenuate as a function of depth
       ! netPen is the netSW as a function of depth
       call sumswoverbands(g, gv, h2d(:,:), optics%opacity_band(:,:,j,:), nsw, j, dt, &
            h_limit_fluxes, .true., pen_sw_bnd_rate, netpen)
       ! Density derivatives
       call calculate_density_derivs(t2d(:,1), tv%S(:,j,1), surfpressure, &
            drhodt, drhods, start, npts, tv%eqn_of_state)
       ! 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.
       skinbuoyflux(g%isc:g%iec,j) = - gorho * ( drhods(g%isc:g%iec) * (netsalt_rate(g%isc:g%iec) &
            - tv%S(g%isc:g%iec,j,1) * netmassinout_rate(g%isc:g%iec)* gv%H_to_m )&
            + drhodt(g%isc:g%iec) * ( netheat_rate(g%isc:g%iec) +        &
            netpen(g%isc:g%iec,1))) * gv%H_to_m ! m^2/s^3
    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)
      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

subroutine diabatic_aux_init(Time, G, GV, param_file, diag, CS, useALEalgorithm, use_ePBL)
  type(time_type),         intent(in)    :: Time
  type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time parameters
  type(diag_ctrl), target, intent(inout) :: diag
  type(diabatic_aux_cs),   pointer       :: CS
  logical,                 intent(in)    :: useALEalgorithm
  logical,                 intent(in)    :: use_ePBL

! Arguments:
!  (in)     Time       = current model time
!  (in)     G          = ocean grid structure
!  (in)     GV - The ocean's vertical grid structure.
!  (in)     param_file = structure indicating the open file to parse for parameter values
!  (in)     diag       = structure used to regulate diagnostic output
!  (in/out) CS         = pointer set to point to the control structure for this module
!  (in)     use_ePBL   = If true, use the implicit energetics planetary boundary
!                        layer scheme to determine the diffusivity in the
!                        surface boundary layer.
  type(vardesc) :: vd

! This "include" declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mod  = "MOM_diabatic_aux" ! This module's name.
  character(len=48)  :: thickness_units
  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

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

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

  if (use_epbl) then
    call get_param(param_file, mod, "IGNORE_FLUXES_OVER_LAND", cs%ignore_fluxes_over_land,&
         "If true, the model does not check if fluxes are being applied\n"//&
         "over land points. This is needed when the ocean is coupled \n"//&
         "with ice shelves and sea ice, since the sea ice mask needs to \n"//&
         "be different than the ocean mask to avoid sea ice formation \n"//&
         "under ice shelves. This flag only works when use_ePBL = True.", default=.false.)
    call get_param(param_file, mod, "DO_RIVERMIX", cs%do_rivermix, &
                 "If true, apply additional mixing whereever there is \n"//&
                 "runoff, so that it is mixed down to RIVERMIX_DEPTH \n"//&
                 "if the ocean is that deep.", default=.false.)
    if (cs%do_rivermix) &
      call get_param(param_file, mod, "RIVERMIX_DEPTH", cs%rivermix_depth, &
                 "The depth to which rivers are mixed if DO_RIVERMIX is \n"//&
                 "defined.", units="m", default=0.0)
  else ; cs%do_rivermix = .false. ; cs%rivermix_depth = 0.0 ; endif
  if (gv%nkml == 0) then
    call get_param(param_file, mod, "USE_RIVER_HEAT_CONTENT", cs%use_river_heat_content, &
                   "If true, use the fluxes%runoff_Hflx field to set the \n"//&
                   "heat carried by runoff, instead of using SST*CP*liq_runoff.", &
                   default=.false.)
    call get_param(param_file, mod, "USE_CALVING_HEAT_CONTENT", cs%use_calving_heat_content, &
                   "If true, use the fluxes%calving_Hflx field to set the \n"//&
                   "heat carried by runoff, instead of using SST*CP*froz_runoff.", &
                   default=.false.)
  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", &
        "meter second-1")
    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',&
          'Watt meter-2', 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',&
          'Watt meter-2', 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)',&
          'Watt meter-2', 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

  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


subroutine diabatic_aux_end(CS)
  type(diabatic_aux_cs), pointer :: CS

  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

end module mom_diabatic_aux