Kelvin_initialization.F90

Go to the documentation of this file.

module kelvin_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_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
use mom_open_boundary,  only : obc_segment_type, register_obc
use mom_open_boundary,  only : obc_direction_n, obc_direction_e
use mom_open_boundary,  only : obc_direction_s, obc_direction_w
use mom_open_boundary,  only : obc_registry_type
use mom_verticalgrid,   only : verticalgrid_type
use mom_time_manager,   only : time_type, set_time, time_type_to_real

implicit none ; private

#include <MOM_memory.h>

public kelvin_set_obc_data, kelvin_initialize_topography
public register_kelvin_obc, kelvin_obc_end

!> Control structure for Kelvin wave open boundaries.
type, public :: kelvin_obc_cs ; private
  integer :: mode = 0         !< Vertical mode
  real    :: coast_angle = 0  !< Angle of coastline
  real    :: coast_offset = 0 !< Longshore distance to coastal angle
  real    :: n0 = 0           !< Brunt-Vaisala frequency
  real    :: h0 = 0           !< Bottom depth
  real    :: f_0              !< Coriolis parameter
  real    :: plx = 0          !< Longshore wave parameter
  real    :: pmz = 0          !< Vertical wave parameter
  real    :: lambda = 0       !< Vertical wave parameter
  real    :: omega            !< Frequency
  real    :: rho_range        !< Density range
  real    :: rho_0            !< Mean density
end type kelvin_obc_cs

! This include declares and sets the variable "version".
#include "version_variable.h"

contains

!> Add Kelvin wave to OBC registry.
function register_kelvin_obc(param_file, CS, OBC_Reg)
  type(param_file_type),    intent(in) :: param_file !< parameter file.
  type(kelvin_obc_cs),      pointer    :: CS         !< Kelvin wave control structure.
  type(obc_registry_type),  pointer    :: OBC_Reg    !< OBC registry.
  logical                              :: register_Kelvin_OBC
  character(len=40)  :: mdl = "register_Kelvin_OBC"  !< This subroutine's name.
  character(len=32)  :: casename = "Kelvin wave"     !< This case's name.
  character(len=200) :: config

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

  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "KELVIN_WAVE_MODE", cs%mode, &
                 "Vertical Kelvin wave mode imposed at upstream open boundary.", &
                 default=0)
  call get_param(param_file, mdl, "F_0", cs%F_0, &
                 default=0.0, do_not_log=.true.)
  call get_param(param_file, mdl, "TOPO_CONFIG", config, do_not_log=.true.)
  if (trim(config) == "Kelvin") then
    call get_param(param_file, mdl, "KELVIN_COAST_OFFSET", cs%coast_offset, &
                   "The distance along the southern and northern boundaries \n"//&
                   "at which the coasts angle in.", &
                   units="km", default=100.0)
    call get_param(param_file, mdl, "KELVIN_COAST_ANGLE", cs%coast_angle, &
                   "The angle of the southern bondary beyond X=KELVIN_COAST_OFFSET.", &
                   units="degrees", default=11.3)
    cs%coast_angle = cs%coast_angle * (atan(1.0)/45.) ! Convert to radians
    cs%coast_offset = cs%coast_offset * 1.e3          ! Convert to m
  endif
  if (cs%mode /= 0) then
    call get_param(param_file, mdl, "DENSITY_RANGE", cs%rho_range, &
                   default=2.0, do_not_log=.true.)
    call get_param(param_file, mdl, "RHO_0", cs%rho_0, &
                   default=1035.0, do_not_log=.true.)
    call get_param(param_file, mdl, "MAXIMUM_DEPTH", cs%H0, &
                   default=1000.0, do_not_log=.true.)
  endif

  ! Register the Kelvin open boundary.
  call register_obc(casename, param_file, obc_reg)
  register_kelvin_obc = .true.

end function register_kelvin_obc

!> Clean up the Kelvin wave OBC from registry.
subroutine kelvin_obc_end(CS)
  type(kelvin_obc_cs), pointer    :: CS         !< Kelvin wave control structure.

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

! -----------------------------------------------------------------------------
!> This subroutine sets up the Kelvin topography and land mask
subroutine kelvin_initialize_topography(D, G, param_file, max_depth)
  type(dyn_horgrid_type),           intent(in)  :: G !< The dynamic horizontal grid type
  real, dimension(SZI_(G),SZJ_(G)), 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
  character(len=40)  :: mdl = "Kelvin_initialize_topography" ! This subroutine's name.
  real :: min_depth ! The minimum and maximum depths in m.
  real :: PI ! 3.1415...
  real :: coast_offset, coast_angle, right_angle
  integer :: i, j

  call mom_mesg("  Kelvin_initialization.F90, Kelvin_initialize_topography: setting topography", 5)

  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0)
  call get_param(param_file, mdl, "KELVIN_COAST_OFFSET", coast_offset, &
                 default=100.0, do_not_log=.true.)
  call get_param(param_file, mdl, "KELVIN_COAST_ANGLE", coast_angle, &
                 default=11.3, do_not_log=.true.)

  coast_angle = coast_angle * (atan(1.0)/45.) ! Convert to radians
  right_angle = 2 * atan(1.0)

  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    d(i,j)=max_depth
    ! Southern side
    if ((g%geoLonT(i,j) > coast_offset) .AND. &
        (atan2(g%geoLatT(i,j),g%geoLonT(i,j) - coast_offset) < &
        coast_angle)) d(i,j)=0.5*min_depth
    ! Northern side
    if ((g%geoLonT(i,j) < g%len_lon - coast_offset) .AND. &
        (atan2(g%len_lat - g%geoLatT(i,j),g%len_lon - coast_offset - g%geoLonT(i,j)) < &
        coast_angle)) d(i,j)=0.5*min_depth
    ! Western side
    if ((g%geoLonT(i,j) < coast_offset) .AND. &
        (atan2(g%geoLatT(i,j),g%geoLonT(i,j) - coast_offset) > &
        right_angle + coast_angle)) d(i,j)=0.5*min_depth
    ! Eastern side
    if ((g%geoLonT(i,j) > g%len_lon - coast_offset) .AND. &
        (atan2(g%len_lat - g%geoLatT(i,j),g%len_lon - coast_offset - g%geoLonT(i,j)) > &
              right_angle + coast_angle)) d(i,j)=0.5*min_depth

    if (d(i,j) > max_depth) d(i,j) = max_depth
    if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
  enddo ; enddo

