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

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

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

revision 1.2 by jmc, Tue Jun 21 18:00:48 2011 UTC revision 1.8 by jmc, Wed Aug 24 02:24:12 2011 UTC
# Line 10  C     !ROUTINE: DIAGNOSTICS_CALC_PHIVEL Line 10  C     !ROUTINE: DIAGNOSTICS_CALC_PHIVEL
10  C     !INTERFACE:  C     !INTERFACE:
11        SUBROUTINE DIAGNOSTICS_CALC_PHIVEL(        SUBROUTINE DIAGNOSTICS_CALC_PHIVEL(
12       I                        listId, md, ndId, ip, im, lm,       I                        listId, md, ndId, ip, im, lm,
13         I                        NrMax,
14       U                        qtmp1, qtmp2,       U                        qtmp1, qtmp2,
15       I                        myTime, myIter, myThid )       I                        myTime, myIter, myThid )
16    
17  C     !DESCRIPTION:  C     !DESCRIPTION:
18  C     Compute Velocity Potential  C     Compute Velocity Potential and Velocity Stream-Function
19    
20  C     !USES:  C     !USES:
21        IMPLICIT NONE        IMPLICIT NONE
# Line 24  C     !USES: Line 25  C     !USES:
25  #include "GRID.h"  #include "GRID.h"
26  #include "DIAGNOSTICS_SIZE.h"  #include "DIAGNOSTICS_SIZE.h"
27  #include "DIAGNOSTICS.h"  #include "DIAGNOSTICS.h"
28    #include "DIAGNOSTICS_CALC.h"
       INTEGER NrMax  
       PARAMETER( NrMax = numLevels )  
   
