pkg/monitor/mon_calc_stats_rs.F

0d404f33b0 Jean*0001 #include "MONITOR_OPTIONS.h"
                0002
                0003 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0004 CBOP
                0005 C     !ROUTINE: MON_CALC_STATS_RS
                0006
                0007 C     !INTERFACE:
                0008       SUBROUTINE MON_CALC_STATS_RS(
                0009      I               myNr, arr, arrhFac, arrMask, arrArea, arrDr,
                0010      O               theMin, theMax, theMean, theSD, theDel2, theVol,
                0011      I               myThid )
                0012
                0013 C     Calculate statistics of global array ``\_RS arr''.
                0014 C     account for volume and mask
                0015
                0016 C     !USES:
                0017       IMPLICIT NONE
                0018 #include "SIZE.h"
                0019 #include "EEPARAMS.h"
                0020
                0021 C     !INPUT/OUTPUT PARAMETERS:
                0022       INTEGER myNr
                0023       _RS arr    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,myNr,nSx,nSy)
                0024       _RS arrhFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy,myNr,nSx,nSy)
                0025       _RS arrMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0026       _RS arrArea(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0027       _RS arrDr(myNr)
                0028       _RL theMin, theMax, theMean, theSD, theDel2, theVol
                0029       INTEGER myThid
                0030 CEOP
                0031
                0032 C     !LOCAL VARIABLES:
                0033       INTEGER bi,bj,i,j,k
                0034       LOGICAL noPnts
                0035       _RL tmpVal
                0036       _RL tmpMask
                0037       _RL tmpVol
                0038       _RL ddx, ddy
                0039       _RL theVar
d3ae19070d Jean*0040       _RL theNbPt
0d404f33b0 Jean*0041       _RL tileMean(nSx,nSy)
                0042       _RL tileVar (nSx,nSy)
                0043       _RL tileSD  (nSx,nSy)
                0044       _RL tileDel2(nSx,nSy)
                0045       _RL tileVol (nSx,nSy)
d3ae19070d Jean*0046       _RL tileNbPt(nSx,nSy)
0d404f33b0 Jean*0047
                0048       theMin = 0.
                0049       theMax = 0.
                0050       theMean= 0.
                0051       theSD  = 0.
                0052       theVar = 0.
                0053       theDel2= 0.
                0054       theVol = 0.
d3ae19070d Jean*0055       theNbPt= 0.
0d404f33b0 Jean*0056       noPnts = .TRUE.
                0057
                0058       DO bj=myByLo(myThid),myByHi(myThid)
                0059        DO bi=myBxLo(myThid),myBxHi(myThid)
d3ae19070d Jean*0060         tileNbPt(bi,bj) = 0.
0d404f33b0 Jean*0061         tileDel2(bi,bj) = 0.
                0062         tileVol (bi,bj) = 0.
                0063         tileMean(bi,bj) = 0.
                0064         tileVar (bi,bj) = 0.
                0065         DO k=1,myNr
                0066          DO j=1,sNy
                0067           DO i=1,sNx
                0068            tmpVal  = arr(i,j,k,bi,bj)
                0069            tmpMask = arrMask(i,j,bi,bj)*arrhFac(i,j,k,bi,bj)
                0070            IF ( tmpMask.GT.0. _d 0 .AND. noPnts ) THEN
                0071             theMin=tmpVal
                0072             theMax=tmpVal
                0073             noPnts=.FALSE.
                0074            ENDIF
                0075            IF ( tmpMask.GT.0. _d 0 ) THEN
                0076             theMin = MIN(theMin,tmpVal)
                0077             theMax = MAX(theMax,tmpVal)
                0078 C--   like old code (but using hFac instead of mask): identical if no partial cell
                0079 c           tileDel2(bi,bj) = tileDel2(bi,bj)
                0080 c    &       + 0.25*ABS(
                0081 c    &          (arr(i+1,j,k,bi,bj)-tmpVal)*arrhFac(i+1,j,k,bi,bj)
                0082 c    &         +(arr(i-1,j,k,bi,bj)-tmpVal)*arrhFac(i-1,j,k,bi,bj)
                0083 c    &         +(arr(i,j+1,k,bi,bj)-tmpVal)*arrhFac(i,j+1,k,bi,bj)
                0084 c    &         +(arr(i,j-1,k,bi,bj)-tmpVal)*arrhFac(i,j-1,k,bi,bj)
                0085 c    &                 )
                0086 C--   New form:
                0087             ddx = arrhFac(i+1,j,k,bi,bj)*arrhFac(i-1,j,k,bi,bj)
                0088             IF ( ddx.GT.0. _d 0 ) THEN
                0089              ddx = (arr(i+1,j,k,bi,bj)-tmpVal)
                0090      &           + (arr(i-1,j,k,bi,bj)-tmpVal)
                0091             ENDIF
                0092             ddy = arrhFac(i,j+1,k,bi,bj)*arrhFac(i,j-1,k,bi,bj)
                0093             IF ( ddy.GT.0. _d 0 ) THEN
                0094              ddy = (arr(i,j+1,k,bi,bj)-tmpVal)
                0095      &           + (arr(i,j-1,k,bi,bj)-tmpVal)
                0096             ENDIF
                0097             tileDel2(bi,bj) = tileDel2(bi,bj) + ddx*ddx + ddy*ddy
                0098
d3ae19070d Jean*0099             tileNbPt(bi,bj) = tileNbPt(bi,bj) + oneRL
0d404f33b0 Jean*0100             tmpVol = arrArea(i,j,bi,bj)*arrDr(k)*tmpMask
                0101             tileVol (bi,bj) = tileVol (bi,bj) + tmpVol
                0102             tileMean(bi,bj) = tileMean(bi,bj) + tmpVol*tmpVal
                0103             tileVar (bi,bj) = tileVar (bi,bj) + tmpVol*tmpVal*tmpVal
                0104            ENDIF
                0105           ENDDO
                0106          ENDDO
                0107         ENDDO
                0108        ENDDO
                0109       ENDDO
                0110
d3ae19070d Jean*0111       CALL GLOBAL_SUM_TILE_RL( tileNbPt, theNbPt, myThid )
0d404f33b0 Jean*0112       CALL GLOBAL_SUM_TILE_RL( tileDel2, theDel2, myThid )
                0113       CALL GLOBAL_SUM_TILE_RL( tileVol , theVol , myThid )
                0114       CALL GLOBAL_SUM_TILE_RL( tileMean, theMean, myThid )
                0115 c     CALL GLOBAL_SUM_TILE_RL( tileVar , theVar , myThid )
                0116
d3ae19070d Jean*0117       IF ( theNbPt.GT.zeroRL ) THEN
                0118 c      theDel2  = theDel2/theNbPt
                0119        theDel2  = SQRT(theDel2)/theNbPt
0d404f33b0 Jean*0120       ENDIF
                0121
                0122       IF ( theVol.GT.0. _d 0 ) THEN
                0123        theMean= theMean/theVol
                0124        theVar = theVar/theVol
                0125        IF ( noPnts ) theMin = theMean
                0126        theMin = -theMin
                0127        _GLOBAL_MAX_RL(theMin,myThid)
                0128        theMin = -theMin
                0129        IF ( noPnts ) theMax = theMean
                0130        _GLOBAL_MAX_RL(theMax,myThid)
                0131
                0132        DO bj=myByLo(myThid),myByHi(myThid)
                0133         DO bi=myBxLo(myThid),myBxHi(myThid)
                0134          tileSD(bi,bj)=0.
                0135          DO k=1,myNr
                0136           DO j=1,sNy
                0137            DO i=1,sNx
                0138             tmpVal=arr(i,j,k,bi,bj)
                0139             tmpMask = arrMask(i,j,bi,bj)*arrhFac(i,j,k,bi,bj)
                0140             IF ( tmpMask.GT.0. _d 0 ) THEN
                0141              tmpVol = arrArea(i,j,bi,bj)*arrDr(k)*tmpMask
                0142              tileSD(bi,bj) = tileSD(bi,bj)
                0143      &                     + tmpVol*(tmpVal-theMean)*(tmpVal-theMean)
                0144             ENDIF
                0145            ENDDO
                0146           ENDDO
                0147          ENDDO
                0148         ENDDO
                0149        ENDDO
                0150
                0151        CALL GLOBAL_SUM_TILE_RL( tileSD, theSD, myThid )
                0152
                0153        theSD = SQRT(theSD/theVol)
                0154 c      theSD = SQRT(theVar-theMean*theMean)
                0155       ENDIF
                0156
                0157       RETURN
                0158       END