MOM6
Kelvin_initialization.F90
Go to the documentation of this file.
1 !> Configures the model for the Kelvin wave experiment.
2 !!
3 !! Kelvin = coastally-trapped Kelvin waves from the ROMS examples.
4 !! Initialize with level surfaces and drive the wave in at the west,
5 !! radiate out at the east.
7 
8 ! This file is part of MOM6. See LICENSE.md for the license.
9 
11 use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
13 use mom_grid, only : ocean_grid_type
21 use mom_time_manager, only : time_type, time_type_to_real
22 
23 implicit none ; private
24 
25 #include <MOM_memory.h>
26 
29 
30 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
31 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
32 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
33 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
34 
35 !> Control structure for Kelvin wave open boundaries.
36 type, public :: kelvin_obc_cs ; private
37  integer :: mode = 0 !< Vertical mode
38  real :: coast_angle = 0 !< Angle of coastline
39  real :: coast_offset1 = 0 !< Longshore distance to coastal angle
40  real :: coast_offset2 = 0 !< Longshore distance to coastal angle
41  real :: h0 = 0 !< Bottom depth
42  real :: f_0 !< Coriolis parameter
43  real :: rho_range !< Density range
44  real :: rho_0 !< Mean density
45  logical :: answers_2018 !< If true, use the order of arithmetic and expressions that recover the
46  !! answers from the end of 2018. Otherwise, use expressions that give
47  !! rotational symmetry and eliminate apparent bugs.
48 end type kelvin_obc_cs
49 
50 ! This include declares and sets the variable "version".
51 #include "version_variable.h"
52 
53 contains
54 
55 !> Add Kelvin wave to OBC registry.
56 function register_kelvin_obc(param_file, CS, OBC_Reg)
57  type(param_file_type), intent(in) :: param_file !< parameter file.
58  type(kelvin_obc_cs), pointer :: cs !< Kelvin wave control structure.
59  type(obc_registry_type), pointer :: obc_reg !< OBC registry.
60 
61  ! Local variables
62  logical :: register_kelvin_obc
63  logical :: default_2018_answers
64  character(len=40) :: mdl = "register_Kelvin_OBC" !< This subroutine's name.
65  character(len=32) :: casename = "Kelvin wave" !< This case's name.
66  character(len=200) :: config
67 
68  if (associated(cs)) then
69  call mom_error(warning, "register_Kelvin_OBC called with an "// &
70  "associated control structure.")
71  return
72  endif
73  allocate(cs)
74 
75  call log_version(param_file, mdl, version, "")
76  call get_param(param_file, mdl, "KELVIN_WAVE_MODE", cs%mode, &
77  "Vertical Kelvin wave mode imposed at upstream open boundary.", &
78  default=0)
79  call get_param(param_file, mdl, "F_0", cs%F_0, &
80  default=0.0, do_not_log=.true.)
81  call get_param(param_file, mdl, "TOPO_CONFIG", config, do_not_log=.true.)
82  if (trim(config) == "Kelvin") then
83  call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_1", cs%coast_offset1, &
84  "The distance along the southern and northern boundaries "//&
85  "at which the coasts angle in.", &
86  units="km", default=100.0)
87  call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", cs%coast_offset2, &
88  "The distance from the southern and northern boundaries "//&
89  "at which the coasts angle in.", &
90  units="km", default=10.0)
91  call get_param(param_file, mdl, "ROTATED_COAST_ANGLE", cs%coast_angle, &
92  "The angle of the southern bondary beyond X=ROTATED_COAST_OFFSET.", &
93  units="degrees", default=11.3)
94  cs%coast_angle = cs%coast_angle * (atan(1.0)/45.) ! Convert to radians
95  cs%coast_offset1 = cs%coast_offset1 * 1.e3 ! Convert to m
96  cs%coast_offset2 = cs%coast_offset2 * 1.e3 ! Convert to m
97  endif
98  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
99  "This sets the default value for the various _2018_ANSWERS parameters.", &
100  default=.true.)
101  call get_param(param_file, mdl, "KELVIN_WAVE_2018_ANSWERS", cs%answers_2018, &
102  "If true, use the order of arithmetic and expressions that recover the "//&
103  "answers from the end of 2018. Otherwise, use expressions that give rotational "//&
104  "symmetry and eliminate apparent bugs.", default=default_2018_answers)
105  if (cs%mode /= 0) then
106  call get_param(param_file, mdl, "DENSITY_RANGE", cs%rho_range, &
107  default=2.0, do_not_log=.true.)
108  call get_param(param_file, mdl, "RHO_0", cs%rho_0, &
109  default=1035.0, do_not_log=.true.)
110  call get_param(param_file, mdl, "MAXIMUM_DEPTH", cs%H0, &
111  default=1000.0, do_not_log=.true.)
112  endif
113 
114  ! Register the Kelvin open boundary.
115  call register_obc(casename, param_file, obc_reg)
116  register_kelvin_obc = .true.
117 
118 end function register_kelvin_obc
119 
120 !> Clean up the Kelvin wave OBC from registry.
121 subroutine kelvin_obc_end(CS)
122  type(kelvin_obc_cs), pointer :: cs !< Kelvin wave control structure.
123 
124  if (associated(cs)) then
125  deallocate(cs)
126  endif
127 end subroutine kelvin_obc_end
128 
129 ! -----------------------------------------------------------------------------
130 !> This subroutine sets up the Kelvin topography and land mask
131 subroutine kelvin_initialize_topography(D, G, param_file, max_depth, US)
132  type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
133  real, dimension(G%isd:G%ied,G%jsd:G%jed), &
134  intent(out) :: d !< Ocean bottom depth in m or Z if US is present
135  type(param_file_type), intent(in) :: param_file !< Parameter file structure
136  real, intent(in) :: max_depth !< Maximum model depth in the units of D
137  type(unit_scale_type), optional, intent(in) :: us !< A dimensional unit scaling type
138 
139  ! Local variables
140  character(len=40) :: mdl = "Kelvin_initialize_topography" ! This subroutine's name.
141  real :: m_to_z ! A dimensional rescaling factor.
142  real :: min_depth ! The minimum and maximum depths [Z ~> m].
143  real :: pi ! 3.1415...
144  real :: coast_offset1, coast_offset2, coast_angle, right_angle
145  integer :: i, j
146 
147  call mom_mesg(" Kelvin_initialization.F90, Kelvin_initialize_topography: setting topography", 5)
148 
149  m_to_z = 1.0 ; if (present(us)) m_to_z = us%m_to_Z
150 
151  call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
152  "The minimum depth of the ocean.", units="m", default=0.0, scale=m_to_z)
153  call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_1", coast_offset1, &
154  default=100.0, do_not_log=.true.)
155  call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", coast_offset2, &
156  default=10.0, do_not_log=.true.)
157  call get_param(param_file, mdl, "ROTATED_COAST_ANGLE", coast_angle, &
158  default=11.3, do_not_log=.true.)
159 
160  coast_angle = coast_angle * (atan(1.0)/45.) ! Convert to radians
161  right_angle = 2 * atan(1.0)
162 
163  do j=g%jsc,g%jec ; do i=g%isc,g%iec
164  d(i,j) = max_depth
165  ! Southern side
166  if ((g%geoLonT(i,j) - g%west_lon > coast_offset1) .AND. &
167  (atan2(g%geoLatT(i,j) - g%south_lat + coast_offset2, &
168  g%geoLonT(i,j) - g%west_lon - coast_offset1) < coast_angle)) &
169  d(i,j) = 0.5*min_depth
170  ! Northern side
171  if ((g%geoLonT(i,j) - g%west_lon < g%len_lon - coast_offset1) .AND. &
172  (atan2(g%len_lat + g%south_lat + coast_offset2 - g%geoLatT(i,j), &
173  g%len_lon + g%west_lon - coast_offset1 - g%geoLonT(i,j)) < coast_angle)) &
174  d(i,j) = 0.5*min_depth
175 
176  if (d(i,j) > max_depth) d(i,j) = max_depth
177  if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
178  enddo ; enddo
179 
180 end subroutine kelvin_initialize_topography
181 
182 !> This subroutine sets the properties of flow at open boundary conditions.
183 subroutine kelvin_set_obc_data(OBC, CS, G, GV, US, h, Time)
184  type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
185  !! whether, where, and what open boundary
186  !! conditions are used.
187  type(kelvin_obc_cs), pointer :: cs !< Kelvin wave control structure.
188  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
189  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
190  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
191  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< layer thickness [H ~> m or kg m-2].
192  type(time_type), intent(in) :: time !< model time.
193 
194  ! The following variables are used to set up the transport in the Kelvin example.
195  real :: time_sec, cff
196  real :: n0 ! Brunt-Vaisala frequency [s-1]
197  real :: plx !< Longshore wave parameter
198  real :: pmz !< Vertical wave parameter
199  real :: lambda !< Offshore decay scale
200  real :: omega !< Wave frequency [s-1]
201  real :: pi
202  integer :: i, j, k, n, is, ie, js, je, isd, ied, jsd, jed, nz
203  integer :: isdb, iedb, jsdb, jedb
204  real :: fac, x, y, x1, y1
205  real :: val1, val2, sina, cosa
206  type(obc_segment_type), pointer :: segment => null()
207 
208  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
209  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
210  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
211 
212  if (.not.associated(obc)) call mom_error(fatal, 'Kelvin_initialization.F90: '// &
213  'Kelvin_set_OBC_data() was called but OBC type was not initialized!')
214 
215  time_sec = time_type_to_real(time)
216  pi = 4.0*atan(1.0)
217  fac = 1.0
218 
219  if (cs%mode == 0) then
220  omega = 2.0 * pi / (12.42 * 3600.0) ! M2 Tide period
221  val1 = us%m_to_Z * sin(omega * time_sec)
222  else
223  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))
224  ! Two wavelengths in domain
225  plx = 4.0 * pi / g%len_lon
226  pmz = pi * cs%mode / cs%H0
227  lambda = pmz * cs%F_0 / n0
228  omega = cs%F_0 * plx / lambda
229 
230  ! lambda = PI * CS%mode * CS%F_0 / (CS%H0 * N0)
231  ! omega = (4.0 * CS%H0 * N0) / (CS%mode * G%len_lon)
232  endif
233 
234  sina = sin(cs%coast_angle)
235  cosa = cos(cs%coast_angle)
236  do n=1,obc%number_of_segments
237  segment => obc%segment(n)
238  if (.not. segment%on_pe) cycle
239  ! Apply values to the inflow end only.
240  if (segment%direction == obc_direction_e) cycle
241  if (segment%direction == obc_direction_n) cycle
242 
243  ! This should be somewhere else...
244  segment%Velocity_nudging_timescale_in = 1.0/(0.3*86400)
245 
246  if (segment%direction == obc_direction_w) then
247  isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
248  jsd = segment%HI%jsd ; jed = segment%HI%jed
249  jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
250  do j=jsd,jed ; do i=isdb,iedb
251  x1 = 1000. * g%geoLonCu(i,j)
252  y1 = 1000. * g%geoLatCu(i,j)
253  x = (x1 - cs%coast_offset1) * cosa + y1 * sina
254  y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
255  if (cs%mode == 0) then
256  cff = sqrt(gv%g_Earth * 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j)))
257  val2 = fac * exp(- us%T_to_s*cs%F_0 * us%m_to_L*y / cff)
258  segment%eta(i,j) = val2 * cos(omega * time_sec)
259  segment%normal_vel_bt(i,j) = (val2 * (val1 * cff * cosa / &
260  (0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j)))) )
261  else
262  ! Not rotated yet
263  segment%eta(i,j) = 0.0
264  segment%normal_vel_bt(i,j) = 0.0
265  if (segment%nudged) then
266  do k=1,nz
267  segment%nudged_normal_vel(i,j,k) = us%m_s_to_L_T * fac * lambda / cs%F_0 * &
268  exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * &
269  cos(omega * time_sec)
270  enddo
271  elseif (segment%specified) then
272  do k=1,nz
273  segment%normal_vel(i,j,k) = us%m_s_to_L_T * fac * lambda / cs%F_0 * &
274  exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * &
275  cos(omega * time_sec)
276  segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i+1,j,k) * g%dyCu(i,j)
277  enddo
278  endif
279  endif
280  enddo ; enddo
281  if (associated(segment%tangential_vel)) then
282  do j=jsdb+1,jedb-1 ; do i=isdb,iedb
283  x1 = 1000. * g%geoLonBu(i,j)
284  y1 = 1000. * g%geoLatBu(i,j)
285  x = (x1 - cs%coast_offset1) * cosa + y1 * sina
286  y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
287  if (cs%answers_2018) then
288  ! Problem: val2 & cff could be functions of space, but are not set in this loop.
289  if (cs%mode == 0) then ; do k=1,nz
290  segment%tangential_vel(i,j,k) = (val2 * (val1 * cff * sina / &
291  (0.25 * (g%bathyT(i+1,j) + g%bathyT(i,j) + g%bathyT(i+1,j+1) + g%bathyT(i,j+1))) ))
292  enddo ; endif
293  else
294  cff =sqrt(gv%g_Earth * 0.5 * (g%bathyT(i+1,j) + g%bathyT(i,j)))
295  val2 = fac * exp(- us%T_to_s*cs%F_0 * us%m_to_L*y / cff)
296  if (cs%mode == 0) then ; do k=1,nz
297  segment%tangential_vel(i,j,k) = (val1 * val2 * cff * sina) / &
298  ( 0.25*((g%bathyT(i,j) + g%bathyT(i+1,j+1)) + (g%bathyT(i+1,j) + g%bathyT(i,j+1))) )
299 
300  enddo ; endif
301  endif
302  enddo ; enddo
303  endif
304  else
305  isd = segment%HI%isd ; ied = segment%HI%ied
306  jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
307  do j=jsdb,jedb ; do i=isd,ied
308  x1 = 1000. * g%geoLonCv(i,j)
309  y1 = 1000. * g%geoLatCv(i,j)
310  x = (x1 - cs%coast_offset1) * cosa + y1 * sina
311  y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
312  if (cs%mode == 0) then
313  cff = sqrt(gv%g_Earth * 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j)))
314  val2 = fac * exp(- 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) * us%m_to_L*y / cff)
315  segment%eta(i,j) = val2 * cos(omega * time_sec)
316  segment%normal_vel_bt(i,j) = us%L_T_to_m_s * (val1 * cff * sina / &
317  (0.5*(g%bathyT(i+1,j) + g%bathyT(i,j)))) * val2
318  else
319  ! Not rotated yet
320  segment%eta(i,j) = 0.0
321  segment%normal_vel_bt(i,j) = 0.0
322  if (segment%nudged) then
323  do k=1,nz
324  segment%nudged_normal_vel(i,j,k) = us%m_s_to_L_T*fac * lambda / cs%F_0 * &
325  exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * cosa
326  enddo
327  elseif (segment%specified) then
328  do k=1,nz
329  segment%normal_vel(i,j,k) = us%m_s_to_L_T*fac * lambda / cs%F_0 * &
330  exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * cosa
331  segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i,j+1,k) * g%dxCv(i,j)
332  enddo
333  endif
334  endif
335  enddo ; enddo
336  if (associated(segment%tangential_vel)) then
337  do j=jsdb,jedb ; do i=isdb+1,iedb-1
338  x1 = 1000. * g%geoLonBu(i,j)
339  y1 = 1000. * g%geoLatBu(i,j)
340  x = (x1 - cs%coast_offset1) * cosa + y1 * sina
341  y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
342  if (cs%answers_2018) then
343  ! Problem: val2 & cff could be functions of space, but are not set in this loop.
344  if (cs%mode == 0) then ; do k=1,nz
345  segment%tangential_vel(i,j,k) = (val2 * (val1 * cff * sina / &
346  (0.25*(g%bathyT(i+1,j) + g%bathyT(i,j) + g%bathyT(i+1,j+1) + g%bathyT(i,j+1)))))
347  enddo ; endif
348  else
349  cff = sqrt(gv%g_Earth * 0.5 * (g%bathyT(i,j+1) + g%bathyT(i,j)))
350  val2 = fac * exp(- 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) * us%m_to_L*y / cff)
351  if (cs%mode == 0) then ; do k=1,nz
352  segment%tangential_vel(i,j,k) = ((val1 * val2 * cff * sina) / &
353  ( 0.25*((g%bathyT(i,j) + g%bathyT(i+1,j+1)) + (g%bathyT(i+1,j) + g%bathyT(i,j+1))) ))
354  enddo ; endif
355  endif
356  enddo ; enddo
357  endif
358  endif
359  enddo
360 
361 end subroutine kelvin_set_obc_data
362 
363 end module kelvin_initialization
mom_time_manager
Wraps the FMS time manager functions.
Definition: MOM_time_manager.F90:2
kelvin_initialization::kelvin_set_obc_data
subroutine, public kelvin_set_obc_data(OBC, CS, G, GV, US, h, Time)
This subroutine sets the properties of flow at open boundary conditions.
Definition: Kelvin_initialization.F90:184
mom_open_boundary::obc_direction_n
integer, parameter, public obc_direction_n
Indicates the boundary is an effective northern boundary.
Definition: MOM_open_boundary.F90:63
mom_open_boundary::obc_direction_s
integer, parameter, public obc_direction_s
Indicates the boundary is an effective southern boundary.
Definition: MOM_open_boundary.F90:64
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_open_boundary::obc_direction_e
integer, parameter, public obc_direction_e
Indicates the boundary is an effective eastern boundary.
Definition: MOM_open_boundary.F90:65
mom_open_boundary::obc_direction_w
integer, parameter, public obc_direction_w
Indicates the boundary is an effective western boundary.
Definition: MOM_open_boundary.F90:66
mom_dyn_horgrid
Contains a shareable dynamic type for describing horizontal grids and metric data and utilty routines...
Definition: MOM_dyn_horgrid.F90:3
kelvin_initialization
Configures the model for the Kelvin wave experiment.
Definition: Kelvin_initialization.F90:6
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
kelvin_initialization::kelvin_obc_cs
Control structure for Kelvin wave open boundaries.
Definition: Kelvin_initialization.F90:36
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
kelvin_initialization::register_kelvin_obc
logical function, public register_kelvin_obc(param_file, CS, OBC_Reg)
Add Kelvin wave to OBC registry.
Definition: Kelvin_initialization.F90:57
mom_open_boundary::obc_registry_type
Type to carry basic OBC information needed for updating values.
Definition: MOM_open_boundary.F90:294
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_open_boundary::register_obc
subroutine, public register_obc(name, param_file, Reg)
register open boundary objects for boundary updates.
Definition: MOM_open_boundary.F90:3899
mom_open_boundary
Controls where open boundary conditions are applied.
Definition: MOM_open_boundary.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
kelvin_initialization::kelvin_initialize_topography
subroutine, public kelvin_initialize_topography(D, G, param_file, max_depth, US)
This subroutine sets up the Kelvin topography and land mask.
Definition: Kelvin_initialization.F90:132
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_open_boundary::ocean_obc_type
Open-boundary data.
Definition: MOM_open_boundary.F90:195
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_open_boundary::obc_none
integer, parameter, public obc_none
Indicates the use of no open boundary.
Definition: MOM_open_boundary.F90:58
mom_open_boundary::obc_segment_type
Open boundary segment data structure.
Definition: MOM_open_boundary.F90:103
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
kelvin_initialization::kelvin_obc_end
subroutine, public kelvin_obc_end(CS)
Clean up the Kelvin wave OBC from registry.
Definition: Kelvin_initialization.F90:122
mom_dyn_horgrid::dyn_horgrid_type
Describes the horizontal ocean grid with only dynamic memory arrays.
Definition: MOM_dyn_horgrid.F90:23
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26