C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/diagnostics/diagstats_global.F,v 1.5 2009/04/28 18:10:47 jmc Exp $
C $Name:  $

#include "DIAG_OPTIONS.h"

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

C     !INTERFACE:
      SUBROUTINE DIAGSTATS_GLOBAL(
     O                    qtmp1, qtmp2,
     I                    undef, nLev, jReg,
     I                    ndId, mate, iSp, iSm, myThid )

C     !DESCRIPTION:
C     Retrieve averaged model diagnostic

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

C     !INPUT PARAMETERS:
C     undef     :: Undefined value
C     nLev      :: 2nd Dimension (max Nb of levels) of qtmp1,2 arrays
C     jReg      :: region Index to be process.
C     ndId      :: diagnostic Id number (in available diagnostics list)
C     mate      :: counter mate Id number if any ; 0 otherwise
C     iSp       :: diagnostics  pointer to storage array
C     iSm       :: counter-mate pointer to storage array
C     myThid    :: my thread Id number
      _RL undef
      INTEGER nLev, jReg, ndId, mate, iSp, iSm
      INTEGER myThid

C     !OUTPUT PARAMETERS:
C     qtmp1    ..... AVERAGED DIAGNOSTIC QUANTITY
C     qtmp2    ..... working array (used for counter mate statistics)
      _RL qtmp1(0:nStats,0:nLev)
      _RL qtmp2(0:nStats,0:nLev)
CEOP

C     !LOCAL VARIABLES:
      INTEGER im, ix, iv
      PARAMETER ( iv = nStats - 2 , im = nStats - 1 , ix = nStats )
      INTEGER bi, bj
      INTEGER i, k, kd, kCnt, klev, kMlev
      _RL     tmpMin, tmpMax, tmpVol

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

C--   Initialize to zero :
      DO k=0,nLev
        DO i=0,nStats
          qtmp1(i,k) = 0.
          qtmp2(i,k) = 0.
        ENDDO
      ENDDO

      klev = kdiag(ndId)
      IF ( mate.GT.0 ) kMlev = kdiag(mate)

      IF (klev.LE.nLev) THEN
C---    Compute global statistics :

C--     Retrieve tile statistics first
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)

          DO k=1,klev
            kd = iSp + k - 1
            IF ( qSdiag(0,jReg,kd,bi,bj).GT.0. ) THEN
             IF ( qtmp1(0,k).LE.0. ) THEN
              DO i=0,nStats
                qtmp1(i,k) = qSdiag(i,jReg,kd,bi,bj)
              ENDDO
             ELSE
              DO i=0,iv
                qtmp1(i,k) = qtmp1(i,k) + qSdiag(i,jReg,kd,bi,bj)
              ENDDO
              qtmp1(im,k) = MIN( qtmp1(im,k),qSdiag(im,jReg,kd,bi,bj) )
              qtmp1(ix,k) = MAX( qtmp1(ix,k),qSdiag(ix,jReg,kd,bi,bj) )
             ENDIF
            ENDIF
          ENDDO
          IF ( mate.GT.0 ) THEN
           DO k=1,kMlev
            kd = iSm + k - 1
            IF ( qSdiag(0,jReg,kd,bi,bj).GT.0. ) THEN
             IF ( qtmp2(0,k).LE.0. ) THEN
              DO i=0,1
                qtmp2(i,k) = qSdiag(i,jReg,kd,bi,bj)
              ENDDO
             ELSE
              DO i=0,1
                qtmp2(i,k) = qtmp2(i,k) + qSdiag(i,jReg,kd,bi,bj)
              ENDDO
             ENDIF
            ENDIF
           ENDDO
          ENDIF

C-       end tile index loops
         ENDDO
        ENDDO

C--     Global min,max & sum (at each level) over all thread & processors :
        DO k=1,klev
           tmpVol = qtmp1(0,k)
           DO i=0,iv
            _GLOBAL_SUM_RL(qtmp1(i,k),myThid)
           ENDDO
           IF ( qtmp1(0,k).GT.0. .AND. tmpVol.LE.0. ) THEN
