benchmark_initialization.F90

!> Initialization for the "bench mark" configuration
module benchmark_initialization

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

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
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_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

implicit none ; private

#include <MOM_memory.h>

public benchmark_initialize_topography
public benchmark_initialize_thickness
public benchmark_init_temperature_salinity

! 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

!> This subroutine sets up the benchmark test case topography.
subroutine benchmark_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 ~> m] 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 :: pi                   ! 3.1415926... calculated as 4*atan(1)
  real :: d0                   ! A constant to make the maximum     !
                               ! basin depth MAXIMUM_DEPTH.         !
  real :: m_to_z  ! A dimensional rescaling factor.
  real :: x, y
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "benchmark_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("  benchmark_initialization.F90, benchmark_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)

  pi = 4.0*atan(1.0)
  d0 = max_depth / 0.5

!  Calculate the depth of the bottom.
  do j=js,je ; do i=is,ie
    x = (g%geoLonT(i,j)-g%west_lon) / g%len_lon
    y = (g%geoLatT(i,j)-g%south_lat) / g%len_lat
!  This sets topography that has a reentrant channel to the south.
    d(i,j) = -d0 * ( y*(1.0 + 0.6*cos(4.0*pi*x)) &
                   + 0.75*exp(-6.0*y) &
                   + 0.05*cos(10.0*pi*x) - 0.7 )
    if (d(i,j) > max_depth) d(i,j) = max_depth
    if (d(i,j) < min_depth) d(i,j) = 0.
  enddo ; enddo

end subroutine benchmark_initialize_topography

!> Initializes layer thicknesses for the benchmark test case,
!! by finding the depths of interfaces in a specified latitude-dependent
!! temperature profile with an exponentially decaying thermocline on top of a
!! linear stratification.
subroutine benchmark_initialize_thickness(h, G, GV, US, param_file, eqn_of_state, &
                                          P_Ref, 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(eos_type),          pointer     :: eqn_of_state !< Equation of state structure
  real,                    intent(in)  :: p_ref       !< The coordinate-density
                                                      !! reference pressure [R L2 T-2 ~> Pa].
  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_(gv)+1)     ! The resting interface heights, in depth units [Z ~> m],
                             ! usually negative because it is positive upward.
  real :: e_pert(szk_(gv)+1) ! Interface height perturbations, positive upward,
                             ! in depth units [Z ~> m].
  real :: eta1d(szk_(gv)+1)  ! Interface height relative to the sea surface
                             ! positive upward, in depth units [Z ~> m].
  real :: sst       !  The initial sea surface temperature [degC].
  real :: t_int     !  The initial temperature of an interface [degC].
  real :: ml_depth  !  The specified initial mixed layer depth, in depth units [Z ~> m].
  real :: thermocline_scale ! The e-folding scale of the thermocline, in depth units [Z ~> m].
  real, dimension(SZK_(GV)) :: &
    t0, s0, &       ! Profiles of temperature [degC] and salinity [ppt]
    rho_guess, &    ! Potential density at T0 & S0 [R ~> kg m-3].
    drho_dt, &      ! Derivative of density with temperature [R degC-1 ~> kg m-3 degC-1].
    drho_ds         ! Derivative of density with salinity [R ppt-1 ~> kg m-3 ppt-1].
  real :: pres(szk_(gv))  ! Reference pressure [R L2 T-2 ~> Pa].
  real :: a_exp     ! The fraction of the overall stratification that is exponential.
  real :: i_ts, i_md ! Inverse lengthscales [Z-1 ~> m-1].
  real :: t_frac    ! A ratio of the interface temperature to the range
                    ! between SST and the bottom temperature.
  real :: err, derr_dz  ! The error between the profile's temperature and the
                    ! interface temperature for a given z and its derivative.
  real :: pi, z
  logical :: just_read
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "benchmark_initialize_thickness" ! This subroutine's name.
  integer :: i, j, k, k1, is, ie, js, je, nz, itt

  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 log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "BENCHMARK_ML_DEPTH_IC", ml_depth, &
                 "Initial mixed layer depth in the benchmark test case.", &
                 units='m', default=50.0, scale=us%m_to_Z, do_not_log=just_read)
  call get_param(param_file, mdl, "BENCHMARK_THERMOCLINE_SCALE", thermocline_scale, &
                 "Initial thermocline depth scale in the benchmark test case.", &
                 default=500.0, units="m", scale=us%m_to_Z, do_not_log=just_read)

  if (just_read) return ! This subroutine has no run-time parameters.

  call mom_mesg("  benchmark_initialization.F90, benchmark_initialize_thickness: setting thickness", 5)

  k1 = gv%nk_rho_varies + 1

  a_exp = 0.9

! This block calculates T0(k) for the purpose of diagnosing where the
! interfaces will be found.
  do k=1,nz
    pres(k) = p_ref ; s0(k) = 35.0
  enddo
  t0(k1) = 29.0
  call calculate_density(t0(k1), s0(k1), pres(k1), rho_guess(k1), eqn_of_state)
  call calculate_density_derivs(t0(k1), s0(k1), pres(k1), drho_dt(k1), drho_ds(k1), eqn_of_state)

