Rossby_front_2d_initialization.F90

!> Initial conditions for the 2D Rossby front test
module rossby_front_2d_initialization

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

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

implicit none ; private

#include <MOM_memory.h>

! Private (module-wise) parameters
character(len=40) :: mdl = "Rossby_front_2d_initialization" !< This module's name.
! This include declares and sets the variable "version".
#include "version_variable.h"

public rossby_front_initialize_thickness
public rossby_front_initialize_temperature_salinity
public rossby_front_initialize_velocity

! Parameters defining the initial conditions of this test case
real, parameter :: frontfractionalwidth = 0.5 !< Width of front as fraction of domain [nondim]
real, parameter :: hmlmin = 0.25 !< Shallowest ML as fractional depth of ocean [nondim]
real, parameter :: hmlmax = 0.75 !< Deepest ML as fractional depth of ocean [nondim]

contains

!> Initialization of thicknesses in 2D Rossby front test
subroutine rossby_front_initialize_thickness(h, G, GV, US, param_file, just_read_params)
  type(ocean_grid_type),   intent(in)  :: g           !< 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_(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.

  integer :: i, j, k, is, ie, js, je, nz
  real    :: tz, dml, eta, stretch, h0
  real    :: min_thickness, t_range
  real    :: drho_dt      ! The partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
  logical :: just_read    ! If true, just read parameters but set nothing.
  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("Rossby_front_2d_initialization.F90, Rossby_front_initialize_thickness: setting thickness")

  if (.not.just_read) call log_version(param_file, mdl, version, "")
  ! Read parameters needed to set 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)
  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, "T_RANGE", t_range, 'Initial temperature range', &
                 units='C', default=0.0, do_not_log=just_read)
  call get_param(param_file, mdl, "DRHO_DT", drho_dt, default=-0.2, scale=us%kg_m3_to_R, do_not_log=.true.)

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

  tz = t_range / g%max_depth

  select case ( coordinatemode(verticalcoordinate) )

    case (regridding_layer, regridding_rho)
      do j = g%jsc,g%jec ; do i = g%isc,g%iec
        dml = hml( g, g%geoLatT(i,j) )
        eta = -( -drho_dt / gv%Rho0 ) * tz * 0.5 * ( dml * dml )
        stretch = ( ( g%max_depth + eta ) / g%max_depth )
        h0 = ( g%max_depth / real(nz) ) * stretch
        do k = 1, nz
          h(i,j,k) = h0 * gv%Z_to_H
        enddo
      enddo ; enddo

    case (regridding_zstar, regridding_sigma)
      do j = g%jsc,g%jec ; do i = g%isc,g%iec
        dml = hml( g, g%geoLatT(i,j) )
        eta = -( -drho_dt / gv%Rho0 ) * tz * 0.5 * ( dml * dml )
        stretch = ( ( g%max_depth + eta ) / g%max_depth )
        h0 = ( g%max_depth / real(nz) ) * stretch
        do k = 1, nz
          h(i,j,k) = h0 * gv%Z_to_H
        enddo
      enddo ; enddo

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

  end select

end subroutine rossby_front_initialize_thickness


!> Initialization of temperature and salinity in the Rossby front test
subroutine rossby_front_initialize_temperature_salinity(T, S, h, G, GV, &
                   param_file, eqn_of_state, just_read_params)
  type(ocean_grid_type),                     intent(in)  :: g  !< 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  !< 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  !< Thickness [H ~> m or kg m-2]
  type(param_file_type),                     intent(in)  :: param_file   !< Parameter file handle
  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.

  integer   :: i, j, k, is, ie, js, je, nz
  real      :: t_ref, s_ref ! Reference salinity and temerature within surface layer
  real      :: t_range      ! Range of salinities and temperatures over the vertical
  real      :: y, zc, zi, dtdz
  logical :: just_read    ! If true, just read parameters but set nothing.
  character(len=40) :: verticalcoordinate
  real      :: pi                   ! 3.1415926... calculated as 4*atan(1)

  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

  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, "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='C',&
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"T_RANGE",t_range,'Initial temperature range',&
                 units='C', default=0.0, do_not_log=just_read)

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

  t(:,:,:) = 0.0
  s(:,:,:) = s_ref
  dtdz = t_range / g%max_depth

  do j = g%jsc,g%jec ; do i = g%isc,g%iec
    zi = 0.
    do k = 1, nz
      zi = zi - h(i,j,k)              ! Bottom interface position
      zc = gv%H_to_Z * (zi - 0.5*h(i,j,k))    ! Position of middle of cell
      zc = min( zc, -hml(g, g%geoLatT(i,j)) ) ! Bound by depth of mixed layer
      t(i,j,k) = t_ref + dtdz * zc ! Linear temperature profile
    enddo
  enddo ; enddo

end subroutine rossby_front_initialize_temperature_salinity


