Ticket #211: ana_m2obc.h

File ana_m2obc.h, 12.5 KB (added by m.hadfield, 17 years ago)
Line 
1 SUBROUTINE ana_m2obc (ng, tile, model)
2!
3!! svn $Id: ana_m2obc.h 216 2008-08-08 18:35:19Z arango $
4!!======================================================================
5!! Copyright (c) 2002-2008 The ROMS/TOMS Group !
6!! Licensed under a MIT/X style license !
7!! See License_ROMS.txt !
8!=======================================================================
9! !
10! This routine sets 2D momentum open boundary conditions using !
11! analytical expressions. !
12! !
13!=======================================================================
14!
15 USE mod_param
16 USE mod_grid
17 USE mod_ncparam
18 USE mod_ocean
19 USE mod_stepping
20!
21! Imported variable declarations.
22!
23 integer, intent(in) :: ng, tile, model
24
25#include "tile.h"
26!
27 CALL ana_m2obc_tile (ng, tile, model, &
28 & LBi, UBi, LBj, UBj, &
29 & IminS, ImaxS, JminS, JmaxS, &
30 & knew(ng), &
31 & GRID(ng) % angler, &
32 & GRID(ng) % h, &
33 & GRID(ng) % pm, &
34 & GRID(ng) % pn, &
35 & GRID(ng) % on_u, &
36#ifdef MASKING
37 & GRID(ng) % umask, &
38#endif
39 & OCEAN(ng) % zeta)
40!
41! Set analytical header file name used.
42!
43 IF (Lanafile) THEN
44 ANANAME(12)=__FILE__
45 END IF
46
47 RETURN
48 END SUBROUTINE ana_m2obc
49!
50!***********************************************************************
51 SUBROUTINE ana_m2obc_tile (ng, tile, model, &
52 & LBi, UBi, LBj, UBj, &
53 & IminS, ImaxS, JminS, JmaxS, &
54 & knew, &
55 & angler, h, pm, pn, on_u, &
56#ifdef MASKING
57 & umask, &
58#endif
59 & zeta)
60!***********************************************************************
61!
62 USE mod_param
63 USE mod_boundary
64 USE mod_grid
65 USE mod_scalars
66!
67! Imported variable declarations.
68!
69 integer, intent(in) :: ng, tile, model
70 integer, intent(in) :: LBi, UBi, LBj, UBj
71 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
72 integer, intent(in) :: knew
73!
74#ifdef ASSUMED_SHAPE
75 real(r8), intent(in) :: angler(LBi:,LBj:)
76 real(r8), intent(in) :: h(LBi:,LBj:)
77 real(r8), intent(in) :: pm(LBi:,LBj:)
78 real(r8), intent(in) :: pn(LBi:,LBj:)
79 real(r8), intent(in) :: on_u(LBi:,LBj:)
80# ifdef MASKING
81 real(r8), intent(in) :: umask(LBi:,LBj:)
82# endif
83 real(r8), intent(in) :: zeta(LBi:,LBj:,:)
84#else
85 real(r8), intent(in) :: angler(LBi:UBi,LBj:UBj)
86 real(r8), intent(in) :: h(LBi:UBi,LBj:UBj)
87 real(r8), intent(in) :: pm(LBi:UBi,LBj:UBj)
88 real(r8), intent(in) :: pn(LBi:UBi,LBj:UBj)
89 real(r8), intent(in) :: on_u(LBi:UBi,LBj:UBj)
90# ifdef MASKING
91 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
92# endif
93 real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,3)
94#endif
95!
96! Local variable declarations.
97!
98 integer :: i, j
99 real(r8) :: angle, cff, fac, major, minor, omega, phase, val
100 real(r8) :: ramp
101#if defined ESTUARY_TEST || defined INLET_TEST
102 real(r8) :: my_area, my_flux, tid_flow, riv_flow, cff1, cff2, &
103 & model_flux
104#endif
105#if defined TEST_CHAN
106 real(r8) :: my_area, my_width
107#endif
108
109#include "set_bounds.h"
110!
111!-----------------------------------------------------------------------
112! 2D momentum open boundary conditions.
113!-----------------------------------------------------------------------
114!
115#if defined ESTUARY_TEST
116 cff1=0.40_r8 ! west end
117 cff2=0.08_r8
118 riv_flow=cff2*300.0_r8*5.0_r8
119 tid_flow=cff1*300.0_r8*10.0_r8
120 IF (WESTERN_EDGE) THEN
121 my_area=0.0_r8
122 my_flux=0.0_r8
123 DO j=Jstr,Jend
124 cff=0.5_r8*(zeta(Istr ,j,knew)+h(Istr ,j)+ &
125 & zeta(Istr-1,j,knew)+h(Istr-1,j))/pn(Istr,j)
126 my_area=my_area+cff
127 END DO
128 my_flux=-tid_flow*SIN(2.0_r8*pi*time(ng)/ &
129 & (12.0_r8*3600.0_r8))-riv_flow
130 DO j=Jstr,Jend
131 BOUNDARY(ng)%ubar_west(j)=my_flux/my_area
132 BOUNDARY(ng)%vbar_west(j)=0.0_r8
133 END DO
134 END IF
135 cff2=0.08_r8 ! east end
136 riv_flow=cff2*300.0_r8*5.0_r8
137 IF (EASTERN_EDGE) THEN
138 my_area=0.0_r8
139 my_flux=0.0_r8
140 DO j=Jstr,Jend
141 cff=0.5_r8*(zeta(Iend ,j,knew)+h(Iend ,j)+ &
142 & zeta(Iend+1,j,knew)+h(Iend+1,j))/pn(Iend,j)
143 my_area=my_area+cff
144 END DO
145 my_flux=-riv_flow
146 DO j=Jstr,Jend
147 BOUNDARY(ng)%ubar_east(j)=my_flux/my_area
148 BOUNDARY(ng)%vbar_east(j)=0.0_r8
149 END DO
150 END IF
151#elif defined KELVIN
152 fac=1.0_r8 ! zeta0
153 omega=2.0_r8*pi/(12.42_r8*3600.0_r8) ! M2 Tide period
154 val=fac*SIN(omega*time(ng))
155 IF (WESTERN_EDGE) THEN
156 DO j=JstrR,JendR
157 cff=SQRT(g/GRID(ng)%h(Istr-1,j))
158 BOUNDARY(ng)%ubar_west(j)=val*cff*EXP(-GRID(ng)%f(Istr-1,j)* &
159 & GRID(ng)%yp(Istr-1,j)/ &
160 & cff)
161 END DO
162 DO j=Jstr,JendR
163 BOUNDARY(ng)%vbar_west(j)=0.0_r8
164 END DO
165 END IF
166 IF (EASTERN_EDGE) THEN
167 DO j=JstrR,JendR
168 cff=1.0_r8/SQRT(g*GRID(ng)%h(Iend,j))
169 val=fac*EXP(-GRID(ng)%f(Iend,j)*GRID(ng)%yp(Istr-1,j)*cff)
170 BOUNDARY(ng)%ubar_east(j)=val*SIN(omega*GRID(ng)%xp(Iend,j)* &
171 & cff-omega*time(ng))
172 END DO
173 DO j=Jstr,JendR
174 BOUNDARY(ng)%vbar_east(j)=0.0_r8
175 END DO
176 END IF
177#elif defined SED_TEST1
178 IF (WESTERN_EDGE) THEN
179 DO j=JstrR,JendR
180 val=0.5_r8*(zeta(Istr-1,j,knew)+h(Istr-1,j)+ &
181 & zeta(Istr ,j,knew)+h(Istr ,j))
182 BOUNDARY(ng)%ubar_west(j)=-10.0_r8/val
183 END DO
184 DO j=Jstr,JendR
185 BOUNDARY(ng)%vbar_west(j)=0.0_r8
186 END DO
187 END IF
188 IF (EASTERN_EDGE) THEN
189 DO j=JstrR,JendR
190 val=0.5_r8*(zeta(Iend ,j,knew)+h(Iend ,j)+ &
191 & zeta(Iend+1,j,knew)+h(Iend+1,j))
192 BOUNDARY(ng)%ubar_east(j)=-10.0_r8/val
193 END DO
194 DO j=Jstr,JendR
195 BOUNDARY(ng)%vbar_east(j)=0.0_r8
196 END DO
197 END IF
198#elif defined TEST_CHAN
199 ramp=MIN(time(ng)/150000.0_r8,1.0_r8)
200 IF (WESTERN_EDGE) THEN
201 my_area =0.0_r8
202 my_width=0.0_r8
203 DO j=Jstr,Jend
204 my_area=my_area+0.5_r8*(zeta(Istr-1,j,knew)+h(Istr-1,j)+ &
205 & zeta(Istr ,j,knew)+h(Istr ,j))* &
206 & on_u(Istr,j)
207 my_width=my_width+on_u(Istr,j)
208 END DO
209 fac=my_width*10.0_r8*1.0_r8*ramp !(width depth ubar)
210 DO j=Jstr,Jend
211 BOUNDARY(ng)%ubar_west(j)=fac/my_area
212 END DO
213 END IF
214 IF (EASTERN_EDGE) THEN
215 my_area =0.0_r8
216 my_width=0.0_r8
217 DO j=Jstr,Jend
218 my_area=my_area+0.5_r8*(zeta(Iend+1,j,knew)+h(Iend+1,j)+ &
219 & zeta(Iend ,j,knew)+h(Iend ,j))* &
220 & on_u(Iend,j)
221 my_width=my_width+on_u(Iend,j)
222 END DO
223 fac=my_width*10.0_r8*1.0_r8*ramp !(width depth ubar)
224 DO j=Jstr,Jend
225 BOUNDARY(ng)%ubar_east(j)=fac/my_area
226 END DO
227 END IF
228# elif defined TRENCH
229 IF (WESTERN_EDGE) THEN
230 my_area=0.0_r8
231 my_width=0.0_r8
232 DO j=Jstr,Jend
233 my_area=my_area+0.5_r8*(zeta(Istr-1,j,knew)+h(Istr-1,j)+ &
234 & zeta(Istr ,j,knew)+h(Istr ,j))* &
235 & on_u(Istr,j)
236 my_width=my_width+on_u(Istr,j)
237 END DO
238 fac=my_width*0.39_r8*0.51_r8 !(width depth ubar)
239 DO j=Jstr,Jend
240 BOUNDARY(ng)%ubar_west(j)=fac/my_area
241 END DO
242 END IF
243 IF (EASTERN_EDGE) THEN
244 my_area=0.0_r8
245 my_width=0.0_r8
246 DO j=Jstr,Jend
247 my_area=my_area+0.5_r8*(zeta(Iend+1,j,knew)+h(Iend+1,j)+ &
248 & zeta(Iend ,j,knew)+h(Iend ,j))* &
249 & on_u(Iend,j)
250 my_width=my_width+on_u(Iend,j)
251 END DO
252 fac=my_width*0.39_r8*0.51_r8 !(width depth ubar)
253 DO j=Jstr,Jend
254 BOUNDARY(ng)%ubar_east(j)=fac/my_area
255 END DO
256 END IF
257#elif defined WEDDELL
258 IF (WESTERN_EDGE) THEN
259 fac=TANH((tdays(ng)-dstart)/1.0_r8)
260 omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8) ! M2 Tide period
261 minor=0.0143_r8+(0.0143_r8+0.010_r8)/REAL(Iend+1,r8)
262 major=0.1144_r8+(0.1144_r8-0.013_r8)/REAL(Iend+1,r8)
263 phase=(318.0_r8+(318.0_r8-355.0_r8)/REAL(Iend+1,r8))*deg2rad
264 angle=(125.0_r8+(125.0_r8- 25.0_r8)/REAL(Iend+1,r8))*deg2rad
265 DO j=JstrR,JendR
266 val=0.5_r8*(angler(Istr-1,j)+angler(Istr,j))
267 BOUNDARY(ng)%ubar_west(j)=fac*(major*COS(angle-val)* &
268 & COS(omega-phase)- &
269 & minor*SIN(angle-val)* &
270 & SIN(omega-phase))
271 END DO
272 DO j=Jstr,JendR
273 val=0.5_r8*(angler(Istr-1,j-1)+angler(Istr-1,j))
274 BOUNDARY(ng)%vbar_west(j)=fac*(major*SIN(angle-val)* &
275 & COS(omega-phase)- &
276 & minor*SIN(angle-val)* &
277 & COS(omega-phase))
278 END DO
279 END IF
280 IF (EASTERN_EDGE) THEN
281 fac=TANH((tdays(ng)-dstart)/1.0_r8)
282 omega=2.0_r8*pi*time(ng)/(12.42_r8*3600.0_r8) ! M2 Tide period
283 minor=0.0143_r8+(0.0143_r8+0.010_r8)
284 major=0.1144_r8+(0.1144_r8-0.013_r8)
285 phase=(318.0_r8+(318.0_r8-355.0_r8))*deg2rad
286 angle=(125.0_r8+(125.0_r8- 25.0_r8))*deg2rad
287 DO j=JstrR,JendR
288 val=0.5_r8*(angler(Iend,j)+angler(Iend+1,j))
289 BOUNDARY(ng)%ubar_east(j)=fac*(major*COS(angle-val)* &
290 & COS(omega-phase)- &
291 & minor*SIN(angle-val)* &
292 & SIN(omega-phase))
293 END DO
294 DO j=Jstr,JendR
295 val=0.5_r8*(angler(Iend+1,j-1)+angler(Iend+1,j))
296 BOUNDARY(ng)%vbar_east(j)=fac*(major*SIN(angle-val)* &
297 & COS(omega-phase)- &
298 & minor*SIN(angle-val)* &
299 & COS(omega-phase))
300 END DO
301 END IF
302#else
303 IF (EASTERN_EDGE) THEN
304 DO j=JstrR,JendR
305 BOUNDARY(ng)%ubar_east(j)=0.0_r8
306 END DO
307 DO j=Jstr,JendR
308 BOUNDARY(ng)%vbar_east(j)=0.0_r8
309 END DO
310 END IF
311 IF (WESTERN_EDGE) THEN
312 DO j=JstrR,JendR
313 BOUNDARY(ng)%ubar_west(j)=0.0_r8
314 END DO
315 DO j=Jstr,JendR
316 BOUNDARY(ng)%vbar_west(j)=0.0_r8
317 END DO
318 END IF
319 IF (SOUTHERN_EDGE) THEN
320 DO i=Istr,IendR
321 BOUNDARY(ng)%ubar_south(i)=0.0_r8
322 END DO
323 DO i=IstrR,IendR
324 BOUNDARY(ng)%vbar_south(i)=0.0_r8
325 END DO
326 END IF
327 IF (NORTHERN_EDGE) THEN
328 DO i=Istr,IendR
329 BOUNDARY(ng)%ubar_north(i)=0.0_r8
330 END DO
331 DO i=IstrR,IendR
332 BOUNDARY(ng)%vbar_north(i)=0.0_r8
333 END DO
334 END IF
335#endif
336 RETURN
337 END SUBROUTINE ana_m2obc_tile