kelvin_initialization Module Reference

Detailed Description

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.

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, US)
	This subroutine sets up the Kelvin topography and land mask. More...
subroutine, public	kelvin_set_obc_data (OBC, CS, G, GV, US, 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,
		type(unit_scale_type), intent(in), optional	US
	)

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 or Z if US is present
[in]	param_file	Parameter file structure
[in]	max_depth	Maximum model depth in the units of D
[in]	us	A dimensional unit scaling type

Definition at line 132 of file Kelvin_initialization.F90.

  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
 

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 122 of file Kelvin_initialization.F90.

  type(Kelvin_OBC_CS), pointer    :: CS         !< Kelvin wave control structure.
 
  if (associated(cs)) then
    deallocate(cs)
  endif

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,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		real, dimension(szi_(g),szj_(g),szk_(g)), 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]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	h	layer thickness [H ~> m or kg m-2].
[in]	time	model time.

Definition at line 184 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.
  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
 

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 57 of file Kelvin_initialization.F90.

  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.