Back to home page

MITgcm
 
 

    

File indexing completed on 2018-03-02 18:42:22 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
91672e10e3 Alis*0001 #include "MONITOR_OPTIONS.h"
0a272ce360 Jean*0002 
2741539ec0 Ed H*0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: MON_VORT3
                0006 
                0007 C     !INTERFACE:
0a272ce360 Jean*0008       SUBROUTINE MON_VORT3(
2741539ec0 Ed H*0009      I     myIter, myThid )
0a272ce360 Jean*0010 
2741539ec0 Ed H*0011 C     !DESCRIPTION:
                0012 C     Calculates stats for Vorticity (z-component).
                0013 
                0014 C     !USES:
                0015       IMPLICIT NONE
0a272ce360 Jean*0016 #include "SIZE.h"
                0017 #include "EEPARAMS.h"
                0018 #include "PARAMS.h"
                0019 #include "DYNVARS.h"
                0020 #include "MONITOR.h"
                0021 #include "GRID.h"
4749c74143 Alis*0022 #ifdef ALLOW_EXCH2
f9f661930b Jean*0023 #include "W2_EXCH2_SIZE.h"
404487e0d3 Andr*0024 #include "W2_EXCH2_TOPOLOGY.h"
4749c74143 Alis*0025 #endif /* ALLOW_EXCH2 */
0a272ce360 Jean*0026 
2741539ec0 Ed H*0027 C     !INPUT PARAMETERS:
0a272ce360 Jean*0028       INTEGER myIter, myThid
2741539ec0 Ed H*0029 CEOP
0a272ce360 Jean*0030 
2741539ec0 Ed H*0031 C     !LOCAL VARIABLES:
0a272ce360 Jean*0032       INTEGER bi,bj,i,j,k
                0033       INTEGER iMax,jMax
c22d41690c Jean*0034       _RL theVol, theArea, tmpVal, tmpAre, tmpVol
0a272ce360 Jean*0035       _RL theMin, theMax, theMean, theVar, volMean, volVar, theSD
4d2b0c1389 Jean*0036 c     _RL areaTile, volTile, sumTile, sqsTile, vSumTile, vSqsTile
                0037       _RL tileArea(nSx,nSy)
                0038       _RL tileVol (nSx,nSy)
                0039       _RL tileSum (nSx,nSy)
                0040       _RL tileVar (nSx,nSy)
                0041       _RL tileVSum(nSx,nSy)
                0042       _RL tileVSq (nSx,nSy)
