pkg/diagnostics/diagnostics_calc_phivel.F

d5ada4102b Jean*0001 #include "DIAG_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP 0
                0005 C     !ROUTINE: DIAGNOSTICS_CALC_PHIVEL
                0006
                0007 C     !INTERFACE:
                0008       SUBROUTINE DIAGNOSTICS_CALC_PHIVEL(
                0009      I                        listId, md, ndId, ip, im, lm,
bc85e7ab92 Jean*0010      I                        NrMax,
48a533dac6 Jean*0011      U                        qtmp1, qtmp2,
d5ada4102b Jean*0012      I                        myTime, myIter, myThid )
                0013
                0014 C     !DESCRIPTION:
bc85e7ab92 Jean*0015 C     Compute Velocity Potential and Velocity Stream-Function
d5ada4102b Jean*0016
                0017 C     !USES:
                0018       IMPLICIT NONE
                0019 #include "SIZE.h"
                0020 #include "EEPARAMS.h"
                0021 #include "PARAMS.h"
                0022 #include "GRID.h"
                0023 #include "DIAGNOSTICS_SIZE.h"
                0024 #include "DIAGNOSTICS.h"
49d8041f88 Jean*0025 #include "DIAGNOSTICS_CALC.h"
d5ada4102b Jean*0026
                0027 C     !INPUT PARAMETERS:
                0028 C     listId  :: Diagnostics list number being written
                0029 C     md      :: field number in the list "listId".
                0030 C     ndId    :: diagnostics  Id number (in available diagnostics list)
                0031 C     ip      :: diagnostics  pointer to storage array
                0032 C     im      :: counter-mate pointer to storage array
                0033 C     lm      :: index in the averageCycle
bc85e7ab92 Jean*0034 C     NrMax   :: 3rd dimension of input/output arrays
48a533dac6 Jean*0035 C     qtmp1   :: horizontal velocity input diag., u-component
                0036 C     qtmp2   :: horizontal velocity input diag., v-component
d5ada4102b Jean*0037 C     myTime  :: current time of simulation (s)
                0038 C     myIter  :: current iteration number
                0039 C     myThid  :: my Thread Id number
                0040       INTEGER listId, md, ndId, ip, im, lm
bc85e7ab92 Jean*0041       INTEGER NrMax
d5ada4102b Jean*0042       _RL     qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
                0043       _RL     qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
                0044       _RL     myTime
                0045       INTEGER myIter, myThid
48a533dac6 Jean*0046
                0047 C     !OUTPUT PARAMETERS:
                0048 C     qtmp1   :: horizontal-velocity potential
                0049 C     qtmp2   :: horizontal-velocity stream-function
d5ada4102b Jean*0050 CEOP
                0051
                0052 C     !FUNCTIONS:
bc85e7ab92 Jean*0053       INTEGER  ILNBLNK
                0054       EXTERNAL ILNBLNK
d5ada4102b Jean*0055
                0056 C     !LOCAL VARIABLES:
bc85e7ab92 Jean*0057 C     bi, bj  :: tile indices
                0058 C     i,j,k   :: loop indices
                0059 C     uTrans  :: horizontal transport, u-component
                0060 C     vTrans  :: horizontal transport, u-component
                0061 C     psiVel  :: horizontal stream-function
                0062 C     psiLoc  :: horizontal stream-function at special location
d5ada4102b Jean*0063       INTEGER bi, bj
bc85e7ab92 Jean*0064       INTEGER i, j, k
d5ada4102b Jean*0065
48a533dac6 Jean*0066       INTEGER ks
d53cbeab8e Jean*0067       INTEGER numIters, nIterMin
d5ada4102b Jean*0068       LOGICAL normaliseMatrice, diagNormaliseRHS
d53cbeab8e Jean*0069       _RL  residCriter, firstResidual, minResidual, lastResidual
d5ada4102b Jean*0070       _RL  a2dMax, a2dNorm
bdbdc9dd8f Jean*0071       _RL  rhsMax, b2dNorm, x2dNorm
d5ada4102b Jean*0072       _RS  aW2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0073       _RS  aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0074       _RL  b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0075       _RL  x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
bc85e7ab92 Jean*0076       _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0077       _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0078       _RL psiVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0079       _RL psiLoc(2)
                0080       INTEGER iL
                0081       CHARACTER*(MAX_LEN_FNAM) dataFName
                0082       CHARACTER*(MAX_LEN_MBUF) msgBuf
