pkg/diagnostics/diagnostics_calc_phivel.F

C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_calc_phivel.F,v 1.10 2014/11/21 20:59:06 jmc Exp $
C $Name:  $

#include "DIAG_OPTIONS.h"

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP 0
C     !ROUTINE: DIAGNOSTICS_CALC_PHIVEL

C     !INTERFACE:
      SUBROUTINE DIAGNOSTICS_CALC_PHIVEL(
     I                        listId, md, ndId, ip, im, lm,
     I                        NrMax,
     U                        qtmp1, qtmp2,
     I                        myTime, myIter, myThid )

C     !DESCRIPTION:
C     Compute Velocity Potential and Velocity Stream-Function

C     !USES:
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "DIAGNOSTICS_SIZE.h"
#include "DIAGNOSTICS.h"
#include "DIAGNOSTICS_CALC.h"

C     !INPUT PARAMETERS:
C     listId  :: Diagnostics list number being written
C     md      :: field number in the list "listId".
C     ndId    :: diagnostics  Id number (in available diagnostics list)
C     ip      :: diagnostics  pointer to storage array
C     im      :: counter-mate pointer to storage array
C     lm      :: index in the averageCycle
C     NrMax   :: 3rd dimension of input/output arrays
C     qtmp1   :: horizontal velocity input diag., u-component
C     qtmp2   :: horizontal velocity input diag., v-component
C     myTime  :: current time of simulation (s)
C     myIter  :: current iteration number
C     myThid  :: my Thread Id number
      INTEGER listId, md, ndId, ip, im, lm
      INTEGER NrMax
      _RL     qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
      _RL     qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
      _RL     myTime
      INTEGER myIter, myThid

C     !OUTPUT PARAMETERS:
C     qtmp1   :: horizontal-velocity potential
C     qtmp2   :: horizontal-velocity stream-function
CEOP

C     !FUNCTIONS:
      INTEGER  ILNBLNK
      EXTERNAL ILNBLNK

