pkg/diagnostics/diagstats_global.F

C $Header:  $
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
1	jmc	1.1	C $Header: $
2			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			I ndId, mate, myThid )
15
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			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 diagnostics long list)
31			C mate ..... counter diagnostic Id number if any ; 0 otherwise
32			C myThid ..... my thread Id number
33			_RL undef
34			INTEGER nLev, jReg, ndId, mate
35			INTEGER myThid
36
37			C !OUTPUT PARAMETERS:
38			C qtmp1 ..... AVERAGED DIAGNOSTIC QUANTITY
39			C qtmp2 ..... working array (used for counter mate statistics)
40			_RL qtmp1(0:nStats,0:nLev)
41			_RL qtmp2(0:nStats,0:nLev)
42			CEOP
43
44			C !LOCAL VARIABLES:
45			INTEGER im, ix, iv
46			PARAMETER ( iv = nStats - 2 , im = nStats - 1 , ix = nStats )
47			INTEGER bi, bj
48			INTEGER i, k, kd, kCnt, klev, kMlev
49
50			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
51
52			C-- Initialize to zero :
53			DO k=0,nLev
54			DO i=0,nStats
55			qtmp1(i,k) = 0.
56			qtmp2(i,k) = 0.
57			ENDDO
58			ENDDO
59
60			klev = kdiag(ndId)
61			IF ( mate.GT.0 ) kMlev = kdiag(mate)
62
63			IF (klev.LE.nLev) THEN
64			C--- Compute global statistics :
65
66			C-- Retrieve tile statistics first
67			DO bj=myByLo(myThid),myByHi(myThid)
68			DO bi=myBxLo(myThid),myBxHi(myThid)
69
70			DO k=1,klev
71			kd = iSdiag(ndId) + k - 1
72			IF ( qSdiag(0,jReg,kd,bi,bj).GT.0. ) THEN
73			IF ( qtmp1(0,k).LE.0. ) THEN
74			DO i=0,nStats
75			qtmp1(i,k) = qSdiag(i,jReg,kd,bi,bj)
76			ENDDO
77			ELSE
78			DO i=0,iv
79			qtmp1(i,k) = qtmp1(i,k) + qSdiag(i,jReg,kd,bi,bj)
80			ENDDO
81			qtmp1(im,k) = MIN( qtmp1(im,k),qSdiag(im,jReg,kd,bi,bj) )
82			qtmp1(ix,k) = MAX( qtmp1(ix,k),qSdiag(ix,jReg,kd,bi,bj) )
83			ENDIF
84			ENDIF
85			ENDDO
86			IF ( mate.GT.0 ) THEN
87			DO k=1,kMlev
88			kd = iSdiag(mate) + k - 1
89			IF ( qSdiag(0,jReg,kd,bi,bj).GT.0. ) THEN
90			IF ( qtmp1(0,k).LE.0. ) THEN
91			DO i=0,nStats
92			qtmp2(i,k) = qSdiag(i,jReg,kd,bi,bj)
93			ENDDO
94			ELSE
95			DO i=0,iv
96			qtmp2(i,k) = qtmp2(i,k) + qSdiag(i,jReg,kd,bi,bj)
97			ENDDO
98			qtmp2(im,k) = MIN( qtmp2(im,k),qSdiag(im,jReg,kd,bi,bj) )
99			qtmp2(ix,k) = MAX( qtmp2(ix,k),qSdiag(ix,jReg,kd,bi,bj) )
100			ENDIF
101			ENDIF
102			ENDDO
103			ENDIF
104
105			C- end tile index loops
106			ENDDO
107			ENDDO
108
109			C-- Global min,max & sum (at each level) over all thread & processors :
110			DO k=1,klev
111			DO i=0,iv
112			_GLOBAL_SUM_R8(qtmp1(i,k),myThid)
113			ENDDO
114			qtmp1(im,k) = -qtmp1(im,k)
115			_GLOBAL_MAX_R8(qtmp1(im,k),myThid)
116			qtmp1(im,k) = -qtmp1(im,k)
117			_GLOBAL_MAX_R8(qtmp1(ix,k),myThid)
118			ENDDO
119			IF ( mate.GT.0 ) THEN
120			DO k=1,kMlev
121			DO i=0,iv
122			_GLOBAL_SUM_R8(qtmp2(i,k),myThid)
123			ENDDO
124			qtmp2(im,k) = -qtmp2(im,k)
125			_GLOBAL_MAX_R8(qtmp2(im,k),myThid)
126			qtmp2(im,k) = -qtmp2(im,k)
127			_GLOBAL_MAX_R8(qtmp2(ix,k),myThid)
128			ENDDO
129			ENDIF
130
131			C-- Vertical integral, min & max :
132			DO k=1,klev
133			IF ( qtmp1(0,0).LE.0. ) THEN
134			DO i=0,nStats
135			qtmp1(i,0) = qtmp1(i,k)
136			ENDDO
137			ELSE
138			DO i=0,iv
139			qtmp1(i,0) = qtmp1(i,0) + qtmp1(i,k)
140			ENDDO
141			qtmp1(im,0) = MIN(qtmp1(im,0),qtmp1(im,k))
142			qtmp1(ix,0) = MAX(qtmp1(ix,0),qtmp1(ix,k))
143			ENDIF
144			ENDDO
145			IF ( mate.GT.0 ) THEN
146			DO k=1,kMlev
147			IF ( qtmp2(0,0).LE.0. ) THEN
148			DO i=0,nStats
149			qtmp2(i,0) = qtmp2(i,k)
150			ENDDO
151			ELSE
152			DO i=0,iv
153			qtmp2(i,0) = qtmp2(i,0) + qtmp2(i,k)
154			ENDDO
155			qtmp2(im,0) = MIN(qtmp2(im,0),qtmp2(im,k))
156			qtmp2(ix,0) = MAX(qtmp2(ix,0),qtmp2(ix,k))
157			ENDIF
158			ENDDO
159			ENDIF
160
161			C-- Average, Standard.Dev.:
162			C- no counter diagnostics => average = Sum / vol :
163			IF ( mate.EQ.0 ) THEN
164			DO k=0,klev
165			IF ( qtmp1(0,k).LE.0. ) THEN
166			DO i=1,nStats
167			qtmp1(i,k) = undef
168			ENDDO
169			ELSE
170			DO i=1,iv
171			qtmp1(i,k) = qtmp1(i,k) / qtmp1(0,k)
172			ENDDO
173			C Variance :
174			qtmp1(iv,k) = qtmp1(iv,k) - qtmp1(1,k)*qtmp1(1,k)
175			C Standard deviation :
176			IF (qtmp1(iv,k).GT.0.) qtmp1(iv,k) = SQRT(qtmp1(iv,k))
177			ENDIF
178			ENDDO
179			ELSE
180			C With counter diagnostics => average = Sum / Sum(counter) :
181			DO k=0,klev
182			kCnt = min(k,kMlev)
183			IF ( qtmp2(0,kCnt).LE.0. ) THEN
184			DO i=1,nStats
185			qtmp1(i,k) = undef
186			ENDDO
187			ELSEIF ( qtmp2(1,kCnt).LE.0. ) THEN
188			DO i=1,iv
189			qtmp1(i,k) = undef
190			ENDDO
191			ELSE
192			DO i=1,iv
193			qtmp1(i,k) = qtmp1(i,k) / qtmp2(1,kCnt)
194			ENDDO
195			C jmc: looks like there is a Pb with how Variance is computed
196			C Variance :
197			qtmp1(iv,k) = qtmp1(iv,k) - qtmp1(1,k)*qtmp1(1,k)
198			C Standard deviation :
199			IF (qtmp1(iv,k).GT.0.) qtmp1(iv,k) = SQRT(qtmp1(iv,k))
200			ENDIF
201			ENDDO
202			ENDIF
203
204			ENDIF
205
206			C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-\|--+----\|
207
208			RETURN
209			END