d5ada4102b Jean*0083
                0084 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0085
5ce5806717 Jean*0086 #ifdef ALLOW_DEBUG
                0087       IF (debugMode) CALL DEBUG_ENTER('DIAGNOSTICS_CALC_PHIVEL',myThid)
                0088 #endif
                0089
bdbdc9dd8f Jean*0090       DO ks = 1,nlevels(listId)
48a533dac6 Jean*0091         k = NINT(levs(ks,listId))
d5ada4102b Jean*0092 C--   Solve for velocity potential for each level:
                0093
                0094         a2dMax = 0. _d 0
                0095         rhsMax = 0. _d 0
                0096         DO bj = myByLo(myThid), myByHi(myThid)
                0097          DO bi = myBxLo(myThid), myBxHi(myThid)
                0098 C-    Initialise fist guess & RHS
a63b8f5615 Jean*0099            DO j = 1-OLy,sNy+OLy
                0100             DO i = 1-OLx,sNx+OLx
d5ada4102b Jean*0101               b2d(i,j,bi,bj) = 0.
                0102               x2d(i,j,bi,bj) = 0.
                0103             ENDDO
                0104            ENDDO
                0105 C-    calculate cg2d matrix:
c8934cd63f Jean*0106 C     Note: Here, at Open-Boundary location, we keep non-zero aW & aS (using
                0107 C     maskInW & maskInS) whereas in S/R CG2D they are zero (product of maskInC)
d5ada4102b Jean*0108            DO j = 1,sNy+1
                0109             DO i = 1,sNx+1
                0110               aW2d(i,j,bi,bj) = dyG(i,j,bi,bj)*recip_dxC(i,j,bi,bj)
                0111      &                         *drF(k)*hFacW(i,j,k,bi,bj)
c8934cd63f Jean*0112      &                         *maskInW(i,j,bi,bj)
d5ada4102b Jean*0113               aS2d(i,j,bi,bj) = dxG(i,j,bi,bj)*recip_dyC(i,j,bi,bj)
                0114      &                         *drF(k)*hFacS(i,j,k,bi,bj)
c8934cd63f Jean*0115      &                         *maskInS(i,j,bi,bj)
d5ada4102b Jean*0116               a2dMax = MAX(a2dMax,aW2d(i,j,bi,bj))
                0117               a2dMax = MAX(a2dMax,aS2d(i,j,bi,bj))
                0118             ENDDO
                0119            ENDDO
                0120
bdbdc9dd8f Jean*0121 C-    calculate horizontal transport
d5ada4102b Jean*0122            DO j = 1,sNy+1
                0123             DO i = 1,sNx+1
bc85e7ab92 Jean*0124               uTrans(i,j,bi,bj) = dyG(i,j,bi,bj)*drF(k)
                0125      &                          *qtmp1(i,j,ks,bi,bj)*maskInW(i,j,bi,bj)
                0126               vTrans(i,j,bi,bj) = dxG(i,j,bi,bj)*drF(k)
                0127      &                          *qtmp2(i,j,ks,bi,bj)*maskInS(i,j,bi,bj)
d5ada4102b Jean*0128             ENDDO
                0129            ENDDO
bdbdc9dd8f Jean*0130 C-   end bi,bj loops
                0131          ENDDO
                0132         ENDDO
                0133
                0134 C-    calculate RHS = rAc*Div(uVel,vVel):
                0135         DO bj = myByLo(myThid), myByHi(myThid)
                0136          DO bi = myBxLo(myThid), myBxHi(myThid)
d5ada4102b Jean*0137            DO j = 1,sNy
                0138             DO i = 1,sNx