0a272ce360 Jean*0043       _RL vort3(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0044       _RS hFacZ(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0045       _RL AZcorner
                0046 #ifdef MONITOR_TEST_HFACZ
                0047       _RL tmpFac
                0048       _RL etaFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0049 #endif
4d2b0c1389 Jean*0050       LOGICAL northWestCorner, northEastCorner
                0051       LOGICAL southWestCorner, southEastCorner
8996cf5a3c Jean*0052       INTEGER iG
                0053 #ifdef ALLOW_EXCH2
                0054       INTEGER myTile
                0055 #endif
0a272ce360 Jean*0056 
                0057       theMin = 1. _d 20
                0058       theMax =-1. _d 20
                0059       theArea= 0. _d 0
                0060       theMean= 0. _d 0
                0061       theVar = 0. _d 0
                0062       theVol = 0. _d 0
                0063       volMean= 0. _d 0
                0064       volVar = 0. _d 0
                0065       theSD  = 0. _d 0
                0066       AZcorner = 1. _d 0
                0067 
                0068       DO bj=myByLo(myThid),myByHi(myThid)
                0069        DO bi=myBxLo(myThid),myBxHi(myThid)
4d2b0c1389 Jean*0070          tileArea(bi,bj)= 0. _d 0
                0071          tileVol(bi,bj) = 0. _d 0
                0072          tileSum(bi,bj) = 0. _d 0
                0073          tileVar(bi,bj) = 0. _d 0
                0074          tileVSum(bi,bj)= 0. _d 0
                0075          tileVSq(bi,bj) = 0. _d 0
0a272ce360 Jean*0076 #ifdef MONITOR_TEST_HFACZ
                0077          tmpFac = 0.
4d2b0c1389 Jean*0078          IF( implicDiv2Dflow.GT.0 .AND. abEps.GT.-0.5 )
0a272ce360 Jean*0079      &    tmpFac = (0.5+abEps) / implicDiv2Dflow
                0080          DO j=1-Oly,sNy+Oly
                0081           DO i=1-Olx,sNx+Olx
                0082             etaFld(i,j) = etaH(i,j,bi,bj)
                0083      &          + tmpFac*(etaN(i,j,bi,bj)-etaH(i,j,bi,bj))
                0084           ENDDO
                0085          ENDDO
                0086 #endif
                0087         DO k=1,Nr
                0088 
                0089          iMax = sNx
                0090          jMax = sNy
                0091          DO j=1,sNy
                0092           DO i=1,sNx
                0093 #ifdef MONITOR_TEST_HFACZ
                0094 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0095 C-    Test various definitions of hFacZ (for 1 layer, flat bottom ocean):
                0096 c          hFacZ(i,j) = 1. +
4d2b0c1389 Jean*0097 c    &           0.25 _d 0*( etaFld(i-1,j-1)
0a272ce360 Jean*0098 c    &                     + etaFld( i ,j-1)
4d2b0c1389 Jean*0099 c    &                     + etaFld(i-1, j )
                0100 c    &                     + etaFld( i , j )
0a272ce360 Jean*0101 c    &                     )*recip_drF(k)
                0102 c          hFacZ(i,j) = 1. +
                0103 c    &           0.25 _d 0*( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
4d2b0c1389 Jean*0104 c    &                     + etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
0a272ce360 Jean*0105 c    &                     + etaFld(i-1, j )*rA(i-1, j ,bi,bj)
4d2b0c1389 Jean*0106 c    &                     + etaFld( i , j )*rA( i , j ,bi,bj)
0a272ce360 Jean*0107 c    &                     )*recip_drF(k)*recip_rAz(i,j,bi,bj)
                0108            hFacZ(i,j) = 1. + 0.125 _d 0*
                0109      &                   ( ( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
                0110      &                      +etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
4d2b0c1389 Jean*0111      &                     )*recip_rAw(i,j-1,bi,bj)
0a272ce360 Jean*0112      &                   + ( etaFld(i-1, j )*rA(i-1, j ,bi,bj)
4d2b0c1389 Jean*0113      &                      +etaFld( i , j )*rA( i , j ,bi,bj)
0a272ce360 Jean*0114      &                     )*recip_rAw(i, j ,bi,bj)
                0115      &                   + ( etaFld(i-1,j-1)*rA(i-1,j-1,bi,bj)
                0116      &                      +etaFld(i-1, j )*rA(i-1, j ,bi,bj)
4d2b0c1389 Jean*0117      &                     )*recip_rAs(i-1,j,bi,bj)
0a272ce360 Jean*0118      &                   + ( etaFld( i ,j-1)*rA( i ,j-1,bi,bj)
4d2b0c1389 Jean*0119      &                     + etaFld( i , j )*rA( i , j ,bi,bj)
0a272ce360 Jean*0120      &                     )*recip_rAs( i ,j,bi,bj)
                0121      &                   )*recip_drF(k)
                0122 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0123 #else
                0124 C-    Standard definition of hFac at vorticity point:
                0125            hFacZ(i,j) =
4d2b0c1389 Jean*0126      &           0.25 _d 0*( _hFacW(i,j-1,k,bi,bj)
                0127      &                     + _hFacW(i, j ,k,bi,bj)
                0128      &                     + _hFacS(i-1,j,k,bi,bj)
0a272ce360 Jean*0129      &                     + _hFacS( i ,j,k,bi,bj)
                0130      &                     )
4d2b0c1389 Jean*0131 #endif /* MONITOR_TEST_HFACZ */
0a272ce360 Jean*0132            vort3(i,j) = recip_rAz(i,j,bi,bj)*(
                0133      &       vVel( i ,j,k,bi,bj)*dyC( i ,j,bi,bj)
                0134      &      -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
                0135      &      -uVel(i, j ,k,bi,bj)*dxC(i, j ,bi,bj)
                0136      &      +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
4d2b0c1389 Jean*0137      &                                       )
0a272ce360 Jean*0138           ENDDO
                0139          ENDDO
                0140 
                0141 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0142 C     Special stuff for Cubed Sphere:
                0143          IF (useCubedSphereExchange) THEN
                0144 c          AZcorner = 0.75 _d 0
                0145            iMax = sNx+1
                0146            jMax = sNy+1
                0147            DO i=1,iMax
                0148             hFacZ(i,jMax)=0.
                0149             vort3(i,jMax)=0.
                0150            ENDDO
                0151            DO j=1,jMax
                0152             hFacZ(iMax,j)=0.
                0153             vort3(iMax,j)=0.
                0154            ENDDO
e98f4d9780 Jean*0155 
                0156            southWestCorner = .TRUE.
                0157            southEastCorner = .TRUE.
                0158            northWestCorner = .TRUE.
                0159            northEastCorner = .TRUE.
                0160            iG = bi+(myXGlobalLo-1)/sNx
4749c74143 Alis*0161 #ifdef ALLOW_EXCH2
c424ee7cc7 Jean*0162            myTile = W2_myTileList(bi,bj)
e98f4d9780 Jean*0163            iG = exch2_myFace(myTile)
4d2b0c1389 Jean*0164            southWestCorner = exch2_isWedge(myTile).EQ.1
e98f4d9780 Jean*0165      &                 .AND. exch2_isSedge(myTile).EQ.1
                0166            southEastCorner = exch2_isEedge(myTile).EQ.1
                0167      &                 .AND. exch2_isSedge(myTile).EQ.1
                0168            northEastCorner = exch2_isEedge(myTile).EQ.1
                0169      &                 .AND. exch2_isNedge(myTile).EQ.1
                0170            northWestCorner = exch2_isWedge(myTile).EQ.1
                0171      &                 .AND. exch2_isNedge(myTile).EQ.1
4749c74143 Alis*0172 #endif /* ALLOW_EXCH2 */
e98f4d9780 Jean*0173 
                0174 C--        avoid to count 3 times the same corner:
                0175            southEastCorner = southEastCorner .AND. iG.EQ.2
                0176            northWestCorner = northWestCorner .AND. iG.EQ.1
                0177            northEastCorner = .FALSE.
                0178 
404487e0d3 Andr*0179 C--       S.W. corner:
                0180           IF ( southWestCorner ) THEN
0a272ce360 Jean*0181            i=1
                0182            j=1
                0183            vort3(i,j)=
                0184      &       +recip_rAz(i,j,bi,bj)/AZcorner*(
                0185      &        vVel(i,j,k,bi,bj)*dyC(i,j,bi,bj)
                0186      &       -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
                0187      &       +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
                0188      &       )
4d2b0c1389 Jean*0189            hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
                0190      &                  + _hFacW(i, j ,k,bi,bj)
0a272ce360 Jean*0191      &                  + _hFacS( i ,j,k,bi,bj)
                0192      &                  )/3. _d 0
404487e0d3 Andr*0193           ENDIF
                0194           IF ( southEastCorner ) THEN
0a272ce360 Jean*0195 C--        S.E. corner:
                0196            i=iMax
                0197            j=1
4d2b0c1389 Jean*0198            vort3(i,j)=
                0199      &       +recip_rAz(i,j,bi,bj)/AZcorner*(
0a272ce360 Jean*0200      &       -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
                0201      &       -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
                0202      &       +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
                0203      &       )
4d2b0c1389 Jean*0204            hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
                0205      &                  + _hFacW(i, j ,k,bi,bj)
                0206      &                  + _hFacS(i-1,j,k,bi,bj)
0a272ce360 Jean*0207      &                  )/3. _d 0
404487e0d3 Andr*0208           ENDIF
                0209           IF ( northWestCorner ) THEN
