baroclinic_zone_initialization.F90

!> Initial conditions for an idealized baroclinic zone
module baroclinic_zone_initialization
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_file_parser,   only : openparameterblock, closeparameterblock
use mom_grid, only : ocean_grid_type
use mom_unit_scaling, only : unit_scale_type
use mom_verticalgrid, only : verticalgrid_type
 
implicit none ; private
 
#include <MOM_memory.h>
#include "version_variable.h"
 
! Private (module-wise) parameters
character(len=40) :: mdl = "baroclinic_zone_initialization" !< This module's name.
 
public baroclinic_zone_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
 
!> Reads the parameters unique to this module
subroutine bcz_params(G, GV, US, param_file, S_ref, dSdz, delta_S, dSdx, T_ref, dTdz, &
                      delta_T, dTdx, L_zone, just_read_params)
  type(ocean_grid_type),   intent(in)  :: G          !< 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(param_file_type),   intent(in)  :: param_file !< Parameter file handle
  real,                    intent(out) :: S_ref      !< Reference salinity [ppt]
  real,                    intent(out) :: dSdz       !< Salinity stratification [ppt Z-1 ~> ppt m-1]
  real,                    intent(out) :: delta_S    !< Salinity difference across baroclinic zone [ppt]
  real,                    intent(out) :: dSdx       !< Linear salinity gradient [ppt m-1]
  real,                    intent(out) :: T_ref      !< Reference temperature [degC]
  real,                    intent(out) :: dTdz       !< Temperature stratification [degC Z-1 ~> degC m-1]
  real,                    intent(out) :: delta_T    !< Temperature difference across baroclinic zone [degC]
  real,                    intent(out) :: dTdx       !< Linear temperature gradient [degC m-1]
  real,                    intent(out) :: L_zone     !< Width of baroclinic zone [m]
  logical,       optional, intent(in)  :: just_read_params !< If present and true, this call will
                                                     !! only read parameters without changing h.
 
  logical :: just_read    ! If true, just read parameters but set nothing.
 
  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params
 
  if (.not.just_read) &
    call log_version(param_file, mdl, version, 'Initialization of an analytic baroclinic zone')
  call openparameterblock(param_file,'BCZIC')
  call get_param(param_file, mdl, "S_REF", s_ref, 'Reference salinity', units='ppt', &
                 default=35., do_not_log=just_read)
  call get_param(param_file, mdl, "DSDZ", dsdz, 'Salinity stratification', &
                 units='ppt/m', default=0.0, scale=us%Z_to_m, do_not_log=just_read)
  call get_param(param_file, mdl,"DELTA_S",delta_s,'Salinity difference across baroclinic zone', &
                 units='ppt', default=0.0, do_not_log=just_read)
  call get_param(param_file, mdl,"DSDX",dsdx,'Meridional salinity difference', &
                 units='ppt/'//trim(g%x_axis_units), default=0.0, do_not_log=just_read)
  call get_param(param_file, mdl,"T_REF",t_ref,'Reference temperature',units='C', &
                 default=10., do_not_log=just_read)
  call get_param(param_file, mdl, "DTDZ", dtdz, 'Temperature stratification', &
                 units='C/m', default=0.0, scale=us%Z_to_m, do_not_log=just_read)
  call get_param(param_file, mdl,"DELTA_T",delta_t,'Temperature difference across baroclinic zone', &
                 units='C', default=0.0, do_not_log=just_read)
  call get_param(param_file, mdl,"DTDX",dtdx,'Meridional temperature difference', &
                 units='C/'//trim(g%x_axis_units), default=0.0, do_not_log=just_read)
  call get_param(param_file, mdl,"L_ZONE",l_zone,'Width of baroclinic zone', &
                 units=g%x_axis_units, default=0.5*g%len_lat, do_not_log=just_read)
  call closeparameterblock(param_file)
 
end subroutine bcz_params
 
!> Initialization of temperature and salinity with the baroclinic zone initial conditions
subroutine baroclinic_zone_init_temperature_salinity(T, S, h, G, GV, US, param_file, &
                                                     just_read_params)
  type(ocean_grid_type),                     intent(in)  :: g  !< 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_(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  !< The model thicknesses [H ~> m or kg m-2]
  type(param_file_type),                     intent(in)  :: param_file  !< Parameter file handle
  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, dtdz, dtdx, delta_t ! Parameters describing temperature distribution
  real      :: s_ref, dsdz, dsdx, delta_s ! Parameters describing salinity distribution
  real      :: l_zone ! Width of baroclinic zone
  real      :: zc, zi ! Depths in depth units [Z ~> m]
  real      :: x, xd, xs, y, yd, fn
  real      :: pi                   ! 3.1415926... calculated as 4*atan(1)
  logical :: just_read    ! If true, just read parameters but set nothing.
 
  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 bcz_params(g, gv, us, param_file, s_ref, dsdz, delta_s, dsdx, t_ref, dtdz, &
                  delta_t, dtdx, l_zone, just_read_params)
 
  if (just_read) return ! All run-time parameters have been read, so return.
 
  t(:,:,:) = 0.
  s(:,:,:) = 0.
  pi = 4.*atan(1.)
 
  do j = g%jsc,g%jec ; do i = g%isc,g%iec
    zi = -g%bathyT(i,j)
    x = g%geoLonT(i,j) - (g%west_lon + 0.5*g%len_lon) ! Relative to center of domain
    xd = x / g%len_lon ! -1/2 < xd 1/2
    y = g%geoLatT(i,j) - (g%south_lat + 0.5*g%len_lat) ! Relative to center of domain
    yd = y / g%len_lat ! -1/2 < yd 1/2
    if (l_zone/=0.) then
      xs = min(1., max(-1., x/l_zone)) ! -1 < ys < 1
      fn = sin((0.5*pi)*xs)
    else
      xs = sign(1., x) ! +/- 1
      fn = xs
    endif
    do k = nz, 1, -1
      zc = zi + 0.5*h(i,j,k)*gv%H_to_Z ! Position of middle of cell
      zi = zi + h(i,j,k)*gv%H_to_Z    ! Top interface position
      t(i,j,k) = t_ref + dtdz * zc  & ! Linear temperature stratification
                 + dtdx * x         & ! Linear gradient
                 + delta_t * fn       ! Smooth fn of width L_zone
      s(i,j,k) = s_ref + dsdz * zc  & ! Linear temperature stratification
                 + dsdx * x         & ! Linear gradient
                 + delta_s * fn       ! Smooth fn of width L_zone
    enddo
  enddo ; enddo
 
end subroutine baroclinic_zone_init_temperature_salinity
 
!> \namespace baroclinic_zone_initialization
!!
!! \section section_baroclinic_zone Description of the baroclinic zone initial conditions
!!
!! yada yada yada
 
end module baroclinic_zone_initialization