Kelvin_initialization.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Configures the model for the Kelvin wave experiment.
6!!
7!! Kelvin = coastally-trapped Kelvin waves from the ROMS examples.
8!! Initialize with level surfaces and drive the wave in at the west,
9!! radiate out at the east.
11
13use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
14use mom_file_parser, only : get_param, log_version, param_file_type
15use mom_grid, only : ocean_grid_type
16use mom_open_boundary, only : ocean_obc_type, obc_none
18use mom_open_boundary, only : obc_direction_n, obc_direction_e, obc_direction_s, obc_direction_w
22use mom_time_manager, only : time_type, time_to_real
23
24implicit none ; private
25
26#include <MOM_memory.h>
27
30
31! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
32! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
33! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
34! vary with the Boussinesq approximation, the Boussinesq variant is given first.
35
36!> Control structure for Kelvin wave open boundaries.
37type, public :: kelvin_obc_cs ; private
38 integer :: mode = 0 !< Vertical mode
39 real :: coast_angle = 0 !< Angle of coastline [rad]
40 real :: coast_offset1 = 0 !< Longshore distance to coastal angle [L ~> m]
41 real :: coast_offset2 = 0 !< Offshore distance to coastal angle [L ~> m]
42 real :: h0 = 0 !< Bottom depth [Z ~> m]
43 real :: f_0 !< Coriolis parameter [T-1 ~> s-1]
44 real :: rho_range !< Density range [R ~> kg m-3]
45 real :: rho_0 !< Mean density [R ~> kg m-3]
46 real :: wave_period !< Period of the mode-0 waves [T ~> s]
47 real :: ssh_amp !< Amplitude of the sea surface height forcing for mode-0 waves [Z ~> m]
48 real :: inflow_amp !< Amplitude of the boundary velocity forcing for internal waves [L T-1 ~> m s-1]
49 real :: obc_nudging_time !< The timescale with which the inflowing open boundary velocities are nudged toward
50 !! their intended values with the Kelvin wave test case [T ~> s], or a negative
51 !! value to retain the value that is set when the OBC segments are initialized.
52 logical :: indexing_bugs !< If true, retain several horizontal indexing bugs that were in the
53 !! original version of Kelvin_set_OBC_data.
54end type kelvin_obc_cs
55
56! This include declares and sets the variable "version".
57#include "version_variable.h"
58
59contains
60
61!> Add Kelvin wave to OBC registry.
62logical function register_kelvin_obc(param_file, CS, US, OBC_Reg)
63 type(param_file_type), intent(in) :: param_file !< parameter file.
64 type(kelvin_obc_cs), pointer :: cs !< Kelvin wave control structure.
65 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
66 type(obc_registry_type), pointer :: obc_reg !< OBC registry.
67
68 ! Local variables
69 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
70 ! recreate the bugs, or if false bugs are only used if actively selected.
71 character(len=40) :: mdl = "register_Kelvin_OBC" !< This subroutine's name.
72 character(len=32) :: casename = "Kelvin wave" !< This case's name.
73 character(len=200) :: config
74
75 if (associated(cs)) then
76 call mom_error(warning, "register_Kelvin_OBC called with an "// &
77 "associated control structure.")
78 return
79 endif
80 allocate(cs)
81
82 call log_version(param_file, mdl, version, "")
83 call get_param(param_file, mdl, "KELVIN_WAVE_MODE", cs%mode, &
84 "Vertical Kelvin wave mode imposed at upstream open boundary.", &
85 default=0)
86 call get_param(param_file, mdl, "F_0", cs%F_0, &
87 default=0.0, units="s-1", scale=us%T_to_s, do_not_log=.true.)
88 call get_param(param_file, mdl, "TOPO_CONFIG", config, fail_if_missing=.true., do_not_log=.true.)
89 if (trim(config) == "Kelvin") then
90 call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_1", cs%coast_offset1, &
91 "The distance along the southern and northern boundaries "//&
92 "at which the coasts angle in.", &
93 units="km", default=100.0, scale=1.0e3*us%m_to_L)
94 call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", cs%coast_offset2, &
95 "The distance from the southern and northern boundaries "//&
96 "at which the coasts angle in.", &
97 units="km", default=10.0, scale=1.0e3*us%m_to_L)
98 call get_param(param_file, mdl, "ROTATED_COAST_ANGLE", cs%coast_angle, &
99 "The angle of the southern bondary beyond X=ROTATED_COAST_OFFSET.", &
100 units="degrees", default=11.3, scale=atan(1.0)/45.) ! Convert to radians
101 else
102 cs%coast_offset1 = 0.0 ; cs%coast_offset2 = 0.0 ; cs%coast_angle = 0.0
103 endif
104 if (cs%mode == 0) then
105 call get_param(param_file, mdl, "KELVIN_WAVE_PERIOD", cs%wave_period, &
106 "The period of the Kelvin wave forcing at the open boundaries. "//&
107 "The default value is the M2 tide period.", &
108 units="s", default=12.42*3600.0, scale=us%s_to_T)
109 call get_param(param_file, mdl, "KELVIN_WAVE_SSH_AMP", cs%ssh_amp, &
110 "The amplitude of the Kelvin wave sea surface height anomaly forcing "//&
111 "at the open boundaries.", units="m", default=1.0, scale=us%m_to_Z)
112 else
113 call get_param(param_file, mdl, "DENSITY_RANGE", cs%rho_range, &
114 units="kg m-3", default=2.0, scale=us%kg_m3_to_R, do_not_log=.true.)
115 call get_param(param_file, mdl, "RHO_0", cs%rho_0, &
116 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R, do_not_log=.true.)
117 call get_param(param_file, mdl, "MAXIMUM_DEPTH", cs%H0, &
118 units="m", default=1000.0, scale=us%m_to_Z, do_not_log=.true.)
119 call get_param(param_file, mdl, "KELVIN_WAVE_INFLOW_AMP", cs%inflow_amp, &
120 "The amplitude of the Kelvin wave sea surface inflow velocity forcing "//&
121 "at the open boundaries.", units="m s-1", default=1.0, scale=us%m_s_to_L_T)
122 endif
123
124 call get_param(param_file, mdl, "KELVIN_WAVE_VEL_NUDGING_TIMESCALE", cs%OBC_nudging_time, &
125 "The timescale with which the inflowing open boundary velocities are nudged toward "//&
126 "their intended values with the Kelvin wave test case, or a negative value to keep "//&
127 "the value that is set when the OBC segments are initialized.", &
128 units="s", default=-1.0, scale=us%s_to_T)
129 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
130 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
131 call get_param(param_file, mdl, "KELVIN_SET_OBC_INDEXING_BUGS", cs%indexing_bugs, &
132 "If true, retain several horizontal indexing bugs that were in the original "//&
133 "version of Kelvin_set_OBC_data.", default=enable_bugs)
134
135 ! Register the Kelvin open boundary.
136 call register_obc(casename, param_file, obc_reg)
137 register_kelvin_obc = .true.
138
139 ! TODO: Revisit and correct the internal Kelvin wave test case.
140 ! Specifically, using wave_speed() and investigating adding eta_anom
141 ! noted in the comments below.
142
143end function register_kelvin_obc
144
145!> Clean up the Kelvin wave OBC from registry.
146subroutine kelvin_obc_end(CS)
147 type(kelvin_obc_cs), pointer :: cs !< Kelvin wave control structure.
148
149 if (associated(cs)) then
150 deallocate(cs)
151 endif
152end subroutine kelvin_obc_end
153
154! -----------------------------------------------------------------------------
155!> This subroutine sets up the Kelvin topography and land mask
156subroutine kelvin_initialize_topography(D, G, param_file, max_depth, US)
157 type(dyn_horgrid_type), intent(in) :: g !< The dynamic horizontal grid type
158 real, dimension(G%isd:G%ied,G%jsd:G%jed), &
159 intent(out) :: d !< Ocean bottom depth [Z ~> m]
160 type(param_file_type), intent(in) :: param_file !< Parameter file structure
161 real, intent(in) :: max_depth !< Maximum model depth [Z ~> m]
162 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
163
164 ! Local variables
165 character(len=40) :: mdl = "Kelvin_initialize_topography" ! This subroutine's name.
166 real :: min_depth ! The minimum and maximum depths [Z ~> m].
167 real :: coast_angle ! Angle of coastline [rad]
168 real :: coast_offset1 ! Longshore distance to coastal angle [L ~> m]
169 real :: coast_offset2 ! Offshore distance to coastal angle [L ~> m]
170 integer :: i, j
171
172 call mom_mesg(" Kelvin_initialization.F90, Kelvin_initialize_topography: setting topography", 5)
173
174 call get_param(param_file, mdl, "MINIMUM_DEPTH", min_depth, &
175 "The minimum depth of the ocean.", units="m", default=0.0, scale=us%m_to_Z)
176 call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_1", coast_offset1, &
177 units="km", default=100.0, do_not_log=.true.)
178 call get_param(param_file, mdl, "ROTATED_COAST_OFFSET_2", coast_offset2, &
179 units="km", default=10.0, do_not_log=.true.)
180 call get_param(param_file, mdl, "ROTATED_COAST_ANGLE", coast_angle, &
181 units="degrees", default=11.3, scale=(atan(1.0)/45.), do_not_log=.true.) ! Convert to radians
182
183 do j=g%jsc,g%jec ; do i=g%isc,g%iec
184 d(i,j) = max_depth
185 ! Southern side
186 if ((g%geoLonT(i,j) - g%west_lon > coast_offset1) .AND. &
187 (atan2(g%geoLatT(i,j) - g%south_lat + coast_offset2, &
188 g%geoLonT(i,j) - g%west_lon - coast_offset1) < coast_angle)) &
189 d(i,j) = 0.5*min_depth
190 ! Northern side
191 if ((g%geoLonT(i,j) - g%west_lon < g%len_lon - coast_offset1) .AND. &
192 (atan2(g%len_lat + g%south_lat + coast_offset2 - g%geoLatT(i,j), &
193 g%len_lon + g%west_lon - coast_offset1 - g%geoLonT(i,j)) < coast_angle)) &
194 d(i,j) = 0.5*min_depth
195
196 if (d(i,j) > max_depth) d(i,j) = max_depth
197 if (d(i,j) < min_depth) d(i,j) = 0.5*min_depth
198 enddo ; enddo
199
200end subroutine kelvin_initialize_topography
201
202!> This subroutine sets the properties of flow at open boundary conditions.
203subroutine kelvin_set_obc_data(OBC, CS, G, GV, US, h, Time)
204 type(ocean_obc_type), pointer :: obc !< This open boundary condition type specifies
205 !! whether, where, and what open boundary
206 !! conditions are used.
207 type(kelvin_obc_cs), pointer :: cs !< Kelvin wave control structure.
208 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
209 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
210 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
211 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< layer thickness [H ~> m or kg m-2].
212 type(time_type), intent(in) :: time !< model time.
213
214 ! The following variables are used to set up the transport in the Kelvin example.
215 real :: time_sec ! The time in the run [T ~> s]
216 real :: cff ! The wave speed [L T-1 ~> m s-1]
217 real :: n0 ! Brunt-Vaisala frequency times a rescaling of slopes [L Z-1 T-1 ~> s-1]
218 real :: lambda ! Offshore decay scale, i.e. the inverse of the deformation radius of a mode [L-1 ~> m-1]
219 real :: omega ! Wave frequency [T-1 ~> s-1]
220 real :: pi ! The ratio of the circumference of a circle to its diameter [nondim]
221 real :: depth_tot(szi_(g),szj_(g)) ! The total depth of the ocean [Z ~> m]
222 real :: depth_tot_vel ! The total depth of the ocean at a velocity point [Z ~> m]
223 real :: depth_tot_corner ! The total depth of the ocean at a vorticity point [Z ~> m]
224 real :: cor_vel ! The Coriolis parameter interpolated to a velocity point [T-1 ~> s-1]
225 real :: mag_ssh ! An overall magnitude of the external wave sea surface height at the coastline [Z ~> m]
226 real :: mag_int ! An overall magnitude of the internal wave at the coastline [L T-1 ~> m s-1]
227 real :: x1, y1 ! Various positions [L ~> m]
228 real :: x, y ! Various positions [L ~> m]
229 real :: sin_wt ! The sine-based periodicity factor [nondim]
230 real :: cos_wt ! The cosine-based periodicity factor [nondim]
231 real :: val2 ! The local wave amplitude [Z ~> m]
232 real :: km_to_l_scale ! A scaling factor from longitudes in km to L [L km-1 ~> 1e3]
233 real :: sina, cosa ! The sine and cosine of the coast angle [nondim]
234 real :: normal_sign ! A variable that corrects the sign of normal velocities for rotation [nondim]
235 real :: trans_sign ! A variable that corrects the sign of transverse velocities for rotation [nondim]
236 type(obc_segment_type), pointer :: segment => null()
237 integer :: unrot_dir ! The unrotated direction of the segment
238 integer :: turns ! Number of index quarter turns
239 integer :: i, j, k, n, is, ie, js, je, isd, ied, jsd, jed, nz
240 integer :: isdb, iedb, jsdb, jedb, isq, ieq, jsq, jeq, is_vel, ie_vel, js_vel, je_vel
241
242 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
243 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
244 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
245
246 if (.not.associated(obc)) call mom_error(fatal, 'Kelvin_initialization.F90: '// &
247 'Kelvin_set_OBC_data() was called but OBC type was not initialized!')
248 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, 'Kelvin_initialization.F90: '// &
249 "Kelvin_set_OBC_data() is only set to work with Cartesian axis units.")
250
251 time_sec = time_to_real(time, scale=us%s_to_T)
252 pi = 4.0*atan(1.0)
253
254 turns = modulo(g%HI%turns, 4)
255
256 if (cs%indexing_bugs .and. (turns /= 0)) call mom_error(fatal, &
257 "Kelvin_set_OBC_data does not support grid rotation when KELVIN_SET_OBC_INDEXING_BUGS is true.")
258
259 do j=jsd,jed ; do i=isd,ied
260 depth_tot(i,j) = 0.0
261 enddo ; enddo
262 do k=1,nz ; do j=jsd,jed ; do i=isd,ied
263 depth_tot(i,j) = depth_tot(i,j) + gv%H_to_Z * h(i,j,k)
264 enddo ; enddo ; enddo
265
266 if (cs%mode == 0) then
267 mag_ssh = cs%ssh_amp
268 omega = 2.0 * pi / cs%wave_period
269 sin_wt = sin(omega * time_sec)
270 else
271 mag_int = cs%inflow_amp
272 n0 = sqrt((cs%rho_range / cs%rho_0) * (gv%g_Earth / cs%H0))
273 lambda = pi * cs%mode * cs%F_0 / (cs%H0 * n0)
274 ! Two wavelengths in domain
275 omega = (4.0 * cs%H0 * n0) / (cs%mode * (g%grid_unit_to_L*g%len_lon))
276 ! If the modal wave speed were calculated via wave_speeds(), we should have
277 ! lambda = CS%F_0 / CS%cg_mode
278 ! omega = (4.0 * PI / (G%grid_unit_to_L*G%len_lon)) * CS%cg_mode
279 endif
280 cos_wt = cos(omega * time_sec)
281
282 sina = sin(cs%coast_angle)
283 cosa = cos(cs%coast_angle)
284 do n=1,obc%number_of_segments
285 segment => obc%segment(n)
286 if (.not. segment%on_pe) cycle
287
288 unrot_dir = segment%direction
289 if (turns /= 0) unrot_dir = rotate_obc_segment_direction(segment%direction, -turns)
290
291 ! Apply values to the inflow end only.
292 if ((unrot_dir == obc_direction_e) .or. (unrot_dir == obc_direction_n)) cycle
293
294 ! Set variables that correct for sign changes during rotation.
295 normal_sign = 1.0
296 if ( (segment%is_E_or_W .and. ((turns == 1) .or. (turns == 2))) .or. &
297 (segment%is_N_or_S .and. ((turns == 2) .or. (turns == 3))) ) normal_sign = -1.0
298
299 ! If OBC_nudging_time is negative, the value of Velocity_nudging_timescale_in that was set
300 ! when the segments are initialized is retained.
301 if (cs%OBC_nudging_time >= 0.0) segment%Velocity_nudging_timescale_in = cs%OBC_nudging_time
302
303 isd = segment%HI%isd ; ied = segment%HI%ied ; isdb = segment%HI%IsdB ; iedb = segment%HI%IedB
304 jsd = segment%HI%jsd ; jed = segment%HI%jed ; jsdb = segment%HI%JsdB ; jedb = segment%HI%JedB
305
306 if (unrot_dir == obc_direction_w) then
307 if (segment%is_E_or_W) then
308 is_vel = isdb ; ie_vel = iedb ; js_vel = jsd ; je_vel = jed
309 else
310 is_vel = isd ; ie_vel = ied ; js_vel = jsdb ; je_vel = jedb
311 endif
312 do j=js_vel,je_vel ; do i=is_vel,ie_vel
313 if (segment%is_E_or_W) then
314 x1 = g%grid_unit_to_L * g%geoLonCu(i,j)
315 y1 = g%grid_unit_to_L * g%geoLatCu(i,j)
316 else
317 x1 = g%grid_unit_to_L * g%geoLonCv(i,j)
318 y1 = g%grid_unit_to_L * g%geoLatCv(i,j)
319 endif
320 x = (x1 - cs%coast_offset1) * cosa + y1 * sina
321 y = -(x1 - cs%coast_offset1) * sina + y1 * cosa
322 if (cs%mode == 0) then
323 ! Use inside bathymetry
324 if (segment%direction == obc_direction_w) then
325 depth_tot_vel = depth_tot(i+1,j)
326 cor_vel = 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1))
327 elseif (segment%direction == obc_direction_s) then
328 depth_tot_vel = depth_tot(i,j+1)
329 cor_vel = 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j))
330 elseif (segment%direction == obc_direction_e) then
331 depth_tot_vel = depth_tot(i,j)
332 cor_vel = 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1))
333 elseif (segment%direction == obc_direction_n) then
334 depth_tot_vel = depth_tot(i,j)
335 cor_vel = 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j))
336 endif
337 cff = sqrt(gv%g_Earth * depth_tot_vel )
338 val2 = mag_ssh * exp(- cor_vel * y / cff)
339 segment%SSH(i,j) = val2 * cos_wt
340 segment%normal_vel_bt(i,j) = (normal_sign*val2) * (sin_wt * cff * cosa / depth_tot_vel )
341 if (segment%nudged) then
342 do k=1,nz
343 segment%nudged_normal_vel(i,j,k) = (normal_sign*val2) * (sin_wt * cff * cosa / depth_tot_vel )
344 enddo
345 elseif (segment%specified) then
346 do k=1,nz
347 segment%normal_vel(i,j,k) = (normal_sign*val2) * (sin_wt * cff * cosa / depth_tot_vel )
348 enddo
349 endif
350 else
351 ! Baroclinic, not rotated yet
352 segment%SSH(i,j) = 0.0
353 segment%normal_vel_bt(i,j) = 0.0
354 ! Use inside bathymetry
355 if (segment%direction == obc_direction_w) then
356 depth_tot_vel = depth_tot(i+1,j)
357 elseif (segment%direction == obc_direction_s) then
358 depth_tot_vel = depth_tot(i,j+1)
359 elseif (segment%direction == obc_direction_e) then
360 depth_tot_vel = depth_tot(i,j)
361 elseif (segment%direction == obc_direction_n) then
362 depth_tot_vel = depth_tot(i,j)
363 endif
364 ! I suspect that the velocities in both of the following loops should instead be
365 ! normal_vel(I,j,k) = CS%inflow_amp * CS%u_struct(k) * exp(-lambda * y) * cos_wt
366 ! In addition, there should be a specification of the interface-height anomalies at the
367 ! open boundaries that are specified as something like
368 ! eta_anom(I,j,K) = (CS%inflow_amp*depth_tot/CS%cg_mode) * CS%w_struct(K) * &
369 ! exp(-lambda * y) * cos_wt
370 ! In these expressions CS%u_struct and CS%w_struct could be returned from the subroutine wave_speeds
371 ! in MOM_wave_speed() based on the horizontally uniform initial state.
372 if (segment%nudged) then
373 do k=1,nz
374 segment%nudged_normal_vel(i,j,k) = (normal_sign*mag_int) * &
375 exp(-lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * cos_wt
376 enddo
377 elseif (segment%specified) then
378 do k=1,nz
379 segment%normal_vel(i,j,k) = (normal_sign*mag_int) * &
380 exp(-lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * cos_wt
381 enddo
382 endif
383 if (associated(segment%h_Reg)) then
384 if (allocated(segment%h_Reg%h)) then
385 do k=1,nz
386 segment%h_Reg%h(i,j,k) = depth_tot_vel / nz + &
387 ((cs%mode * pi) * cs%inflow_amp / (n0 * nz)) * &
388 cos(((pi * k) * cs%mode) / nz) * &
389 exp(-lambda * y) * cos_wt
390 enddo
391 endif
392 endif
393 endif
394 enddo ; enddo
395 endif
396
397 if (unrot_dir == obc_direction_s) then
398 if (segment%is_E_or_W) then
399 is_vel = isdb ; ie_vel = iedb ; js_vel = jsd ; je_vel = jed
400 else
401 is_vel = isd ; ie_vel = ied ; js_vel = jsdb ; je_vel = jedb
402 endif
403 do j=js_vel,je_vel ; do i=is_vel,ie_vel
404 if (segment%is_E_or_W) then
405 x1 = g%grid_unit_to_L * g%geoLonCu(i,j)
406 y1 = g%grid_unit_to_L * g%geoLatCu(i,j)
407 else
408 x1 = g%grid_unit_to_L * g%geoLonCv(i,j)
409 y1 = g%grid_unit_to_L * g%geoLatCv(i,j)
410 endif
411 x = (x1 - cs%coast_offset1) * cosa + y1 * sina
412 y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
413 if (cs%mode == 0) then
414 if (segment%direction == obc_direction_w) then
415 depth_tot_vel = depth_tot(i+1,j)
416 cor_vel = 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1))
417 elseif (segment%direction == obc_direction_s) then
418 depth_tot_vel = depth_tot(i,j+1)
419 cor_vel = 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j))
420 elseif (segment%direction == obc_direction_e) then
421 depth_tot_vel = depth_tot(i,j)
422 cor_vel = 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i,j-1))
423 elseif (segment%direction == obc_direction_n) then
424 depth_tot_vel = depth_tot(i,j)
425 cor_vel = 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j))
426 endif
427 cff = sqrt(gv%g_Earth * depth_tot_vel )
428 val2 = mag_ssh * exp(- cor_vel * y / cff)
429 segment%SSH(i,j) = val2 * cos_wt
430 segment%normal_vel_bt(i,j) = (sin_wt * cff * sina / depth_tot_vel ) * (normal_sign*val2)
431 if (segment%nudged) then
432 do k=1,nz
433 segment%nudged_normal_vel(i,j,k) = (sin_wt * cff * sina / depth_tot_vel) * (normal_sign*val2)
434 enddo
435 elseif (segment%specified) then
436 do k=1,nz
437 segment%normal_vel(i,j,k) = (sin_wt * cff * sina / depth_tot_vel ) * (normal_sign*val2)
438 enddo
439 endif
440 else
441 ! Not rotated yet (also see the notes above on how this case might be improved)
442 segment%SSH(i,j) = 0.0
443 segment%normal_vel_bt(i,j) = 0.0
444 if (segment%nudged) then
445 do k=1,nz
446 segment%nudged_normal_vel(i,j,k) = (normal_sign*mag_int) * &
447 exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * cosa
448 ! This is missing cos_wt
449 enddo
450 elseif (segment%specified) then
451 do k=1,nz
452 segment%normal_vel(i,j,k) = (normal_sign*mag_int) * &
453 exp(- lambda * y) * cos(pi * cs%mode * (k - 0.5) / nz) * cosa
454 ! This is missing cos_wt
455 enddo
456 endif
457 endif
458 enddo ; enddo
459 endif
460
461 if (allocated(segment%tangential_vel)) then
462 trans_sign = 1.0
463 if (segment%is_E_or_W) then
464 isq = isdb ; ieq = iedb ; jsq = jsdb+1 ; jeq = jedb-1
465 if ((turns == 2) .or. (turns == 3)) trans_sign = -1.0
466 else
467 isq = isdb+1 ; ieq = iedb-1 ; jsq = jsdb ; jeq = jedb
468 if ((turns == 1) .or. (turns == 2)) trans_sign = -1.0
469 endif
470
471 if ((unrot_dir == obc_direction_w) .or. (unrot_dir == obc_direction_s)) then
472 do j=jsq,jeq ; do i=isq,ieq
473 if (segment%direction == obc_direction_w) then
474 depth_tot_corner = 0.5*(depth_tot(i+1,j+1) + depth_tot(i+1,j))
475 elseif (segment%direction == obc_direction_e) then
476 depth_tot_corner = 0.5*(depth_tot(i,j+1) + depth_tot(i,j))
477 elseif (segment%direction == obc_direction_s) then
478 depth_tot_corner = 0.5*(depth_tot(i+1,j+1) + depth_tot(i,j+1))
479 elseif (segment%direction == obc_direction_n) then
480 depth_tot_corner = 0.5*(depth_tot(i+1,j) + depth_tot(i,j))
481 endif
482 x1 = g%grid_unit_to_L * g%geoLonBu(i,j)
483 y1 = g%grid_unit_to_L * g%geoLatBu(i,j)
484 x = (x1 - cs%coast_offset1) * cosa + y1 * sina
485 y = - (x1 - cs%coast_offset1) * sina + y1 * cosa
486 cff = sqrt(gv%g_Earth * depth_tot_corner )
487 val2 = (trans_sign*mag_ssh) * exp(- g%CoriolisBu(i,j) * y / cff)
488 if (cs%indexing_bugs) then
489 if (unrot_dir == obc_direction_w) then
490 cff = sqrt(gv%g_Earth * depth_tot(i+1,j) )
491 val2 = (trans_sign*mag_ssh) * exp(- g%CoriolisBu(i,j) * y / cff)
492 endif
493 if (unrot_dir == obc_direction_s) then
494 cff = sqrt(gv%g_Earth * depth_tot(i,j+1) )
495 val2 = (trans_sign*mag_ssh) * exp(- 0.5 * (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j)) * y / cff)
496 endif
497 endif
498 if (cs%mode == 0) then ; do k=1,nz
499 segment%tangential_vel(i,j,k) = (sin_wt * val2 * cff * sina) / depth_tot_corner
500 enddo ; endif
501 enddo ; enddo
502 endif
503 endif
504
505 if (segment%specified .and. (.not.segment%nudged) .and. &
506 ((unrot_dir == obc_direction_s) .or. (unrot_dir == obc_direction_w))) then
507 if (segment%direction == obc_direction_w) then
508 do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
509 segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i+1,j,k) * g%dyCu(i,j)
510 enddo ; enddo ; enddo
511 elseif (segment%direction == obc_direction_e) then
512 do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
513 segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i,j,k) * g%dyCu(i,j)
514 enddo ; enddo ; enddo
515 elseif (segment%direction == obc_direction_s) then
516 do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
517 segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i,j+1,k) * g%dxCv(i,j)
518 enddo ; enddo ; enddo
519 elseif (segment%direction == obc_direction_n) then
520 do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
521 segment%normal_trans(i,j,k) = segment%normal_vel(i,j,k) * h(i,j,k) * g%dxCv(i,j)
522 enddo ; enddo ; enddo
523 endif
524 endif
525
526 enddo
527
528end subroutine kelvin_set_obc_data
529
530end module kelvin_initialization