MOM_full_convection.F90

!> Does full convective adjustment of unstable regions via a strong diffusivity.
module mom_full_convection

! 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, eos_domain

implicit none ; private

#include <MOM_memory.h>

public full_convection

contains

!> Calculate new temperatures and salinities that have been subject to full convective mixing.
subroutine full_convection(G, GV, US, h, tv, T_adj, S_adj, p_surf, Kddt_smooth, &
                           Kddt_convect, halo)
  type(ocean_grid_type),   intent(in)    :: g     !< The ocean's grid structure
  type(verticalgrid_type), intent(in)    :: gv    !< The ocean's vertical grid structure
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: h     !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),   intent(in)    :: tv    !< A structure pointing to various
                                                  !! thermodynamic variables
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(out)   :: t_adj !< Adjusted potential temperature [degC].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(out)   :: s_adj !< Adjusted salinity [ppt].
  real, dimension(:,:),    pointer       :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa] (or NULL).
  real,                    intent(in)    :: kddt_smooth  !< A smoothing vertical
                                                  !! diffusivity times a timestep [H2 ~> m2 or kg2 m-4].
  real,          optional, intent(in)    :: kddt_convect !< A large convecting vertical
                                                  !! diffusivity times a timestep [H2 ~> m2 or kg2 m-4].
  integer,       optional, intent(in)    :: halo  !< Halo width over which to compute

  ! Local variables
  real, dimension(SZI_(G),SZK_(G)+1) :: &
    drho_dt, &  ! The derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
    drho_ds     ! The derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
  real :: h_neglect, h0 ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected [H ~> m or kg m-2].
