MOM_isopycnal_slopes.F90

!> Calculations of isoneutral slopes and stratification.
module mom_isopycnal_slopes

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

use mom_grid, only : ocean_grid_type
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_eos, only : calculate_density_derivs
use mom_open_boundary, only : ocean_obc_type, obc_none
use mom_open_boundary, only : obc_direction_e, obc_direction_w, obc_direction_n, obc_direction_s

implicit none ; private

#include <MOM_memory.h>

public calc_isoneutral_slopes, vert_fill_ts

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

contains

!> Calculate isopycnal slopes, and optionally return N2 used in calculation.
subroutine calc_isoneutral_slopes(G, GV, US, h, e, tv, dt_kappa_smooth, &
                                  slope_x, slope_y, N2_u, N2_v, halo, OBC) !, eta_to_m)
  type(ocean_grid_type),                       intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),                     intent(in)    :: gv   !< The ocean's vertical grid structure
  type(unit_scale_type),                       intent(in)    :: us   !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),    intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1),  intent(in)    :: e    !< Interface heights [Z ~> m] or units
                                                                     !! given by 1/eta_to_m)
  type(thermo_var_ptrs),                       intent(in)    :: tv   !< A structure pointing to various
                                                                     !! thermodynamic variables
  real,                                        intent(in)    :: dt_kappa_smooth !< A smoothing vertical diffusivity
                                                                     !! times a smoothing timescale [Z2 ~> m2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: slope_x !< Isopycnal slope in i-direction [nondim]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: slope_y !< Isopycnal slope in j-direction [nondim]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), &
                                     optional, intent(inout) :: n2_u !< Brunt-Vaisala frequency squared at
                                                                     !! interfaces between u-points [T-2 ~> s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), &
                                     optional, intent(inout) :: n2_v !< Brunt-Vaisala frequency squared at
                                                                     !! interfaces between u-points [T-2 ~> s-2]
  integer,                           optional, intent(in)    :: halo !< Halo width over which to compute
  type(ocean_obc_type),              optional, pointer       :: obc  !< Open boundaries control structure.

  ! real,                              optional, intent(in)    :: eta_to_m !< The conversion factor from the units
  !  (This argument has been tested but for now serves no purpose.)  !! of eta to m; US%Z_to_m by default.
  ! Local variables
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
    t, &          ! The temperature [degC], with the values in
                  ! in massless layers filled vertically by diffusion.
    s !, &          ! The filled salinity [ppt], with the values in
                  ! in massless layers filled vertically by diffusion.
