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

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

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

revision 1.6 by jmc, Fri Feb 8 22:16:09 2002 UTC revision 1.20 by jmc, Tue Apr 28 18:18:29 2009 UTC
# Line 6  C $Name$ Line 6  C $Name$
6        SUBROUTINE OBCS_CALC( bi, bj, futureTime, futureIter,        SUBROUTINE OBCS_CALC( bi, bj, futureTime, futureIter,
7       &                      uVel, vVel, wVel, theta, salt,       &                      uVel, vVel, wVel, theta, salt,
8       &                      myThid )       &                      myThid )
9  C     /==========================================================\  C     |==========================================================|
10  C     | SUBROUTINE OBCS_CALC                                     |  C     | SUBROUTINE OBCS_CALC                                     |
11  C     | o Calculate future boundary data at open boundaries      |  C     | o Calculate future boundary data at open boundaries      |
12  C     |   at time = futureTime                                   |  C     |   at time = futureTime                                   |
13  C     |==========================================================|  C     |==========================================================|
14  C     |                                                          |  C     |                                                          |
15  C     \==========================================================/  C     |==========================================================|
16        IMPLICIT NONE        IMPLICIT NONE
17    
18  C     === Global variables ===  C     === Global variables ===
19  #include "SIZE.h"  #include "SIZE.h"
20  #include "EEPARAMS.h"  #include "EEPARAMS.h"
21  #include "PARAMS.h"  #include "PARAMS.h"
22    #include "GRID.h"
23  #include "OBCS.h"  #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 ==  C     == Routine arguments ==
32        INTEGER bi, bj        INTEGER bi, bj
# Line 35  C     == Routine arguments == Line 42  C     == Routine arguments ==
42  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
43    
44  C     == Local variables ==  C     == Local variables ==
45    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        INTEGER I, J , K, I_obc, J_obc        INTEGER I, J , K, I_obc, J_obc
49          CHARACTER*(MAX_LEN_MBUF) msgBuf
50    #ifdef ALLOW_OBCS_BALANCE
51          _RL Tr_T, Ar_T, Tr, Ar
52    #endif /* ALLOW_OBCS_BALANCE */
53    #ifdef ALLOW_PTRACERS
54          INTEGER iTracer
55    #endif /* ALLOW_PTRACERS */
56    
57    
58    #ifdef ALLOW_DEBUG
59          IF ( debugLevel .GE. debLevB )
60         &     CALL DEBUG_ENTER('OBCS_CALC',myThid)
61    #endif
62    
63    #ifdef ALLOW_OBCS_EAST
64  C     Eastern OB  C     Eastern OB
65    #ifdef ALLOW_DEBUG
66          IF ( debugLevel .GE. debLevB )
67         &     CALL DEBUG_MSG('OBCS_CALC: East',myThid)
68    #endif
69        IF (useOrlanskiEast) THEN        IF (useOrlanskiEast) THEN
70    #ifdef ALLOW_ORLANSKI
71          CALL ORLANSKI_EAST(          CALL ORLANSKI_EAST(
72       &          bi, bj, futureTime,       &          bi, bj, futureTime,
73       &          uVel, vVel, wVel, theta, salt,       &          uVel, vVel, wVel, theta, salt,
74       &          myThid )       &          myThid )
75    #endif
76        ELSE        ELSE
77          DO K=1,Nr          DO K=1,Nr
78            DO J=1-Oly,sNy+Oly            DO J=1-Oly,sNy+Oly
# Line 62  C     Eastern OB Line 92  C     Eastern OB
92            ENDDO            ENDDO
93          ENDDO          ENDDO
94        ENDIF        ENDIF
95    #endif /* ALLOW_OBCS_EAST */
96    
97    C ------------------------------------------------------------------------------
98    
99    #ifdef ALLOW_OBCS_WEST
100  C     Western OB  C     Western OB
101    #ifdef ALLOW_DEBUG
102          IF ( debugLevel .GE. debLevB )
103         &     CALL DEBUG_MSG('OBCS_CALC: West',myThid)
104    #endif
105        IF (useOrlanskiWest) THEN        IF (useOrlanskiWest) THEN
106    #ifdef ALLOW_ORLANSKI
107          CALL ORLANSKI_WEST(          CALL ORLANSKI_WEST(
108       &          bi, bj, futureTime,       &          bi, bj, futureTime,
109       &          uVel, vVel, wVel, theta, salt,       &          uVel, vVel, wVel, theta, salt,
110       &          myThid )       &          myThid )
111    #endif
112        ELSE        ELSE
113          DO K=1,Nr          DO K=1,Nr
114            DO J=1-Oly,sNy+Oly            DO J=1-Oly,sNy+Oly
# Line 80  C     Western OB Line 120  C     Western OB
120                OBWs(J,K,bi,bj)=sRef(K)                OBWs(J,K,bi,bj)=sRef(K)
121  #ifdef ALLOW_NONHYDROSTATIC  #ifdef ALLOW_NONHYDROSTATIC
122                OBWw(J,K,bi,bj)=0.                OBWw(J,K,bi,bj)=0.
123  #endif  #endif
124  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
125                OBWeta(J,bi,bj)=0.                OBWeta(J,bi,bj)=0.
126  #endif  #endif
127              ENDIF             ENDIF
128            ENDDO            ENDDO
129          ENDDO          ENDDO
130        ENDIF        ENDIF
131    #endif /* ALLOW_OBCS_WEST */
132    
133    C ------------------------------------------------------------------------------
134    
135    #ifdef ALLOW_OBCS_NORTH
136  C         Northern OB  C         Northern OB
137    #ifdef ALLOW_DEBUG
138          IF ( debugLevel .GE. debLevB )
139         &     CALL DEBUG_MSG('OBCS_CALC: North',myThid)
140    #endif
141        IF (useOrlanskiNorth) THEN        IF (useOrlanskiNorth) THEN
142    #ifdef ALLOW_ORLANSKI
143          CALL ORLANSKI_NORTH(          CALL ORLANSKI_NORTH(
144       &          bi, bj, futureTime,       &          bi, bj, futureTime,
145       &          uVel, vVel, wVel, theta, salt,       &          uVel, vVel, wVel, theta, salt,
146       &          myThid )       &          myThid )
147    #endif
148        ELSE        ELSE
149          DO K=1,Nr          DO K=1,Nr
150            DO I=1-Olx,sNx+Olx            DO I=1-Olx,sNx+Olx
# Line 114  C         Northern OB Line 164  C         Northern OB
164            ENDDO            ENDDO
165          ENDDO          ENDDO
166        ENDIF        ENDIF
167    #endif /* ALLOW_OBCS_NORTH */
168    
169    C ------------------------------------------------------------------------------
170    
171    #ifdef ALLOW_OBCS_SOUTH
172  C         Southern OB  C         Southern OB
173    #ifdef ALLOW_DEBUG
174          IF ( debugLevel .GE. debLevB )
175         &     CALL DEBUG_MSG('OBCS_CALC: South',myThid)
176    #endif
177        IF (useOrlanskiSouth) THEN          IF (useOrlanskiSouth) THEN  
178    #ifdef ALLOW_ORLANSKI
179          CALL ORLANSKI_SOUTH(          CALL ORLANSKI_SOUTH(
180       &          bi, bj, futureTime,       &          bi, bj, futureTime,
181       &          uVel, vVel, wVel, theta, salt,       &          uVel, vVel, wVel, theta, salt,
182       &          myThid )       &          myThid )
183    #endif
184        ELSE        ELSE
185          DO K=1,Nr          DO K=1,Nr
186            DO I=1-Olx,sNx+Olx            DO I=1-Olx,sNx+Olx
# Line 140  C         Southern OB Line 200  C         Southern OB
200            ENDDO            ENDDO
201          ENDDO          ENDDO
202        ENDIF        ENDIF
203    #endif /* ALLOW_OBCS_SOUTH */
204    
205    #ifdef ALLOW_PTRACERS
206          IF ( usePTRACERS ) THEN
207    C
208    C     Calculate some default open boundary conditions for passive tracers:
209    C     The default is a homogeneous v.Neumann conditions, that is, the
210    C     tracer gradient across the open boundary is nearly zero;
211    C     only nearly, because the boundary conditions are applied throughout
212    C     the time step during which the interior field does change; therefore
213    C     we have to use the values from the previous time step here. If you
214    C     really want exact v.Neumann conditions, you have to modify
215    C     obcs_apply_ptracer directly.
216    C
217    # ifdef ALLOW_OBCS_EAST
218    C     Eastern OB
219    #  ifdef ALLOW_DEBUG
220           IF ( debugLevel .GE. debLevB )
221         &      CALL DEBUG_MSG('OBCS_CALC: East, pTracers',myThid)
222    #  endif
223           IF (useOrlanskiEast) THEN
224    #  ifdef ALLOW_ORLANSKI
225            WRITE(msgBuf,'(A)')
226         &       'OBCS_CALC: ERROR: useOrlanskiEast Rad OBC with'
227            CALL PRINT_ERROR( msgBuf , 1)
228            WRITE(msgBuf,'(A)')
229         &       'OBCS_CALC: ERROR: pTracers not yet implemented'
230            CALL PRINT_ERROR( msgBuf , 1)
231            STOP 'ABNORMAL END: S/R OBCS_CALC'
232    #  endif
233           ELSE
234            DO iTracer=1,PTRACERS_numInUse
235             DO K=1,Nr
236              DO J=1-Oly,sNy+Oly
237               I_obc=OB_Ie(J,bi,bj)
238               IF (I_obc.ne.0) THEN
239                OBEptr(J,K,bi,bj,iTracer) =
240         &           pTracer(I_obc-1,J,K,bi,bj,iTracer)
241         &           *_maskW(I_obc,J,K,bi,bj)
242               ENDIF
243              ENDDO
244             ENDDO
245            ENDDO
246           ENDIF
247    # endif /* ALLOW_OBCS_EAST */
248    
249    C ------------------------------------------------------------------------------
250    
251    # ifdef ALLOW_OBCS_WEST
252    C     Western OB
253    #  ifdef ALLOW_DEBUG
254           IF ( debugLevel .GE. debLevB )
255         &      CALL DEBUG_MSG('OBCS_CALC: West, pTracers',myThid)
256    #  endif
257           IF (useOrlanskiWest) THEN
258    #  ifdef ALLOW_ORLANSKI
259            WRITE(msgBuf,'(A)')
260         &       'OBCS_CALC: ERROR: useOrlanskiWest Rad OBC with'
261            CALL PRINT_ERROR( msgBuf , 1)
262            WRITE(msgBuf,'(A)')
263         &       'OBCS_CALC: ERROR: pTracers not yet implemented'
264            CALL PRINT_ERROR( msgBuf , 1)
265            STOP 'ABNORMAL END: S/R OBCS_CALC'
266    #  endif
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    #  ifdef ALLOW_ORLANSKI
293            WRITE(msgBuf,'(A)')
294         &       'OBCS_CALC: ERROR: useOrlanskiNorth Rad OBC with'
295            CALL PRINT_ERROR( msgBuf , 1)
296            WRITE(msgBuf,'(A)')
297         &       'OBCS_CALC: ERROR: pTracers not yet implemented'
298            CALL PRINT_ERROR( msgBuf , 1)
299            STOP 'ABNORMAL END: S/R OBCS_CALC'
300    #  endif
301           ELSE
302            DO iTracer=1,PTRACERS_numInUse
303             DO K=1,Nr
304              DO I=1-Olx,sNx+Olx
305               J_obc=OB_Jn(I,bi,bj)
306               IF (J_obc.ne.0) THEN
307                OBNptr(I,K,bi,bj,iTracer) =
308         &           pTracer(I,J_obc-1,K,bi,bj,iTracer)
309         &           *_maskS(I,J_obc,K,bi,bj)
310               ENDIF
311              ENDDO
312             ENDDO
313            ENDDO
314           ENDIF
315    # endif /* ALLOW_OBCS_NORTH */
316    
317    C ------------------------------------------------------------------------------
318    
319    # ifdef ALLOW_OBCS_SOUTH
320    C         Southern OB
321    # ifdef ALLOW_DEBUG
322           IF ( debugLevel .GE. debLevB )
323         &      CALL DEBUG_MSG('OBCS_CALC: South, pTracers',myThid)
324    #endif
325           IF (useOrlanskiSouth) THEN  
326    #ifdef ALLOW_ORLANSKI
327            WRITE(msgBuf,'(A)')
328         &       'OBCS_CALC: ERROR: useOrlanskiSouth Rad OBC with'
329            CALL PRINT_ERROR( msgBuf , 1)
330            WRITE(msgBuf,'(A)')
331         &       'OBCS_CALC: ERROR: pTracers not yet implemented'
332            CALL PRINT_ERROR( msgBuf , 1)
333            STOP 'ABNORMAL END: S/R OBCS_CALC'
334    #endif
335           ELSE
336            DO iTracer=1,PTRACERS_numInUse
337             DO K=1,Nr
338              DO I=1-Olx,sNx+Olx
339               J_obc=OB_Js(I,bi,bj)
340               IF (J_obc.ne.0) THEN
341                OBSptr(I,K,bi,bj,iTracer) =
342         &           pTracer(I,J_obc+1,K,bi,bj,iTracer)
343         &           *_maskS(I,J_obc+1,K,bi,bj)
344               ENDIF
345              ENDDO
346             ENDDO
347            ENDDO
348           ENDIF
349    # endif /* ALLOW_OBCS_SOUTH */
350    C     end if (usePTracers)
351          ENDIF    
352    #endif /* ALLOW_PTRACERS */
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 K=1,Nr
380            DO J=1,sNy
381             I_obc=OB_Ie(J,bi,bj)
382             IF (I_obc.ne.0) THEN
383              Ar = drF(k)*hFacC(I_obc,j,k,bi,bj)*dyG(I_obc,j,bi,bj)
384              Ar_T = Ar_T + Ar
385              Tr_T = Tr_T + Ar * OBEu(J,K,bi,bj)
386             ENDIF
387            ENDDO
388           ENDDO
389           _GLOBAL_SUM_RL( Ar_T , myThid )
390           IF ( Ar_T .GT. 0. _d 0 ) THEN
391            _GLOBAL_SUM_RL( Tr_T , myThid )
392            Tr_T = (0. - Tr_T)/Ar_T
393            DO K=1,Nr
394             DO J=1-Oly,sNy+Oly
395              I_obc=OB_Ie(J,bi,bj)
396              IF (I_obc.ne.0) THEN
397               OBEu(J,K,bi,bj) = OBEu(J,K,bi,bj) + Tr_T
398    c          OBEv(J,K,bi,bj) = 0.
399              ENDIF
400             ENDDO
401            ENDDO
402           ENDIF
403    #endif
404    
405    #ifdef ALLOW_OBCS_WEST
406           Tr_T = 0. _d 0
407           Ar_T = 0. _d 0
408           DO K=1,Nr
409            DO J=1,sNy
410             I_obc=OB_Iw(J,bi,bj)
411             IF (I_obc.ne.0) THEN
412              Ar = drF(k)*hFacC(I_obc,j,k,bi,bj)*dyG(I_obc,j,bi,bj)
413              Ar_T = Ar_T + Ar
414              Tr_T = Tr_T + Ar * OBWu(J,K,bi,bj)
415             ENDIF
416            ENDDO
417           ENDDO
418           _GLOBAL_SUM_RL( Ar_T , myThid )
419           IF ( Ar_T .GT. 0. _d 0 ) THEN
420            _GLOBAL_SUM_RL( Tr_T , myThid )
421            Tr_T = (0. - Tr_T)/Ar_T
422            DO K=1,Nr
423             DO J=1-Oly,sNy+Oly
424              I_obc=OB_Iw(J,bi,bj)
425              IF (I_obc.ne.0) THEN
426               OBWu(J,K,bi,bj) = OBWu(J,K,bi,bj) + Tr_T
427    c          OBWv(J,K,bi,bj) = 0.
428              ENDIF
429             ENDDO
430            ENDDO
431           ENDIF
432    #endif
433    
434    #ifdef ALLOW_OBCS_NORTH
435           Tr_T = 0. _d 0
436           Ar_T = 0. _d 0
437           DO K=1,Nr
438            DO I=1,sNx
439             J_obc=OB_Jn(I,bi,bj)
440             IF (J_obc.ne.0) THEN
441              Ar = drF(k)*hFacC(i,J_obc,k,bi,bj)*dxG(i,J_obc,bi,bj)
442              Ar_T = Ar_T + Ar
443              Tr_T = Tr_T + Ar * OBNv(I,K,bi,bj)
444             ENDIF
445            ENDDO
446           ENDDO
447           _GLOBAL_SUM_RL( Ar_T , myThid )
448           IF ( Ar_T .GT. 0. _d 0 ) THEN
449            _GLOBAL_SUM_RL( Tr_T , myThid )
450            Tr_T = (0. - Tr_T)/Ar_T
451            DO K=1,Nr
452             DO I=1-Olx,sNx+Olx
453              J_obc=OB_Jn(I,bi,bj)
454              IF (J_obc.ne.0) THEN
455    c          OBNu(I,K,bi,bj) = 0.
456               OBNv(I,K,bi,bj) = OBNv(I,K,bi,bj) + Tr_T
457              ENDIF
458             ENDDO
459            ENDDO
460           ENDIF
461    #endif
462    
463    #ifdef ALLOW_OBCS_SOUTH
464           Tr_T = 0. _d 0
465           Ar_T = 0. _d 0
466           DO K=1,Nr
467            DO I=1,sNx
468             J_obc=OB_Js(I,bi,bj)
469             IF (J_obc.ne.0) THEN
470              Ar = drF(k)*hFacC(i,J_obc,k,bi,bj)*dxG(i,J_obc,bi,bj)
471              Ar_T = Ar_T + Ar
472              Tr_T = Tr_T + Ar * OBSv(I,K,bi,bj)
473             ENDIF
474            ENDDO
475           ENDDO
476           _GLOBAL_SUM_RL( Ar_T , myThid )
477           IF ( Ar_T .GT. 0. _d 0 ) THEN
478            _GLOBAL_SUM_RL( Tr_T , myThid )
479            Tr_T = (0. - Tr_T)/Ar_T
480            DO K=1,Nr
481             DO I=1-Olx,sNx+Olx
482              J_obc=OB_Js(I,bi,bj)
483              IF (J_obc.ne.0) THEN
484    c          OBSu(I,K,bi,bj) = 0.
485               OBSv(I,K,bi,bj) = OBSv(I,K,bi,bj) + Tr_T
486              ENDIF
487             ENDDO
488            ENDDO
489           ENDIF
490    #endif
491    
492          ENDIF
493    #endif /* ALLOW_OBCS_BALANCE */
494    
495  #endif /* ALLOW_OBCS */  #endif /* ALLOW_OBCS */
496    
497    #ifdef ALLOW_DEBUG
498          IF ( debugLevel .GE. debLevB )
499         &     CALL DEBUG_LEAVE('OBCS_CALC',myThid)
500    #endif
501        RETURN        RETURN
502        END        END
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.20
  ViewVC Help
Powered by ViewVC 1.1.22