ISOMIP_initialization.F90

!> Configures the ISOMIP test case.
module isomip_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_sponge, only : sponge_cs, set_up_sponge_field, initialize_sponge
use mom_dyn_horgrid, only : dyn_horgrid_type
use mom_error_handler, only : mom_mesg, mom_error, fatal, is_root_pe, warning
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_io, only : file_exists
use mom_io, only : mom_read_data
use mom_io, only : slasher
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 regrid_consts, only : regridding_sigma_shelf_zstar
implicit none ; private

#include <MOM_memory.h>

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

! The following routines are visible to the outside world
public isomip_initialize_topography
public isomip_initialize_thickness
public isomip_initialize_temperature_salinity
public isomip_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 for the ISOMIP configuration
subroutine isomip_initialize_topography(D, G, param_file, max_depth, US)
  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 m or Z if US is present
  type(param_file_type),           intent(in)  :: param_file !< Parameter file structure
  real,                            intent(in)  :: max_depth !< Maximum model depth in the units of D
  type(unit_scale_type), optional, intent(in)  :: us !< A dimensional unit scaling type

  ! Local variables
  real :: min_depth ! The minimum and maximum depths [Z ~> m].
  real :: m_to_z  ! A dimensional rescaling factor.
  ! The following variables are used to set up the bathymetry in the ISOMIP example.
  real :: bmax            ! max depth of bedrock topography
  real :: b0,b2,b4,b6     ! first, second, third and fourth bedrock topography coeff
  real :: xbar            ! characteristic along-flow lenght scale of the bedrock
  real :: dc              ! depth of the trough compared with side walls [Z ~> m].
  real :: fc              ! characteristic width of the side walls of the channel
  real :: wc              ! half-width of the trough
  real :: ly              ! domain width (across ice flow)
  real :: bx, by          ! dummy vatiables [Z ~> m].
  real :: xtil            ! dummy vatiable
  logical :: is_2d         ! If true, use 2D setup
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "ISOMIP_initialize_topography" ! This subroutine's name.
  integer :: i, j, is, ie, js, je, isd, ied, jsd, jed
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  call mom_mesg("  ISOMIP_initialization.F90, ISOMIP_initialize_topography: setting topography", 5)

  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z

  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_z)
  call get_param(param_file, mdl, "ISOMIP_2D",is_2d,'If true, use a 2D setup.', default=.false.)

  ! The following variables should be transformed into runtime parameters?
  bmax = 720.0*m_to_z ; dc = 500.0*m_to_z
  b0 = -150.0*m_to_z ; b2 = -728.8*m_to_z ; b4 = 343.91*m_to_z ; b6 = -50.57*m_to_z
  xbar = 300.0e3 ; fc = 4.0e3 ; wc = 24.0e3 ; ly = 80.0e3
  bx = 0.0 ; by = 0.0 ; xtil = 0.0


  if (is_2d) then
    do j=js,je ; do i=is,ie
      ! 2D setup
      xtil = g%geoLonT(i,j)*1.0e3/xbar
      !xtil = 450*1.0e3/xbar
      bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
      !by = (dc/(1.+exp(-2.*(G%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + &
      !        (dc/(1.+exp(2.*(G%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc)))

      ! slice at y = 40 km
      by = (dc / (1.+exp(-2.*(40.0*1.0e3- ly/2. - wc)/fc))) + &
           (dc / (1.+exp(2.*(40.0*1.0e3- ly/2. + wc)/fc)))

      d(i,j) = -max(bx+by, -bmax)
      if (d(i,j) > max_depth) d(i,j) = max_depth
      if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
    enddo ; enddo

  else
    do j=js,je ; do i=is,ie
      ! 3D setup
      ! ===== TEST =====
      !if (G%geoLonT(i,j)<500.) then
      !  xtil = 500.*1.0e3/xbar
      !else
      !  xtil = G%geoLonT(i,j)*1.0e3/xbar
      !endif
      ! ===== TEST =====

      xtil = g%geoLonT(i,j)*1.0e3/xbar

      bx = b0 + b2*xtil**2 + b4*xtil**4 + b6*xtil**6
      by = (dc / (1.+exp(-2.*(g%geoLatT(i,j)*1.0e3- ly/2. - wc)/fc))) + &
           (dc / (1.+exp(2.*(g%geoLatT(i,j)*1.0e3- ly/2. + wc)/fc)))

      d(i,j) = -max(bx+by, -bmax)
      if (d(i,j) > max_depth) d(i,j) = max_depth
      if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
    enddo ; enddo
  endif

end subroutine isomip_initialize_topography

!> Initialization of thicknesses
subroutine isomip_initialize_thickness ( h, G, GV, US, param_file, tv, 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.
  type(thermo_var_ptrs),   intent(in)  :: tv          !< A structure containing pointers to any
                                                      !! available thermodynamic fields, including
                                                      !! the eqn. of state.
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
  ! Local variables
  real :: e0(szk_(g)+1)   ! The resting interface heights, in depth units [Z ~> m],
                          !  usually negative because it is positive upward.
  real :: eta1d(szk_(g)+1)! Interface height relative to the sea surface
                          ! positive upward, in depth units [Z ~> m].
  integer :: i, j, k, is, ie, js, je, nz, tmp1
  real    :: x
  real    :: min_thickness, s_sur, s_bot, t_sur, t_bot
  real    :: rho_sur, rho_bot  ! Surface and bottom densities [R ~> kg m-3]
  real    :: rho_range    ! The range of densities [R ~> kg m-3]
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=256) :: mesg  ! The text of an error message
  character(len=40) :: verticalcoordinate

  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")

  call get_param(param_file, mdl,"MIN_THICKNESS", min_thickness, &
                 'Minimum layer thickness', units='m', default=1.e-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)

  select case ( coordinatemode(verticalcoordinate) )

  case ( regridding_layer, regridding_rho ) ! Initial thicknesses for isopycnal coordinates
    call get_param(param_file, mdl, "ISOMIP_T_SUR", t_sur, &
                   "Temperature at the surface (interface)", units="degC", default=-1.9, do_not_log=just_read)
    call get_param(param_file, mdl, "ISOMIP_S_SUR", s_sur, &
                   "Salinity at the surface (interface)", units="ppt",  default=33.8, do_not_log=just_read)
    call get_param(param_file, mdl, "ISOMIP_T_BOT", t_bot, &
                   "Temperature at the bottom (interface)", units="degC", default=-1.9, do_not_log=just_read)
    call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot,&
                   "Salinity at the bottom (interface)", units="ppt", default=34.55, do_not_log=just_read)

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

    ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT
    call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state)
    ! write(mesg,*) 'Surface density is:', rho_sur
    ! call MOM_mesg(mesg,5)
    call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state)
    ! write(mesg,*) 'Bottom density is:', rho_bot
    ! call MOM_mesg(mesg,5)
    rho_range = rho_bot - rho_sur
    ! write(mesg,*) 'Density range is:', rho_range
    ! call MOM_mesg(mesg,5)

    ! Construct notional interface positions
    e0(1) = 0.
    do k=2,nz
      e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
      e0(k) = min( 0., e0(k) ) ! Bound by surface
      e0(k) = max( -g%max_depth, e0(k) ) ! Bound by possible deepest point in model
      ! write(mesg,*) 'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)', &
      !     G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)
      ! call MOM_mesg(mesg,5)
    enddo
    e0(nz+1) = -g%max_depth

    ! Calculate thicknesses
    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, regridding_sigma_shelf_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

  case default
      call mom_error(fatal,"isomip_initialize: "// &
      "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")

  end select

end subroutine isomip_initialize_thickness

!> Initial values for temperature and salinity
subroutine isomip_initialize_temperature_salinity ( T, S, h, G, GV, US, 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
  type(unit_scale_type),                     intent(in)  :: us !< A dimensional unit scaling type
  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 T & S.
  ! Local variables
  integer   :: i, j, k, is, ie, js, je, nz, itt
  real      :: x, ds, dt
  real      :: rho_sur, rho_bot  ! Surface and bottom densities [R ~> kg m-3]
  real      :: xi0, xi1 ! Heights in depth units [Z ~> m].
  real      :: s_sur, s_bot ! Salinity at the surface and bottom [ppt]
  real      :: t_sur, t_bot ! Temperature at the bottom [degC]
  real      :: dt_dz  ! Vertical gradient of temperature [degC Z-1 ~> degC m-1].
  real      :: ds_dz  ! Vertical gradient of salinity [ppt Z-1 ~> ppt m-1].
  real      :: z            ! vertical position in z space [Z ~> m]
  character(len=256) :: mesg ! The text of an error message
  character(len=40) :: verticalcoordinate, density_profile
  real :: rho_tmp
  logical :: just_read       ! If true, just read parameters but set nothing.
  logical :: fit_salin       ! If true, accept the prescribed temperature and fit the salinity.
  real :: t0(szk_(g)), s0(szk_(g))
  real :: drho_dt(szk_(g))   ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
  real :: drho_ds(szk_(g))   ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
  real :: rho_guess(szk_(g)) ! Potential density at T0 & S0 [R ~> kg m-3].
  real :: pres(szk_(g))      ! An array of the reference pressure [R L2 T-2 ~> Pa]. (zero here)
  real :: drho_dt1           ! A prescribed derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
  real :: drho_ds1           ! A prescribed derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
  real :: t_ref, s_ref
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  pres(:) = 0.0

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

  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, "ISOMIP_T_SUR",t_sur, &
                 "Temperature at the surface (interface)", units="degC", default=-1.9, do_not_log=just_read)
  call get_param(param_file, mdl, "ISOMIP_S_SUR", s_sur, &
                 "Salinity at the surface (interface)", units="ppt", default=33.8, do_not_log=just_read)
  call get_param(param_file, mdl, "ISOMIP_T_BOT", t_bot, &
                 "Temperature at the bottom (interface)", units="degC", default=-1.9, do_not_log=just_read)
  call get_param(param_file, mdl, "ISOMIP_S_BOT", s_bot, &
                 "Salinity at the bottom (interface)", units="ppt", default=34.55, do_not_log=just_read)

  call calculate_density(t_sur, s_sur, 0.0, rho_sur, eqn_of_state)
  ! write(mesg,*) 'Density in the surface layer:', rho_sur
  ! call MOM_mesg(mesg,5)
  call calculate_density(t_bot, s_bot, 0.0, rho_bot, eqn_of_state)
  ! write(mesg,*) 'Density in the bottom layer::', rho_bot
  ! call MOM_mesg(mesg,5)

  select case ( coordinatemode(verticalcoordinate) )

    case (  regridding_rho, regridding_zstar, regridding_sigma_shelf_zstar, regridding_sigma )
      if (just_read) return ! All run-time parameters have been read, so return.

      ds_dz = (s_sur - s_bot) / g%max_depth
      dt_dz = (t_sur - t_bot) / g%max_depth
      do j=js,je ; do i=is,ie
        xi0 = -g%bathyT(i,j)
        do k = nz,1,-1
          xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z ! Depth in middle of layer
          s(i,j,k) = s_sur + ds_dz * xi0
          t(i,j,k) = t_sur + dt_dz * xi0
          xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z ! Depth at top of layer
        enddo
      enddo ; enddo

    case ( regridding_layer )
      call get_param(param_file, mdl, "FIT_SALINITY", fit_salin, &
                  "If true, accept the prescribed temperature and fit the "//&
                  "salinity; otherwise take salinity and fit temperature.", &
                  default=.false., do_not_log=just_read)
      call get_param(param_file, mdl, "DRHO_DS", drho_ds1, &
                  "Partial derivative of density with salinity.", &
                  units="kg m-3 PSU-1", scale=us%kg_m3_to_R, fail_if_missing=.not.just_read, do_not_log=just_read)
      call get_param(param_file, mdl, "DRHO_DT", drho_dt1, &
                  "Partial derivative of density with temperature.", &
                  units="kg m-3 K-1", scale=us%kg_m3_to_R, fail_if_missing=.not.just_read, do_not_log=just_read)
      call get_param(param_file, mdl, "T_REF", t_ref, &
                  "A reference temperature used in initialization.", &
                  units="degC", fail_if_missing=.not.just_read, do_not_log=just_read)
      call get_param(param_file, mdl, "S_REF", s_ref, &
                  "A reference salinity used in initialization.", units="PSU", &
                  default=35.0, do_not_log=just_read)
      if (just_read) return ! All run-time parameters have been read, so return.

      ! write(mesg,*) 'read drho_dS, drho_dT', drho_dS1, drho_dT1
      ! call MOM_mesg(mesg,5)

      ds_dz = (s_sur - s_bot) / g%max_depth
      dt_dz = (t_sur - t_bot) / g%max_depth

      do j=js,je ; do i=is,ie
        xi0 = 0.0
        do k = 1,nz
          !T0(k) = T_Ref; S0(k) = S_Ref
          xi1 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z
          s0(k) = s_sur - ds_dz * xi1
          t0(k) = t_sur - dt_dz * xi1
          xi0 = xi0 + h(i,j,k) * gv%H_to_Z
          ! write(mesg,*) 'S,T,xi0,xi1,k',S0(k),T0(k),xi0,xi1,k
          ! call MOM_mesg(mesg,5)
        enddo

        call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state, (/1,1/) )
        ! write(mesg,*) 'computed drho_dS, drho_dT', drho_dS(1), drho_dT(1)
        ! call MOM_mesg(mesg,5)
        call calculate_density(t0(1), s0(1), pres(1), rho_guess(1), eqn_of_state)

        if (fit_salin) then
          ! A first guess of the layers' salinity.
          do k=nz,1,-1
            s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds1)
          enddo
          ! Refine the guesses for each layer.
          do itt=1,6
            call calculate_density(t0, s0, pres, rho_guess, eqn_of_state)
            call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state)
            do k=1,nz
              s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds1)
            enddo
          enddo

        else
          ! A first guess of the layers' temperatures.
          do k=nz,1,-1
            t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt1
          enddo

          do itt=1,6
            call calculate_density(t0, s0, pres, rho_guess, eqn_of_state)
            call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state)
            do k=1,nz
              t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
            enddo
          enddo
        endif

        do k=1,nz
          t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
        enddo

      enddo ; enddo

    case default
      call mom_error(fatal,"isomip_initialize: "// &
      "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")

  end select

  ! for debugging
  !i=G%iec; j=G%jec
  !do k = 1,nz
  !  call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,eqn_of_state, scale=US%kg_m3_to_R)
  !  write(mesg,*) 'k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k)
  !  call MOM_mesg(mesg,5)
  !enddo

