C $Header: /u/gcmpack/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