dumbbell_initialization.F90

!> Configures the model for the idealized dumbbell test case.
module dumbbell_initialization

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

use mom_domains, only : sum_across_pes
use mom_dyn_horgrid, only : dyn_horgrid_type
use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_get_input, only : directories
use mom_grid, only : ocean_grid_type
use mom_sponge, only : set_up_sponge_field, initialize_sponge, sponge_cs
use mom_tracer_registry, only : tracer_registry_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, calculate_density_derivs, eos_type
use regrid_consts, only : coordinatemode, default_coordinate_mode
use regrid_consts, only : regridding_layer, regridding_zstar
use regrid_consts, only : regridding_rho, regridding_sigma
use mom_ale_sponge,    only : ale_sponge_cs, set_up_ale_sponge_field, initialize_ale_sponge

implicit none ; private

#include <MOM_memory.h>

character(len=40) :: mdl = "dumbbell_initialization" !< This module's name.

public dumbbell_initialize_topography
public dumbbell_initialize_thickness
public dumbbell_initialize_temperature_salinity
public dumbbell_initialize_sponges

! 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

!> Initialization of topography.
subroutine dumbbell_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
  integer   :: i, j
  real      :: x, y, delta, dblen, dbfrac
  logical   :: dbrotate

  call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, &
                'Lateral Length scale for dumbbell.',&
                 units='k', default=600., do_not_log=.false.)
  call get_param(param_file, mdl,"DUMBBELL_FRACTION",dbfrac, &
                'Meridional fraction for narrow part of dumbbell.',&
                 units='nondim', default=0.5, do_not_log=.false.)
  call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
                'Logical for rotation of dumbbell domain.',&
                 units='nondim', default=.false., do_not_log=.false.)

  if (g%x_axis_units == 'm') then
    dblen=dblen*1.e3
  endif

  if (dbrotate) then
    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      ! Compute normalized zonal coordinates (x,y=0 at center of domain)
      x = ( g%geoLonT(i,j) ) / g%len_lon
      y = ( g%geoLatT(i,j)  ) / dblen
      d(i,j) = g%max_depth
      if ((y>=-0.25 .and. y<=0.25) .and. (x <= -0.5*dbfrac .or. x >= 0.5*dbfrac)) then
        d(i,j) = 0.0
      endif
    enddo ; enddo
  else
    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      ! Compute normalized zonal coordinates (x,y=0 at center of domain)
      x = ( g%geoLonT(i,j) ) / dblen
      y = ( g%geoLatT(i,j)  ) / g%len_lat
      d(i,j) = g%max_depth
      if ((x>=-0.25 .and. x<=0.25) .and. (y <= -0.5*dbfrac .or. y >= 0.5*dbfrac)) then
        d(i,j) = 0.0
      endif
    enddo ; enddo
  endif

end subroutine dumbbell_initialize_topography