end subroutine isomip_initialize_temperature_salinity

!> Sets up the the inverse restoration time (Idamp), and
! the values towards which the interface heights and an arbitrary
! number of tracers should be restored within each sponge.
subroutine isomip_initialize_sponges(G, GV, US, tv, PF, use_ALE, CSp, ACSp)
  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 containing pointers
                                            !! to any available thermodynamic
                                            !! fields, potential temperature and
                                            !! salinity or mixed layer density.
                                            !! Absent fields have NULL ptrs.
  type(param_file_type), intent(in) :: pf   !< A structure indicating the
                                            !! open file to parse for model
                                            !! parameter values.
  logical, intent(in) :: use_ale            !< If true, indicates model is in ALE mode
  type(sponge_cs),   pointer    :: csp      !< Layer-mode sponge structure
  type(ale_sponge_cs),   pointer    :: acsp !< ALE-mode sponge structure
  ! Local variables
  real :: t(szi_(g),szj_(g),szk_(g))  ! A temporary array for temp [degC]
  real :: s(szi_(g),szj_(g),szk_(g))  ! A temporary array for salt [ppt]
  ! real :: RHO(SZI_(G),SZJ_(G),SZK_(G))  ! A temporary array for RHO [R ~> kg m-3]
  real :: h(szi_(g),szj_(g),szk_(g))  ! A temporary array for thickness [H ~> m or kg m-2]
  real :: idamp(szi_(g),szj_(g))    ! The inverse damping rate [T-1 ~> s-1].
  real :: tnudg                     ! Nudging time scale [T ~> s]
  real :: s_sur, t_sur              ! Surface salinity and temerature in sponge
  real :: s_bot, t_bot              ! Bottom salinity and temerature in sponge
  real :: t_ref, s_ref              ! reference T and S
  real :: rho_sur, rho_bot          ! Surface and bottom densities [R ~> kg m-3]
  real :: rho_range                 ! The range of densities [R ~> kg m-3]
  real :: dt_dz, ds_dz              ! Gradients of T and S in degC/Z and PPT/Z.

  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 :: eta(szi_(g),szj_(g),szk_(g)+1) ! A temporary array for eta [Z ~> m].
  real :: min_depth, dummy1, z
  real :: rho_dummy, min_thickness, rho_tmp, xi0
  character(len=40) :: verticalcoordinate, filename, state_file
  character(len=40) :: temp_var, salt_var, eta_var, inputdir

  character(len=40)  :: mdl = "ISOMIP_initialize_sponges" ! This subroutine's name.
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  call get_param(pf, mdl, "MIN_THICKNESS", min_thickness, "Minimum layer thickness", &
                 units="m", default=1.e-3, scale=us%m_to_Z)

  call get_param(pf, mdl, "REGRIDDING_COORDINATE_MODE", verticalcoordinate, &
            default=default_coordinate_mode)

  call get_param(pf, mdl, "ISOMIP_TNUDG", tnudg, "Nudging time scale for sponge layers (days)", &
                 default=0.0, scale=86400.0*us%s_to_T)

  call get_param(pf, mdl, "T_REF", t_ref, "Reference temperature", default=10.0, &
                 do_not_log=.true.)

  call get_param(pf, mdl, "S_REF", s_ref, "Reference salinity", default=35.0, &
                 do_not_log=.true.)

  call get_param(pf, mdl, "ISOMIP_S_SUR_SPONGE", s_sur, &
                 "Surface salinity in sponge layer.", units="ppt", default=s_ref) ! units="ppt")

  call get_param(pf, mdl, "ISOMIP_S_BOT_SPONGE", s_bot, &
                 "Bottom salinity in sponge layer.", units="ppt", default=s_ref) ! units="ppt")

  call get_param(pf, mdl, "ISOMIP_T_SUR_SPONGE", t_sur, &
                 "Surface temperature in sponge layer.", units="degC", default=t_ref) ! units="degC")

  call get_param(pf, mdl, "ISOMIP_T_BOT_SPONGE", t_bot, &
                 "Bottom temperature in sponge layer.", units="degC", default=t_ref) ! units="degC")

  t(:,:,:) = 0.0 ; s(:,:,:) = 0.0 ; idamp(:,:) = 0.0 !; RHO(:,:,:) = 0.0