bc85e7ab92 Jean*0139               b2d(i,j,bi,bj) = (
                0140      &                    ( uTrans(i+1,j,bi,bj) - uTrans(i,j,bi,bj) )
                0141      &                  + ( vTrans(i,j+1,bi,bj) - vTrans(i,j,bi,bj) )
d5ada4102b Jean*0142      &                         )*maskInC(i,j,bi,bj)
                0143             ENDDO
                0144            ENDDO
                0145          ENDDO
                0146         ENDDO
                0147
c8934cd63f Jean*0148 #ifdef ALLOW_OBCS
2b093fab14 Jean*0149 C  There is ambiguity in splitting OB cross flow into divergent (grad.Phi)
                0150 C    contribution and rotational (rot.Psi) contribution:
                0151 C  a) In most cases, will interpret most of the OB cross flow (except
                0152 C     the net inflow which has to come from grad.Phi) as non divergent
                0153 C     and only keep the divergence associated with the net inflow
                0154 C     (assuming here uniform distribution along the OB section;
                0155 C     This processing must be done for each domain-connected section
                0156 C     of the OB (using pkg/obcs OB[N,S,E,W]_connect Id) otherwise the
                0157 C     solver will not converge.
                0158 C  b) When setting a null domain-connected Id for some OB section,
                0159 C     we can recover the other extreme where the OB cross flow is
                0160 C     interpreted as a pure divergent (grad.Phi) contribution (-> constant
                0161 C      Psi along this section). This is done by keeping the full RHS just
                0162 C     outside OB (i.e., where tracer OBCS are applied.
c8934cd63f Jean*0163         IF ( useOBCS ) THEN
                0164          CALL OBCS_DIAG_BALANCE(
                0165      U             b2d,
                0166      I             uTrans, vTrans, k,
                0167      I             myTime, myIter, myThid )
                0168         ENDIF
                0169 #endif /* ALLOW_OBCS */
                0170
d5ada4102b Jean*0171 C-    Normalise Matrice & RHS :
d53cbeab8e Jean*0172         diagNormaliseRHS = diagCG_resTarget.GT.0.
d5ada4102b Jean*0173         normaliseMatrice = .TRUE.
                0174         diagNormaliseRHS = .TRUE.
2b093fab14 Jean*0175         IF ( diagNormaliseRHS ) THEN
                0176          DO bj = myByLo(myThid), myByHi(myThid)
                0177           DO bi = myBxLo(myThid), myBxHi(myThid)
                0178            DO j = 1,sNy
                0179             DO i = 1,sNx
41c4545f8f Jean*0180               rhsMax = MAX(ABS(b2d(i,j,bi,bj)),rhsMax)
2b093fab14 Jean*0181             ENDDO
                0182            ENDDO
                0183           ENDDO
                0184          ENDDO
                0185         ENDIF
d5ada4102b Jean*0186         a2dNorm = 1. _d 0
bdbdc9dd8f Jean*0187         b2dNorm = 1. _d 0
                0188         x2dNorm = 1. _d 0
d5ada4102b Jean*0189         IF ( normaliseMatrice ) THEN
                0190           _GLOBAL_MAX_RL( a2dMax, myThid )
49d8041f88 Jean*0191           IF ( a2dMax .GT. 0. _d 0 ) a2dNorm = 1. _d 0/a2dMax
bdbdc9dd8f Jean*0192           b2dNorm = a2dNorm
d5ada4102b Jean*0193         ENDIF
                0194         IF ( diagNormaliseRHS ) THEN
                0195           _GLOBAL_MAX_RL( rhsMax, myThid )
bdbdc9dd8f Jean*0196           IF ( rhsMax .GT. 0. _d 0 ) THEN
                0197             b2dNorm = 1. _d 0/rhsMax
                0198             x2dNorm = a2dNorm*rhsMax
                0199           ENDIF
d53cbeab8e Jean*0200           residCriter = diagCG_resTarget
d5ada4102b Jean*0201         ELSE
d53cbeab8e Jean*0202           residCriter = a2dNorm * ABS(diagCG_resTarget)
a63b8f5615 Jean*0203      &                          * globalArea / deltaTMom
d5ada4102b Jean*0204         ENDIF
                0205         IF ( normaliseMatrice .OR. diagNormaliseRHS ) THEN
                0206           DO bj = myByLo(myThid), myByHi(myThid)
                0207            DO bi = myBxLo(myThid), myBxHi(myThid)
                0208             DO j = 1,sNy+1
                0209              DO i = 1,sNx+1
                0210               aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*a2dNorm
                0211               aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*a2dNorm
bdbdc9dd8f Jean*0212               b2d(i,j,bi,bj)  = b2d(i,j,bi,bj) *b2dNorm
                0213 c             x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) /x2dNorm
