lock_exchange_initialization.F90

!> Initialization of the "lock exchange" experiment.
!! lock_exchange = A 2-d density driven hydraulic exchange flow.
module lock_exchange_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_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

implicit none ; private

#include <MOM_memory.h>

public lock_exchange_initialize_thickness

contains

!> This subroutine initializes layer thicknesses for the lock_exchange experiment.
! -----------------------------------------------------------------------------
subroutine lock_exchange_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_(gv))     ! The resting interface heights [Z ~> m], usually
                           ! negative because it is positive upward.
  real :: e_pert(szk_(gv)) ! Interface height perturbations, positive upward [Z ~> m].
  real :: eta1d(szk_(gv)+1)! Interface height relative to the sea surface
                           ! positive upward [Z ~> m].
  real :: front_displacement ! Vertical displacement acrodd front
  real :: thermocline_thickness ! Thickness of stratified region
  logical :: just_read    ! If true, just read parameters but set nothing.
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "lock_exchange_initialize_thickness" ! This subroutine's name.
  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("  lock_exchange_initialization.F90, lock_exchange_initialize_thickness: setting thickness", 5)

  if (.not.just_read) call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "FRONT_DISPLACEMENT", front_displacement, &
                 "The vertical displacement of interfaces across the front. "//&
                 "A value larger in magnitude that MAX_DEPTH is truncated,", &
                 units="m", fail_if_missing=.not.just_read, do_not_log=just_read, scale=us%m_to_Z)
  call get_param(param_file, mdl, "THERMOCLINE_THICKNESS", thermocline_thickness, &
                 "The thickness of the thermocline in the lock exchange "//&
                 "experiment.  A value of zero creates a two layer system "//&
                 "with vanished layers in between the two inflated layers.", &
                 default=0., units="m", do_not_log=just_read, scale=us%m_to_Z)

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

  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    do k=2,nz
      eta1d(k) = -0.5 * g%max_depth & ! Middle of column
              - thermocline_thickness * ( (real(k-1))/real(nz) -0.5 ) ! Stratification
      if (g%geoLonT(i,j)-g%west_lon < 0.5 * g%len_lon) then
        eta1d(k) = eta1d(k) + 0.5 * front_displacement
      elseif (g%geoLonT(i,j)-g%west_lon > 0.5 * g%len_lon) then
        eta1d(k) = eta1d(k) - 0.5 * front_displacement
      endif
    enddo
    eta1d(nz+1) = -g%max_depth ! Force bottom interface to bottom
    do k=nz,2,-1 ! Make sure interfaces increase upwards
      eta1d(k) = max( eta1d(k), eta1d(k+1) + gv%Angstrom_Z )
    enddo
    eta1d(1) = 0. ! Force bottom interface to bottom
    do k=2,nz ! Make sure interfaces decrease downwards
      eta1d(k) = min( eta1d(k), eta1d(k-1) - gv%Angstrom_Z )
    enddo
    do k=nz,1,-1
      h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
    enddo
  enddo ; enddo

end subroutine lock_exchange_initialize_thickness
! -----------------------------------------------------------------------------

end module lock_exchange_initialization