0a272ce360 Jean*0210 C--        N.W. corner:
                0211            i=1
                0212            j=jMax
                0213            vort3(i,j)=
                0214      &       +recip_rAz(i,j,bi,bj)/AZcorner*(
                0215      &        vVel(i,j,k,bi,bj)*dyC(i,j,bi,bj)
                0216      &       -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
                0217      &       +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
                0218      &       )
4d2b0c1389 Jean*0219            hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
                0220      &                  + _hFacW(i, j ,k,bi,bj)
0a272ce360 Jean*0221      &                  + _hFacS( i ,j,k,bi,bj)
                0222      &                  )/3. _d 0
404487e0d3 Andr*0223           ENDIF
                0224           IF ( northEastCorner ) THEN
0a272ce360 Jean*0225 C--        N.E. corner:
                0226            i=iMax
                0227            j=jMax
                0228            vort3(i,j)=
                0229      &       +recip_rAz(i,j,bi,bj)/AZcorner*(
                0230      &       -vVel(i-1,j,k,bi,bj)*dyC(i-1,j,bi,bj)
                0231      &       -uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
                0232      &       +uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
                0233      &       )
4d2b0c1389 Jean*0234            hFacZ(i,j) = ( _hFacW(i,j-1,k,bi,bj)
                0235      &                  + _hFacW(i, j ,k,bi,bj)
                0236      &                  + _hFacS(i-1,j,k,bi,bj)
0a272ce360 Jean*0237      &                  )/3. _d 0
404487e0d3 Andr*0238           ENDIF
0a272ce360 Jean*0239          ENDIF
                0240 
                0241 C-    Special stuff for North & South Poles, LatLon grid
                0242          IF ( usingSphericalPolarGrid ) THEN
                0243           IF (yG(1,sNy+1,bi,bj).EQ.90.) THEN
                0244            jMax = sNy+1
                0245            j = jMax
                0246            DO i=1,sNx
                0247             vort3(i,j) = 0.
                0248             vort3(1,j) = vort3(1,j)
                0249      &                 + uVel(i,j-1,k,bi,bj)*dxC(i,j-1,bi,bj)
                0250             hFacZ(i,j) = 0.
                0251 #ifndef MONITOR_TEST_HFACZ
