/[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.17 - (hide annotations) (download)
Mon Nov 5 19:19:05 2007 UTC (16 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint59p, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59k, checkpoint59j
Changes since 1.16: +3 -2 lines
split PTRACERS.h in 2 header files: PTRACERS_FIELDS.h & PTRACERS_PARAMS.h

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

  ViewVC Help
Powered by ViewVC 1.1.22