C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/diagnostics/diagstats_global.F,v 1.1 2005/05/20 07:28:52 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, 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 diagnostics long list)
C     mate   ..... counter diagnostic Id number if any ; 0 otherwise
C     myThid ..... my thread Id number
      _RL undef
      INTEGER nLev, jReg, ndId, mate
      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

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 = iSdiag(ndId) + 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 = iSdiag(mate) + k - 1
            IF ( qSdiag(0,jReg,kd,bi,bj).GT.0. ) THEN
             IF ( qtmp1(0,k).LE.0. ) THEN
              DO i=0,nStats
                qtmp2(i,k) = qSdiag(i,jReg,kd,bi,bj)
              ENDDO
             ELSE
              DO i=0,iv
                qtmp2(i,k) = qtmp2(i,k) + qSdiag(i,jReg,kd,bi,bj)
              ENDDO
              qtmp2(im,k) = MIN( qtmp2(im,k),qSdiag(im,jReg,kd,bi,bj) )
              qtmp2(ix,k) = MAX( qtmp2(ix,k),qSdiag(ix,jReg,kd,bi,bj) )
             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
           DO i=0,iv
            _GLOBAL_SUM_R8(qtmp1(i,k),myThid)
           ENDDO
            qtmp1(im,k) = -qtmp1(im,k)
            _GLOBAL_MAX_R8(qtmp1(im,k),myThid)
            qtmp1(im,k) = -qtmp1(im,k)
            _GLOBAL_MAX_R8(qtmp1(ix,k),myThid)
        ENDDO
        IF ( mate.GT.0 ) THEN
         DO k=1,kMlev
           DO i=0,iv
            _GLOBAL_SUM_R8(qtmp2(i,k),myThid)
           ENDDO
            qtmp2(im,k) = -qtmp2(im,k)
            _GLOBAL_MAX_R8(qtmp2(im,k),myThid)
            qtmp2(im,k) = -qtmp2(im,k)
            _GLOBAL_MAX_R8(qtmp2(ix,k),myThid)
         ENDDO
        ENDIF

C--     Vertical integral, min & max :
        DO k=1,klev
           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
        ENDDO
        IF ( mate.GT.0 ) THEN
         DO k=1,kMlev
           IF ( qtmp2(0,0).LE.0. ) THEN
             DO i=0,nStats
              qtmp2(i,0) = qtmp2(i,k)
             ENDDO
           ELSE
             DO i=0,iv
              qtmp2(i,0) = qtmp2(i,0) + qtmp2(i,k)
             ENDDO
              qtmp2(im,0) = MIN(qtmp2(im,0),qtmp2(im,k))
              qtmp2(ix,0) = MAX(qtmp2(ix,0),qtmp2(ix,k))
           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
        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