!> Initializes the layer thicknesses to be uniform in the dumbbell test case
subroutine dumbbell_initialize_thickness ( h, G, GV, US, param_file, just_read_params)
  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_(GV)), &
                           intent(out) :: h           !< The thickness that is being initialized [H ~> m or kg m-2].
  type(param_file_type),   intent(in)  :: param_file  !< A structure indicating the open file
                                                      !! to parse for model parameter values.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

  real :: e0(szk_(g)+1)   ! The resting interface heights [Z ~> m], usually
                          ! negative because it is positive upward.
  real :: eta1d(szk_(g)+1)! Interface height relative to the sea surface
                          ! positive upward [Z ~> m].
  real :: min_thickness   ! The minimum layer thicknesses [Z ~> m].
  real :: s_surf, s_range, s_ref, s_light, s_dense ! Various salinities [ppt].
  real :: eta_ic_quanta   ! The granularity of quantization of intial interface heights [Z-1 ~> m-1].
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=20) :: verticalcoordinate
  logical :: just_read    ! If true, just read parameters but set nothing.
  integer :: i, j, k, is, ie, js, je, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

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

  if (.not.just_read) &
    call mom_mesg("MOM_initialization.F90, initialize_thickness_uniform: setting thickness")

  if (.not.just_read) call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, &
                'Minimum thickness for layer',&
                 units='m', default=1.0e-3, do_not_log=just_read, scale=us%m_to_Z)
  call get_param(param_file, mdl,"REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
                 default=default_coordinate_mode, do_not_log=just_read)

  ! WARNING: this routine specifies the interface heights so that the last layer
  !          is vanished, even at maximum depth. In order to have a uniform
  !          layer distribution, use this line of code within the loop:
  !          e0(k) = -G%max_depth * real(k-1) / real(nz)
  !          To obtain a thickness distribution where the last layer is
  !          vanished and the other thicknesses uniformly distributed, use:
  !          e0(k) = -G%max_depth * real(k-1) / real(nz-1)
  !do k=1,nz+1
  !  e0(k) = -G%max_depth * real(k-1) / real(nz)
  !enddo

  select case ( coordinatemode(verticalcoordinate) )

  case ( regridding_layer, regridding_rho ) ! Initial thicknesses for isopycnal coordinates
    call get_param(param_file, mdl,"INITIAL_SSS", s_surf, default=34., do_not_log=.true.)
    call get_param(param_file, mdl,"INITIAL_S_RANGE", s_range, default=2., 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, "TS_RANGE_S_LIGHT", s_light, default = s_ref, do_not_log=.true.)
    call get_param(param_file, mdl, "TS_RANGE_S_DENSE", s_dense, default = s_ref, do_not_log=.true.)
    call get_param(param_file, mdl, "INTERFACE_IC_QUANTA", eta_ic_quanta, &
                   "The granularity of initial interface height values "//&
                   "per meter, to avoid sensivity to order-of-arithmetic changes.", &
                   default=2048.0, units="m-1", scale=us%Z_to_m, do_not_log=just_read)
    if (just_read) return ! All run-time parameters have been read, so return.

    do k=1,nz+1
      ! Salinity of layer k is S_light + (k-1)/(nz-1) * (S_dense - S_light)
      ! Salinity of interface K is S_light + (K-3/2)/(nz-1) * (S_dense - S_light)
      ! Salinity at depth z should be S(z) = S_surf - S_range * z/max_depth
      ! Equating: S_surf - S_range * z/max_depth = S_light + (K-3/2)/(nz-1) * (S_dense - S_light)
      ! Equating: - S_range * z/max_depth = S_light - S_surf + (K-3/2)/(nz-1) * (S_dense - S_light)
      ! Equating: z/max_depth = - ( S_light - S_surf + (K-3/2)/(nz-1) * (S_dense - S_light) ) / S_range
      e0(k) = - g%max_depth * ( ( s_light  - s_surf ) + ( s_dense - s_light ) * &
                ( (real(k)-1.5) / real(nz-1) ) ) / s_range
      ! Force round numbers ... the above expression has irrational factors ...
      if (eta_ic_quanta > 0.0) &
        e0(k) = nint(eta_ic_quanta*e0(k)) / eta_ic_quanta
      e0(k) = min(real(1-k)*gv%Angstrom_Z, e0(k)) ! Bound by surface
      e0(k) = max(-g%max_depth, e0(k)) ! Bound by bottom
    enddo
    do j=js,je ; do i=is,ie
      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
          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

  case ( regridding_zstar )                       ! Initial thicknesses for z coordinates
    if (just_read) return ! All run-time parameters have been read, so return.
    do j=js,je ; do i=is,ie
      eta1d(nz+1) = -g%bathyT(i,j)
      do k=nz,1,-1
        eta1d(k) = -g%max_depth * real(k-1) / real(nz)
        if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
          eta1d(k) = eta1d(k+1) + min_thickness
          h(i,j,k) = gv%Z_to_H * min_thickness
        else
          h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
        endif
      enddo
    enddo ; enddo

  case ( regridding_sigma )             ! Initial thicknesses for sigma coordinates
    if (just_read) return ! All run-time parameters have been read, so return.
    do j=js,je ; do i=is,ie
      h(i,j,:) = gv%Z_to_H * g%bathyT(i,j) / dfloat(nz)
    enddo ; enddo

end select

end subroutine dumbbell_initialize_thickness

!> Initial values for temperature and salinity for the dumbbell test case
subroutine dumbbell_initialize_temperature_salinity ( T, S, h, G, GV, param_file, &
                                                  eqn_of_state, just_read_params)
  type(ocean_grid_type),                     intent(in)  :: g !< Ocean grid structure
  type(verticalgrid_type),                   intent(in) :: gv !< Vertical grid structure
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: t !< Potential temperature [degC]
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: s !< Salinity [ppt]
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in)  :: h !< Layer thickness [H ~> m or kg m-2]
  type(param_file_type),                     intent(in)  :: param_file !< Parameter file structure
  type(eos_type),                            pointer     :: eqn_of_state !< Equation of state structure
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.

  ! Local variables
  integer :: i, j, k, is, ie, js, je, nz, k_light
  real    :: xi0, xi1, dxi, r, s_surf, t_surf, s_range, t_range
  real    :: x, y, dblen
  real    :: t_ref, t_light, t_dense, s_ref, s_light, s_dense, a1, frac_dense, k_frac, res_rat
  logical :: just_read    ! If true, just read parameters but set nothing.
  logical :: dbrotate     ! If true, rotate the domain.
  character(len=20) :: verticalcoordinate, density_profile

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

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

  t_surf = 20.0

  call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
                 default=default_coordinate_mode, do_not_log=just_read)
  call get_param(param_file, mdl,"INITIAL_DENSITY_PROFILE", density_profile, &
                 'Initial profile shape. Valid values are "linear", "parabolic" '// &
                 'and "exponential".', default='linear', do_not_log=just_read)
  call get_param(param_file, mdl,"DUMBBELL_SREF", s_surf, &
                 'DUMBBELL REFERENCE SALINITY', units='1e-3', default=34., do_not_log=just_read)
  call get_param(param_file, mdl,"DUMBBELL_S_RANGE", s_range, &
                 'DUMBBELL salinity range (right-left)', units='1e-3', &
                 default=2., do_not_log=just_read)
  call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, &
                'Lateral Length scale for dumbbell ',&
                 units='k', default=600., do_not_log=just_read)
  call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
                'Logical for rotation of dumbbell domain.',&
                 units='nondim', default=.false., do_not_log=just_read)

  if (g%x_axis_units == 'm') then
    dblen=dblen*1.e3
  endif

  do j=g%jsc,g%jec
    do i=g%isc,g%iec
    ! Compute normalized zonal coordinates (x,y=0 at center of domain)
      if (dbrotate) then
        ! This is really y in the rotated case
        x = ( g%geoLatT(i,j) ) / dblen
      else
        x = ( g%geoLonT(i,j) ) / dblen
      endif
      do k=1,nz
        t(i,j,k)=t_surf
      enddo
      if (x>=0. ) then
        do k=1,nz
          s(i,j,k)=s_surf + 0.5*s_range
        enddo
      endif
      if (x<0. ) then
        do k=1,nz
          s(i,j,k)=s_surf - 0.5*s_range
        enddo
      endif

    enddo
  enddo

