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

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

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


Revision 1.4 - (show annotations) (download)
Mon Dec 12 15:28:16 2011 UTC (12 years, 5 months ago) by mlosch
Branch: MAIN
CVS Tags: checkpoint63h, checkpoint63g
Changes since 1.3: +10 -4 lines
fix problem with non-initialise halos in tmp-fields for TARGET_NEC_SX
at the cost of extra loops and less vectorisation

1 C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagstats_calc.F,v 1.3 2011/11/21 09:37:10 mlosch Exp $
2 C $Name: $
3
4 #include "DIAG_OPTIONS.h"
5
6 CBOP
7 C !ROUTINE: DIAGSTATS_CALC
8 C !INTERFACE:
9 SUBROUTINE DIAGSTATS_CALC(
10 O statArr,
11 I inpArr, frcArr, scaleFact, power, useFract,
12 I regId, regMskVal,
13 I nStats,sizI1,sizI2,sizJ1,sizJ2, iRun,jRun,
14 I regMask, arrMask, arrhFac, arrArea,
15 I arrDr, specialVal, exclSpVal, useWeight,
16 I myThid )
17
18 C !DESCRIPTION:
19 C Compute statistics for this tile, level, region
20
21 C !USES:
22 IMPLICIT NONE
23
24 #include "EEPARAMS.h"
25 #include "SIZE.h"
26
27 C !INPUT/OUTPUT PARAMETERS:
28 C == Routine Arguments ==
29 C statArr :: output statistics array
30 C inpArr :: input field array to process (compute stats & add to statFld)
31 C frcArr :: fraction used for weighted-average diagnostics
32 C scaleFact :: scaling factor
33 C power :: option to fill-in with the field square (power=2)
34 C useFract :: if True, use fraction-weight
35 C regId :: region number Id
36 C regMskVal :: region-mask identificator value
37 C (point i,j belong to region "regId" <=> regMask(i,j) = regMskVal)
38 C nStats :: size of output array: statArr
39 C sizI1,sizI2 :: size of inpArr array: 1rst index range (min,max)
40 C sizJ1,sizJ2 :: size of inpArr array: 2nd index range (min,max)
41 C iRun,jRun :: range of 1rst & 2nd index to process
42 C regMask :: regional mask
43 C arrMask :: mask for this input array
44 C arrhFac :: weight factor (horizontally varying)
45 C arrArea :: Area weighting factor
46 C arrDr :: uniform weighting factor
47 C specialVal :: special value in input array (to exclude if exclSpVal=T)
48 C exclSpVal :: if T, exclude "specialVal" in input array
49 C useWeight :: use weight factor "arrhFac"
50 Cc k,bi,bj :: level and tile indices used for weighting (mask,area ...)
51 Cc parsFld :: parser field with characteristics of the diagnostics
52 C myThid :: my Thread Id number
53 INTEGER nStats,sizI1,sizI2,sizJ1,sizJ2
54 INTEGER iRun, jRun
55 _RL statArr(0:nStats)
56 _RL inpArr (sizI1:sizI2,sizJ1:sizJ2)
57 _RL frcArr (sizI1:sizI2,sizJ1:sizJ2)
58 _RL scaleFact
59 INTEGER power
60 LOGICAL useFract
61 INTEGER regId
62 _RS regMskVal
63 _RS regMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
64 _RS arrMask(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
65 _RS arrhFac(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
66 _RS arrArea(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
67 _RL arrDr
68 _RL specialVal
69 LOGICAL exclSpVal
70 LOGICAL useWeight
71 INTEGER myThid
72 CEOP
73
74 C !LOCAL VARIABLES:
75 C i,j :: loop indices
76 INTEGER i, j, n
77 INTEGER im, ix
78 #ifndef TARGET_NEC_SX
79 _RL tmpVol
80 _RL tmpFld
81 #else
82 C Extra variables and fields to support vectorization.
83 C This code also uses the intrinsic F90 routines SUM, MINVAL, MAXVAL
84 C and thus will not compile with a F77 compiler.
85 _RL arrMaskL(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
86 _RL tmpFld (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
87 _RL tmpVol (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
88 #endif
89 _RL tmpFac
90
91 im = nStats - 1
92 ix = nStats
93 DO n=0,nStats
94 statArr(n) = 0.
95 ENDDO
96 tmpFac = scaleFact
97 IF ( power.EQ.2) tmpFac = scaleFact*scaleFact
98
99 #ifndef TARGET_NEC_SX
100 IF ( regId.EQ.0 .AND. useFract .AND. exclSpVal ) THEN
101
102 DO j = 1,jRun
103 DO i = 1,iRun
104 IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0.
105 & .AND. inpArr(i,j).NE.specialVal ) THEN
106 IF ( power.EQ.2) THEN
107 tmpFld = tmpFac*inpArr(i,j)*inpArr(i,j)
108 ELSE
109 tmpFld = tmpFac*inpArr(i,j)
110 ENDIF
111 IF ( statArr(0).EQ.0. ) THEN
112 statArr(im) = tmpFld
113 statArr(ix) = tmpFld
114 ELSE
115 statArr(im) = MIN(tmpFld,statArr(im))
116 statArr(ix) = MAX(tmpFld,statArr(ix))
117 ENDIF
118 IF ( useWeight ) THEN
119 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)*frcArr(i,j)
120 ELSE
121 tmpVol = arrDr*arrArea(i,j)*frcArr(i,j)
122 ENDIF
123 statArr(0) = statArr(0) + tmpVol
124 statArr(1) = statArr(1) + tmpVol*tmpFld
125 statArr(2) = statArr(2) + tmpVol*tmpFld*tmpFld
126 ENDIF
127 ENDDO
128 ENDDO
129
130 ELSEIF ( regId.EQ.0 .AND. useFract ) THEN
131
132 DO j = 1,jRun
133 DO i = 1,iRun
134 IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0. ) THEN
135 IF ( power.EQ.2) THEN
136 tmpFld = tmpFac*inpArr(i,j)*inpArr(i,j)
137 ELSE
138 tmpFld = tmpFac*inpArr(i,j)
139 ENDIF
140 IF ( statArr(0).EQ.0. ) THEN
141 statArr(im) = tmpFld
142 statArr(ix) = tmpFld
143 ELSE
144 statArr(im) = MIN(tmpFld,statArr(im))
145 statArr(ix) = MAX(tmpFld,statArr(ix))
146 ENDIF
147 IF ( useWeight ) THEN
148 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)*frcArr(i,j)
149 ELSE
150 tmpVol = arrDr*arrArea(i,j)*frcArr(i,j)
151 ENDIF
152 statArr(0) = statArr(0) + tmpVol
153 statArr(1) = statArr(1) + tmpVol*tmpFld
154 statArr(2) = statArr(2) + tmpVol*tmpFld*tmpFld
155 ENDIF
156 ENDDO
157 ENDDO
158
159 ELSEIF ( regId.EQ.0 .AND. exclSpVal ) THEN
160
161 DO j = 1,jRun
162 DO i = 1,iRun
163 IF ( arrMask(i,j).NE.0.
164 & .AND. inpArr(i,j).NE.specialVal ) THEN
165 IF ( power.EQ.2) THEN
166 tmpFld = tmpFac*inpArr(i,j)*inpArr(i,j)
167 ELSE
168 tmpFld = tmpFac*inpArr(i,j)
169 ENDIF
170 IF ( statArr(0).EQ.0. ) THEN
171 statArr(im) = tmpFld
172 statArr(ix) = tmpFld
173 ELSE
174 statArr(im) = MIN(tmpFld,statArr(im))
175 statArr(ix) = MAX(tmpFld,statArr(ix))
176 ENDIF
177 IF ( useWeight ) THEN
178 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)
179 ELSE
180 tmpVol = arrDr*arrArea(i,j)
181 ENDIF
182 statArr(0) = statArr(0) + tmpVol
183 statArr(1) = statArr(1) + tmpVol*tmpFld
184 statArr(2) = statArr(2) + tmpVol*tmpFld*tmpFld
185 ENDIF
186 ENDDO
187 ENDDO
188
189 ELSEIF ( regId.EQ.0 ) THEN
190
191 DO j = 1,jRun
192 DO i = 1,iRun
193 IF ( arrMask(i,j).NE.0. ) THEN
194 IF ( power.EQ.2) THEN
195 tmpFld = tmpFac*inpArr(i,j)*inpArr(i,j)
196 ELSE
197 tmpFld = tmpFac*inpArr(i,j)
198 ENDIF
199 IF ( statArr(0).EQ.0. ) THEN
200 statArr(im) = tmpFld
201 statArr(ix) = tmpFld
202 ELSE
203 statArr(im) = MIN(tmpFld,statArr(im))
204 statArr(ix) = MAX(tmpFld,statArr(ix))
205 ENDIF
206 IF ( useWeight ) THEN
207 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)
208 ELSE
209 tmpVol = arrDr*arrArea(i,j)
210 ENDIF
211 statArr(0) = statArr(0) + tmpVol
212 statArr(1) = statArr(1) + tmpVol*tmpFld
213 statArr(2) = statArr(2) + tmpVol*tmpFld*tmpFld
214 ENDIF
215 ENDDO
216 ENDDO
217
218 ELSEIF ( useFract .AND. exclSpVal ) THEN
219
220 DO j = 1,jRun
221 DO i = 1,iRun
222 IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0.
223 & .AND. inpArr(i,j).NE.specialVal
224 & .AND. regMask(i,j).EQ.regMskVal ) THEN
225 IF ( power.EQ.2) THEN
226 tmpFld = tmpFac*inpArr(i,j)*inpArr(i,j)
227 ELSE
228 tmpFld = tmpFac*inpArr(i,j)
229 ENDIF
230 IF ( statArr(0).EQ.0. ) THEN
231 statArr(im) = tmpFld
232 statArr(ix) = tmpFld
233 ELSE
234 statArr(im) = MIN(tmpFld,statArr(im))
235 statArr(ix) = MAX(tmpFld,statArr(ix))
236 ENDIF
237 IF ( useWeight ) THEN
238 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)*frcArr(i,j)
239 ELSE
240 tmpVol = arrDr*arrArea(i,j)*frcArr(i,j)
241 ENDIF
242 statArr(0) = statArr(0) + tmpVol
243 statArr(1) = statArr(1) + tmpVol*tmpFld
244 statArr(2) = statArr(2) + tmpVol*tmpFld*tmpFld
245 ENDIF
246 ENDDO
247 ENDDO
248
249 ELSEIF ( useFract ) THEN
250
251 DO j = 1,jRun
252 DO i = 1,iRun
253 IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0.
254 & .AND. regMask(i,j).EQ.regMskVal ) THEN
255 IF ( power.EQ.2) THEN
256 tmpFld = tmpFac*inpArr(i,j)*inpArr(i,j)
257 ELSE
258 tmpFld = tmpFac*inpArr(i,j)
259 ENDIF
260 IF ( statArr(0).EQ.0. ) THEN
261 statArr(im) = tmpFld
262 statArr(ix) = tmpFld
263 ELSE
264 statArr(im) = MIN(tmpFld,statArr(im))
265 statArr(ix) = MAX(tmpFld,statArr(ix))
266 ENDIF
267 IF ( useWeight ) THEN
268 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)*frcArr(i,j)
269 ELSE
270 tmpVol = arrDr*arrArea(i,j)*frcArr(i,j)
271 ENDIF
272 statArr(0) = statArr(0) + tmpVol
273 statArr(1) = statArr(1) + tmpVol*tmpFld
274 statArr(2) = statArr(2) + tmpVol*tmpFld*tmpFld
275 ENDIF
276 ENDDO
277 ENDDO
278
279 ELSEIF ( exclSpVal ) THEN
280
281 DO j = 1,jRun
282 DO i = 1,iRun
283 IF ( arrMask(i,j).NE.0.
284 & .AND. inpArr(i,j).NE.specialVal
285 & .AND. regMask(i,j).EQ.regMskVal ) THEN
286 IF ( power.EQ.2) THEN
287 tmpFld = tmpFac*inpArr(i,j)*inpArr(i,j)
288 ELSE
289 tmpFld = tmpFac*inpArr(i,j)
290 ENDIF
291 IF ( statArr(0).EQ.0. ) THEN
292 statArr(im) = tmpFld
293 statArr(ix) = tmpFld
294 ELSE
295 statArr(im) = MIN(tmpFld,statArr(im))
296 statArr(ix) = MAX(tmpFld,statArr(ix))
297 ENDIF
298 IF ( useWeight ) THEN
299 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)
300 ELSE
301 tmpVol = arrDr*arrArea(i,j)
302 ENDIF
303 statArr(0) = statArr(0) + tmpVol
304 statArr(1) = statArr(1) + tmpVol*tmpFld
305 statArr(2) = statArr(2) + tmpVol*tmpFld*tmpFld
306 ENDIF
307 ENDDO
308 ENDDO
309
310 ELSE
311
312 DO j = 1,jRun
313 DO i = 1,iRun
314 IF ( arrMask(i,j).NE.0.
315 & .AND. regMask(i,j).EQ.regMskVal ) THEN
316 IF ( power.EQ.2) THEN
317 tmpFld = tmpFac*inpArr(i,j)*inpArr(i,j)
318 ELSE
319 tmpFld = tmpFac*inpArr(i,j)
320 ENDIF
321 IF ( statArr(0).EQ.0. ) THEN
322 statArr(im) = tmpFld
323 statArr(ix) = tmpFld
324 ELSE
325 statArr(im) = MIN(tmpFld,statArr(im))
326 statArr(ix) = MAX(tmpFld,statArr(ix))
327 ENDIF
328 IF ( useWeight ) THEN
329 tmpVol = arrDr*arrhFac(i,j)*arrArea(i,j)
330 ELSE
331 tmpVol = arrDr*arrArea(i,j)
332 ENDIF
333 statArr(0) = statArr(0) + tmpVol
334 statArr(1) = statArr(1) + tmpVol*tmpFld
335 statArr(2) = statArr(2) + tmpVol*tmpFld*tmpFld
336 ENDIF
337 ENDDO
338 ENDDO
339
340 ENDIF
341
342 #else /* TARGET_NEC_SX defined */
343
344 arrMaskL = 0. _d 0
345
346 IF ( regId.EQ.0 .AND. useFract .AND. exclSpVal ) THEN
347
348 DO j = 1,jRun
349 DO i = 1,iRun
350 IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0.
351 & .AND. inpArr(i,j).NE.specialVal )
352 & arrMaskL(i,j) = 1. _d 0
353 ENDDO
354 ENDDO
355 IF ( useWeight ) THEN
356 tmpVol = arrhFac*arrArea*frcArr
357 ELSE
358 tmpVol = arrArea*frcArr
359 ENDIF
360
361 ELSEIF ( regId.EQ.0 .AND. useFract ) THEN
362
363 DO j = 1,jRun
364 DO i = 1,iRun
365 IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0.)
366 & arrMaskL(i,j) = 1. _d 0
367 ENDDO
368 ENDDO
369 IF ( useWeight ) THEN
370 tmpVol = arrhFac*arrArea*frcArr
371 ELSE
372 tmpVol = arrArea*frcArr
373 ENDIF
374
375 ELSEIF ( regId.EQ.0 .AND. exclSpVal ) THEN
376
377 DO j = 1,jRun
378 DO i = 1,iRun
379 IF ( arrMask(i,j).NE.0. .AND. inpArr(i,j).NE.specialVal )
380 & arrMaskL(i,j) = 1. _d 0
381 ENDDO
382 ENDDO
383 IF ( useWeight ) THEN
384 tmpVol = arrhFac*arrArea
385 ELSE
386 tmpVol = arrArea
387 ENDIF
388
389 ELSEIF ( regId.EQ.0 ) THEN
390
391 DO j = 1,jRun
392 DO i = 1,iRun
393 IF ( arrMask(i,j).NE.0. ) arrMaskL(i,j) = 1. _d 0
394 ENDDO
395 ENDDO
396 IF ( useWeight ) THEN
397 tmpVol = arrhFac*arrArea
398 ELSE
399 tmpVol = arrArea
400 ENDIF
401
402 ELSEIF ( useFract .AND. exclSpVal ) THEN
403
404 DO j = 1,jRun
405 DO i = 1,iRun
406 IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0.
407 & .AND. inpArr(i,j).NE.specialVal
408 & .AND. regMask(i,j).EQ.regMskVal )
409 & arrMaskL(i,j) = 1. _d 0
410 ENDDO
411 ENDDO
412 IF ( useWeight ) THEN
413 tmpVol = arrhFac*arrArea*frcArr
414 ELSE
415 tmpVol = arrArea*frcArr
416 ENDIF
417
418 ELSEIF ( useFract ) THEN
419
420 DO j = 1,jRun
421 DO i = 1,iRun
422 IF ( arrMask(i,j).NE.0. .AND. frcArr(i,j).NE.0.
423 & .AND. regMask(i,j).EQ.regMskVal )
424 & arrMaskL(i,j) = 1. _d 0
425 ENDDO
426 ENDDO
427 IF ( useWeight ) THEN
428 tmpVol = arrhFac*arrArea*frcArr
429 ELSE
430 tmpVol = arrArea*frcArr
431 ENDIF
432
433 ELSEIF ( exclSpVal ) THEN
434
435 DO j = 1,jRun
436 DO i = 1,iRun
437 IF ( arrMask(i,j).NE.0.
438 & .AND. inpArr(i,j).NE.specialVal
439 & .AND. regMask(i,j).EQ.regMskVal )
440 & arrMaskL(i,j) = 1. _d 0
441 ENDDO
442 ENDDO
443 IF ( useWeight ) THEN
444 tmpVol = arrhFac*arrArea
445 ELSE
446 tmpVol = arrArea
447 ENDIF
448
449 ELSE
450
451 DO j = 1,jRun
452 DO i = 1,iRun
453 IF ( arrMask(i,j).NE.0.
454 & .AND. regMask(i,j).EQ.regMskVal )
455 & arrMaskL(i,j) = 1. _d 0
456 ENDDO
457 ENDDO
458 IF ( useWeight ) THEN
459 tmpVol = arrhFac*arrArea
460 ELSE
461 tmpVol = arrArea
462 ENDIF
463
464 ENDIF
465 C inpArr can be undefined/non-initialised in overlaps, so we need
466 C to clean this fields first by copying the defined range to tmpFld
467 tmpFld = 0. _d 0
468 DO j = 1,jRun
469 DO i = 1,iRun
470 tmpFld(i,j) = inpArr(i,j)
471 ENDDO
472 ENDDO
473 IF ( power.EQ.2) THEN
474 tmpFld = tmpFld*tmpFld
475 ENDIF
476 C sum up the volume
477 tmpVol = tmpVol*arrMaskL
478 statArr(0) = SUM(tmpVol)*arrDr
479 C compute and sum up volume*field
480 tmpVol = tmpVol*tmpFld
481 statArr(1) = SUM(tmpVol)*tmpFac*arrDr
482 C compute and sum up volume*field**2
483 tmpVol = tmpVol*tmpFld
484 statArr(2) = SUM(tmpVol)*tmpFac*tmpFac*arrDr
485 statArr(im) = MINVAL(tmpFld, MASK = arrMaskL>0.)*tmpFac
486 statArr(ix) = MAXVAL(tmpFld, MASK = arrMaskL>0.)*tmpFac
487
488 #endif /* TARGET_NEC_SX */
489
490 RETURN
491 END

  ViewVC Help
Powered by ViewVC 1.1.22