C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_calc_phivel.F,v 1.11 2014/11/25 01:12:11 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 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