Kelvin_initialization.F90

!> 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.
module kelvin_initialization

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

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_unit_scaling, only : unit_scale_type
use mom_verticalgrid,   only : verticalgrid_type
use mom_time_manager,   only : time_type, 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

! 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.

!> 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_offset1 = 0 !< Longshore distance to coastal angle
  real    :: coast_offset2 = 0 !< Longshore distance to coastal angle
  real    :: h0 = 0            !< Bottom depth
  real    :: f_0               !< Coriolis parameter
  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.

  ! Local variables
  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, "ROTATED_COAST_OFFSET_1", cs%coast_offset1, &
                   "The distance along the southern and northern boundaries "//&
                   "at which the coasts angle in.", &
                   units="km", default=100.0)
    call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", cs%coast_offset2, &
                   "The distance from the southern and northern boundaries "//&
                   "at which the coasts angle in.", &
                   units="km", default=10.0)
    call get_param(param_file, mdl, "ROTATED_COAST_ANGLE", cs%coast_angle, &
                   "The angle of the southern bondary beyond X=ROTATED_COAST_OFFSET.", &
                   units="degrees", default=11.3)
    cs%coast_angle = cs%coast_angle * (atan(1.0)/45.) ! Convert to radians
    cs%coast_offset1 = cs%coast_offset1 * 1.e3          ! Convert to m
    cs%coast_offset2 = cs%coast_offset2 * 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, 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
  character(len=40)  :: mdl = "Kelvin_initialize_topography" ! This subroutine's name.
  real :: m_to_Z  ! A dimensional rescaling factor.
  real :: min_depth ! The minimum and maximum depths [Z ~> m].
  real :: PI ! 3.1415...
  real :: coast_offset1, coast_offset2, coast_angle, right_angle
  integer :: i, j

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

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

  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
                 "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_z)
  call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_1", coast_offset1, &
                 default=100.0, do_not_log=.true.)
  call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", coast_offset2, &
                 default=10.0, do_not_log=.true.)
  call get_param(param_file, mdl, "ROTATED_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) - g%west_lon > coast_offset1) .AND. &
        (atan2(g%geoLatT(i,j) - g%south_lat + coast_offset2, &
         g%geoLonT(i,j) - g%west_lon - coast_offset1) < coast_angle)) &
      d(i,j) = 0.5*min_depth
    ! Northern side
    if ((g%geoLonT(i,j) - g%west_lon < g%len_lon - coast_offset1) .AND. &
        (atan2(g%len_lat + g%south_lat + coast_offset2 - g%geoLatT(i,j), &
         g%len_lon + g%west_lon - coast_offset1 - g%geoLonT(i,j)) < 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, GV, US, 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.
  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(in) :: h !< layer thickness [H ~> m or kg m-2].
  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 :: N0           ! Brunt-Vaisala frequency [s-1]
  real :: plx          !< Longshore wave parameter
  real :: pmz          !< Vertical wave parameter
  real :: lambda       !< Offshore decay scale
  real :: omega        !< Wave frequency [s-1]
  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 => null()

  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
    omega = 2.0 * pi / (12.42 * 3600.0)      ! M2 Tide period
    val1 = us%m_to_Z * sin(omega * time_sec)
  else
    n0 = us%L_to_m*us%s_to_T * sqrt((cs%rho_range / cs%rho_0) * gv%g_Earth * (us%m_to_Z * cs%H0))
    ! Two wavelengths in domain
    plx = 4.0 * pi / g%len_lon
    pmz = pi * cs%mode / cs%H0
    lambda = pmz * cs%F_0 / n0
    omega = cs%F_0 * plx / lambda

    ! lambda = PI * CS%mode * CS%F_0 / (CS%H0 * N0)
    ! omega = (4.0 * CS%H0 * N0)  / (CS%mode * G%len_lon)
  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%Velocity_nudging_timescale_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
      jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
      do j=jsd,jed ; do i=isdb,iedb
        x1 = 1000. * g%geoLonCu(i,j)
        y1 = 1000. * g%geoLatCu(i,j)
        x = (x1 - cs%coast_offset1) * cosa + y1 * sina
        y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
        if (cs%mode == 0) then
          ! Use inside bathymetry
          cff = sqrt(gv%g_Earth * g%bathyT(i+1,j) )
          val2 = fac * exp(- us%T_to_s*cs%F_0 * us%m_to_L*y / cff)
          segment%eta(i,j) = val2 * cos(omega * time_sec)
          segment%normal_vel_bt(i,j) = (val2 * (val1 * cff * cosa / &
                 (g%bathyT(i+1,j) )) )
          if (segment%nudged) then
            do k=1,nz
              segment%nudged_normal_vel(i,j,k) = (val2 * (val1 * cff * cosa / &
                     (g%bathyT(i+1,j))) )
            enddo
          elseif (segment%specified) then
            do k=1,nz
              segment%normal_vel(i,j,k) = (val2 * (val1 * cff * cosa / &
                     (g%bathyT(i+1,j) )) )
              segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i+1,j,k) * g%dyCu(i,j)
            enddo
          endif
        else
          ! Not rotated yet
          segment%eta(i,j) = 0.0
          segment%normal_vel_bt(i,j) = 0.0
          if (segment%nudged) then
            do k=1,nz
              segment%nudged_normal_vel(i,j,k) = us%m_s_to_L_T * fac * lambda / cs%F_0 * &
                   exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * &
                   cos(omega * time_sec)
            enddo
          elseif (segment%specified) then
            do k=1,nz
              segment%normal_vel(i,j,k) = us%m_s_to_L_T * fac * lambda / cs%F_0 * &
                   exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * &
                   cos(omega * time_sec)
              segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i+1,j,k) * g%dyCu(i,j)
            enddo
          endif
        endif
      enddo ; enddo
      if (associated(segment%tangential_vel)) then
        do j=jsdb+1,jedb-1 ; do i=isdb,iedb
          x1 = 1000. * g%geoLonBu(i,j)
          y1 = 1000. * g%geoLatBu(i,j)
          x = (x1 - cs%coast_offset1) * cosa + y1 * sina
          y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
          cff =sqrt(gv%g_Earth * g%bathyT(i+1,j) )
          val2 = fac * exp(- us%T_to_s*cs%F_0 * us%m_to_L*y / cff)
          if (cs%mode == 0) then ; do k=1,nz
            segment%tangential_vel(i,j,k) = (val1 * val2 * cff * sina) / &
               ( 0.5*(g%bathyT(i+1,j+1) +  g%bathyT(i+1,j) ) )

          enddo ; endif
        enddo ; enddo
      endif
    else ! Must be south
      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_offset1) * cosa + y1 * sina
        y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
        if (cs%mode == 0) then
          cff = sqrt(gv%g_Earth * g%bathyT(i,j+1) )
          val2 = fac * exp(- 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) * us%m_to_L*y / cff)
          segment%eta(i,j) = val2 * cos(omega * time_sec)
          segment%normal_vel_bt(i,j) = us%L_T_to_m_s * (val1 * cff * sina / &
                 (g%bathyT(i,j+1) )) * val2
          if (segment%nudged) then
            do k=1,nz
              segment%nudged_normal_vel(i,j,k) = us%L_T_to_m_s * (val1 * cff * sina / &
                     (g%bathyT(i,j+1) )) * val2
            enddo
          elseif (segment%specified) then
            do k=1,nz
              segment%normal_vel(i,j,k) = us%L_T_to_m_s * (val1 * cff * sina / &
                     (g%bathyT(i,j+1) )) * val2
              segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i,j+1,k) * g%dxCv(i,j)
            enddo
          endif
        else
          ! Not rotated yet
          segment%eta(i,j) = 0.0
          segment%normal_vel_bt(i,j) = 0.0
          if (segment%nudged) then
            do k=1,nz
              segment%nudged_normal_vel(i,j,k) = us%m_s_to_L_T*fac * lambda / cs%F_0 * &
                   exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * cosa
            enddo
          elseif (segment%specified) then
            do k=1,nz
              segment%normal_vel(i,j,k) = us%m_s_to_L_T*fac * lambda / cs%F_0 * &
                   exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * cosa
              segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i,j+1,k) * g%dxCv(i,j)
            enddo
          endif
        endif
      enddo ; enddo
      if (associated(segment%tangential_vel)) then
        do j=jsdb,jedb ; do i=isdb+1,iedb-1
          x1 = 1000. * g%geoLonBu(i,j)
          y1 = 1000. * g%geoLatBu(i,j)
          x = (x1 - cs%coast_offset1) * cosa + y1 * sina
          y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
          cff = sqrt(gv%g_Earth * g%bathyT(i,j+1) )
          val2 = fac * exp(- 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) * us%m_to_L*y / cff)
          if (cs%mode == 0) then ; do k=1,nz
            segment%tangential_vel(i,j,k) = ((val1 * val2 * cff * sina) / &
                ( 0.5*((g%bathyT(i+1,j+1)) + g%bathyT(i,j+1))) )
          enddo ; endif
        enddo ; enddo
      endif
    endif
  enddo

end subroutine kelvin_set_obc_data

end module kelvin_initialization