ROMS
Loading...
Searching...
No Matches
ana_m2obc.h
Go to the documentation of this file.
1!!
2 SUBROUTINE ana_m2obc (ng, tile, model)
3!
4!! git $Id$
5!!======================================================================
6!! Copyright (c) 2002-2025 The ROMS Group !
7!! Licensed under a MIT/X style license !
8!! See License_ROMS.md !
9!=======================================================================
10! !
11! This routine sets 2D momentum open boundary conditions using !
12! analytical expressions. !
13! !
14!=======================================================================
15!
16 USE mod_param
17 USE mod_grid
18 USE mod_ncparam
19 USE mod_ocean
20 USE mod_stepping
21!
22! Imported variable declarations.
23!
24 integer, intent(in) :: ng, tile, model
25!
26! Local variable declarations.
27!
28 character (len=*), parameter :: MyFile = &
29 & __FILE__
30!
31#include "tile.h"
32!
33 CALL ana_m2obc_tile (ng, tile, model, &
34 & lbi, ubi, lbj, ubj, &
35 & imins, imaxs, jmins, jmaxs, &
36 & knew(ng), &
37 & grid(ng) % angler, &
38 & grid(ng) % h, &
39 & grid(ng) % pm, &
40 & grid(ng) % pn, &
41 & grid(ng) % on_u, &
42#ifdef MASKING
43 & grid(ng) % umask, &
44#endif
45 & ocean(ng) % zeta)
46!
47! Set analytical header file name used.
48!
49#ifdef DISTRIBUTE
50 IF (lanafile) THEN
51#else
52 IF (lanafile.and.(tile.eq.0)) THEN
53#endif
54 ananame(12)=myfile
55 END IF
56!
57 RETURN
58 END SUBROUTINE ana_m2obc
59!
60!***********************************************************************
61 SUBROUTINE ana_m2obc_tile (ng, tile, model, &
62 & LBi, UBi, LBj, UBj, &
63 & IminS, ImaxS, JminS, JmaxS, &
64 & knew, &
65 & angler, h, pm, pn, on_u, &
66#ifdef MASKING
67 & umask, &
68#endif
69 & zeta)
70!***********************************************************************
71!
72 USE mod_param
73 USE mod_boundary
74 USE mod_grid
75 USE mod_ncparam
76 USE mod_scalars
77!
78! Imported variable declarations.
79!
80 integer, intent(in) :: ng, tile, model
81 integer, intent(in) :: LBi, UBi, LBj, UBj
82 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
83 integer, intent(in) :: knew
84!
85#ifdef ASSUMED_SHAPE
86 real(r8), intent(in) :: angler(LBi:,LBj:)
87 real(r8), intent(in) :: h(LBi:,LBj:)
88 real(r8), intent(in) :: pm(LBi:,LBj:)
89 real(r8), intent(in) :: pn(LBi:,LBj:)
90 real(r8), intent(in) :: on_u(LBi:,LBj:)
91# ifdef MASKING
92 real(r8), intent(in) :: umask(LBi:,LBj:)
93# endif
94 real(r8), intent(in) :: zeta(LBi:,LBj:,:)
95#else
96 real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
97 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
98 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
99 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
100 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
101# ifdef MASKING
102 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
103# endif
104 real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,:)
105#endif
106!
107! Local variable declarations.
108!
109 integer :: i, j
110!
111 real(r8) :: angle, cff, fac, major, minor, omega, phase, val
112 real(r8) :: ramp
113#if defined ESTUARY_TEST || defined INLET_TEST
114 real(r8) :: my_area, my_flux, tid_flow, riv_flow, cff1, cff2, &
115 & model_flux
116#endif
117#if defined TEST_CHAN
118 real(r8) :: my_area, my_width
119#endif
120
121#include "set_bounds.h"
122!
123!-----------------------------------------------------------------------
124! 2D momentum open boundary conditions.
125!-----------------------------------------------------------------------
126!
127#if defined ESTUARY_TEST
128 IF (lbc(iwest,isubar,ng)%acquire.and. &
129 & lbc(iwest,isvbar,ng)%acquire.and. &
130 & domain(ng)%Western_Edge(tile)) THEN
131 cff1=0.40_r8 ! west end
132 cff2=0.08_r8
133 riv_flow=cff2*300.0_r8*5.0_r8
134 tid_flow=cff1*300.0_r8*10.0_r8
135 my_area=0.0_r8
136 my_flux=0.0_r8
137 DO j=jstrp,jendp
138 cff=0.5_r8*(zeta(istr ,j,knew)+h(istr ,j)+ &
139 & zeta(istr-1,j,knew)+h(istr-1,j))/pn(istr,j)
140 my_area=my_area+cff
141 END DO
142 my_flux=-tid_flow*sin(2.0_r8*pi*time(ng)/ &
143 & (12.0_r8*3600.0_r8))-riv_flow
144 DO j=jstrp,jendp
145 boundary(ng)%ubar_west(j)=my_flux/my_area
146 boundary(ng)%vbar_west(j)=0.0_r8
147 END DO
148 END IF
149
150 IF (lbc(ieast,isubar,ng)%acquire.and. &
151 & lbc(ieast,isvbar,ng)%acquire.and. &
152 & domain(ng)%Eastern_Edge(tile)) THEN
153 cff2=0.08_r8 ! east end
154 riv_flow=cff2*300.0_r8*5.0_r8
155 my_area=0.0_r8
156 my_flux=0.0_r8
157 DO j=jstrp,jendp
158 cff=0.5_r8*(zeta(iend ,j,knew)+h(iend ,j)+ &
159 & zeta(iend+1,j,knew)+h(iend+1,j))/pn(iend,j)
160 my_area=my_area+cff
161 END DO
162 my_flux=-riv_flow
163 DO j=jstrp,jendp
164 boundary(ng)%ubar_east(j)=my_flux/my_area
165 boundary(ng)%vbar_east(j)=0.0_r8
166 END DO
167 END IF
168
169#elif defined KELVIN
170 fac=1.0_r8 ! zeta0
171 omega=2.0_r8*pi/(12.42_r8*3600.0_r8) ! M2 Tide period
172 val=fac*sin(omega*time(ng))
173 IF (lbc(iwest,isubar,ng)%acquire.and. &
174 & lbc(iwest,isvbar,ng)%acquire.and. &
175 & domain(ng)%Western_Edge(tile)) THEN
176 DO j=jstrt,jendt
177 cff=sqrt(g*grid(ng)%h(istr-1,j))
178 boundary(ng)%ubar_west(j)=(val*cff/grid(ng)%h(istr-1,j))* &
179 & exp(-grid(ng)%f(istr-1,j)* &
180 & grid(ng)%yp(istr-1,j)/cff)
181 END DO
182 DO j=jstrp,jendt
183 boundary(ng)%vbar_west(j)=0.0_r8
184 END DO
185 END IF
186
187 IF (lbc(ieast,isubar,ng)%acquire.and. &
188 & lbc(ieast,isvbar,ng)%acquire.and. &
189 & domain(ng)%Eastern_Edge(tile)) THEN
190 DO j=jstrt,jendt
191 cff=sqrt(g*grid(ng)%h(iend,j))
192 val=fac*exp(-grid(ng)%f(iend,j)*grid(ng)%yp(istr-1,j)/cff)
193 boundary(ng)%ubar_east(j)=(val*cff/grid(ng)%h(iend,j))* &
194 & sin(omega*grid(ng)%xp(iend,j)/cff- &
195 & omega*time(ng))
196 END DO
197 DO j=jstrp,jendt
198 boundary(ng)%vbar_east(j)=0.0_r8
199 END DO
200 END IF
201
202#elif defined SED_TEST1
203 IF (lbc(iwest,isubar,ng)%acquire.and. &
204 & lbc(iwest,isvbar,ng)%acquire.and. &
205 & domain(ng)%Western_Edge(tile)) THEN
206 DO j=jstrt,jendt
207 val=0.5_r8*(zeta(istr-1,j,knew)+h(istr-1,j)+ &
208 & zeta(istr ,j,knew)+h(istr ,j))
209 boundary(ng)%ubar_west(j)=-10.0_r8/val
210 END DO
211 DO j=jstrp,jendt
212 boundary(ng)%vbar_west(j)=0.0_r8
213 END DO
214 END IF
215
216 IF (lbc(ieast,isubar,ng)%acquire.and. &
217 & lbc(ieast,isvbar,ng)%acquire.and. &
218 & domain(ng)%Eastern_Edge(tile)) THEN
219 DO j=jstrt,jendt
220 val=0.5_r8*(zeta(iend ,j,knew)+h(iend ,j)+ &
221 & zeta(iend+1,j,knew)+h(iend+1,j))
222 boundary(ng)%ubar_east(j)=-10.0_r8/val
223 END DO
224 DO j=jstrp,jendt
225 boundary(ng)%vbar_east(j)=0.0_r8
226 END DO
227 END IF
228
229#elif defined TEST_CHAN
230 ramp=min(time(ng)/150000.0_r8,1.0_r8)
231 IF (lbc(iwest,isubar,ng)%acquire.and. &
232 & lbc(iwest,isvbar,ng)%acquire.and. &
233 & domain(ng)%Western_Edge(tile)) THEN
234 my_area =0.0_r8
235 my_width=0.0_r8
236 DO j=jstr,jend
237 my_area=my_area+0.5_r8*(zeta(istr-1,j,knew)+h(istr-1,j)+ &
238 & zeta(istr ,j,knew)+h(istr ,j))* &
239 & on_u(istr,j)
240 my_width=my_width+on_u(istr,j)
241 END DO
242 fac=my_width*10.0_r8*1.0_r8*ramp !(width depth ubar)
243 DO j=jstr,jend
244 boundary(ng)%ubar_west(j)=fac/my_area
245 END DO
246 END IF
247
248 IF (lbc(ieast,isubar,ng)%acquire.and. &
249 & lbc(ieast,isvbar,ng)%acquire.and. &
250 & domain(ng)%Eastern_Edge(tile)) THEN
251 my_area =0.0_r8
252 my_width=0.0_r8
253 DO j=jstr,jend
254 my_area=my_area+0.5_r8*(zeta(iend+1,j,knew)+h(iend+1,j)+ &
255 & zeta(iend ,j,knew)+h(iend ,j))* &
256 & on_u(iend,j)
257 my_width=my_width+on_u(iend,j)
258 END DO
259 fac=my_width*10.0_r8*1.0_r8*ramp !(width depth ubar)
260 DO j=jstr,jend
261 boundary(ng)%ubar_east(j)=fac/my_area
262 END DO
263 END IF
264
265#elif defined TRENCH
266 IF (lbc(iwest,isubar,ng)%acquire.and. &
267 & lbc(iwest,isvbar,ng)%acquire.and. &
268 & domain(ng)%Western_Edge(tile)) THEN
269 my_area=0.0_r8
270 my_width=0.0_r8
271 DO j=jstr,jend
272 my_area=my_area+0.5_r8*(zeta(istr-1,j,knew)+h(istr-1,j)+ &
273 & zeta(istr ,j,knew)+h(istr ,j))* &
274 & on_u(istr,j)
275 my_width=my_width+on_u(istr,j)
276 END DO
277 fac=my_width*0.39_r8*0.51_r8 !(width depth ubar)
278 DO j=jstr,jend
279 boundary(ng)%ubar_west(j)=fac/my_area
280 END DO
281 END IF
282
283 IF (lbc(ieast,isubar,ng)%acquire.and. &
284 & lbc(ieast,isvbar,ng)%acquire.and. &
285 & domain(ng)%Eastern_Edge(tile)) THEN
286 my_area=0.0_r8
287 my_width=0.0_r8
288 DO j=jstr,jend
289 my_area=my_area+0.5_r8*(zeta(iend+1,j,knew)+h(iend+1,j)+ &
290 & zeta(iend ,j,knew)+h(iend ,j))* &
291 & on_u(iend,j)
292 my_width=my_width+on_u(iend,j)
293 END DO
294 fac=my_width*0.39_r8*0.51_r8 !(width depth ubar)
295 DO j=jstr,jend
296 boundary(ng)%ubar_east(j)=fac/my_area
297 END DO
298 END IF
299
300#elif defined WEDDELL
301 IF (lbc(iwest,isubar,ng)%acquire.and. &
302 & lbc(iwest,isvbar,ng)%acquire.and. &
303 & domain(ng)%Western_Edge(tile)) THEN
304 fac=tanh((tdays(ng)-dstart)/1.0_r8)
305 omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8) ! M2 Tide period
306 minor=0.0143_r8+(0.0143_r8+0.010_r8)/real(iend+1,r8)
307 major=0.1144_r8+(0.1144_r8-0.013_r8)/real(iend+1,r8)
308 phase=(318.0_r8+(318.0_r8-355.0_r8)/real(iend+1,r8))*deg2rad
309 angle=(125.0_r8+(125.0_r8- 25.0_r8)/real(iend+1,r8))*deg2rad
310 DO j=jstrt,jendt
311 val=0.5_r8*(angler(istr-1,j)+angler(istr,j))
312 boundary(ng)%ubar_west(j)=fac*(major*cos(angle-val)* &
313 & cos(omega-phase)- &
314 & minor*sin(angle-val)* &
315 & sin(omega-phase))
316 END DO
317 DO j=jstrp,jendt
318 val=0.5_r8*(angler(istr-1,j-1)+angler(istr-1,j))
319 boundary(ng)%vbar_west(j)=fac*(major*sin(angle-val)* &
320 & cos(omega-phase)- &
321 & minor*sin(angle-val)* &
322 & cos(omega-phase))
323 END DO
324 END IF
325
326 IF (lbc(ieast,isubar,ng)%acquire.and. &
327 & lbc(ieast,isvbar,ng)%acquire.and. &
328 & domain(ng)%Eastern_Edge(tile)) THEN
329 fac=tanh((tdays(ng)-dstart)/1.0_r8)
330 omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8) ! M2 Tide period
331 minor=0.0143_r8+(0.0143_r8+0.010_r8)
332 major=0.1144_r8+(0.1144_r8-0.013_r8)
333 phase=(318.0_r8+(318.0_r8-355.0_r8))*deg2rad
334 angle=(125.0_r8+(125.0_r8- 25.0_r8))*deg2rad
335 DO j=jstrt,jendt
336 val=0.5_r8*(angler(iend,j)+angler(iend+1,j))
337 boundary(ng)%ubar_east(j)=fac*(major*cos(angle-val)* &
338 & cos(omega-phase)- &
339 & minor*sin(angle-val)* &
340 & sin(omega-phase))
341 END DO
342 DO j=jstrp,jendt
343 val=0.5_r8*(angler(iend+1,j-1)+angler(iend+1,j))
344 boundary(ng)%vbar_east(j)=fac*(major*sin(angle-val)* &
345 & cos(omega-phase)- &
346 & minor*sin(angle-val)* &
347 & cos(omega-phase))
348 END DO
349 END IF
350#else
351 IF (lbc(ieast,isubar,ng)%acquire.and. &
352 & lbc(ieast,isvbar,ng)%acquire.and. &
353 & domain(ng)%Eastern_Edge(tile)) THEN
354 DO j=jstrt,jendt
355 boundary(ng)%ubar_east(j)=0.0_r8
356 END DO
357 DO j=jstrp,jendt
358 boundary(ng)%vbar_east(j)=0.0_r8
359 END DO
360 END IF
361
362 IF (lbc(iwest,isubar,ng)%acquire.and. &
363 & lbc(iwest,isvbar,ng)%acquire.and. &
364 & domain(ng)%Western_Edge(tile)) THEN
365 DO j=jstrt,jendt
366 boundary(ng)%ubar_west(j)=0.0_r8
367 END DO
368 DO j=jstrp,jendt
369 boundary(ng)%vbar_west(j)=0.0_r8
370 END DO
371 END IF
372
373 IF (lbc(isouth,isubar,ng)%acquire.and. &
374 & lbc(isouth,isvbar,ng)%acquire.and. &
375 & domain(ng)%Southern_Edge(tile)) THEN
376 DO i=istrp,iendt
377 boundary(ng)%ubar_south(i)=0.0_r8
378 END DO
379 DO i=istrt,iendt
380 boundary(ng)%vbar_south(i)=0.0_r8
381 END DO
382 END IF
383
384 IF (lbc(inorth,isubar,ng)%acquire.and. &
385 & lbc(inorth,isvbar,ng)%acquire.and. &
386 & domain(ng)%Northern_Edge(tile)) THEN
387 DO i=istrp,iendt
388 boundary(ng)%ubar_north(i)=0.0_r8
389 END DO
390 DO i=istrt,iendt
391 boundary(ng)%vbar_north(i)=0.0_r8
392 END DO
393 END IF
394#endif
395!
396 RETURN
397 END SUBROUTINE ana_m2obc_tile
subroutine ana_m2obc(ng, tile, model)
Definition ana_m2obc.h:3
subroutine ana_m2obc_tile(ng, tile, model, lbi, ubi, lbj, ubj, imins, imaxs, jmins, jmaxs, knew, angler, h, pm, pn, on_u, umask, zeta)
Definition ana_m2obc.h:70
type(t_boundary), dimension(:), allocatable boundary
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
integer isvbar
logical lanafile
integer isubar
character(len=256), dimension(39) ananame
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
type(t_lbc), dimension(:,:,:), allocatable lbc
Definition mod_param.F:375
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
integer, parameter iwest
real(dp), dimension(:), allocatable tdays
real(dp) dstart
real(dp), parameter deg2rad
integer, parameter isouth
integer, parameter ieast
real(dp) g
real(dp), dimension(:), allocatable time
integer, parameter inorth
real(dp), parameter pi
integer, dimension(:), allocatable knew