MOM_diabatic_driver.F90

!> This routine drives the diabatic/dianeutral physics for MOM
module mom_diabatic_driver

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


use mom_bulk_mixed_layer,    only : bulkmixedlayer, bulkmixedlayer_init, bulkmixedlayer_cs
use mom_debugging,           only : hchksum
use mom_checksum_packages,   only : mom_state_chksum, mom_state_stats
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_cvmix_shear,         only : cvmix_shear_is_used
use mom_cvmix_ddiff,         only : cvmix_ddiff_is_used
use mom_diabatic_aux,        only : diabatic_aux_init, diabatic_aux_end, diabatic_aux_cs
use mom_diabatic_aux,        only : make_frazil, adjust_salt, differential_diffuse_t_s, tridiagts
use mom_diabatic_aux,        only : find_uv_at_h, diagnosemldbydensitydifference, applyboundaryfluxesinout
use mom_diabatic_aux,        only : diagnosemldbyenergy
use mom_diabatic_aux,        only : set_pen_shortwave
use mom_diag_mediator,       only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator,       only : diag_ctrl, time_type, diag_update_remap_grids
use mom_diag_mediator,       only : diag_ctrl, query_averaging_enabled, enable_averages, disable_averaging
use mom_diag_mediator,       only : diag_grid_storage, diag_grid_storage_init, diag_grid_storage_end
use mom_diag_mediator,       only : diag_copy_diag_to_storage, diag_copy_storage_to_diag
use mom_diag_mediator,       only : diag_save_grids, diag_restore_grids
use mom_diapyc_energy_req,   only : diapyc_energy_req_init, diapyc_energy_req_end
use mom_diapyc_energy_req,   only : diapyc_energy_req_calc, diapyc_energy_req_test, diapyc_energy_req_cs
use mom_cvmix_conv,          only : cvmix_conv_init, cvmix_conv_cs
use mom_cvmix_conv,          only : cvmix_conv_end, calculate_cvmix_conv
use mom_domains,             only : pass_var, to_west, to_south, to_all, omit_corners
use mom_domains,             only : create_group_pass, do_group_pass, group_pass_type
use mom_energetic_pbl,       only : energetic_pbl, energetic_pbl_init
use mom_energetic_pbl,       only : energetic_pbl_end, energetic_pbl_cs
use mom_energetic_pbl,       only : energetic_pbl_get_mld
use mom_entrain_diffusive,   only : entrainment_diffusive, entrain_diffusive_init
use mom_entrain_diffusive,   only : entrain_diffusive_end, entrain_diffusive_cs
use mom_eos,                 only : calculate_density, calculate_tfreeze, eos_domain
use mom_error_handler,       only : mom_error, fatal, warning, calltree_showquery,mom_mesg
use mom_error_handler,       only : calltree_enter, calltree_leave, calltree_waypoint
use mom_file_parser,         only : get_param, log_version, param_file_type, read_param
use mom_forcing_type,        only : forcing, mom_forcing_chksum
use mom_forcing_type,        only : calculatebuoyancyflux2d, forcing_singlepointprint
use mom_geothermal,          only : geothermal_entraining, geothermal_in_place
use mom_geothermal,          only : geothermal_init, geothermal_end, geothermal_cs
use mom_grid,                only : ocean_grid_type
use mom_int_tide_input,      only : set_int_tide_input, int_tide_input_init
use mom_int_tide_input,      only : int_tide_input_end, int_tide_input_cs, int_tide_input_type
use mom_interface_heights,   only : find_eta
use mom_internal_tides,      only : propagate_int_tide
use mom_internal_tides,      only : internal_tides_init, internal_tides_end, int_tide_cs
use mom_kappa_shear,         only : kappa_shear_is_used
use mom_cvmix_kpp,           only : kpp_cs, kpp_init, kpp_compute_bld, kpp_calculate
use mom_cvmix_kpp,           only : kpp_end, kpp_get_bld
use mom_cvmix_kpp,           only : kpp_nonlocaltransport_temp, kpp_nonlocaltransport_saln
use mom_opacity,             only : opacity_init, opacity_end, opacity_cs
use mom_opacity,             only : absorbremainingsw, optics_type, optics_nbands
use mom_open_boundary,       only : ocean_obc_type
use mom_regularize_layers,   only : regularize_layers, regularize_layers_init, regularize_layers_cs
use mom_set_diffusivity,     only : set_diffusivity, set_bbl_tke
use mom_set_diffusivity,     only : set_diffusivity_init, set_diffusivity_end
use mom_set_diffusivity,     only : set_diffusivity_cs
use mom_sponge,              only : apply_sponge, sponge_cs
use mom_ale_sponge,          only : apply_ale_sponge, ale_sponge_cs
use mom_time_manager,        only : time_type, real_to_time, operator(-), operator(<=)
use mom_tracer_flow_control, only : call_tracer_column_fns, tracer_flow_control_cs
use mom_tracer_diabatic,     only : tracer_vertdiff
use mom_unit_scaling,        only : unit_scale_type
use mom_variables,           only : thermo_var_ptrs, vertvisc_type, accel_diag_ptrs
use mom_variables,           only : cont_diag_ptrs, mom_thermovar_chksum, p3d
use mom_verticalgrid,        only : verticalgrid_type, get_thickness_units
use mom_wave_speed,          only : wave_speeds
use mom_wave_interface,      only : wave_parameters_cs


implicit none ; private

#include <MOM_memory.h>

public diabatic
public diabatic_driver_init
public diabatic_driver_end
public extract_diabatic_member
public adiabatic
public adiabatic_driver_init
! public legacy_diabatic

! 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 this module
type, public:: diabatic_cs; private

  logical :: use_legacy_diabatic     !< If true (default), use the a legacy version of the diabatic
                                     !! algorithm. This is temporary and is needed to avoid change
                                     !! in answers.
  logical :: bulkmixedlayer          !< If true, a refined bulk mixed layer is used with
                                     !! nkml sublayers (and additional buffer layers).
  logical :: use_energetic_pbl       !< If true, use the implicit energetics planetary
                                     !! boundary layer scheme to determine the diffusivity
                                     !! in the surface boundary layer.
  logical :: use_kpp                 !< If true, use CVMix/KPP boundary layer scheme to determine the
                                     !! OBLD and the diffusivities within this layer.
  logical :: use_kappa_shear         !< If true, use the kappa_shear module to find the
                                     !! shear-driven diapycnal diffusivity.
  logical :: use_cvmix_shear         !< If true, use the CVMix module to find the
                                     !! shear-driven diapycnal diffusivity.
  logical :: use_cvmix_ddiff         !< If true, use the CVMix double diffusion module.
  logical :: use_cvmix_conv          !< If true, use the CVMix module to get enhanced
                                     !! mixing due to convection.
  logical :: double_diffuse          !< If true, some form of double-diffusive mixing is used.
  logical :: use_sponge              !< If true, sponges may be applied anywhere in the
                                     !! domain.  The exact location and properties of
                                     !! those sponges are set by calls to
                                     !! initialize_sponge and set_up_sponge_field.
  logical :: use_geothermal          !< If true, apply geothermal heating.
  logical :: use_int_tides           !< If true, use the code that advances a separate set
                                     !! of equations for the internal tide energy density.
  logical :: epbl_is_additive        !< If true, the diffusivity from ePBL is added to all
                                     !! other diffusivities. Otherwise, the larger of kappa-
                                     !! shear and ePBL diffusivities are used.
  real    :: epbl_prandtl            !< The Prandtl number used by ePBL to convert vertical
                                     !! diffusivities into viscosities.
  integer :: nmode = 1               !< Number of baroclinic modes to consider
  real    :: uniform_test_cg         !< Uniform group velocity of internal tide
                                     !! for testing internal tides [L T-1 ~> m s-1]
  logical :: usealealgorithm         !< If true, use the ALE algorithm rather than layered
                                     !! isopycnal/stacked shallow water mode. This logical
                                     !! passed by argument to diabatic_driver_init.
  logical :: aggregate_fw_forcing    !< Determines whether net incoming/outgoing surface
                                     !! FW fluxes are applied separately or combined before
                                     !! being applied.
  real :: ml_mix_first               !< The nondimensional fraction of the mixed layer
                                     !! algorithm that is applied before diffusive mixing.
                                     !! The default is 0, while 0.5 gives Strang splitting
                                     !! and 1 is a sensible value too.  Note that if there
                                     !! are convective instabilities in the initial state,
                                     !! the first call may do much more than the second.
  integer :: nkbl                    !< The number of buffer layers (if bulk_mixed_layer)
  logical :: massless_match_targets  !< If true (the default), keep the T & S
                                     !! consistent with the target values.
  logical :: mix_boundary_tracers    !< If true, mix the passive tracers in massless
                                     !! layers at the bottom into the interior as though
                                     !! a diffusivity of Kd_min_tr (see below) were
                                     !! operating.
  real    :: kd_bbl_tr               !< A bottom boundary layer tracer diffusivity that
                                     !! will allow for explicitly specified bottom fluxes
                                     !! [Z2 T-1 ~> m2 s-1].  The entrainment at the bottom is at
                                     !! least sqrt(Kd_BBL_tr*dt) over the same distance.
  real    :: kd_min_tr               !< A minimal diffusivity that should always be
                                     !! applied to tracers, especially in massless layers
                                     !! near the bottom [Z2 T-1 ~> m2 s-1].
  real    :: minimum_forcing_depth   !< The smallest depth over which heat and freshwater
                                     !! fluxes are applied [H ~> m or kg m-2].
  real    :: evap_cfl_limit = 0.8    !< The largest fraction of a layer that can be
                                     !! evaporated in one time-step [nondim].
  integer :: halo_ts_diff = 0        !< The temperature, salinity and thickness halo size that
                                     !! must be valid for the diffusivity calculations.
  logical :: usekpp = .false.        !< use CVMix/KPP diffusivities and non-local transport
  logical :: kppispassive            !< If true, KPP is in passive mode, not changing answers.
  logical :: debug                   !< If true, write verbose checksums for debugging purposes.
  logical :: debugconservation       !< If true, monitor conservation and extrema.
  logical :: tracer_tridiag          !< If true, use tracer_vertdiff instead of tridiagTS for
                                     !< vertical diffusion of T and S
  logical :: debug_energy_req        !< If true, test the mixing energy requirement code.
  type(diag_ctrl), pointer :: diag   !< structure used to regulate timing of diagnostic output
  real :: mlddensitydifference       !< Density difference used to determine MLD_user [R ~> kg m-3]
  real :: dz_subml_n2                !< The distance over which to calculate a diagnostic of the
                                     !! average stratification at the base of the mixed layer [Z ~> m].
  real :: mld_en_vals(3)             !< Energy values for energy mixed layer diagnostics

  !>@{ Diagnostic IDs
  integer :: id_cg1      = -1                 ! diag handle for mode-1 speed
  integer, allocatable, dimension(:) :: id_cn ! diag handle for all mode speeds
  integer :: id_wd       = -1, id_ea       = -1, id_eb           = -1 ! used by layer diabatic
  integer :: id_dudt_dia = -1, id_dvdt_dia = -1, id_ea_s         = -1, id_eb_s     = -1
  integer :: id_ea_t     = -1, id_eb_t     = -1
  integer :: id_kd_heat  = -1, id_kd_salt  = -1, id_kd_interface = -1, id_kd_epbl  = -1
  integer :: id_tdif     = -1, id_tadv     = -1, id_sdif         = -1, id_sadv     = -1
  integer :: id_mld_003  = -1, id_mld_0125  = -1, id_mld_user     = -1, id_mlotstsq = -1
  integer :: id_mld_en1 = -1, id_mld_en2 = -1, id_mld_en3= -1
  integer :: id_submln2  = -1
  integer :: id_hf_dudt_dia_2d = -1, id_hf_dvdt_dia_2d = -1

  ! diagnostic for fields prior to applying diapycnal physics
  integer :: id_u_predia = -1, id_v_predia = -1, id_h_predia = -1
  integer :: id_t_predia = -1, id_s_predia = -1, id_e_predia = -1

  integer :: id_diabatic_diff_temp_tend     = -1
  integer :: id_diabatic_diff_saln_tend     = -1
  integer :: id_diabatic_diff_heat_tend     = -1
  integer :: id_diabatic_diff_salt_tend     = -1
  integer :: id_diabatic_diff_heat_tend_2d  = -1
  integer :: id_diabatic_diff_salt_tend_2d  = -1
  integer :: id_diabatic_diff_h= -1

  integer :: id_boundary_forcing_h       = -1
  integer :: id_boundary_forcing_h_tendency   = -1
  integer :: id_boundary_forcing_temp_tend    = -1
  integer :: id_boundary_forcing_saln_tend    = -1
  integer :: id_boundary_forcing_heat_tend    = -1
  integer :: id_boundary_forcing_salt_tend    = -1
  integer :: id_boundary_forcing_heat_tend_2d = -1
  integer :: id_boundary_forcing_salt_tend_2d = -1

  integer :: id_frazil_h    = -1
  integer :: id_frazil_temp_tend    = -1
  integer :: id_frazil_heat_tend    = -1
  integer :: id_frazil_heat_tend_2d = -1
  !>@}

  logical :: diabatic_diff_tendency_diag = .false. !< If true calculate diffusive tendency diagnostics
  logical :: boundary_forcing_tendency_diag = .false. !< If true calculate frazil diagnostics
  logical :: frazil_tendency_diag = .false. !< If true calculate frazil tendency diagnostics
  real, allocatable, dimension(:,:,:) :: frazil_heat_diag !< diagnose 3d heat tendency from frazil
  real, allocatable, dimension(:,:,:) :: frazil_temp_diag !< diagnose 3d temp tendency from frazil

  type(diabatic_aux_cs),        pointer :: diabatic_aux_csp      => null() !< Control structure for a child module
  type(entrain_diffusive_cs),   pointer :: entrain_diffusive_csp => null() !< Control structure for a child module
  type(bulkmixedlayer_cs),      pointer :: bulkmixedlayer_csp    => null() !< Control structure for a child module
  type(energetic_pbl_cs),       pointer :: energetic_pbl_csp     => null() !< Control structure for a child module
  type(regularize_layers_cs),   pointer :: regularize_layers_csp => null() !< Control structure for a child module
  type(geothermal_cs),          pointer :: geothermal_csp        => null() !< Control structure for a child module
  type(int_tide_cs),            pointer :: int_tide_csp          => null() !< Control structure for a child module
  type(int_tide_input_cs),      pointer :: int_tide_input_csp    => null() !< Control structure for a child module
  type(int_tide_input_type),    pointer :: int_tide_input        => null() !< Control structure for a child module
  type(opacity_cs),             pointer :: opacity_csp           => null() !< Control structure for a child module
  type(set_diffusivity_cs),     pointer :: set_diff_csp          => null() !< Control structure for a child module
  type(sponge_cs),              pointer :: sponge_csp            => null() !< Control structure for a child module
  type(ale_sponge_cs),          pointer :: ale_sponge_csp        => null() !< Control structure for a child module
  type(tracer_flow_control_cs), pointer :: tracer_flow_csp       => null() !< Control structure for a child module
  type(optics_type),            pointer :: optics                => null() !< Control structure for a child module
  type(kpp_cs),                 pointer :: kpp_csp               => null() !< Control structure for a child module
  type(cvmix_conv_cs),          pointer :: cvmix_conv_csp        => null() !< Control structure for a child module
  type(diapyc_energy_req_cs),   pointer :: diapyc_en_rec_csp     => null() !< Control structure for a child module

  type(group_pass_type) :: pass_hold_eb_ea !< For group halo pass
  type(group_pass_type) :: pass_kv         !< For group halo pass
  type(diag_grid_storage) :: diag_grids_prev!< Stores diagnostic grids at some previous point in the algorithm
  ! Data arrays for communicating between components
  real, allocatable, dimension(:,:,:) :: kpp_nltheat    !< KPP non-local transport for heat [m s-1]
  real, allocatable, dimension(:,:,:) :: kpp_nltscalar  !< KPP non-local transport for scalars [m s-1]
  real, allocatable, dimension(:,:,:) :: kpp_buoy_flux  !< KPP forcing buoyancy flux [L2 T-3 ~> m2 s-3]
  real, allocatable, dimension(:,:)   :: kpp_temp_flux  !< KPP effective temperature flux [degC m s-1]
  real, allocatable, dimension(:,:)   :: kpp_salt_flux  !< KPP effective salt flux [ppt m s-1]

  type(time_type), pointer :: time !< Pointer to model time (needed for sponges)
end type diabatic_cs

!>@{ clock ids
integer :: id_clock_entrain, id_clock_mixedlayer, id_clock_set_diffusivity
integer :: id_clock_tracers, id_clock_tridiag, id_clock_pass, id_clock_sponge
integer :: id_clock_geothermal, id_clock_differential_diff, id_clock_remap
integer :: id_clock_kpp
!>@}

contains

!>  This subroutine imposes the diapycnal mass fluxes and the
!!  accompanying diapycnal advection of momentum and tracers.
subroutine diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
                    G, GV, US, CS, OBC, WAVES)
  type(ocean_grid_type),                     intent(inout) :: G         !< ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: GV        !< ocean vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u         !< zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v         !< meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h         !< thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),                     intent(inout) :: tv        !< points to thermodynamic fields
                                                                        !! unused have NULL ptrs
  real, dimension(:,:),                      pointer       :: Hml       !< Active mixed layer depth [Z ~> m]
  type(forcing),                             intent(inout) :: fluxes    !< points to forcing fields
                                                                        !! unused fields have NULL ptrs
  type(vertvisc_type),                       intent(inout) :: visc      !< vertical viscosities, BBL properies, and
  type(accel_diag_ptrs),                     intent(inout) :: ADp       !< related points to accelerations in momentum
                                                                        !! equations, to enable the later derived
                                                                        !! diagnostics, like energy budgets
  type(cont_diag_ptrs),                      intent(inout) :: CDp       !< points to terms in continuity equations
  real,                                      intent(in)    :: dt        !< time increment [T ~> s]
  type(time_type),                           intent(in)    :: Time_end  !< Time at the end of the interval
  type(unit_scale_type),                     intent(in)    :: US        !< A dimensional unit scaling type
  type(diabatic_cs),                         pointer       :: CS        !< module control structure
  type(ocean_obc_type),            optional, pointer       :: OBC       !< Open boundaries control structure.
  type(wave_parameters_cs),        optional, pointer       :: Waves     !< Surface gravity waves

  ! local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    eta      ! Interface heights before diapycnal mixing [m].
  real, dimension(SZI_(G),SZJ_(G),CS%nMode) :: &
    cn_IGW   ! baroclinic internal gravity wave speeds
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag             ! diagnostic array for temp
  integer :: i, j, k, m, is, ie, js, je, nz
  logical :: showCallTree ! If true, show the call tree

  if (g%ke == 1) return

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

  if (.not. associated(cs)) call mom_error(fatal, "MOM_diabatic_driver: "// &
         "Module must be initialized before it is used.")
  if (dt == 0.0) call mom_error(fatal, "MOM_diabatic_driver: "// &
        "diabatic was called with a zero length timestep.")
  if (dt < 0.0) call mom_error(fatal, "MOM_diabatic_driver: "// &
        "diabatic was called with a negative timestep.")

  showcalltree = calltree_showquery()

  ! Offer diagnostics of various state varables at the start of diabatic
  ! these are mostly for debugging purposes.
  if (cs%id_u_predia > 0) call post_data(cs%id_u_predia, u, cs%diag)
  if (cs%id_v_predia > 0) call post_data(cs%id_v_predia, v, cs%diag)
  if (cs%id_h_predia > 0) call post_data(cs%id_h_predia, h, cs%diag)
  if (cs%id_T_predia > 0) call post_data(cs%id_T_predia, tv%T, cs%diag)
  if (cs%id_S_predia > 0) call post_data(cs%id_S_predia, tv%S, cs%diag)
  if (cs%id_e_predia > 0) then
    call find_eta(h, tv, g, gv, us, eta, eta_to_m=1.0)
    call post_data(cs%id_e_predia, eta, cs%diag)
  endif

  if (cs%debug) then
    call mom_state_chksum("Start of diabatic ", u, v, h, g, gv, us, haloshift=0)
    call mom_forcing_chksum("Start of diabatic", fluxes, g, us, haloshift=0)
  endif
  if (cs%debugConservation) call mom_state_stats('Start of diabatic', u, v, h, tv%T, tv%S, g, gv, us)

  if (cs%debug_energy_req) &
    call diapyc_energy_req_test(h, dt, tv, g, gv, us, cs%diapyc_en_rec_CSp)

  call cpu_clock_begin(id_clock_set_diffusivity)
  call set_bbl_tke(u, v, h, fluxes, visc, g, gv, us, cs%set_diff_CSp, obc=obc)
  call cpu_clock_end(id_clock_set_diffusivity)

  ! Frazil formation keeps the temperature above the freezing point.
  ! make_frazil is deliberately called at both the beginning and at
  ! the end of the diabatic processes.
  if (associated(tv%T) .AND. associated(tv%frazil)) then
    ! For frazil diagnostic, the first call covers the first half of the time step
    call enable_averages(0.5*dt, time_end - real_to_time(0.5*us%T_to_s*dt), cs%diag)
    if (cs%frazil_tendency_diag) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        temp_diag(i,j,k) = tv%T(i,j,k)
      enddo ; enddo ; enddo
    endif

    if (associated(fluxes%p_surf_full)) then
      call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp, fluxes%p_surf_full, halo=cs%halo_TS_diff)
    else
      call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp, halo=cs%halo_TS_diff)
    endif
    if (showcalltree) call calltree_waypoint("done with 1st make_frazil (diabatic)")

    if (cs%frazil_tendency_diag) then
      call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, g, gv, us, cs)
      if (cs%id_frazil_h > 0) call post_data(cs%id_frazil_h, h, cs%diag)
    endif
    call disable_averaging(cs%diag)
  endif ! associated(tv%T) .AND. associated(tv%frazil)
  if (cs%debugConservation) call mom_state_stats('1st make_frazil', u, v, h, tv%T, tv%S, g, gv, us)

  if (cs%use_int_tides) then
    ! This block provides an interface for the unresolved low-mode internal tide module.
    call set_int_tide_input(u, v, h, tv, fluxes, cs%int_tide_input, dt, g, gv, us, &
                            cs%int_tide_input_CSp)
    cn_igw(:,:,:) = 0.0
    if (cs%uniform_test_cg > 0.0) then
      do m=1,cs%nMode ; cn_igw(:,:,m) = cs%uniform_test_cg ; enddo
    else
      call wave_speeds(h, tv, g, gv, us, cs%nMode, cn_igw, full_halos=.true.)
    endif

    call propagate_int_tide(h, tv, cn_igw, cs%int_tide_input%TKE_itidal_input, cs%int_tide_input%tideamp, &
                            cs%int_tide_input%Nb, dt, g, gv, us, cs%int_tide_CSp)
    if (showcalltree) call calltree_waypoint("done with propagate_int_tide (diabatic)")
  endif ! end CS%use_int_tides

  if (cs%useALEalgorithm .and. cs%use_legacy_diabatic) then
    call diabatic_ale_legacy(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
                      g, gv, us, cs, waves)
  elseif (cs%useALEalgorithm) then
    call diabatic_ale(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
                      g, gv, us, cs, waves)
  else
    call layered_diabatic(u, v, h, tv, hml, fluxes, visc, adp, cdp, dt, time_end, &
                          g, gv, us, cs, waves)
  endif


  call cpu_clock_begin(id_clock_pass)
  if (associated(visc%Kv_shear)) &
    call pass_var(visc%Kv_shear, g%Domain, to_all+omit_corners, halo=1)
  call cpu_clock_end(id_clock_pass)

  call disable_averaging(cs%diag)
  ! Frazil formation keeps temperature above the freezing point.
  ! make_frazil is deliberately called at both the beginning and at
  ! the end of the diabatic processes.
  if (associated(tv%T) .AND. associated(tv%frazil)) then
    call enable_averages(0.5*dt, time_end, cs%diag)
    if (cs%frazil_tendency_diag) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        temp_diag(i,j,k) = tv%T(i,j,k)
      enddo ; enddo ; enddo
    endif

    if (associated(fluxes%p_surf_full)) then
      call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp, fluxes%p_surf_full)
    else
      call make_frazil(h, tv, g, gv, us, cs%diabatic_aux_CSp)
    endif

    if (cs%frazil_tendency_diag) then
      call diagnose_frazil_tendency(tv, h, temp_diag, 0.5*dt, g, gv, us, cs)
      if (cs%id_frazil_h > 0 ) call post_data(cs%id_frazil_h, h, cs%diag)
    endif

    if (showcalltree) call calltree_waypoint("done with 2nd make_frazil (diabatic)")
    if (cs%debugConservation) call mom_state_stats('2nd make_frazil', u, v, h, tv%T, tv%S, g, gv, us)
    call disable_averaging(cs%diag)

  endif  ! endif for frazil


  ! Diagnose mixed layer depths.
  call enable_averages(dt, time_end, cs%diag)
  if (cs%id_MLD_003 > 0 .or. cs%id_subMLN2 > 0 .or. cs%id_mlotstsq > 0) then
    call diagnosemldbydensitydifference(cs%id_MLD_003, h, tv, 0.03*us%kg_m3_to_R, g, gv, us, cs%diag, &
                                        id_n2subml=cs%id_subMLN2, id_mldsq=cs%id_mlotstsq, dz_subml=cs%dz_subML_N2)
  endif
  if (cs%id_MLD_0125 > 0) then
    call diagnosemldbydensitydifference(cs%id_MLD_0125, h, tv, 0.125*us%kg_m3_to_R, g, gv, us, cs%diag)
  endif
  if (cs%id_MLD_user > 0) then
    call diagnosemldbydensitydifference(cs%id_MLD_user, h, tv, cs%MLDdensityDifference, g, gv, us, cs%diag)
  endif
  if ((cs%id_MLD_EN1 > 0) .or. (cs%id_MLD_EN2 > 0) .or. (cs%id_MLD_EN3 > 0)) then
    call diagnosemldbyenergy((/cs%id_MLD_EN1, cs%id_MLD_EN2, cs%id_MLD_EN3/),&
         h, tv, g, gv, us, cs%MLD_EN_VALS, cs%diag)
  endif
  if (cs%use_int_tides) then
    if (cs%id_cg1 > 0) call post_data(cs%id_cg1, cn_igw(:,:,1),cs%diag)
    do m=1,cs%nMode ; if (cs%id_cn(m) > 0) call post_data(cs%id_cn(m), cn_igw(:,:,m), cs%diag) ; enddo
  endif
  call disable_averaging(cs%diag)

  if (cs%debugConservation) call mom_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, g, gv, us)

