coord_adapt.F90

!> Regrid columns for the adaptive coordinate
module coord_adapt

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

use mom_eos,           only : calculate_density_derivs
use mom_error_handler, only : mom_error, fatal
use mom_unit_scaling,  only : unit_scale_type
use mom_variables,     only : ocean_grid_type, thermo_var_ptrs
use mom_verticalgrid,  only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

!> Control structure for adaptive coordinates (coord_adapt).
type, public :: adapt_cs ; private

  !> Number of layers/levels
  integer :: nk

  !> Nominal near-surface resolution [H ~> m or kg m-2]
  real, allocatable, dimension(:) :: coordinateresolution

  !> Ratio of optimisation and diffusion timescales
  real :: adapttimeratio

  !> Nondimensional coefficient determining how much optimisation to apply
  real :: adaptalpha

  !> Near-surface zooming depth [H ~> m or kg m-2]
  real :: adaptzoom

  !> Near-surface zooming coefficient
  real :: adaptzoomcoeff

  !> Stratification-dependent diffusion coefficient
  real :: adaptbuoycoeff

  !> Reference density difference for stratification-dependent diffusion [R ~> kg m-3]
  real :: adaptdrho0

  !> If true, form a HYCOM1-like mixed layet by preventing interfaces
  !! from becoming shallower than the depths set by coordinateResolution
  logical :: adaptdomin  = .false.
end type adapt_cs

public init_coord_adapt, set_adapt_params, build_adapt_column, end_coord_adapt

contains

!> Initialise an adapt_CS with parameters
subroutine init_coord_adapt(CS, nk, coordinateResolution, m_to_H, kg_m3_to_R)
  type(adapt_cs),     pointer    :: cs !< Unassociated pointer to hold the control structure
  integer,            intent(in) :: nk !< Number of layers in the grid
  real, dimension(:), intent(in) :: coordinateresolution !< Nominal near-surface resolution [m] or
                                       !! other units specified with m_to_H
  real,               intent(in) :: m_to_h !< A conversion factor from m to the units of thicknesses
  real,               intent(in) :: kg_m3_to_r !< A conversion factor from kg m-3 to the units of density

  if (associated(cs)) call mom_error(fatal, "init_coord_adapt: CS already associated")
  allocate(cs)
  allocate(cs%coordinateResolution(nk))

  cs%nk = nk
  cs%coordinateResolution(:) = coordinateresolution(:)

  ! Set real parameter default values
  cs%adaptTimeRatio = 1e-1 ! Nondim.
  cs%adaptAlpha     = 1.0  ! Nondim.
  cs%adaptZoom      = 200.0 * m_to_h    ! [H ~> m or kg m-2]
  cs%adaptZoomCoeff = 0.0  ! Nondim.
  cs%adaptBuoyCoeff = 0.0  ! Nondim.
  cs%adaptDrho0     = 0.5 * kg_m3_to_r  ! [R ~> kg m-3]

end subroutine init_coord_adapt

!> Clean up the coordinate control structure
subroutine end_coord_adapt(CS)
  type(adapt_cs), pointer :: cs  !< The control structure for this module

  ! nothing to do
  if (.not. associated(cs)) return
  deallocate(cs%coordinateResolution)
  deallocate(cs)
end subroutine end_coord_adapt

!> This subtroutine can be used to set the parameters for coord_adapt module
subroutine set_adapt_params(CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoomCoeff, &
                            adaptBuoyCoeff, adaptDrho0, adaptDoMin)
  type(adapt_cs),    pointer    :: cs  !< The control structure for this module
  real,    optional, intent(in) :: adapttimeratio !< Ratio of optimisation and diffusion timescales
  real,    optional, intent(in) :: adaptalpha     !< Nondimensional coefficient determining
                                                  !! how much optimisation to apply
  real,    optional, intent(in) :: adaptzoom      !< Near-surface zooming depth [H ~> m or kg m-2]
  real,    optional, intent(in) :: adaptzoomcoeff !< Near-surface zooming coefficient
  real,    optional, intent(in) :: adaptbuoycoeff !< Stratification-dependent diffusion coefficient
  real,    optional, intent(in) :: adaptdrho0  !< Reference density difference for
                                               !! stratification-dependent diffusion [R ~> kg m-3]
  logical, optional, intent(in) :: adaptdomin  !< If true, form a HYCOM1-like mixed layer by
                                               !! preventing interfaces from becoming shallower than
                                               !! the depths set by coordinateResolution

  if (.not. associated(cs)) call mom_error(fatal, "set_adapt_params: CS not associated")

  if (present(adapttimeratio)) cs%adaptTimeRatio = adapttimeratio
  if (present(adaptalpha)) cs%adaptAlpha = adaptalpha
  if (present(adaptzoom)) cs%adaptZoom = adaptzoom
  if (present(adaptzoomcoeff)) cs%adaptZoomCoeff = adaptzoomcoeff
  if (present(adaptbuoycoeff)) cs%adaptBuoyCoeff = adaptbuoycoeff
  if (present(adaptdrho0)) cs%adaptDrho0 = adaptdrho0
  if (present(adaptdomin)) cs%adaptDoMin = adaptdomin
