/[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.3 by jmc, Fri Jun 24 22:15:48 2011 UTC revision 1.4 by jmc, Fri Jul 1 18:50:00 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    c#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    C     i,jPsi0 :: indices of grid-point location where Psi == 0
67        INTEGER bi, bj        INTEGER bi, bj
68  c     CHARACTER*(MAX_LEN_MBUF) msgBuf        INTEGER i, j, k
69    
70        INTEGER ks        INTEGER ks
71        INTEGER numIters        INTEGER numIters
# Line 70  c     CHARACTER*(MAX_LEN_MBUF) msgBuf Line 77  c     CHARACTER*(MAX_LEN_MBUF) msgBuf
77        _RS  aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RS  aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
78        _RL  b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
79        _RL  x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        _RL  x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
80        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
81        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
82          _RL psiVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
83          _RL psiLoc(2)
84          INTEGER iPsi0(nSx,nSy)
85          INTEGER jPsi0(nSx,nSy)
86          INTEGER iL
87          CHARACTER*(MAX_LEN_FNAM) dataFName
88          CHARACTER*(MAX_LEN_MBUF) msgBuf
89    
90  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
91    
92          DO bj=myByLo(myThid),myByHi(myThid)
93           DO bi=myBxLo(myThid),myBxHi(myThid)
94             iPsi0(bi,bj) = -1
95             jPsi0(bi,bj) =  0
96           ENDDO
97          ENDDO
98    
99        DO ks = 1,kdiag(ndId)        DO ks = 1,kdiag(ndId)
100          k = NINT(levs(ks,listId))          k = NINT(levs(ks,listId))
101  C--   Solve for velocity potential for each level:  C--   Solve for velocity potential for each level:
# Line 107  C-    calculate cg2d matrix: Line 128  C-    calculate cg2d matrix:
128  C-    calculate RHS = Div(uVel,vVel):  C-    calculate RHS = Div(uVel,vVel):
129             DO j = 1,sNy+1             DO j = 1,sNy+1
130              DO i = 1,sNx+1              DO i = 1,sNx+1
131                uTrans(i,j) = dyG(i,j,bi,bj)*drF(k)                uTrans(i,j,bi,bj) = dyG(i,j,bi,bj)*drF(k)
132       &                     *qtmp1(i,j,ks,bi,bj)*maskInW(i,j,bi,bj)       &                          *qtmp1(i,j,ks,bi,bj)*maskInW(i,j,bi,bj)
133                vTrans(i,j) = dxG(i,j,bi,bj)*drF(k)                vTrans(i,j,bi,bj) = dxG(i,j,bi,bj)*drF(k)
134       &                     *qtmp2(i,j,ks,bi,bj)*maskInS(i,j,bi,bj)       &                          *qtmp2(i,j,ks,bi,bj)*maskInS(i,j,bi,bj)
135              ENDDO              ENDDO
136             ENDDO             ENDDO
137             DO j = 1,sNy             DO j = 1,sNy
138              DO i = 1,sNx              DO i = 1,sNx
139                b2d(i,j,bi,bj) = ( ( uTrans(i+1,j) - uTrans(i,j) )                b2d(i,j,bi,bj) = (
140       &                         + ( vTrans(i,j+1) - vTrans(i,j) )       &                    ( uTrans(i+1,j,bi,bj) - uTrans(i,j,bi,bj) )
141         &                  + ( vTrans(i,j+1,bi,bj) - vTrans(i,j,bi,bj) )
142       &                         )*maskInC(i,j,bi,bj)       &                         )*maskInC(i,j,bi,bj)
143                rhsMax = MAX(ABS(b2d(I,J,bi,bj)),rhsMax)                rhsMax = MAX(ABS(b2d(I,J,bi,bj)),rhsMax)
144              ENDDO              ENDDO
# Line 183  c             x2d(i,j,bi,bj) =  x2d(i,j, Line 205  c             x2d(i,j,bi,bj) =  x2d(i,j,
205            _END_MASTER( myThid )            _END_MASTER( myThid )
206          ENDIF          ENDIF
207    
208  c       _EXCH_XY_RL( x2d, myThid )          _EXCH_XY_RL( x2d, myThid )
209    
210  C-    Un-normalise the answer  C-    Un-normalise the answer
211          IF (diagNormaliseRHS) THEN          IF (diagNormaliseRHS) THEN
212            DO bj = myByLo(myThid), myByHi(myThid)            DO bj = myByLo(myThid), myByHi(myThid)
213             DO bi = myBxLo(myThid), myBxHi(myThid)             DO bi = myBxLo(myThid), myBxHi(myThid)
214              DO j = 1,sNy              DO j = 1-Oly,sNy+Oly
215               DO i = 1,sNx               DO i = 1-Olx,sNx+Olx
216  c             x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) /rhsNorm                x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) /rhsNorm
               qtmp1(i,j,ks,bi,bj) =  x2d(i,j,bi,bj)/rhsNorm  