!   Set up sponges for ISOMIP configuration
  call get_param(pf, mdl, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0, scale=us%m_to_Z)

  if (associated(csp)) call mom_error(fatal, &
        "ISOMIP_initialize_sponges called with an associated control structure.")
  if (associated(acsp)) call mom_error(fatal, &
        "ISOMIP_initialize_sponges called with an associated ALE-sponge control structure.")

  !  Here the inverse damping time [T-1 ~> s-1], is set. Set Idamp to 0
  !  wherever there is no sponge, and the subroutines that are called
  !  will automatically set up the sponges only where Idamp is positive
  !  and mask2dT is 1.

  do i=is,ie; do j=js,je
    if (g%bathyT(i,j) <= min_depth) then
      idamp(i,j) = 0.0
    elseif (g%geoLonT(i,j) >= 790.0 .AND. g%geoLonT(i,j) <= 800.0) then
      dummy1 = (g%geoLonT(i,j)-790.0)/(800.0-790.0)
      idamp(i,j) = (1.0/tnudg) * max(0.0,dummy1)
    else
      idamp(i,j) = 0.0
    endif

  enddo ; enddo

  ! Compute min/max density using T_SUR/S_SUR and T_BOT/S_BOT
  call calculate_density(t_sur, s_sur, 0.0, rho_sur, tv%eqn_of_state)
  !write (mesg,*) 'Surface density in sponge:', rho_sur
  ! call MOM_mesg(mesg,5)
  call calculate_density(t_bot, s_bot, 0.0, rho_bot, tv%eqn_of_state)
  !write (mesg,*) 'Bottom density in sponge:', rho_bot
  ! call MOM_mesg(mesg,5)
  rho_range = rho_bot - rho_sur
  !write (mesg,*) 'Density range in sponge:', rho_range
  ! call MOM_mesg(mesg,5)

  if (use_ale) then

    select case ( coordinatemode(verticalcoordinate) )

     case ( regridding_rho )
       ! Construct notional interface positions
       e0(1) = 0.
       do k=2,nz
         e0(k) = -g%max_depth * ( 0.5 * ( gv%Rlay(k-1) + gv%Rlay(k) ) - rho_sur ) / rho_range
         e0(k) = min( 0., e0(k) ) ! Bound by surface
         e0(k) = max( -g%max_depth, e0(k) ) ! Bound by possible deepest point in model
         ! write(mesg,*) 'G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)',&
         !       G%max_depth,GV%Rlay(k-1),GV%Rlay(k),e0(k)
         ! call MOM_mesg(mesg,5)
       enddo
       e0(nz+1) = -g%max_depth

       ! Calculate thicknesses
       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, regridding_sigma_shelf_zstar )   ! Initial thicknesses for z coordinates
       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) = min_thickness * gv%Z_to_H
           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
        do j=js,je ; do i=is,ie
          h(i,j,:) = gv%Z_to_H * (g%bathyT(i,j) / dfloat(nz))
        enddo ; enddo

      case default
         call mom_error(fatal,"ISOMIP_initialize_sponges: "// &
         "Unrecognized i.c. setup - set REGRIDDING_COORDINATE_MODE")

    end select

    !  This call sets up the damping rates and interface heights.
    !  This sets the inverse damping timescale fields in the sponges.
    call initialize_ale_sponge(idamp, g, pf, acsp, h, nz)

    ds_dz = (s_sur - s_bot) / g%max_depth
    dt_dz = (t_sur - t_bot) / g%max_depth
    do j=js,je ; do i=is,ie
      xi0 = -g%bathyT(i,j)
      do k = nz,1,-1
        xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z ! Depth in middle of layer
        s(i,j,k) = s_sur + ds_dz * xi0
        t(i,j,k) = t_sur + dt_dz * xi0
        xi0 = xi0 + 0.5 * h(i,j,k) * gv%H_to_Z ! Depth at top of layer
      enddo
    enddo ; enddo
    ! for debugging
    !i=G%iec; j=G%jec
    !do k = 1,nz
    !  call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state, scale=US%kg_m3_to_R)
    !  write(mesg,*) 'Sponge - k,h,T,S,rho,Rlay',k,h(i,j,k),T(i,j,k),S(i,j,k),rho_tmp,GV%Rlay(k)
    !  call MOM_mesg(mesg,5)
    !enddo

    !   Now register all of the fields which are damped in the sponge.   !
    ! By default, momentum is advected vertically within the sponge, but !
    ! momentum is typically not damped within the sponge.                !

    !  The remaining calls to set_up_sponge_field can be in any order. !
    if ( associated(tv%T) ) then
      call set_up_ale_sponge_field(t, g, tv%T, acsp)
    endif
    if ( associated(tv%S) ) then
      call set_up_ale_sponge_field(s, g, tv%S, acsp)
    endif

  else ! layer mode
    ! 1) Read eta, salt and temp from IC file
    call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
    inputdir = slasher(inputdir)
    ! GM: get two different files, one with temp and one with salt values
    ! this is work around to avoid having wrong values near the surface
    ! because of the FIT_SALINITY option. To get salt values right in the
    ! sponge, FIT_SALINITY=False. The oposite is true for temp. One can
    ! combined the *correct* temp and salt values in one file instead.
    call get_param(pf, mdl, "ISOMIP_SPONGE_FILE", state_file, &
              "The name of the file with temps., salts. and interfaces to "//&
              "damp toward.", fail_if_missing=.true.)
    call get_param(pf, mdl, "SPONGE_PTEMP_VAR", temp_var, &
              "The name of the potential temperature variable in "//&
              "SPONGE_STATE_FILE.", default="Temp")
    call get_param(pf, mdl, "SPONGE_SALT_VAR", salt_var, &
              "The name of the salinity variable in "//&
              "SPONGE_STATE_FILE.", default="Salt")
    call get_param(pf, mdl, "SPONGE_ETA_VAR", eta_var, &
              "The name of the interface height variable in "//&
              "SPONGE_STATE_FILE.", default="eta")

    !read temp and eta
    filename = trim(inputdir)//trim(state_file)
    if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
          "ISOMIP_initialize_sponges: Unable to open "//trim(filename))
    call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
    call mom_read_data(filename, temp_var, t(:,:,:), g%Domain)
    call mom_read_data(filename, salt_var, s(:,:,:), g%Domain)

    ! for debugging
    !i=G%iec; j=G%jec
    !do k = 1,nz
    !  call calculate_density(T(i,j,k),S(i,j,k),0.0,rho_tmp,tv%eqn_of_state, scale=US%kg_m3_to_R)
    !  write(mesg,*) 'Sponge - k,eta,T,S,rho,Rlay',k,eta(i,j,k),T(i,j,k),&
    !              S(i,j,k),rho_tmp,GV%Rlay(k)
    !  call MOM_mesg(mesg,5)
    !enddo

    ! Set the inverse damping rates so that the model will know where to
    ! apply the sponges, along with the interface heights.
    call initialize_sponge(idamp, eta, g, pf, csp, gv)
    ! Apply sponge in tracer fields
    call set_up_sponge_field(t, tv%T, g, nz, csp)
    call set_up_sponge_field(s, tv%S, g, nz, csp)

  endif

end subroutine isomip_initialize_sponges

!> \namespace isomip_initialization
!!
!! See this paper for details: http://www.geosci-model-dev-discuss.net/8/9859/2015/gmdd-8-9859-2015.pdf
end module isomip_initialization