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

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

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


Revision 1.14 - (hide annotations) (download)
Wed May 30 09:21:42 2007 UTC (17 years ago) by mlosch
Branch: MAIN
Changes since 1.13: +91 -83 lines
bug fix as reported by Mark Hadfield

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

  ViewVC Help
Powered by ViewVC 1.1.22