29    
30  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
31  C     listId  :: Diagnostics list number being written  C     listId  :: Diagnostics list number being written
# Line 36  C     ndId    :: diagnostics  Id number Line 34  C     ndId    :: diagnostics  Id number
34  C     ip      :: diagnostics  pointer to storage array  C     ip      :: diagnostics  pointer to storage array
35  C     im      :: counter-mate pointer to storage array  C     im      :: counter-mate pointer to storage array
36  C     lm      :: index in the averageCycle  C     lm      :: index in the averageCycle
37    C     NrMax   :: 3rd dimension of input/output arrays
38  C     qtmp1   :: horizontal velocity input diag., u-component  C     qtmp1   :: horizontal velocity input diag., u-component
39  C     qtmp2   :: horizontal velocity input diag., v-component  C     qtmp2   :: horizontal velocity input diag., v-component
40  C     myTime  :: current time of simulation (s)  C     myTime  :: current time of simulation (s)
41  C     myIter  :: current iteration number  C     myIter  :: current iteration number
42  C     myThid  :: my Thread Id number  C     myThid  :: my Thread Id number
43        INTEGER listId, md, ndId, ip, im, lm        INTEGER listId, md, ndId, ip, im, lm
44          INTEGER NrMax
45        _RL     qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)        _RL     qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
46        _RL     qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)        _RL     qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
47        _RL     myTime        _RL     myTime
# Line 53  C     qtmp2   :: horizontal-velocity str Line 53  C     qtmp2   :: horizontal-velocity str
53  CEOP  CEOP
54    
55  C     !FUNCTIONS:  C     !FUNCTIONS:
56          INTEGER  ILNBLNK
57          EXTERNAL ILNBLNK
58    
59  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
60  C     i,j,k :: loop indices  C     bi, bj  :: tile indices
61        INTEGER i, j, k  C     i,j,k   :: loop indices
62    C     uTrans  :: horizontal transport, u-component
63    C     vTrans  :: horizontal transport, u-component
64    C     psiVel  :: horizontal stream-function
65    C     psiLoc  :: horizontal stream-function at special location
66        INTEGER bi, bj        INTEGER bi, bj
67        CHARACTER*(MAX_LEN_MBUF) msgBuf        INTEGER i, j, k
68    
69        INTEGER ks        INTEGER ks
70        INTEGER numIters        INTEGER numIters, nIterMin
71        LOGICAL normaliseMatrice, diagNormaliseRHS        LOGICAL normaliseMatrice, diagNormaliseRHS
72        _RL  residCriter, firstResidual, lastResidual        _RL  residCriter, firstResidual, minResidual, lastResidual
73        _RL  a2dMax, a2dNorm        _RL  a2dMax, a2dNorm
74        _RL  rhsMax, rhsNorm        _RL  rhsMax, b2dNorm, x2dNorm
75        _RS  aW2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS  aW2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
76        _RS  aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS  aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
77        _RL  b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
78        _RL  x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
79        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
80        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
81          _RL psiVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
82          _RL psiLoc(2)
83          INTEGER iL
84          CHARACTER*(MAX_LEN_FNAM) dataFName
85          CHARACTER*(MAX_LEN_MBUF) msgBuf
86    
87  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
88    
89        DO ks = 1,kdiag(ndId)  #ifdef ALLOW_DEBUG
90          IF (debugMode) CALL DEBUG_ENTER('DIAGNOSTICS_CALC_PHIVEL',myThid)
91    #endif
92    
93          DO ks = 1,nlevels(listId)
94          k = NINT(levs(ks,listId))          k = NINT(levs(ks,listId))
95  C--   Solve for velocity potential for each level:  C--   Solve for velocity potential for each level:
96    
# Line 95  C-    calculate cg2d matrix: Line 110  C-    calculate cg2d matrix:
110              DO i = 1,sNx+1              DO i = 1,sNx+1
111                aW2d(i,j,bi,bj) = dyG(i,j,bi,bj)*recip_dxC(i,j,bi,bj)                aW2d(i,j,bi,bj) = dyG(i,j,bi,bj)*recip_dxC(i,j,bi,bj)
112       &                         *drF(k)*hFacW(i,j,k,bi,bj)       &                         *drF(k)*hFacW(i,j,k,bi,bj)
113       &                         *maskInW(i,j,bi,bj)       &                         *maskInC(i-1,j,bi,bj)*maskInC(i,j,bi,bj)
114                aS2d(i,j,bi,bj) = dxG(i,j,bi,bj)*recip_dyC(i,j,bi,bj)                aS2d(i,j,bi,bj) = dxG(i,j,bi,bj)*recip_dyC(i,j,bi,bj)
115       &                         *drF(k)*hFacS(i,j,k,bi,bj)       &                         *drF(k)*hFacS(i,j,k,bi,bj)
116       &                         *maskInS(i,j,bi,bj)       &                         *maskInC(i,j-1,bi,bj)*maskInC(i,j,bi,bj)
117                a2dMax = MAX(a2dMax,aW2d(i,j,bi,bj))                a2dMax = MAX(a2dMax,aW2d(i,j,bi,bj))
118                a2dMax = MAX(a2dMax,aS2d(i,j,bi,bj))                a2dMax = MAX(a2dMax,aS2d(i,j,bi,bj))
119              ENDDO              ENDDO
120             ENDDO             ENDDO
121    
122  C-    calculate RHS = Div(uVel,vVel):  C-    calculate horizontal transport
123             DO j = 1,sNy+1             DO j = 1,sNy+1
124              DO i = 1,sNx+1              DO i = 1,sNx+1
125                uTrans(i,j) = dyG(i,j,bi,bj)*drF(k)                uTrans(i,j,bi,bj) = dyG(i,j,bi,bj)*drF(k)
126       &                     *qtmp1(i,j,ks,bi,bj)*maskInW(i,j,bi,bj)       &                          *qtmp1(i,j,ks,bi,bj)*maskInW(i,j,bi,bj)
127                vTrans(i,j) = dxG(i,j,bi,bj)*drF(k)                vTrans(i,j,bi,bj) = dxG(i,j,bi,bj)*drF(k)
128       &                     *qtmp2(i,j,ks,bi,bj)*maskInS(i,j,bi,bj)       &                          *qtmp2(i,j,ks,bi,bj)*maskInS(i,j,bi,bj)
129              ENDDO              ENDDO
130             ENDDO             ENDDO
131    C-   end bi,bj loops
132             ENDDO
133            ENDDO
134    
135    #ifdef ALLOW_OBCS
136            IF ( useOBCS ) THEN
137             CALL OBCS_DIAG_BALANCE(
138         I             k,
139         U             uTrans, vTrans,
140         I             myTime, myIter, myThid )
141            ENDIF
142    #endif /* ALLOW_OBCS */
143    
144    C-    calculate RHS = rAc*Div(uVel,vVel):
145            DO bj = myByLo(myThid), myByHi(myThid)
146             DO bi = myBxLo(myThid), myBxHi(myThid)
147             DO j = 1,sNy             DO j = 1,sNy
148              DO i = 1,sNx              DO i = 1,sNx
149                b2d(i,j,bi,bj) = ( ( uTrans(i+1,j) - uTrans(i,j) )                b2d(i,j,bi,bj) = (
150       &                         + ( vTrans(i,j+1) - vTrans(i,j) )       &                    ( uTrans(i+1,j,bi,bj) - uTrans(i,j,bi,bj) )
151         &                  + ( vTrans(i,j+1,bi,bj) - vTrans(i,j,bi,bj) )
152       &                         )*maskInC(i,j,bi,bj)       &                         )*maskInC(i,j,bi,bj)
153                rhsMax = MAX(ABS(b2d(I,J,bi,bj)),rhsMax)                rhsMax = MAX(ABS(b2d(I,J,bi,bj)),rhsMax)
154              ENDDO              ENDDO
155             ENDDO             ENDDO
   
 C-   end bi,bj loops  