end subroutine diabatic


!> Applies diabatic forcing and diapycnal mixing of temperature, salinity and other tracers for use
!! with an ALE algorithm.  This version uses an older set of algorithms compared with diabatic_ALE.
subroutine diabatic_ale_legacy(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
                           G, GV, US, CS, Waves)
  type(ocean_grid_type),                     intent(inout) :: G         !< ocean 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, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u         !< zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v         !< meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h         !< thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),                     intent(inout) :: tv        !< points to thermodynamic fields
                                                                        !! unused have NULL ptrs
  real, dimension(:,:),                      pointer       :: Hml       !< Active mixed layer depth [Z ~> m]
  type(forcing),                             intent(inout) :: fluxes    !< points to forcing fields
                                                                        !! unused fields have NULL ptrs
  type(vertvisc_type),                       intent(inout) :: visc      !< vertical viscosities, BBL properies, and
  type(accel_diag_ptrs),                     intent(inout) :: ADp       !< related points to accelerations in momentum
                                                                        !! equations, to enable the later derived
                                                                        !! diagnostics, like energy budgets
  type(cont_diag_ptrs),                      intent(inout) :: CDp       !< points to terms in continuity equations
  real,                                      intent(in)    :: dt        !< time increment [T ~> s]
  type(time_type),                           intent(in)    :: Time_end  !< Time at the end of the interval
  type(diabatic_cs),                         pointer       :: CS        !< module control structure
  type(wave_parameters_cs),        optional, pointer       :: Waves     !< Surface gravity waves

  ! local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    ea_s,     &  ! amount of fluid entrained from the layer above within
                 ! one time step  [H ~> m or kg m-2]
    eb_s,     &  ! amount of fluid entrained from the layer below within
                 ! one time step  [H ~> m or kg m-2]
    ea_t,     &  ! amount of fluid entrained from the layer above within
                 ! one time step  [H ~> m or kg m-2]
    eb_t,     &  ! amount of fluid entrained from the layer below within
                 ! one time step  [H ~> m or kg m-2]
    h_orig, &    ! initial layer thicknesses [H ~> m or kg m-2]
    h_prebound, & ! initial layer thicknesses [H ~> m or kg m-2]
    hold,   &    ! layer thickness before diapycnal entrainment, and later
                 ! the initial layer thicknesses (if a mixed layer is used),
                 ! [H ~> m or kg m-2]
    dsv_dt, &    ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
    dsv_ds, &    ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
    ctke,   &    ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
    u_h,    &    ! zonal and meridional velocities at thickness points after
    v_h          ! entrainment [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G)) :: &
    SkinBuoyFlux! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag                ! diagnostic array for thickness
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag             ! diagnostic array for temp
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag             ! diagnostic array for salinity
  real, dimension(SZI_(G),SZJ_(G))      :: tendency_2d           ! depth integrated content tendency for diagn

  real :: net_ent  ! The net of ea-eb at an interface.

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    eatr, &  ! The equivalent of ea and eb for tracers, which differ from ea and
    ebtr     ! eb in that they tend to homogenize tracers in massless layers
             ! near the boundaries [H ~> m or kg m-2] (for Bous or non-Bouss)

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    Kd_int,   & ! diapycnal diffusivity of interfaces [Z2 T-1 ~> m2 s-1]
    Kd_heat,  & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1]
    Kd_salt,  & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1]
    Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to
                ! Kd_int [Z2 T-1 ~> m2 s-1].
    kd_extra_s , & !  The extra diffusivity of salinity due to double diffusion relative to
                ! Kd_int [Z2 T-1 ~> m2 s-1].
    kd_epbl,  & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]
    tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
    tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
    sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
    sadv_flx    ! advective diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

  real, allocatable, dimension(:,:) :: &
    hf_dudt_dia_2d, hf_dvdt_dia_2d ! Depth sum of diapycnal mixing accelaration * fract. thickness [L T-2 ~> m s-2].

  logical :: in_boundary(szi_(g)) ! True if there are no massive layers below,
                                  ! where massive is defined as sufficiently thick that
                                  ! the no-flux boundary conditions have not restricted
                                  ! the entrainment - usually sqrt(Kd*dt).

  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 :: h_neglect2   ! h_neglect^2 [H2 ~> m2 or kg2 m-4]
  real :: add_ent      ! Entrainment that needs to be added when mixing tracers
                       ! [H ~> m or kg m-2]
  real :: eaval        ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
  real :: hval         ! hval is 2*h at velocity grid points [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 :: Tr_ea_BBL    ! The diffusive tracer thickness in the BBL that is
                       ! coupled to the bottom within a timestep [H ~> m or kg m-2]

  real :: htot(szib_(g))             ! The summed thickness from the bottom [H ~> m or kg m-2].
  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 :: Ent_int ! The diffusive entrainment rate at an interface [H ~> m or kg m-2]
  real :: Idt     ! The inverse time step [T-1 ~> s-1]

  integer :: dir_flag     ! An integer encoding the directions in which to do halo updates.
  logical :: showCallTree ! If true, show the call tree
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m

  integer :: ig, jg      ! global indices for testing testing itide point source (BDM)
  real :: Kd_add_here    ! An added diffusivity [Z2 T-1 ~> m2 s-1].

  is   = g%isc  ; ie  = g%iec  ; js  = g%jsc  ; je  = g%jec ; nz = g%ke
  isq  = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
  kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0

  showcalltree = calltree_showquery()
  if (showcalltree) call calltree_enter("diabatic_ALE_legacy(), MOM_diabatic_driver.F90")

  ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
  call enable_averages(dt, time_end, cs%diag)

  if (cs%use_geothermal) then
    call cpu_clock_begin(id_clock_geothermal)
    call geothermal_in_place(h, tv, dt, g, gv, us, cs%geothermal_CSp, halo=cs%halo_TS_diff)
    call cpu_clock_end(id_clock_geothermal)
    if (showcalltree) call calltree_waypoint("geothermal (diabatic)")
    if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)
  endif

  ! Whenever thickness changes let the diag manager know, target grids
  ! for vertical remapping may need to be regenerated.
  call diag_update_remap_grids(cs%diag)

  ! Set_pen_shortwave estimates the optical properties of the water column.
  ! It will need to be modified later to include information about the
  ! biological properties and layer thicknesses.
  if (associated(cs%optics)) &
    call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)

  if (cs%debug) call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)

  if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then
    if (cs%use_geothermal) then
      ! The presence of eatr and ebtr causes find_uv_at_h to use a tridiagonal solver,
      ! which changes answers at the level of roundoff because ((A*B / A) /= B).
      eatr(:,:,:) = 0.0 ; ebtr(:,:,:) = 0.0
      call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us, eatr, ebtr)
    else
      call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
    endif
    if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
  endif

  call cpu_clock_begin(id_clock_set_diffusivity)
  ! Sets: Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb
  ! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow
  if (cs%debug) &
    call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)
  if (cs%double_diffuse) then
    call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, cs%set_diff_CSp, &
                         kd_int=kd_int, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)
  else
    call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, &
                         cs%set_diff_CSp, kd_int=kd_int)
  endif
  call cpu_clock_end(id_clock_set_diffusivity)
  if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")

  ! Set diffusivities for heat and salt separately

  if (cs%useKPP) then
    ! Add contribution from double diffusion
    if (cs%double_diffuse) then
      !$OMP parallel do default(shared)
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_salt(i,j,k) = kd_int(i,j,k) + kd_extra_s(i,j,k)
        kd_heat(i,j,k) = kd_int(i,j,k) + kd_extra_t(i,j,k)
      enddo ; enddo ; enddo
    else
      !$OMP parallel do default(shared)
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_salt(i,j,k) = kd_int(i,j,k)
        kd_heat(i,j,k) = kd_int(i,j,k)
      enddo ; enddo ; enddo
    endif
  endif

  if (cs%debug) then
    call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)
    call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)
    call mom_thermovar_chksum("after set_diffusivity ", tv, g)
    call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
    call hchksum(kd_heat, "after set_diffusivity Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
    call hchksum(kd_salt, "after set_diffusivity Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
  endif

  if (cs%useKPP) then
    call cpu_clock_begin(id_clock_kpp)
    ! total vertical viscosity in the interior is represented via visc%Kv_shear

    ! KPP needs the surface buoyancy flux but does not update state variables.
    ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
    ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux
    ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second)
    ! unlike other instances where the fluxes are integrated in time over a time-step.
    call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
                                 cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
    ! The KPP scheme calculates boundary layer diffusivities and non-local transport.

    call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &
                         fluxes%ustar, cs%KPP_buoy_flux, waves=waves)

    call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
                       kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)

    if (associated(hml)) then
      call kpp_get_bld(cs%KPP_CSp, hml(:,:), g, us)
      call pass_var(hml, g%domain, halo=1)
      ! If visc%MLD exists, copy KPP's BLD into it
      if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
    endif

    if (.not.cs%KPPisPassive) then
      !$OMP parallel do default(shared)
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_int(i,j,k) = min( kd_salt(i,j,k),  kd_heat(i,j,k) )
      enddo ; enddo ; enddo
      if (cs%double_diffuse) then
        !$OMP parallel do default(shared)
        do k=1,nz+1 ; do j=js,je ; do i=is,ie
          kd_extra_s(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))
          kd_extra_t(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))
        enddo ; enddo ; enddo
      endif
    endif ! not passive

    call cpu_clock_end(id_clock_kpp)
    if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")
    if (cs%debug) then
      call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)
      call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)
      call mom_thermovar_chksum("after KPP", tv, g)
      call hchksum(kd_heat, "after KPP Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
      call hchksum(kd_salt, "after KPP Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
    endif

  endif  ! endif for KPP


  if (cs%useKPP) then
    call cpu_clock_begin(id_clock_kpp)
    if (cs%debug) then
      call hchksum(cs%KPP_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(cs%KPP_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(cs%KPP_NLTheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)
      call hchksum(cs%KPP_NLTscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
    endif
    ! Apply non-local transport of heat and salt
    ! Changes: tv%T, tv%S
    call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat,   cs%KPP_temp_flux, &
                                    us%T_to_s*dt, tv%T, us%Q_to_J_kg*tv%C_p)
    call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, &
                                    us%T_to_s*dt, tv%S)
    call cpu_clock_end(id_clock_kpp)
    if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
    if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)

    if (cs%debug) then
      call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)
      call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)
      call mom_thermovar_chksum("after KPP_applyNLT ", tv, g)
    endif
  endif ! endif for KPP

  ! This is the "old" method for applying differential diffusion.
  ! Changes: tv%T, tv%S
  if (cs%double_diffuse .and. associated(tv%T)) then

    call cpu_clock_begin(id_clock_differential_diff)
    call differential_diffuse_t_s(h, tv%T, tv%S, kd_extra_t, kd_extra_s, dt, g, gv)
    call cpu_clock_end(id_clock_differential_diff)

    if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")
    if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g, gv, us)

    ! increment heat and salt diffusivity.
    ! CS%useKPP==.true. already has extra_T and extra_S included
    if (.not. cs%useKPP) then
      !$OMP parallel do default(shared)
      do k=2,nz ; do j=js,je ; do i=is,ie
        kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)
        kd_salt(i,j,k) = kd_salt(i,j,k) + kd_extra_s(i,j,k)
      enddo ; enddo ; enddo
    endif

  endif

  ! Calculate vertical mixing due to convection (computed via CVMix)
  if (cs%use_CVMix_conv) then
    ! Increment vertical diffusion and viscosity due to convection
    call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml, kd=kd_int, kv=visc%Kv_slow)
  endif

  ! This block sets ea, eb from h and Kd_int.
  do j=js,je ; do i=is,ie
    ea_s(i,j,1) = 0.0
  enddo ; enddo
  !$OMP parallel do default(shared)  private(hval)
  do k=2,nz ; do j=js,je ; do i=is,ie
    hval=1.0/(h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
    ea_s(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_int(i,j,k)
    eb_s(i,j,k-1) = ea_s(i,j,k)
    ea_t(i,j,k-1) = ea_s(i,j,k-1) ; eb_t(i,j,k-1) = eb_s(i,j,k-1)
  enddo ; enddo ; enddo
  do j=js,je ; do i=is,ie
    eb_s(i,j,nz) = 0.0
    ea_t(i,j,nz) = ea_s(i,j,nz) ; eb_t(i,j,nz) = eb_s(i,j,nz)
  enddo ; enddo
  if (showcalltree) call calltree_waypoint("done setting ea,eb from Kd_int (diabatic)")

  if (cs%debug) then
    call mom_forcing_chksum("after calc_entrain ", fluxes, g, us, haloshift=0)
    call mom_thermovar_chksum("after calc_entrain ", tv, g)
    call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, us, haloshift=0)
    call hchksum(ea_s, "after calc_entrain ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
    call hchksum(eb_s, "after calc_entrain eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
  endif

  ! Save fields before boundary forcing is applied for tendency diagnostics
  if (cs%boundary_forcing_tendency_diag) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      h_diag(i,j,k)    = h(i,j,k)
      temp_diag(i,j,k) = tv%T(i,j,k)
      saln_diag(i,j,k) = tv%S(i,j,k)
    enddo ; enddo ; enddo
  endif

  ! Apply forcing
  call cpu_clock_begin(id_clock_remap)

  ! Changes made to following fields:  h, tv%T and tv%S.
  do k=1,nz ; do j=js,je ; do i=is,ie
    h_prebound(i,j,k) = h(i,j,k)
  enddo ; enddo ; enddo
  if (cs%use_energetic_PBL) then

    skinbuoyflux(:,:) = 0.0
    call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
            optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, cs%evap_CFL_limit,                         &
            cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux)

    if (cs%debug) then
      call hchksum(ea_t, "after applyBoundaryFluxes ea_t", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_t, "after applyBoundaryFluxes eb_t", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(ea_s, "after applyBoundaryFluxes ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_s, "after applyBoundaryFluxes eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(ctke, "after applyBoundaryFluxes cTKE", g%HI, haloshift=0, &
                   scale=us%RZ3_T3_to_W_m2*us%T_to_s)
      call hchksum(dsv_dt, "after applyBoundaryFluxes dSV_dT", g%HI, haloshift=0, scale=us%kg_m3_to_R)
      call hchksum(dsv_ds, "after applyBoundaryFluxes dSV_dS", g%HI, haloshift=0, scale=us%kg_m3_to_R)
    endif

    call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
    call energetic_pbl(h, u_h, v_h, tv, fluxes, dt, kd_epbl, g, gv, us, &
                       cs%energetic_PBL_CSp, dsv_dt, dsv_ds, ctke, skinbuoyflux, waves=waves)

    if (associated(hml)) then
      call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hml(:,:), g, us)
      call pass_var(hml, g%domain, halo=1)
      ! If visc%MLD exists, copy ePBL's MLD into it
      if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
    elseif (associated(visc%MLD)) then
      call energetic_pbl_get_mld(cs%energetic_PBL_CSp, visc%MLD, g, us)
      call pass_var(visc%MLD, g%domain, halo=1)
    endif

    ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
    do k=2,nz ; do j=js,je ; do i=is,ie
      if (cs%ePBL_is_additive) then
        kd_add_here = kd_epbl(i,j,k)
        visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%ePBL_Prandtl*kd_epbl(i,j,k)
      else
        kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)
        visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), cs%ePBL_Prandtl*kd_epbl(i,j,k))
      endif

      ent_int = kd_add_here * (gv%Z_to_H**2 * dt) / &
                  (0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect)
      eb_s(i,j,k-1) = eb_s(i,j,k-1) + ent_int
      ea_s(i,j,k) = ea_s(i,j,k) + ent_int
      kd_int(i,j,k) = kd_int(i,j,k) + kd_add_here

      ! for diagnostics
      kd_heat(i,j,k) = kd_heat(i,j,k) + kd_int(i,j,k)
      kd_salt(i,j,k) = kd_salt(i,j,k) + kd_int(i,j,k)

    enddo ; enddo ; enddo

    if (cs%debug) then
      call hchksum(ea_t, "after ePBL ea_t", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_t, "after ePBL eb_t", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(ea_s, "after ePBL ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_s, "after ePBL eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(kd_epbl, "after ePBL Kd_ePBL", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
    endif

  else
    call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
                                  optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, &
                                  cs%evap_CFL_limit, cs%minimum_forcing_depth)

  endif   ! endif for CS%use_energetic_PBL

  ! diagnose the tendencies due to boundary forcing
  ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme
  !  so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards
  if (cs%boundary_forcing_tendency_diag) then
    call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, g, gv, us, cs)
    if (cs%id_boundary_forcing_h > 0) call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h = h_diag)
  endif
  ! Boundary fluxes may have changed T, S, and h
  call diag_update_remap_grids(cs%diag)
  call cpu_clock_end(id_clock_remap)
  if (cs%debug) then
    call mom_forcing_chksum("after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)
    call mom_thermovar_chksum("after applyBoundaryFluxes ", tv, g)
    call mom_state_chksum("after applyBoundaryFluxes ", u, v, h, g, gv, us, haloshift=0)
  endif
  if (showcalltree) call calltree_waypoint("done with applyBoundaryFluxes (diabatic)")
  if (cs%debugConservation)  call mom_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, g, gv, us)

  ! Update h according to divergence of the difference between
  ! ea and eb. We keep a record of the original h in hold.
  ! In the following, the checks for negative values are to guard against
  ! instances where entrainment drives a layer to negative thickness.
  !### This code may be unnecessary, but the negative-thickness checks do appear to change
  !    answers slightly in some cases.
  !$OMP parallel do default(shared)
  do j=js,je
    do i=is,ie
      hold(i,j,1) = h(i,j,1)
      ! Does nothing with ALE:  h(i,j,1) = h(i,j,1) + (eb_s(i,j,1) - ea_s(i,j,2))
      hold(i,j,nz) = h(i,j,nz)
      ! Does nothing with ALE:  h(i,j,nz) = h(i,j,nz) + (ea_s(i,j,nz) - eb_s(i,j,nz-1))
      if (h(i,j,1) <= 0.0) h(i,j,1) = gv%Angstrom_H
      if (h(i,j,nz) <= 0.0) h(i,j,nz) = gv%Angstrom_H
    enddo
    do k=2,nz-1 ; do i=is,ie
      hold(i,j,k) = h(i,j,k)
      ! Does nothing with ALE:  h(i,j,k) = h(i,j,k) + ((ea_s(i,j,k) - eb_s(i,j,k-1)) + &
      !                                                (eb_s(i,j,k) - ea_s(i,j,k+1)))
      if (h(i,j,k) <= 0.0) h(i,j,k) = gv%Angstrom_H
    enddo ; enddo
  enddo
  ! Checks for negative thickness may have changed layer thicknesses
  call diag_update_remap_grids(cs%diag)

  if (cs%debug) then
    call mom_state_chksum("after negative check ", u, v, h, g, gv, us, haloshift=0)
    call mom_forcing_chksum("after negative check ", fluxes, g, us, haloshift=0)
    call mom_thermovar_chksum("after negative check ", tv, g)
  endif
  if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")
  if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)

  ! calculate change in temperature & salinity due to dia-coordinate surface diffusion
  if (associated(tv%T)) then

    if (cs%debug) then
      call hchksum(ea_t, "before triDiagTS ea_t ", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_t, "before triDiagTS eb_t ", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(ea_s, "before triDiagTS ea_s ", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_s, "before triDiagTS eb_s ", g%HI, haloshift=0, scale=gv%H_to_m)
    endif

    call cpu_clock_begin(id_clock_tridiag)
    !  Keep salinity from falling below a small but positive threshold.
    !  This constraint is needed for SIS1 ice model, which can extract
    !  more salt than is present in the ocean. SIS2 does not suffer
    !  from this limitation, in which case we can let salinity=0 and still
    !  have salt conserved with SIS2 ice. So for SIS2, we can run with
    !  BOUND_SALINITY=False in MOM.F90.
    if (associated(tv%S) .and. associated(tv%salt_deficit)) &
      call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

    if (cs%diabatic_diff_tendency_diag) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        temp_diag(i,j,k) = tv%T(i,j,k)
        saln_diag(i,j,k) = tv%S(i,j,k)
      enddo ; enddo ; enddo
    endif

    ! Changes T and S via the tridiagonal solver; no change to h
    do k=1,nz ; do j=js,je ; do i=is,ie
      ea_t(i,j,k) = ea_s(i,j,k) ; eb_t(i,j,k) = eb_s(i,j,k)
    enddo ; enddo ; enddo
    if (cs%tracer_tridiag) then
      call tracer_vertdiff(hold, ea_t, eb_t, dt, tv%T, g, gv)
      call tracer_vertdiff(hold, ea_s, eb_s, dt, tv%S, g, gv)
    else
      call tridiagts(g, gv, is, ie, js, je, hold, ea_s, eb_s, tv%T, tv%S)
    endif

    ! diagnose temperature, salinity, heat, and salt tendencies
    ! Note: hold here refers to the thicknesses from before the dual-entraintment when using
    ! the bulk mixed layer scheme. Otherwise in ALE-mode, layer thicknesses will (not?) have changed
    ! In either case, tendencies should be posted on hold
    if (cs%diabatic_diff_tendency_diag) then
      call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, us, cs)
      if (cs%id_diabatic_diff_h > 0) call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h = hold)
    endif

    call cpu_clock_end(id_clock_tridiag)

    if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")

  endif  ! endif corresponding to if (associated(tv%T))

  if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)

  if (cs%debug) then
    call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)
    call mom_thermovar_chksum("after mixed layer ", tv, g)
  endif

  ! Whenever thickness changes let the diag manager know, as the
  ! target grids for vertical remapping may need to be regenerated.
  if (associated(adp%du_dt_dia) .or. associated(adp%dv_dt_dia)) &
    ! Remapped d[uv]dt_dia require east/north halo updates of h
    call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
  call diag_update_remap_grids(cs%diag)

  ! diagnostics
  idt = 1.0 / dt
  if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0)) then
    do j=js,je ; do i=is,ie
      tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
      tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
    enddo ; enddo
    !$OMP parallel do default(shared)
    do k=2,nz ; do j=js,je ; do i=is,ie
      tdif_flx(i,j,k) = (idt * 0.5*(ea_t(i,j,k) + eb_t(i,j,k-1))) * &
                        (tv%T(i,j,k-1) - tv%T(i,j,k))
      tadv_flx(i,j,k) = (idt * (ea_t(i,j,k) - eb_t(i,j,k-1))) * &
                    0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
    enddo ; enddo ; enddo
  endif
  if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0)) then
    do j=js,je ; do i=is,ie
      sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
      sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
    enddo ; enddo
    !$OMP parallel do default(shared)
    do k=2,nz ; do j=js,je ; do i=is,ie
      sdif_flx(i,j,k) = (idt * 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))) * &
                        (tv%S(i,j,k-1) - tv%S(i,j,k))
      sadv_flx(i,j,k) = (idt * (ea_s(i,j,k) - eb_s(i,j,k-1))) * &
                    0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
    enddo ; enddo ; enddo
  endif

  ! mixing of passive tracers from massless boundary layers to interior
  call cpu_clock_begin(id_clock_tracers)

  if (cs%mix_boundary_tracers) then
    tr_ea_bbl = gv%Z_to_H * sqrt(dt*cs%Kd_BBL_tr)
    !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
    do j=js,je
      do i=is,ie
        ebtr(i,j,nz) = eb_s(i,j,nz)
        htot(i) = 0.0
        in_boundary(i) = (g%mask2dT(i,j) > 0.0)
      enddo
      do k=nz,2,-1 ; do i=is,ie
        if (in_boundary(i)) then
          htot(i) = htot(i) + h(i,j,k)
          !   If diapycnal mixing has been suppressed because this is a massless
          ! layer near the bottom, add some mixing of tracers between these
          ! layers.  This flux is based on the harmonic mean of the two
          ! thicknesses, as this corresponds pretty closely (to within
          ! differences in the density jumps between layers) with what is done
          ! in the calculation of the fluxes in the first place.  Kd_min_tr
          ! should be much less than the values that have been set in Kd_int,
          ! perhaps a molecular diffusivity.
          add_ent = ((dt * cs%Kd_min_tr) * gv%Z_to_H**2) * &
                    ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
                     (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
                    0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))
          if (htot(i) < tr_ea_bbl) then
            add_ent = max(0.0, add_ent, &
                          (tr_ea_bbl - htot(i)) - min(ea_s(i,j,k),eb_s(i,j,k-1)))
          elseif (add_ent < 0.0) then
            add_ent = 0.0 ; in_boundary(i) = .false.
          endif

          ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
          eatr(i,j,k) = ea_s(i,j,k) + add_ent
        else
          ebtr(i,j,k-1) = eb_s(i,j,k-1) ; eatr(i,j,k) = ea_s(i,j,k)
        endif

        if (cs%double_diffuse) then ; if (kd_extra_s(i,j,k) > 0.0) then
          add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
             (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) +  h_neglect)
          ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
          eatr(i,j,k) = eatr(i,j,k) + add_ent
        endif ; endif
      enddo ; enddo
      do i=is,ie ; eatr(i,j,1) = ea_s(i,j,1) ; enddo

    enddo

    ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
    ! so h_prebound is used for the old thickness.
    !### I think that in the following, ea_s and eb_s should be eatr and ebtr.
    call call_tracer_column_fns(h_prebound, h, ea_s, eb_s, fluxes, hml, dt, g, gv, us, tv, &
                              cs%optics, cs%tracer_flow_CSp, cs%debug, &
                              evap_cfl_limit = cs%evap_CFL_limit, &
                              minimum_forcing_depth=cs%minimum_forcing_depth)

  elseif (cs%double_diffuse) then  ! extra diffusivity for passive tracers

    do j=js,je ; do i=is,ie
      ebtr(i,j,nz) = eb_s(i,j,nz) ; eatr(i,j,1) = ea_s(i,j,1)
    enddo ; enddo
    !$OMP parallel do default(shared) private(add_ent)
    do k=nz,2,-1 ; do j=js,je ; do i=is,ie
      if (kd_extra_s(i,j,k) > 0.0) then
        add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
           (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + h_neglect)
      else
        add_ent = 0.0
      endif
      ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
      eatr(i,j,k) = ea_s(i,j,k) + add_ent
    enddo ; enddo ; enddo

    ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
    call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
                                cs%optics, cs%tracer_flow_CSp, cs%debug,&
                                evap_cfl_limit = cs%evap_CFL_limit, &
                                minimum_forcing_depth=cs%minimum_forcing_depth)
  else
    ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
    !### eatr and ebtr may not be initialized or may be 0, depending on CS%geothermal.
    call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
                                cs%optics, cs%tracer_flow_CSp, cs%debug, &
                                evap_cfl_limit = cs%evap_CFL_limit, &
                                minimum_forcing_depth=cs%minimum_forcing_depth)
  endif  ! (CS%mix_boundary_tracers)

  call cpu_clock_end(id_clock_tracers)

  ! Apply ALE sponge
  if (cs%use_sponge) then
    call cpu_clock_begin(id_clock_sponge)
    if (associated(cs%ALE_sponge_CSp)) then
      call apply_ale_sponge(h, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)
    endif

    call cpu_clock_end(id_clock_sponge)
    if (cs%debug) then
      call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)
      call mom_thermovar_chksum("apply_sponge ", tv, g)
    endif
  endif ! CS%use_sponge

  call disable_averaging(cs%diag)
  ! Diagnose the diapycnal diffusivities and other related quantities.
  call enable_averages(dt, time_end, cs%diag)

  if (cs%id_Kd_interface > 0) call post_data(cs%id_Kd_interface, kd_int,  cs%diag)
  if (cs%id_Kd_heat      > 0) call post_data(cs%id_Kd_heat,      kd_heat, cs%diag)
  if (cs%id_Kd_salt      > 0) call post_data(cs%id_Kd_salt,      kd_salt, cs%diag)
  if (cs%id_Kd_ePBL      > 0) call post_data(cs%id_Kd_ePBL,      kd_epbl, cs%diag)

  if (cs%id_ea         > 0) call post_data(cs%id_ea,         ea_s, cs%diag)
  if (cs%id_eb         > 0) call post_data(cs%id_eb,         eb_s, cs%diag)
  if (cs%id_ea_t       > 0) call post_data(cs%id_ea_t,       ea_t, cs%diag)
  if (cs%id_eb_t       > 0) call post_data(cs%id_eb_t,       eb_t, cs%diag)
  if (cs%id_ea_s       > 0) call post_data(cs%id_ea_s,       ea_s, cs%diag)
  if (cs%id_eb_s       > 0) call post_data(cs%id_eb_s,       eb_s, cs%diag)

  if (cs%id_dudt_dia > 0) call post_data(cs%id_dudt_dia, adp%du_dt_dia,  cs%diag)
  if (cs%id_dvdt_dia > 0) call post_data(cs%id_dvdt_dia, adp%dv_dt_dia,  cs%diag)
  if (cs%id_wd       > 0) call post_data(cs%id_wd,       cdp%diapyc_vel, cs%diag)

  if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)
  if (cs%id_Tadv > 0) call post_data(cs%id_Tadv, tadv_flx, cs%diag)
  if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)
  if (cs%id_Sadv > 0) call post_data(cs%id_Sadv, sadv_flx, cs%diag)

  !! Diagnostics for terms multiplied by fractional thicknesses
  if (cs%id_hf_dudt_dia_2d > 0) then
    allocate(hf_dudt_dia_2d(g%IsdB:g%IedB,g%jsd:g%jed))
    hf_dudt_dia_2d(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=isq,ieq
      hf_dudt_dia_2d(i,j) = hf_dudt_dia_2d(i,j) + adp%du_dt_dia(i,j,k) * adp%diag_hfrac_u(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_hf_dudt_dia_2d, hf_dudt_dia_2d, cs%diag)
    deallocate(hf_dudt_dia_2d)
  endif

  if (cs%id_hf_dvdt_dia_2d > 0) then
    allocate(hf_dvdt_dia_2d(g%isd:g%ied,g%JsdB:g%JedB))
    hf_dvdt_dia_2d(:,:) = 0.0
    do k=1,nz ; do j=jsq,jeq ; do i=is,ie
      hf_dvdt_dia_2d(i,j) = hf_dvdt_dia_2d(i,j) + adp%dv_dt_dia(i,j,k) * adp%diag_hfrac_v(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_hf_dvdt_dia_2d, hf_dvdt_dia_2d, cs%diag)
    deallocate(hf_dvdt_dia_2d)
  endif

  call disable_averaging(cs%diag)

  if (showcalltree) call calltree_leave("diabatic_ALE_legacy()")

end subroutine diabatic_ale_legacy


!>  This subroutine imposes the diapycnal mass fluxes and the
!!  accompanying diapycnal advection of momentum and tracers.
subroutine diabatic_ale(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
                        G, GV, US, CS, Waves)
  type(ocean_grid_type),                     intent(inout) :: G         !< ocean 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, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u         !< zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v         !< meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h         !< thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),                     intent(inout) :: tv        !< points to thermodynamic fields
                                                                        !! unused have NULL ptrs
  real, dimension(:,:),                      pointer       :: Hml       !< Active mixed layer depth [Z ~> m]
  type(forcing),                             intent(inout) :: fluxes    !< points to forcing fields
                                                                        !! unused fields have NULL ptrs
  type(vertvisc_type),                       intent(inout) :: visc      !< vertical viscosities, BBL properies, and
  type(accel_diag_ptrs),                     intent(inout) :: ADp       !< related points to accelerations in momentum
                                                                        !! equations, to enable the later derived
                                                                        !! diagnostics, like energy budgets
  type(cont_diag_ptrs),                      intent(inout) :: CDp       !< points to terms in continuity equations
  real,                                      intent(in)    :: dt        !< time increment [T ~> s]
  type(time_type),                           intent(in)    :: Time_end  !< Time at the end of the interval
  type(diabatic_cs),                         pointer       :: CS        !< module control structure
  type(wave_parameters_cs),        optional, pointer       :: Waves     !< Surface gravity waves

  ! local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    ea_s,     &  ! amount of fluid entrained from the layer above within
                 ! one time step  [H ~> m or kg m-2]
    eb_s,     &  ! amount of fluid entrained from the layer below within
                 ! one time step  [H ~> m or kg m-2]
    ea_t,     &  ! amount of fluid entrained from the layer above within
                 ! one time step  [H ~> m or kg m-2]
    eb_t,     &  ! amount of fluid entrained from the layer below within
                 ! one time step  [H ~> m or kg m-2]
    h_prebound, & ! initial layer thicknesses [H ~> m or kg m-2]
    dsv_dt, &    ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
    dsv_ds, &    ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
    ctke,   &    ! convective TKE requirements for each layer [R Z3 T-2 ~> J m-2].
    u_h,    &    ! zonal and meridional velocities at thickness points after
    v_h          ! entrainment [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G)) :: &
    SkinBuoyFlux! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag                ! diagnostic array for thickness
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag             ! diagnostic array for temp
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag             ! diagnostic array for salinity
  real, dimension(SZI_(G),SZJ_(G))      :: tendency_2d           ! depth integrated content tendency for diagn

  real :: net_ent  ! The net of ea-eb at an interface.

  real, dimension(SZI_(G),SZJ_(G),SZK_(G))  :: &
    eatr, &  ! The equivalent of ea and eb for tracers, which differ from ea and
    ebtr     ! eb in that they tend to homogenize tracers in massless layers
             ! near the boundaries [H ~> m or kg m-2] (for Bous or non-Bouss)

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    Kd_int,   & ! diapycnal diffusivity of interfaces [Z2 T-1 ~> m2 s-1]
    Kd_heat,  & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1]
    Kd_salt,  & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1]
    Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to
                ! Kd_int [Z2 T-1 ~> m2 s-1].
    kd_extra_s , & !  The extra diffusivity of salinity due to double diffusion relative to
                ! Kd_int [Z2 T-1 ~> m2 s-1].
    kd_epbl,  & ! boundary layer or convective diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]
    tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
    tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
    sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
    sadv_flx    ! advective diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

  real, allocatable, dimension(:,:) :: &
    hf_dudt_dia_2d, hf_dvdt_dia_2d ! Depth sum of diapycnal mixing accelaration * fract. thickness [L T-2 ~> m s-2].

  logical :: in_boundary(szi_(g)) ! True if there are no massive layers below,
                                  ! where massive is defined as sufficiently thick that
                                  ! the no-flux boundary conditions have not restricted
                                  ! the entrainment - usually sqrt(Kd*dt).

  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 :: h_neglect2   ! h_neglect^2 [H2 ~> m2 or kg2 m-4]
  real :: add_ent      ! Entrainment that needs to be added when mixing tracers
                       ! [H ~> m or kg m-2]
  real :: eaval        ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
  real :: hval         ! hval is 2*h at velocity grid points [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 :: Tr_ea_BBL    ! The diffusive tracer thickness in the BBL that is
                       ! coupled to the bottom within a timestep [H ~> m or kg m-2]

  real :: htot(szib_(g))             ! The summed thickness from the bottom [H ~> m or kg m-2].
  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 :: Kd_add_here    ! An added diffusivity [Z2 T-1 ~> m2 s-1].
  real :: Idt     ! The inverse time step [T-1 ~> s-1]

  integer :: dir_flag     ! An integer encoding the directions in which to do halo updates.
  logical :: showCallTree ! If true, show the call tree
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m

  is   = g%isc  ; ie  = g%iec  ; js  = g%jsc  ; je  = g%jec ; nz = g%ke
  isq  = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
  kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0
  ea_s(:,:,:) = 0.0; eb_s(:,:,:) = 0.0; ea_t(:,:,:) = 0.0; eb_t(:,:,:) = 0.0

  showcalltree = calltree_showquery()
  if (showcalltree) call calltree_enter("diabatic_ALE(), MOM_diabatic_driver.F90")

  if (.not. (cs%useALEalgorithm)) call mom_error(fatal, "MOM_diabatic_driver: "// &
         "The ALE algorithm must be enabled when using MOM_diabatic_driver.")

  ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
  call enable_averages(dt, time_end, cs%diag)

  if (cs%use_geothermal) then
    call cpu_clock_begin(id_clock_geothermal)
    call geothermal_in_place(h, tv, dt, g, gv, us, cs%geothermal_CSp, halo=cs%halo_TS_diff)
    call cpu_clock_end(id_clock_geothermal)
    if (showcalltree) call calltree_waypoint("geothermal (diabatic)")
    if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)
  endif

  ! Whenever thickness changes let the diag manager know, target grids
  ! for vertical remapping may need to be regenerated.
  call diag_update_remap_grids(cs%diag)

  ! Set_pen_shortwave estimates the optical properties of the water column.
  ! It will need to be modified later to include information about the
  ! biological properties and layer thicknesses.
  if (associated(cs%optics)) &
    call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)

  if (cs%debug) call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)

  if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then
    if (cs%use_geothermal) then
      ! The presence of eatr and ebtr causes find_uv_at_h to use a tridiagonal solver,
      ! which changes answers at the level of roundoff because ((A*B / A) /= B).
      eatr(:,:,:) = 0.0 ; ebtr(:,:,:) = 0.0
      call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us, eatr, ebtr)
    else
      call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
    endif
    if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
  endif

  call cpu_clock_begin(id_clock_set_diffusivity)
  ! Sets: Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb
  ! Also changes: visc%Kd_shear, visc%Kv_shear and visc%Kv_slow
  if (cs%debug) &
    call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)
  if (cs%double_diffuse) then
    call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, cs%set_diff_CSp, &
                         kd_int=kd_int, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)
  else
    call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, &
                         cs%set_diff_CSp, kd_int=kd_int)
  endif
  call cpu_clock_end(id_clock_set_diffusivity)
  if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")

  if (cs%debug) then
    call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)
    call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)
    call mom_thermovar_chksum("after set_diffusivity ", tv, g)
    call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
  endif

  ! Set diffusivities for heat and salt separately

  ! Add contribution from double diffusion
  if (cs%double_diffuse) then
    !$OMP parallel do default(shared)
    do k=1,nz+1 ; do j=js,je ; do i=is,ie
      kd_salt(i,j,k) = kd_int(i,j,k) + kd_extra_s(i,j,k)
      kd_heat(i,j,k) = kd_int(i,j,k) + kd_extra_t(i,j,k)
    enddo ; enddo ; enddo
  else
    !$OMP parallel do default(shared)
    do k=1,nz+1 ; do j=js,je ; do i=is,ie
      kd_salt(i,j,k) = kd_int(i,j,k)
      kd_heat(i,j,k) = kd_int(i,j,k)
    enddo ; enddo ; enddo
  endif

  if (cs%debug) then
    call hchksum(kd_heat, "after set_diffusivity Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
    call hchksum(kd_salt, "after set_diffusivity Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
  endif

  if (cs%useKPP) then
    call cpu_clock_begin(id_clock_kpp)
    ! total vertical viscosity in the interior is represented via visc%Kv_shear
    do k=1,nz+1 ; do j=js,je ; do i=is,ie
      visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + visc%Kv_slow(i,j,k)
    enddo ; enddo ; enddo

    ! KPP needs the surface buoyancy flux but does not update state variables.
    ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
    ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux
    ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second)
    ! unlike other instances where the fluxes are integrated in time over a time-step.
    call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
                                 cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)

    ! The KPP scheme calculates boundary layer diffusivities and non-local transport.
    call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &
                         fluxes%ustar, cs%KPP_buoy_flux, waves=waves)

    call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
                       kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)

    if (associated(hml)) then
      call kpp_get_bld(cs%KPP_CSp, hml(:,:), g, us)
      call pass_var(hml, g%domain, halo=1)
      ! If visc%MLD exists, copy KPP's BLD into it
      if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
    endif

    call cpu_clock_end(id_clock_kpp)
    if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")
    if (cs%debug) then
      call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)
      call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)
      call mom_thermovar_chksum("after KPP", tv, g)
      call hchksum(kd_heat, "after KPP Kd_heat", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
      call hchksum(kd_salt, "after KPP Kd_salt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
    endif

  endif  ! endif for KPP


  if (cs%useKPP) then
    call cpu_clock_begin(id_clock_kpp)
    if (cs%debug) then
      call hchksum(cs%KPP_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(cs%KPP_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(cs%KPP_NLTheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)
      call hchksum(cs%KPP_NLTscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
    endif
    ! Apply non-local transport of heat and salt
    ! Changes: tv%T, tv%S
    call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat,   cs%KPP_temp_flux, &
                                    us%T_to_s*dt, tv%T, us%Q_to_J_kg*tv%C_p)
    call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, &
                                    us%T_to_s*dt, tv%S)
    call cpu_clock_end(id_clock_kpp)
    if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
    if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)

    if (cs%debug) then
      call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)
      call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)
      call mom_thermovar_chksum("after KPP_applyNLT ", tv, g)
    endif
  endif ! endif for KPP

  ! This is the "old" method for applying differential diffusion.
  ! Changes: tv%T, tv%S
  if (cs%double_diffuse .and. associated(tv%T) .and. (.not.cs%use_CVMix_ddiff)) then

    call cpu_clock_begin(id_clock_differential_diff)
    call differential_diffuse_t_s(h, tv%T, tv%S, kd_extra_t, kd_extra_s, dt, g, gv)
    call cpu_clock_end(id_clock_differential_diff)

    if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")
    if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g, gv, us)

    ! increment heat and salt diffusivity.
    ! CS%useKPP==.true. already has extra_T and extra_S included
    if (.not. cs%useKPP) then
      !$OMP parallel do default(shared)
      do k=2,nz ; do j=js,je ; do i=is,ie
        kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)
        kd_salt(i,j,k) = kd_salt(i,j,k) + kd_extra_s(i,j,k)
      enddo ; enddo ; enddo
    endif

  endif

  ! Calculate vertical mixing due to convection (computed via CVMix)
  if (cs%use_CVMix_conv) then
    ! Increment vertical diffusion and viscosity due to convection
    if (cs%useKPP) then
      call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml, kd=kd_heat, kv=visc%Kv_shear, kd_aux=kd_salt)
    else
      call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml, kd=kd_heat, kv=visc%Kv_slow, kd_aux=kd_salt)
    endif
  endif

  ! Save fields before boundary forcing is applied for tendency diagnostics
  if (cs%boundary_forcing_tendency_diag) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      h_diag(i,j,k)    = h(i,j,k)
      temp_diag(i,j,k) = tv%T(i,j,k)
      saln_diag(i,j,k) = tv%S(i,j,k)
    enddo ; enddo ; enddo
  endif

  ! Apply forcing
  call cpu_clock_begin(id_clock_remap)

  ! Changes made to following fields:  h, tv%T and tv%S.
  do k=1,nz ; do j=js,je ; do i=is,ie
    h_prebound(i,j,k) = h(i,j,k)
  enddo ; enddo ; enddo
  if (cs%use_energetic_PBL) then

    skinbuoyflux(:,:) = 0.0
    call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
            optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, cs%evap_CFL_limit, &
            cs%minimum_forcing_depth, ctke, dsv_dt, dsv_ds, skinbuoyflux=skinbuoyflux)

    if (cs%debug) then
      call hchksum(ea_t, "after applyBoundaryFluxes ea_t", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_t, "after applyBoundaryFluxes eb_t", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(ea_s, "after applyBoundaryFluxes ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_s, "after applyBoundaryFluxes eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(ctke, "after applyBoundaryFluxes cTKE", g%HI, haloshift=0, &
                   scale=us%RZ3_T3_to_W_m2*us%T_to_s)
      call hchksum(dsv_dt, "after applyBoundaryFluxes dSV_dT", g%HI, haloshift=0, scale=us%kg_m3_to_R)
      call hchksum(dsv_ds, "after applyBoundaryFluxes dSV_dS", g%HI, haloshift=0, scale=us%kg_m3_to_R)
    endif

    call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
    call energetic_pbl(h, u_h, v_h, tv, fluxes, dt, kd_epbl, g, gv, us, &
                       cs%energetic_PBL_CSp, dsv_dt, dsv_ds, ctke, skinbuoyflux, waves=waves)

    if (associated(hml)) then
      call energetic_pbl_get_mld(cs%energetic_PBL_CSp, hml(:,:), g, us)
      call pass_var(hml, g%domain, halo=1)
      ! If visc%MLD exists, copy ePBL's MLD into it
      if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
    elseif (associated(visc%MLD)) then
      call energetic_pbl_get_mld(cs%energetic_PBL_CSp, visc%MLD, g, us)
      call pass_var(visc%MLD, g%domain, halo=1)
    endif

    ! Augment the diffusivities and viscosity due to those diagnosed in energetic_PBL.
    do k=2,nz ; do j=js,je ; do i=is,ie
      if (cs%ePBL_is_additive) then
        kd_add_here = kd_epbl(i,j,k)
        visc%Kv_shear(i,j,k) = visc%Kv_shear(i,j,k) + cs%ePBL_Prandtl*kd_epbl(i,j,k)
      else
        kd_add_here = max(kd_epbl(i,j,k) - visc%Kd_shear(i,j,k), 0.0)
        visc%Kv_shear(i,j,k) = max(visc%Kv_shear(i,j,k), cs%ePBL_Prandtl*kd_epbl(i,j,k))
      endif

      kd_heat(i,j,k) = kd_heat(i,j,k) + kd_add_here
      kd_salt(i,j,k) = kd_salt(i,j,k) + kd_add_here

    enddo ; enddo ; enddo

    if (cs%debug) then
      call hchksum(ea_t, "after ePBL ea_t", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_t, "after ePBL eb_t", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(ea_s, "after ePBL ea_s", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_s, "after ePBL eb_s", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(kd_epbl, "after ePBL Kd_ePBL", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
    endif

  else
    call applyboundaryfluxesinout(cs%diabatic_aux_CSp, g, gv, us, dt, fluxes, cs%optics, &
                                  optics_nbands(cs%optics), h, tv, cs%aggregate_FW_forcing, &
                                  cs%evap_CFL_limit, cs%minimum_forcing_depth)

  endif   ! endif for CS%use_energetic_PBL

  ! diagnose the tendencies due to boundary forcing
  ! At this point, the diagnostic grids have not been updated since the call to the boundary layer scheme
  !  so all tendency diagnostics need to be posted on h_diag, and grids rebuilt afterwards
  if (cs%boundary_forcing_tendency_diag) then
    call diagnose_boundary_forcing_tendency(tv, h, temp_diag, saln_diag, h_diag, dt, g, gv, us, cs)
    if (cs%id_boundary_forcing_h > 0) call post_data(cs%id_boundary_forcing_h, h, cs%diag, alt_h = h_diag)
  endif
  ! Boundary fluxes may have changed T, S, and h
  call diag_update_remap_grids(cs%diag)
  call cpu_clock_end(id_clock_remap)
  if (cs%debug) then
    call mom_forcing_chksum("after applyBoundaryFluxes ", fluxes, g, us, haloshift=0)
    call mom_thermovar_chksum("after applyBoundaryFluxes ", tv, g)
    call mom_state_chksum("after applyBoundaryFluxes ", u, v, h, g, gv, us, haloshift=0)
  endif
  if (showcalltree) call calltree_waypoint("done with applyBoundaryFluxes (diabatic)")
  if (cs%debugConservation)  call mom_state_stats('applyBoundaryFluxes', u, v, h, tv%T, tv%S, g, gv, us)

  if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")
  if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)

  ! calculate change in temperature & salinity due to dia-coordinate surface diffusion
  if (associated(tv%T)) then

    if (cs%debug) then
      call hchksum(ea_t, "before triDiagTS ea_t ", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_t, "before triDiagTS eb_t ", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(ea_s, "before triDiagTS ea_s ", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(eb_s, "before triDiagTS eb_s ", g%HI, haloshift=0, scale=gv%H_to_m)
    endif

    call cpu_clock_begin(id_clock_tridiag)
    !  Keep salinity from falling below a small but positive threshold.
    !  This constraint is needed for SIS1 ice model, which can extract
    !  more salt than is present in the ocean. SIS2 does not suffer
    !  from this limitation, in which case we can let salinity=0 and still
    !  have salt conserved with SIS2 ice. So for SIS2, we can run with
    !  BOUND_SALINITY=False in MOM.F90.
    if (associated(tv%S) .and. associated(tv%salt_deficit)) &
      call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

    if (cs%diabatic_diff_tendency_diag) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        temp_diag(i,j,k) = tv%T(i,j,k)
        saln_diag(i,j,k) = tv%S(i,j,k)
      enddo ; enddo ; enddo
    endif

    ! set ea_t=eb_t=Kd_heat and ea_s=eb_s=Kd_salt on interfaces for use in the tridiagonal solver.
    do j=js,je ; do i=is,ie
      ea_t(i,j,1) = 0.; ea_s(i,j,1) = 0.
    enddo ; enddo

    !$OMP parallel do default(shared) private(hval)
    do k=2,nz ; do j=js,je ; do i=is,ie
      hval = 1.0 / (h_neglect + 0.5*(h(i,j,k-1) + h(i,j,k)))
      ea_t(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_heat(i,j,k)
      eb_t(i,j,k-1) = ea_t(i,j,k)
      ea_s(i,j,k) = (gv%Z_to_H**2) * dt * hval * kd_salt(i,j,k)
      eb_s(i,j,k-1) = ea_s(i,j,k)
    enddo ; enddo ; enddo
    do j=js,je ; do i=is,ie
      eb_t(i,j,nz) = 0. ; eb_s(i,j,nz) = 0.
    enddo ; enddo
    if (showcalltree) call calltree_waypoint("done setting ea_t,ea_s,eb_t,eb_s from Kd_heat" //&
       "and Kd_salt (diabatic)")

    ! Initialize halo regions of ea_t, eb_t, ea_s and eb_s to default values.
    !$OMP parallel do default(shared)
    do k=1,nz
      do i=is-1,ie+1
        ea_t(i,js-1,k) = 0.0 ; eb_t(i,js-1,k) = 0.0
        ea_s(i,js-1,k) = 0.0 ; eb_s(i,js-1,k) = 0.0
        ea_t(i,je+1,k) = 0.0 ; eb_t(i,je+1,k) = 0.0
        ea_s(i,je+1,k) = 0.0 ; eb_s(i,je+1,k) = 0.0
      enddo
      do j=js,je
        ea_t(is-1,j,k) = 0.0 ; eb_t(is-1,j,k) = 0.0
        ea_s(is-1,j,k) = 0.0 ; eb_s(is-1,j,k) = 0.0
        ea_t(ie+1,j,k) = 0.0 ; eb_t(ie+1,j,k) = 0.0
        ea_s(ie+1,j,k) = 0.0 ; eb_s(ie+1,j,k) = 0.0
      enddo
    enddo

  ! Changes T and S via the tridiagonal solver; no change to h
    call tracer_vertdiff(h, ea_t, eb_t, dt, tv%T, g, gv)
    call tracer_vertdiff(h, ea_s, eb_s, dt, tv%S, g, gv)


    ! In ALE-mode, layer thicknesses do not change. Therefore, we can use h below
    if (cs%diabatic_diff_tendency_diag) then
      call diagnose_diabatic_diff_tendency(tv, h, temp_diag, saln_diag, dt, g, gv, us, cs)
    endif
    call cpu_clock_end(id_clock_tridiag)

    if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")

  endif  ! endif corresponding to if (associated(tv%T))

  if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)

  if (cs%debug) then
    call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)
    call mom_thermovar_chksum("after mixed layer ", tv, g)
  endif

  ! Whenever thickness changes let the diag manager know, as the
  ! target grids for vertical remapping may need to be regenerated.
  call diag_update_remap_grids(cs%diag)

  ! diagnostics
  idt = 1.0 / dt
  if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0)) then
    do j=js,je ; do i=is,ie
      tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
      tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
    enddo ; enddo
    !$OMP parallel do default(shared)
    do k=2,nz ; do j=js,je ; do i=is,ie
      tdif_flx(i,j,k) = (idt * 0.5*(ea_t(i,j,k) + eb_t(i,j,k-1))) * &
                        (tv%T(i,j,k-1) - tv%T(i,j,k))
      tadv_flx(i,j,k) = (idt * (ea_t(i,j,k) - eb_t(i,j,k-1))) * &
                    0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
    enddo ; enddo ; enddo
  endif
  if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0)) then
    do j=js,je ; do i=is,ie
      sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
      sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
    enddo ; enddo
    !$OMP parallel do default(shared)
    do k=2,nz ; do j=js,je ; do i=is,ie
      sdif_flx(i,j,k) = (idt * 0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))) * &
                        (tv%S(i,j,k-1) - tv%S(i,j,k))
      sadv_flx(i,j,k) = (idt * (ea_s(i,j,k) - eb_s(i,j,k-1))) * &
                    0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
    enddo ; enddo ; enddo
  endif

  ! mixing of passive tracers from massless boundary layers to interior
  call cpu_clock_begin(id_clock_tracers)

  if (cs%mix_boundary_tracers) then
    tr_ea_bbl = gv%Z_to_H * sqrt(dt*cs%Kd_BBL_tr)
    !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
    do j=js,je
      do i=is,ie
        ebtr(i,j,nz) = eb_s(i,j,nz)
        htot(i) = 0.0
        in_boundary(i) = (g%mask2dT(i,j) > 0.0)
      enddo
      do k=nz,2,-1 ; do i=is,ie
        if (in_boundary(i)) then
          htot(i) = htot(i) + h(i,j,k)
          !   If diapycnal mixing has been suppressed because this is a massless
          ! layer near the bottom, add some mixing of tracers between these
          ! layers.  This flux is based on the harmonic mean of the two
          ! thicknesses, as this corresponds pretty closely (to within
          ! differences in the density jumps between layers) with what is done
          ! in the calculation of the fluxes in the first place.  Kd_min_tr
          ! should be much less than the values that have been set in Kd_int,
          ! perhaps a molecular diffusivity.
          add_ent = ((dt * cs%Kd_min_tr) * gv%Z_to_H**2) * &
                    ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
                     (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
                    0.5*(ea_s(i,j,k) + eb_s(i,j,k-1))
          if (htot(i) < tr_ea_bbl) then
            add_ent = max(0.0, add_ent, &
                          (tr_ea_bbl - htot(i)) - min(ea_s(i,j,k),eb_s(i,j,k-1)))
          elseif (add_ent < 0.0) then
            add_ent = 0.0 ; in_boundary(i) = .false.
          endif

          ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
          eatr(i,j,k) = ea_s(i,j,k) + add_ent
        else
          ebtr(i,j,k-1) = eb_s(i,j,k-1) ; eatr(i,j,k) = ea_s(i,j,k)
        endif

        if (cs%double_diffuse) then ; if (kd_extra_s(i,j,k) > 0.0) then
          add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
             (0.5 * (h(i,j,k-1) + h(i,j,k)) + &
              h_neglect)
          ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
          eatr(i,j,k) = eatr(i,j,k) + add_ent
        endif ; endif
      enddo ; enddo
      do i=is,ie ; eatr(i,j,1) = ea_s(i,j,1) ; enddo

    enddo

    ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
    ! so use h_prebound as the old value.
    !### I think that in the following, ea_s and eb_s should be eatr and ebtr.
    call call_tracer_column_fns(h_prebound, h, ea_s, eb_s, fluxes, hml, dt, g, gv, us, tv, &
                              cs%optics, cs%tracer_flow_CSp, cs%debug, &
                              evap_cfl_limit = cs%evap_CFL_limit, &
                              minimum_forcing_depth=cs%minimum_forcing_depth)

  elseif (cs%double_diffuse) then  ! extra diffusivity for passive tracers

    do j=js,je ; do i=is,ie
      ebtr(i,j,nz) = eb_s(i,j,nz) ; eatr(i,j,1) = ea_s(i,j,1)
    enddo ; enddo
    !$OMP parallel do default(shared) private(add_ent)
    do k=nz,2,-1 ; do j=js,je ; do i=is,ie
      if (kd_extra_s(i,j,k) > 0.0) then
        add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
           (0.5 * (h(i,j,k-1) + h(i,j,k))  + &
            h_neglect)
      else
        add_ent = 0.0
      endif
      ebtr(i,j,k-1) = eb_s(i,j,k-1) + add_ent
      eatr(i,j,k) = ea_s(i,j,k) + add_ent
    enddo ; enddo ; enddo

    ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
    call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
                                cs%optics, cs%tracer_flow_CSp, cs%debug,&
                                evap_cfl_limit = cs%evap_CFL_limit, &
                                minimum_forcing_depth=cs%minimum_forcing_depth)
  else
    ! For passive tracers, the changes in thickness due to boundary fluxes has yet to be applied
    !### eatr and ebtr may not be initialized or may be 0, depending on CS%geothermal.
    call call_tracer_column_fns(h_prebound, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
                                cs%optics, cs%tracer_flow_CSp, cs%debug, &
                                evap_cfl_limit = cs%evap_CFL_limit, &
                                minimum_forcing_depth=cs%minimum_forcing_depth)
  endif  ! (CS%mix_boundary_tracers)

  call cpu_clock_end(id_clock_tracers)

  ! Apply ALE sponge
  if (cs%use_sponge) then
    call cpu_clock_begin(id_clock_sponge)
    if (associated(cs%ALE_sponge_CSp)) then
      call apply_ale_sponge(h, dt, g, gv, us, cs%ALE_sponge_CSp, cs%Time)
    endif

    call cpu_clock_end(id_clock_sponge)
    if (cs%debug) then
      call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)
      call mom_thermovar_chksum("apply_sponge ", tv, g)
    endif
  endif ! CS%use_sponge

  call cpu_clock_begin(id_clock_pass)
  if (g%symmetric) then ; dir_flag = to_all+omit_corners
  else ; dir_flag = to_west+to_south+omit_corners ; endif
  call create_group_pass(cs%pass_hold_eb_ea, eb_t, g%Domain, dir_flag, halo=1)
  call create_group_pass(cs%pass_hold_eb_ea, eb_s, g%Domain, dir_flag, halo=1)
  call create_group_pass(cs%pass_hold_eb_ea, ea_t, g%Domain, dir_flag, halo=1)
  call create_group_pass(cs%pass_hold_eb_ea, ea_s, g%Domain, dir_flag, halo=1)
  call do_group_pass(cs%pass_hold_eb_ea, g%Domain)
  ! visc%Kv_slow is not in the group pass because it has larger vertical extent.
  if (associated(visc%Kv_slow)) &
    call pass_var(visc%Kv_slow, g%Domain, to_all+omit_corners, halo=1)
  call cpu_clock_end(id_clock_pass)

  call disable_averaging(cs%diag)

  ! Diagnose the diapycnal diffusivities and other related quantities.
  call enable_averages(dt, time_end, cs%diag)

  if (cs%id_Kd_interface > 0) call post_data(cs%id_Kd_interface, kd_int,  cs%diag)
  if (cs%id_Kd_heat      > 0) call post_data(cs%id_Kd_heat,      kd_heat, cs%diag)
  if (cs%id_Kd_salt      > 0) call post_data(cs%id_Kd_salt,      kd_salt, cs%diag)
  if (cs%id_Kd_ePBL      > 0) call post_data(cs%id_Kd_ePBL,      kd_epbl, cs%diag)

  if (cs%id_ea_t       > 0) call post_data(cs%id_ea_t,       ea_t, cs%diag)
  if (cs%id_eb_t       > 0) call post_data(cs%id_eb_t,       eb_t, cs%diag)
  if (cs%id_ea_s       > 0) call post_data(cs%id_ea_s,       ea_s, cs%diag)
  if (cs%id_eb_s       > 0) call post_data(cs%id_eb_s,       eb_s, cs%diag)

  if (cs%id_dudt_dia > 0) call post_data(cs%id_dudt_dia, adp%du_dt_dia,  cs%diag)
  if (cs%id_dvdt_dia > 0) call post_data(cs%id_dvdt_dia, adp%dv_dt_dia,  cs%diag)

  if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)
  if (cs%id_Tadv > 0) call post_data(cs%id_Tadv, tadv_flx, cs%diag)
  if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)
  if (cs%id_Sadv > 0) call post_data(cs%id_Sadv, sadv_flx, cs%diag)

  !! Diagnostics for terms multiplied by fractional thicknesses
  if (cs%id_hf_dudt_dia_2d > 0) then
    allocate(hf_dudt_dia_2d(g%IsdB:g%IedB,g%jsd:g%jed))
    hf_dudt_dia_2d(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=isq,ieq
      hf_dudt_dia_2d(i,j) = hf_dudt_dia_2d(i,j) + adp%du_dt_dia(i,j,k) * adp%diag_hfrac_u(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_hf_dudt_dia_2d, hf_dudt_dia_2d, cs%diag)
    deallocate(hf_dudt_dia_2d)
  endif

  if (cs%id_hf_dvdt_dia_2d > 0) then
    allocate(hf_dvdt_dia_2d(g%isd:g%ied,g%JsdB:g%JedB))
    hf_dvdt_dia_2d(:,:) = 0.0
    do k=1,nz ; do j=jsq,jeq ; do i=is,ie
      hf_dvdt_dia_2d(i,j) = hf_dvdt_dia_2d(i,j) + adp%dv_dt_dia(i,j,k) * adp%diag_hfrac_v(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_hf_dvdt_dia_2d, hf_dvdt_dia_2d, cs%diag)
    deallocate(hf_dvdt_dia_2d)
  endif

  call disable_averaging(cs%diag)

  if (showcalltree) call calltree_leave("diabatic_ALE()")

end subroutine diabatic_ale

!> Imposes the diapycnal mass fluxes and the accompanying diapycnal advection of momentum and tracers
!! using the original MOM6 algorithms.
subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end, &
                            G, GV, US, CS, WAVES)
  type(ocean_grid_type),                     intent(inout) :: G         !< ocean 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, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: u         !< zonal velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: v         !< meridional velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: h         !< thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),                     intent(inout) :: tv        !< points to thermodynamic fields
                                                                        !! unused have NULL ptrs
  real, dimension(:,:),                      pointer       :: Hml       !< Active mixed layer depth [Z ~> m]
  type(forcing),                             intent(inout) :: fluxes    !< points to forcing fields
                                                                        !! unused fields have NULL ptrs
  type(vertvisc_type),                       intent(inout) :: visc      !< vertical viscosities, BBL properies, and
  type(accel_diag_ptrs),                     intent(inout) :: ADp       !< related points to accelerations in momentum
                                                                        !! equations, to enable the later derived
                                                                        !! diagnostics, like energy budgets
  type(cont_diag_ptrs),                      intent(inout) :: CDp       !< points to terms in continuity equations
  real,                                      intent(in)    :: dt        !< time increment [T ~> s]
  type(time_type),                           intent(in)    :: Time_end  !< Time at the end of the interval
  type(diabatic_cs),                         pointer       :: CS        !< module control structure
  type(wave_parameters_cs),        optional, pointer       :: Waves     !< Surface gravity waves

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    ea,     &    ! amount of fluid entrained from the layer above within
                 ! one time step  [H ~> m or kg m-2]
    eb,     &    ! amount of fluid entrained from the layer below within
                 ! one time step  [H ~> m or kg m-2]
    kd_lay, &    ! diapycnal diffusivity of layers [Z2 T-1 ~> m2 s-1]
    h_orig, &    ! initial layer thicknesses [H ~> m or kg m-2]
    h_prebound, & ! initial layer thicknesses [H ~> m or kg m-2]
    hold,   &    ! layer thickness before diapycnal entrainment, and later
                 ! the initial layer thicknesses (if a mixed layer is used),
                 ! [H ~> m or kg m-2]
    dsv_dt, &    ! The partial derivative of specific volume with temperature [R-1 degC-1 ~> m3 kg-1 degC-1]
    dsv_ds, &    ! The partial derivative of specific volume with salinity [R-1 ppt-1 ~> m3 kg-1 ppt-1].
    u_h,    &    ! zonal and meridional velocities at thickness points after
    v_h          ! entrainment [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G)) :: &
    Rcv_ml, &   ! coordinate density of mixed layer, used for applying sponges
    SkinBuoyFlux! 2d surface buoyancy flux [Z2 T-3 ~> m2 s-3], used by ePBL
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: h_diag                ! diagnostic array for thickness
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: temp_diag             ! diagnostic array for temp
  real, dimension(SZI_(G),SZJ_(G),G%ke) :: saln_diag             ! diagnostic array for salinity
  real, dimension(SZI_(G),SZJ_(G))      :: tendency_2d           ! depth integrated content tendency for diagn

  real :: net_ent  ! The net of ea-eb at an interface.

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), target :: &
             ! These are targets so that the space can be shared with eaml & ebml.
    eatr, &  ! The equivalent of ea and eb for tracers, which differ from ea and
    ebtr     ! eb in that they tend to homogenize tracers in massless layers
             ! near the boundaries [H ~> m or kg m-2] (for Bous or non-Bouss)

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    Kd_int,   & ! diapycnal diffusivity of interfaces [Z2 T-1 ~> m2 s-1]
    Kd_heat,  & ! diapycnal diffusivity of heat [Z2 T-1 ~> m2 s-1]
    Kd_salt,  & ! diapycnal diffusivity of salt and passive tracers [Z2 T-1 ~> m2 s-1]
    Kd_extra_T , & ! The extra diffusivity of temperature due to double diffusion relative to
                ! Kd_int [Z2 T-1 ~> m2 s-1].
    kd_extra_s , & !  The extra diffusivity of salinity due to double diffusion relative to
                ! Kd_int [Z2 T-1 ~> m2 s-1].
    kd_epbl,  & ! test array of diapycnal diffusivities at interfaces [Z2 T-1 ~> m2 s-1]
    tdif_flx, & ! diffusive diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
    tadv_flx, & ! advective diapycnal heat flux across interfaces [degC H T-1 ~> degC m s-1 or degC kg m-2 s-1]
    sdif_flx, & ! diffusive diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
    sadv_flx    ! advective diapycnal salt flux across interfaces [ppt H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]

  real, allocatable, dimension(:,:) :: &
    hf_dudt_dia_2d, hf_dvdt_dia_2d ! Depth sum of diapycnal mixing accelaration * fract. thickness [L T-2 ~> m s-2].

  ! The following 5 variables are only used with a bulk mixed layer.
  real, pointer, dimension(:,:,:) :: &
    eaml, &  ! The equivalent of ea due to mixed layer processes [H ~> m or kg m-2].
    ebml     ! The equivalent of eb due to mixed layer processes [H ~> m or kg m-2].
             ! eaml and ebml are pointers to eatr and ebtr so as to reuse the memory as
             ! the arrays are not needed at the same time.

  integer :: kb(szi_(g),szj_(g)) ! index of the lightest layer denser
                                 ! than the buffer layer [nondim]

  real :: p_ref_cv(szi_(g))      ! Reference pressure for the potential density that defines the
                                 ! coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].

  logical :: in_boundary(szi_(g)) ! True if there are no massive layers below,
                                  ! where massive is defined as sufficiently thick that
                                  ! the no-flux boundary conditions have not restricted
                                  ! the entrainment - usually sqrt(Kd*dt).

  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 :: h_neglect2   ! h_neglect^2 [H2 ~> m2 or kg2 m-4]
  real :: add_ent      ! Entrainment that needs to be added when mixing tracers
                       ! [H ~> m or kg m-2]
  real :: eaval        ! eaval is 2*ea at velocity grid points [H ~> m or kg m-2]
  real :: hval         ! hval is 2*h at velocity grid points [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 :: Tr_ea_BBL    ! The diffusive tracer thickness in the BBL that is
                       ! coupled to the bottom within a timestep [H ~> m or kg m-2]

  real :: htot(szib_(g))             ! The summed thickness from the bottom [H ~> m or kg m-2].
  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 :: dt_mix  ! The amount of time over which to apply mixing [T ~> s]
  real :: Idt     ! The inverse time step [T-1 ~> s-1]
  real :: Idt_accel  ! The inverse time step times rescaling factors [T-1 ~> s-1]

  integer :: dir_flag     ! An integer encoding the directions in which to do halo updates.
  logical :: showCallTree ! If true, show the call tree
  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb, m, halo

  is   = g%isc  ; ie  = g%iec  ; js  = g%jsc  ; je  = g%jec ; nz = g%ke
  isq  = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  nkmb = gv%nk_rho_varies
  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect*h_neglect
  kd_heat(:,:,:) = 0.0 ; kd_salt(:,:,:) = 0.0


  showcalltree = calltree_showquery()
  if (showcalltree) call calltree_enter("layered_diabatic(), MOM_diabatic_driver.F90")

  ! set equivalence between the same bits of memory for these arrays
  eaml => eatr ; ebml => ebtr

  ! For all other diabatic subroutines, the averaging window should be the entire diabatic timestep
  call enable_averages(dt, time_end, cs%diag)

  if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
    halo = cs%halo_TS_diff
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=js-halo,je+halo ; do i=is-halo,ie+halo
      h_orig(i,j,k) = h(i,j,k) ; eaml(i,j,k) = 0.0 ; ebml(i,j,k) = 0.0
    enddo ; enddo ; enddo
  endif

  if (cs%use_geothermal) then
    call cpu_clock_begin(id_clock_geothermal)
    call geothermal_entraining(h, tv, dt, eaml, ebml, g, gv, us, cs%geothermal_CSp, halo=cs%halo_TS_diff)
    call cpu_clock_end(id_clock_geothermal)
    if (showcalltree) call calltree_waypoint("geothermal (diabatic)")
    if (cs%debugConservation) call mom_state_stats('geothermal', u, v, h, tv%T, tv%S, g, gv, us)
  endif

  ! Whenever thickness changes let the diag manager know, target grids
  ! for vertical remapping may need to be regenerated.
  call diag_update_remap_grids(cs%diag)

  ! Set_pen_shortwave estimates the optical properties of the water column.
  ! It will need to be modified later to include information about the
  ! biological properties and layer thicknesses.
  if (associated(cs%optics)) &
    call set_pen_shortwave(cs%optics, fluxes, g, gv, us, cs%diabatic_aux_CSp, cs%opacity_CSp, cs%tracer_flow_CSp)

  if (cs%bulkmixedlayer) then
    if (cs%debug) call mom_forcing_chksum("Before mixedlayer", fluxes, g, us, haloshift=0)

    if (cs%ML_mix_first > 0.0) then
!  This subroutine
!    (1) Cools the mixed layer.
!    (2) Performs convective adjustment by mixed layer entrainment.
!    (3) Heats the mixed layer and causes it to detrain to
!        Monin-Obukhov depth or minimum mixed layer depth.
!    (4) Uses any remaining TKE to drive mixed layer entrainment.
!    (5) Possibly splits buffer layer into two isopycnal layers (when using isopycnal coordinate)
      call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)

      call cpu_clock_begin(id_clock_mixedlayer)
      if (cs%ML_mix_first < 1.0) then
        ! Changes: h, tv%T, tv%S, eaml and ebml  (G is also inout???)
        call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt*cs%ML_mix_first, &
                            eaml,ebml, g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
                            hml, cs%aggregate_FW_forcing, dt, last_call=.false.)
        ! Changes: h, tv%T, tv%S, eaml and ebml  (G is also inout???)
        call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt, eaml, ebml, &
                            g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
                            hml, cs%aggregate_FW_forcing, dt, last_call=.true.)
      endif

      !  Keep salinity from falling below a small but positive threshold.
      !  This constraint is needed for SIS1 ice model, which can extract
      !  more salt than is present in the ocean. SIS2 does not suffer
      !  from this limitation, in which case we can let salinity=0 and still
      !  have salt conserved with SIS2 ice. So for SIS2, we can run with
      !  BOUND_SALINITY=False in MOM.F90.
      if (associated(tv%S) .and. associated(tv%salt_deficit)) &
        call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)
      call cpu_clock_end(id_clock_mixedlayer)
      if (cs%debug) then
        call mom_state_chksum("After mixedlayer ", u, v, h, g, gv, us, haloshift=0)
        call mom_forcing_chksum("After mixedlayer", fluxes, g, us, haloshift=0)
      endif
      if (showcalltree) call calltree_waypoint("done with 1st bulkmixedlayer (diabatic)")
      if (cs%debugConservation) call mom_state_stats('1st bulkmixedlayer', u, v, h, tv%T, tv%S, g, gv, us)
    endif
  endif

  if (cs%debug) &
    call mom_state_chksum("before find_uv_at_h", u, v, h, g, gv, us, haloshift=0)
  if (cs%use_kappa_shear .or. cs%use_CVMix_shear) then
    if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
      call find_uv_at_h(u, v, h_orig, u_h, v_h, g, gv, us, eaml, ebml)
      if (cs%debug) then
        call hchksum(eaml, "after find_uv_at_h eaml", g%HI, scale=gv%H_to_m)
        call hchksum(ebml, "after find_uv_at_h ebml", g%HI, scale=gv%H_to_m)
      endif
    else
      call find_uv_at_h(u, v, h, u_h, v_h, g, gv, us)
    endif
    if (showcalltree) call calltree_waypoint("done with find_uv_at_h (diabatic)")
  endif

  call cpu_clock_begin(id_clock_set_diffusivity)
  ! Sets: Kd_lay, Kd_int, Kd_extra_T, Kd_extra_S and visc%TKE_turb
  ! Also changes: visc%Kd_shear and visc%Kv_shear
  if ((cs%halo_TS_diff > 0) .and. (cs%ML_mix_first > 0.0)) then
    if (associated(tv%T)) call pass_var(tv%T, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
    if (associated(tv%T)) call pass_var(tv%S, g%Domain, halo=cs%halo_TS_diff, complete=.false.)
    call pass_var(h, g%domain, halo=cs%halo_TS_diff, complete=.true.)
  endif
  if (cs%debug) &
    call mom_state_chksum("before set_diffusivity", u, v, h, g, gv, us, haloshift=cs%halo_TS_diff)
  if (cs%double_diffuse) then
    call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, cs%set_diff_CSp, &
                         kd_lay=kd_lay, kd_int=kd_int, kd_extra_t=kd_extra_t, kd_extra_s=kd_extra_s)
  else
    call set_diffusivity(u, v, h, u_h, v_h, tv, fluxes, cs%optics, visc, dt, g, gv, us, &
                         cs%set_diff_CSp, kd_lay=kd_lay, kd_int=kd_int)
  endif
  call cpu_clock_end(id_clock_set_diffusivity)
  if (showcalltree) call calltree_waypoint("done with set_diffusivity (diabatic)")

  if (cs%debug) then
    call mom_state_chksum("after set_diffusivity ", u, v, h, g, gv, us, haloshift=0)
    call mom_forcing_chksum("after set_diffusivity ", fluxes, g, us, haloshift=0)
    call mom_thermovar_chksum("after set_diffusivity ", tv, g)
    call hchksum(kd_lay, "after set_diffusivity Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
    call hchksum(kd_int, "after set_diffusivity Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
  endif


  if (cs%useKPP) then
    call cpu_clock_begin(id_clock_kpp)
    ! KPP needs the surface buoyancy flux but does not update state variables.
    ! We could make this call higher up to avoid a repeat unpacking of the surface fluxes.
    ! Sets: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux
    ! NOTE: CS%KPP_buoy_flux, CS%KPP_temp_flux, CS%KPP_salt_flux are returned as rates (i.e. stuff per second)
    ! unlike other instances where the fluxes are integrated in time over a time-step.
    call calculatebuoyancyflux2d(g, gv, us, fluxes, cs%optics, h, tv%T, tv%S, tv, &
                                 cs%KPP_buoy_flux, cs%KPP_temp_flux, cs%KPP_salt_flux)
    ! The KPP scheme calculates boundary layer diffusivities and non-local transport.

    ! Set diffusivities for heat and salt separately

    if (cs%double_diffuse) then
      ! Add contribution from double diffusion
      !$OMP parallel do default(shared)
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_salt(i,j,k) = kd_int(i,j,k) + kd_extra_s(i,j,k)
        kd_heat(i,j,k) = kd_int(i,j,k) + kd_extra_t(i,j,k)
      enddo ; enddo ; enddo
    else
      !$OMP parallel do default(shared)
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_salt(i,j,k) = kd_int(i,j,k)
        kd_heat(i,j,k) = kd_int(i,j,k)
      enddo ; enddo ; enddo
    endif

    call kpp_compute_bld(cs%KPP_CSp, g, gv, us, h, tv%T, tv%S, u, v, tv, &
                         fluxes%ustar, cs%KPP_buoy_flux, waves=waves)

    call kpp_calculate(cs%KPP_CSp, g, gv, us, h, fluxes%ustar, cs%KPP_buoy_flux, kd_heat, &
                       kd_salt, visc%Kv_shear, cs%KPP_NLTheat, cs%KPP_NLTscalar, waves=waves)

    if (associated(hml)) then
      call kpp_get_bld(cs%KPP_CSp, hml(:,:), g, us)
      call pass_var(hml, g%domain, halo=1)
      ! If visc%MLD exists, copy KPP's BLD into it
      if (associated(visc%MLD)) visc%MLD(:,:) = hml(:,:)
    endif

    if (.not. cs%KPPisPassive) then
      !$OMP parallel do default(shared)
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        kd_int(i,j,k) = min( kd_salt(i,j,k),  kd_heat(i,j,k) )
      enddo ; enddo ; enddo
      if (cs%double_diffuse) then
        !$OMP parallel do default(shared)
        do k=1,nz+1 ; do j=js,je ; do i=is,ie
          kd_extra_s(i,j,k) = (kd_salt(i,j,k) - kd_int(i,j,k))
          kd_extra_t(i,j,k) = (kd_heat(i,j,k) - kd_int(i,j,k))
        enddo ; enddo ; enddo
      endif
    endif ! not passive

    call cpu_clock_end(id_clock_kpp)
    if (showcalltree) call calltree_waypoint("done with KPP_calculate (diabatic)")
    if (cs%debug) then
      call mom_state_chksum("after KPP", u, v, h, g, gv, us, haloshift=0)
      call mom_forcing_chksum("after KPP", fluxes, g, us, haloshift=0)
      call mom_thermovar_chksum("after KPP", tv, g)
      call hchksum(kd_lay, "after KPP Kd_lay", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
      call hchksum(kd_int, "after KPP Kd_Int", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
    endif

  endif  ! endif for KPP

  ! Add vertical diff./visc. due to convection (computed via CVMix)
  if (cs%use_CVMix_conv) then
    call calculate_cvmix_conv(h, tv, g, gv, us, cs%CVMix_conv_csp, hml, kd=kd_int, kv=visc%Kv_slow)
  endif

  if (cs%useKPP) then
    call cpu_clock_begin(id_clock_kpp)
    if (cs%debug) then
      call hchksum(cs%KPP_temp_flux, "before KPP_applyNLT netHeat", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(cs%KPP_salt_flux, "before KPP_applyNLT netSalt", g%HI, haloshift=0, scale=gv%H_to_m)
      call hchksum(cs%KPP_NLTheat, "before KPP_applyNLT NLTheat", g%HI, haloshift=0)
      call hchksum(cs%KPP_NLTscalar, "before KPP_applyNLT NLTscalar", g%HI, haloshift=0)
    endif
    ! Apply non-local transport of heat and salt
    ! Changes: tv%T, tv%S
    call kpp_nonlocaltransport_temp(cs%KPP_CSp, g, gv, h, cs%KPP_NLTheat,   cs%KPP_temp_flux, &
                                    us%T_to_s*dt, tv%T, us%Q_to_J_kg*tv%C_p)
    call kpp_nonlocaltransport_saln(cs%KPP_CSp, g, gv, h, cs%KPP_NLTscalar, cs%KPP_salt_flux, &
                                    us%T_to_s*dt, tv%S)
    call cpu_clock_end(id_clock_kpp)
    if (showcalltree) call calltree_waypoint("done with KPP_applyNonLocalTransport (diabatic)")
    if (cs%debugConservation) call mom_state_stats('KPP_applyNonLocalTransport', u, v, h, tv%T, tv%S, g, gv, us)

    if (cs%debug) then
      call mom_state_chksum("after KPP_applyNLT ", u, v, h, g, gv, us, haloshift=0)
      call mom_forcing_chksum("after KPP_applyNLT ", fluxes, g, us, haloshift=0)
      call mom_thermovar_chksum("after KPP_applyNLT ", tv, g)
    endif
  endif ! endif for KPP

  ! Differential diffusion done here.
  ! Changes: tv%T, tv%S
  if (cs%double_diffuse .and. associated(tv%T)) then

    call cpu_clock_begin(id_clock_differential_diff)
    call differential_diffuse_t_s(h, tv%T, tv%S, kd_extra_t, kd_extra_s, dt, g, gv)
    call cpu_clock_end(id_clock_differential_diff)
    if (showcalltree) call calltree_waypoint("done with differential_diffuse_T_S (diabatic)")
    if (cs%debugConservation) call mom_state_stats('differential_diffuse_T_S', u, v, h, tv%T, tv%S, g, gv, us)

    ! increment heat and salt diffusivity.
    ! CS%useKPP==.true. already has extra_T and extra_S included
    if (.not. cs%useKPP) then
      !$OMP parallel do default(shared)
      do k=2,nz ; do j=js,je ; do i=is,ie
        kd_heat(i,j,k) = kd_heat(i,j,k) + kd_extra_t(i,j,k)
        kd_salt(i,j,k) = kd_salt(i,j,k) + kd_extra_s(i,j,k)
      enddo ; enddo ; enddo
    endif

  endif

  ! Calculate layer entrainments and detrainments from diffusivities and differences between
  ! layer and target densities (i.e. do remapping as well as diffusion).
  call cpu_clock_begin(id_clock_entrain)
  ! Calculate appropriately limited diapycnal mass fluxes to account
  ! for diapycnal diffusion and advection.  Sets: ea, eb. Changes: kb
  call entrainment_diffusive(h, tv, fluxes, dt, g, gv, us, cs%entrain_diffusive_CSp, &
                             ea, eb, kb, kd_lay=kd_lay, kd_int=kd_int)
  call cpu_clock_end(id_clock_entrain)
  if (showcalltree) call calltree_waypoint("done with Entrainment_diffusive (diabatic)")

  if (cs%debug) then
    call mom_forcing_chksum("after calc_entrain ", fluxes, g, us, haloshift=0)
    call mom_thermovar_chksum("after calc_entrain ", tv, g)
    call mom_state_chksum("after calc_entrain ", u, v, h, g, gv, us, haloshift=0)
    call hchksum(ea, "after calc_entrain ea", g%HI, haloshift=0, scale=gv%H_to_m)
    call hchksum(eb, "after calc_entrain eb", g%HI, haloshift=0, scale=gv%H_to_m)
  endif

  ! Save fields before boundary forcing is applied for tendency diagnostics
  if (cs%boundary_forcing_tendency_diag) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      h_diag(i,j,k)    = h(i,j,k)
      temp_diag(i,j,k) = tv%T(i,j,k)
      saln_diag(i,j,k) = tv%S(i,j,k)
    enddo ; enddo ; enddo
  endif

  ! Update h according to divergence of the difference between
  ! ea and eb. We keep a record of the original h in hold.
  ! In the following, the checks for negative values are to guard
  ! against instances where entrainment drives a layer to
  ! negative thickness.  This situation will never happen if
  ! enough iterations are permitted in Calculate_Entrainment.
  ! Even if too few iterations are allowed, it is still guarded
  ! against.  In other words the checks are probably unnecessary.
  !$OMP parallel do default(shared)
  do j=js,je
    do i=is,ie
      hold(i,j,1) = h(i,j,1)
      h(i,j,1) = h(i,j,1) + (eb(i,j,1) - ea(i,j,2))
      hold(i,j,nz) = h(i,j,nz)
      h(i,j,nz) = h(i,j,nz) + (ea(i,j,nz) - eb(i,j,nz-1))
      if (h(i,j,1) <= 0.0) then
        h(i,j,1) = gv%Angstrom_H
      endif
      if (h(i,j,nz) <= 0.0) then
        h(i,j,nz) = gv%Angstrom_H
      endif
    enddo
    do k=2,nz-1 ; do i=is,ie
      hold(i,j,k) = h(i,j,k)
      h(i,j,k) = h(i,j,k) + ((ea(i,j,k) - eb(i,j,k-1)) + &
                    (eb(i,j,k) - ea(i,j,k+1)))
      if (h(i,j,k) <= 0.0) then
        h(i,j,k) = gv%Angstrom_H
      endif
    enddo ; enddo
  enddo
  ! Checks for negative thickness may have changed layer thicknesses
  call diag_update_remap_grids(cs%diag)

  if (cs%debug) then
    call mom_state_chksum("after negative check ", u, v, h, g, gv, us, haloshift=0)
    call mom_forcing_chksum("after negative check ", fluxes, g, us, haloshift=0)
    call mom_thermovar_chksum("after negative check ", tv, g)
  endif
  if (showcalltree) call calltree_waypoint("done with h=ea-eb (diabatic)")
  if (cs%debugConservation) call mom_state_stats('h=ea-eb', u, v, h, tv%T, tv%S, g, gv, us)

  ! Here, T and S are updated according to ea and eb.
  ! If using the bulk mixed layer, T and S are also updated
  ! by surface fluxes (in fluxes%*).
  ! This is a very long block.
  if (cs%bulkmixedlayer) then

    if (associated(tv%T)) then
      call cpu_clock_begin(id_clock_tridiag)
      ! Temperature and salinity (as state variables) are treated
      ! differently from other tracers to insure massless layers that
      ! are lighter than the mixed layer have temperatures and salinities
      ! that correspond to their prescribed densities.
      if (cs%massless_match_targets) then
        !$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) + h_neglect
            b1(i) = 1.0 / (h_tr + eb(i,j,1))
            d1(i) = h_tr * b1(i)
            tv%T(i,j,1) = b1(i) * (h_tr*tv%T(i,j,1))
            tv%S(i,j,1) = b1(i) * (h_tr*tv%S(i,j,1))
          enddo
          do k=2,nkmb ; do i=is,ie
            c1(i,k) = eb(i,j,k-1) * b1(i)
            h_tr = hold(i,j,k) + h_neglect
            b_denom_1 = h_tr + d1(i)*ea(i,j,k)
            b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
            if (k<nkmb) d1(i) = b_denom_1 * b1(i)
            tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
            tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
          enddo ; enddo

          do k=nkmb+1,nz ; do i=is,ie
            if (k == kb(i,j)) then
              c1(i,k) = eb(i,j,k-1) * b1(i)
              d1(i) = (((eb(i,j,nkmb)-eb(i,j,k-1)) + hold(i,j,nkmb) + h_neglect) + &
                       d1(i)*ea(i,j,nkmb)) * b1(i)
              h_tr = hold(i,j,k) + h_neglect
              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)
              tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,nkmb))
              tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,nkmb))
            elseif (k > kb(i,j)) then
              c1(i,k) = eb(i,j,k-1) * b1(i)
              h_tr = hold(i,j,k) + h_neglect
              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)
              tv%T(i,j,k) = b1(i) * (h_tr*tv%T(i,j,k) + ea(i,j,k)*tv%T(i,j,k-1))
              tv%S(i,j,k) = b1(i) * (h_tr*tv%S(i,j,k) + ea(i,j,k)*tv%S(i,j,k-1))
            elseif (eb(i,j,k) < eb(i,j,k-1)) then ! (note that k < kb(i,j))
              !   The bottommost buffer layer might entrain all the mass from some
              ! of the interior layers that are thin and lighter in the coordinate
              ! density than that buffer layer.  The T and S of these newly
              ! massless interior layers are unchanged.
              tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%T(i,j,k)
              tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + b1(i) * (eb(i,j,k-1) - eb(i,j,k)) * tv%S(i,j,k)
            endif
          enddo ; enddo

          do k=nz-1,nkmb,-1 ; do i=is,ie
            if (k >= kb(i,j)) then
              tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
              tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
            endif
          enddo ; enddo
          do i=is,ie ; if (kb(i,j) <= nz) then
            tv%T(i,j,nkmb) = tv%T(i,j,nkmb) + c1(i,kb(i,j))*tv%T(i,j,kb(i,j))
            tv%S(i,j,nkmb) = tv%S(i,j,nkmb) + c1(i,kb(i,j))*tv%S(i,j,kb(i,j))
          endif ; enddo
          do k=nkmb-1,1,-1 ; do i=is,ie
            tv%T(i,j,k) = tv%T(i,j,k) + c1(i,k+1)*tv%T(i,j,k+1)
            tv%S(i,j,k) = tv%S(i,j,k) + c1(i,k+1)*tv%S(i,j,k+1)
          enddo ; enddo
        enddo ! end of j loop
      else ! .not. massless_match_targets
        ! This simpler form allows T & S to be too dense for the layers
        ! between the buffer layers and the interior.
        ! Changes: T, S
        if (cs%tracer_tridiag) then
          call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
          call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
        else
          call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
        endif
      endif ! massless_match_targets
      call cpu_clock_end(id_clock_tridiag)

    endif ! endif for associated(T)
    if (cs%debugConservation) call mom_state_stats('BML tridiag', u, v, h, tv%T, tv%S, g, gv, us)

    if ((cs%ML_mix_first > 0.0) .or. cs%use_geothermal) then
      ! The mixed layer code has already been called, but there is some needed
      ! bookkeeping.
      !$OMP parallel do default(shared)
      do k=1,nz ; do j=js,je ; do i=is,ie
        hold(i,j,k) = h_orig(i,j,k)
        ea(i,j,k) = ea(i,j,k) + eaml(i,j,k)
        eb(i,j,k) = eb(i,j,k) + ebml(i,j,k)
      enddo ; enddo ; enddo
      if (cs%debug) then
        call hchksum(ea, "after ea = ea + eaml", g%HI, haloshift=0, scale=gv%H_to_m)
        call hchksum(eb, "after eb = eb + ebml", g%HI, haloshift=0, scale=gv%H_to_m)
      endif
    endif

    if (cs%ML_mix_first < 1.0) then
    !  Call the mixed layer code now, perhaps for a second time.
    !  This subroutine (1)  Cools the mixed layer.
    !    (2) Performs convective adjustment by mixed layer entrainment.
    !    (3) Heats the mixed layer and causes it to detrain to
    !        Monin-Obukhov depth or minimum mixed layer depth.
    !    (4) Uses any remaining TKE to drive mixed layer entrainment.
    !    (5) Possibly splits the buffer layer into two isopycnal layers.

      call find_uv_at_h(u, v, hold, u_h, v_h, g, gv, us, ea, eb)
      if (cs%debug) call mom_state_chksum("find_uv_at_h1 ", u, v, h, g, gv, us, haloshift=0)

      dt_mix = min(dt, dt*(1.0 - cs%ML_mix_first))
      call cpu_clock_begin(id_clock_mixedlayer)
      ! Changes: h, tv%T, tv%S, ea and eb  (G is also inout???)
      call bulkmixedlayer(h, u_h, v_h, tv, fluxes, dt_mix, ea, eb, &
                          g, gv, us, cs%bulkmixedlayer_CSp, cs%optics, &
                          hml, cs%aggregate_FW_forcing, dt, last_call=.true.)

      !  Keep salinity from falling below a small but positive threshold.
      !  This constraint is needed for SIS1 ice model, which can extract
      !  more salt than is present in the ocean. SIS2 does not suffer
      !  from this limitation, in which case we can let salinity=0 and still
      !  have salt conserved with SIS2 ice. So for SIS2, we can run with
      !  BOUND_SALINITY=False in MOM.F90.
      if (associated(tv%S) .and. associated(tv%salt_deficit)) &
        call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

      call cpu_clock_end(id_clock_mixedlayer)
      if (showcalltree) call calltree_waypoint("done with 2nd bulkmixedlayer (diabatic)")
      if (cs%debugConservation) call mom_state_stats('2nd bulkmixedlayer', u, v, h, tv%T, tv%S, g, gv, us)
    endif

  else  ! following block for when NOT using BULKMIXEDLAYER

    ! calculate change in temperature & salinity due to dia-coordinate surface diffusion
    if (associated(tv%T)) then

      if (cs%debug) then
        call hchksum(ea, "before triDiagTS ea ", g%HI, haloshift=0, scale=gv%H_to_m)
        call hchksum(eb, "before triDiagTS eb ", g%HI, haloshift=0, scale=gv%H_to_m)
      endif
      call cpu_clock_begin(id_clock_tridiag)

      !  Keep salinity from falling below a small but positive threshold.
      !  This constraint is needed for SIS1 ice model, which can extract
      !  more salt than is present in the ocean. SIS2 does not suffer
      !  from this limitation, in which case we can let salinity=0 and still
      !  have salt conserved with SIS2 ice. So for SIS2, we can run with
      !  BOUND_SALINITY=False in MOM.F90.
      if (associated(tv%S) .and. associated(tv%salt_deficit)) &
        call adjust_salt(h, tv, g, gv, cs%diabatic_aux_CSp)

      if (cs%diabatic_diff_tendency_diag) then
        do k=1,nz ; do j=js,je ; do i=is,ie
          temp_diag(i,j,k) = tv%T(i,j,k)
          saln_diag(i,j,k) = tv%S(i,j,k)
        enddo ; enddo ; enddo
      endif

      ! Changes T and S via the tridiagonal solver; no change to h
      if (cs%tracer_tridiag) then
        call tracer_vertdiff(hold, ea, eb, dt, tv%T, g, gv)
        call tracer_vertdiff(hold, ea, eb, dt, tv%S, g, gv)
      else
        call tridiagts(g, gv, is, ie, js, je, hold, ea, eb, tv%T, tv%S)
      endif

      ! diagnose temperature, salinity, heat, and salt tendencies
      ! Note: hold here refers to the thicknesses from before the dual-entraintment when using
      ! the bulk mixed layer scheme, so tendencies should be posted on hold.
      if (cs%diabatic_diff_tendency_diag) then
        call diagnose_diabatic_diff_tendency(tv, hold, temp_diag, saln_diag, dt, g, gv, us, cs)
        if (cs%id_diabatic_diff_h > 0) call post_data(cs%id_diabatic_diff_h, hold, cs%diag, alt_h=hold)
      endif

      call cpu_clock_end(id_clock_tridiag)
      if (showcalltree) call calltree_waypoint("done with triDiagTS (diabatic)")

    endif  ! endif corresponding to if (associated(tv%T))
    if (cs%debugConservation) call mom_state_stats('triDiagTS', u, v, h, tv%T, tv%S, g, gv, us)

  endif  ! endif for the BULKMIXEDLAYER block

  if (cs%debug) then
    call mom_state_chksum("after mixed layer ", u, v, h, g, gv, us, haloshift=0)
    call mom_thermovar_chksum("after mixed layer ", tv, g)
    call hchksum(ea, "after mixed layer ea", g%HI, scale=gv%H_to_m)
    call hchksum(eb, "after mixed layer eb", g%HI, scale=gv%H_to_m)
  endif

  call cpu_clock_begin(id_clock_remap)
  call regularize_layers(h, tv, dt, ea, eb, g, gv, us, cs%regularize_layers_CSp)
  call cpu_clock_end(id_clock_remap)
  if (showcalltree) call calltree_waypoint("done with regularize_layers (diabatic)")
  if (cs%debugConservation) call mom_state_stats('regularize_layers', u, v, h, tv%T, tv%S, g, gv, us)

  ! Whenever thickness changes let the diag manager know, as the
  ! target grids for vertical remapping may need to be regenerated.
  if (associated(adp%du_dt_dia) .or. associated(adp%dv_dt_dia)) &
    ! Remapped d[uv]dt_dia require east/north halo updates of h
    call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
  call diag_update_remap_grids(cs%diag)

  ! diagnostics
  idt = 1.0 / dt
  if ((cs%id_Tdif > 0) .or. (cs%id_Tadv > 0)) then
    do j=js,je ; do i=is,ie
      tdif_flx(i,j,1) = 0.0 ; tdif_flx(i,j,nz+1) = 0.0
      tadv_flx(i,j,1) = 0.0 ; tadv_flx(i,j,nz+1) = 0.0
    enddo ; enddo
    !$OMP parallel do default(shared)
    do k=2,nz ; do j=js,je ; do i=is,ie
      tdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
                        (tv%T(i,j,k-1) - tv%T(i,j,k))
      tadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
                    0.5*(tv%T(i,j,k-1) + tv%T(i,j,k))
    enddo ; enddo ; enddo
  endif
  if ((cs%id_Sdif > 0) .or. (cs%id_Sadv > 0)) then
    do j=js,je ; do i=is,ie
      sdif_flx(i,j,1) = 0.0 ; sdif_flx(i,j,nz+1) = 0.0
      sadv_flx(i,j,1) = 0.0 ; sadv_flx(i,j,nz+1) = 0.0
    enddo ; enddo
    !$OMP parallel do default(shared)
    do k=2,nz ; do j=js,je ; do i=is,ie
      sdif_flx(i,j,k) = (idt * 0.5*(ea(i,j,k) + eb(i,j,k-1))) * &
                        (tv%S(i,j,k-1) - tv%S(i,j,k))
      sadv_flx(i,j,k) = (idt * (ea(i,j,k) - eb(i,j,k-1))) * &
                    0.5*(tv%S(i,j,k-1) + tv%S(i,j,k))
    enddo ; enddo ; enddo
  endif

  ! mixing of passive tracers from massless boundary layers to interior
  call cpu_clock_begin(id_clock_tracers)
  if (cs%mix_boundary_tracers) then
    tr_ea_bbl = gv%Z_to_H * sqrt(dt*cs%Kd_BBL_tr)
    !$OMP parallel do default(shared) private(htot,in_boundary,add_ent)
    do j=js,je
      do i=is,ie
        ebtr(i,j,nz) = eb(i,j,nz)
        htot(i) = 0.0
        in_boundary(i) = (g%mask2dT(i,j) > 0.0)
      enddo
      do k=nz,2,-1 ; do i=is,ie
        if (in_boundary(i)) then
          htot(i) = htot(i) + h(i,j,k)
          !   If diapycnal mixing has been suppressed because this is a massless
          ! layer near the bottom, add some mixing of tracers between these
          ! layers.  This flux is based on the harmonic mean of the two
          ! thicknesses, as this corresponds pretty closely (to within
          ! differences in the density jumps between layers) with what is done
          ! in the calculation of the fluxes in the first place.  Kd_min_tr
          ! should be much less than the values that have been set in Kd_lay,
          ! perhaps a molecular diffusivity.
          add_ent = ((dt * cs%Kd_min_tr) * gv%Z_to_H**2) * &
                    ((h(i,j,k-1)+h(i,j,k)+h_neglect) / &
                     (h(i,j,k-1)*h(i,j,k)+h_neglect2)) - &
                    0.5*(ea(i,j,k) + eb(i,j,k-1))
          if (htot(i) < tr_ea_bbl) then
            add_ent = max(0.0, add_ent, &
                          (tr_ea_bbl - htot(i)) - min(ea(i,j,k), eb(i,j,k-1)))
          elseif (add_ent < 0.0) then
            add_ent = 0.0 ; in_boundary(i) = .false.
          endif

          ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
          eatr(i,j,k) = ea(i,j,k) + add_ent
        else
          ebtr(i,j,k-1) = eb(i,j,k-1) ; eatr(i,j,k) = ea(i,j,k)
        endif
        if (cs%double_diffuse) then ; if (kd_extra_s(i,j,k) > 0.0) then
          add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
             (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
              h_neglect)
          ebtr(i,j,k-1) = ebtr(i,j,k-1) + add_ent
          eatr(i,j,k) = eatr(i,j,k) + add_ent
        endif ; endif
      enddo ; enddo
      do i=is,ie ; eatr(i,j,1) = ea(i,j,1) ; enddo

    enddo

    call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
                              cs%optics, cs%tracer_flow_CSp, cs%debug)

  elseif (cs%double_diffuse) then  ! extra diffusivity for passive tracers

    do j=js,je ; do i=is,ie
      ebtr(i,j,nz) = eb(i,j,nz) ; eatr(i,j,1) = ea(i,j,1)
    enddo ; enddo
    !$OMP parallel do default(shared) private(add_ent)
    do k=nz,2,-1 ; do j=js,je ; do i=is,ie
      if (kd_extra_s(i,j,k) > 0.0) then
        add_ent = ((dt * kd_extra_s(i,j,k)) * gv%Z_to_H**2) / &
           (0.25 * ((h(i,j,k-1) + h(i,j,k)) + (hold(i,j,k-1) + hold(i,j,k))) + &
            h_neglect)
      else
        add_ent = 0.0
      endif
      ebtr(i,j,k-1) = eb(i,j,k-1) + add_ent
      eatr(i,j,k) = ea(i,j,k) + add_ent
    enddo ; enddo ; enddo

    call call_tracer_column_fns(hold, h, eatr, ebtr, fluxes, hml, dt, g, gv, us, tv, &
                                cs%optics, cs%tracer_flow_CSp, cs%debug)

  else
    call call_tracer_column_fns(hold, h, ea, eb, fluxes, hml, dt, g, gv, us, tv, &
                                cs%optics, cs%tracer_flow_CSp, cs%debug)

  endif  ! (CS%mix_boundary_tracers)

  call cpu_clock_end(id_clock_tracers)

  ! sponges
  if (cs%use_sponge) then
    call cpu_clock_begin(id_clock_sponge)
    ! Layer mode sponge
    if (cs%bulkmixedlayer .and. associated(tv%eqn_of_state)) then
      do i=is,ie ; p_ref_cv(i) = tv%P_Ref ; enddo
      eosdom(:) = eos_domain(g%HI)
      !$OMP parallel do default(shared)
      do j=js,je
        call calculate_density(tv%T(:,j,1), tv%S(:,j,1), p_ref_cv, rcv_ml(:,j), &
                               tv%eqn_of_state, eosdom)
      enddo
      call apply_sponge(h, dt, g, gv, us, ea, eb, cs%sponge_CSp, rcv_ml)
    else
      call apply_sponge(h, dt, g, gv, us, ea, eb, cs%sponge_CSp)
    endif
    call cpu_clock_end(id_clock_sponge)
    if (cs%debug) then
      call mom_state_chksum("apply_sponge ", u, v, h, g, gv, us, haloshift=0)
      call mom_thermovar_chksum("apply_sponge ", tv, g)
    endif
  endif ! CS%use_sponge