!    Rho           ! Density itself, when a nonlinear equation of state is not in use [R ~> kg m-3].
  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
    pres          ! The pressure at an interface [R L2 T-2 ~> Pa].
  real, dimension(SZIB_(G)) :: &
    drho_dt_u, &  ! The derivative of density with temperature at u points [R degC-1 ~> kg m-3 degC-1].
    drho_ds_u     ! The derivative of density with salinity at u points [R ppt-1 ~> kg m-3 ppt-1].
  real, dimension(SZI_(G)) :: &
    drho_dt_v, &  ! The derivative of density with temperature at v points [R degC-1 ~> kg m-3 degC-1].
    drho_ds_v     ! The derivative of density with salinity at v points [R ppt-1 ~> kg m-3 ppt-1].
  real, dimension(SZIB_(G)) :: &
    t_u, &        ! Temperature on the interface at the u-point [degC].
    s_u, &        ! Salinity on the interface at the u-point [ppt].
    pres_u        ! Pressure on the interface at the u-point [R L2 T-2 ~> Pa].
  real, dimension(SZI_(G)) :: &
    t_v, &        ! Temperature on the interface at the v-point [degC].
    s_v, &        ! Salinity on the interface at the v-point [ppt].
    pres_v        ! Pressure on the interface at the v-point [R L2 T-2 ~> Pa].
  real :: drdia, drdib  ! Along layer zonal- and meridional- potential density
  real :: drdja, drdjb  ! gradients in the layers above (A) and below (B) the
                        ! interface times the grid spacing [R ~> kg m-3].
  real :: drdkl, drdkr  ! Vertical density differences across an interface [R ~> kg m-3].
  real :: hg2a, hg2b    ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
  real :: hg2l, hg2r    ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
  real :: haa, hab, hal, har  ! Arithmetic mean thicknesses [H ~> m or kg m-2].
  real :: dzal, dzar    ! Temporary thicknesses in eta units [Z ~> m].
  real :: wta, wtb, wtl, wtr  ! Unscaled weights, with various units.
  real :: drdx, drdy    ! Zonal and meridional density gradients [R L-1 ~> kg m-4].
  real :: drdz          ! Vertical density gradient [R Z-1 ~> kg m-4].
  real :: slope         ! The slope of density surfaces, calculated in a way
                        ! that is always between -1 and 1.
  real :: mag_grad2     ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8].
  real :: slope2_ratio  ! The ratio of the slope squared to slope_max squared.
  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 :: dz_neglect    ! A change in interface heighs that is so small it is usually lost
                        ! in roundoff and can be neglected [Z ~> m].
  logical :: use_eos    ! If true, density is calculated from T & S using an equation of state.
  real :: g_rho0        ! The gravitational acceleration divided by density [Z2 T-2 R-1 ~> m5 kg-2 s-2]
  real :: z_to_l        ! A conversion factor between from units for e to the
                        ! units for lateral distances.
  real :: l_to_z        ! A conversion factor between from units for lateral distances
                        ! to the units for e.
  real :: h_to_z        ! A conversion factor from thickness units to the units of e.

  logical :: present_n2_u, present_n2_v
  integer, dimension(2) :: eosdom_u, eosdom_v ! Domains for the equation of state calculations at u and v points
  integer :: is, ie, js, je, nz, isdb
  integer :: i, j, k
  integer :: l_seg
  logical :: local_open_u_bc, local_open_v_bc

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

  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect**2
  z_to_l = us%Z_to_L ; h_to_z = gv%H_to_Z
  ! if (present(eta_to_m)) then
  !   Z_to_L = eta_to_m*US%m_to_L ; H_to_Z = GV%H_to_m / eta_to_m
  ! endif
  l_to_z = 1.0 / z_to_l
  dz_neglect = gv%H_subroundoff * h_to_z

  local_open_u_bc = .false.
  local_open_v_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_u_bc = obc%open_u_BCs_exist_globally
    local_open_v_bc = obc%open_v_BCs_exist_globally
  endif ; endif

  use_eos = associated(tv%eqn_of_state)

  present_n2_u = PRESENT(n2_u)
  present_n2_v = PRESENT(n2_v)
  g_rho0 = (us%L_to_Z*l_to_z*gv%g_Earth) / gv%Rho0
  if (present_n2_u) then
    do j=js,je ; do i=is-1,ie
      n2_u(i,j,1) = 0.
      n2_u(i,j,nz+1) = 0.
    enddo ; enddo
  endif
  if (present_n2_v) then
    do j=js-1,je ; do i=is,ie
      n2_v(i,j,1) = 0.
      n2_v(i,j,nz+1) = 0.
    enddo ; enddo
  endif

  if (use_eos) then
    if (present(halo)) then
      call vert_fill_ts(h, tv%T, tv%S, dt_kappa_smooth, t, s, g, gv, halo+1)
    else
      call vert_fill_ts(h, tv%T, tv%S, dt_kappa_smooth, t, s, g, gv, 1)
    endif
  endif

  ! Find the maximum and minimum permitted streamfunction.
  if (associated(tv%p_surf)) then
    !$OMP parallel do default(shared)
    do j=js-1,je+1 ; do i=is-1,ie+1
      pres(i,j,1) = tv%p_surf(i,j)
    enddo ; enddo
  else
    !$OMP parallel do default(shared)
    do j=js-1,je+1 ; do i=is-1,ie+1
      pres(i,j,1) = 0.0
    enddo ; enddo
  endif
  !$OMP parallel do default(shared)
  do j=js-1,je+1
    do k=1,nz ; do i=is-1,ie+1
      pres(i,j,k+1) = pres(i,j,k) + gv%g_Earth * gv%H_to_RZ * h(i,j,k)
    enddo ; enddo
  enddo

  eosdom_u(1) = is-1 - (g%IsdB-1) ; eosdom_u(2) = ie - (g%IsdB-1)

  !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv,h,e, &
  !$OMP                                  h_neglect,dz_neglect,Z_to_L,L_to_Z,H_to_Z,h_neglect2, &
  !$OMP                                  present_N2_u,G_Rho0,N2_u,slope_x,EOSdom_u,local_open_u_BC, &
  !$OMP                                  OBC) &
  !$OMP                          private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u,      &
  !$OMP                                  drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
  !$OMP                                  haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz,  &
  !$OMP                                  drdx,mag_grad2,Slope,slope2_Ratio,l_seg)
  do j=js,je ; do k=nz,2,-1
    if (.not.(use_eos)) then
      drdia = 0.0 ; drdib = 0.0
      drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = gv%Rlay(k)-gv%Rlay(k-1)
    endif

    ! Calculate the zonal isopycnal slope.
    if (use_eos) then
      do i=is-1,ie
        pres_u(i) = 0.5*(pres(i,j,k) + pres(i+1,j,k))
        t_u(i) = 0.25*((t(i,j,k) + t(i+1,j,k)) + (t(i,j,k-1) + t(i+1,j,k-1)))
        s_u(i) = 0.25*((s(i,j,k) + s(i+1,j,k)) + (s(i,j,k-1) + s(i+1,j,k-1)))
      enddo
      call calculate_density_derivs(t_u, s_u, pres_u, drho_dt_u, drho_ds_u, &
                                    tv%eqn_of_state, eosdom_u)
    endif

    do i=is-1,ie
      if (use_eos) then
        ! Estimate the horizontal density gradients along layers.
        drdia = drho_dt_u(i) * (t(i+1,j,k-1)-t(i,j,k-1)) + &
                drho_ds_u(i) * (s(i+1,j,k-1)-s(i,j,k-1))
        drdib = drho_dt_u(i) * (t(i+1,j,k)-t(i,j,k)) + &
                drho_ds_u(i) * (s(i+1,j,k)-s(i,j,k))

        ! Estimate the vertical density gradients times the grid spacing.
        drdkl = (drho_dt_u(i) * (t(i,j,k)-t(i,j,k-1)) + &
                 drho_ds_u(i) * (s(i,j,k)-s(i,j,k-1)))
        drdkr = (drho_dt_u(i) * (t(i+1,j,k)-t(i+1,j,k-1)) + &
                 drho_ds_u(i) * (s(i+1,j,k)-s(i+1,j,k-1)))
      endif


      if (use_eos) then
        hg2a = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
        hg2b = h(i,j,k)*h(i+1,j,k) + h_neglect2
        hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
        hg2r = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
        haa = 0.5*(h(i,j,k-1) + h(i+1,j,k-1))
        hab = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect
        hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
        har = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
        if (gv%Boussinesq) then
          dzal = hal * h_to_z ; dzar = har * h_to_z
        else
          dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
          dzar = 0.5*(e(i+1,j,k-1) - e(i+1,j,k+1)) + dz_neglect
        endif
        ! Use the harmonic mean thicknesses to weight the horizontal gradients.
        ! These unnormalized weights have been rearranged to minimize divisions.
        wta = hg2a*hab ; wtb = hg2b*haa
        wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)

        drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
        ! The expression for drdz above is mathematically equivalent to:
        !   drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
        !          ((hg2L/haL) + (hg2R/haR))
        ! This is the gradient of density along geopotentials.
        drdx = ((wta * drdia + wtb * drdib) / (wta + wtb) - &
                drdz * (e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)

        ! This estimate of slope is accurate for small slopes, but bounded
        ! to be between -1 and 1.
        mag_grad2 = drdx**2 + (l_to_z*drdz)**2
        if (mag_grad2 > 0.0) then
          slope_x(i,j,k) = drdx / sqrt(mag_grad2)
        else ! Just in case mag_grad2 = 0 ever.
          slope_x(i,j,k) = 0.0
        endif

        if (present_n2_u) n2_u(i,j,k) = g_rho0 * drdz * g%mask2dCu(i,j) ! Square of buoyancy frequency [T-2 ~> s-2]

      else ! With .not.use_EOS, the layers are constant density.
        slope_x(i,j,k) = (z_to_l*(e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
      endif
      if (local_open_u_bc) then
        l_seg = obc%segnum_u(i,j)
        if (l_seg /= obc_none) then
          if (obc%segment(l_seg)%open) then
            slope_x(i,j,k) = 0.
            ! This and/or the masking code below is to make slopes match inside
            ! land mask. Might not be necessary except for DEBUG output.
!           if (OBC%segment(OBC%segnum_u(I,j))%direction == OBC_DIRECTION_E) then
!             slope_x(I+1,j,K) = 0.
!           else
!             slope_x(I-1,j,K) = 0.
!           endif
          endif
        endif
        slope_x(i,j,k) = slope_x(i,j,k) * max(g%mask2dT(i,j),g%mask2dT(i+1,j))
      endif

    enddo ! I
  enddo ; enddo ! end of j-loop

  eosdom_v(1) = is - (g%isd-1) ; eosdom_v(2) = ie - (g%isd-1)

  ! Calculate the meridional isopycnal slope.
  !$OMP parallel do default(none) shared(nz,is,ie,js,je,IsdB,use_EOS,G,GV,US,pres,T,S,tv, &
  !$OMP                                  h,h_neglect,e,dz_neglect,Z_to_L,L_to_Z,H_to_Z, &
  !$OMP                                  h_neglect2,present_N2_v,G_Rho0,N2_v,slope_y,EOSdom_v, &
  !$OMP                                  local_open_v_BC,OBC) &
  !$OMP                          private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v,      &
  !$OMP                                  drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
  !$OMP                                  haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz,  &
  !$OMP                                  drdy,mag_grad2,Slope,slope2_Ratio,l_seg)
  do j=js-1,je ; do k=nz,2,-1
    if (.not.(use_eos)) then
      drdja = 0.0 ; drdjb = 0.0
      drdkl = gv%Rlay(k)-gv%Rlay(k-1) ; drdkr = gv%Rlay(k)-gv%Rlay(k-1)
    endif

    if (use_eos) then
      do i=is,ie
        pres_v(i) = 0.5*(pres(i,j,k) + pres(i,j+1,k))
        t_v(i) = 0.25*((t(i,j,k) + t(i,j+1,k)) + (t(i,j,k-1) + t(i,j+1,k-1)))
        s_v(i) = 0.25*((s(i,j,k) + s(i,j+1,k)) + (s(i,j,k-1) + s(i,j+1,k-1)))
      enddo
      call calculate_density_derivs(t_v, s_v, pres_v, drho_dt_v, drho_ds_v, tv%eqn_of_state, &
                                    eosdom_v)
    endif
    do i=is,ie
      if (use_eos) then
        ! Estimate the horizontal density gradients along layers.
        drdja = drho_dt_v(i) * (t(i,j+1,k-1)-t(i,j,k-1)) + &
                drho_ds_v(i) * (s(i,j+1,k-1)-s(i,j,k-1))
        drdjb = drho_dt_v(i) * (t(i,j+1,k)-t(i,j,k)) + &
                drho_ds_v(i) * (s(i,j+1,k)-s(i,j,k))

        ! Estimate the vertical density gradients times the grid spacing.
        drdkl = (drho_dt_v(i) * (t(i,j,k)-t(i,j,k-1)) + &
                 drho_ds_v(i) * (s(i,j,k)-s(i,j,k-1)))
        drdkr = (drho_dt_v(i) * (t(i,j+1,k)-t(i,j+1,k-1)) + &
                 drho_ds_v(i) * (s(i,j+1,k)-s(i,j+1,k-1)))
      endif

      if (use_eos) then
        hg2a = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
        hg2b = h(i,j,k)*h(i,j+1,k) + h_neglect2
        hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
        hg2r = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
        haa = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
        hab = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect
        hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
        har = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
        if (gv%Boussinesq) then
          dzal = hal * h_to_z ; dzar = har * h_to_z
        else
          dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
          dzar = 0.5*(e(i,j+1,k-1) - e(i,j+1,k+1)) + dz_neglect
        endif
        ! Use the harmonic mean thicknesses to weight the horizontal gradients.
        ! These unnormalized weights have been rearranged to minimize divisions.
        wta = hg2a*hab ; wtb = hg2b*haa
        wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)

        drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
        ! The expression for drdz above is mathematically equivalent to:
        !   drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
        !          ((hg2L/haL) + (hg2R/haR))
        ! This is the gradient of density along geopotentials.
        drdy = ((wta * drdja + wtb * drdjb) / (wta + wtb) - &
                drdz * (e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)

        ! This estimate of slope is accurate for small slopes, but bounded
        ! to be between -1 and 1.
        mag_grad2 = drdy**2 + (l_to_z*drdz)**2
        if (mag_grad2 > 0.0) then
          slope_y(i,j,k) = drdy / sqrt(mag_grad2)
        else ! Just in case mag_grad2 = 0 ever.
          slope_y(i,j,k) = 0.0
        endif

        if (present_n2_v) n2_v(i,j,k) = g_rho0 * drdz * g%mask2dCv(i,j) ! Square of buoyancy frequency [T-2 ~> s-2]

      else ! With .not.use_EOS, the layers are constant density.
        slope_y(i,j,k) = (z_to_l*(e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
      endif
      if (local_open_v_bc) then
        l_seg = obc%segnum_v(i,j)
        if (l_seg /= obc_none) then
          if (obc%segment(l_seg)%open) then
            slope_y(i,j,k) = 0.
            ! This and/or the masking code below is to make slopes match inside
            ! land mask. Might not be necessary except for DEBUG output.
!           if (OBC%segment(OBC%segnum_v(i,J))%direction == OBC_DIRECTION_N) then
!             slope_y(i,J+1,K) = 0.
!           else
!             slope_y(i,J-1,K) = 0.
!           endif
          endif
        endif
        slope_y(i,j,k) = slope_y(i,j,k) * max(g%mask2dT(i,j),g%mask2dT(i,j+1))
      endif

    enddo ! i
  enddo ; enddo ! end of j-loop

end subroutine calc_isoneutral_slopes

!> Returns tracer arrays (nominally T and S) with massless layers filled with
!! sensible values, by diffusing vertically with a small but constant diffusivity.
subroutine vert_fill_ts(h, T_in, S_in, kappa_dt, T_f, S_f, G, GV, halo_here, larger_h_denom)
  type(ocean_grid_type),                    intent(in)  :: g    !< The ocean's grid structure
  type(verticalgrid_type),                  intent(in)  :: gv   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)  :: h    !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)  :: t_in !< Input temperature [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in)  :: s_in !< Input salinity [ppt]
  real,                                     intent(in)  :: kappa_dt !< A vertical diffusivity to use for smoothing
                                                                !! times a smoothing timescale [Z2 ~> m2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: t_f  !< Filled temperature [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: s_f  !< Filled salinity [ppt]
  integer,                        optional, intent(in)  :: halo_here !< Number of halo points to work on,
                                                                !! 0 by default
  logical,                        optional, intent(in)  :: larger_h_denom !< Present and true, add a large
                                                                !! enough minimal thickness in the denominator of
                                                                !! the flux calculations so that the fluxes are
                                                                !! never so large as eliminate the transmission
                                                                !! of information across groups of massless layers.
  ! Local variables
  real :: ent(szi_(g),szk_(g)+1)   ! The diffusive entrainment (kappa*dt)/dz
                                   ! between layers in a timestep [H ~> m or kg m-2].
  real :: b1(szi_(g)), d1(szi_(g)) ! b1, c1, and d1 are variables used by the
  real :: c1(szi_(g),szk_(g))      ! tridiagonal solver.
  real :: kap_dt_x2                ! The 2*kappa_dt converted to H units [H2 ~> m2 or kg2 m-4].
  real :: h_neglect                ! A negligible thickness [H ~> m or kg m-2], to allow for zero thicknesses.
  real :: h0                       ! A negligible thickness to allow for zero thickness layers without
                                   ! completely decouping groups of layers [H ~> m or kg m-2].
                                   ! Often 0 < h_neglect << h0.
  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].
  integer :: i, j, k, is, ie, js, je, nz, halo

  halo=0 ; if (present(halo_here)) halo = max(halo_here,0)

  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo ; nz = gv%ke

  h_neglect = gv%H_subroundoff
  kap_dt_x2 = (2.0*kappa_dt)*gv%Z_to_H**2
  h0 = h_neglect
  if (present(larger_h_denom)) then
    if (larger_h_denom) h0 = 1.0e-16*sqrt(kappa_dt)*gv%Z_to_H
  endif

  if (kap_dt_x2 <= 0.0) then
    !$OMP parallel do default(shared)
    do k=1,nz ; do j=js,je ; do i=is,ie
      t_f(i,j,k) = t_in(i,j,k) ; s_f(i,j,k) = s_in(i,j,k)
    enddo ; enddo ; enddo
  else
   !$OMP parallel do default(shared) private(ent,b1,d1,c1,h_tr)
    do j=js,je
      do i=is,ie
        ent(i,2) = kap_dt_x2 / ((h(i,j,1)+h(i,j,2)) + h0)
        h_tr = h(i,j,1) + h_neglect
        b1(i) = 1.0 / (h_tr + ent(i,2))
        d1(i) = b1(i) * h_tr
        t_f(i,j,1) = (b1(i)*h_tr)*t_in(i,j,1)
        s_f(i,j,1) = (b1(i)*h_tr)*s_in(i,j,1)
      enddo
      do k=2,nz-1 ; do i=is,ie
        ent(i,k+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0)
        h_tr = h(i,j,k) + h_neglect
        c1(i,k) = ent(i,k) * b1(i)
        b1(i) = 1.0 / ((h_tr + d1(i)*ent(i,k)) + ent(i,k+1))
        d1(i) = b1(i) * (h_tr + d1(i)*ent(i,k))
        t_f(i,j,k) = b1(i) * (h_tr*t_in(i,j,k) + ent(i,k)*t_f(i,j,k-1))
        s_f(i,j,k) = b1(i) * (h_tr*s_in(i,j,k) + ent(i,k)*s_f(i,j,k-1))
      enddo ; enddo
      do i=is,ie
        c1(i,nz) = ent(i,nz) * b1(i)
        h_tr = h(i,j,nz) + h_neglect
        b1(i) = 1.0 / (h_tr + d1(i)*ent(i,nz))
        t_f(i,j,nz) = b1(i) * (h_tr*t_in(i,j,nz) + ent(i,nz)*t_f(i,j,nz-1))
        s_f(i,j,nz) = b1(i) * (h_tr*s_in(i,j,nz) + ent(i,nz)*s_f(i,j,nz-1))
      enddo
      do k=nz-1,1,-1 ; do i=is,ie
        t_f(i,j,k) = t_f(i,j,k) + c1(i,k+1)*t_f(i,j,k+1)
        s_f(i,j,k) = s_f(i,j,k) + c1(i,k+1)*s_f(i,j,k+1)
      enddo ; enddo
    enddo
  endif

end subroutine vert_fill_ts

end module mom_isopycnal_slopes