4d2b0c1389 Jean*0252             hFacZ(1,j) = hFacZ(1,j) + _hFacW(i,j-1,k,bi,bj)
0a272ce360 Jean*0253            ENDDO
                0254 #else
4d2b0c1389 Jean*0255             hFacZ(1,j) = hFacZ(1,j) + etaFld(i,j-1)
0a272ce360 Jean*0256            ENDDO
                0257             hFacZ(1,j) = sNx + hFacZ(1,j)*recip_drF(k)
                0258 #endif
                0259             hFacZ(1,j) = hFacZ(1,j) / FLOAT(sNx)
                0260             vort3(1,j) = vort3(1,j)*recip_rAz(1,j,bi,bj)
                0261           ENDIF
                0262           IF (yG(1,1,bi,bj).EQ.-90.) THEN
                0263            j = 1
                0264            DO i=1,sNx
                0265             vort3(i,j) = 0.
                0266             vort3(1,j) = vort3(1,j)
                0267      &                 - uVel(i,j,k,bi,bj)*dxC(i,j,bi,bj)
                0268             hFacZ(i,j) = 0.
                0269 #ifndef MONITOR_TEST_HFACZ
4d2b0c1389 Jean*0270             hFacZ(1,j) = hFacZ(1,j) + _hFacW(i,j,k,bi,bj)
0a272ce360 Jean*0271            ENDDO
                0272 #else
4d2b0c1389 Jean*0273             hFacZ(1,j) = hFacZ(1,j) + etaFld(i,j)
0a272ce360 Jean*0274            ENDDO
                0275             hFacZ(1,j) = sNx + hFacZ(1,j)*recip_drF(k)
                0276 #endif
                0277             hFacZ(1,j) = hFacZ(1,j) / FLOAT(sNx)
                0278             vort3(1,j) = vort3(1,j)*recip_rAz(1,j,bi,bj)
                0279           ENDIF
                0280          ENDIF
                0281 
                0282 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0283 
4d2b0c1389 Jean*0284          DO j=1,jMax
                0285           DO i=1,iMax
0a272ce360 Jean*0286            IF (hFacZ(i,j).GT.0. _d 0) THEN
                0287             tmpVal = vort3(i,j)
                0288             tmpAre = rAz(i,j,bi,bj)*drF(k)
                0289             tmpVol = rAz(i,j,bi,bj)*drF(k)*hFacZ(i,j)
4d2b0c1389 Jean*0290             tileArea(bi,bj)  = tileArea(bi,bj)  + tmpAre
0a272ce360 Jean*0291 C-       min,max of relative vorticity ("r")
                0292             theMin = MIN(theMin,tmpVal)
                0293             theMax = MAX(theMax,tmpVal)
                0294 C-       average & std.dev of absolute vorticity ("a")
                0295             tmpVal = tmpVal + fCoriG(i,j,bi,bj)