217               ENDDO               ENDDO
218              ENDDO              ENDDO
219             ENDDO             ENDDO
220            ENDDO            ENDDO
221          ENDIF          ENDIF
222    
223    C-    Compte divergence-free transport:
224            DO bj = myByLo(myThid), myByHi(myThid)
225             DO bi = myBxLo(myThid), myBxHi(myThid)
226              DO j = 1,sNy+1
227               DO i = 1,sNx+1
228                uTrans(i,j,bi,bj) = uTrans(i,j,bi,bj)
229         &                        - ( x2d(i,j,bi,bj) - x2d(i-1,j,bi,bj) )
230         &                         *recip_dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
231         &                         *drF(k)*hFacW(i,j,k,bi,bj)
232         &                         *maskInW(i,j,bi,bj)
233                vTrans(i,j,bi,bj) = vTrans(i,j,bi,bj)
234         &                        - ( x2d(i,j,bi,bj) - x2d(i,j-1,bi,bj) )
235         &                         *recip_dyC(i,j,bi,bj)*dxG(i,j,bi,bj)
236         &                         *drF(k)*hFacS(i,j,k,bi,bj)
237         &                         *maskInS(i,j,bi,bj)
238               ENDDO
239              ENDDO
240             ENDDO
241            ENDDO
242            CALL DIAG_CALC_PSIVEL(
243         I                         k, iPsi0, jPsi0, uTrans, vTrans,
244         O                         psiVel, psiLoc,
245         I                         myTime, myIter, myThid )
246    
247            IF ( useCubedSphereExchange .AND.
248         &       diag_mdsio .AND. myProcId.EQ.0 ) THEN
249    C-      Missing-corner value are not written in MDS output file
250    C       Write separately these 2 values (should be part of DIAGNOSTICS_OUT)
251             _BEGIN_MASTER( myThid)
252             IF ( diagLoc_ioUnit.EQ.0 ) THEN
253              CALL MDSFINDUNIT( diagLoc_ioUnit, myThid )
254              WRITE(dataFName,'(2A,I10.10,A)')
255         &         'diags_CScorners', '.', nIter0, '.txt'
256              OPEN( diagLoc_ioUnit, FILE=dataFName, STATUS='unknown' )
257              iL = ILNBLNK(dataFName)
258              WRITE(msgBuf,'(2A,I6,2A)') 'DIAGNOSTICS_CALC_PHIVEL: ',
259         &         'open unit=',diagLoc_ioUnit, ', file: ',dataFName(1:iL)
260              CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
261         &                        SQUEEZE_RIGHT, myThid )
262             ENDIF
263             IF ( diagLoc_ioUnit.GT.0 ) THEN
264              WRITE(diagLoc_ioUnit,'(1P2E18.10,A,2I4,I8,A,2I4,I6,2A)')
265         &      psiLoc, ' #', k, lm, myIter,
266         &      ' :',listId, md, ndId, ' ', cdiag(ndId)
267    C-       check accuracy (f1.SW-corner = f6.NW-corner = f5-NE-corner)
268    c         WRITE(0,'(1P2E18.10,A,2I4,I8)')
269    c    &         psiVel(1,1+sNy,nSx,nSy)- psiVel(1,1,1,1),
270    c    &     psiVel(1+sNx,1+sNy,nSx-1,nSy)-psiVel(1,1,1,1),
271    c    &      ' #', k, lm, myIter
272             ENDIF
273             _END_MASTER( myThid)
274            ENDIF
275    
276    C-    Put the results back in qtmp[1,2]
277            DO bj = myByLo(myThid), myByHi(myThid)
278             DO bi = myBxLo(myThid), myBxHi(myThid)
279              DO j = 1,sNy+1
280               DO i = 1,sNx+1
281                  qtmp1(i,j,ks,bi,bj) =  x2d(i,j,bi,bj)
282                  qtmp2(i,j,ks,bi,bj) =  psiVel(i,j,bi,bj)
283               ENDDO
284              ENDDO
285             ENDDO
286            ENDDO
287    
288        ENDDO        ENDDO
289    
290        RETURN        RETURN
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.22