C     !LOCAL VARIABLES:
C     bi, bj  :: tile indices
C     i,j,k   :: loop indices
C     uTrans  :: horizontal transport, u-component
C     vTrans  :: horizontal transport, u-component
C     psiVel  :: horizontal stream-function
C     psiLoc  :: horizontal stream-function at special location
      INTEGER bi, bj
      INTEGER i, j, k

      INTEGER ks
      INTEGER numIters, nIterMin
      LOGICAL normaliseMatrice, diagNormaliseRHS
      _RL  residCriter, firstResidual, minResidual, lastResidual
      _RL  a2dMax, a2dNorm
      _RL  rhsMax, b2dNorm, x2dNorm
      _RS  aW2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS  aS2d(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL  b2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL  x2d (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL psiVel(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL psiLoc(2)
      INTEGER iL
      CHARACTER*(MAX_LEN_FNAM) dataFName
      CHARACTER*(MAX_LEN_MBUF) msgBuf

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_ENTER('DIAGNOSTICS_CALC_PHIVEL',myThid)
#endif

      DO ks = 1,nlevels(listId)
        k = NINT(levs(ks,listId))
C--   Solve for velocity potential for each level:

        a2dMax = 0. _d 0
        rhsMax = 0. _d 0
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
C-    Initialise fist guess & RHS
           DO j = 1-OLy,sNy+OLy
            DO i = 1-OLx,sNx+OLx
              b2d(i,j,bi,bj) = 0.
              x2d(i,j,bi,bj) = 0.
            ENDDO
           ENDDO
C-    calculate cg2d matrix:
C     Note: Here, at Open-Boundary location, we keep non-zero aW & aS (using
C     maskInW & maskInS) whereas in S/R CG2D they are zero (product of maskInC)
           DO j = 1,sNy+1
            DO i = 1,sNx+1
              aW2d(i,j,bi,bj) = dyG(i,j,bi,bj)*recip_dxC(i,j,bi,bj)
     &                         *drF(k)*hFacW(i,j,k,bi,bj)
     &                         *maskInW(i,j,bi,bj)
              aS2d(i,j,bi,bj) = dxG(i,j,bi,bj)*recip_dyC(i,j,bi,bj)
     &                         *drF(k)*hFacS(i,j,k,bi,bj)
     &                         *maskInS(i,j,bi,bj)
              a2dMax = MAX(a2dMax,aW2d(i,j,bi,bj))
              a2dMax = MAX(a2dMax,aS2d(i,j,bi,bj))
            ENDDO
           ENDDO

C-    calculate horizontal transport
           DO j = 1,sNy+1
            DO i = 1,sNx+1
              uTrans(i,j,bi,bj) = dyG(i,j,bi,bj)*drF(k)
     &                          *qtmp1(i,j,ks,bi,bj)*maskInW(i,j,bi,bj)
              vTrans(i,j,bi,bj) = dxG(i,j,bi,bj)*drF(k)
     &                          *qtmp2(i,j,ks,bi,bj)*maskInS(i,j,bi,bj)
            ENDDO
           ENDDO
C-   end bi,bj loops
         ENDDO
        ENDDO

C-    calculate RHS = rAc*Div(uVel,vVel):
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
           DO j = 1,sNy
            DO i = 1,sNx
              b2d(i,j,bi,bj) = (
     &                    ( uTrans(i+1,j,bi,bj) - uTrans(i,j,bi,bj) )
     &                  + ( vTrans(i,j+1,bi,bj) - vTrans(i,j,bi,bj) )
     &                         )*maskInC(i,j,bi,bj)
            ENDDO
           ENDDO
         ENDDO
        ENDDO

#ifdef ALLOW_OBCS
C  There is ambiguity in splitting OB cross flow into divergent (grad.Phi)
C    contribution and rotational (rot.Psi) contribution:
C  a) In most cases, will interpret most of the OB cross flow (except
C     the net inflow which has to come from grad.Phi) as non divergent
C     and only keep the divergence associated with the net inflow
C     (assuming here uniform distribution along the OB section;
C     This processing must be done for each domain-connected section
C     of the OB (using pkg/obcs OB[N,S,E,W]_connect Id) otherwise the
C     solver will not converge.
C  b) When setting a null domain-connected Id for some OB section,
C     we can recover the other extreme where the OB cross flow is
C     interpreted as a pure divergent (grad.Phi) contribution (-> constant
C      Psi along this section). This is done by keeping the full RHS just
C     outside OB (i.e., where tracer OBCS are applied.
        IF ( useOBCS ) THEN
         CALL OBCS_DIAG_BALANCE(
     U             b2d,
     I             uTrans, vTrans, k,
     I             myTime, myIter, myThid )
        ENDIF
#endif /* ALLOW_OBCS */

C-    Normalise Matrice & RHS :
        diagNormaliseRHS = diagCG_resTarget.GT.0.
        normaliseMatrice = .TRUE.
        diagNormaliseRHS = .TRUE.
        IF ( diagNormaliseRHS ) THEN
         DO bj = myByLo(myThid), myByHi(myThid)
          DO bi = myBxLo(myThid), myBxHi(myThid)
           DO j = 1,sNy
            DO i = 1,sNx
              rhsMax = MAX(ABS(b2d(I,J,bi,bj)),rhsMax)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDIF
        a2dNorm = 1. _d 0
        b2dNorm = 1. _d 0
        x2dNorm = 1. _d 0
        IF ( normaliseMatrice ) THEN
          _GLOBAL_MAX_RL( a2dMax, myThid )
          IF ( a2dMax .GT. 0. _d 0 ) a2dNorm = 1. _d 0/a2dMax
          b2dNorm = a2dNorm
        ENDIF
        IF ( diagNormaliseRHS ) THEN
          _GLOBAL_MAX_RL( rhsMax, myThid )
          IF ( rhsMax .GT. 0. _d 0 ) THEN
            b2dNorm = 1. _d 0/rhsMax
            x2dNorm = a2dNorm*rhsMax
          ENDIF
          residCriter = diagCG_resTarget
        ELSE
          residCriter = a2dNorm * ABS(diagCG_resTarget)
     &                          * globalArea / deltaTMom
        ENDIF
        IF ( normaliseMatrice .OR. diagNormaliseRHS ) THEN
          DO bj = myByLo(myThid), myByHi(myThid)
           DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j = 1,sNy+1
             DO i = 1,sNx+1
              aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*a2dNorm
              aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*a2dNorm
              b2d(i,j,bi,bj)  = b2d(i,j,bi,bj) *b2dNorm
c             x2d(i,j,bi,bj) =  x2d(i,j,bi,bj) /x2dNorm
             ENDDO
            ENDDO
           ENDDO
          ENDDO
        ENDIF

        IF ( debugLevel.GE.debLevA .AND. k.EQ.1 ) THEN
          _BEGIN_MASTER( myThid )
          WRITE(standardMessageUnit,'(A,I9,2(A,1P1E13.6),A,1P1E9.2)')
     &     ' diag_cg2d (it=', myIter,') a2dNorm,x2dNorm=', a2dNorm,
     &     ' ,', x2dNorm, ' ; Criter=', residCriter
          _END_MASTER( myThid )
        ENDIF

        numIters = diagCG_maxIters
        CALL DIAG_CG2D(
     I                aW2d, aS2d, b2d,
     I                diagCG_pcOffDFac, residCriter,
     O                firstResidual, minResidual, lastResidual,
     U                x2d, numIters,
     O                nIterMin,
     I                diagCG_prtResFrq, myThid )

        IF ( debugLevel.GE.debLevA ) THEN
          _BEGIN_MASTER( myThid )
          WRITE(standardMessageUnit,'(A,I4,A,2I6,A,1P3E14.7)')
     &    ' diag_cg2d : k=', k, ' , it=', nIterMin, numIters,
     &    ' ; ini,min,last_Res=',firstResidual,minResidual,lastResidual
          _END_MASTER( myThid )
        ENDIF

        _EXCH_XY_RL( x2d, myThid )

C-    Un-normalise the answer
        IF (diagNormaliseRHS) THEN
          DO bj = myByLo(myThid), myByHi(myThid)
           DO bi = myBxLo(myThid), myBxHi(myThid)
            DO j = 1-OLy,sNy+OLy
             DO i = 1-OLx,sNx+OLx
              x2d(i,j,bi,bj) =  x2d(i,j,bi,bj)*x2dNorm
             ENDDO
            ENDDO
           ENDDO
          ENDDO
        ENDIF

C-    Compte divergence-free transport:
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j = 1,sNy+1
           DO i = 1,sNx+1
            uTrans(i,j,bi,bj) = uTrans(i,j,bi,bj)
     &                        - ( x2d(i,j,bi,bj) - x2d(i-1,j,bi,bj) )
     &                         *recip_dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
     &                         *drF(k)*hFacW(i,j,k,bi,bj)
     &                         *maskInW(i,j,bi,bj)
            vTrans(i,j,bi,bj) = vTrans(i,j,bi,bj)
     &                        - ( x2d(i,j,bi,bj) - x2d(i,j-1,bi,bj) )
     &                         *recip_dyC(i,j,bi,bj)*dxG(i,j,bi,bj)
     &                         *drF(k)*hFacS(i,j,k,bi,bj)
     &                         *maskInS(i,j,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDDO
        CALL DIAG_CALC_PSIVEL(
     I                         k, iPsi0, jPsi0, uTrans, vTrans,
     O                         psiVel, psiLoc,
     I                         prtFirstCall, myTime, myIter, myThid )

        _BEGIN_MASTER( myThid)
        IF ( useCubedSphereExchange .AND.
     &       diag_mdsio .AND. myProcId.EQ.0 ) THEN
C-      Missing-corner value are not written in MDS output file
C       Write separately these 2 values (should be part of DIAGNOSTICS_OUT)
         IF ( diagLoc_ioUnit.EQ.0 ) THEN
          CALL MDSFINDUNIT( diagLoc_ioUnit, myThid )
          WRITE(dataFName,'(2A,I10.10,A)')
     &         'diags_CScorners', '.', nIter0, '.txt'
          OPEN( diagLoc_ioUnit, FILE=dataFName, STATUS='unknown' )
          iL = ILNBLNK(dataFName)
          WRITE(msgBuf,'(2A,I6,2A)') 'DIAGNOSTICS_CALC_PHIVEL: ',
     &         'open unit=',diagLoc_ioUnit, ', file: ',dataFName(1:iL)
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
     &                        SQUEEZE_RIGHT, myThid )
         ENDIF
         IF ( diagLoc_ioUnit.GT.0 ) THEN
          WRITE(diagLoc_ioUnit,'(1P2E18.10,A,2I4,I8,A,2I4,I6,2A)')
     &      psiLoc, ' #', k, lm, myIter,
     &      ' :',listId, md, ndId, ' ', cdiag(ndId)
C-       check accuracy (f1.SW-corner = f6.NW-corner = f5-NE-corner)
c         WRITE(0,'(1P2E18.10,A,2I4,I8)')
c    &         psiVel(1,1+sNy,nSx,nSy)- psiVel(1,1,1,1),
c    &     psiVel(1+sNx,1+sNy,nSx-1,nSy)-psiVel(1,1,1,1),
c    &      ' #', k, lm, myIter
         ENDIF
        ENDIF
        IF ( prtFirstCall ) prtFirstCall = .FALSE.
        _END_MASTER( myThid)

C-    Put the results back in qtmp[1,2]
        DO bj = myByLo(myThid), myByHi(myThid)
         DO bi = myBxLo(myThid), myBxHi(myThid)
          DO j = 1,sNy+1
           DO i = 1,sNx+1
              qtmp1(i,j,ks,bi,bj) =  x2d(i,j,bi,bj)
              qtmp2(i,j,ks,bi,bj) =  psiVel(i,j,bi,bj)
           ENDDO
          ENDDO
         ENDDO
        ENDDO

      ENDDO

#ifdef ALLOW_DEBUG
      IF (debugMode) CALL DEBUG_LEAVE('DIAGNOSTICS_CALC_PHIVEL',myThid)
#endif

      RETURN
      END
1	jmc	1.11	C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_calc_phivel.F,v 1.10 2014/11/21 20:59:06 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "DIAG_OPTIONS.h"
5
6			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
7			CBOP 0
8			C !ROUTINE: DIAGNOSTICS_CALC_PHIVEL
9
10			C !INTERFACE:
11			SUBROUTINE DIAGNOSTICS_CALC_PHIVEL(
12			I listId, md, ndId, ip, im, lm,
13	jmc	1.4	I NrMax,
14	jmc	1.2	U qtmp1, qtmp2,
15	jmc	1.1	I myTime, myIter, myThid )
16
17			C !DESCRIPTION:
18	jmc	1.4	C Compute Velocity Potential and Velocity Stream-Function
19	jmc	1.1
20			C !USES:
21			IMPLICIT NONE
22			#include "SIZE.h"
23			#include "EEPARAMS.h"
24			#include "PARAMS.h"
25			#include "GRID.h"
26			#include "DIAGNOSTICS_SIZE.h"
27			#include "DIAGNOSTICS.h"
28	jmc	1.5	#include "DIAGNOSTICS_CALC.h"
29	jmc	1.1
30			C !INPUT PARAMETERS:
31			C listId :: Diagnostics list number being written
32			C md :: field number in the list "listId".
33			C ndId :: diagnostics Id number (in available diagnostics list)
34			C ip :: diagnostics pointer to storage array
35			C im :: counter-mate pointer to storage array
36			C lm :: index in the averageCycle
37	jmc	1.4	C NrMax :: 3rd dimension of input/output arrays
38	jmc	1.2	C qtmp1 :: horizontal velocity input diag., u-component
39			C qtmp2 :: horizontal velocity input diag., v-component
40	jmc	1.1	C myTime :: current time of simulation (s)
41			C myIter :: current iteration number
42			C myThid :: my Thread Id number
43			INTEGER listId, md, ndId, ip, im, lm
44	jmc	1.4	INTEGER NrMax
45	jmc	1.1	_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)
47			_RL myTime
48			INTEGER myIter, myThid
49	jmc	1.2
50			C !OUTPUT PARAMETERS:
51			C qtmp1 :: horizontal-velocity potential
52			C qtmp2 :: horizontal-velocity stream-function
53	jmc	1.1	CEOP
54
55			C !FUNCTIONS:
56	jmc	1.4	INTEGER ILNBLNK
57			EXTERNAL ILNBLNK
58	jmc	1.1
59			C !LOCAL VARIABLES:
60	jmc	1.4	C bi, bj :: tile indices
61			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
67	jmc	1.1	INTEGER i, j, k
68
69	jmc	1.2	INTEGER ks
70	jmc	1.6	INTEGER numIters, nIterMin
71	jmc	1.1	LOGICAL normaliseMatrice, diagNormaliseRHS
72	jmc	1.6	_RL residCriter, firstResidual, minResidual, lastResidual
73	jmc	1.1	_RL a2dMax, a2dNorm
74	jmc	1.7	_RL rhsMax, b2dNorm, x2dNorm
75	jmc	1.1	_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)
77			_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)
79	jmc	1.4	_RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
80			_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	jmc	1.1
87			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
88
89	jmc	1.8	#ifdef ALLOW_DEBUG
90			IF (debugMode) CALL DEBUG_ENTER('DIAGNOSTICS_CALC_PHIVEL',myThid)
91			#endif
92
93	jmc	1.7	DO ks = 1,nlevels(listId)
94	jmc	1.2	k = NINT(levs(ks,listId))
95	jmc	1.1	C-- Solve for velocity potential for each level:
96
97			a2dMax = 0. _d 0
98			rhsMax = 0. _d 0
99			DO bj = myByLo(myThid), myByHi(myThid)
100			DO bi = myBxLo(myThid), myBxHi(myThid)
101			C- Initialise fist guess & RHS
102	jmc	1.9	DO j = 1-OLy,sNy+OLy
103			DO i = 1-OLx,sNx+OLx
104	jmc	1.1	b2d(i,j,bi,bj) = 0.
105			x2d(i,j,bi,bj) = 0.
106			ENDDO
107			ENDDO
108			C- calculate cg2d matrix:
109	jmc	1.10	C Note: Here, at Open-Boundary location, we keep non-zero aW & aS (using
110			C maskInW & maskInS) whereas in S/R CG2D they are zero (product of maskInC)
111	jmc	1.1	DO j = 1,sNy+1
112			DO i = 1,sNx+1
113			aW2d(i,j,bi,bj) = dyG(i,j,bi,bj)*recip_dxC(i,j,bi,bj)
114			& drF(k)hFacW(i,j,k,bi,bj)
115	jmc	1.10	& *maskInW(i,j,bi,bj)
116	jmc	1.1	aS2d(i,j,bi,bj) = dxG(i,j,bi,bj)*recip_dyC(i,j,bi,bj)
117			& drF(k)hFacS(i,j,k,bi,bj)
118	jmc	1.10	& *maskInS(i,j,bi,bj)
119	jmc	1.1	a2dMax = MAX(a2dMax,aW2d(i,j,bi,bj))
120			a2dMax = MAX(a2dMax,aS2d(i,j,bi,bj))
121			ENDDO
122			ENDDO
123
124	jmc	1.7	C- calculate horizontal transport
125	jmc	1.1	DO j = 1,sNy+1
126			DO i = 1,sNx+1
127	jmc	1.4	uTrans(i,j,bi,bj) = dyG(i,j,bi,bj)*drF(k)
128			& qtmp1(i,j,ks,bi,bj)maskInW(i,j,bi,bj)
129			vTrans(i,j,bi,bj) = dxG(i,j,bi,bj)*drF(k)
130			& qtmp2(i,j,ks,bi,bj)maskInS(i,j,bi,bj)
131	jmc	1.1	ENDDO
132			ENDDO
133	jmc	1.7	C- end bi,bj loops
134			ENDDO
135			ENDDO
136
137			C- calculate RHS = rAc*Div(uVel,vVel):
138			DO bj = myByLo(myThid), myByHi(myThid)
139			DO bi = myBxLo(myThid), myBxHi(myThid)
140	jmc	1.1	DO j = 1,sNy
141			DO i = 1,sNx
142	jmc	1.4	b2d(i,j,bi,bj) = (
143			& ( uTrans(i+1,j,bi,bj) - uTrans(i,j,bi,bj) )
144			& + ( vTrans(i,j+1,bi,bj) - vTrans(i,j,bi,bj) )
145	jmc	1.1	& )*maskInC(i,j,bi,bj)
146			ENDDO
147			ENDDO
148			ENDDO
149			ENDDO
150
151	jmc	1.10	#ifdef ALLOW_OBCS
152	jmc	1.11	C There is ambiguity in splitting OB cross flow into divergent (grad.Phi)
153			C contribution and rotational (rot.Psi) contribution:
154			C a) In most cases, will interpret most of the OB cross flow (except
155			C the net inflow which has to come from grad.Phi) as non divergent
156			C and only keep the divergence associated with the net inflow
157			C (assuming here uniform distribution along the OB section;
158			C This processing must be done for each domain-connected section
159			C of the OB (using pkg/obcs OB[N,S,E,W]_connect Id) otherwise the
160			C solver will not converge.
161			C b) When setting a null domain-connected Id for some OB section,
162			C we can recover the other extreme where the OB cross flow is
163			C interpreted as a pure divergent (grad.Phi) contribution (-> constant
164			C Psi along this section). This is done by keeping the full RHS just
165			C outside OB (i.e., where tracer OBCS are applied.
166	jmc	1.10	IF ( useOBCS ) THEN
167			CALL OBCS_DIAG_BALANCE(
168			U b2d,
169			I uTrans, vTrans, k,
170			I myTime, myIter, myThid )
171			ENDIF
172			#endif /* ALLOW_OBCS */
173
174	jmc	1.1	C- Normalise Matrice & RHS :
175	jmc	1.6	diagNormaliseRHS = diagCG_resTarget.GT.0.
176	jmc	1.1	normaliseMatrice = .TRUE.
177			diagNormaliseRHS = .TRUE.
178	jmc	1.11	IF ( diagNormaliseRHS ) THEN
179			DO bj = myByLo(myThid), myByHi(myThid)
180			DO bi = myBxLo(myThid), myBxHi(myThid)
181			DO j = 1,sNy
182			DO i = 1,sNx
183			rhsMax = MAX(ABS(b2d(I,J,bi,bj)),rhsMax)
184			ENDDO
185			ENDDO
186			ENDDO
187			ENDDO
188			ENDIF
189	jmc	1.1	a2dNorm = 1. _d 0
190	jmc	1.7	b2dNorm = 1. _d 0
191			x2dNorm = 1. _d 0
192	jmc	1.1	IF ( normaliseMatrice ) THEN
193			_GLOBAL_MAX_RL( a2dMax, myThid )
194	jmc	1.5	IF ( a2dMax .GT. 0. _d 0 ) a2dNorm = 1. _d 0/a2dMax
195	jmc	1.7	b2dNorm = a2dNorm
196	jmc	1.1	ENDIF
197			IF ( diagNormaliseRHS ) THEN
198			_GLOBAL_MAX_RL( rhsMax, myThid )
199	jmc	1.7	IF ( rhsMax .GT. 0. _d 0 ) THEN
200			b2dNorm = 1. _d 0/rhsMax
201			x2dNorm = a2dNorm*rhsMax
202			ENDIF
203	jmc	1.6	residCriter = diagCG_resTarget
204	jmc	1.1	ELSE
205	jmc	1.6	residCriter = a2dNorm * ABS(diagCG_resTarget)
206	jmc	1.9	& * globalArea / deltaTMom
207	jmc	1.1	ENDIF
208			IF ( normaliseMatrice .OR. diagNormaliseRHS ) THEN
209			DO bj = myByLo(myThid), myByHi(myThid)
210			DO bi = myBxLo(myThid), myBxHi(myThid)
211			DO j = 1,sNy+1
212			DO i = 1,sNx+1
213			aW2d(i,j,bi,bj) = aW2d(i,j,bi,bj)*a2dNorm
214			aS2d(i,j,bi,bj) = aS2d(i,j,bi,bj)*a2dNorm
215	jmc	1.7	b2d(i,j,bi,bj) = b2d(i,j,bi,bj) *b2dNorm
216			c x2d(i,j,bi,bj) = x2d(i,j,bi,bj) /x2dNorm
217	jmc	1.1	ENDDO
218			ENDDO
219			ENDDO
220			ENDDO
221			ENDIF
222
223			IF ( debugLevel.GE.debLevA .AND. k.EQ.1 ) THEN
224			_BEGIN_MASTER( myThid )
225			WRITE(standardMessageUnit,'(A,I9,2(A,1P1E13.6),A,1P1E9.2)')
226	jmc	1.7	& ' diag_cg2d (it=', myIter,') a2dNorm,x2dNorm=', a2dNorm,
227			& ' ,', x2dNorm, ' ; Criter=', residCriter
228	jmc	1.1	_END_MASTER( myThid )
229			ENDIF
230
231	jmc	1.6	numIters = diagCG_maxIters
232	jmc	1.1	CALL DIAG_CG2D(
233			I aW2d, aS2d, b2d,
234	jmc	1.9	I diagCG_pcOffDFac, residCriter,
235	jmc	1.6	O firstResidual, minResidual, lastResidual,
236	jmc	1.1	U x2d, numIters,
237	jmc	1.6	O nIterMin,
238			I diagCG_prtResFrq, myThid )
239	jmc	1.1
240			IF ( debugLevel.GE.debLevA ) THEN
241			_BEGIN_MASTER( myThid )
242	jmc	1.6	WRITE(standardMessageUnit,'(A,I4,A,2I6,A,1P3E14.7)')
243			& ' diag_cg2d : k=', k, ' , it=', nIterMin, numIters,
244			& ' ; ini,min,last_Res=',firstResidual,minResidual,lastResidual
245	jmc	1.1	_END_MASTER( myThid )
246			ENDIF
247
248	jmc	1.4	_EXCH_XY_RL( x2d, myThid )
249	jmc	1.1
250			C- Un-normalise the answer
251			IF (diagNormaliseRHS) THEN
252			DO bj = myByLo(myThid), myByHi(myThid)
253			DO bi = myBxLo(myThid), myBxHi(myThid)
254	jmc	1.9	DO j = 1-OLy,sNy+OLy
255			DO i = 1-OLx,sNx+OLx
256	jmc	1.7	x2d(i,j,bi,bj) = x2d(i,j,bi,bj)*x2dNorm
257	jmc	1.1	ENDDO
258			ENDDO
259			ENDDO
260			ENDDO
261			ENDIF
262
263	jmc	1.4	C- Compte divergence-free transport:
264			DO bj = myByLo(myThid), myByHi(myThid)
265			DO bi = myBxLo(myThid), myBxHi(myThid)
266			DO j = 1,sNy+1
267			DO i = 1,sNx+1
268			uTrans(i,j,bi,bj) = uTrans(i,j,bi,bj)
269			& - ( x2d(i,j,bi,bj) - x2d(i-1,j,bi,bj) )
270			& recip_dxC(i,j,bi,bj)dyG(i,j,bi,bj)
271			& drF(k)hFacW(i,j,k,bi,bj)
272	jmc	1.10	& *maskInW(i,j,bi,bj)
273	jmc	1.4	vTrans(i,j,bi,bj) = vTrans(i,j,bi,bj)
274			& - ( x2d(i,j,bi,bj) - x2d(i,j-1,bi,bj) )
275			& recip_dyC(i,j,bi,bj)dxG(i,j,bi,bj)
276			& drF(k)hFacS(i,j,k,bi,bj)
277	jmc	1.10	& *maskInS(i,j,bi,bj)
278	jmc	1.4	ENDDO
279			ENDDO
280			ENDDO
281			ENDDO
282			CALL DIAG_CALC_PSIVEL(
283			I k, iPsi0, jPsi0, uTrans, vTrans,
284			O psiVel, psiLoc,
285	jmc	1.7	I prtFirstCall, myTime, myIter, myThid )
286	jmc	1.4
287	jmc	1.7	_BEGIN_MASTER( myThid)
288	jmc	1.4	IF ( useCubedSphereExchange .AND.
289			& diag_mdsio .AND. myProcId.EQ.0 ) THEN
290			C- Missing-corner value are not written in MDS output file
291			C Write separately these 2 values (should be part of DIAGNOSTICS_OUT)
292			IF ( diagLoc_ioUnit.EQ.0 ) THEN
293			CALL MDSFINDUNIT( diagLoc_ioUnit, myThid )
294			WRITE(dataFName,'(2A,I10.10,A)')
295			& 'diags_CScorners', '.', nIter0, '.txt'
296			OPEN( diagLoc_ioUnit, FILE=dataFName, STATUS='unknown' )
297			iL = ILNBLNK(dataFName)
298			WRITE(msgBuf,'(2A,I6,2A)') 'DIAGNOSTICS_CALC_PHIVEL: ',
299			& 'open unit=',diagLoc_ioUnit, ', file: ',dataFName(1:iL)
300			CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
301			& SQUEEZE_RIGHT, myThid )
302			ENDIF
303			IF ( diagLoc_ioUnit.GT.0 ) THEN
304			WRITE(diagLoc_ioUnit,'(1P2E18.10,A,2I4,I8,A,2I4,I6,2A)')
305			& psiLoc, ' #', k, lm, myIter,
306			& ' :',listId, md, ndId, ' ', cdiag(ndId)
307			C- check accuracy (f1.SW-corner = f6.NW-corner = f5-NE-corner)
308			c WRITE(0,'(1P2E18.10,A,2I4,I8)')
309			c & psiVel(1,1+sNy,nSx,nSy)- psiVel(1,1,1,1),
310			c & psiVel(1+sNx,1+sNy,nSx-1,nSy)-psiVel(1,1,1,1),
311			c & ' #', k, lm, myIter
312			ENDIF
313			ENDIF
314	jmc	1.7	IF ( prtFirstCall ) prtFirstCall = .FALSE.
315			_END_MASTER( myThid)
316	jmc	1.4
317			C- Put the results back in qtmp[1,2]
318			DO bj = myByLo(myThid), myByHi(myThid)
319			DO bi = myBxLo(myThid), myBxHi(myThid)
320			DO j = 1,sNy+1
321			DO i = 1,sNx+1
322			qtmp1(i,j,ks,bi,bj) = x2d(i,j,bi,bj)
323			qtmp2(i,j,ks,bi,bj) = psiVel(i,j,bi,bj)
324			ENDDO
325			ENDDO
326			ENDDO
327			ENDDO
328
329	jmc	1.1	ENDDO
330
331	jmc	1.8	#ifdef ALLOW_DEBUG
332			IF (debugMode) CALL DEBUG_LEAVE('DIAGNOSTICS_CALC_PHIVEL',myThid)
333			#endif
334
335	jmc	1.1	RETURN
336			END