/[MITgcm]/MITgcm/pkg/diagnostics/diag_calc_psivel.F
ViewVC logotype

Diff of /MITgcm/pkg/diagnostics/diag_calc_psivel.F

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

revision 1.2 by jmc, Wed Jul 6 01:45:32 2011 UTC revision 1.3 by jmc, Mon Aug 1 20:40:14 2011 UTC
# Line 9  C     !INTERFACE: Line 9  C     !INTERFACE:
9        SUBROUTINE DIAG_CALC_PSIVEL(        SUBROUTINE DIAG_CALC_PSIVEL(
10       I                k, iPsi0, jPsi0, uTrans, vTrans,       I                k, iPsi0, jPsi0, uTrans, vTrans,
11       O                psiVel, psiLoc,       O                psiVel, psiLoc,
12       U                myTIme, myIter, myThid )       I                prtMsg, myTime, myIter, myThid )
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
15  C     | SUBROUTINE DIAG_CALC_PSIVEL  C     | SUBROUTINE DIAG_CALC_PSIVEL
# Line 25  C     === Global data === Line 25  C     === Global data ===
25  #include "SIZE.h"  #include "SIZE.h"
26  #include "EEPARAMS.h"  #include "EEPARAMS.h"
27  #include "PARAMS.h"  #include "PARAMS.h"
28  c#include "GRID.h"  #ifdef ALLOW_OBCS
29  c#include "SURFACE.h"  # include "GRID.h"
30    # include "OBCS_GRID.h"
31    #endif /* ALLOW_OBCS */
32    
33  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
34  C     === Routine arguments ===  C     === Routine arguments ===
# Line 36  C     uTrans  :: horizontal transport, u Line 38  C     uTrans  :: horizontal transport, u
38  C     vTrans  :: horizontal transport, u-component  C     vTrans  :: horizontal transport, u-component
39  C     psiVel  :: horizontal stream-function  C     psiVel  :: horizontal stream-function
40  C     psiLoc  :: horizontal stream-function at special location  C     psiLoc  :: horizontal stream-function at special location
41    C     prtMsg  :: do print message to standard-output
42  C     myTime  :: current time of simulation (s)  C     myTime  :: current time of simulation (s)
43  C     myIter  :: current iteration number  C     myIter  :: current iteration number
44  C     myThid  :: my Thread Id number  C     myThid  :: my Thread Id number
# Line 46  C     myThid  :: my Thread Id number Line 49  C     myThid  :: my Thread Id number
49        _RL  vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
50        _RL  psiVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  psiVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
51        _RL  psiLoc(2)        _RL  psiLoc(2)
52          LOGICAL prtMsg
53        _RL  myTime        _RL  myTime
54        INTEGER myIter        INTEGER myIter
55        INTEGER myThid        INTEGER myThid
# Line 65  C     psiOri  :: stream-function value a Line 69  C     psiOri  :: stream-function value a
69        _RL    offSet        _RL    offSet
70        LOGICAL zeroPsi        LOGICAL zeroPsi
71  c     CHARACTER*(MAX_LEN_MBUF) msgBuf  c     CHARACTER*(MAX_LEN_MBUF) msgBuf
72    #ifdef ALLOW_OBCS
73          INTEGER is, js, ix, jx, iy, jy, ijCnt
74          INTEGER npass, nPts, prev_nPts
75          LOGICAL kPsi(1:sNx+1,1:sNy+1)
76    #endif /* ALLOW_OBCS */
77  CEOP  CEOP
78    
79  C--   Initialise  C--   Initialise
# Line 75  C--   Initialise Line 84  C--   Initialise
84  C--   step.1 : compute Psi over each tile separately  C--   step.1 : compute Psi over each tile separately
85        DO bj=myByLo(myThid),myByHi(myThid)        DO bj=myByLo(myThid),myByHi(myThid)
86         DO bi=myBxLo(myThid),myBxHi(myThid)         DO bi=myBxLo(myThid),myBxHi(myThid)
87           dPsiX (bi,bj) = 0.          dPsiX (bi,bj) = 0.
88           dPsiY (bi,bj) = 0.          dPsiY (bi,bj) = 0.
89           psiOri(bi,bj) = 0.          psiOri(bi,bj) = 0.
90            IF ( useOBCS ) THEN
91    C-    Case with OBC:
92    #ifdef ALLOW_OBCS
93    C- Note: OB may introduce discontinuity in domain & tile stream-function map;
94    C     within a tile: define a local "is-set" mask (=kPsi) and
95    C                    propagate stream-function value without assumption.
96    C     between tiles: present code is not "general", likely to work with
97    C                    simple OB setting and/or simple tile connection (no exch2).
98    C     A truly general algorithm requires to change CUMULSUM_Z_TILE (adding 1
99    C     more input dPsi/tile) and to account for disabled tile-connection due
100    C     to OB when setting cumsum tile-mapping (W2_SET_MAP_CUMSUM).
101             DO j=1,sNy+1
102              DO i=1,sNx+1
103                kPsi(i,j) = .FALSE.
104                psiVel(i,j,bi,bj) = 0.
105              ENDDO
106             ENDDO
107    C-    select starting point
108             ijCnt = sNx+sNy+1
109             is = 0
110             js = 0
111             DO j=1,sNy
112              DO i=1,sNx
113               IF ( OBCS_insideMask(i,j,bi,bj).EQ.1.
114         &                      .AND. (i+j).LE.ijCnt ) THEN
115                 is = i
116                 js = j
117                 ijCnt = i+j
118               ENDIF
119              ENDDO
120             ENDDO
121             IF ( is.EQ.0 ) THEN
122               nPts = 0
123             ELSE
124               kPsi(is,js) = .TRUE.
125               nPts = 1
126             ENDIF
127             npass = 0
128             prev_nPts = 0
129             DO WHILE ( nPts.GT.prev_nPts )
130               prev_nPts = nPts
131               npass = npass + 1
132               DO j=1,sNy+1
133                DO i=1,sNx
134                 IF ( OBCS_insideMask(i,j-1,bi,bj).EQ.1. .OR.
135         &            OBCS_insideMask(i, j ,bi,bj).EQ.1. ) THEN
136                   IF ( kPsi(i,j) .AND. .NOT.kPsi(i+1,j) ) THEN
137                     nPts = nPts + 1
138                     kPsi(i+1,j) = .TRUE.
139                     psiVel(i+1,j,bi,bj) = psiVel(i,j,bi,bj)
140         &                               + vTrans(i,j,bi,bj)
141                   ENDIF
142                   IF ( .NOT.kPsi(i,j) .AND. kPsi(i+1,j) ) THEN
143                     nPts = nPts + 1
144                     kPsi(i,j) = .TRUE.
145                     psiVel(i,j,bi,bj) = psiVel(i+1,j,bi,bj)
146         &                               - vTrans(i,j,bi,bj)
147                   ENDIF
148                 ENDIF
149                ENDDO
150               ENDDO
151               DO j=1,sNy
152                DO i=1,sNx+1
153                 IF ( OBCS_insideMask(i-1,j,bi,bj).EQ.1. .OR.
154         &            OBCS_insideMask( i ,j,bi,bj).EQ.1. ) THEN
155                   IF ( kPsi(i,j) .AND. .NOT.kPsi(i,j+1) ) THEN
156                     nPts = nPts + 1
157                     kPsi(i,j+1) = .TRUE.
158                     psiVel(i,j+1,bi,bj) = psiVel(i,j,bi,bj)
159         &                               - uTrans(i,j,bi,bj)
160                   ENDIF
161                   IF ( .NOT.kPsi(i,j) .AND. kPsi(i,j+1) ) THEN
162                     nPts = nPts + 1
163                     kPsi(i,j) = .TRUE.
164                     psiVel(i,j,bi,bj) = psiVel(i,j+1,bi,bj)
165         &                               + uTrans(i,j,bi,bj)
166                   ENDIF
167                 ENDIF
168                ENDDO
169               ENDDO
170               IF ( prtMsg .AND. nPts.GT.prev_nPts ) THEN
171                 _BEGIN_MASTER( myThid )
172                 WRITE(standardMessageUnit,'(A,2I4,A,I6,I8)')
173         &          ' diag_calc_psivel: bi,bj=', bi, bj,
174         &          ' : npass,nPts=', npass, nPts
175                 _END_MASTER( myThid )
176               ENDIF
177    C-    end do while (npass count)
178             ENDDO
179    C-    set tile increments
180             ix = 0
181             jx = 0
182             DO i=sNx+1,1,-1
183              DO j=1,sNy
184               IF ( kPsi(i,j) .AND. jx.EQ.0 ) THEN
185                 ix = i
186                 jx = j
187               ENDIF
188              ENDDO
189             ENDDO
190             IF ( jx.NE.0 ) dPsiX (bi,bj) = psiVel(ix,jx,bi,bj)
191             iy = 0
192             jy = 0
193             DO j=sNy+1,1,-1
194              DO i=1,sNx
195               IF ( kPsi(i,j) .AND. iy.EQ.0 ) THEN
196                 iy = i
197                 jy = j
198               ENDIF
199              ENDDO
200             ENDDO
201             IF ( iy.NE.0 ) dPsiY (bi,bj) = psiVel(iy,jy,bi,bj)
202             IF ( prtMsg ) THEN
203               _BEGIN_MASTER( myThid )
204               WRITE(standardMessageUnit,'(3(A,2I4))')
205         &          ' diag_calc_psivel:            is,js=', is,js,
206         &                 ' ; ix,jx=', ix,jx, ' ; iy,jy=', iy,jy
207               _END_MASTER( myThid )
208             ENDIF
209    #endif /* ALLOW_OBCS */
210            ELSE
211    C-    Case without OBC:
212           psiVel(1,1,bi,bj) = psiOri(bi,bj)           psiVel(1,1,bi,bj) = psiOri(bi,bj)
213           j = 1           j = 1
214           DO i=1,sNx           DO i=1,sNx
# Line 93  C-     note: can vectorise inner loop Line 224  C-     note: can vectorise inner loop
224           ENDDO           ENDDO
225           dPsiX (bi,bj) = psiVel(1+sNx,1,bi,bj)           dPsiX (bi,bj) = psiVel(1+sNx,1,bi,bj)
226           dPsiY (bi,bj) = psiVel(1,1+sNy,bi,bj)           dPsiY (bi,bj) = psiVel(1,1+sNy,bi,bj)
227    C-    end with/without OBC cases
228            ENDIF
229         ENDDO         ENDDO
230        ENDDO        ENDDO
231    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3
  ViewVC Help
Powered by ViewVC 1.1.22