tidal_bay_initialization.F90

!> Configures the model for the "tidal_bay" experiment.
!! tidal_bay = Tidally resonant bay from Zygmunt Kowalik's class on tides.
module tidal_bay_initialization

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

use mom_coms,           only : reproducing_sum
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_registry_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 tidal_bay_set_obc_data, tidal_bay_obc_end
public register_tidal_bay_obc

!> Control structure for tidal bay open boundaries.
type, public :: tidal_bay_obc_cs ; private
  real :: tide_flow = 3.0e6         !< Maximum tidal flux.
end type tidal_bay_obc_cs

contains

!> Add tidal bay to OBC registry.
function register_tidal_bay_obc(param_file, CS, OBC_Reg)
  type(param_file_type),    intent(in) :: param_file !< parameter file.
  type(tidal_bay_obc_cs),   pointer    :: cs         !< tidal bay control structure.
  type(obc_registry_type),  pointer    :: obc_reg    !< OBC registry.
  logical                              :: register_tidal_bay_obc
  character(len=32)  :: casename = "tidal bay"       !< This case's name.

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

  ! Register the open boundaries.
  call register_obc(casename, param_file, obc_reg)
  register_tidal_bay_obc = .true.

end function register_tidal_bay_obc

!> Clean up the tidal bay OBC from registry.
subroutine tidal_bay_obc_end(CS)
  type(tidal_bay_obc_cs), pointer    :: cs         !< tidal bay control structure.

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

!> This subroutine sets the properties of flow at open boundary conditions.
subroutine tidal_bay_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(tidal_bay_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 :: time_sec, cff
  real :: my_flux, total_area
  real :: pi
  real, allocatable :: my_area(:,:)
  character(len=40)  :: mdl = "tidal_bay_set_OBC_data" ! This subroutine's name.
  integer :: i, j, k, itt, is, ie, js, je, isd, ied, jsd, jed, nz, 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 ; 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

  pi = 4.0*atan(1.0)

  if (.not.associated(obc)) return

  allocate(my_area(1:1,js:je))

  time_sec = time_type_to_real(time)
  cff = 0.1*sin(2.0*pi*time_sec/(12.0*3600.0))
  my_area=0.0
  my_flux=0.0
  segment => obc%segment(1)

  do j=segment%HI%jsc,segment%HI%jec ; do i=segment%HI%IscB,segment%HI%IecB
    if (obc%segnum_u(i,j) /= obc_none) then
      do k=1,nz
        my_area(1,j) = my_area(1,j) + h(i,j,k)*g%US%L_to_m*g%dyCu(i,j)
      enddo
    endif
  enddo ; enddo
  total_area = reproducing_sum(my_area)
  my_flux = - cs%tide_flow*sin(2.0*pi*time_sec/(12.0*3600.0))

  do n = 1, obc%number_of_segments
    segment => obc%segment(n)

    if (.not. segment%on_pe) cycle

    segment%normal_vel_bt(:,:) = g%US%m_s_to_L_T*my_flux/total_area
    segment%eta(:,:) = cff

  enddo ! end segment loop

end subroutine tidal_bay_set_obc_data

end module tidal_bay_initialization