! A first guess of the layers' temperatures.
  do k=1,nz
    t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
  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
      t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
    enddo
  enddo

  pi = 4.0*atan(1.0)
  i_ts = 1.0 / thermocline_scale
  i_md = 1.0 / g%max_depth
  do j=js,je ; do i=is,ie
    sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
                               cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))

    do k=1,nz ; e_pert(k) = 0.0 ; enddo

    !   This sets the initial thickness (in [H ~> m or kg m-2]) of the layers.  The thicknesses
    ! are set to insure that: 1. each layer is at least  Gv%Angstrom_m thick, and
    ! 2. the interfaces are where they should be based on the resting depths and interface
    ! height perturbations, as long at this doesn't interfere with 1.
    eta1d(nz+1) = -g%bathyT(i,j)

    do k=nz,2,-1
      t_int = 0.5*(t0(k) + t0(k-1))
      t_frac = (t_int - t0(nz)) / (sst - t0(nz))
      ! Find the z such that T_frac = a exp(z/thermocline_scale) + (1-a) (z+D)/D
      z = 0.0
      do itt=1,6
        err = a_exp * exp(z*i_ts) + (1.0 - a_exp) * (z*i_md + 1.0) - t_frac
        derr_dz = a_exp * i_ts * exp(z*i_ts) + (1.0 - a_exp) * i_md
        z = z - err / derr_dz
      enddo
      e0(k) = z
!       e0(K) = -ML_depth + thermocline_scale * log((T_int - T0(nz)) / (SST - T0(nz)))

      eta1d(k) = e0(k) + e_pert(k)

      if (eta1d(k) > -ml_depth) eta1d(k) = -ml_depth

      if (eta1d(k) < eta1d(k+1) + gv%Angstrom_Z) &
        eta1d(k) = eta1d(k+1) + gv%Angstrom_Z

      h(i,j,k) = max(gv%Z_to_H * (eta1d(k) - eta1d(k+1)), gv%Angstrom_H)
    enddo
    h(i,j,1) = max(gv%Z_to_H * (0.0 - eta1d(2)), gv%Angstrom_H)

  enddo ; enddo

end subroutine benchmark_initialize_thickness

!> Initializes layer temperatures and salinities for benchmark
subroutine benchmark_init_temperature_salinity(T, S, G, GV, US, param_file, &
               eqn_of_state, P_Ref, 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.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: t       !< The potential temperature
                                                                   !! that is being initialized.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(out) :: s       !< The salinity that is being
                                                                   !! initialized.
  type(unit_scale_type),                     intent(in)  :: us     !< A dimensional unit scaling type
  type(param_file_type),               intent(in)  :: param_file   !< A structure indicating the
                                                                   !! open file to parse for
                                                                   !! model parameter values.
  type(eos_type),                      pointer     :: eqn_of_state !< Equation of state structure
  real,                                intent(in)  :: p_ref        !< The coordinate-density
                                                                   !! reference pressure [R L2 T-2 ~> Pa].
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                      !! only read parameters without changing h.
  ! Local variables
  real :: t0(szk_(g)), s0(szk_(g))
  real :: pres(szk_(g))      ! Reference pressure [R L2 T-2 ~> Pa].
  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 :: pi        ! 3.1415926... calculated as 4*atan(1)
  real :: sst       !  The initial sea surface temperature [degC].
  real :: lat
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=40)  :: mdl = "benchmark_init_temperature_salinity" ! This subroutine's name.
  integer :: i, j, k, k1, is, ie, js, je, nz, itt

  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 (just_read) return ! All run-time parameters have been read, so return.

  k1 = gv%nk_rho_varies + 1

  do k=1,nz
    pres(k) = p_ref ; s0(k) = 35.0
  enddo

  t0(k1) = 29.0
  call calculate_density(t0(k1), s0(k1), pres(k1), rho_guess(k1), eqn_of_state)
  call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state, (/k1,k1/) )

! A first guess of the layers' temperatures.                         !
  do k=1,nz
    t0(k) = t0(k1) + (gv%Rlay(k) - rho_guess(k1)) / drho_dt(k1)
  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
      t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
    enddo
  enddo

  do k=1,nz ; do i=is,ie ; do j=js,je
    t(i,j,k) = t0(k)
    s(i,j,k) = s0(k)
  enddo ; enddo ; enddo
  pi = 4.0*atan(1.0)
  do i=is,ie ; do j=js,je
    sst = 0.5*(t0(k1)+t0(nz)) - 0.9*0.5*(t0(k1)-t0(nz)) * &
                               cos(pi*(g%geoLatT(i,j)-g%south_lat)/(g%len_lat))
    do k=1,k1-1
      t(i,j,k) = sst
    enddo
  enddo ; enddo

end subroutine benchmark_init_temperature_salinity

end module benchmark_initialization