!   Save the diapycnal mass fluxes as a diagnostic field.
  if (associated(cdp%diapyc_vel)) then
    !$OMP parallel do default(shared)
    do j=js,je
      do k=2,nz ; do i=is,ie
        cdp%diapyc_vel(i,j,k) = us%s_to_T*idt * (ea(i,j,k) - eb(i,j,k-1))
      enddo ; enddo
      do i=is,ie
        cdp%diapyc_vel(i,j,1) = 0.0
        cdp%diapyc_vel(i,j,nz+1) = 0.0
      enddo
    enddo
  endif

! For momentum, it is only the net flux that homogenizes within
! the mixed layer.  Vertical viscosity that is proportional to the
! mixed layer turbulence is applied elsewhere.
  if (cs%bulkmixedlayer) then
    if (cs%debug) then
      call hchksum(ea, "before net flux rearrangement ea", g%HI, scale=gv%H_to_m)
      call hchksum(eb, "before net flux rearrangement eb", g%HI, scale=gv%H_to_m)
    endif
    !$OMP parallel do default(shared) private(net_ent)
    do j=js,je
      do k=2,gv%nkml ; do i=is,ie
        net_ent = ea(i,j,k) - eb(i,j,k-1)
        ea(i,j,k) = max(net_ent, 0.0)
        eb(i,j,k-1) = max(-net_ent, 0.0)
      enddo ; enddo
    enddo
    if (cs%debug) then
      call hchksum(ea, "after net flux rearrangement ea", g%HI, scale=gv%H_to_m)
      call hchksum(eb, "after net flux rearrangement eb", g%HI, scale=gv%H_to_m)
    endif
  endif