! logical :: use_EOS    ! If true, density is calculated from T & S using an equation of state.
  real, dimension(SZI_(G),SZK0_(G)) :: &
    te_a, & ! A partially updated temperature estimate including the influnce from
            ! mixing with layers above rescaled by a factor of d_a [degC].
            ! This array is discreted on tracer cells, but contains an extra
            ! layer at the top for algorithmic convenience.
    se_a    ! A partially updated salinity estimate including the influnce from
            ! mixing with layers above rescaled by a factor of d_a [ppt].
            ! This array is discreted on tracer cells, but contains an extra
            ! layer at the top for algorithmic convenience.
  real, dimension(SZI_(G),SZK_(G)+1) :: &
    te_b, & ! A partially updated temperature estimate including the influnce from
            ! mixing with layers below rescaled by a factor of d_b [degC].
            ! This array is discreted on tracer cells, but contains an extra
            ! layer at the bottom for algorithmic convenience.
    se_b    ! A partially updated salinity estimate including the influnce from
            ! mixing with layers below rescaled by a factor of d_b [ppt].
            ! This array is discreted on tracer cells, but contains an extra
            ! layer at the bottom for algorithmic convenience.
  real, dimension(SZI_(G),SZK_(G)+1) :: &
    c_a, &  ! The fractional influence of the properties of the layer below
            ! in the final properies with a downward-first solver, nondim.
    d_a, &  ! The fractional influence of the properties of the layer in question
            ! and layers above in the final properies with a downward-first solver, nondim.
            ! d_a = 1.0 - c_a
    c_b, &  ! The fractional influence of the properties of the layer above
            ! in the final properies with a upward-first solver, nondim.
    d_b     ! The fractional influence of the properties of the layer in question
            ! and layers below in the final properies with a upward-first solver, nondim.
            ! d_b = 1.0 - c_b
  real, dimension(SZI_(G),SZK_(G)+1) :: &
    mix     !< The amount of mixing across the interface between layers [H ~> m or kg m-2].
  real :: mix_len  ! The length-scale of mixing, when it is active [H ~> m or kg m-2]
  real :: h_b, h_a ! The thicknessses of the layers above and below an interface [H ~> m or kg m-2]
  real :: b_b, b_a ! Inverse pivots used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].

  real :: kap_dt_x2 ! The product of 2*kappa*dt [H2 ~> m2 or kg2 m-4].

  logical, dimension(SZI_(G)) :: do_i ! Do more work on this column.
  logical, dimension(SZI_(G)) :: last_down ! The last setup pass was downward.
  integer, dimension(SZI_(G)) :: change_ct ! The number of interfaces where the
                         ! mixing has changed this iteration.
  integer :: changed_col ! The number of colums whose mixing changed.
  integer :: i, j, k, is, ie, js, je, nz, itt

  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

  if (.not.associated(tv%eqn_of_state)) return

  h_neglect = gv%H_subroundoff
  kap_dt_x2 = 0.0
  if (present(kddt_convect)) kap_dt_x2 = 2.0*kddt_convect
  mix_len = (1.0e20 * nz) * (g%max_depth * gv%Z_to_H)
  h0 = 1.0e-16*sqrt(kddt_smooth) + h_neglect

  do j=js,je
    mix(:,:) = 0.0 ; d_b(:,:) = 1.0
    ! These would be Te_b(:,:) = tv%T(:,j,:), etc., but the values are not used
    te_b(:,:) = 0.0 ; se_b(:,:) = 0.0

    call smoothed_drdt_drds(h, tv, kddt_smooth, drho_dt, drho_ds, g, gv, us, j, p_surf, halo)

    do i=is,ie
      do_i(i) = (g%mask2dT(i,j) > 0.0)

      d_a(i,1) = 1.0
      last_down(i) = .true. ! This is set for debuggers.
      ! These are extra values are used for convenience in the stability test
      te_a(i,0) = 0.0 ; se_a(i,0) = 0.0
    enddo

    do itt=1,nz ! At least 2 interfaces will change with each full pass, or the
                ! iterations stop, so the maximum count of nz is very conservative.

      do i=is,ie ; change_ct(i) = 0 ; enddo
      ! Move down the water column, finding unstable interfaces, and building up the
      ! temporary arrays for the tridiagonal solver.
      do k=2,nz ; do i=is,ie ; if (do_i(i)) then

        h_a = h(i,j,k-1) + h_neglect ; h_b = h(i,j,k) + h_neglect
        if (mix(i,k) <= 0.0) then
          if (is_unstable(drho_dt(i,k), drho_ds(i,k), h_a, h_b, mix(i,k-1), mix(i,k+1), &
                          tv%T(i,j,k-1), tv%T(i,j,k), tv%S(i,j,k-1), tv%S(i,j,k), &
                          te_a(i,k-2), te_b(i,k+1), se_a(i,k-2), se_b(i,k+1), &
                          d_a(i,k-1), d_b(i,k+1))) then
            mix(i,k) = mix_len
            if (kap_dt_x2 > 0.0) mix(i,k) = kap_dt_x2 / ((h(i,j,k-1)+h(i,j,k)) + h0)
            change_ct(i) = change_ct(i) + 1
          endif
        endif

        b_a = 1.0 / ((h_a + d_a(i,k-1)*mix(i,k-1)) + mix(i,k))
        if (mix(i,k) <= 0.0) then
          c_a(i,k) = 0.0 ; d_a(i,k) = 1.0
        else
          d_a(i,k) = b_a * (h_a + d_a(i,k-1)*mix(i,k-1)) ! = 1.0-c_a(i,K)
          c_a(i,k) = 1.0 ; if (d_a(i,k) > epsilon(b_a)) c_a(i,k) = b_a * mix(i,k)
        endif

        if (k>2) then
          te_a(i,k-1) = b_a * (h_a*tv%T(i,j,k-1) + mix(i,k-1)*te_a(i,k-2))
          se_a(i,k-1) = b_a * (h_a*tv%S(i,j,k-1) + mix(i,k-1)*se_a(i,k-2))
        else
          te_a(i,k-1) = b_a * (h_a*tv%T(i,j,k-1))
          se_a(i,k-1) = b_a * (h_a*tv%S(i,j,k-1))
        endif
      endif ; enddo ; enddo

      ! Determine which columns might have further instabilities.
      changed_col = 0
      do i=is,ie ; if (do_i(i)) then
        if (change_ct(i) == 0) then
          last_down(i) = .true. ; do_i(i) = .false.
        else
          changed_col = changed_col + 1 ; change_ct(i) = 0
        endif
      endif ; enddo
      if (changed_col == 0) exit ! No more columns are unstable.

      ! This is the same as above, but with the direction reversed (bottom to top)
      do k=nz,2,-1 ; do i=is,ie ; if (do_i(i)) then

        h_a = h(i,j,k-1) + h_neglect ; h_b = h(i,j,k) + h_neglect
        if (mix(i,k) <= 0.0) then
          if (is_unstable(drho_dt(i,k), drho_ds(i,k), h_a, h_b, mix(i,k-1), mix(i,k+1), &
                          tv%T(i,j,k-1), tv%T(i,j,k), tv%S(i,j,k-1), tv%S(i,j,k), &
                          te_a(i,k-2), te_b(i,k+1), se_a(i,k-2), se_b(i,k+1), &
                          d_a(i,k-1), d_b(i,k+1))) then
            mix(i,k) = mix_len
            if (kap_dt_x2 > 0.0) mix(i,k) = kap_dt_x2 / ((h(i,j,k-1)+h(i,j,k)) + h0)
            change_ct(i) = change_ct(i) + 1
          endif
        endif

        b_b = 1.0 / ((h_b + d_b(i,k+1)*mix(i,k+1)) + mix(i,k))
        if (mix(i,k) <= 0.0) then
          c_b(i,k) = 0.0 ; d_b(i,k) = 1.0
        else
          d_b(i,k) = b_b * (h_b + d_b(i,k+1)*mix(i,k+1)) ! = 1.0-c_b(i,K)
          c_b(i,k) = 1.0 ; if (d_b(i,k) > epsilon(b_b)) c_b(i,k) = b_b * mix(i,k)
        endif

        if (k<nz) then
          te_b(i,k) = b_b * (h_b*tv%T(i,j,k) + mix(i,k+1)*te_b(i,k+1))
          se_b(i,k) = b_b * (h_b*tv%S(i,j,k) + mix(i,k+1)*se_b(i,k+1))
        else
          te_b(i,k) = b_b * (h_b*tv%T(i,j,k))
          se_b(i,k) = b_b * (h_b*tv%S(i,j,k))
        endif
      endif ; enddo ; enddo

      ! Determine which columns might have further instabilities.
      changed_col = 0
      do i=is,ie ; if (do_i(i)) then
        if (change_ct(i) == 0) then
          last_down(i) = .false. ; do_i(i) = .false.
        else
          changed_col = changed_col + 1 ; change_ct(i) = 0
        endif
      endif ; enddo
      if (changed_col == 0) exit ! No more columns are unstable.

    enddo  ! End of iterations, all columns are now stable.

    ! Do the final return pass on the columns where the penultimate pass was downward.
    do i=is,ie ; do_i(i) = ((g%mask2dT(i,j) > 0.0) .and. last_down(i)) ; enddo
    do i=is,ie ; if (do_i(i)) then
      h_a = h(i,j,nz) + h_neglect
      b_a = 1.0 / (h_a + d_a(i,nz)*mix(i,nz))
      t_adj(i,j,nz) = b_a * (h_a*tv%T(i,j,nz) + mix(i,nz)*te_a(i,nz-1))
      s_adj(i,j,nz) = b_a * (h_a*tv%S(i,j,nz) + mix(i,nz)*se_a(i,nz-1))
    endif ; enddo
    do k=nz-1,1,-1 ; do i=is,ie ; if (do_i(i)) then
      t_adj(i,j,k) = te_a(i,k) + c_a(i,k+1)*t_adj(i,j,k+1)
      s_adj(i,j,k) = se_a(i,k) + c_a(i,k+1)*s_adj(i,j,k+1)
    endif ; enddo ; enddo

    do i=is,ie ; if (do_i(i)) then
      k = 1 ! A hook for debugging.
    endif ; enddo

    ! Do the final return pass on the columns where the penultimate pass was upward.
    ! Also do a simple copy of T & S values on land points.
    do i=is,ie
      do_i(i) = ((g%mask2dT(i,j) > 0.0) .and. .not.last_down(i))
      if (do_i(i)) then
        h_b = h(i,j,1) + h_neglect
        b_b = 1.0 / (h_b + d_b(i,2)*mix(i,2))
        t_adj(i,j,1) = b_b * (h_b*tv%T(i,j,1) + mix(i,2)*te_b(i,2))
        s_adj(i,j,1) = b_b * (h_b*tv%S(i,j,1) + mix(i,2)*se_b(i,2))
      elseif (g%mask2dT(i,j) <= 0.0) then
        t_adj(i,j,1) = tv%T(i,j,1) ; s_adj(i,j,1) = tv%S(i,j,1)
      endif
    enddo
    do k=2,nz ; do i=is,ie
      if (do_i(i)) then
        t_adj(i,j,k) = te_b(i,k) + c_b(i,k)*t_adj(i,j,k-1)
        s_adj(i,j,k) = se_b(i,k) + c_b(i,k)*s_adj(i,j,k-1)
      elseif (g%mask2dT(i,j) <= 0.0) then
        t_adj(i,j,k) = tv%T(i,j,k) ; s_adj(i,j,k) = tv%S(i,j,k)
      endif
    enddo ; enddo

    do i=is,ie ; if (do_i(i)) then
      k = 1 ! A hook for debugging.
    endif ; enddo

  enddo ! j-loop

  k = 1 ! A hook for debugging.

  ! The following set of expressions for the final values are derived from the the partial
  ! updates for the estimated temperatures and salinities around an interface, then directly
  ! solving for the final temperatures and salinities.  They are here for later reference
  ! and to document an intermediate step in the stability calculation.
    ! hp_a = (h_a + d_a(i,K-1)*mix(i,K-1))
    ! hp_b = (h_b + d_b(i,K+1)*mix(i,K+1))
    ! b2_c = 1.0 / (hp_a*hp_b + (hp_a + hp_b) * mix(i,K))
    ! Th_a = h_a*tv%T(i,j,k-1) + mix(i,K-1)*Te_a(i,k-2)
    ! Th_b = h_b*tv%T(i,j,k)   + mix(i,K+1)*Te_b(i,k+1)
    ! T_fin(i,k)   = ( (hp_a + mix(i,K)) * Th_b  + Th_a * mix(i,K) ) * b2_c
    ! T_fin(i,k-1) = ( (hp_b + mix(i,K)) * Th_a  + Th_b * mix(i,K) ) * b2_c
    ! Sh_a = h_a*tv%S(i,j,k-1) + mix(i,K-1)*Se_a(i,k-2)
    ! Sh_b = h_b*tv%S(i,j,k)   + mix(i,K+1)*Se_b(i,k+1)
    ! S_fin(i,k)   = ( (hp_a + mix(i,K)) * Sh_b  + Sh_a * mix(i,K) ) * b2_c
    ! S_fin(i,k-1) = ( (hp_b + mix(i,K)) * Sh_a  + Sh_b * mix(i,K) ) * b2_c

