dense_water_initialization.F90

!> Initialization routines for the dense water formation
!! and overflow experiment.
module dense_water_initialization

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

use mom_ale_sponge,    only : ale_sponge_cs, set_up_ale_sponge_field, initialize_ale_sponge
use mom_dyn_horgrid,   only : dyn_horgrid_type
use mom_eos,           only : eos_type
use mom_error_handler, only : mom_error, fatal
use mom_file_parser,   only : get_param, param_file_type
use mom_grid,          only : ocean_grid_type
use mom_sponge,        only : sponge_cs
use mom_unit_scaling,  only : unit_scale_type
use mom_variables,     only : thermo_var_ptrs
use mom_verticalgrid,  only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

public dense_water_initialize_topography
public dense_water_initialize_ts
public dense_water_initialize_sponges

character(len=40) :: mdl = "dense_water_initialization" !< Module name

real, parameter :: default_sill  = 0.2  !< Default depth of the sill [nondim]
real, parameter :: default_shelf = 0.4  !< Default depth of the shelf [nondim]
real, parameter :: default_mld   = 0.25 !< Default depth of the mixed layer [nondim]

contains

!> Initialize the topography field for the dense water experiment
subroutine dense_water_initialize_topography(D, G, param_file, max_depth)
  type(dyn_horgrid_type),  intent(in)  :: g !< The dynamic horizontal grid type
  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
                           intent(out) :: d !< Ocean bottom depth in the units of depth_max
  type(param_file_type),   intent(in)  :: param_file !< Parameter file structure
  real,                    intent(in)  :: max_depth !< Maximum ocean depth in arbitrary units

  ! Local variables
  real, dimension(5) :: domain_params ! nondimensional widths of all domain sections
  real :: sill_frac, shelf_frac
  integer :: i, j
  real :: x

  call get_param(param_file, mdl, "DENSE_WATER_DOMAIN_PARAMS", domain_params, &
       "Fractional widths of all the domain sections for the dense water experiment.\n"//&
       "As a 5-element vector:\n"//&
       "  - open ocean, the section at maximum depth\n"//&
       "  - downslope, the downward overflow slope\n"//&
       "  - sill separating downslope from upslope\n"//&
       "  - upslope, the upward slope accumulating dense water\n"//&
       "  - the shelf in the dense formation region.", &
       units="nondim", fail_if_missing=.true.)
  call get_param(param_file, mdl, "DENSE_WATER_SILL_DEPTH", sill_frac, &
       "Depth of the sill separating downslope from upslope, as fraction of basin depth.", &
       units="nondim", default=default_sill)
  call get_param(param_file, mdl, "DENSE_WATER_SHELF_DEPTH", shelf_frac, &
       "Depth of the shelf region accumulating dense water for overflow, as fraction of basin depth.", &
       units="nondim", default=default_shelf)

  do i = 2, 5
    ! turn widths into positions
    domain_params(i) = domain_params(i-1) + domain_params(i)
  enddo

  do j = g%jsc,g%jec
    do i = g%isc,g%iec
      ! compute normalised zonal coordinate
      x = (g%geoLonT(i,j) - g%west_lon) / g%len_lon

      if (x <= domain_params(1)) then
        ! open ocean region
        d(i,j) = max_depth
      elseif (x <= domain_params(2)) then
        ! downslope region, linear
        d(i,j) = max_depth - (1.0 - sill_frac) * max_depth * &
             (x - domain_params(1)) / (domain_params(2) - domain_params(1))
      elseif (x <= domain_params(3)) then
        ! sill region
        d(i,j) = sill_frac * max_depth
      elseif (x <= domain_params(4)) then
        ! upslope region
        d(i,j) = sill_frac * max_depth + (shelf_frac - sill_frac) * max_depth * &
             (x - domain_params(3)) / (domain_params(4) - domain_params(3))
      else
        ! shelf region
        d(i,j) = shelf_frac * max_depth
      endif
    enddo
  enddo

end subroutine dense_water_initialize_topography