! Initialize halo regions of ea, eb, and hold to default values.
  !$OMP parallel do default(shared)
  do k=1,nz
    do i=is-1,ie+1
      hold(i,js-1,k) = gv%Angstrom_H ; ea(i,js-1,k) = 0.0 ; eb(i,js-1,k) = 0.0
      hold(i,je+1,k) = gv%Angstrom_H ; ea(i,je+1,k) = 0.0 ; eb(i,je+1,k) = 0.0
    enddo
    do j=js,je
      hold(is-1,j,k) = gv%Angstrom_H ; ea(is-1,j,k) = 0.0 ; eb(is-1,j,k) = 0.0
      hold(ie+1,j,k) = gv%Angstrom_H ; ea(ie+1,j,k) = 0.0 ; eb(ie+1,j,k) = 0.0
    enddo
  enddo

  call cpu_clock_begin(id_clock_pass)
  if (g%symmetric) then ; dir_flag = to_all+omit_corners
  else ; dir_flag = to_west+to_south+omit_corners ; endif
  call create_group_pass(cs%pass_hold_eb_ea, hold, g%Domain, dir_flag, halo=1)
  call create_group_pass(cs%pass_hold_eb_ea, eb, g%Domain, dir_flag, halo=1)
  call create_group_pass(cs%pass_hold_eb_ea, ea, g%Domain, dir_flag, halo=1)
  call do_group_pass(cs%pass_hold_eb_ea, g%Domain)
  call cpu_clock_end(id_clock_pass)

  !  Use a tridiagonal solver to determine effect of the diapycnal
  !  advection on velocity field. It is assumed that water leaves
  !  or enters the ocean with the surface velocity.
  if (cs%debug) then
    call mom_state_chksum("before u/v tridiag ", u, v, h, g, gv, us, haloshift=0)
    call hchksum(ea, "before u/v tridiag ea", g%HI, scale=gv%H_to_m)
    call hchksum(eb, "before u/v tridiag eb", g%HI, scale=gv%H_to_m)
    call hchksum(hold, "before u/v tridiag hold", g%HI, scale=gv%H_to_m)
  endif
  call cpu_clock_begin(id_clock_tridiag)
  idt_accel = 1.0 / dt
  !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)
  do j=js,je
    do i=isq,ieq
      if (associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,1) = u(i,j,1)
      hval = (hold(i,j,1) + hold(i+1,j,1)) + (ea(i,j,1) + ea(i+1,j,1)) + h_neglect
      b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i+1,j,1)))
      d1(i) = hval * b1(i)
      u(i,j,1) = b1(i) * (hval * u(i,j,1))
    enddo
    do k=2,nz ; do i=isq,ieq
      if (associated(adp%du_dt_dia)) adp%du_dt_dia(i,j,k) = u(i,j,k)
      c1(i,k) = (eb(i,j,k-1)+eb(i+1,j,k-1)) * b1(i)
      eaval = ea(i,j,k) + ea(i+1,j,k)
      hval = hold(i,j,k) + hold(i+1,j,k) + h_neglect
      b1(i) = 1.0 / ((eb(i,j,k) + eb(i+1,j,k)) + (hval + d1(i)*eaval))
      d1(i) = (hval + d1(i)*eaval) * b1(i)
      u(i,j,k) = (hval*u(i,j,k) + eaval*u(i,j,k-1))*b1(i)
    enddo ; enddo
    do k=nz-1,1,-1 ; do i=isq,ieq
      u(i,j,k) = u(i,j,k) + c1(i,k+1)*u(i,j,k+1)
      if (associated(adp%du_dt_dia)) &
        adp%du_dt_dia(i,j,k) = (u(i,j,k) - adp%du_dt_dia(i,j,k)) * idt_accel
    enddo ; enddo
    if (associated(adp%du_dt_dia)) then
      do i=isq,ieq
        adp%du_dt_dia(i,j,nz) = (u(i,j,nz)-adp%du_dt_dia(i,j,nz)) * idt_accel
      enddo
    endif
  enddo
  if (cs%debug) then
    call mom_state_chksum("aft 1st loop tridiag ", u, v, h, g, gv, us, haloshift=0)
  endif
  !$OMP parallel do default(shared) private(hval,b1,d1,c1,eaval)
  do j=jsq,jeq
    do i=is,ie
      if (associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,1) = v(i,j,1)
      hval = (hold(i,j,1) + hold(i,j+1,1)) + (ea(i,j,1) + ea(i,j+1,1)) + h_neglect
      b1(i) = 1.0 / (hval + (eb(i,j,1) + eb(i,j+1,1)))
      d1(i) = hval * b1(i)
      v(i,j,1) = b1(i) * (hval * v(i,j,1))
    enddo
    do k=2,nz ; do i=is,ie
      if (associated(adp%dv_dt_dia)) adp%dv_dt_dia(i,j,k) = v(i,j,k)
      c1(i,k) = (eb(i,j,k-1)+eb(i,j+1,k-1)) * b1(i)
      eaval = ea(i,j,k) + ea(i,j+1,k)
      hval = hold(i,j,k) + hold(i,j+1,k) + h_neglect
      b1(i) = 1.0 / ((eb(i,j,k) + eb(i,j+1,k)) + (hval + d1(i)*eaval))
      d1(i) = (hval + d1(i)*eaval) * b1(i)
      v(i,j,k) = (hval*v(i,j,k) + eaval*v(i,j,k-1))*b1(i)
    enddo ; enddo
    do k=nz-1,1,-1 ; do i=is,ie
      v(i,j,k) = v(i,j,k) + c1(i,k+1)*v(i,j,k+1)
      if (associated(adp%dv_dt_dia)) &
        adp%dv_dt_dia(i,j,k) = (v(i,j,k) - adp%dv_dt_dia(i,j,k)) * idt_accel
    enddo ; enddo
    if (associated(adp%dv_dt_dia)) then
      do i=is,ie
        adp%dv_dt_dia(i,j,nz) = (v(i,j,nz)-adp%dv_dt_dia(i,j,nz)) * idt_accel
      enddo
    endif
  enddo
  call cpu_clock_end(id_clock_tridiag)
  if (cs%debug) then
    call mom_state_chksum("after u/v tridiag ", u, v, h, g, gv, us, haloshift=0)
  endif

  call disable_averaging(cs%diag)
  ! Diagnose the diapycnal diffusivities and other related quantities.
  call enable_averages(dt, time_end, cs%diag)

  if (cs%id_Kd_interface > 0) call post_data(cs%id_Kd_interface, kd_int,  cs%diag)
  if (cs%id_Kd_heat      > 0) call post_data(cs%id_Kd_heat,      kd_heat, cs%diag)
  if (cs%id_Kd_salt      > 0) call post_data(cs%id_Kd_salt,      kd_salt, cs%diag)
  if (cs%id_Kd_ePBL      > 0) call post_data(cs%id_Kd_ePBL,      kd_epbl, cs%diag)

  if (cs%id_ea       > 0) call post_data(cs%id_ea,       ea, cs%diag)
  if (cs%id_eb       > 0) call post_data(cs%id_eb,       eb, cs%diag)

  if (cs%id_dudt_dia > 0) call post_data(cs%id_dudt_dia, adp%du_dt_dia,  cs%diag)
  if (cs%id_dvdt_dia > 0) call post_data(cs%id_dvdt_dia, adp%dv_dt_dia,  cs%diag)
  if (cs%id_wd       > 0) call post_data(cs%id_wd,       cdp%diapyc_vel, cs%diag)

  if (cs%id_Tdif > 0) call post_data(cs%id_Tdif, tdif_flx, cs%diag)
  if (cs%id_Tadv > 0) call post_data(cs%id_Tadv, tadv_flx, cs%diag)
  if (cs%id_Sdif > 0) call post_data(cs%id_Sdif, sdif_flx, cs%diag)
  if (cs%id_Sadv > 0) call post_data(cs%id_Sadv, sadv_flx, cs%diag)

  !! Diagnostics for terms multiplied by fractional thicknesses
  if (cs%id_hf_dudt_dia_2d > 0) then
    allocate(hf_dudt_dia_2d(g%IsdB:g%IedB,g%jsd:g%jed))
    hf_dudt_dia_2d(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=isq,ieq
      hf_dudt_dia_2d(i,j) = hf_dudt_dia_2d(i,j) + adp%du_dt_dia(i,j,k) * adp%diag_hfrac_u(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_hf_dudt_dia_2d, hf_dudt_dia_2d, cs%diag)
    deallocate(hf_dudt_dia_2d)
  endif

  if (cs%id_hf_dvdt_dia_2d > 0) then
    allocate(hf_dvdt_dia_2d(g%isd:g%ied,g%JsdB:g%JedB))
    hf_dvdt_dia_2d(:,:) = 0.0
    do k=1,nz ; do j=jsq,jeq ; do i=is,ie
      hf_dvdt_dia_2d(i,j) = hf_dvdt_dia_2d(i,j) + adp%dv_dt_dia(i,j,k) * adp%diag_hfrac_v(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_hf_dvdt_dia_2d, hf_dvdt_dia_2d, cs%diag)
    deallocate(hf_dvdt_dia_2d)
  endif

  call disable_averaging(cs%diag)

  if (showcalltree) call calltree_leave("layered_diabatic()")

end subroutine layered_diabatic

!> Returns pointers or values of members within the diabatic_CS type. For extensibility,
!! each returned argument is an optional argument
subroutine extract_diabatic_member(CS, opacity_CSp, optics_CSp, evap_CFL_limit, minimum_forcing_depth, &
                                   KPP_CSp, energetic_PBL_CSp, diabatic_aux_CSp, diabatic_halo)
  type(diabatic_cs), intent(in   )           :: CS !< module control structure
  ! All output arguments are optional
  type(opacity_cs),  optional, pointer       :: opacity_CSp !< A pointer to be set to the opacity control structure
  type(optics_type), optional, pointer       :: optics_CSp  !< A pointer to be set to the optics control structure
  type(kpp_cs),      optional, pointer       :: KPP_CSp     !< A pointer to be set to the KPP CS
  type(energetic_pbl_cs), optional, pointer  :: energetic_PBL_CSp !< A pointer to be set to the ePBL CS
  real,              optional, intent(  out) :: evap_CFL_limit !<The largest fraction of a layer that can be
                                                            !! evaporated in one time-step [nondim].
  real,              optional, intent(  out) :: minimum_forcing_depth !< The smallest depth over which heat
                                                            !! and freshwater fluxes are applied [H ~> m or kg m-2].
  type(diabatic_aux_cs), optional, pointer   :: diabatic_aux_CSp !< A pointer to be set to the diabatic_aux
                                                            !! control structure
  integer,           optional, intent(  out) :: diabatic_halo !< The halo size where the diabatic algorithms
                                                            !! assume thermodynamics properties are valid.

  ! Pointers to control structures
  if (present(opacity_csp))       opacity_csp => cs%opacity_CSp
  if (present(optics_csp))        optics_csp  => cs%optics
  if (present(kpp_csp))           kpp_csp     => cs%KPP_CSp
  if (present(energetic_pbl_csp)) energetic_pbl_csp => cs%energetic_PBL_CSp

  ! Constants within diabatic_CS
  if (present(evap_cfl_limit))        evap_cfl_limit = cs%evap_CFL_limit
  if (present(minimum_forcing_depth)) minimum_forcing_depth = cs%minimum_forcing_depth
  if (present(diabatic_halo)) diabatic_halo = cs%halo_TS_diff

end subroutine extract_diabatic_member

!> Routine called for adiabatic physics
subroutine adiabatic(h, tv, fluxes, dt, G, GV, US, CS)
  type(ocean_grid_type),   intent(inout) :: G      !< ocean grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(inout) :: h      !< thickness [H ~> m or kg m-2]
  type(thermo_var_ptrs),   intent(inout) :: tv     !< points to thermodynamic fields
  type(forcing),           intent(inout) :: fluxes !< boundary fluxes
  real,                    intent(in)    :: dt     !< time step [T ~> s]
  type(verticalgrid_type), intent(in)    :: GV     !< ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: US     !< A dimensional unit scaling type
  type(diabatic_cs),       pointer       :: CS     !< module control structure

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: zeros  ! An array of zeros.

  zeros(:,:,:) = 0.0

  call call_tracer_column_fns(h, h, zeros, zeros, fluxes, zeros(:,:,1), dt, g, gv, us, tv, &
                              cs%optics, cs%tracer_flow_CSp, cs%debug)

end subroutine adiabatic


!> This routine diagnoses tendencies from application of diabatic diffusion
!! using ALE algorithm. Note that layer thickness is not altered by
!! diabatic diffusion.
subroutine diagnose_diabatic_diff_tendency(tv, h, temp_old, saln_old, dt, G, GV, US, CS)
  type(ocean_grid_type),                     intent(in) :: G        !< ocean grid structure
  type(verticalgrid_type),                   intent(in) :: GV       !< ocean vertical grid structure
  type(thermo_var_ptrs),                     intent(in) :: tv       !< points to updated thermodynamic fields
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h        !< thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: temp_old !< temperature prior to diabatic physics
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: saln_old !< salinity prior to diabatic physics [ppt]
  real,                                      intent(in) :: dt       !< time step [T ~> s]
  type(unit_scale_type),                     intent(in) :: US       !< A dimensional unit scaling type
  type(diabatic_cs),                         pointer    :: CS       !< module control structure

  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
  real, dimension(SZI_(G),SZJ_(G))         :: work_2d
  real :: Idt  ! The inverse of the timestep [T-1 ~> s-1]
  real :: ppt2mks = 0.001  ! Conversion factor from g/kg to kg/kg.
  integer :: i, j, k, is, ie, js, je, nz
  logical :: do_saln_tend   ! Calculate salinity-based tendency diagnosics

  is  = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  idt = 0.0 ; if (dt > 0.0) idt = 1. / dt
  work_3d(:,:,:) = 0.0
  work_2d(:,:)   = 0.0


  ! temperature tendency
  do k=1,nz ; do j=js,je ; do i=is,ie
    work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt
  enddo ; enddo ; enddo
  if (cs%id_diabatic_diff_temp_tend > 0) then
    call post_data(cs%id_diabatic_diff_temp_tend, work_3d, cs%diag, alt_h=h)
  endif

  ! heat tendency
  if (cs%id_diabatic_diff_heat_tend > 0 .or. cs%id_diabatic_diff_heat_tend_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = h(i,j,k)*gv%H_to_RZ * tv%C_p * work_3d(i,j,k)
    enddo ; enddo ; enddo
    if (cs%id_diabatic_diff_heat_tend > 0) then
      call post_data(cs%id_diabatic_diff_heat_tend, work_3d, cs%diag, alt_h=h)
    endif
    if (cs%id_diabatic_diff_heat_tend_2d > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = 0.0
        do k=1,nz
          work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_diabatic_diff_heat_tend_2d, work_2d, cs%diag)
    endif
  endif

  ! salinity tendency
  do_saln_tend = cs%id_diabatic_diff_saln_tend > 0 &
    .or. cs%id_diabatic_diff_salt_tend > 0 &
    .or. cs%id_diabatic_diff_salt_tend_2d > 0

  if (do_saln_tend) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = (tv%S(i,j,k) - saln_old(i,j,k)) * idt
    enddo ; enddo ; enddo

    if (cs%id_diabatic_diff_saln_tend > 0) &
      call post_data(cs%id_diabatic_diff_saln_tend, work_3d, cs%diag, alt_h=h)

    ! salt tendency
    if (cs%id_diabatic_diff_salt_tend > 0 .or. cs%id_diabatic_diff_salt_tend_2d > 0) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        work_3d(i,j,k) = h(i,j,k)*gv%H_to_RZ * ppt2mks * work_3d(i,j,k)
      enddo ; enddo ; enddo
      if (cs%id_diabatic_diff_salt_tend > 0) then
        call post_data(cs%id_diabatic_diff_salt_tend, work_3d, cs%diag, alt_h=h)
      endif
      if (cs%id_diabatic_diff_salt_tend_2d > 0) then
        do j=js,je ; do i=is,ie
          work_2d(i,j) = 0.0
          do k=1,nz
            work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
          enddo
        enddo ; enddo
        call post_data(cs%id_diabatic_diff_salt_tend_2d, work_2d, cs%diag)
      endif
    endif
  endif

end subroutine diagnose_diabatic_diff_tendency


!> This routine diagnoses tendencies from application of boundary fluxes.
!! These impacts are generally 3d, in particular for penetrative shortwave.
!! Other fluxes contribute 3d in cases when the layers vanish or are very thin,
!! in which case we distribute the flux into k > 1 layers.
subroutine diagnose_boundary_forcing_tendency(tv, h, temp_old, saln_old, h_old, &
                                              dt, G, GV, US, CS)
  type(ocean_grid_type),   intent(in) :: G        !< ocean grid structure
  type(verticalgrid_type), intent(in) :: GV       !< ocean vertical grid structure
  type(thermo_var_ptrs),   intent(in) :: tv       !< points to updated thermodynamic fields
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: h        !< thickness after boundary flux application [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: temp_old !< temperature prior to boundary flux application [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: saln_old !< salinity prior to boundary flux application [ppt]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: h_old    !< thickness prior to boundary flux application [H ~> m or kg m-2]
  real,                    intent(in) :: dt       !< time step [T ~> s]
  type(unit_scale_type),   intent(in) :: US       !< A dimensional unit scaling type
  type(diabatic_cs),       pointer    :: CS       !< module control structure

  ! Local variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: work_3d
  real, dimension(SZI_(G),SZJ_(G))         :: work_2d
  real :: Idt  ! The inverse of the timestep [T-1 ~> s-1]
  real :: ppt2mks = 0.001  ! Conversion factor from g/kg to kg/kg.
  integer :: i, j, k, is, ie, js, je, nz

  is  = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  idt = 0.0 ; if (dt > 0.0) idt = 1. / dt
  work_3d(:,:,:) = 0.0
  work_2d(:,:)   = 0.0

  ! Thickness tendency
  if (cs%id_boundary_forcing_h_tendency > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = (h(i,j,k) - h_old(i,j,k))*idt
    enddo ; enddo ; enddo
    call post_data(cs%id_boundary_forcing_h_tendency, work_3d, cs%diag, alt_h=h_old)
  endif

  ! temperature tendency
  if (cs%id_boundary_forcing_temp_tend > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = (tv%T(i,j,k)-temp_old(i,j,k))*idt
    enddo ; enddo ; enddo
    call post_data(cs%id_boundary_forcing_temp_tend, work_3d, cs%diag, alt_h=h_old)
  endif

  ! heat tendency
  if (cs%id_boundary_forcing_heat_tend > 0 .or. cs%id_boundary_forcing_heat_tend_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = gv%H_to_RZ * tv%C_p * idt * (h(i,j,k) * tv%T(i,j,k) - h_old(i,j,k) * temp_old(i,j,k))
    enddo ; enddo ; enddo
    if (cs%id_boundary_forcing_heat_tend > 0) then
      call post_data(cs%id_boundary_forcing_heat_tend, work_3d, cs%diag, alt_h = h_old)
    endif
    if (cs%id_boundary_forcing_heat_tend_2d > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = 0.0
        do k=1,nz
          work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_boundary_forcing_heat_tend_2d, work_2d, cs%diag)
    endif
  endif

  ! salinity tendency
  if (cs%id_boundary_forcing_saln_tend > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = (tv%S(i,j,k)-saln_old(i,j,k))*idt
    enddo ; enddo ; enddo
    call post_data(cs%id_boundary_forcing_saln_tend, work_3d, cs%diag, alt_h = h_old)
  endif

  ! salt tendency
  if (cs%id_boundary_forcing_salt_tend > 0 .or. cs%id_boundary_forcing_salt_tend_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_3d(i,j,k) = gv%H_to_RZ * ppt2mks * idt * (h(i,j,k) * tv%S(i,j,k) - h_old(i,j,k) * saln_old(i,j,k))
    enddo ; enddo ; enddo
    if (cs%id_boundary_forcing_salt_tend > 0) then
      call post_data(cs%id_boundary_forcing_salt_tend, work_3d, cs%diag, alt_h = h_old)
    endif
    if (cs%id_boundary_forcing_salt_tend_2d > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = 0.0
        do k=1,nz
          work_2d(i,j) = work_2d(i,j) + work_3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_boundary_forcing_salt_tend_2d, work_2d, cs%diag)
    endif
  endif

end subroutine diagnose_boundary_forcing_tendency


!> This routine diagnoses tendencies for temperature and heat from frazil formation.
!! This routine is called twice from within subroutine diabatic; at start and at
!! end of the diabatic processes. The impacts from frazil are generally a function
!! of depth.  Hence, when checking heat budget, be sure to remove HFSIFRAZIL from HFDS in k=1.
subroutine diagnose_frazil_tendency(tv, h, temp_old, dt, G, GV, US, CS)
  type(ocean_grid_type),                    intent(in) :: G        !< ocean grid structure
  type(verticalgrid_type),                  intent(in) :: GV       !< ocean vertical grid structure
  type(thermo_var_ptrs),                    intent(in) :: tv       !< points to updated thermodynamic fields
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h        !< thickness [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: temp_old !< temperature prior to frazil formation [degC]
  real,                                     intent(in) :: dt       !< time step [T ~> s]
  type(unit_scale_type),                    intent(in) :: US       !< A dimensional unit scaling type
  type(diabatic_cs),                        pointer    :: CS       !< module control structure

  real, dimension(SZI_(G),SZJ_(G))         :: work_2d
  real    :: Idt ! The inverse of the timestep [T-1 ~> s-1]
  integer :: i, j, k, is, ie, js, je, nz

  is  = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  idt = 0.0 ; if (dt > 0.0) idt = 1. / dt

  ! temperature tendency
  if (cs%id_frazil_temp_tend > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%frazil_temp_diag(i,j,k) = idt * (tv%T(i,j,k)-temp_old(i,j,k))
    enddo ; enddo ; enddo
    call post_data(cs%id_frazil_temp_tend, cs%frazil_temp_diag(:,:,:), cs%diag)
  endif

  ! heat tendency
  if (cs%id_frazil_heat_tend > 0 .or. cs%id_frazil_heat_tend_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%frazil_heat_diag(i,j,k) = gv%H_to_RZ * tv%C_p * h(i,j,k) * idt * (tv%T(i,j,k)-temp_old(i,j,k))
    enddo ; enddo ; enddo
    if (cs%id_frazil_heat_tend > 0) call post_data(cs%id_frazil_heat_tend, cs%frazil_heat_diag(:,:,:), cs%diag)

    ! As a consistency check, we must have
    ! FRAZIL_HEAT_TENDENCY_2d = HFSIFRAZIL
    if (cs%id_frazil_heat_tend_2d > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = 0.0
        do k=1,nz
          work_2d(i,j) = work_2d(i,j) + cs%frazil_heat_diag(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_frazil_heat_tend_2d, work_2d, cs%diag)
    endif
  endif

end subroutine diagnose_frazil_tendency


!> A simplified version of diabatic_driver_init that will allow
!! tracer column functions to be called without allowing any
!! of the diabatic processes to be used.
subroutine adiabatic_driver_init(Time, G, param_file, diag, CS, &
                                tracer_flow_CSp)
  type(time_type),         intent(in)    :: Time             !< current model time
  type(ocean_grid_type),   intent(in)    :: G                !< model grid structure
  type(param_file_type),   intent(in)    :: param_file       !< the file to parse for parameter values
  type(diag_ctrl), target, intent(inout) :: diag             !< regulates diagnostic output
  type(diabatic_cs),       pointer       :: CS               !< module control structure
  type(tracer_flow_control_cs), pointer  :: tracer_flow_CSp  !< pointer to control structure of the
                                                             !! tracer flow control module

! This "include" declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_diabatic_driver" ! This module's name.

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

  cs%diag => diag
  if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp

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

end subroutine adiabatic_driver_init


!> This routine initializes the diabatic driver module.
subroutine diabatic_driver_init(Time, G, GV, US, param_file, useALEalgorithm, diag, &
                                ADp, CDp, CS, tracer_flow_CSp, sponge_CSp, &
                                ALE_sponge_CSp)
  type(time_type), target                :: Time             !< model time
  type(ocean_grid_type),   intent(inout) :: G                !< model grid structure
  type(verticalgrid_type), intent(in)    :: GV               !< model vertical grid structure
  type(unit_scale_type),   intent(in)    :: US               !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file       !< file to parse for parameter values
  logical,                 intent(in)    :: useALEalgorithm  !< logical for whether to use ALE remapping
  type(diag_ctrl), target, intent(inout) :: diag             !< structure to regulate diagnostic output
  type(accel_diag_ptrs),   intent(inout) :: ADp              !< pointers to accelerations in momentum equations,
                                                             !! to enable diagnostics, like energy budgets
  type(cont_diag_ptrs),    intent(inout) :: CDp              !< pointers to terms in continuity equations
  type(diabatic_cs),       pointer       :: CS               !< module control structure
  type(tracer_flow_control_cs), pointer  :: tracer_flow_CSp  !< pointer to control structure of the
                                                             !! tracer flow control module
  type(sponge_cs),         pointer       :: sponge_CSp       !< pointer to the sponge module control structure
  type(ale_sponge_cs),     pointer       :: ALE_sponge_CSp   !< pointer to the ALE sponge module control structure

  real    :: Kd  ! A diffusivity used in the default for other tracer diffusivities, in MKS units [m2 s-1]
  integer :: num_mode
  logical :: use_temperature
  character(len=20) :: EN1, EN2, EN3

! This "include" declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_diabatic_driver" ! This module's name.
  character(len=48)  :: thickness_units
  character(len=40)  :: var_name
  character(len=160) :: var_descript
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nbands, m
  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_driver_init called with an "// &
                            "associated control structure.")
    return
  else
    allocate(cs)
  endif

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

  if (associated(tracer_flow_csp)) cs%tracer_flow_CSp => tracer_flow_csp
  if (associated(sponge_csp))      cs%sponge_CSp      => sponge_csp
  if (associated(ale_sponge_csp))  cs%ALE_sponge_CSp  => ale_sponge_csp

  cs%useALEalgorithm = usealealgorithm
  cs%bulkmixedlayer = (gv%nkml > 0)

  ! Set default, read and log parameters
  call log_version(param_file, mdl, version, &
                   "The following parameters are used for diabatic processes.", &
                   log_to_all=.true., debugging=.true.)
  call get_param(param_file, mdl, "USE_LEGACY_DIABATIC_DRIVER", cs%use_legacy_diabatic, &
                 "If true, use a legacy version of the diabatic subroutine. "//&
                 "This is temporary and is needed to avoid change in answers.", &
                 default=.true.)
  call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &
                 "If true, sponges may be applied anywhere in the domain. "//&
                 "The exact location and properties of those sponges are "//&
                 "specified via calls to initialize_sponge and possibly "//&
                 "set_up_sponge_field.", default=.false.)
  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, "ENERGETICS_SFC_PBL", cs%use_energetic_PBL, &
                 "If true, use an implied energetics planetary boundary "//&
                 "layer scheme to determine the diffusivity and viscosity "//&
                 "in the surface boundary layer.", default=.false.)
  call get_param(param_file, mdl, "EPBL_IS_ADDITIVE", cs%ePBL_is_additive, &
                 "If true, the diffusivity from ePBL is added to all "//&
                 "other diffusivities. Otherwise, the larger of kappa-shear "//&
                 "and ePBL diffusivities are used.", default=.true.)
  call get_param(param_file, mdl, "PRANDTL_EPBL", cs%ePBL_Prandtl, &
                 "The Prandtl number used by ePBL to convert vertical diffusivities into "//&
                 "viscosities.", default=1.0, units="nondim", do_not_log=.not.cs%use_energetic_PBL)
  call get_param(param_file, mdl, "USE_KPP", cs%use_KPP, &
                 "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
                 "to calculate diffusivities and non-local transport in the OBL.", &
                 default=.false., do_not_log=.true.)
  cs%use_CVMix_ddiff = cvmix_ddiff_is_used(param_file)

  cs%use_kappa_shear = kappa_shear_is_used(param_file)
  cs%use_CVMix_shear = cvmix_shear_is_used(param_file)

  if (cs%bulkmixedlayer) then
    call get_param(param_file, mdl, "ML_MIX_FIRST", cs%ML_mix_first, &
                 "The fraction of the mixed layer mixing that is applied "//&
                 "before interior diapycnal mixing.  0 by default.", &
                 units="nondim", default=0.0)
    call get_param(param_file, mdl, "NKBL", cs%nkbl, default=2, do_not_log=.true.)
  else
    cs%ML_mix_first = 0.0
  endif
  if (use_temperature) then
    call get_param(param_file, mdl, "DO_GEOTHERMAL", cs%use_geothermal, &
                 "If true, apply geothermal heating.", default=.false.)
  else
    cs%use_geothermal = .false.
  endif
  call get_param(param_file, mdl, "INTERNAL_TIDES", cs%use_int_tides, &
                 "If true, use the code that advances a separate set of "//&
                 "equations for the internal tide energy density.", default=.false.)
  cs%nMode = 1
  if (cs%use_int_tides) then
    call get_param(param_file, mdl, "INTERNAL_TIDE_MODES", cs%nMode, &
                 "The number of distinct internal tide modes "//&
                 "that will be calculated.", default=1, do_not_log=.true.)
    call get_param(param_file, mdl, "UNIFORM_TEST_CG", cs%uniform_test_cg, &
                 "If positive, a uniform group velocity of internal tide for test case", &
                 default=-1., units="m s-1", scale=us%m_s_to_L_T)
  endif

  call get_param(param_file, mdl, "MASSLESS_MATCH_TARGETS", &
                                cs%massless_match_targets, &
                 "If true, the temperature and salinity of massless layers "//&
                 "are kept consistent with their target densities. "//&
                 "Otherwise the properties of massless layers evolve "//&
                 "diffusively to match massive neighboring layers.", &
                 default=.true.)

  call get_param(param_file, mdl, "AGGREGATE_FW_FORCING", cs%aggregate_FW_forcing, &
                 "If true, the net incoming and outgoing fresh water fluxes are combined "//&
                 "and applied as either incoming or outgoing depending on the sign of the net. "//&
                 "If false, the net incoming fresh water flux is added to the model and "//&
                 "thereafter the net outgoing is removed from the topmost non-vanished "//&
                 "layers of the updated state.", default=.true.)

  call get_param(param_file, mdl, "DEBUG", cs%debug, &
                 "If true, write out verbose debugging data.", &
                 default=.false., debuggingparam=.true.)
  call get_param(param_file, mdl, "DEBUG_CONSERVATION", cs%debugConservation, &
                 "If true, monitor conservation and extrema.", &
                 default=.false., debuggingparam=.true.)

  call get_param(param_file, mdl, "DEBUG_ENERGY_REQ", cs%debug_energy_req, &
                 "If true, debug the energy requirements.", default=.false., do_not_log=.true.)
  call get_param(param_file, mdl, "MIX_BOUNDARY_TRACERS", cs%mix_boundary_tracers, &
                 "If true, mix the passive tracers in massless layers at "//&
                 "the bottom into the interior as though a diffusivity of "//&
                 "KD_MIN_TR were operating.", default=.true.)

  if (cs%mix_boundary_tracers) then
    call get_param(param_file, mdl, "KD", kd, default=0.0)
    call get_param(param_file, mdl, "KD_MIN_TR", cs%Kd_min_tr, &
                 "A minimal diffusivity that should always be applied to "//&
                 "tracers, especially in massless layers near the bottom. "//&
                 "The default is 0.1*KD.", units="m2 s-1", default=0.1*kd, scale=us%m2_s_to_Z2_T)
    call get_param(param_file, mdl, "KD_BBL_TR", cs%Kd_BBL_tr, &
                 "A bottom boundary layer tracer diffusivity that will "//&
                 "allow for explicitly specified bottom fluxes. The "//&
                 "entrainment at the bottom is at least sqrt(Kd_BBL_tr*dt) "//&
                 "over the same distance.", units="m2 s-1", default=0., scale=us%m2_s_to_Z2_T)
  endif

  call get_param(param_file, mdl, "TRACER_TRIDIAG", cs%tracer_tridiag, &
                 "If true, use the passive tracer tridiagonal solver for T and S", &
                 default=.false.)

  call get_param(param_file, mdl, "MINIMUM_FORCING_DEPTH", cs%minimum_forcing_depth, &
                 "The smallest depth over which forcing can be applied. This "//&
                 "only takes effect when near-surface layers become thin "//&
                 "relative to this scale, in which case the forcing tendencies "//&
                 "scaled down by distributing the forcing over this depth scale.", &
                 units="m", default=0.001, scale=gv%m_to_H)
  call get_param(param_file, mdl, "EVAP_CFL_LIMIT", cs%evap_CFL_limit, &
                 "The largest fraction of a layer than can be lost to forcing "//&
                 "(e.g. evaporation, sea-ice formation) in one time-step. The unused "//&
                 "mass loss is passed down through the column.", &
                 units="nondim", default=0.8)


  ! Register all available diagnostics for this module.
  thickness_units = get_thickness_units(gv)

  cs%id_ea_t = register_diag_field('ocean_model', 'ea_t', diag%axesTL, time, &
      'Layer (heat) entrainment from above per timestep', 'm', conversion=gv%H_to_m)
  cs%id_eb_t = register_diag_field('ocean_model', 'eb_t', diag%axesTL, time, &
      'Layer (heat) entrainment from below per timestep', 'm', conversion=gv%H_to_m)
  cs%id_ea_s = register_diag_field('ocean_model', 'ea_s', diag%axesTL, time, &
      'Layer (salt) entrainment from above per timestep', 'm', conversion=gv%H_to_m)
  cs%id_eb_s = register_diag_field('ocean_model', 'eb_s', diag%axesTL, time, &
      'Layer (salt) entrainment from below per timestep', 'm', conversion=gv%H_to_m)
  ! used by layer diabatic
  cs%id_ea = register_diag_field('ocean_model', 'ea', diag%axesTL, time, &
      'Layer entrainment from above per timestep', 'm', conversion=gv%H_to_m)
  cs%id_eb = register_diag_field('ocean_model', 'eb', diag%axesTL, time, &
      'Layer entrainment from below per timestep', 'm', conversion=gv%H_to_m)
  cs%id_wd = register_diag_field('ocean_model', 'wd', diag%axesTi, time, &
    'Diapycnal velocity', 'm s-1', conversion=gv%H_to_m)
  if (cs%id_wd > 0) call safe_alloc_ptr(cdp%diapyc_vel,isd,ied,jsd,jed,nz+1)

  cs%id_dudt_dia = register_diag_field('ocean_model', 'dudt_dia', diag%axesCuL, time, &
      'Zonal Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)
  cs%id_dvdt_dia = register_diag_field('ocean_model', 'dvdt_dia', diag%axesCvL, time, &
      'Meridional Acceleration from Diapycnal Mixing', 'm s-2', conversion=us%L_T2_to_m_s2)

  cs%id_hf_dudt_dia_2d = register_diag_field('ocean_model', 'hf_dudt_dia_2d', diag%axesCu1, time, &
      'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Diapycnal Mixing', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
  if (cs%id_hf_dudt_dia_2d > 0) then
    call safe_alloc_ptr(adp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
  endif

  cs%id_hf_dvdt_dia_2d = register_diag_field('ocean_model', 'hf_dvdt_dia_2d', diag%axesCv1, time, &
      'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Diapycnal Mixing', 'm s-2', &
      conversion=us%L_T2_to_m_s2)
  if (cs%id_hf_dvdt_dia_2d > 0) then
    call safe_alloc_ptr(adp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
    call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
  endif

  if (cs%use_int_tides) then
    cs%id_cg1 = register_diag_field('ocean_model', 'cn1', diag%axesT1, &
                 time, 'First baroclinic mode (eigen) speed', 'm s-1', conversion=us%L_T_to_m_s)
    allocate(cs%id_cn(cs%nMode)) ; cs%id_cn(:) = -1
    do m=1,cs%nMode
      write(var_name, '("cn_mode",i1)') m
      write(var_descript, '("Baroclinic (eigen) speed of mode ",i1)') m
      cs%id_cn(m) = register_diag_field('ocean_model',var_name, diag%axesT1, &
                   time, var_descript, 'm s-1', conversion=us%L_T_to_m_s)
      call mom_mesg("Registering "//trim(var_name)//", Described as: "//var_descript, 5)
    enddo
  endif

  if (use_temperature) then
    cs%id_Tdif = register_diag_field('ocean_model',"Tflx_dia_diff", diag%axesTi, &
        time, "Diffusive diapycnal temperature flux across interfaces", &
        "degC m s-1", conversion=gv%H_to_m*us%s_to_T)
    cs%id_Tadv = register_diag_field('ocean_model',"Tflx_dia_adv", diag%axesTi, &
        time, "Advective diapycnal temperature flux across interfaces", &
        "degC m s-1", conversion=gv%H_to_m*us%s_to_T)
    cs%id_Sdif = register_diag_field('ocean_model',"Sflx_dia_diff", diag%axesTi, &
        time, "Diffusive diapycnal salnity flux across interfaces", &
        "psu m s-1", conversion=gv%H_to_m*us%s_to_T)
    cs%id_Sadv = register_diag_field('ocean_model',"Sflx_dia_adv", diag%axesTi, &
        time, "Advective diapycnal salnity flux across interfaces", &
        "psu m s-1", conversion=gv%H_to_m*us%s_to_T)
    cs%id_MLD_003 = register_diag_field('ocean_model', 'MLD_003', diag%axesT1, time, &
        'Mixed layer depth (delta rho = 0.03)', 'm', conversion=us%Z_to_m, &
        cmor_field_name='mlotst', cmor_long_name='Ocean Mixed Layer Thickness Defined by Sigma T', &
        cmor_standard_name='ocean_mixed_layer_thickness_defined_by_sigma_t')
    cs%id_mlotstsq = register_diag_field('ocean_model', 'mlotstsq', diag%axesT1, time, &
        long_name='Square of Ocean Mixed Layer Thickness Defined by Sigma T', &
        standard_name='square_of_ocean_mixed_layer_thickness_defined_by_sigma_t', &
        units='m2', conversion=us%Z_to_m**2)
    cs%id_MLD_0125 = register_diag_field('ocean_model', 'MLD_0125', diag%axesT1, time, &
        'Mixed layer depth (delta rho = 0.125)', 'm', conversion=us%Z_to_m)
    call get_param(param_file, mdl, "MLD_EN_VALS", cs%MLD_EN_VALS, &
         "The energy values used to compute MLDs.  If not set (or all set to 0.), the "//&
         "default will overwrite to 25., 2500., 250000.",units='J/m2', default=0., &
         scale=us%kg_m3_to_R*us%m_to_Z**3*us%T_to_s**2)
    if ((cs%MLD_EN_VALS(1)==0.).and.(cs%MLD_EN_VALS(2)==0.).and.(cs%MLD_EN_VALS(3)==0.)) then
      cs%MLD_EN_VALS = (/25.*us%kg_m3_to_R*us%m_to_Z*us%m_to_L**2*us%T_to_s**2,&
           2500.*us%kg_m3_to_R*us%m_to_Z*us%m_to_L**2*us%T_to_s**2,&
           250000.*us%kg_m3_to_R*us%m_to_Z*us%m_to_L**2*us%T_to_s**2/)
    endif
    write(en1,'(F10.2)') cs%MLD_EN_VALS(1)*us%R_to_kg_m3*us%Z_to_m*us%L_to_m**2*us%s_to_T**2
    write(en2,'(F10.2)') cs%MLD_EN_VALS(2)*us%R_to_kg_m3*us%Z_to_m*us%L_to_m**2*us%s_to_T**2
    write(en3,'(F10.2)') cs%MLD_EN_VALS(3)*us%R_to_kg_m3*us%Z_to_m*us%L_to_m**2*us%s_to_T**2
    cs%id_MLD_EN1 = register_diag_field('ocean_model', 'MLD_EN1', diag%axesT1, time, &
         'Mixed layer depth for energy value set to '//trim(en1)//' J/m2 (Energy set by 1st MLD_EN_VALS)', &
         'm', conversion=us%Z_to_m)
    cs%id_MLD_EN2 = register_diag_field('ocean_model', 'MLD_EN2', diag%axesT1, time, &
         'Mixed layer depth for energy value set to '//trim(en2)//' J/m2 (Energy set by 2nd MLD_EN_VALS)', &
         'm', conversion=us%Z_to_m)
    cs%id_MLD_EN3 = register_diag_field('ocean_model', 'MLD_EN3', diag%axesT1, time, &
         'Mixed layer depth for energy value set to '//trim(en3)//' J/m2 (Energy set by 3rd MLD_EN_VALS)', &
         'm', conversion=us%Z_to_m)
    cs%id_subMLN2  = register_diag_field('ocean_model', 'subML_N2', diag%axesT1, time, &
        'Squared buoyancy frequency below mixed layer', 's-2', conversion=us%s_to_T**2)
    cs%id_MLD_user = register_diag_field('ocean_model', 'MLD_user', diag%axesT1, time, &
        'Mixed layer depth (used defined)', 'm', conversion=us%Z_to_m)
  endif
  call get_param(param_file, mdl, "DIAG_MLD_DENSITY_DIFF", cs%MLDdensityDifference, &
                 "The density difference used to determine a diagnostic mixed "//&
                 "layer depth, MLD_user, following the definition of Levitus 1982. "//&
                 "The MLD is the depth at which the density is larger than the "//&
                 "surface density by the specified amount.", &
                 units='kg/m3', default=0.1, scale=us%kg_m3_to_R)
  call get_param(param_file, mdl, "DIAG_DEPTH_SUBML_N2", cs%dz_subML_N2, &
                 "The distance over which to calculate a diagnostic of the "//&
                 "stratification at the base of the mixed layer.", &
                 units='m', default=50.0, scale=us%m_to_Z)

  if (cs%id_dudt_dia > 0) call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
  if (cs%id_dvdt_dia > 0) call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)

  ! diagnostics for values prior to diabatic and prior to ALE
  cs%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, time, &
      'Zonal velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, time, &
      'Meridional velocity before diabatic forcing', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, time, &
      'Layer Thickness before diabatic forcing', &
      trim(thickness_units), conversion=gv%H_to_MKS, v_extensive=.true.)
  cs%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, time, &
      'Interface Heights before diabatic forcing', 'm')
  if (use_temperature) then
    cs%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, time, &
        'Potential Temperature', 'degC')
    cs%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, time, &
        'Salinity', 'PSU')
  endif


  !call set_diffusivity_init(Time, G, param_file, diag, CS%set_diff_CSp, CS%int_tide_CSp)
  cs%id_Kd_interface = register_diag_field('ocean_model', 'Kd_interface', diag%axesTi, time, &
      'Total diapycnal diffusivity at interfaces', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
  if (cs%use_energetic_PBL) then
      cs%id_Kd_ePBL = register_diag_field('ocean_model', 'Kd_ePBL', diag%axesTi, time, &
          'ePBL diapycnal diffusivity at interfaces', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
  endif

  cs%id_Kd_heat = register_diag_field('ocean_model', 'Kd_heat', diag%axesTi, time, &
      'Total diapycnal diffusivity for heat at interfaces', 'm2 s-1', conversion=us%Z2_T_to_m2_s, &
       cmor_field_name='difvho',                                                   &
       cmor_standard_name='ocean_vertical_heat_diffusivity',                       &
       cmor_long_name='Ocean vertical heat diffusivity')
  cs%id_Kd_salt = register_diag_field('ocean_model', 'Kd_salt', diag%axesTi, time, &
      'Total diapycnal diffusivity for salt at interfaces', 'm2 s-1',  conversion=us%Z2_T_to_m2_s, &
       cmor_field_name='difvso',                                                   &
       cmor_standard_name='ocean_vertical_salt_diffusivity',                       &
       cmor_long_name='Ocean vertical salt diffusivity')

  ! CS%useKPP is set to True if KPP-scheme is to be used, False otherwise.
  ! KPP_init() allocated CS%KPP_Csp and also sets CS%KPPisPassive
  cs%useKPP = kpp_init(param_file, g, gv, us, diag, time, cs%KPP_CSp, passive=cs%KPPisPassive)
  if (cs%useKPP) then
    allocate( cs%KPP_NLTheat(isd:ied,jsd:jed,nz+1) )   ; cs%KPP_NLTheat(:,:,:)   = 0.
    allocate( cs%KPP_NLTscalar(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_NLTscalar(:,:,:) = 0.
  endif
  if (cs%useKPP) then
    allocate( cs%KPP_buoy_flux(isd:ied,jsd:jed,nz+1) ) ; cs%KPP_buoy_flux(:,:,:) = 0.
    allocate( cs%KPP_temp_flux(isd:ied,jsd:jed) )      ; cs%KPP_temp_flux(:,:)   = 0.
    allocate( cs%KPP_salt_flux(isd:ied,jsd:jed) )      ; cs%KPP_salt_flux(:,:)   = 0.
  endif


  ! diagnostics for tendencies of temp and saln due to diabatic processes
  ! available only for ALE algorithm.
  ! diagnostics for tendencies of temp and heat due to frazil
  cs%id_diabatic_diff_h = register_diag_field('ocean_model', 'diabatic_diff_h', diag%axesTL, time, &
      long_name='Cell thickness used during diabatic diffusion', &
      units='m', conversion=gv%H_to_m, v_extensive=.true.)
  if (cs%useALEalgorithm) then
    cs%id_diabatic_diff_temp_tend = register_diag_field('ocean_model', &
        'diabatic_diff_temp_tendency', diag%axesTL, time,              &
        'Diabatic diffusion temperature tendency', 'degC s-1', conversion=us%s_to_T)
    if (cs%id_diabatic_diff_temp_tend > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    cs%id_diabatic_diff_saln_tend = register_diag_field('ocean_model',&
        'diabatic_diff_saln_tendency', diag%axesTL, time,             &
        'Diabatic diffusion salinity tendency', 'psu s-1', conversion=us%s_to_T)
    if (cs%id_diabatic_diff_saln_tend > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    cs%id_diabatic_diff_heat_tend = register_diag_field('ocean_model',                             &
        'diabatic_heat_tendency', diag%axesTL, time,                                               &
        'Diabatic diffusion heat tendency',                                                        &
        'W m-2', conversion=us%QRZ_T_to_W_m2, cmor_field_name='opottempdiff',            &
        cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'// &
                           'due_to_parameterized_dianeutral_mixing',                               &
        cmor_long_name='Tendency of sea water potential temperature expressed as heat content '//  &
                       'due to parameterized dianeutral mixing', &
        v_extensive=.true.)
    if (cs%id_diabatic_diff_heat_tend > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    cs%id_diabatic_diff_salt_tend = register_diag_field('ocean_model',                   &
        'diabatic_salt_tendency', diag%axesTL, time,                                     &
        'Diabatic diffusion of salt tendency',                                           &
        'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s, cmor_field_name='osaltdiff', &
        cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &
                           'due_to_parameterized_dianeutral_mixing',                     &
        cmor_long_name='Tendency of sea water salinity expressed as salt content '//     &
                       'due to parameterized dianeutral mixing', &
        v_extensive=.true.)
    if (cs%id_diabatic_diff_salt_tend > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    ! This diagnostic should equal to roundoff if all is working well.
    cs%id_diabatic_diff_heat_tend_2d = register_diag_field('ocean_model',                        &
        'diabatic_heat_tendency_2d', diag%axesT1, time,                                          &
        'Depth integrated diabatic diffusion heat tendency',                                     &
        'W m-2', conversion=us%QRZ_T_to_W_m2, cmor_field_name='opottempdiff_2d',      &
        cmor_standard_name='tendency_of_sea_water_potential_temperature_expressed_as_heat_content_'//&
                           'due_to_parameterized_dianeutral_mixing_depth_integrated',            &
        cmor_long_name='Tendency of sea water potential temperature expressed as heat content '//&
                       'due to parameterized dianeutral mixing depth integrated')
    if (cs%id_diabatic_diff_heat_tend_2d > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    ! This diagnostic should equal to roundoff if all is working well.
    cs%id_diabatic_diff_salt_tend_2d = register_diag_field('ocean_model',                &
        'diabatic_salt_tendency_2d', diag%axesT1, time,                                  &
        'Depth integrated diabatic diffusion salt tendency',                             &
        'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s, cmor_field_name='osaltdiff_2d',              &
        cmor_standard_name='tendency_of_sea_water_salinity_expressed_as_salt_content_'// &
                           'due_to_parameterized_dianeutral_mixing_depth_integrated',    &
        cmor_long_name='Tendency of sea water salinity expressed as salt content '//     &
                       'due to parameterized dianeutral mixing depth integrated')
    if (cs%id_diabatic_diff_salt_tend_2d > 0) then
      cs%diabatic_diff_tendency_diag = .true.
    endif

    ! diagnostics for tendencies of thickness temp and saln due to boundary forcing
    ! available only for ALE algorithm.
  ! diagnostics for tendencies of temp and heat due to frazil
    cs%id_boundary_forcing_h = register_diag_field('ocean_model', 'boundary_forcing_h', diag%axesTL, time, &
        long_name='Cell thickness after applying boundary forcing', &
        units='m', conversion=gv%H_to_m, v_extensive=.true.)
    cs%id_boundary_forcing_h_tendency = register_diag_field('ocean_model',   &
        'boundary_forcing_h_tendency', diag%axesTL, time,                &
        'Cell thickness tendency due to boundary forcing', 'm s-1', &
        conversion=gv%H_to_m*us%s_to_T, v_extensive=.true.)
    if (cs%id_boundary_forcing_h_tendency > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    cs%id_boundary_forcing_temp_tend = register_diag_field('ocean_model',&
        'boundary_forcing_temp_tendency', diag%axesTL, time,             &
        'Boundary forcing temperature tendency', 'degC s-1', conversion=us%s_to_T)
    if (cs%id_boundary_forcing_temp_tend > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    cs%id_boundary_forcing_saln_tend = register_diag_field('ocean_model',&
        'boundary_forcing_saln_tendency', diag%axesTL, time,             &
        'Boundary forcing saln tendency', 'psu s-1', conversion=us%s_to_T)
    if (cs%id_boundary_forcing_saln_tend > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    cs%id_boundary_forcing_heat_tend = register_diag_field('ocean_model',&
        'boundary_forcing_heat_tendency', diag%axesTL, time,             &
        'Boundary forcing heat tendency', &
        'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive = .true.)
    if (cs%id_boundary_forcing_heat_tend > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    cs%id_boundary_forcing_salt_tend = register_diag_field('ocean_model',&
        'boundary_forcing_salt_tendency', diag%axesTL, time,             &
        'Boundary forcing salt tendency', 'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s, &
        v_extensive = .true.)
    if (cs%id_boundary_forcing_salt_tend > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    ! This diagnostic should equal to surface heat flux if all is working well.
    cs%id_boundary_forcing_heat_tend_2d = register_diag_field('ocean_model',&
        'boundary_forcing_heat_tendency_2d', diag%axesT1, time,             &
        'Depth integrated boundary forcing of ocean heat', &
        'W m-2', conversion=us%QRZ_T_to_W_m2)
    if (cs%id_boundary_forcing_heat_tend_2d > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif

    ! This diagnostic should equal to surface salt flux if all is working well.
    cs%id_boundary_forcing_salt_tend_2d = register_diag_field('ocean_model',&
        'boundary_forcing_salt_tendency_2d', diag%axesT1, time,             &
        'Depth integrated boundary forcing of ocean salt', &
        'kg m-2 s-1', conversion=us%RZ_T_to_kg_m2s)
    if (cs%id_boundary_forcing_salt_tend_2d > 0) then
      cs%boundary_forcing_tendency_diag = .true.
    endif
  endif

  ! diagnostics for tendencies of temp and heat due to frazil
  cs%id_frazil_h = register_diag_field('ocean_model', 'frazil_h', diag%axesTL, time, &
      long_name='Cell Thickness', standard_name='cell_thickness', &
      units='m', conversion=gv%H_to_m, v_extensive=.true.)

  ! diagnostic for tendency of temp due to frazil
  cs%id_frazil_temp_tend = register_diag_field('ocean_model',&
      'frazil_temp_tendency', diag%axesTL, time,             &
      'Temperature tendency due to frazil formation', 'degC s-1', conversion=us%s_to_T)
  if (cs%id_frazil_temp_tend > 0) then
    cs%frazil_tendency_diag = .true.
  endif

  ! diagnostic for tendency of heat due to frazil
  cs%id_frazil_heat_tend = register_diag_field('ocean_model',&
      'frazil_heat_tendency', diag%axesTL, time,             &
      'Heat tendency due to frazil formation', &
      'W m-2', conversion=us%QRZ_T_to_W_m2, v_extensive=.true.)
  if (cs%id_frazil_heat_tend > 0) then
    cs%frazil_tendency_diag = .true.
  endif

  ! if all is working propertly, this diagnostic should equal to hfsifrazil
  cs%id_frazil_heat_tend_2d = register_diag_field('ocean_model',&
      'frazil_heat_tendency_2d', diag%axesT1, time,             &
      'Depth integrated heat tendency due to frazil formation', &
      'W m-2', conversion=us%QRZ_T_to_W_m2)
  if (cs%id_frazil_heat_tend_2d > 0) then
    cs%frazil_tendency_diag = .true.
  endif

  if (cs%frazil_tendency_diag) then
    allocate(cs%frazil_temp_diag(isd:ied,jsd:jed,nz) ) ; cs%frazil_temp_diag(:,:,:) = 0.
    allocate(cs%frazil_heat_diag(isd:ied,jsd:jed,nz) ) ; cs%frazil_heat_diag(:,:,:) = 0.
  endif

  ! CS%use_CVMix_conv is set to True if CVMix convection will be used, otherwise it is False.
  cs%use_CVMix_conv = cvmix_conv_init(time, g, gv, us, param_file, diag, cs%CVMix_conv_csp)

  call entrain_diffusive_init(time, g, gv, us, param_file, diag, cs%entrain_diffusive_CSp, &
                              just_read_params=cs%useALEalgorithm)

  ! initialize the geothermal heating module
  if (cs%use_geothermal) &
    call geothermal_init(time, g, gv, us, param_file, diag, cs%geothermal_CSp, usealealgorithm)

  ! initialize module for internal tide induced mixing
  if (cs%use_int_tides) then
    call int_tide_input_init(time, g, gv, us, param_file, diag, cs%int_tide_input_CSp, &
                             cs%int_tide_input)
    call internal_tides_init(time, g, gv, us, param_file, diag, cs%int_tide_CSp)
  endif

  ! initialize module for setting diffusivities
  call set_diffusivity_init(time, g, gv, us, param_file, diag, cs%set_diff_CSp, cs%int_tide_CSp, &
                            halo_ts=cs%halo_TS_diff, double_diffuse=cs%double_diffuse)

  if (cs%useKPP .and. (cs%double_diffuse .and. .not.cs%use_CVMix_ddiff)) &
    call mom_error(fatal, 'diabatic_driver_init: DOUBLE_DIFFUSION (old method) does not work '//&
                          'with KPP. Please set DOUBLE_DIFFUSION=False and USE_CVMIX_DDIFF=True.')

  ! set up the clocks for this module
  id_clock_entrain = cpu_clock_id('(Ocean diabatic entrain)', grain=clock_module)
  if (cs%bulkmixedlayer) &
    id_clock_mixedlayer = cpu_clock_id('(Ocean mixed layer)', grain=clock_module)
  id_clock_remap = cpu_clock_id('(Ocean vert remap)', grain=clock_module)
  if (cs%use_geothermal) &
    id_clock_geothermal = cpu_clock_id('(Ocean geothermal)', grain=clock_routine)
  id_clock_set_diffusivity = cpu_clock_id('(Ocean set_diffusivity)', grain=clock_module)
  id_clock_kpp = cpu_clock_id('(Ocean KPP)', grain=clock_module)
  id_clock_tracers = cpu_clock_id('(Ocean tracer_columns)', grain=clock_module_driver+5)
  if (cs%use_sponge) &
    id_clock_sponge = cpu_clock_id('(Ocean sponges)', grain=clock_module)
  id_clock_tridiag = cpu_clock_id('(Ocean diabatic tridiag)', grain=clock_routine)
  id_clock_pass = cpu_clock_id('(Ocean diabatic message passing)', grain=clock_routine)
  id_clock_differential_diff = -1 ; if (cs%double_diffuse .and. .not.cs%use_CVMix_ddiff) &
    id_clock_differential_diff = cpu_clock_id('(Ocean differential diffusion)', grain=clock_routine)

  ! initialize the auxiliary diabatic driver module
  call diabatic_aux_init(time, g, gv, us, param_file, diag, cs%diabatic_aux_CSp, &
                         cs%useALEalgorithm, cs%use_energetic_PBL)

  ! initialize the boundary layer modules
  if (cs%bulkmixedlayer) &
    call bulkmixedlayer_init(time, g, gv, us, param_file, diag, cs%bulkmixedlayer_CSp)
  if (cs%use_energetic_PBL) &
    call energetic_pbl_init(time, g, gv, us, param_file, diag, cs%energetic_PBL_CSp)

  call regularize_layers_init(time, g, gv, param_file, diag, cs%regularize_layers_CSp)

  if (cs%debug_energy_req) &
    call diapyc_energy_req_init(time, g, gv, us, param_file, diag, cs%diapyc_en_rec_CSp)

  ! obtain information about the number of bands for penetrative shortwave
  if (use_temperature) then
    call get_param(param_file, mdl, "PEN_SW_NBANDS", nbands, default=1)
    if (nbands > 0) then
      allocate(cs%optics)
      call opacity_init(time, g, gv, us, param_file, diag, cs%opacity_CSp, cs%optics)
    endif
  endif

  ! Initialize the diagnostic grid storage
  call diag_grid_storage_init(cs%diag_grids_prev, g, diag)

end subroutine diabatic_driver_init


!> Routine to close the diabatic driver module
subroutine diabatic_driver_end(CS)
  type(diabatic_cs), pointer :: CS    !< module control structure

  if (.not.associated(cs)) return

  call diabatic_aux_end(cs%diabatic_aux_CSp)

  call entrain_diffusive_end(cs%entrain_diffusive_CSp)
  call set_diffusivity_end(cs%set_diff_CSp)

  if (cs%useKPP) then
    deallocate( cs%KPP_buoy_flux )
    deallocate( cs%KPP_temp_flux )
    deallocate( cs%KPP_salt_flux )
    deallocate( cs%KPP_NLTheat )
    deallocate( cs%KPP_NLTscalar )
    call kpp_end(cs%KPP_CSp)
  endif

  if (cs%use_CVMix_conv) call cvmix_conv_end(cs%CVMix_conv_csp)

  if (cs%use_energetic_PBL) &
    call energetic_pbl_end(cs%energetic_PBL_CSp)
  if (cs%debug_energy_req) &
    call diapyc_energy_req_end(cs%diapyc_en_rec_CSp)

  if (associated(cs%optics)) then
    call opacity_end(cs%opacity_CSp, cs%optics)
    deallocate(cs%optics)
  endif

  ! GMM, the following is commented out because arrays in
  ! CS%diag_grids_prev are neither pointers or allocatables
  ! and, therefore, cannot be deallocated.

  !call diag_grid_storage_end(CS%diag_grids_prev)

  deallocate(cs)

end subroutine diabatic_driver_end


!> \namespace mom_diabatic_driver
!!
!!  By Robert Hallberg, 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.
!!
!!  \section section_diabatic Outline of MOM diabatic
!!
!!  * 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 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 treated
!!  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_driver