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

Contents of /MITgcm/pkg/diagnostics/diagstats_local.F

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


Revision 1.1 - (show annotations) (download)
Fri May 20 07:28:52 2005 UTC (18 years, 11 months ago) by jmc
Branch: MAIN
Add new capability: compute & write Global/Regional & per level statistics

1 C $Header: $
2 C $Name: $
3
4 #include "DIAG_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: DIAGSTATS_LOCAL
8 C !INTERFACE:
9 SUBROUTINE DIAGSTATS_LOCAL(
10 U statFld,
11 I inpFld,
12 I sizI1,sizI2,sizJ1,sizJ2,sizK,sizTx,sizTy,
13 I iRun,jRun,kIn,biIn,bjIn,
14 I k,bi,bj, region2fill, ndId, parsFld,
15 I myThid)
16
17 C !DESCRIPTION:
18 C Update array statFld
19 C by adding statistics over the range [1:iRun],[1:jRun]
20 C from input field array inpFld
21 C- note:
22 C a) this S/R should not see DIAGNOSTICS pkg commons blocks (in DIAGNOSTICS.h)
23 C b) for main grid variables, get area & weigting factors (to compute global mean)
24 C from the main common blocks.
25 C c) for other type of grids, call a specifics S/R which include its own
26 C grid common blocks
27
28 C !USES:
29 IMPLICIT NONE
30
31 #include "EEPARAMS.h"
32 #include "SIZE.h"
33 #include "DIAGNOSTICS_SIZE.h"
34 #include "PARAMS.h"
35 #include "GRID.h"
36 #include "SURFACE.h"
37 #ifdef ALLOW_FIZHI
38 #include "fizhi_SIZE.h"
39 #include "gridalt_mapping.h"
40 #endif
41
42 C !INPUT/OUTPUT PARAMETERS:
43 C == Routine Arguments ==
44 C statFld :: cumulative statistics array (updated)
45 C inpFld :: input field array to process (compute stats & add to statFld)
46 C sizI1,sizI2 :: size of inpFld array: 1rst index range (min,max)
47 C sizJ1,sizJ2 :: size of inpFld array: 2nd index range (min,max)
48 C sizK :: size of inpFld array: 3rd dimension
49 C sizTx,sizTy :: size of inpFld array: tile dimensions
50 C iRun,jRun :: range of 1rst & 2nd index
51 C kIn :: level index of inFld array to porcess
52 C biIn,bjIn :: tile indices of inFld array to process
53 C k,bi,bj :: level and tile indices used for weighting (mask,area ...)
54 C region2fill :: indicates whether to compute statistics over this region
55 C ndId :: Diagnostics Id Number (in available diag. list)
56 C parsFld :: parser field with characteristics of the diagnostics
57 C myThid :: my Thread Id number
58 _RL statFld(0:nStats,0:nRegions)
59 INTEGER sizI1,sizI2,sizJ1,sizJ2
60 INTEGER sizK,sizTx,sizTy
61 _RL inpFld(sizI1:sizI2,sizJ1:sizJ2,sizK,sizTx,sizTy)
62 INTEGER iRun, jRun, kIn, biIn, bjIn
63 INTEGER k, bi, bj, ndId
64 INTEGER region2fill(0:nRegions)
65 CHARACTER*16 parsFld
66 INTEGER myThid
67 CEOP
68
69 C !LOCAL VARIABLES:
70 C i,j :: loop indices
71 INTEGER i, n, km
72 INTEGER im, ix, iv
73 PARAMETER ( iv = nStats - 2 , im = nStats - 1 , ix = nStats )
74 LOGICAL exclSpVal
75 LOGICAL useWeight
76 _RL statLoc(0:nStats)
77 _RL drLoc
78 _RL specialVal
79 _RL getcon
80 EXTERNAL getcon
81
82 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
83
84 useWeight = .FALSE.
85 exclSpVal = .FALSE.
86 specialVal = 0.
87 IF ( useFIZHI ) THEN
88 exclSpVal = .TRUE.
89 specialVal = getcon('UNDEF')
90 ENDIF
91
92 DO n=0,nRegions
93 IF (region2fill(n).NE.0) THEN
94 C--- Compute statistics for this tile, level and region:
95
96 C- case of a special region (no specific regional mask)
97 IF ( n.EQ.0 ) THEN
98
99 IF ( parsFld(10:10) .EQ. 'R' ) THEN
100
101 drLoc = drF(k)
102 IF ( parsFld(9:9).EQ.'L') drLoc = drC(k)
103 IF ( parsFld(9:9).EQ.'U') drLoc = drC(MIN(k+1,Nr))
104 IF ( parsFld(9:9).EQ.'M') useWeight = .TRUE.
105
106 IF ( parsFld(2:2).EQ.'U' ) THEN
107 CALL DIAGSTATS_R_CALC(
108 O statLoc,
109 I inpFld(sizI1,sizJ1,kIn,biIn,bjIn),
110 I sizI1,sizI2,sizJ1,sizJ2,iRun,jRun,
111 I maskH(1-Olx,1-Oly,bi,bj), maskW(1-Olx,1-Oly,k,bi,bj),
112 I hFacW(1-Olx,1-Oly,k,bi,bj), rAw(1-Olx,1-Oly,bi,bj),
113 I drLoc, specialVal, exclSpVal, useWeight, myThid )
114 c I drLoc, k,bi,bj, parsFld, myThid )
115 ELSEIF ( parsFld(2:2).EQ.'V' ) THEN
116 CALL DIAGSTATS_R_CALC(
117 O statLoc,
118 I inpFld(sizI1,sizJ1,kIn,biIn,bjIn),
119 I sizI1,sizI2,sizJ1,sizJ2,iRun,jRun,
120 I maskH(1-Olx,1-Oly,bi,bj), maskS(1-Olx,1-Oly,k,bi,bj),
121 I hFacS(1-Olx,1-Oly,k,bi,bj), rAs(1-Olx,1-Oly,bi,bj),
122 I drLoc, specialVal, exclSpVal, useWeight, myThid )
123 ELSE
124 CALL DIAGSTATS_R_CALC(
125 O statLoc,
126 I inpFld(sizI1,sizJ1,kIn,biIn,bjIn),
127 I sizI1,sizI2,sizJ1,sizJ2,iRun,jRun,
128 I maskH(1-Olx,1-Oly,bi,bj), maskC(1-Olx,1-Oly,k,bi,bj),
129 I hFacC(1-Olx,1-Oly,k,bi,bj), rA(1-Olx,1-Oly,bi,bj),
130 I drLoc, specialVal, exclSpVal, useWeight, myThid )
131 ENDIF
132
133 #ifdef ALLOW_FIZHI
134 c ELSEIF ( parsFld(10:10) .EQ. 'L' ) THEN
135 ELSEIF ( parsFld(10:10) .EQ. 'M' ) THEN
136 drLoc = 1. _d 0
137 km = 1 + Nrphys - k
138 CALL DIAGSTATS_R_CALC(
139 O statLoc,
140 I inpFld(sizI1,sizJ1,kIn,biIn,bjIn),
141 I sizI1,sizI2,sizJ1,sizJ2,iRun,jRun,
142 I maskH(1-Olx,1-Oly,bi,bj), maskH(1-Olx,1-Oly,bi,bj),
143 I dpphys(1-Olx,1-Oly,km,bi,bj), rA(1-Olx,1-Oly,bi,bj),
144 I drLoc, specialVal, exclSpVal, useWeight, myThid )
145 #endif
146 #ifdef ALLOW_LAND
147 c ELSEIF ( parsFld(10:10) .EQ. 'G' ) THEN
148 #endif
149 c ELSEIF ( parsFld(10:10) .EQ. 'I' ) THEN
150 c ELSEIF ( parsFld(10:10) .EQ. '1' ) THEN
151 ELSE
152
153 km = 1
154 IF ( usingPCoords ) km = Nr
155 drLoc = 1. _d 0
156 IF ( parsFld(2:2).EQ.'U' ) THEN
157 CALL DIAGSTATS_R_CALC(
158 O statLoc,
159 I inpFld(sizI1,sizJ1,kIn,biIn,bjIn),
160 I sizI1,sizI2,sizJ1,sizJ2,iRun,jRun,
161 I maskH(1-Olx,1-Oly,bi,bj), maskW(1-Olx,1-Oly,km,bi,bj),
162 I maskW(1-Olx,1-Oly,km,bi,bj),rAw(1-Olx,1-Oly,bi,bj),
163 I drLoc, specialVal, exclSpVal, useWeight, myThid )
164 ELSEIF ( parsFld(2:2).EQ.'V' ) THEN
165 CALL DIAGSTATS_R_CALC(
166 O statLoc,
167 I inpFld(sizI1,sizJ1,kIn,biIn,bjIn),
168 I sizI1,sizI2,sizJ1,sizJ2,iRun,jRun,
169 I maskH(1-Olx,1-Oly,bi,bj), maskS(1-Olx,1-Oly,km,bi,bj),
170 I maskS(1-Olx,1-Oly,km,bi,bj),rAs(1-Olx,1-Oly,bi,bj),
171 I drLoc, specialVal, exclSpVal, useWeight, myThid )
172 ELSE
173 CALL DIAGSTATS_R_CALC(
174 O statLoc,
175 I inpFld(sizI1,sizJ1,kIn,biIn,bjIn),
176 I sizI1,sizI2,sizJ1,sizJ2,iRun,jRun,
177 I maskH(1-Olx,1-Oly,bi,bj), maskH(1-Olx,1-Oly,bi,bj),
178 I maskH(1-Olx,1-Oly,bi,bj), rA(1-Olx,1-Oly,bi,bj),
179 I drLoc, specialVal, exclSpVal, useWeight, myThid )
180 ENDIF
181
182 ENDIF
183
184 C Update cumulative statistics array
185 IF ( statLoc(0).GT.0. ) THEN
186 IF ( statFld(0,n).LE.0. ) THEN
187 statFld(im,n) = statLoc(im)
188 statFld(ix,n) = statLoc(ix)
189 ELSE
190 statFld(im,n) = MIN( statFld(im,n), statLoc(im) )
191 statFld(ix,n) = MAX( statFld(ix,n), statLoc(ix) )
192 ENDIF
193 DO i=0,iv
194 statFld(i,n) = statFld(i,n) + statLoc(i)
195 ENDDO
196 ENDIF
197
198 ENDIF
199
200 C--- processing region "n" ends here.
201 ENDIF
202 ENDDO
203
204 RETURN
205 END
206
207 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
208
209 CBOP
210 C !ROUTINE: DIAGSTATS_CALC
211 C !INTERFACE:
212 SUBROUTINE DIAGSTATS_R_CALC(
213 U statArr,
214 I inpArr,
215 I sizI1,sizI2,sizJ1,sizJ2, iRun,jRun,
216 I regMask, arrMask, arrhFac, arrArea,
217 I arrDr, specialVal, exclSpVal, useWeight,
218 I myThid )
219 c I arrDr, k,bi,bj, parsFld, myThid )
220
221 C !DESCRIPTION:
222 C Compute statistics for this tile, level, region
223
224 C !USES:
225 IMPLICIT NONE
226
227 #include "EEPARAMS.h"
228 #include "SIZE.h"
229 #include "DIAGNOSTICS_SIZE.h"
230 c #include "PARAMS.h"
231 c #include "GRID.h"
232 c #include "SURFACE.h"
233
234 C !INPUT/OUTPUT PARAMETERS:
235 C == Routine Arguments ==
236 C statArr :: cumulative statistics array (updated)
237 C inpArr :: input field array to process (compute stats & add to statFld)
238 C sizI1,sizI2 :: size of inpArr array: 1rst index range (min,max)
239 C sizJ1,sizJ2 :: size of inpArr array: 2nd index range (min,max)
240 C iRun,jRun :: range of 1rst & 2nd index to process
241 C regMask :: regional mask
242 C arrMask :: mask for this input array
243 C arrhFac :: weight factor (horizontally varying)
244 C arrArea :: Area weighting factor
245 C arrDr :: uniform weighting factor
246 C specialVal :: special value in input array (to exclude if exclSpVal=T)
247 C exclSpVal :: if T, exclude "specialVal" in input array
248 C useWeight :: use weight factor "arrhFac"
249 Cc k,bi,bj :: level and tile indices used for weighting (mask,area ...)
250 Cc parsFld :: parser field with characteristics of the diagnostics
251 C myThid :: my Thread Id number
252 _RL statArr(0:nStats)
253 INTEGER sizI1,sizI2,sizJ1,sizJ2
254 INTEGER iRun, jRun
255 _RL inpArr (sizI1:sizI2,sizJ1:sizJ2)
256 _RS regMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
257 _RS arrMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
258 _RS arrhFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
259 _RS arrArea(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
260 _RL arrDr
261 _RL specialVal
262 LOGICAL exclSpVal
263 LOGICAL useWeight
264 c INTEGER k, bi, bj
265 c CHARACTER*16 parsFld
266 INTEGER myThid
267 CEOP
268
269 C !LOCAL VARIABLES:
270 C i,j :: loop indices
271 INTEGER i, j, n
272 INTEGER im, ix, iv
273 PARAMETER ( iv = nStats - 2 , im = nStats - 1 , ix = nStats )
274 _RL tmpVol
275
276 DO n=0,nStats
277 statArr(n) = 0.
278 ENDDO
279
280 IF ( exclSpVal ) THEN
281
282 DO j = 1,jRun
283 DO i = 1,iRun
284 IF (arrMask(i,j).NE.0. .AND. inpArr(i,j).NE.specialVal) THEN
285 IF ( statArr(0).EQ.0. ) THEN
286 statArr(im) = inpArr(i,j)
287 statArr(ix) = inpArr(i,j)
288 ELSE
289 statArr(im) = MIN(inpArr(i,j),statArr(im))
290 statArr(ix) = MAX(inpArr(i,j),statArr(ix))
291 ENDIF
292 IF ( useWeight ) THEN
293 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)
294 ELSE
295 tmpVol = arrDr*arrArea(i,j)
296 ENDIF
297 statArr(0) = statArr(0) + tmpVol
298 statArr(1) = statArr(1) + tmpVol*inpArr(i,j)
299 statArr(2) = statArr(2) + tmpVol*inpArr(i,j)*inpArr(i,j)
300 ENDIF
301 ENDDO
302 ENDDO
303
304 ELSE
305
306 DO j = 1,jRun
307 DO i = 1,iRun
308 c IF ( regMask(i,j).NE.0. .AND. arrMask(i,j).NE.0. ) THEN
309 IF ( arrMask(i,j).NE.0. ) THEN
310 IF ( statArr(0).EQ.0. ) THEN
311 statArr(im) = inpArr(i,j)
312 statArr(ix) = inpArr(i,j)
313 ELSE
314 statArr(im) = MIN(inpArr(i,j),statArr(im))
315 statArr(ix) = MAX(inpArr(i,j),statArr(ix))
316 ENDIF
317 IF ( useWeight ) THEN
318 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)
319 ELSE
320 tmpVol = arrDr*arrArea(i,j)
321 ENDIF
322 statArr(0) = statArr(0) + tmpVol
323 statArr(1) = statArr(1) + tmpVol*inpArr(i,j)
324 statArr(2) = statArr(2) + tmpVol*inpArr(i,j)*inpArr(i,j)
325 ENDIF
326 ENDDO
327 ENDDO
328
329 ENDIF
330
331 RETURN
332 END

  ViewVC Help
Powered by ViewVC 1.1.22