end subroutine kelvin_initialize_topography

!> This subroutine sets the properties of flow at open boundary conditions.
subroutine kelvin_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(kelvin_obc_cs),    pointer    :: CS   !< Kelvin wave 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 Kelvin example.
  real :: time_sec, cff
  real :: PI
  integer :: i, j, k, n, is, ie, js, je, isd, ied, jsd, jed, nz
  integer :: IsdB, IedB, JsdB, JedB
  real    :: fac, x, y, x1, y1
  real    :: val1, val2, sina, cosa
  type(obc_segment_type), pointer :: segment

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  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)) call mom_error(fatal, 'Kelvin_initialization.F90: '// &
        'Kelvin_set_OBC_data() was called but OBC type was not initialized!')

  time_sec = time_type_to_real(time)
  pi = 4.0*atan(1.0)
  fac = 1.0

  if (cs%mode == 0) then
    cs%omega = 2.0 * pi / (12.42 * 3600.0)      ! M2 Tide period
    val1 = sin(cs%omega * time_sec)
  else
    cs%N0 = sqrt(cs%rho_range / cs%rho_0 * g%g_Earth * cs%H0)
    ! Two wavelengths in domain
    cs%plx = 4.0 * pi / g%len_lon
    cs%pmz = pi * cs%mode / cs%H0
    cs%lambda = cs%pmz * cs%F_0 / cs%N0
    cs%omega = cs%F_0 * cs%plx / cs%lambda
  endif

  sina = sin(cs%coast_angle)
  cosa = cos(cs%coast_angle)
  do n=1,obc%number_of_segments
    segment => obc%segment(n)
    if (.not. segment%on_pe) cycle
    ! Apply values to the inflow end only.
    if (segment%direction == obc_direction_e) cycle
    if (segment%direction == obc_direction_n) cycle

    ! This should be somewhere else...
    segment%Tnudge_in = 1.0/(0.3*86400)

    if (segment%direction == obc_direction_w) then
      isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
      jsd = segment%HI%jsd ; jed = segment%HI%jed
      do j=jsd,jed ; do i=isdb,iedb
        x1 = 1000. * g%geoLonCu(i,j)
        y1 = 1000. * g%geoLatCu(i,j)
        x = (x1 - cs%coast_offset) * cosa + y1 * sina
        y = - (x1 - cs%coast_offset) * sina + y1 * cosa
        if (cs%mode == 0) then
          cff = sqrt(g%g_Earth * 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j)))
          val2 = fac * exp(- cs%F_0 * y / cff)
          segment%eta(i,j) = val2 * cos(cs%omega * time_sec)
          segment%normal_vel_bt(i,j) = val1 * cff * cosa /         &
                 (0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j))) * val2
        else
          segment%eta(i,j) = 0.0
          segment%normal_vel_bt(i,j) = 0.0
          do k=1,nz
            segment%nudged_normal_vel(i,j,k) = fac * cs%lambda / cs%F_0 * &
                   exp(- cs%lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * &
                   cos(cs%omega * time_sec)
          enddo
        endif
      enddo ; enddo
    else
      isd = segment%HI%isd ; ied = segment%HI%ied
      jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
      do j=jsdb,jedb ; do i=isd,ied
        x1 = 1000. * g%geoLonCv(i,j)
        y1 = 1000. * g%geoLatCv(i,j)
        x = (x1 - cs%coast_offset) * cosa + y1 * sina
        y = - (x1 - cs%coast_offset) * sina + y1 * cosa
        if (cs%mode == 0) then
          cff = sqrt(g%g_Earth * 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j)))
          val2 = fac * exp(- 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) * y / cff)
          segment%eta(i,j) = val2 * cos(cs%omega * time_sec)
          segment%normal_vel_bt(i,j) = val1 * cff * sina /       &
                 (0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j))) * val2
        else
          segment%eta(i,j) = 0.0
          segment%normal_vel_bt(i,j) = 0.0
          do k=1,nz
            segment%nudged_normal_vel(i,j,k) = fac * cs%lambda / cs%F_0 * &
                   exp(- cs%lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * cosa
          enddo
        endif
      enddo ; enddo
    endif
  enddo

end subroutine kelvin_set_obc_data

!> \class Kelvin_Initialization
!!
!! The module configures the model for the Kelvin wave experiment.
!! Kelvin = coastally-trapped Kelvin waves from the ROMS examples.
!! Initialize with level surfaces and drive the wave in at the west,
!! radiate out at the east.
end module kelvin_initialization