end subroutine full_convection

!> This function returns True if the profiles around the given interface will be
!! statically unstable after mixing above below.  The arguments are the ocean state
!! above and below, including partial calculations from a tridiagonal solver.
function is_unstable(dRho_dT, dRho_dS, h_a, h_b, mix_A, mix_B, T_a, T_b, S_a, S_b, &
                     Te_aa, Te_bb, Se_aa, Se_bb, d_A, d_B)
  real, intent(in) :: drho_dt !< The derivative of in situ density with temperature [R degC-1 ~> kg m-3 degC-1]
  real, intent(in) :: drho_ds !< The derivative of in situ density with salinity [R ppt-1 ~> kg m-3 ppt-1]
  real, intent(in) :: h_a     !< The thickness of the layer above [H ~> m or kg m-2]
  real, intent(in) :: h_b     !< The thickness of the layer below [H ~> m or kg m-2]
  real, intent(in) :: mix_a   !< The time integrated mixing rate of the interface above [H ~> m or kg m-2]
  real, intent(in) :: mix_b   !< The time integrated mixing rate of the interface below [H ~> m or kg m-2]
  real, intent(in) :: t_a     !< The initial temperature of the layer above [degC]
  real, intent(in) :: t_b     !< The initial temperature of the layer below [degC]
  real, intent(in) :: s_a     !< The initial salinity of the layer below [ppt]
  real, intent(in) :: s_b     !< The initial salinity of the layer below [ppt]
  real, intent(in) :: te_aa   !< The estimated temperature two layers above rescaled by d_A [degC]
  real, intent(in) :: te_bb   !< The estimated temperature two layers below rescaled by d_B [degC]
  real, intent(in) :: se_aa   !< The estimated salinity two layers above rescaled by d_A [ppt]
  real, intent(in) :: se_bb   !< The estimated salinity two layers below rescaled by d_B [ppt]
  real, intent(in) :: d_a     !< The rescaling dependency across the interface above, nondim.
  real, intent(in) :: d_b     !< The rescaling dependency across the interface below, nondim.
  logical :: is_unstable !< The return value, true if the profile is statically unstable
                         !! around the interface in question.

  ! These expressions for the local stability are long, but they have been carefully
  ! grouped for accuracy even when the mixing rates are huge or tiny, and common
  ! positive definite factors that would appear in the final expression for the
  ! locally referenced potential density difference across an interface have been omitted.
  is_unstable = (drho_dt * ((h_a * h_b * (t_b - t_a) + &
                             mix_a*mix_b * (d_a*te_bb - d_b*te_aa)) + &
                            (h_a*mix_b * (te_bb - d_b*t_a) + &
                             h_b*mix_a * (d_a*t_b - te_aa)) ) + &
                 drho_ds * ((h_a * h_b * (s_b - s_a) + &
                             mix_a*mix_b * (d_a*se_bb - d_b*se_aa)) + &
                            (h_a*mix_b * (se_bb - d_b*s_a) + &
                             h_b*mix_a * (d_a*s_b - se_aa)) ) < 0.0)
