/[MITgcm]/MITgcm/pkg/obcs/obcs_calc.F
ViewVC logotype

Contents of /MITgcm/pkg/obcs/obcs_calc.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.25 - (show annotations) (download)
Thu Dec 17 23:52:11 2009 UTC (14 years, 5 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62a, checkpoint62, checkpoint62b
Changes since 1.24: +8 -16 lines
avoid unused variable

1 C $Header: /u/gcmpack/MITgcm/pkg/obcs/obcs_calc.F,v 1.24 2009/12/15 17:03:28 jahn Exp $
2 C $Name: $
3
4 #include "OBCS_OPTIONS.h"
5
6 SUBROUTINE OBCS_CALC( futureTime, futureIter,
7 & uVel, vVel, wVel, theta, salt,
8 & myThid )
9 C |==========================================================|
10 C | SUBROUTINE OBCS_CALC |
11 C | o Calculate future boundary data at open boundaries |
12 C | at time = futureTime |
13 C |==========================================================|
14 C | |
15 C |==========================================================|
16 IMPLICIT NONE
17
18 C === Global variables ===
19 #include "SIZE.h"
20 #include "EEPARAMS.h"
21 #include "PARAMS.h"
22 #include "GRID.h"
23 #include "OBCS.h"
24 #ifdef ALLOW_PTRACERS
25 #include "PTRACERS_SIZE.h"
26 #include "PTRACERS_PARAMS.h"
27 #include "PTRACERS_FIELDS.h"
28 #include "OBCS_PTRACERS.h"
29 #endif /* ALLOW_PTRACERS */
30
31 C == Routine arguments ==
32 INTEGER futureIter
33 _RL futureTime
34 _RL uVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
35 _RL vVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
36 _RL wVel (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
37 _RL theta(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
38 _RL salt (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
39 INTEGER myThid
40
41 #ifdef ALLOW_OBCS
42
43 C == Local variables ==
44 C bi, bj :: tile indices
45 C I,J,K :: loop indices
46 C I_obc, J_obc :: local index of open boundary
47 C msgBuf :: Informational/error message buffer
48 INTEGER bi, bj
49 INTEGER I, J, K, I_obc, J_obc
50 #ifdef ALLOW_OBCS_BALANCE
51 _RL Tr_T, Ar_T, Tr, Ar
52 #endif /* ALLOW_OBCS_BALANCE */
53 #ifdef ALLOW_PTRACERS
54 CHARACTER*(MAX_LEN_MBUF) msgBuf
55 INTEGER iTracer
56 #endif /* ALLOW_PTRACERS */
57
58
59 #ifdef ALLOW_DEBUG
60 IF ( debugLevel .GE. debLevB )
61 & CALL DEBUG_ENTER('OBCS_CALC',myThid)
62 #endif
63
64 DO bj=myByLo(myThid),myByHi(myThid)
65 DO bi=myBxLo(myThid),myBxHi(myThid)
66
67 #ifdef ALLOW_OBCS_EAST
68 C Eastern OB
69 #ifdef ALLOW_DEBUG
70 IF ( debugLevel .GE. debLevB )
71 & CALL DEBUG_MSG('OBCS_CALC: East',myThid)
72 #endif
73 IF (useOrlanskiEast) THEN
74 #ifdef ALLOW_ORLANSKI
75 CALL ORLANSKI_EAST(
76 & bi, bj, futureTime,
77 & uVel, vVel, wVel, theta, salt,
78 & myThid )
79 #endif
80 ELSE
81 DO K=1,Nr
82 DO J=1-Oly,sNy+Oly
83 I_obc=OB_Ie(J,bi,bj)
84 IF (I_obc.ne.0) THEN
85 OBEu(J,K,bi,bj)=0.
86 OBEv(J,K,bi,bj)=0.
87 OBEt(J,K,bi,bj)=tRef(K)
88 OBEs(J,K,bi,bj)=sRef(K)
89 #ifdef ALLOW_NONHYDROSTATIC
90 OBEw(J,K,bi,bj)=0.
91 #endif
92 #ifdef NONLIN_FRSURF
93 OBEeta(J,bi,bj)=0.
94 #endif
95 ENDIF
96 ENDDO
97 ENDDO
98 ENDIF
99 #endif /* ALLOW_OBCS_EAST */
100
101 C ------------------------------------------------------------------------------
102
103 #ifdef ALLOW_OBCS_WEST
104 C Western OB
105 #ifdef ALLOW_DEBUG
106 IF ( debugLevel .GE. debLevB )
107 & CALL DEBUG_MSG('OBCS_CALC: West',myThid)
108 #endif
109 IF (useOrlanskiWest) THEN
110 #ifdef ALLOW_ORLANSKI
111 CALL ORLANSKI_WEST(
112 & bi, bj, futureTime,
113 & uVel, vVel, wVel, theta, salt,
114 & myThid )
115 #endif
116 ELSE
117 DO K=1,Nr
118 DO J=1-Oly,sNy+Oly
119 I_obc=OB_Iw(J,bi,bj)
120 IF (I_obc.ne.0) THEN
121 OBWu(J,K,bi,bj)=0.
122 OBWv(J,K,bi,bj)=0.
123 OBWt(J,K,bi,bj)=tRef(K)
124 OBWs(J,K,bi,bj)=sRef(K)
125 #ifdef ALLOW_NONHYDROSTATIC
126 OBWw(J,K,bi,bj)=0.
127 #endif
128 #ifdef NONLIN_FRSURF
129 OBWeta(J,bi,bj)=0.
130 #endif
131 ENDIF
132 ENDDO
133 ENDDO
134 ENDIF
135 #endif /* ALLOW_OBCS_WEST */
136
137 C ------------------------------------------------------------------------------
138
139 #ifdef ALLOW_OBCS_NORTH
140 C Northern OB
141 #ifdef ALLOW_DEBUG
142 IF ( debugLevel .GE. debLevB )
143 & CALL DEBUG_MSG('OBCS_CALC: North',myThid)
144 #endif
145 IF (useOrlanskiNorth) THEN
146 #ifdef ALLOW_ORLANSKI
147 CALL ORLANSKI_NORTH(
148 & bi, bj, futureTime,
149 & uVel, vVel, wVel, theta, salt,
150 & myThid )
151 #endif
152 ELSE
153 DO K=1,Nr
154 DO I=1-Olx,sNx+Olx
155 J_obc=OB_Jn(I,bi,bj)
156 IF (J_obc.ne.0) THEN
157 OBNv(I,K,bi,bj)=0.
158 OBNu(I,K,bi,bj)=0.
159 OBNt(I,K,bi,bj)=tRef(K)
160 OBNs(I,K,bi,bj)=sRef(K)
161 #ifdef ALLOW_NONHYDROSTATIC
162 OBNw(I,K,bi,bj)=0.
163 #endif
164 #ifdef NONLIN_FRSURF
165 OBNeta(I,bi,bj)=0.
166 #endif
167 ENDIF
168 ENDDO
169 ENDDO
170 ENDIF
171 #endif /* ALLOW_OBCS_NORTH */
172
173 C ------------------------------------------------------------------------------
174
175 #ifdef ALLOW_OBCS_SOUTH
176 C Southern OB
177 #ifdef ALLOW_DEBUG
178 IF ( debugLevel .GE. debLevB )
179 & CALL DEBUG_MSG('OBCS_CALC: South',myThid)
180 #endif
181 IF (useOrlanskiSouth) THEN
182 #ifdef ALLOW_ORLANSKI
183 CALL ORLANSKI_SOUTH(
184 & bi, bj, futureTime,
185 & uVel, vVel, wVel, theta, salt,
186 & myThid )
187 #endif
188 ELSE
189 DO K=1,Nr
190 DO I=1-Olx,sNx+Olx
191 J_obc=OB_Js(I,bi,bj)
192 IF (J_obc.ne.0) THEN
193 OBSu(I,K,bi,bj)=0.
194 OBSv(I,K,bi,bj)=0.
195 OBSt(I,K,bi,bj)=tRef(K)
196 OBSs(I,K,bi,bj)=sRef(K)
197 #ifdef ALLOW_NONHYDROSTATIC
198 OBSw(I,K,bi,bj)=0.
199 #endif
200 #ifdef NONLIN_FRSURF
201 OBSeta(I,bi,bj)=0.
202 #endif
203 ENDIF
204 ENDDO
205 ENDDO
206 ENDIF
207 #endif /* ALLOW_OBCS_SOUTH */
208
209 #ifdef ALLOW_PTRACERS
210 IF ( usePTRACERS ) THEN
211 C
212 C Calculate some default open boundary conditions for passive tracers:
213 C The default is a homogeneous v.Neumann conditions, that is, the
214 C tracer gradient across the open boundary is nearly zero;
215 C only nearly, because the boundary conditions are applied throughout
216 C the time step during which the interior field does change; therefore
217 C we have to use the values from the previous time step here. If you
218 C really want exact v.Neumann conditions, you have to modify
219 C obcs_apply_ptracer directly.
220 C
221 # ifdef ALLOW_OBCS_EAST
222 C Eastern OB
223 # ifdef ALLOW_DEBUG
224 IF ( debugLevel .GE. debLevB )
225 & CALL DEBUG_MSG('OBCS_CALC: East, pTracers',myThid)
226 # endif
227 IF (useOrlanskiEast) THEN
228 WRITE(msgBuf,'(A)')
229 & 'OBCS_CALC: ERROR: useOrlanskiEast Rad OBC with'
230 CALL PRINT_ERROR( msgBuf, myThid )
231 WRITE(msgBuf,'(A)')
232 & 'OBCS_CALC: ERROR: pTracers not yet implemented'
233 CALL PRINT_ERROR( msgBuf, myThid )
234 STOP 'ABNORMAL END: S/R OBCS_CALC'
235 ELSE
236 DO iTracer=1,PTRACERS_numInUse
237 DO K=1,Nr
238 DO J=1-Oly,sNy+Oly
239 I_obc=OB_Ie(J,bi,bj)
240 IF (I_obc.ne.0) THEN
241 OBEptr(J,K,bi,bj,iTracer) =
242 & pTracer(I_obc-1,J,K,bi,bj,iTracer)
243 & *_maskW(I_obc,J,K,bi,bj)
244 ENDIF
245 ENDDO
246 ENDDO
247 ENDDO
248 ENDIF
249 # endif /* ALLOW_OBCS_EAST */
250
251 C ------------------------------------------------------------------------------
252
253 # ifdef ALLOW_OBCS_WEST
254 C Western OB
255 # ifdef ALLOW_DEBUG
256 IF ( debugLevel .GE. debLevB )
257 & CALL DEBUG_MSG('OBCS_CALC: West, pTracers',myThid)
258 # endif
259 IF (useOrlanskiWest) THEN
260 WRITE(msgBuf,'(A)')
261 & 'OBCS_CALC: ERROR: useOrlanskiWest Rad OBC with'
262 CALL PRINT_ERROR( msgBuf, myThid )
263 WRITE(msgBuf,'(A)')
264 & 'OBCS_CALC: ERROR: pTracers not yet implemented'
265 CALL PRINT_ERROR( msgBuf, myThid )
266 STOP 'ABNORMAL END: S/R OBCS_CALC'
267 ELSE
268 DO iTracer=1,PTRACERS_numInUse
269 DO K=1,Nr
270 DO J=1-Oly,sNy+Oly
271 I_obc=OB_Iw(J,bi,bj)
272 IF (I_obc.ne.0) THEN
273 OBWptr(J,K,bi,bj,iTracer) =
274 & pTracer(I_obc+1,J,K,bi,bj,iTracer)
275 & *_maskW(I_obc+1,J,K,bi,bj)
276 ENDIF
277 ENDDO
278 ENDDO
279 ENDDO
280 ENDIF
281 # endif /* ALLOW_OBCS_WEST */
282
283 C ------------------------------------------------------------------------------
284
285 # ifdef ALLOW_OBCS_NORTH
286 C Northern OB
287 # ifdef ALLOW_DEBUG
288 IF ( debugLevel .GE. debLevB )
289 & CALL DEBUG_MSG('OBCS_CALC: North, pTracers',myThid)
290 # endif
291 IF (useOrlanskiNorth) THEN
292 WRITE(msgBuf,'(A)')
293 & 'OBCS_CALC: ERROR: useOrlanskiNorth Rad OBC with'
294 CALL PRINT_ERROR( msgBuf, myThid )
295 WRITE(msgBuf,'(A)')
296 & 'OBCS_CALC: ERROR: pTracers not yet implemented'
297 CALL PRINT_ERROR( msgBuf, myThid )
298 STOP 'ABNORMAL END: S/R OBCS_CALC'
299 ELSE
300 DO iTracer=1,PTRACERS_numInUse
301 DO K=1,Nr
302 DO I=1-Olx,sNx+Olx
303 J_obc=OB_Jn(I,bi,bj)
304 IF (J_obc.ne.0) THEN
305 OBNptr(I,K,bi,bj,iTracer) =
306 & pTracer(I,J_obc-1,K,bi,bj,iTracer)
307 & *_maskS(I,J_obc,K,bi,bj)
308 ENDIF
309 ENDDO
310 ENDDO
311 ENDDO
312 ENDIF
313 # endif /* ALLOW_OBCS_NORTH */
314
315 C ------------------------------------------------------------------------------
316
317 # ifdef ALLOW_OBCS_SOUTH
318 C Southern OB
319 # ifdef ALLOW_DEBUG
320 IF ( debugLevel .GE. debLevB )
321 & CALL DEBUG_MSG('OBCS_CALC: South, pTracers',myThid)
322 #endif
323 IF (useOrlanskiSouth) THEN
324 WRITE(msgBuf,'(A)')
325 & 'OBCS_CALC: ERROR: useOrlanskiSouth Rad OBC with'
326 CALL PRINT_ERROR( msgBuf, myThid )
327 WRITE(msgBuf,'(A)')
328 & 'OBCS_CALC: ERROR: pTracers not yet implemented'
329 CALL PRINT_ERROR( msgBuf, myThid )
330 STOP 'ABNORMAL END: S/R OBCS_CALC'
331 ELSE
332 DO iTracer=1,PTRACERS_numInUse
333 DO K=1,Nr
334 DO I=1-Olx,sNx+Olx
335 J_obc=OB_Js(I,bi,bj)
336 IF (J_obc.ne.0) THEN
337 OBSptr(I,K,bi,bj,iTracer) =
338 & pTracer(I,J_obc+1,K,bi,bj,iTracer)
339 & *_maskS(I,J_obc+1,K,bi,bj)
340 ENDIF
341 ENDDO
342 ENDDO
343 ENDDO
344 ENDIF
345 # endif /* ALLOW_OBCS_SOUTH */
346 C end if (usePTracers)
347 ENDIF
348 #endif /* ALLOW_PTRACERS */
349
350 C-- end bi,bj loops.
351 ENDDO
352 ENDDO
353
354 C ------------------------------------------------------------------------------
355
356 #ifdef ALLOW_OBCS_PRESCRIBE
357 IF (useOBCSprescribe) THEN
358 C-- Calculate future values on open boundaries
359 #ifdef ALLOW_DEBUG
360 IF ( debugLevel .GE. debLevB )
361 & CALL DEBUG_CALL('OBCS_PRESCRIBE_READ',myThid)
362 #endif
363 CALL OBCS_PRESCRIBE_READ( futureTime, futureIter, myThid )
364 ENDIF
365 #endif /* ALLOW_OBCS_PRESCRIBE */
366
367 C ------------------------------------------------------------------------------
368
369 #ifdef ALLOW_OBCS_BALANCE
370 IF ( useOBCSbalance) THEN
371 #ifdef ALLOW_DEBUG
372 IF ( debugLevel .GE. debLevB )
373 & CALL DEBUG_MSG('useOBCSbalance=.TRUE.',myThid)
374 #endif
375
376 #ifdef ALLOW_OBCS_EAST
377 Tr_T = 0. _d 0
378 Ar_T = 0. _d 0
379 DO bj=myByLo(myThid),myByHi(myThid)
380 DO bi=myBxLo(myThid),myBxHi(myThid)
381 DO K=1,Nr
382 DO J=1,sNy
383 I_obc=OB_Ie(J,bi,bj)
384 IF (I_obc.ne.0) THEN
385 Ar = drF(k)*hFacC(I_obc,j,k,bi,bj)*dyG(I_obc,j,bi,bj)
386 Ar_T = Ar_T + Ar
387 Tr_T = Tr_T + Ar * OBEu(J,K,bi,bj)
388 ENDIF
389 ENDDO
390 ENDDO
391 ENDDO
392 ENDDO
393 _GLOBAL_SUM_RL( Ar_T , myThid )
394 IF ( Ar_T .GT. 0. _d 0 ) THEN
395 _GLOBAL_SUM_RL( Tr_T , myThid )
396 Tr_T = (0. - Tr_T)/Ar_T
397 DO bj=myByLo(myThid),myByHi(myThid)
398 DO bi=myBxLo(myThid),myBxHi(myThid)
399 DO K=1,Nr
400 DO J=1-Oly,sNy+Oly
401 I_obc=OB_Ie(J,bi,bj)
402 IF (I_obc.ne.0) THEN
403 OBEu(J,K,bi,bj) = OBEu(J,K,bi,bj) + Tr_T
404 c OBEv(J,K,bi,bj) = 0.
405 ENDIF
406 ENDDO
407 ENDDO
408 ENDDO
409 ENDDO
410 ENDIF
411 #endif
412
413 #ifdef ALLOW_OBCS_WEST
414 Tr_T = 0. _d 0
415 Ar_T = 0. _d 0
416 DO bj=myByLo(myThid),myByHi(myThid)
417 DO bi=myBxLo(myThid),myBxHi(myThid)
418 DO K=1,Nr
419 DO J=1,sNy
420 I_obc=OB_Iw(J,bi,bj)
421 IF (I_obc.ne.0) THEN
422 Ar = drF(k)*hFacC(I_obc,j,k,bi,bj)*dyG(I_obc,j,bi,bj)
423 Ar_T = Ar_T + Ar
424 Tr_T = Tr_T + Ar * OBWu(J,K,bi,bj)
425 ENDIF
426 ENDDO
427 ENDDO
428 ENDDO
429 ENDDO
430 _GLOBAL_SUM_RL( Ar_T , myThid )
431 IF ( Ar_T .GT. 0. _d 0 ) THEN
432 _GLOBAL_SUM_RL( Tr_T , myThid )
433 Tr_T = (0. - Tr_T)/Ar_T
434 DO bj=myByLo(myThid),myByHi(myThid)
435 DO bi=myBxLo(myThid),myBxHi(myThid)
436 DO K=1,Nr
437 DO J=1-Oly,sNy+Oly
438 I_obc=OB_Iw(J,bi,bj)
439 IF (I_obc.ne.0) THEN
440 OBWu(J,K,bi,bj) = OBWu(J,K,bi,bj) + Tr_T
441 c OBWv(J,K,bi,bj) = 0.
442 ENDIF
443 ENDDO
444 ENDDO
445 ENDDO
446 ENDDO
447 ENDIF
448 #endif
449
450 #ifdef ALLOW_OBCS_NORTH
451 Tr_T = 0. _d 0
452 Ar_T = 0. _d 0
453 DO bj=myByLo(myThid),myByHi(myThid)
454 DO bi=myBxLo(myThid),myBxHi(myThid)
455 DO K=1,Nr
456 DO I=1,sNx
457 J_obc=OB_Jn(I,bi,bj)
458 IF (J_obc.ne.0) THEN
459 Ar = drF(k)*hFacC(i,J_obc,k,bi,bj)*dxG(i,J_obc,bi,bj)
460 Ar_T = Ar_T + Ar
461 Tr_T = Tr_T + Ar * OBNv(I,K,bi,bj)
462 ENDIF
463 ENDDO
464 ENDDO
465 ENDDO
466 ENDDO
467 _GLOBAL_SUM_RL( Ar_T , myThid )
468 IF ( Ar_T .GT. 0. _d 0 ) THEN
469 _GLOBAL_SUM_RL( Tr_T , myThid )
470 Tr_T = (0. - Tr_T)/Ar_T
471 DO bj=myByLo(myThid),myByHi(myThid)
472 DO bi=myBxLo(myThid),myBxHi(myThid)
473 DO K=1,Nr
474 DO I=1-Olx,sNx+Olx
475 J_obc=OB_Jn(I,bi,bj)
476 IF (J_obc.ne.0) THEN
477 c OBNu(I,K,bi,bj) = 0.
478 OBNv(I,K,bi,bj) = OBNv(I,K,bi,bj) + Tr_T
479 ENDIF
480 ENDDO
481 ENDDO
482 ENDDO
483 ENDDO
484 ENDIF
485 #endif
486
487 #ifdef ALLOW_OBCS_SOUTH
488 Tr_T = 0. _d 0
489 Ar_T = 0. _d 0
490 DO bj=myByLo(myThid),myByHi(myThid)
491 DO bi=myBxLo(myThid),myBxHi(myThid)
492 DO K=1,Nr
493 DO I=1,sNx
494 J_obc=OB_Js(I,bi,bj)
495 IF (J_obc.ne.0) THEN
496 Ar = drF(k)*hFacC(i,J_obc,k,bi,bj)*dxG(i,J_obc,bi,bj)
497 Ar_T = Ar_T + Ar
498 Tr_T = Tr_T + Ar * OBSv(I,K,bi,bj)
499 ENDIF
500 ENDDO
501 ENDDO
502 ENDDO
503 ENDDO
504 _GLOBAL_SUM_RL( Ar_T , myThid )
505 IF ( Ar_T .GT. 0. _d 0 ) THEN
506 _GLOBAL_SUM_RL( Tr_T , myThid )
507 Tr_T = (0. - Tr_T)/Ar_T
508 DO bj=myByLo(myThid),myByHi(myThid)
509 DO bi=myBxLo(myThid),myBxHi(myThid)
510 DO K=1,Nr
511 DO I=1-Olx,sNx+Olx
512 J_obc=OB_Js(I,bi,bj)
513 IF (J_obc.ne.0) THEN
514 c OBSu(I,K,bi,bj) = 0.
515 OBSv(I,K,bi,bj) = OBSv(I,K,bi,bj) + Tr_T
516 ENDIF
517 ENDDO
518 ENDDO
519 ENDDO
520 ENDDO
521 ENDIF
522 #endif
523
524 ENDIF
525 #endif /* ALLOW_OBCS_BALANCE */
526
527 #endif /* ALLOW_OBCS */
528
529 #ifdef ALLOW_DEBUG
530 IF ( debugLevel .GE. debLevB )
531 & CALL DEBUG_LEAVE('OBCS_CALC',myThid)
532 #endif
533 RETURN
534 END

  ViewVC Help
Powered by ViewVC 1.1.22