!> Initialize the temperature and salinity for the dense water experiment
subroutine dense_water_initialize_ts(G, GV, param_file, eqn_of_state, T, S, h, just_read_params)
  type(ocean_grid_type),                     intent(in)  :: g !< Horizontal grid control structure
  type(verticalgrid_type),                   intent(in)  :: gv !< Vertical grid control structure
  type(param_file_type),                     intent(in)  :: param_file !< Parameter file structure
  type(eos_type),                            pointer     :: eqn_of_state !< EOS structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: t !< Output temperature [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s !< Output salinity [ppt]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h !< Layer thicknesses [H ~> m or kg m-2]
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
  ! Local variables
  real :: mld, s_ref, s_range, t_ref
  real :: zi, zmid
  logical :: just_read    ! If true, just read parameters but set nothing.
  integer :: i, j, k, nz

  nz = gv%ke

  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  call get_param(param_file, mdl, "DENSE_WATER_MLD", mld, &
       "Depth of unstratified mixed layer as a fraction of the water column.", &
       units="nondim", default=default_mld, do_not_log=just_read)
  call get_param(param_file, mdl, "S_REF", s_ref, 'Reference salinity', &
                 default=35.0, units='1e-3', do_not_log=just_read)
  call get_param(param_file, mdl,"T_REF", t_ref, 'Reference temperature', units='degC', &
                fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"S_RANGE", s_range, 'Initial salinity range', &
                units='1e-3', default=2.0, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  ! uniform temperature everywhere
  t(:,:,:) = t_ref

  do j = g%jsc,g%jec
    do i = g%isc,g%iec
      zi = 0.
      do k = 1,nz
        ! nondimensional middle of layer
        zmid = zi + 0.5 * h(i,j,k) / (gv%Z_to_H * g%max_depth)

        if (zmid < mld) then
          ! use reference salinity in the mixed layer
          s(i,j,k) = s_ref
        else
          ! linear between bottom of mixed layer and bottom
          s(i,j,k) = s_ref + s_range * (zmid - mld) / (1.0 - mld)
        endif

        zi = zi + h(i,j,k) / (gv%Z_to_H * g%max_depth)
      enddo
    enddo
  enddo
end subroutine dense_water_initialize_ts

!> Initialize the restoring sponges for the dense water experiment
subroutine dense_water_initialize_sponges(G, GV, US, tv, param_file, use_ALE, CSp, ACSp)
  type(ocean_grid_type),   intent(in) :: g !< Horizontal grid control structure
  type(verticalgrid_type), intent(in) :: gv !< Vertical grid control structure
  type(unit_scale_type),   intent(in) :: us !< A dimensional unit scaling type
  type(thermo_var_ptrs),   intent(in) :: tv !< Thermodynamic variables
  type(param_file_type),   intent(in) :: param_file !< Parameter file structure
  logical,                 intent(in) :: use_ale !< ALE flag
  type(sponge_cs),         pointer    :: csp !< Layered sponge control structure pointer
  type(ale_sponge_cs),     pointer    :: acsp !< ALE sponge control structure pointer
  ! Local variables
  real :: west_sponge_time_scale, east_sponge_time_scale ! Sponge timescales [T ~> s]
  real :: west_sponge_width, east_sponge_width

  real, dimension(SZI_(G),SZJ_(G)) :: idamp ! inverse damping timescale [T-1 ~> s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h  ! sponge thicknesses [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: t  ! sponge temperature [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: s  ! sponge salinity [ppt]
  real, dimension(SZK_(GV)+1) :: e0, eta1d ! interface positions for ALE sponge [Z ~> m]

  integer :: i, j, k, nz
  real :: x, zi, zmid, dist
  real :: mld, s_ref, s_range, s_dense, t_ref, sill_height

  nz = gv%ke

  call get_param(param_file, mdl, "DENSE_WATER_WEST_SPONGE_TIME_SCALE", west_sponge_time_scale, &
       "The time scale on the west (outflow) of the domain for restoring. If zero, the sponge is disabled.", &
       units="s", default=0., scale=us%s_to_T)
  call get_param(param_file, mdl, "DENSE_WATER_WEST_SPONGE_WIDTH", west_sponge_width, &
       "The fraction of the domain in which the western (outflow) sponge is active.", &
       units="nondim", default=0.1)
  call get_param(param_file, mdl, "DENSE_WATER_EAST_SPONGE_TIME_SCALE", east_sponge_time_scale, &
       "The time scale on the east (outflow) of the domain for restoring. If zero, the sponge is disabled.", &
       units="s", default=0., scale=us%s_to_T)
  call get_param(param_file, mdl, "DENSE_WATER_EAST_SPONGE_WIDTH", east_sponge_width, &
       "The fraction of the domain in which the eastern (outflow) sponge is active.", &
       units="nondim", default=0.1)

  call get_param(param_file, mdl, "DENSE_WATER_EAST_SPONGE_SALT", s_dense, &
       "Salt anomaly of the dense water being formed in the overflow region.", &
       units="1e-3", default=4.0)

  call get_param(param_file, mdl, "DENSE_WATER_MLD", mld, default=default_mld, do_not_log=.true.)
  call get_param(param_file, mdl, "DENSE_WATER_SILL_HEIGHT", sill_height, default=default_sill, do_not_log=.true.)

  call get_param(param_file, mdl, "S_REF", s_ref, default=35.0, do_not_log=.true.)
  call get_param(param_file, mdl, "S_RANGE", s_range, do_not_log=.true.)
  call get_param(param_file, mdl, "T_REF", t_ref, do_not_log=.true.)

  ! no active sponges
  if (west_sponge_time_scale <= 0. .and. east_sponge_time_scale <= 0.) return

  ! everywhere is initially unsponged
  idamp(:,:) = 0.0

  do j = g%jsc, g%jec
    do i = g%isc,g%iec
      if (g%mask2dT(i,j) > 0.) then
        ! nondimensional x position
        x = (g%geoLonT(i,j) - g%west_lon) / g%len_lon

        if (west_sponge_time_scale > 0. .and. x < west_sponge_width) then
          dist = 1. - x / west_sponge_width
          ! scale restoring by depth into sponge
          idamp(i,j) = 1. / west_sponge_time_scale * max(0., min(1., dist))
        elseif (east_sponge_time_scale > 0. .and. x > (1. - east_sponge_width)) then
          dist = 1. - (1. - x) / east_sponge_width
          idamp(i,j) = 1. / east_sponge_time_scale * max(0., min(1., dist))
        endif
      endif
    enddo
  enddo

  if (use_ale) then
    ! construct a uniform grid for the sponge
    do k = 1,nz
      e0(k) = -g%max_depth * (real(k - 1) / real(nz))
    enddo
    e0(nz+1) = -g%max_depth

    do j = g%jsc,g%jec
      do i = g%isc,g%iec
        eta1d(nz+1) = -g%bathyT(i,j)
        do k = nz,1,-1
          eta1d(k) = e0(k)

          if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
            ! is this layer vanished?
            eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
            h(i,j,k) = gv%Angstrom_H
          else
            h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
          endif
        enddo
      enddo
    enddo

    call initialize_ale_sponge(idamp, g, param_file, acsp, h, nz)

    ! construct temperature and salinity for the sponge
    ! start with initial condition
    t(:,:,:) = t_ref
    s(:,:,:) = s_ref

    do j = g%jsc,g%jec
      do i = g%isc,g%iec
        zi = 0.
        x = (g%geoLonT(i,j) - g%west_lon) / g%len_lon
        do k = 1,nz
          ! nondimensional middle of layer
          zmid = zi + 0.5 * h(i,j,k) / (gv%Z_to_H * g%max_depth)

          if (x > (1. - east_sponge_width)) then
            !if (zmid >= 0.9 * sill_height) &
                 s(i,j,k) = s_ref + s_dense
          else
            ! linear between bottom of mixed layer and bottom
            if (zmid >= mld) &
                 s(i,j,k) = s_ref + s_range * (zmid - mld) / (1.0 - mld)
          endif

          zi = zi + h(i,j,k) / (gv%Z_to_H * g%max_depth)
        enddo
      enddo
    enddo

    if (associated(tv%T)) call set_up_ale_sponge_field(t, g, tv%T, acsp)
    if (associated(tv%S)) call set_up_ale_sponge_field(s, g, tv%S, acsp)
  else
    call mom_error(fatal, "dense_water_initialize_sponges: trying to use non ALE sponge")
  endif
end subroutine dense_water_initialize_sponges

end module dense_water_initialization

!> \namespace dense_water_initialization
!!
!! This experiment consists of a shelf accumulating dense water, which spills
!! over an upward slope and a sill, before flowing down a slope into an open
!! ocean region. It's intended as a test of one of the motivating situations
!! for the adaptive coordinate.
!!
!! The nondimensional widths of the 5 regions are controlled by the
!! <code>DENSE_WATER_DOMAIN_PARAMS</code>, and the heights of the sill and shelf
!! as a fraction of the total domain depth are controlled by
!! <code>DENSE_WATER_SILL_HEIGHT</code> and <code>DENSE_WATER_SHELF_HEIGHT</code>.
!!
!! The density in the domain is governed by a linear equation of state, and
!! is set up with a mixed layer of non-dimensional depth <code>DENSE_WATER_MLD</code>
!! below which there is a linear stratification from <code>S_REF</code>, increasing
!! by <code>S_RANGE</code> to the bottom.
!!
!! To force the experiment, there are sponges on the inflow and outflow of the
!! domain. The inflow sponge has a salinity anomaly of
!! <code>DENSE_WATER_EAST_SPONGE_SALT</code> through the entire depth. The outflow
!! sponge simply restores to the initial condition. Both sponges have controllable
!! widths and restoring timescales.