156           ENDDO           ENDDO
157          ENDDO          ENDDO
158    
159  C-    Normalise Matrice & RHS :  C-    Normalise Matrice & RHS :
160          diagNormaliseRHS = cg2dTargetResWunit.LE.0.          diagNormaliseRHS = diagCG_resTarget.GT.0.
161          normaliseMatrice = .TRUE.          normaliseMatrice = .TRUE.
162          diagNormaliseRHS = .TRUE.          diagNormaliseRHS = .TRUE.
163          a2dNorm = 1. _d 0          a2dNorm = 1. _d 0
164          rhsNorm = 1. _d 0          b2dNorm = 1. _d 0
165            x2dNorm = 1. _d 0
166          IF ( normaliseMatrice ) THEN          IF ( normaliseMatrice ) THEN
167            _GLOBAL_MAX_RL( a2dMax, myThid )            _GLOBAL_MAX_RL( a2dMax, myThid )
168            IF ( a2dMax .GE. 0. _d 0 ) a2dNorm = 1. _d 0/a2dMax            IF ( a2dMax .GT. 0. _d 0 ) a2dNorm = 1. _d 0/a2dMax
169              b2dNorm = a2dNorm
170          ENDIF          ENDIF
171          IF ( diagNormaliseRHS ) THEN          IF ( diagNormaliseRHS ) THEN
172            _GLOBAL_MAX_RL( rhsMax, myThid )            _GLOBAL_MAX_RL( rhsMax, myThid )
173            IF ( rhsMax .GE. 0. _d 0 ) rhsNorm = 1. _d 0/(a2dNorm*rhsMax)            IF ( rhsMax .GT. 0. _d 0 ) THEN
174            residCriter = cg2dTargetResidual              b2dNorm = 1. _d 0/rhsMax
175                x2dNorm = a2dNorm*rhsMax
176              ENDIF
177              residCriter = diagCG_resTarget
178          ELSE          ELSE
179            residCriter = a2dNorm * cg2dTargetResWunit            residCriter = a2dNorm * ABS(diagCG_resTarget)
180       &                          * globalArea / deltaTmom       &                          * globalArea / deltaTmom
181          ENDIF          ENDIF
182          IF ( normaliseMatrice .OR. diagNormaliseRHS ) THEN          IF ( normaliseMatrice .OR. diagNormaliseRHS ) THEN
# Line 151  C-    Normalise Matrice & RHS : Line 186  C-    Normalise Matrice & RHS :
186               DO i = 1,sNx+1               DO i = 1,sNx+1
187                aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*a2dNorm                aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*a2dNorm
188                aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*a2dNorm                aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*a2dNorm
189                b2d(i,j,bi,bj)  = b2d(i,j,bi,bj) *a2dNorm*rhsNorm                b2d(i,j,bi,bj)  = b2d(i,j,bi,bj) *b2dNorm
190  c             x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) *rhsNorm  c             x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) /x2dNorm
191               ENDDO               ENDDO
192              ENDDO              ENDDO
193             ENDDO             ENDDO
# Line 162  c             x2d(i,j,bi,bj) =  x2d(i,j, Line 197  c             x2d(i,j,bi,bj) =  x2d(i,j,
197          IF ( debugLevel.GE.debLevA .AND. k.EQ.1 ) THEN          IF ( debugLevel.GE.debLevA .AND. k.EQ.1 ) THEN
198            _BEGIN_MASTER( myThid )            _BEGIN_MASTER( myThid )
199            WRITE(standardMessageUnit,'(A,I9,2(A,1P1E13.6),A,1P1E9.2)')            WRITE(standardMessageUnit,'(A,I9,2(A,1P1E13.6),A,1P1E9.2)')
200       &     ' diag_cg2d (it=', myIter,') a2dNorm,rhsNorm=', a2dNorm,       &     ' diag_cg2d (it=', myIter,') a2dNorm,x2dNorm=', a2dNorm,
201       &     ' ,', rhsNorm, ' ; Criter=', residCriter       &     ' ,', x2dNorm, ' ; Criter=', residCriter
202            _END_MASTER( myThid )            _END_MASTER( myThid )
203          ENDIF          ENDIF
204    
205          numIters = cg2dMaxIters          numIters = diagCG_maxIters
206          CALL DIAG_CG2D(          CALL DIAG_CG2D(
207       I                aW2d, aS2d, b2d,       I                aW2d, aS2d, b2d,
208       I                residCriter,       I                residCriter,
209       O                firstResidual, lastResidual,       O                firstResidual, minResidual, lastResidual,
210       U                x2d, numIters,       U                x2d, numIters,
211       I                myThid )       O                nIterMin,
212         I                diagCG_prtResFrq, myThid )
213    
214          IF ( debugLevel.GE.debLevA ) THEN          IF ( debugLevel.GE.debLevA ) THEN
215            _BEGIN_MASTER( myThid )            _BEGIN_MASTER( myThid )
216            WRITE(standardMessageUnit,'(A,I4,A,I6,A,1P2E17.9)')            WRITE(standardMessageUnit,'(A,I4,A,2I6,A,1P3E14.7)')
217       &     ' diag_cg2d : k=', k, ' , iters=', numIters,       &    ' diag_cg2d : k=', k, ' , it=', nIterMin, numIters,
218       &     ' ; ini,last_Res=', firstResidual, lastResidual       &    ' ; ini,min,last_Res=',firstResidual,minResidual,lastResidual
219            _END_MASTER( myThid )            _END_MASTER( myThid )
220          ENDIF          ENDIF
221    
222  c       _EXCH_XY_RL( x2d, myThid )          _EXCH_XY_RL( x2d, myThid )
223    
224  C-    Un-normalise the answer  C-    Un-normalise the answer
225          IF (diagNormaliseRHS) THEN          IF (diagNormaliseRHS) THEN
226            DO bj = myByLo(myThid), myByHi(myThid)            DO bj = myByLo(myThid), myByHi(myThid)
227             DO bi = myBxLo(myThid), myBxHi(myThid)             DO bi = myBxLo(myThid), myBxHi(myThid)
228              DO j = 1,sNy              DO j = 1-Oly,sNy+Oly
229               DO i = 1,sNx               DO i = 1-Olx,sNx+Olx
230  c             x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) /rhsNorm                x2d(i,j,bi,bj) =  x2d(i,j,bi,bj)*x2dNorm
               qtmp1(i,j,ks,bi,bj) =  x2d(i,j,bi,bj)/rhsNorm  
231               ENDDO               ENDDO
232              ENDDO              ENDDO
233             ENDDO             ENDDO
234            ENDDO            ENDDO
235          ENDIF          ENDIF
236    
237    C-    Compte divergence-free transport:
238            DO bj = myByLo(myThid), myByHi(myThid)
239             DO bi = myBxLo(myThid), myBxHi(myThid)
240              DO j = 1,sNy+1
241               DO i = 1,sNx+1
242                uTrans(i,j,bi,bj) = uTrans(i,j,bi,bj)
243         &                        - ( x2d(i,j,bi,bj) - x2d(i-1,j,bi,bj) )
244         &                         *recip_dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
245         &                         *drF(k)*hFacW(i,j,k,bi,bj)
246         &                         *maskInC(i-1,j,bi,bj)*maskInC(i,j,bi,bj)
247                vTrans(i,j,bi,bj) = vTrans(i,j,bi,bj)
248         &                        - ( x2d(i,j,bi,bj) - x2d(i,j-1,bi,bj) )
249         &                         *recip_dyC(i,j,bi,bj)*dxG(i,j,bi,bj)
250         &                         *drF(k)*hFacS(i,j,k,bi,bj)
251         &                         *maskInC(i,j-1,bi,bj)*maskInC(i,j,bi,bj)
252               ENDDO
253              ENDDO
254             ENDDO
255            ENDDO
256            CALL DIAG_CALC_PSIVEL(
257         I                         k, iPsi0, jPsi0, uTrans, vTrans,
258         O                         psiVel, psiLoc,
259         I                         prtFirstCall, myTime, myIter, myThid )
260    
261            _BEGIN_MASTER( myThid)
262            IF ( useCubedSphereExchange .AND.
263         &       diag_mdsio .AND. myProcId.EQ.0 ) THEN
264    C-      Missing-corner value are not written in MDS output file
265    C       Write separately these 2 values (should be part of DIAGNOSTICS_OUT)
266             IF ( diagLoc_ioUnit.EQ.0 ) THEN
267              CALL MDSFINDUNIT( diagLoc_ioUnit, myThid )
268              WRITE(dataFName,'(2A,I10.10,A)')
269         &         'diags_CScorners', '.', nIter0, '.txt'
270              OPEN( diagLoc_ioUnit, FILE=dataFName, STATUS='unknown' )
271              iL = ILNBLNK(dataFName)
272              WRITE(msgBuf,'(2A,I6,2A)') 'DIAGNOSTICS_CALC_PHIVEL: ',
273         &         'open unit=',diagLoc_ioUnit, ', file: ',dataFName(1:iL)
274              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
275         &                        SQUEEZE_RIGHT, myThid )
276             ENDIF
277             IF ( diagLoc_ioUnit.GT.0 ) THEN
278              WRITE(diagLoc_ioUnit,'(1P2E18.10,A,2I4,I8,A,2I4,I6,2A)')
279         &      psiLoc, ' #', k, lm, myIter,
280         &      ' :',listId, md, ndId, ' ', cdiag(ndId)
281    C-       check accuracy (f1.SW-corner = f6.NW-corner = f5-NE-corner)
282    c         WRITE(0,'(1P2E18.10,A,2I4,I8)')
283    c    &         psiVel(1,1+sNy,nSx,nSy)- psiVel(1,1,1,1),
284    c    &     psiVel(1+sNx,1+sNy,nSx-1,nSy)-psiVel(1,1,1,1),
285    c    &      ' #', k, lm, myIter
286             ENDIF
287            ENDIF
288            IF ( prtFirstCall ) prtFirstCall = .FALSE.
289            _END_MASTER( myThid)
290    
291    C-    Put the results back in qtmp[1,2]
292            DO bj = myByLo(myThid), myByHi(myThid)
293             DO bi = myBxLo(myThid), myBxHi(myThid)
294              DO j = 1,sNy+1
295               DO i = 1,sNx+1
296                  qtmp1(i,j,ks,bi,bj) =  x2d(i,j,bi,bj)
297                  qtmp2(i,j,ks,bi,bj) =  psiVel(i,j,bi,bj)
298               ENDDO
299              ENDDO
300             ENDDO
301            ENDDO
302    
303        ENDDO        ENDDO
304    
305    #ifdef ALLOW_DEBUG
306          IF (debugMode) CALL DEBUG_LEAVE('DIAGNOSTICS_CALC_PHIVEL',myThid)
307    #endif
308    
309        RETURN        RETURN
310        END        END
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22