coord_adapt.F90

Go to the documentation of this file.

!> Regrid columns for the adaptive coordinate
module coord_adapt

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

implicit none ; private

#include <MOM_memory.h>

type, public :: adapt_cs
  private

  !> Number of layers/levels
  integer :: nk

  !> Nominal near-surface resolution
  real, allocatable, dimension(:) :: coordinateresolution

  !> Ratio of optimisation and diffusion timescales
  real :: adapttimeratio = 1e-1

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

  !> Near-surface zooming depth
  real :: adaptzoom      = 200.0

  !> Near-surface zooming coefficient
  real :: adaptzoomcoeff = 0.0

  !> Stratification-dependent diffusion coefficient
  real :: adaptbuoycoeff = 0.0

  !> Reference density difference for stratification-dependent diffusion
  real :: adaptdrho0     = 0.5

  !> 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)
  type(adapt_cs),     pointer    :: CS !< Unassociated pointer to hold the control structure
  integer,            intent(in) :: nk
  real, dimension(:), intent(in) :: coordinateResolution

  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(:)
end subroutine init_coord_adapt

!> Clean up the coordinate control structure
subroutine end_coord_adapt(CS)
  type(adapt_cs), pointer :: CS

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

subroutine set_adapt_params(CS, adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoomCoeff, &
     adaptBuoyCoeff, adaptDrho0, adaptDoMin)
  type(adapt_cs),    pointer    :: CS
  real, optional,    intent(in) :: adaptTimeRatio, adaptAlpha, adaptZoom, adaptZoomCoeff
  real, optional,    intent(in) :: adaptBuoyCoeff, adaptDrho0
  logical, optional, intent(in) :: adaptDoMin

  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, tv, i, j, zInt, tInt, sInt, h, zNext)
  type(adapt_cs),                              intent(in)    :: CS
  type(ocean_grid_type),                       intent(in)    :: G    !< The ocean's grid structure
  type(verticalgrid_type),                     intent(in)    :: GV   !< The ocean's vertical grid structure
  type(thermo_var_ptrs),                       intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
  integer,                                     intent(in)    :: i, j
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in)    :: zInt, tInt, sInt
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),   intent(in)    :: h    !< Layer thicknesses, in H (usually 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, drdz, nominal_z, stretching
  real, dimension(SZK_(GV)+1) :: alpha, beta, del2sigma ! drho/dT and drho/dS
  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%m_to_H

  ! initialize del2sigma to zero
  del2sigma(:) = 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_Pa, &
         alpha, beta, 2, nz - 1, tv%eqn_of_state)

    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_Pa, &
         alpha, beta, 2, nz - 1, tv%eqn_of_state)

    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_Pa, &
         alpha, beta, 2, nz - 1, tv%eqn_of_state)

    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_Pa, &
         alpha, beta, 2, nz - 1, tv%eqn_of_state)

    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_Pa, &
       alpha, beta, 1, nz + 1, tv%eqn_of_state)
  do k = 2, nz
    ! TODO make lower bound here configurable
    del2sigma(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)

    ! 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), del2sigma(k) > 0.)
    del2sigma(k) = 0.5 * cs%adaptAlpha * &
         sign(min(abs(del2sigma(k)), 0.5 * h_up), del2sigma(k))

    ! update interface positions so we can diffuse them
    znext(k) = zint(i,j,k) + del2sigma(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 * gv%m_to_H + 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