4d2b0c1389 Jean*0296             tileSum(bi,bj) = tileSum(bi,bj) + tmpAre*tmpVal
                0297             tileVar(bi,bj) = tileVar(bi,bj) + tmpAre*tmpVal*tmpVal
0a272ce360 Jean*0298 C-       average & std.dev of potential vorticity ("p")
                0299             tmpVal = tmpVal / hFacZ(i,j)
4d2b0c1389 Jean*0300             tileVol(bi,bj) = tileVol(bi,bj) + tmpVol
                0301             tileVSum(bi,bj)= tileVSum(bi,bj)+ tmpVol*tmpVal
                0302             tileVSq(bi,bj) = tileVSq(bi,bj) + tmpVol*tmpVal*tmpVal
0a272ce360 Jean*0303            ENDIF
                0304           ENDDO
                0305          ENDDO
                0306 
                0307         ENDDO
4d2b0c1389 Jean*0308 c       theArea= theArea + tileArea(bi,bj)
                0309 c       theVol = theVol  + tileVol (bi,bj)
                0310 c       theMean= theMean + tileSum(bi,bj)
                0311 c       theVar = theVar  + tileVar(bi,bj)
                0312 c       volMean= volMean + tileVSum(bi,bj)
                0313 c       volVar = volVar  + tileVSq(bi,bj)
0a272ce360 Jean*0314        ENDDO
                0315       ENDDO
                0316 
                0317       theMin = -theMin
7163a40534 Jean*0318       _GLOBAL_MAX_RL(theMin, myThid)
                0319       _GLOBAL_MAX_RL(theMax, myThid)
0a272ce360 Jean*0320       theMin = -theMin
7163a40534 Jean*0321 c     _GLOBAL_SUM_RL(theArea,myThid)
                0322 c     _GLOBAL_SUM_RL(theVol, myThid)
                0323 c     _GLOBAL_SUM_RL(theMean,myThid)
                0324 c     _GLOBAL_SUM_RL(theVar, myThid)
                0325 c     _GLOBAL_SUM_RL(volMean,myThid)
                0326 c     _GLOBAL_SUM_RL(volVar ,myThid)
4d2b0c1389 Jean*0327       CALL GLOBAL_SUM_TILE_RL( tileArea, theArea, myThid )
                0328       CALL GLOBAL_SUM_TILE_RL( tileVol,  theVol,  myThid )
                0329       CALL GLOBAL_SUM_TILE_RL( tileSum,  theMean, myThid )
                0330       CALL GLOBAL_SUM_TILE_RL( tileVar,  theVar,  myThid )
                0331       CALL GLOBAL_SUM_TILE_RL( tileVSum, volMean, myThid )
                0332       CALL GLOBAL_SUM_TILE_RL( tileVSq,  volVar,  myThid )
                0333       IF (theArea.GT.0.) THEN
0a272ce360 Jean*0334         theMean= theMean/theArea
                0335         theVar = theVar /theArea
                0336         theVar = theVar - theMean*theMean
                0337 c       IF (theVar.GT.0.) theSD = SQRT(theVar)
                0338         IF (theVar.GT.0.) theVar = SQRT(theVar)
                0339       ENDIF
                0340       IF (theVol.GT.0.) THEN
                0341         volMean= volMean/theVol
                0342         volVar = volVar /theVol
                0343         volVar = volVar - volMean*volMean
                0344         IF (volVar.GT.0.) theSD = SQRT(volVar)
                0345       ENDIF
                0346 
                0347 C-    Print stats for (relative/absolute) Vorticity AND Pot.Vort.
                0348       CALL MON_SET_PREF('vort',myThid)
                0349       CALL MON_OUT_RL(mon_string_none,theMin, '_r_min',  myThid)
                0350       CALL MON_OUT_RL(mon_string_none,theMax, '_r_max',  myThid)
                0351       CALL MON_OUT_RL(mon_string_none,theMean,'_a_mean', myThid)
                0352       CALL MON_OUT_RL(mon_string_none,theVar, '_a_sd',   myThid)
                0353       CALL MON_OUT_RL(mon_string_none,volMean,'_p_mean', myThid)
                0354       CALL MON_OUT_RL(mon_string_none,theSD,  '_p_sd',   myThid)
                0355 c     CALL MON_OUT_RL(mon_string_none,theVol,mon_foot_vol,myThid)
                0356 
                0357       RETURN
                0358       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!