shelfwave_initialization.F90

!> Configures the model for the idealized shelfwave test case.
module shelfwave_initialization

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

use mom_domains,        only : sum_across_pes
use mom_dyn_horgrid,    only : dyn_horgrid_type
use mom_error_handler,  only : mom_mesg, mom_error, fatal, warning, is_root_pe
use mom_file_parser,    only : get_param, log_version, param_file_type
use mom_grid,           only : ocean_grid_type
use mom_open_boundary,  only : ocean_obc_type, obc_none, obc_direction_w
use mom_open_boundary,  only : obc_segment_type, register_obc
use mom_open_boundary,  only : obc_registry_type
use mom_time_manager,   only : time_type, time_type_to_real
use mom_unit_scaling,   only : unit_scale_type

implicit none ; private

#include <MOM_memory.h>

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

! The following routines are visible to the outside world
public shelfwave_initialize_topography
public shelfwave_set_obc_data
public register_shelfwave_obc, shelfwave_obc_end

!> Control structure for shelfwave open boundaries.
type, public :: shelfwave_obc_cs ; private
  real :: lx = 100.0        !< Long-shore length scale of bathymetry.
  real :: ly = 50.0         !< Cross-shore length scale.
  real :: f0 = 1.e-4        !< Coriolis parameter.
  real :: jj = 1            !< Cross-shore wave mode.
  real :: kk                !< Parameter.
  real :: ll                !< Longshore wavenumber.
  real :: alpha             !< 1/Ly.
  real :: omega             !< Frequency.
end type shelfwave_obc_cs

contains

!> Add shelfwave to OBC registry.
function register_shelfwave_obc(param_file, CS, OBC_Reg)
  type(param_file_type),    intent(in) :: param_file !< parameter file.
  type(shelfwave_obc_cs),   pointer    :: cs         !< shelfwave control structure.
  type(obc_registry_type),  pointer    :: obc_reg    !< OBC registry.
  logical                              :: register_shelfwave_obc
  ! Local variables
  real :: kk, ll, pi, len_lat

  character(len=32)  :: casename = "shelfwave"       !< This case's name.

  pi = 4.0*atan(1.0)

  if (associated(cs)) then
    call mom_error(warning, "register_shelfwave_OBC called with an "// &
                            "associated control structure.")
    return
  endif
  allocate(cs)

  ! Register the tracer for horizontal advection & diffusion.
  call register_obc(casename, param_file, obc_reg)
  call get_param(param_file, mdl,"F_0",cs%f0, &
                 do_not_log=.true.)
  call get_param(param_file, mdl,"LENLAT",len_lat, &
                 do_not_log=.true.)
  call get_param(param_file, mdl,"SHELFWAVE_X_WAVELENGTH",cs%Lx, &
                 "Length scale of shelfwave in x-direction.",&
                 units="Same as x,y", default=100.)
  call get_param(param_file, mdl,"SHELFWAVE_Y_LENGTH_SCALE",cs%Ly, &
                 "Length scale of exponential dropoff of topography "//&
                 "in the y-direction.", &
                 units="Same as x,y", default=50.)
  call get_param(param_file, mdl,"SHELFWAVE_Y_MODE",cs%jj, &
                 "Cross-shore wave mode.",               &
                 units="nondim", default=1.)
  cs%alpha = 1. / cs%Ly
  cs%ll = 2. * pi / cs%Lx
  cs%kk = cs%jj * pi / len_lat
  cs%omega = 2 * cs%alpha * cs%f0 * cs%ll / &
             (cs%kk*cs%kk + cs%alpha*cs%alpha + cs%ll*cs%ll)
  register_shelfwave_obc = .true.

end function register_shelfwave_obc

!> Clean up the shelfwave OBC from registry.
subroutine shelfwave_obc_end(CS)
  type(shelfwave_obc_cs), pointer    :: cs         !< shelfwave control structure.

  if (associated(cs)) then
    deallocate(cs)
  endif
end subroutine shelfwave_obc_end

!> Initialization of topography.
subroutine shelfwave_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 :: m_to_z  ! A dimensional rescaling factor.
  integer   :: i, j
  real      :: y, rly, ly, h0

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

  call get_param(param_file, mdl,"SHELFWAVE_Y_LENGTH_SCALE",ly, &
                 default=50., do_not_log=.true.)
  call get_param(param_file, mdl,"MINIMUM_DEPTH", h0, &
                 default=10., units="m", scale=m_to_z, do_not_log=.true.)

  rly = 0. ; if (ly>0.) rly = 1. / ly

  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    ! Compute normalized zonal coordinates (x,y=0 at center of domain)
    y = ( g%geoLatT(i,j) - g%south_lat )
    d(i,j) = h0 * exp(2 * rly * y)
  enddo ; enddo

end subroutine shelfwave_initialize_topography

!> This subroutine sets the properties of flow at open boundary conditions.
subroutine shelfwave_set_obc_data(OBC, CS, G, h, Time)
  type(ocean_obc_type),   pointer    :: obc  !< This open boundary condition type specifies
                                             !! whether, where, and what open boundary
                                             !! conditions are used.
  type(shelfwave_obc_cs), pointer    :: cs   !< tidal bay control structure.
  type(ocean_grid_type),  intent(in) :: g    !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in) :: h !< layer thickness.
  type(time_type),        intent(in) :: time !< model time.

  ! The following variables are used to set up the transport in the shelfwave example.
  real :: my_amp, time_sec
  real :: cos_wt, cos_ky, sin_wt, sin_ky, omega, alpha
  real :: x, y, jj, kk, ll
  character(len=40)  :: mdl = "shelfwave_set_OBC_data" ! This subroutine's name.
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, n
  integer :: isdb, iedb, jsdb, jedb
  type(obc_segment_type), pointer :: segment => null()

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (.not.associated(obc)) return

  time_sec = time_type_to_real(time)
  omega = cs%omega
  alpha = cs%alpha
  my_amp = 1.0
  jj = cs%jj
  kk = cs%kk
  ll = cs%ll
  do n = 1, obc%number_of_segments
    segment => obc%segment(n)
    if (.not. segment%on_pe) cycle
    if (segment%direction /= obc_direction_w) cycle

    isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
    jsd = segment%HI%jsd ; jed = segment%HI%jed
    do j=jsd,jed ; do i=isdb,iedb
      x = g%geoLonCu(i,j) - g%west_lon
      y = g%geoLatCu(i,j) - g%south_lat
      sin_wt = sin(ll*x - omega*time_sec)
      cos_wt = cos(ll*x - omega*time_sec)
      sin_ky = sin(kk * y)
      cos_ky = cos(kk * y)
      segment%normal_vel_bt(i,j) = g%US%m_s_to_L_T*my_amp * exp(- alpha * y) * cos_wt * &
           (alpha * sin_ky + kk * cos_ky)
!     segment%tangential_vel_bt(I,j) = G%US%m_s_to_L_T*my_amp * ll * exp(- alpha * y) * sin_wt * sin_ky
!     segment%vorticity_bt(I,j) = my_amp * exp(- alpha * y) * cos_wt * sin_ky&
!           (ll*ll + kk*kk + alpha*alpha)
    enddo ; enddo
  enddo

end subroutine shelfwave_set_obc_data

end module shelfwave_initialization