pkg/diagnostics/diagstats_global.F

C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagstats_global.F,v 1.3 2005/07/10 00:52:12 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_R8(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_R8(tmpMin,myThid)
           _GLOBAL_MAX_R8(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_R8(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
1	jmc	1.4	C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagstats_global.F,v 1.3 2005/07/10 00:52:12 jmc Exp $
2	jmc	1.1	C $Name: $
3
4			#include "DIAG_OPTIONS.h"
5
6			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
7			CBOP 0
8			C !ROUTINE: DIAGSTATS_GLOBAL
9
10			C !INTERFACE:
11			SUBROUTINE DIAGSTATS_GLOBAL(
12			O qtmp1, qtmp2,
13			I undef, nLev, jReg,
14	jmc	1.2	I ndId, mate, iSp, iSm, myThid )
15	jmc	1.1
16			C !DESCRIPTION:
17			C Retrieve averaged model diagnostic
18
19			C !USES:
20			IMPLICIT NONE
21			#include "EEPARAMS.h"
22			#include "SIZE.h"
23			#include "DIAGNOSTICS_SIZE.h"
24			#include "DIAGNOSTICS.h"
25
26			C !INPUT PARAMETERS:
27	jmc	1.2	C undef :: Undefined value
28			C nLev :: 2nd Dimension (max Nb of levels) of qtmp1,2 arrays
29			C jReg :: region Index to be process.
30			C ndId :: diagnostic Id number (in available diagnostics list)
31			C mate :: counter mate Id number if any ; 0 otherwise
32			C iSp :: diagnostics pointer to storage array
33			C iSm :: counter-mate pointer to storage array
34			C myThid :: my thread Id number
35	jmc	1.1	_RL undef
36	jmc	1.2	INTEGER nLev, jReg, ndId, mate, iSp, iSm
37	jmc	1.1	INTEGER myThid
38
39			C !OUTPUT PARAMETERS:
40			C qtmp1 ..... AVERAGED DIAGNOSTIC QUANTITY
41			C qtmp2 ..... working array (used for counter mate statistics)
42			_RL qtmp1(0:nStats,0:nLev)
43			_RL qtmp2(0:nStats,0:nLev)
44			CEOP
45
46			C !LOCAL VARIABLES:
47			INTEGER im, ix, iv
48			PARAMETER ( iv = nStats - 2 , im = nStats - 1 , ix = nStats )
49			INTEGER bi, bj
50			INTEGER i, k, kd, kCnt, klev, kMlev
51	jmc	1.4	_RL tmpMin, tmpMax, tmpVol
52	jmc	1.1
53			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
54
55			C-- Initialize to zero :
56			DO k=0,nLev
57			DO i=0,nStats
58			qtmp1(i,k) = 0.
59			qtmp2(i,k) = 0.
60			ENDDO
61			ENDDO
62
63			klev = kdiag(ndId)
64			IF ( mate.GT.0 ) kMlev = kdiag(mate)
65
66			IF (klev.LE.nLev) THEN
67			C--- Compute global statistics :
68
69			C-- Retrieve tile statistics first
70			DO bj=myByLo(myThid),myByHi(myThid)
71			DO bi=myBxLo(myThid),myBxHi(myThid)
72
73			DO k=1,klev
74	jmc	1.2	kd = iSp + k - 1
75	jmc	1.1	IF ( qSdiag(0,jReg,kd,bi,bj).GT.0. ) THEN
76			IF ( qtmp1(0,k).LE.0. ) THEN
77			DO i=0,nStats
78			qtmp1(i,k) = qSdiag(i,jReg,kd,bi,bj)
79			ENDDO
80			ELSE
81			DO i=0,iv
82			qtmp1(i,k) = qtmp1(i,k) + qSdiag(i,jReg,kd,bi,bj)
83			ENDDO
84			qtmp1(im,k) = MIN( qtmp1(im,k),qSdiag(im,jReg,kd,bi,bj) )
85			qtmp1(ix,k) = MAX( qtmp1(ix,k),qSdiag(ix,jReg,kd,bi,bj) )
86			ENDIF
87			ENDIF
88			ENDDO
89			IF ( mate.GT.0 ) THEN
90			DO k=1,kMlev
91	jmc	1.2	kd = iSm + k - 1
92	jmc	1.1	IF ( qSdiag(0,jReg,kd,bi,bj).GT.0. ) THEN
93	jmc	1.3	IF ( qtmp2(0,k).LE.0. ) THEN
94	jmc	1.4	DO i=0,1
95	jmc	1.1	qtmp2(i,k) = qSdiag(i,jReg,kd,bi,bj)
96			ENDDO
97			ELSE
98	jmc	1.4	DO i=0,1
99	jmc	1.1	qtmp2(i,k) = qtmp2(i,k) + qSdiag(i,jReg,kd,bi,bj)
100			ENDDO
101			ENDIF
102			ENDIF
103			ENDDO
104			ENDIF
105
106			C- end tile index loops
107			ENDDO
108			ENDDO
109
110			C-- Global min,max & sum (at each level) over all thread & processors :
111			DO k=1,klev
112	jmc	1.4	tmpVol = qtmp1(0,k)
113	jmc	1.1	DO i=0,iv
114			_GLOBAL_SUM_R8(qtmp1(i,k),myThid)
115			ENDDO
116	jmc	1.4	IF ( qtmp1(0,k).GT.0. .AND. tmpVol.LE.0. ) THEN
117			C- In case 1 processor has only empty tiles:
118			tmpMax = qtmp1(1,k)/qtmp1(0,k)
119			tmpmin = -tmpMax
120			ELSE
121			tmpMin = -qtmp1(im,k)
122			tmpMax = qtmp1(ix,k)
123			ENDIF
124			_GLOBAL_MAX_R8(tmpMin,myThid)
125			_GLOBAL_MAX_R8(tmpMax,myThid)
126			qtmp1(im,k) = -tmpMin
127			qtmp1(ix,k) = tmpMax
128	jmc	1.1	ENDDO
129			IF ( mate.GT.0 ) THEN
130			DO k=1,kMlev
131	jmc	1.4	DO i=0,1
132	jmc	1.1	_GLOBAL_SUM_R8(qtmp2(i,k),myThid)
133			ENDDO
134			ENDDO
135			ENDIF
136
137			C-- Vertical integral, min & max :
138			DO k=1,klev
139	jmc	1.4	IF ( qtmp1(0,k).GT.0. ) THEN
140	jmc	1.1	IF ( qtmp1(0,0).LE.0. ) THEN
141			DO i=0,nStats
142			qtmp1(i,0) = qtmp1(i,k)
143			ENDDO
144			ELSE
145			DO i=0,iv
146			qtmp1(i,0) = qtmp1(i,0) + qtmp1(i,k)
147			ENDDO
148			qtmp1(im,0) = MIN(qtmp1(im,0),qtmp1(im,k))
149			qtmp1(ix,0) = MAX(qtmp1(ix,0),qtmp1(ix,k))
150			ENDIF
151	jmc	1.4	ENDIF
152	jmc	1.1	ENDDO
153			IF ( mate.GT.0 ) THEN
154			DO k=1,kMlev
155	jmc	1.4	IF ( qtmp2(0,k).GT.0. ) THEN
156	jmc	1.1	IF ( qtmp2(0,0).LE.0. ) THEN
157	jmc	1.4	DO i=0,1
158	jmc	1.1	qtmp2(i,0) = qtmp2(i,k)
159			ENDDO
160			ELSE
161	jmc	1.4	DO i=0,1
162	jmc	1.1	qtmp2(i,0) = qtmp2(i,0) + qtmp2(i,k)
163			ENDDO
164			ENDIF
165	jmc	1.4	ENDIF
166	jmc	1.1	ENDDO
167			ENDIF
168
169			C-- Average, Standard.Dev.:
170			C- no counter diagnostics => average = Sum / vol :
171			IF ( mate.EQ.0 ) THEN
172			DO k=0,klev
173			IF ( qtmp1(0,k).LE.0. ) THEN
174			DO i=1,nStats
175			qtmp1(i,k) = undef
176			ENDDO
177			ELSE
178			DO i=1,iv
179			qtmp1(i,k) = qtmp1(i,k) / qtmp1(0,k)
180			ENDDO
181			C Variance :
182			qtmp1(iv,k) = qtmp1(iv,k) - qtmp1(1,k)*qtmp1(1,k)
183			C Standard deviation :
184			IF (qtmp1(iv,k).GT.0.) qtmp1(iv,k) = SQRT(qtmp1(iv,k))
185			ENDIF
186			ENDDO
187	jmc	1.3	C return global (& vertically integrated) volume in qtmp2(0,0):
188			qtmp2(0,0) = qtmp1(0,0)
189	jmc	1.1	ELSE
190			C With counter diagnostics => average = Sum / Sum(counter) :
191			DO k=0,klev
192			kCnt = min(k,kMlev)
193			IF ( qtmp2(0,kCnt).LE.0. ) THEN
194			DO i=1,nStats
195			qtmp1(i,k) = undef
196			ENDDO
197			ELSEIF ( qtmp2(1,kCnt).LE.0. ) THEN
198			DO i=1,iv
199			qtmp1(i,k) = undef
200			ENDDO
201			ELSE
202			DO i=1,iv
203			qtmp1(i,k) = qtmp1(i,k) / qtmp2(1,kCnt)
204			ENDDO
205			C jmc: looks like there is a Pb with how Variance is computed
206			C Variance :
207			qtmp1(iv,k) = qtmp1(iv,k) - qtmp1(1,k)*qtmp1(1,k)
208			C Standard deviation :
209			IF (qtmp1(iv,k).GT.0.) qtmp1(iv,k) = SQRT(qtmp1(iv,k))
210			ENDIF
211			ENDDO
212			ENDIF
213
214			ENDIF
215
216			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
217
218			RETURN
219			END