d5ada4102b Jean*0214              ENDDO
                0215             ENDDO
                0216            ENDDO
                0217           ENDDO
                0218         ENDIF
                0219
41c4545f8f Jean*0220         IF ( diag_dBugLevel.GE.debLevA .AND. k.EQ.1 ) THEN
d5ada4102b Jean*0221           _BEGIN_MASTER( myThid )
                0222           WRITE(standardMessageUnit,'(A,I9,2(A,1P1E13.6),A,1P1E9.2)')
bdbdc9dd8f Jean*0223      &     ' diag_cg2d (it=', myIter,') a2dNorm,x2dNorm=', a2dNorm,
                0224      &     ' ,', x2dNorm, ' ; Criter=', residCriter
d5ada4102b Jean*0225           _END_MASTER( myThid )
                0226         ENDIF
                0227
d53cbeab8e Jean*0228         numIters = diagCG_maxIters
d5ada4102b Jean*0229         CALL DIAG_CG2D(
                0230      I                aW2d, aS2d, b2d,
a63b8f5615 Jean*0231      I                diagCG_pcOffDFac, residCriter,
d53cbeab8e Jean*0232      O                firstResidual, minResidual, lastResidual,
d5ada4102b Jean*0233      U                x2d, numIters,
d53cbeab8e Jean*0234      O                nIterMin,
                0235      I                diagCG_prtResFrq, myThid )
d5ada4102b Jean*0236
41c4545f8f Jean*0237         IF ( diag_dBugLevel.GE.debLevA ) THEN
d5ada4102b Jean*0238           _BEGIN_MASTER( myThid )
d53cbeab8e Jean*0239           WRITE(standardMessageUnit,'(A,I4,A,2I6,A,1P3E14.7)')
                0240      &    ' diag_cg2d : k=', k, ' , it=', nIterMin, numIters,
                0241      &    ' ; ini,min,last_Res=',firstResidual,minResidual,lastResidual
d5ada4102b Jean*0242           _END_MASTER( myThid )
                0243         ENDIF
                0244
bc85e7ab92 Jean*0245         _EXCH_XY_RL( x2d, myThid )
d5ada4102b Jean*0246
                0247 C-    Un-normalise the answer
                0248         IF (diagNormaliseRHS) THEN
                0249           DO bj = myByLo(myThid), myByHi(myThid)
                0250            DO bi = myBxLo(myThid), myBxHi(myThid)
a63b8f5615 Jean*0251             DO j = 1-OLy,sNy+OLy
                0252              DO i = 1-OLx,sNx+OLx
bdbdc9dd8f Jean*0253               x2d(i,j,bi,bj) =  x2d(i,j,bi,bj)*x2dNorm
d5ada4102b Jean*0254              ENDDO
                0255             ENDDO
                0256            ENDDO
                0257           ENDDO
                0258         ENDIF
                0259
bc85e7ab92 Jean*0260 C-    Compte divergence-free transport:
                0261         DO bj = myByLo(myThid), myByHi(myThid)
                0262          DO bi = myBxLo(myThid), myBxHi(myThid)
                0263           DO j = 1,sNy+1
                0264            DO i = 1,sNx+1
                0265             uTrans(i,j,bi,bj) = uTrans(i,j,bi,bj)
                0266      &                        - ( x2d(i,j,bi,bj) - x2d(i-1,j,bi,bj) )
                0267      &                         *recip_dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
                0268      &                         *drF(k)*hFacW(i,j,k,bi,bj)