!> Initialization of u and v in the Rossby front test
subroutine rossby_front_initialize_velocity(u, v, h, G, GV, US, param_file, just_read_params)
  type(ocean_grid_type),      intent(in)  :: g  !< Grid structure
  type(verticalgrid_type),    intent(in)  :: gv !< Vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                              intent(out) :: u  !< i-component of velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                              intent(out) :: v  !< j-component of velocity [L T-1 ~> m s-1]
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), &
                              intent(in)  :: h  !< Thickness [H ~> m or kg m-2]
  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.
  logical,          optional, intent(in)  :: just_read_params !< If present and true, this call
                                                !! will only read parameters without setting u & v.

  real    :: y            ! Non-dimensional coordinate across channel, 0..pi
  real    :: t_range      ! Range of salinities and temperatures over the vertical
  real    :: dudt         ! Factor to convert dT/dy into dU/dz, g*alpha/f [L2 Z-1 T-1 degC-1 ~> m s-1 degC-1]
  real    :: drho_dt      ! The partial derivative of density with temperature [R degC-1 ~> kg m-3 degC-1]
  real    :: dml, zi, zc, zm ! Depths [Z ~> m].
  real    :: f            ! The local Coriolis parameter [T-1 ~> s-1]
  real    :: ty           ! The meridional temperature gradient [degC L-1 ~> degC m-1]
  real    :: hatu         ! Interpolated layer thickness [Z ~> m].
  integer :: i, j, k, is, ie, js, je, nz
  logical :: just_read    ! If true, just read parameters but set nothing.
  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

  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, "T_RANGE", t_range, 'Initial temperature range', &
                 units='C', default=0.0, do_not_log=just_read)
  call get_param(param_file, mdl, "DRHO_DT", drho_dt, default=-0.2, scale=us%kg_m3_to_R, do_not_log=.true.)

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

  v(:,:,:) = 0.0
  u(:,:,:) = 0.0

  do j = g%jsc,g%jec ; do i = g%isc-1,g%iec+1
    f = 0.5* (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1) )
    dudt = 0.0 ; if (abs(f) > 0.0) &
      dudt = ( gv%g_Earth*drho_dt ) / ( f * gv%Rho0 )
    dml = hml( g, g%geoLatT(i,j) )
    ty = us%L_to_m*dtdy( g, t_range, g%geoLatT(i,j) )
    zi = 0.
    do k = 1, nz
      hatu = 0.5*(h(i,j,k)+h(i+1,j,k)) * gv%H_to_Z
      zi = zi - hatu              ! Bottom interface position
      zc = zi - 0.5*hatu          ! Position of middle of cell
      zm = max( zc + dml, 0. )    ! Height above bottom of mixed layer
      u(i,j,k) = dudt * ty * zm   ! Thermal wind starting at base of ML
    enddo
  enddo ; enddo

end subroutine rossby_front_initialize_velocity

!> Pseudo coordinate across domain used by Hml() and dTdy()
!! returns a coordinate from -PI/2 .. PI/2 squashed towards the
!! center of the domain.
real function ypseudo( G, lat )
  type(ocean_grid_type), intent(in) :: g   !< Grid structure
  real,                  intent(in) :: lat !< Latitude
  ! Local
  real :: y, pi

  pi = 4.0 * atan(1.0)
  ypseudo = ( ( lat - g%south_lat ) / g%len_lat ) - 0.5 ! -1/2 .. 1/.2
  ypseudo = pi * max(-0.5, min(0.5, ypseudo / frontfractionalwidth))
end function ypseudo


!> Analytic prescription of mixed layer depth in 2d Rossby front test,
!! in the same units as G%max_depth
real function hml( G, lat )
  type(ocean_grid_type), intent(in) :: g   !< Grid structure
  real,                  intent(in) :: lat !< Latitude
  ! Local
  real :: dhml, hmlmean

  dhml = 0.5 * ( hmlmax - hmlmin ) * g%max_depth
  hmlmean = 0.5 * ( hmlmin + hmlmax ) * g%max_depth
  hml = hmlmean + dhml * sin( ypseudo(g, lat) )
end function hml


!> Analytic prescription of mixed layer temperature gradient in 2d Rossby front test
real function dtdy( G, dT, lat )
  type(ocean_grid_type), intent(in) :: g     !< Grid structure
  real,                  intent(in) :: dt    !< Top to bottom temperature difference
  real,                  intent(in) :: lat   !< Latitude
  ! Local
  real :: pi, dhml, dhdy
  real :: km = 1.e3 ! AXIS_UNITS = 'k' (1000 m)

  pi = 4.0 * atan(1.0)
  dhml = 0.5 * ( hmlmax - hmlmin ) * g%max_depth
  dhdy = dhml * ( pi / ( frontfractionalwidth * g%len_lat * km ) ) * cos( ypseudo(g, lat) )
  dtdy = -( dt / g%max_depth ) * dhdy

end function dtdy


!> \namespace rossby_front_2d_initialization
!!
!! \section section_Rossby_front_2d Description of the 2d Rossby front initial conditions
!!
!! Consistent with a linear equation of state, the system has a constant stratification
!! below the mixed layer, stratified in temperature only. Isotherms are flat below the
!! mixed layer and vertical within. Salinity is constant. The mixed layer has a half sine
!! form so that there are no mixed layer or temperature gradients at the side walls.
!!
!! Below the mixed layer the potential temperature, \f$\theta(z)\f$, is given by
!! \f[ \theta(z) = \theta_0 - \Delta \theta \left( z + h_{ML} \right) \f]
!! where \f$ \theta_0 \f$ and \f$ \Delta \theta \f$ are external model parameters.
!!
!! The depth of the mixed layer, \f$H_{ML}\f$ is
!! \f[ h_{ML}(y) = h_{min} + \left( h_{max} - h_{min} \right) \cos{\pi y/L} \f].
!! The temperature in mixed layer is given by the reference temperature at \f$z=h_{ML}\f$
!! so that
!! \f{eqnarray} \theta(y,z) =
!!     \theta_0 - \Delta \theta \left( z + h_{ML} \right) & \forall & z < h_{ML}(y) T.B.D.
!! \f}

end module rossby_front_2d_initialization