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
  logical :: answers_2018    !< If true, use the order of arithmetic and expressions that recover the
                             !! answers from the end of 2018.  Otherwise, use expressions that give
                             !! rotational symmetry and eliminate apparent bugs.
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
  logical :: default_2018_answers
  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
  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.true.)
  call get_param(param_file, mdl, "KELVIN_WAVE_2018_ANSWERS", cs%answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use expressions that give rotational "//&
                 "symmetry and eliminate apparent bugs.", default=default_2018_answers)
  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
          cff = sqrt(gv%g_Earth * 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,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) = us%L_T_to_m_s * (val2 * (val1 * cff * cosa / &
                 (0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j)))) )
        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) = 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) = 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
          if (cs%answers_2018) then
            ! Problem: val2 & cff could be functions of space, but are not set in this loop.
            if (cs%mode == 0) then ; do k=1,nz
              segment%tangential_vel(i,j,k) = us%L_T_to_m_s * (val2 * (val1 * cff * sina / &
                 (0.25 * (g%bathyT(i+1,j) + g%bathyT(i,j) + g%bathyT(i+1,j+1) + g%bathyT(i,j+1))) ))
            enddo ; endif
          else
            cff =sqrt(gv%g_Earth * 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,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) = us%L_T_to_m_s * (val1 * val2 * cff * sina) / &
                 ( 0.25*((g%bathyT(i,j) + g%bathyT(i+1,j+1)) +  (g%bathyT(i+1,j) + g%bathyT(i,j+1))) )
 
            enddo ; endif
          endif
        enddo ; enddo
      endif
    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_offset1) * cosa + y1 * sina
        y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
        if (cs%mode == 0) then
          cff = sqrt(gv%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)) * 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 / &
                 (0.5*(g%bathyT(i+1,j) + g%bathyT(i,j)))) * val2
        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) = 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) = 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
          if (cs%answers_2018) then
            ! Problem: val2 & cff could be functions of space, but are not set in this loop.
            if (cs%mode == 0) then ; do k=1,nz
              segment%tangential_vel(i,j,k) = us%L_T_to_m_s * (val2 * (val1 * cff * sina / &
                 (0.25*(g%bathyT(i+1,j) + g%bathyT(i,j) +  g%bathyT(i+1,j+1) + g%bathyT(i,j+1)))))
            enddo ; endif
          else
            cff = sqrt(gv%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)) * us%m_to_L*y / cff)
            if (cs%mode == 0) then ; do k=1,nz
              segment%tangential_vel(i,j,k) = us%L_T_to_m_s * ((val1 * val2 * cff * sina) / &
                  ( 0.25*((g%bathyT(i,j) + g%bathyT(i+1,j+1)) + (g%bathyT(i+1,j) +  g%bathyT(i,j+1))) ))
            enddo ; endif
          endif
        enddo ; enddo
      endif
    endif
  enddo
 
end subroutine kelvin_set_obc_data
 
end module kelvin_initialization