end subroutine set_adapt_params

subroutine build_adapt_column(CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, zNext)
  type(adapt_cs),                              intent(in)    :: cs   !< The control structure for this module
  type(ocean_grid_type),                       intent(in)    :: g    !< The ocean's grid structure
  type(verticalgrid_type),                     intent(in)    :: gv   !< The ocean's vertical grid structure
  type(unit_scale_type),                       intent(in)    :: us   !< A dimensional unit scaling type
  type(thermo_var_ptrs),                       intent(in)    :: tv   !< A structure pointing to various
                                                                     !! thermodynamic variables
  integer,                                     intent(in)    :: i    !< The i-index of the column to work on
  integer,                                     intent(in)    :: j    !< The j-index of the column to work on
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: zint !< Interface heights [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: tint !< Interface temperatures [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: sint !< Interface salinities [ppt]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2]
  real, dimension(SZK_(GV)+1),                 intent(inout) :: znext !< updated interface positions

  ! Local variables
  integer :: k, nz
  real :: h_up, b1, b_denom_1, d1, depth, nominal_z, stretching
  real :: drdz  ! The vertical density gradient [R H-1 ~> kg m-4 or m-1]
  real, dimension(SZK_(GV)+1) :: alpha ! drho/dT [R degC-1 ~> kg m-3 degC-1]
  real, dimension(SZK_(GV)+1) :: beta  ! drho/dS [R ppt-1 ~> kg m-3 ppt-1]
  real, dimension(SZK_(GV)+1) :: del2sigma ! Laplacian of in situ density times grid spacing [R ~> kg m-3]
  real, dimension(SZK_(GV)+1) :: dh_d2s ! Thickness change in response to del2sigma [H ~> m or kg m-2]
  real, dimension(SZK_(GV)) :: kgrid, c1 ! grid diffusivity on layers, and tridiagonal work array

  nz = cs%nk

  ! set bottom and surface of zNext
  znext(1) = 0.
  znext(nz+1) = zint(i,j,nz+1)

  ! local depth for scaling diffusivity
  depth = g%bathyT(i,j) * gv%Z_to_H

  ! initialize del2sigma and the thickness change response to it zero
  del2sigma(:) = 0.0 ; dh_d2s(:) = 0.0

  ! calculate del-squared of neutral density by a
  ! stencilled finite difference
  ! TODO: this needs to be adjusted to account for vanished layers near topography

  ! up (j-1)
  if (g%mask2dT(i,j-1) > 0.) then
    call calculate_density_derivs( &
         0.5 * (tint(i,j,2:nz) + tint(i,j-1,2:nz)), &
         0.5 * (sint(i,j,2:nz) + sint(i,j-1,2:nz)), &
         0.5 * (zint(i,j,2:nz) + zint(i,j-1,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
         alpha, beta, tv%eqn_of_state, (/2,nz/) )

    del2sigma(2:nz) = del2sigma(2:nz) + &
         (alpha(2:nz) * (tint(i,j-1,2:nz) - tint(i,j,2:nz)) + &
          beta(2:nz)  * (sint(i,j-1,2:nz) - sint(i,j,2:nz)))
  endif
  ! down (j+1)
  if (g%mask2dT(i,j+1) > 0.) then
    call calculate_density_derivs( &
         0.5 * (tint(i,j,2:nz) + tint(i,j+1,2:nz)), &
         0.5 * (sint(i,j,2:nz) + sint(i,j+1,2:nz)), &
         0.5 * (zint(i,j,2:nz) + zint(i,j+1,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
         alpha, beta, tv%eqn_of_state, (/2,nz/) )

    del2sigma(2:nz) = del2sigma(2:nz) + &
         (alpha(2:nz) * (tint(i,j+1,2:nz) - tint(i,j,2:nz)) + &
          beta(2:nz)  * (sint(i,j+1,2:nz) - sint(i,j,2:nz)))
  endif
  ! left (i-1)
  if (g%mask2dT(i-1,j) > 0.) then
    call calculate_density_derivs( &
         0.5 * (tint(i,j,2:nz) + tint(i-1,j,2:nz)), &
         0.5 * (sint(i,j,2:nz) + sint(i-1,j,2:nz)), &
         0.5 * (zint(i,j,2:nz) + zint(i-1,j,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
         alpha, beta, tv%eqn_of_state, (/2,nz/) )

    del2sigma(2:nz) = del2sigma(2:nz) + &
         (alpha(2:nz) * (tint(i-1,j,2:nz) - tint(i,j,2:nz)) + &
          beta(2:nz)  * (sint(i-1,j,2:nz) - sint(i,j,2:nz)))
  endif
  ! right (i+1)
  if (g%mask2dT(i+1,j) > 0.) then
    call calculate_density_derivs( &
         0.5 * (tint(i,j,2:nz) + tint(i+1,j,2:nz)), &
         0.5 * (sint(i,j,2:nz) + sint(i+1,j,2:nz)), &
         0.5 * (zint(i,j,2:nz) + zint(i+1,j,2:nz)) * (gv%H_to_RZ * gv%g_Earth), &
         alpha, beta, tv%eqn_of_state, (/2,nz/) )

    del2sigma(2:nz) = del2sigma(2:nz) + &
         (alpha(2:nz) * (tint(i+1,j,2:nz) - tint(i,j,2:nz)) + &
          beta(2:nz)  * (sint(i+1,j,2:nz) - sint(i,j,2:nz)))
  endif

  ! at this point, del2sigma contains the local neutral density curvature at
  ! h-points, on interfaces
  ! we need to divide by drho/dz to give an interfacial displacement
  !
  ! a positive curvature means we're too light relative to adjacent columns,
  ! so del2sigma needs to be positive too (push the interface deeper)
  call calculate_density_derivs(tint(i,j,:), sint(i,j,:), zint(i,j,:) * (gv%H_to_RZ * gv%g_Earth), &
       alpha, beta, tv%eqn_of_state, (/1,nz+1/) )
  do k = 2, nz
    ! TODO make lower bound here configurable
    dh_d2s(k) = del2sigma(k) * (0.5 * (h(i,j,k-1) + h(i,j,k))) / &
         max(alpha(k) * (tv%T(i,j,k) - tv%T(i,j,k-1)) + &
             beta(k)  * (tv%S(i,j,k) - tv%S(i,j,k-1)), 1e-20*us%kg_m3_to_R)

    ! don't move the interface so far that it would tangle with another
    ! interface in the direction we're moving (or exceed a Nyquist limit
    ! that could cause oscillations of the interface)
    h_up = merge(h(i,j,k), h(i,j,k-1), dh_d2s(k) > 0.)
    dh_d2s(k) = 0.5 * cs%adaptAlpha * &
         sign(min(abs(del2sigma(k)), 0.5 * h_up), dh_d2s(k))

    ! update interface positions so we can diffuse them
    znext(k) = zint(i,j,k) + dh_d2s(k)
  enddo

  ! solve diffusivity equation to smooth grid
  ! upper diagonal coefficients: -kGrid(2:nz)
  ! lower diagonal coefficients: -kGrid(1:nz-1)
  ! diagonal coefficients:       1 + (kGrid(1:nz-1) + kGrid(2:nz))
  !
  ! first, calculate the diffusivities within layers
  do k = 1, nz
    ! calculate the dr bit of drdz
    drdz = 0.5 * (alpha(k) + alpha(k+1)) * (tint(i,j,k+1) - tint(i,j,k)) + &
           0.5 * (beta(k)  + beta(k+1))  * (sint(i,j,k+1) - sint(i,j,k))
    ! divide by dz from the new interface positions
    drdz = drdz / (znext(k) - znext(k+1) + gv%H_subroundoff)
    ! don't do weird stuff in unstably-stratified regions
    drdz = max(drdz, 0.)

    ! set vertical grid diffusivity
    kgrid(k) = (cs%adaptTimeRatio * nz**2 * depth) * &
         (cs%adaptZoomCoeff / (cs%adaptZoom + 0.5*(znext(k) + znext(k+1))) + &
         (cs%adaptBuoyCoeff * drdz / cs%adaptDrho0) + &
         max(1.0 - cs%adaptZoomCoeff - cs%adaptBuoyCoeff, 0.0) / depth)
  enddo

  ! initial denominator (first diagonal element)
  b1 = 1.0
  ! initial Q_1 = 1 - q_1 = 1 - 0/1
  d1 = 1.0
  ! work on all interior interfaces
  do k = 2, nz
    ! calculate numerator of Q_k
    b_denom_1 = 1. + d1 * kgrid(k-1)
    ! update denominator for k
    b1 = 1.0 / (b_denom_1 + kgrid(k))

    ! calculate q_k
    c1(k) = kgrid(k) * b1
    ! update Q_k = 1 - q_k
    d1 = b_denom_1 * b1

    ! update RHS
    znext(k) = b1 * (znext(k) + kgrid(k-1)*znext(k-1))
  enddo
  ! final substitution
  do k = nz, 2, -1
    znext(k) = znext(k) + c1(k)*znext(k+1)
  enddo

  if (cs%adaptDoMin) then
    nominal_z = 0.
    stretching = zint(i,j,nz+1) / depth

    do k = 2, nz+1
      nominal_z = nominal_z + cs%coordinateResolution(k-1) * stretching
      ! take the deeper of the calculated and nominal positions
      znext(k) = max(znext(k), nominal_z)
      ! interface can't go below topography
      znext(k) = min(znext(k), zint(i,j,nz+1))
    enddo
  endif
end subroutine build_adapt_column

end module coord_adapt