shelfwave_initialization.F90

Go to the documentation of this file.

module shelfwave_initialization
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

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, set_time, time_type_to_real

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\n"//&
                 "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 )
  ! Arguments
  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
  type(param_file_type),              intent(in)  :: param_file !< Parameter file structure
  real,                               intent(in)  :: max_depth  !< Maximum depth of model in m

  ! Local variables
  integer   :: i, j
  real      :: y, rLy, Ly, H0

  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., do_not_log=.true.)

  rly = 0. ; if (ly>0.) rly = 1. / ly
  do i=g%isc,g%iec
    do j=g%jsc,g%jec
      ! 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 tidal_bay 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

  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) = my_amp * exp(- alpha * y) * cos_wt * &
           (alpha * sin_ky + kk * cos_ky)
!     segment%tangential_vel_bt(I,j) = 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

!> \namespace shelfwave_initialization
!!
!! The module configures the model for the idealized shelfwave
!! test case.
end module shelfwave_initialization