/[MITgcm]/MITgcm/pkg/diagnostics/diagstats_global.F
ViewVC logotype

Annotation of /MITgcm/pkg/diagnostics/diagstats_global.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.3 - (hide annotations) (download)
Sun Jul 10 00:52:12 2005 UTC (18 years, 10 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint58l_post, checkpoint57t_post, checkpoint57o_post, checkpoint58e_post, mitgcm_mapl_00, checkpoint58u_post, checkpoint58w_post, checkpoint57m_post, checkpoint57s_post, checkpoint58r_post, checkpoint57y_post, checkpoint58n_post, checkpoint58x_post, checkpoint58t_post, checkpoint58h_post, checkpoint57y_pre, checkpoint58q_post, checkpoint57v_post, checkpoint58j_post, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint59j, checkpoint57r_post, checkpoint59, checkpoint58, checkpoint58f_post, checkpoint57x_post, checkpoint57n_post, checkpoint58d_post, checkpoint58c_post, checkpoint57w_post, checkpoint57p_post, checkpint57u_post, checkpoint58a_post, checkpoint58i_post, checkpoint57q_post, checkpoint58g_post, checkpoint58o_post, checkpoint57z_post, checkpoint58y_post, checkpoint58k_post, checkpoint58v_post, checkpoint58s_post, checkpoint58p_post, checkpoint58b_post, checkpoint58m_post, checkpoint57l_post
Changes since 1.2: +4 -2 lines 
refine cases where WARNING is issued for an empty diagnostics.
1 jmc 1.3 C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagstats_global.F,v 1.2 2005/06/26 16:51:49 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    
52     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
53    
54     C-- Initialize to zero :
55     DO k=0,nLev
56     DO i=0,nStats
57     qtmp1(i,k) = 0.
58     qtmp2(i,k) = 0.
59     ENDDO
60     ENDDO
61    
62     klev = kdiag(ndId)
63     IF ( mate.GT.0 ) kMlev = kdiag(mate)
64    
65     IF (klev.LE.nLev) THEN
66     C--- Compute global statistics :
67    
68     C-- Retrieve tile statistics first
69     DO bj=myByLo(myThid),myByHi(myThid)
70     DO bi=myBxLo(myThid),myBxHi(myThid)
71    
72     DO k=1,klev
73 jmc 1.2 kd = iSp + k - 1
74 jmc 1.1 IF ( qSdiag(0,jReg,kd,bi,bj).GT.0. ) THEN
75     IF ( qtmp1(0,k).LE.0. ) THEN
76     DO i=0,nStats
77     qtmp1(i,k) = qSdiag(i,jReg,kd,bi,bj)
78     ENDDO
79     ELSE
80     DO i=0,iv
81     qtmp1(i,k) = qtmp1(i,k) + qSdiag(i,jReg,kd,bi,bj)
82     ENDDO
83     qtmp1(im,k) = MIN( qtmp1(im,k),qSdiag(im,jReg,kd,bi,bj) )
84     qtmp1(ix,k) = MAX( qtmp1(ix,k),qSdiag(ix,jReg,kd,bi,bj) )
85     ENDIF
86     ENDIF
87     ENDDO
88     IF ( mate.GT.0 ) THEN
89     DO k=1,kMlev
90 jmc 1.2 kd = iSm + k - 1
91 jmc 1.1 IF ( qSdiag(0,jReg,kd,bi,bj).GT.0. ) THEN
92 jmc 1.3 IF ( qtmp2(0,k).LE.0. ) THEN
93 jmc 1.1 DO i=0,nStats
94     qtmp2(i,k) = qSdiag(i,jReg,kd,bi,bj)
95     ENDDO
96     ELSE
97     DO i=0,iv
98     qtmp2(i,k) = qtmp2(i,k) + qSdiag(i,jReg,kd,bi,bj)
99     ENDDO
100     qtmp2(im,k) = MIN( qtmp2(im,k),qSdiag(im,jReg,kd,bi,bj) )
101     qtmp2(ix,k) = MAX( qtmp2(ix,k),qSdiag(ix,jReg,kd,bi,bj) )
102     ENDIF
103     ENDIF
104     ENDDO
105     ENDIF
106    
107     C- end tile index loops
108     ENDDO
109     ENDDO
110    
111     C-- Global min,max & sum (at each level) over all thread & processors :
112     DO k=1,klev
113     DO i=0,iv
114     _GLOBAL_SUM_R8(qtmp1(i,k),myThid)
115     ENDDO
116     qtmp1(im,k) = -qtmp1(im,k)
117     _GLOBAL_MAX_R8(qtmp1(im,k),myThid)
118     qtmp1(im,k) = -qtmp1(im,k)
119     _GLOBAL_MAX_R8(qtmp1(ix,k),myThid)
120     ENDDO
121     IF ( mate.GT.0 ) THEN
122     DO k=1,kMlev
123     DO i=0,iv
124     _GLOBAL_SUM_R8(qtmp2(i,k),myThid)
125     ENDDO
126     qtmp2(im,k) = -qtmp2(im,k)
127     _GLOBAL_MAX_R8(qtmp2(im,k),myThid)
128     qtmp2(im,k) = -qtmp2(im,k)
129     _GLOBAL_MAX_R8(qtmp2(ix,k),myThid)
130     ENDDO
131     ENDIF
132    
133     C-- Vertical integral, min & max :
134     DO k=1,klev
135     IF ( qtmp1(0,0).LE.0. ) THEN
136     DO i=0,nStats
137     qtmp1(i,0) = qtmp1(i,k)
138     ENDDO
139     ELSE
140     DO i=0,iv
141     qtmp1(i,0) = qtmp1(i,0) + qtmp1(i,k)
142     ENDDO
143     qtmp1(im,0) = MIN(qtmp1(im,0),qtmp1(im,k))
144     qtmp1(ix,0) = MAX(qtmp1(ix,0),qtmp1(ix,k))
145     ENDIF
146     ENDDO
147     IF ( mate.GT.0 ) THEN
148     DO k=1,kMlev
149     IF ( qtmp2(0,0).LE.0. ) THEN
150     DO i=0,nStats
151     qtmp2(i,0) = qtmp2(i,k)
152     ENDDO
153     ELSE
154     DO i=0,iv
155     qtmp2(i,0) = qtmp2(i,0) + qtmp2(i,k)
156     ENDDO
157     qtmp2(im,0) = MIN(qtmp2(im,0),qtmp2(im,k))
158     qtmp2(ix,0) = MAX(qtmp2(ix,0),qtmp2(ix,k))
159     ENDIF
160     ENDDO
161     ENDIF
162    
163     C-- Average, Standard.Dev.:
164     C- no counter diagnostics => average = Sum / vol :
165     IF ( mate.EQ.0 ) THEN
166     DO k=0,klev
167     IF ( qtmp1(0,k).LE.0. ) THEN
168     DO i=1,nStats
169     qtmp1(i,k) = undef
170     ENDDO
171     ELSE
172     DO i=1,iv
173     qtmp1(i,k) = qtmp1(i,k) / qtmp1(0,k)
174     ENDDO
175     C Variance :
176     qtmp1(iv,k) = qtmp1(iv,k) - qtmp1(1,k)*qtmp1(1,k)
177     C Standard deviation :
178     IF (qtmp1(iv,k).GT.0.) qtmp1(iv,k) = SQRT(qtmp1(iv,k))
179     ENDIF
180     ENDDO
181 jmc 1.3 C return global (& vertically integrated) volume in qtmp2(0,0):
182     qtmp2(0,0) = qtmp1(0,0)
183 jmc 1.1 ELSE
184     C With counter diagnostics => average = Sum / Sum(counter) :
185     DO k=0,klev
186     kCnt = min(k,kMlev)
187     IF ( qtmp2(0,kCnt).LE.0. ) THEN
188     DO i=1,nStats
189     qtmp1(i,k) = undef
190     ENDDO
191     ELSEIF ( qtmp2(1,kCnt).LE.0. ) THEN
192     DO i=1,iv
193     qtmp1(i,k) = undef
194     ENDDO
195     ELSE
196     DO i=1,iv
197     qtmp1(i,k) = qtmp1(i,k) / qtmp2(1,kCnt)
198     ENDDO
199     C jmc: looks like there is a Pb with how Variance is computed
200     C Variance :
201     qtmp1(iv,k) = qtmp1(iv,k) - qtmp1(1,k)*qtmp1(1,k)
202     C Standard deviation :
203     IF (qtmp1(iv,k).GT.0.) qtmp1(iv,k) = SQRT(qtmp1(iv,k))
204     ENDIF
205     ENDDO
206     ENDIF
207    
208     ENDIF
209    
210     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
211    
212     RETURN
213     END
  ViewVC Help
Powered by ViewVC 1.1.22