end function is_unstable

!> Returns the partial derivatives of locally referenced potential density with
!! temperature and salinity after the properties have been smoothed with a small
!! constant diffusivity.
subroutine smoothed_drdt_drds(h, tv, Kddt, dR_dT, dR_dS, G, GV, US, j, p_surf, halo)
  type(ocean_grid_type),   intent(in)  :: G    !< The ocean's grid structure
  type(verticalGrid_type), intent(in)  :: GV   !< The ocean's vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)  :: h    !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),   intent(in)  :: tv   !< A structure pointing to various
                                               !! thermodynamic variables
  real,                    intent(in)  :: Kddt !< A diffusivity times a time increment [H2 ~> m2 or kg2 m-4].
  real, dimension(SZI_(G),SZK_(G)+1), &
                           intent(out) :: dR_dT !< Derivative of locally referenced
                                               !! potential density with temperature [R degC-1 ~> kg m-3 degC-1]
  real, dimension(SZI_(G),SZK_(G)+1), &
                           intent(out) :: dR_dS !< Derivative of locally referenced
                                               !! potential density with salinity [R degC-1 ~> kg m-3 ppt-1]
  type(unit_scale_type),   intent(in)  :: US   !< A dimensional unit scaling type
  integer,                 intent(in)  :: j    !< The j-point to work on.
  real, dimension(:,:),    pointer     :: p_surf !< The pressure at the ocean surface [R L2 T-2 ~> Pa].
  integer,       optional, intent(in)  :: halo !< Halo width over which to compute

  ! Local variables
  real :: mix(SZI_(G),SZK_(G)+1)   ! The diffusive mixing length (kappa*dt)/dz
                                   ! between layers within 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 :: T_f(SZI_(G),SZK_(G))     ! Filtered temperatures [degC]
  real :: S_f(SZI_(G),SZK_(G))     ! Filtered salinities [ppt]
  real :: pres(SZI_(G))            ! Interface pressures [R L2 T-2 ~> Pa].
  real :: T_EOS(SZI_(G))           ! Filtered and vertically averaged temperatures [degC]
  real :: S_EOS(SZI_(G))           ! Filtered and vertically averaged salinities [ppt]
  real :: kap_dt_x2                ! The product of 2*kappa*dt [H2 ~> m2 or kg2 m-4].
  real :: h_neglect, h0            ! Negligible thicknesses to allow for zero thicknesses,
                                   ! [H ~> m or kg m-2].
  real :: h_tr                     ! The thickness at tracer points, plus h_neglect [H ~> m or kg m-2].
  integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
  integer :: i, k, is, ie, nz

  if (present(halo)) then
    is = g%isc-halo ; ie = g%iec+halo
  else
    is = g%isc ; ie = g%iec
  endif
  nz = g%ke

  h_neglect = gv%H_subroundoff
  kap_dt_x2 = 2.0*kddt

  if (kddt <= 0.0) then
    do k=1,nz ; do i=is,ie
      t_f(i,k) = tv%T(i,j,k) ; s_f(i,k) = tv%S(i,j,k)
    enddo ; enddo
  else
    h0 = 1.0e-16*sqrt(kddt) + h_neglect
    do i=is,ie
      mix(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 + mix(i,2))
      d1(i) = b1(i) * h(i,j,1)
      t_f(i,1) = (b1(i)*h_tr)*tv%T(i,j,1)
      s_f(i,1) = (b1(i)*h_tr)*tv%S(i,j,1)
    enddo
    do k=2,nz-1 ; do i=is,ie
      mix(i,k+1) = kap_dt_x2 / ((h(i,j,k)+h(i,j,k+1)) + h0)

      c1(i,k) = mix(i,k) * b1(i)
      h_tr = h(i,j,k) + h_neglect
      b1(i) = 1.0 / ((h_tr + d1(i)*mix(i,k)) + mix(i,k+1))
      d1(i) = b1(i) * (h_tr + d1(i)*mix(i,k))
      t_f(i,k) = b1(i) * (h_tr*tv%T(i,j,k) + mix(i,k)*t_f(i,k-1))
      s_f(i,k) = b1(i) * (h_tr*tv%S(i,j,k) + mix(i,k)*s_f(i,k-1))
    enddo ; enddo
    do i=is,ie
      c1(i,nz) = mix(i,nz) * b1(i)
      h_tr = h(i,j,nz) + h_neglect
      b1(i) = 1.0 / (h_tr + d1(i)*mix(i,nz))
      t_f(i,nz) = b1(i) * (h_tr*tv%T(i,j,nz) + mix(i,nz)*t_f(i,nz-1))
      s_f(i,nz) = b1(i) * (h_tr*tv%S(i,j,nz) + mix(i,nz)*s_f(i,nz-1))
    enddo
    do k=nz-1,1,-1 ; do i=is,ie
      t_f(i,k) = t_f(i,k) + c1(i,k+1)*t_f(i,k+1)
      s_f(i,k) = s_f(i,k) + c1(i,k+1)*s_f(i,k+1)
    enddo ; enddo
  endif

  if (associated(p_surf)) then
    do i=is,ie ; pres(i) = p_surf(i,j) ; enddo
  else
    do i=is,ie ; pres(i) = 0.0 ; enddo
  endif
  eosdom(:) = eos_domain(g%HI, halo)
  call calculate_density_derivs(t_f(:,1), s_f(:,1), pres, dr_dt(:,1), dr_ds(:,1), tv%eqn_of_state, eosdom)
  do i=is,ie ; pres(i) = pres(i) + h(i,j,1)*(gv%H_to_RZ*gv%g_Earth) ; enddo
  do k=2,nz
    do i=is,ie
      t_eos(i) = 0.5*(t_f(i,k-1) + t_f(i,k))
      s_eos(i) = 0.5*(s_f(i,k-1) + s_f(i,k))
    enddo
    call calculate_density_derivs(t_eos, s_eos, pres, dr_dt(:,k), dr_ds(:,k), tv%eqn_of_state, eosdom)
    do i=is,ie ; pres(i) = pres(i) + h(i,j,k)*(gv%H_to_RZ*gv%g_Earth) ; enddo
  enddo
  call calculate_density_derivs(t_f(:,nz), s_f(:,nz), pres, dr_dt(:,nz+1), dr_ds(:,nz+1), &
                                tv%eqn_of_state, eosdom)

end subroutine smoothed_drdt_drds

end module mom_full_convection