c8934cd63f Jean*0269      &                         *maskInW(i,j,bi,bj)
bc85e7ab92 Jean*0270             vTrans(i,j,bi,bj) = vTrans(i,j,bi,bj)
                0271      &                        - ( x2d(i,j,bi,bj) - x2d(i,j-1,bi,bj) )
                0272      &                         *recip_dyC(i,j,bi,bj)*dxG(i,j,bi,bj)
                0273      &                         *drF(k)*hFacS(i,j,k,bi,bj)
c8934cd63f Jean*0274      &                         *maskInS(i,j,bi,bj)
bc85e7ab92 Jean*0275            ENDDO
                0276           ENDDO
                0277          ENDDO
                0278         ENDDO
                0279         CALL DIAG_CALC_PSIVEL(
                0280      I                         k, iPsi0, jPsi0, uTrans, vTrans,
                0281      O                         psiVel, psiLoc,
bdbdc9dd8f Jean*0282      I                         prtFirstCall, myTime, myIter, myThid )
bc85e7ab92 Jean*0283
bdbdc9dd8f Jean*0284         _BEGIN_MASTER( myThid)
bc85e7ab92 Jean*0285         IF ( useCubedSphereExchange .AND.
                0286      &       diag_mdsio .AND. myProcId.EQ.0 ) THEN
                0287 C-      Missing-corner value are not written in MDS output file
                0288 C       Write separately these 2 values (should be part of DIAGNOSTICS_OUT)
                0289          IF ( diagLoc_ioUnit.EQ.0 ) THEN
                0290           CALL MDSFINDUNIT( diagLoc_ioUnit, myThid )
                0291           WRITE(dataFName,'(2A,I10.10,A)')
                0292      &         'diags_CScorners', '.', nIter0, '.txt'
                0293           OPEN( diagLoc_ioUnit, FILE=dataFName, STATUS='unknown' )
                0294           iL = ILNBLNK(dataFName)
                0295           WRITE(msgBuf,'(2A,I6,2A)') 'DIAGNOSTICS_CALC_PHIVEL: ',
                0296      &         'open unit=',diagLoc_ioUnit, ', file: ',dataFName(1:iL)
                0297           CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
                0298      &                        SQUEEZE_RIGHT, myThid )
                0299          ENDIF
                0300          IF ( diagLoc_ioUnit.GT.0 ) THEN
                0301           WRITE(diagLoc_ioUnit,'(1P2E18.10,A,2I4,I8,A,2I4,I6,2A)')
                0302      &      psiLoc, ' #', k, lm, myIter,
                0303      &      ' :',listId, md, ndId, ' ', cdiag(ndId)
                0304 C-       check accuracy (f1.SW-corner = f6.NW-corner = f5-NE-corner)
                0305 c         WRITE(0,'(1P2E18.10,A,2I4,I8)')
                0306 c    &         psiVel(1,1+sNy,nSx,nSy)- psiVel(1,1,1,1),
                0307 c    &     psiVel(1+sNx,1+sNy,nSx-1,nSy)-psiVel(1,1,1,1),
                0308 c    &      ' #', k, lm, myIter
                0309          ENDIF
                0310         ENDIF
bdbdc9dd8f Jean*0311         IF ( prtFirstCall ) prtFirstCall = .FALSE.
                0312         _END_MASTER( myThid)
bc85e7ab92 Jean*0313
                0314 C-    Put the results back in qtmp[1,2]
                0315         DO bj = myByLo(myThid), myByHi(myThid)
                0316          DO bi = myBxLo(myThid), myBxHi(myThid)
                0317           DO j = 1,sNy+1
                0318            DO i = 1,sNx+1
                0319               qtmp1(i,j,ks,bi,bj) =  x2d(i,j,bi,bj)
                0320               qtmp2(i,j,ks,bi,bj) =  psiVel(i,j,bi,bj)
                0321            ENDDO
                0322           ENDDO
                0323          ENDDO
                0324         ENDDO
                0325
d5ada4102b Jean*0326       ENDDO
                0327
5ce5806717 Jean*0328 #ifdef ALLOW_DEBUG
                0329       IF (debugMode) CALL DEBUG_LEAVE('DIAGNOSTICS_CALC_PHIVEL',myThid)
                0330 #endif
                0331
d5ada4102b Jean*0332       RETURN
                0333       END