end subroutine dumbbell_initialize_temperature_salinity

!> Initialize the restoring sponges for the dumbbell test case
subroutine dumbbell_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

  real :: sponge_time_scale  ! The damping time scale [T ~> s]

  real, dimension(SZI_(G),SZJ_(G)) :: idamp ! inverse damping timescale [T-1 ~> s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h, t, s ! sponge thicknesses, temp and salt
  real, dimension(SZK_(GV)+1) :: e0, eta1d ! interface positions for ALE sponge

  integer :: i, j, k, nz
  real :: x, zi, zmid, dist, min_thickness, dblen
  real :: mld, s_ref, s_range, s_dense, t_ref, sill_height
  logical :: dbrotate    ! If true, rotate the domain.

  call get_param(param_file, mdl,"DUMBBELL_LEN",dblen, &
                'Lateral Length scale for dumbbell ',&
                 units='k', default=600., do_not_log=.true.)
  call get_param(param_file, mdl, "DUMBBELL_ROTATION", dbrotate, &
                'Logical for rotation of dumbbell domain.',&
                 units='nondim', default=.false., do_not_log=.true.)

  if (g%x_axis_units == 'm') then
    dblen=dblen*1.e3
  endif

  nz = gv%ke

  call get_param(param_file, mdl, "DUMBBELL_SPONGE_TIME_SCALE", sponge_time_scale, &
       "The time scale in the reservoir for restoring. If zero, the sponge is disabled.", &
       units="s", default=0., scale=us%s_to_T)
  call get_param(param_file, mdl, "DUMBBELL_SREF", s_ref, do_not_log=.true.)
  call get_param(param_file, mdl, "DUMBBELL_S_RANGE", s_range, do_not_log=.true.)
  call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, &
                'Minimum thickness for layer',&
                 units='m', default=1.0e-3, do_not_log=.true., scale=us%m_to_Z)

  ! no active sponges
  if (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
        if (dbrotate) then
          ! This is really y in the rotated case
          x = ( g%geoLatT(i,j) ) / dblen
        else
          x = ( g%geoLonT(i,j) ) / dblen
        endif
        if (x > 0.25 .or. x < -0.25) then
          ! scale restoring by depth into sponge
          idamp(i,j) = 1. / sponge_time_scale
        endif
      endif
    enddo
  enddo

  if (use_ale) then
    ! construct a uniform grid for the sponge
    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) = -g%max_depth * real(k-1) / real(nz)
        if (eta1d(k) < (eta1d(k+1) + min_thickness)) then
          eta1d(k) = eta1d(k+1) + min_thickness
          h(i,j,k) = gv%Z_to_H * min_thickness
        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
    s(:,:,:) = 0.0

    do j=g%jsc,g%jec ; do i=g%isc,g%iec
      ! Compute normalized zonal coordinates (x,y=0 at center of domain)
      if (dbrotate) then
        ! This is really y in the rotated case
        x = ( g%geoLatT(i,j) ) / dblen
      else
        x = ( g%geoLonT(i,j) ) / dblen
      endif
      if (x>=0.25 ) then
        do k=1,nz
          s(i,j,k)=s_ref + 0.5*s_range
        enddo
      endif
      if (x<=-0.25 ) then
        do k=1,nz
          s(i,j,k)=s_ref - 0.5*s_range
        enddo
      endif
    enddo ; enddo
  endif

  if (associated(tv%S)) call set_up_ale_sponge_field(s, g, tv%S, acsp)

end subroutine dumbbell_initialize_sponges

end module dumbbell_initialization