23 implicit none ;
private
25 #include <MOM_memory.h>
38 real :: coast_angle = 0
39 real :: coast_offset1 = 0
40 real :: coast_offset2 = 0
45 logical :: answers_2018
51 #include "version_variable.h"
63 logical :: default_2018_answers
64 character(len=40) :: mdl =
"register_Kelvin_OBC"
65 character(len=32) :: casename =
"Kelvin wave"
66 character(len=200) :: config
68 if (
associated(cs))
then
69 call mom_error(warning,
"register_Kelvin_OBC called with an "// &
70 "associated control structure.")
76 call get_param(param_file, mdl,
"KELVIN_WAVE_MODE", cs%mode, &
77 "Vertical Kelvin wave mode imposed at upstream open boundary.", &
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.)
95 cs%coast_offset1 = cs%coast_offset1 * 1.e3
96 cs%coast_offset2 = cs%coast_offset2 * 1.e3
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.", &
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.)
115 call register_obc(casename, param_file, obc_reg)
124 if (
associated(cs))
then
133 real,
dimension(G%isd:G%ied,G%jsd:G%jed), &
136 real,
intent(in) :: max_depth
140 character(len=40) :: mdl =
"Kelvin_initialize_topography"
144 real :: coast_offset1, coast_offset2, coast_angle, right_angle
147 call mom_mesg(
" Kelvin_initialization.F90, Kelvin_initialize_topography: setting topography", 5)
149 m_to_z = 1.0 ;
if (
present(us)) m_to_z = us%m_to_Z
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.)
160 coast_angle = coast_angle * (atan(1.0)/45.)
161 right_angle = 2 * atan(1.0)
163 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
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
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
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
191 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
192 type(time_type),
intent(in) :: time
195 real :: time_sec, cff
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
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
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!')
215 time_sec = time_type_to_real(time)
219 if (cs%mode == 0)
then
220 omega = 2.0 * pi / (12.42 * 3600.0)
221 val1 = us%m_to_Z * sin(omega * time_sec)
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))
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
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
240 if (segment%direction == obc_direction_e) cycle
241 if (segment%direction == obc_direction_n) cycle
244 segment%Velocity_nudging_timescale_in = 1.0/(0.3*86400)
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)))) )
263 segment%eta(i,j) = 0.0
264 segment%normal_vel_bt(i,j) = 0.0
265 if (segment%nudged)
then
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)
271 elseif (segment%specified)
then
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)
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
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))) ))
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))) )
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
320 segment%eta(i,j) = 0.0
321 segment%normal_vel_bt(i,j) = 0.0
322 if (segment%nudged)
then
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
327 elseif (segment%specified)
then
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)
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
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)))))
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))) ))