kelvin_initialization Module Reference

Data Types
type	kelvin_obc_cs
	Control structure for Kelvin wave open boundaries. More...

Functions/Subroutines
logical function, public	register_kelvin_obc (param_file, CS, OBC_Reg)
	Add Kelvin wave to OBC registry. More...
subroutine, public	kelvin_obc_end (CS)
	Clean up the Kelvin wave OBC from registry. More...
subroutine, public	kelvin_initialize_topography (D, G, param_file, max_depth)
	This subroutine sets up the Kelvin topography and land mask. More...
subroutine, public	kelvin_set_obc_data (OBC, CS, G, h, Time)
	This subroutine sets the properties of flow at open boundary conditions. More...

Function/Subroutine Documentation

kelvin_initialize_topography()

subroutine, public kelvin_initialization::kelvin_initialize_topography	(	real, dimension( g %isd: g %ied, g %jsd: g %jed), intent(out)	D,
		type(dyn_horgrid_type), intent(in)	G,
		type(param_file_type), intent(in)	param_file,
		real, intent(in)	max_depth
	)

This subroutine sets up the Kelvin topography and land mask.

Parameters

[in]	g	The dynamic horizontal grid type
[out]	d	Ocean bottom depth in m
[in]	param_file	Parameter file structure
[in]	max_depth	Maximum depth of model in m

Definition at line 124 of file Kelvin_initialization.F90.

References mom_error_handler::mom_mesg().

Referenced by mom_fixed_initialization::mom_initialize_topography().

  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

Here is the call graph for this function:

Here is the caller graph for this function:

kelvin_obc_end()

subroutine, public kelvin_initialization::kelvin_obc_end

(

type(kelvin_obc_cs), pointer

CS

)

Clean up the Kelvin wave OBC from registry.

Parameters

cs	Kelvin wave control structure.

Definition at line 114 of file Kelvin_initialization.F90.

Referenced by mom_boundary_update::obc_register_end().

  type(kelvin_obc_cs), pointer    :: cs         !< Kelvin wave control structure.

  if (associated(cs)) then
    deallocate(cs)
  endif

Here is the caller graph for this function:

kelvin_set_obc_data()

subroutine, public kelvin_initialization::kelvin_set_obc_data	(	type(ocean_obc_type), pointer	OBC,
		type(kelvin_obc_cs), pointer	CS,
		type(ocean_grid_type), intent(in)	G,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		type(time_type), intent(in)	Time
	)

This subroutine sets the properties of flow at open boundary conditions.

Parameters

	obc	This open boundary condition type specifies whether, where, and what open boundary conditions are used.
	cs	Kelvin wave control structure.
[in]	g	The ocean's grid structure.
[in]	h	layer thickness.
[in]	time	model time.

Definition at line 174 of file Kelvin_initialization.F90.

  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

register_kelvin_obc()

logical function, public kelvin_initialization::register_kelvin_obc	(	type(param_file_type), intent(in)	param_file,
		type(kelvin_obc_cs), pointer	CS,
		type(obc_registry_type), pointer	OBC_Reg
	)

Add Kelvin wave to OBC registry.

Parameters

[in]	param_file	parameter file.
	cs	Kelvin wave control structure.
	obc_reg	OBC registry.

Definition at line 64 of file Kelvin_initialization.F90.

References mom_error_handler::mom_error().

  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.

Here is the call graph for this function: