/[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.6 by jmc, Fri Jul 22 19:53:40 2011 UTC revision 1.7 by jmc, Mon Aug 1 20:39:47 2011 UTC
# Line 71  C     psiLoc  :: horizontal stream-funct Line 71  C     psiLoc  :: horizontal stream-funct
71        LOGICAL normaliseMatrice, diagNormaliseRHS        LOGICAL normaliseMatrice, diagNormaliseRHS
72        _RL  residCriter, firstResidual, minResidual, 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)
# Line 86  C     psiLoc  :: horizontal stream-funct Line 86  C     psiLoc  :: horizontal stream-funct
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)        DO ks = 1,nlevels(listId)
90          k = NINT(levs(ks,listId))          k = NINT(levs(ks,listId))
91  C--   Solve for velocity potential for each level:  C--   Solve for velocity potential for each level:
92    
# Line 115  C-    calculate cg2d matrix: Line 115  C-    calculate cg2d matrix:
115              ENDDO              ENDDO
116             ENDDO             ENDDO
117    
118  C-    calculate RHS = Div(uVel,vVel):  C-    calculate horizontal transport
119             DO j = 1,sNy+1             DO j = 1,sNy+1
120              DO i = 1,sNx+1              DO i = 1,sNx+1
121                uTrans(i,j,bi,bj) = dyG(i,j,bi,bj)*drF(k)                uTrans(i,j,bi,bj) = dyG(i,j,bi,bj)*drF(k)
# Line 124  C-    calculate RHS = Div(uVel,vVel): Line 124  C-    calculate RHS = Div(uVel,vVel):
124       &                          *qtmp2(i,j,ks,bi,bj)*maskInS(i,j,bi,bj)       &                          *qtmp2(i,j,ks,bi,bj)*maskInS(i,j,bi,bj)
125              ENDDO              ENDDO
126             ENDDO             ENDDO
127    C-   end bi,bj loops
128             ENDDO
129            ENDDO
130    
131    #ifdef ALLOW_OBCS
132            IF ( useOBCS ) THEN
133             CALL OBCS_DIAG_BALANCE(
134         I             k,
135         U             uTrans, vTrans,
136         I             myTime, myIter, myThid )
137            ENDIF
138    #endif /* ALLOW_OBCS */
139    
140    C-    calculate RHS = rAc*Div(uVel,vVel):
141            DO bj = myByLo(myThid), myByHi(myThid)
142             DO bi = myBxLo(myThid), myBxHi(myThid)
143             DO j = 1,sNy             DO j = 1,sNy
144              DO i = 1,sNx              DO i = 1,sNx
145                b2d(i,j,bi,bj) = (                b2d(i,j,bi,bj) = (
# Line 133  C-    calculate RHS = Div(uVel,vVel): Line 149  C-    calculate RHS = Div(uVel,vVel):
149                rhsMax = MAX(ABS(b2d(I,J,bi,bj)),rhsMax)                rhsMax = MAX(ABS(b2d(I,J,bi,bj)),rhsMax)
150              ENDDO              ENDDO
151             ENDDO             ENDDO
   
 C-   end bi,bj loops  
152           ENDDO           ENDDO
153          ENDDO          ENDDO
154    
# Line 143  C-    Normalise Matrice & RHS : Line 157  C-    Normalise Matrice & RHS :
157          normaliseMatrice = .TRUE.          normaliseMatrice = .TRUE.
158          diagNormaliseRHS = .TRUE.          diagNormaliseRHS = .TRUE.
159          a2dNorm = 1. _d 0          a2dNorm = 1. _d 0
160          rhsNorm = 1. _d 0          b2dNorm = 1. _d 0
161            x2dNorm = 1. _d 0
162          IF ( normaliseMatrice ) THEN          IF ( normaliseMatrice ) THEN
163            _GLOBAL_MAX_RL( a2dMax, myThid )            _GLOBAL_MAX_RL( a2dMax, myThid )
164            IF ( a2dMax .GT. 0. _d 0 ) a2dNorm = 1. _d 0/a2dMax            IF ( a2dMax .GT. 0. _d 0 ) a2dNorm = 1. _d 0/a2dMax
165              b2dNorm = a2dNorm
166          ENDIF          ENDIF
167          IF ( diagNormaliseRHS ) THEN          IF ( diagNormaliseRHS ) THEN
168            _GLOBAL_MAX_RL( rhsMax, myThid )            _GLOBAL_MAX_RL( rhsMax, myThid )
169            IF ( rhsMax .GT. 0. _d 0 ) rhsNorm = 1. _d 0/(a2dNorm*rhsMax)            IF ( rhsMax .GT. 0. _d 0 ) THEN
170                b2dNorm = 1. _d 0/rhsMax
171                x2dNorm = a2dNorm*rhsMax
172              ENDIF
173            residCriter = diagCG_resTarget            residCriter = diagCG_resTarget
174          ELSE          ELSE
175            residCriter = a2dNorm * ABS(diagCG_resTarget)            residCriter = a2dNorm * ABS(diagCG_resTarget)
# Line 163  C-    Normalise Matrice & RHS : Line 182  C-    Normalise Matrice & RHS :
182               DO i = 1,sNx+1               DO i = 1,sNx+1
183                aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*a2dNorm                aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*a2dNorm
184                aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*a2dNorm                aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*a2dNorm
185                b2d(i,j,bi,bj)  = b2d(i,j,bi,bj) *a2dNorm*rhsNorm                b2d(i,j,bi,bj)  = b2d(i,j,bi,bj) *b2dNorm
186  c             x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) *rhsNorm  c             x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) /x2dNorm
187               ENDDO               ENDDO
188              ENDDO              ENDDO
189             ENDDO             ENDDO
# Line 174  c             x2d(i,j,bi,bj) =  x2d(i,j, Line 193  c             x2d(i,j,bi,bj) =  x2d(i,j,
193          IF ( debugLevel.GE.debLevA .AND. k.EQ.1 ) THEN          IF ( debugLevel.GE.debLevA .AND. k.EQ.1 ) THEN
194            _BEGIN_MASTER( myThid )            _BEGIN_MASTER( myThid )
195            WRITE(standardMessageUnit,'(A,I9,2(A,1P1E13.6),A,1P1E9.2)')            WRITE(standardMessageUnit,'(A,I9,2(A,1P1E13.6),A,1P1E9.2)')
196       &     ' diag_cg2d (it=', myIter,') a2dNorm,rhsNorm=', a2dNorm,       &     ' diag_cg2d (it=', myIter,') a2dNorm,x2dNorm=', a2dNorm,
197       &     ' ,', rhsNorm, ' ; Criter=', residCriter       &     ' ,', x2dNorm, ' ; Criter=', residCriter
198            _END_MASTER( myThid )            _END_MASTER( myThid )
199          ENDIF          ENDIF
200    
# Line 204  C-    Un-normalise the answer Line 223  C-    Un-normalise the answer
223             DO bi = myBxLo(myThid), myBxHi(myThid)             DO bi = myBxLo(myThid), myBxHi(myThid)
224              DO j = 1-Oly,sNy+Oly              DO j = 1-Oly,sNy+Oly
225               DO i = 1-Olx,sNx+Olx               DO i = 1-Olx,sNx+Olx
226                x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) /rhsNorm                x2d(i,j,bi,bj) =  x2d(i,j,bi,bj)*x2dNorm
227               ENDDO               ENDDO
228              ENDDO              ENDDO
229             ENDDO             ENDDO
# Line 233  C-    Compte divergence-free transport: Line 252  C-    Compte divergence-free transport:
252          CALL DIAG_CALC_PSIVEL(          CALL DIAG_CALC_PSIVEL(
253       I                         k, iPsi0, jPsi0, uTrans, vTrans,       I                         k, iPsi0, jPsi0, uTrans, vTrans,
254       O                         psiVel, psiLoc,       O                         psiVel, psiLoc,
255       I                         myTime, myIter, myThid )       I                         prtFirstCall, myTime, myIter, myThid )
256    
257            _BEGIN_MASTER( myThid)
258          IF ( useCubedSphereExchange .AND.          IF ( useCubedSphereExchange .AND.
259       &       diag_mdsio .AND. myProcId.EQ.0 ) THEN       &       diag_mdsio .AND. myProcId.EQ.0 ) THEN
260  C-      Missing-corner value are not written in MDS output file  C-      Missing-corner value are not written in MDS output file
261  C       Write separately these 2 values (should be part of DIAGNOSTICS_OUT)  C       Write separately these 2 values (should be part of DIAGNOSTICS_OUT)
          _BEGIN_MASTER( myThid)  
262           IF ( diagLoc_ioUnit.EQ.0 ) THEN           IF ( diagLoc_ioUnit.EQ.0 ) THEN
263            CALL MDSFINDUNIT( diagLoc_ioUnit, myThid )            CALL MDSFINDUNIT( diagLoc_ioUnit, myThid )
264            WRITE(dataFName,'(2A,I10.10,A)')            WRITE(dataFName,'(2A,I10.10,A)')
# Line 261  c    &         psiVel(1,1+sNy,nSx,nSy)- Line 280  c    &         psiVel(1,1+sNy,nSx,nSy)-
280  c    &     psiVel(1+sNx,1+sNy,nSx-1,nSy)-psiVel(1,1,1,1),  c    &     psiVel(1+sNx,1+sNy,nSx-1,nSy)-psiVel(1,1,1,1),
281  c    &      ' #', k, lm, myIter  c    &      ' #', k, lm, myIter
282           ENDIF           ENDIF
          _END_MASTER( myThid)  
283          ENDIF          ENDIF
284            IF ( prtFirstCall ) prtFirstCall = .FALSE.
285            _END_MASTER( myThid)
286    
287  C-    Put the results back in qtmp[1,2]  C-    Put the results back in qtmp[1,2]
288          DO bj = myByLo(myThid), myByHi(myThid)          DO bj = myByLo(myThid), myByHi(myThid)
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.7
  ViewVC Help
Powered by ViewVC 1.1.22