C-      In case 1 processor has only empty tiles:
             tmpMax = qtmp1(1,k)/qtmp1(0,k)
             tmpmin = -tmpMax
           ELSE
             tmpMin = -qtmp1(im,k)
             tmpMax =  qtmp1(ix,k)
           ENDIF
           _GLOBAL_MAX_RL(tmpMin,myThid)
           _GLOBAL_MAX_RL(tmpMax,myThid)
           qtmp1(im,k) = -tmpMin
           qtmp1(ix,k) =  tmpMax
        ENDDO
        IF ( mate.GT.0 ) THEN
         DO k=1,kMlev
           DO i=0,1
            _GLOBAL_SUM_RL(qtmp2(i,k),myThid)
           ENDDO
         ENDDO
        ENDIF

C--     Vertical integral, min & max :
        DO k=1,klev
          IF ( qtmp1(0,k).GT.0. ) THEN
           IF ( qtmp1(0,0).LE.0. ) THEN
             DO i=0,nStats
              qtmp1(i,0) = qtmp1(i,k)
             ENDDO
           ELSE
             DO i=0,iv
              qtmp1(i,0) = qtmp1(i,0) + qtmp1(i,k)
             ENDDO
              qtmp1(im,0) = MIN(qtmp1(im,0),qtmp1(im,k))
              qtmp1(ix,0) = MAX(qtmp1(ix,0),qtmp1(ix,k))
           ENDIF
          ENDIF
        ENDDO
        IF ( mate.GT.0 ) THEN
         DO k=1,kMlev
          IF ( qtmp2(0,k).GT.0. ) THEN
           IF ( qtmp2(0,0).LE.0. ) THEN
             DO i=0,1
              qtmp2(i,0) = qtmp2(i,k)
             ENDDO
           ELSE
             DO i=0,1
              qtmp2(i,0) = qtmp2(i,0) + qtmp2(i,k)
             ENDDO
           ENDIF
          ENDIF
         ENDDO
        ENDIF

C--     Average, Standard.Dev.:
C-      no counter diagnostics => average = Sum / vol :
        IF ( mate.EQ.0 ) THEN
          DO k=0,klev
           IF ( qtmp1(0,k).LE.0. ) THEN
             DO i=1,nStats
              qtmp1(i,k) = undef
             ENDDO
           ELSE
             DO i=1,iv
              qtmp1(i,k) = qtmp1(i,k) / qtmp1(0,k)
             ENDDO
C            Variance :
             qtmp1(iv,k) = qtmp1(iv,k) - qtmp1(1,k)*qtmp1(1,k)
C            Standard deviation :
             IF (qtmp1(iv,k).GT.0.) qtmp1(iv,k) = SQRT(qtmp1(iv,k))
           ENDIF
          ENDDO
C         return global (& vertically integrated) volume in qtmp2(0,0):
          qtmp2(0,0) = qtmp1(0,0)
        ELSE
C       With counter diagnostics => average = Sum / Sum(counter) :
          DO k=0,klev
           kCnt = min(k,kMlev)
           IF ( qtmp2(0,kCnt).LE.0. ) THEN
             DO i=1,nStats
              qtmp1(i,k) = undef
             ENDDO
           ELSEIF ( qtmp2(1,kCnt).LE.0. ) THEN
             DO i=1,iv
              qtmp1(i,k) = undef
             ENDDO
           ELSE
             DO i=1,iv
              qtmp1(i,k) = qtmp1(i,k) / qtmp2(1,kCnt)
             ENDDO
C jmc: looks like there is a Pb with how Variance is computed
C            Variance :
             qtmp1(iv,k) = qtmp1(iv,k) - qtmp1(1,k)*qtmp1(1,k)
C            Standard deviation :
             IF (qtmp1(iv,k).GT.0.) qtmp1(iv,k) = SQRT(qtmp1(iv,k))
           ENDIF
          ENDDO
        ENDIF

      ENDIF

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

      RETURN
      END