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

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

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


Revision 1.5 - (hide annotations) (download)
Thu Sep 3 20:41:37 2009 UTC (14 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64y, checkpoint64x, checkpoint64z, checkpoint64q, checkpoint64p, checkpoint64s, checkpoint64r, checkpoint64u, checkpoint64t, checkpoint64w, checkpoint64v, checkpoint64i, checkpoint64h, checkpoint64k, checkpoint64j, checkpoint64m, checkpoint64l, checkpoint64o, checkpoint64n, checkpoint64a, checkpoint64c, checkpoint64b, checkpoint64e, checkpoint64d, checkpoint64g, checkpoint64f, checkpoint63p, checkpoint63q, checkpoint63r, checkpoint63s, checkpoint63l, checkpoint63m, checkpoint63n, checkpoint63o, checkpoint63h, checkpoint63i, checkpoint63j, checkpoint63k, checkpoint63d, checkpoint63e, checkpoint63f, checkpoint63g, checkpoint63a, checkpoint63b, checkpoint63c, checkpoint64, checkpoint65, checkpoint62, checkpoint63, checkpoint65a, checkpoint62c, checkpoint62b, checkpoint62a, checkpoint62g, checkpoint62f, checkpoint62e, checkpoint62d, checkpoint62k, checkpoint62j, checkpoint62i, checkpoint62h, checkpoint62o, checkpoint62n, checkpoint62m, checkpoint62l, checkpoint62s, checkpoint62r, checkpoint62q, checkpoint62p, checkpoint62w, checkpoint62v, checkpoint62u, checkpoint62t, checkpoint62z, checkpoint62y, checkpoint62x, checkpoint61v, checkpoint61w, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.4: +182 -10 lines
hack (with copy to RL array) for "RS" version of input field and fraction.

1 jmc 1.5 C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagstats_fill.F,v 1.4 2008/02/05 15:31:19 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "DIAG_OPTIONS.h"
5    
6 jmc 1.5 C-- File diagstats_fill.F:
7     C-- Contents:
8     C-- o DIAGSTATS_FILL
9     C-- o DIAGSTATS_TO_RL
10    
11     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
12    
13 jmc 1.1 CBOP
14     C !ROUTINE: DIAGSTATS_FILL
15     C !INTERFACE:
16 jmc 1.2 SUBROUTINE DIAGSTATS_FILL(
17 jmc 1.5 I inpFldRL, fracFldRL,
18     #ifndef REAL4_IS_SLOW
19     I inpFldRS, fracFldRS,
20     #endif
21     I scaleFact, power, arrType, nLevFract,
22 jmc 1.3 I ndId, kInQSd, region2fill, kLev, nLevs,
23     I bibjflg, biArg, bjArg, myThid )
24 jmc 1.1
25     C !DESCRIPTION:
26     C***********************************************************************
27 jmc 1.2 C compute statistics over 1 tile
28     C and increment the diagnostics array
29 jmc 1.3 C using a scaling factor & square option (power=2),
30 jmc 1.2 C and with the option to use a fraction-weight (assumed
31     C to be the counter-mate of the current diagnostics)
32 jmc 1.1 C***********************************************************************
33     C !USES:
34     IMPLICIT NONE
35    
36     C == Global variables ===
37     #include "EEPARAMS.h"
38     #include "SIZE.h"
39     #include "DIAGNOSTICS_SIZE.h"
40     #include "DIAGNOSTICS.h"
41    
42     C !INPUT PARAMETERS:
43     C***********************************************************************
44     C Arguments Description
45     C ----------------------
46 jmc 1.5 C inpFldRL :: Field to increment diagnostics array (arrType=0,1)
47     C fracFldRL :: fraction used for weighted average diagnostics (arrType=0,2)
48     C inpFldRS :: Field to increment diagnostics array (arrType=2,3)
49     C fracFldRS :: fraction used for weighted average diagnostics (arrType=1,3)
50 jmc 1.3 C scaleFact :: scaling factor
51     C power :: option to fill-in with the field square (power=2)
52 jmc 1.5 C arrType :: select which array & fraction (RL/RS) to process:
53     C 0: both RL ; 1: inpRL & fracRS ; 2: inpRS,fracRL ; 3: both RS
54     C nLevFract :: number of levels of the fraction field, =0: do not use fraction
55 jmc 1.3 C ndId :: Diagnostics Id Number (in available diag list) of diag to process
56     C kInQSd :: Pointer to the slot in qSdiag to fill
57     C region2fill :: array, indicates whether to compute statistics over region
58 jmc 1.1 C "j" (if region2fill(j)=1) or not (if region2fill(j)=0)
59 jmc 1.3 C kLev :: Integer flag for vertical levels:
60 jmc 1.1 C > 0 (any integer): WHICH single level to increment in qSdiag.
61     C 0,-1 to increment "nLevs" levels in qSdiag,
62 jmc 1.2 C 0 : fill-in in the same order as the input array
63 jmc 1.1 C -1: fill-in in reverse order.
64 jmc 1.3 C nLevs :: indicates Number of levels of the input field array
65 jmc 1.1 C (whether to fill-in all the levels (kLev<1) or just one (kLev>0))
66 jmc 1.3 C bibjflg :: Integer flag to indicate instructions for bi bj loop
67 jmc 1.1 C 0 indicates that the bi-bj loop must be done here
68     C 1 indicates that the bi-bj loop is done OUTSIDE
69     C 2 indicates that the bi-bj loop is done OUTSIDE
70     C AND that we have been sent a local array (with overlap regions)
71     C 3 indicates that the bi-bj loop is done OUTSIDE
72     C AND that we have been sent a local array
73     C AND that the array has no overlap region (interior only)
74     C NOTE - bibjflg can be NEGATIVE to indicate not to increment counter
75 jmc 1.3 C biArg :: X-direction tile number - used for bibjflg=1-3
76     C bjArg :: Y-direction tile number - used for bibjflg=1-3
77     C myThid :: my thread Id number
78 jmc 1.1 C***********************************************************************
79     C NOTE: User beware! If a local (1 tile only) array
80     C is sent here, bibjflg MUST NOT be set to 0
81     C or there will be out of bounds problems!
82     C***********************************************************************
83 jmc 1.5 _RL inpFldRL(*)
84     _RL fracFldRL(*)
85     #ifndef REAL4_IS_SLOW
86     _RS inpFldRS(*)
87     _RS fracFldRS(*)
88     #endif
89 jmc 1.2 _RL scaleFact
90 jmc 1.3 INTEGER power
91 jmc 1.5 INTEGER arrType
92 jmc 1.2 INTEGER nLevFract
93 jmc 1.1 INTEGER ndId, kInQSd
94     INTEGER region2fill(0:nRegions)
95     INTEGER kLev, nLevs, bibjflg, biArg, bjArg
96     INTEGER myThid
97     CEOP
98    
99     C !LOCAL VARIABLES:
100     C ===============
101 jmc 1.3 C useFract :: flag to increment (or not) with fraction-weigted inpFld
102 jmc 1.2 LOGICAL useFract
103     INTEGER sizF
104 jmc 1.1 INTEGER sizI1,sizI2,sizJ1,sizJ2
105     INTEGER sizTx,sizTy
106     INTEGER iRun, jRun, k, bi, bj
107     INTEGER kFirst, kLast
108     INTEGER kd, kd0, ksgn, kStore
109     CHARACTER*8 parms1
110     CHARACTER*(MAX_LEN_MBUF) msgBuf
111     INTEGER km, km0
112 jmc 1.5 #ifndef REAL4_IS_SLOW
113     _RL tmpFldRL( sNx+1,sNy+1)
114     _RL tmpFracRL(sNx+1,sNy+1)
115     #endif
116 jmc 1.1
117     C If-sequence to see if we are a valid and an active diagnostic
118     c IF ( ndId.NE.0 .AND. kInQSd.NE.0 ) THEN
119    
120     C- select range for 1rst & 2nd indices to accumulate
121 jmc 1.2 C depending on variable location on C-grid,
122 jmc 1.1 parms1 = gdiag(ndId)(1:8)
123     IF ( parms1(2:2).EQ.'Z' ) THEN
124     iRun = sNx+1
125     jRun = sNy+1
126     c ELSEIF ( parms1(2:2).EQ.'U' ) THEN
127     c iRun = sNx+1
128     c jRun = sNy
129     c ELSEIF ( parms1(2:2).EQ.'V' ) THEN
130     c iRun = sNx
131     c jRun = sNy+1
132     ELSE
133     iRun = sNx
134     jRun = sNy
135     ENDIF
136    
137     C- Dimension of the input array:
138     IF (ABS(bibjflg).EQ.3) THEN
139     sizI1 = 1
140     sizI2 = sNx
141     sizJ1 = 1
142     sizJ2 = sNy
143     iRun = sNx
144     jRun = sNy
145     ELSE
146     sizI1 = 1-OLx
147     sizI2 = sNx+OLx
148     sizJ1 = 1-OLy
149     sizJ2 = sNy+OLy
150     ENDIF
151     IF (ABS(bibjflg).GE.2) THEN
152     sizTx = 1
153     sizTy = 1
154     ELSE
155     sizTx = nSx
156     sizTy = nSy
157     ENDIF
158     C- Which part of inpFld to add : k = 3rd index,
159     C and do the loop >> do k=kFirst,kLast <<
160     IF (kLev.LE.0) THEN
161     kFirst = 1
162     kLast = nLevs
163     ELSEIF ( nLevs.EQ.1 ) THEN
164     kFirst = 1
165     kLast = 1
166     ELSEIF ( kLev.LE.nLevs ) THEN
167     kFirst = kLev
168     kLast = kLev
169     ELSE
170     STOP 'ABNORMAL END in DIAGSTATS_FILL: kLev > nLevs > 0'
171     ENDIF
172 jmc 1.2 C- Which part of qSdiag to update: kd = 3rd index,
173 jmc 1.1 C and do the loop >> do k=kFirst,kLast ; kd = kd0 + k*ksgn <<
174 jmc 1.2 C 1rst try this: for the mask: km = km0 + k*ksgn so that kd= km + kInQSd - 1
175 jmc 1.1 IF ( kLev.EQ.-1 ) THEN
176     ksgn = -1
177     kd0 = kInQSd + nLevs
178     km0 = 1 + nLevs
179     ELSEIF ( kLev.EQ.0 ) THEN
180     ksgn = 1
181     kd0 = kInQSd - 1
182     km0 = 0
183     ELSE
184     ksgn = 0
185     kd0 = kInQSd + kLev - 1
186     km0 = kLev
187     ENDIF
188 jmc 1.2 C- Set fraction-weight option :
189     useFract = nLevFract.GT.0
190     IF ( useFract ) THEN
191     sizF = nLevFract
192     ELSE
193     sizF = 1
194     ENDIF
195 jmc 1.1
196     C- Check for consistency with Nb of levels reserved in storage array
197     kStore = kd0 + MAX(ksgn*kFirst,ksgn*kLast) - kInQSd + 1
198     IF ( kStore.GT.kdiag(ndId) ) THEN
199     _BEGIN_MASTER(myThid)
200 jmc 1.4 WRITE(msgBuf,'(2A,I4,A)') 'DIAGSTATS_FILL: ',
201 jmc 1.1 & 'exceed Nb of levels(=',kdiag(ndId),' ) reserved '
202     CALL PRINT_ERROR( msgBuf , myThid )
203 jmc 1.4 WRITE(msgBuf,'(2A,I6,2A)') 'DIAGSTATS_FILL: ',
204 jmc 1.1 & 'for Diagnostics #', ndId, ' : ', cdiag(ndId)
205     CALL PRINT_ERROR( msgBuf , myThid )
206     WRITE(msgBuf,'(2A,2I4,I3)') 'calling DIAGSTATS_FILL ',
207     I 'with kLev,nLevs,bibjFlg=', kLev,nLevs,bibjFlg
208     CALL PRINT_ERROR( msgBuf , myThid )
209     WRITE(msgBuf,'(2A,I6,A)') 'DIAGSTATS_FILL: ',
210     I '==> trying to store up to ', kStore, ' levels'
211     CALL PRINT_ERROR( msgBuf , myThid )
212     STOP 'ABNORMAL END: S/R DIAGSTATS_FILL'
213     _END_MASTER(myThid)
214     ENDIF
215    
216 jmc 1.5 #ifndef REAL4_IS_SLOW
217     IF ( arrType.EQ.0 .OR. ( arrType.EQ.1 .AND. .NOT.useFract ) ) THEN
218     #endif
219 jmc 1.1 IF ( bibjflg.EQ.0 ) THEN
220     DO bj=myByLo(myThid), myByHi(myThid)
221     DO bi=myBxLo(myThid), myBxHi(myThid)
222     DO k = kFirst,kLast
223     kd = kd0 + ksgn*k
224     km = km0 + ksgn*k
225     CALL DIAGSTATS_LOCAL(
226     U qSdiag(0,0,kd,bi,bj),
227 jmc 1.5 I inpFldRL, fracFldRL,
228 jmc 1.3 I scaleFact, power, useFract, sizF,
229 jmc 1.1 I sizI1,sizI2,sizJ1,sizJ2,nLevs,sizTx,sizTy,
230     I iRun,jRun,k,bi,bj,
231     I km, bi, bj, region2fill,
232     I ndId, gdiag(ndId), myThid )
233     ENDDO
234     ENDDO
235     ENDDO
236     ELSE
237     bi = MIN(biArg,sizTx)
238     bj = MIN(bjArg,sizTy)
239     DO k = kFirst,kLast
240     kd = kd0 + ksgn*k
241     km = km0 + ksgn*k
242     CALL DIAGSTATS_LOCAL(
243     U qSdiag(0,0,kd,biArg,bjArg),
244 jmc 1.5 I inpFldRL, fracFldRL,
245 jmc 1.3 I scaleFact, power, useFract, sizF,
246 jmc 1.1 I sizI1,sizI2,sizJ1,sizJ2,nLevs,sizTx,sizTy,
247     I iRun,jRun,k,bi,bj,
248     I km, biArg, bjArg, region2fill,
249     I ndId, gdiag(ndId), myThid )
250     ENDDO
251     ENDIF
252    
253 jmc 1.5 #ifndef REAL4_IS_SLOW
254     ELSE
255     IF ( bibjflg.EQ.0 ) THEN
256     DO bj=myByLo(myThid), myByHi(myThid)
257     DO bi=myBxLo(myThid), myBxHi(myThid)
258     DO k = kFirst,kLast
259     kd = kd0 + ksgn*k
260     km = km0 + ksgn*k
261     CALL DIAGSTATS_TO_RL(
262     I inpFldRL, fracFldRL, inpFldRS, fracFldRS,
263     O tmpFldRL, tmpFracRL,
264     I arrType, useFract, sizF,
265     I sizI1,sizI2,sizJ1,sizJ2,nLevs,sizTx,sizTy,
266     I iRun,jRun,k,bi,bj, myThid )
267     CALL DIAGSTATS_LOCAL(
268     U qSdiag(0,0,kd,bi,bj),
269     I tmpFldRL, tmpFracRL,
270     I scaleFact, power, useFract, 1,
271     I 1, iRun, 1, jRun, 1, 1, 1,
272     I iRun, jRun, 1, 1, 1,
273     I km, bi, bj, region2fill,
274     I ndId, gdiag(ndId), myThid )
275     ENDDO
276     ENDDO
277     ENDDO
278     ELSE
279     bi = MIN(biArg,sizTx)
280     bj = MIN(bjArg,sizTy)
281     DO k = kFirst,kLast
282     kd = kd0 + ksgn*k
283     km = km0 + ksgn*k
284     CALL DIAGSTATS_TO_RL(
285     I inpFldRL, fracFldRL, inpFldRS, fracFldRS,
286     O tmpFldRL, tmpFracRL,
287     I arrType, useFract, sizF,
288     I sizI1,sizI2,sizJ1,sizJ2,nLevs,sizTx,sizTy,
289     I iRun,jRun,k,bi,bj, myThid )
290     CALL DIAGSTATS_LOCAL(
291     U qSdiag(0,0,kd,biArg,bjArg),
292     I tmpFldRL, tmpFracRL,
293     I scaleFact, power, useFract, 1,
294     I 1, iRun, 1, jRun, 1, 1, 1,
295     I iRun, jRun, 1, 1, 1,
296     I km, biArg, bjArg, region2fill,
297     I ndId, gdiag(ndId), myThid )
298     ENDDO
299     ENDIF
300     ENDIF
301     #endif /* ndef REAL4_IS_SLOW */
302    
303 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
304     c ELSE
305    
306     c ENDIF
307    
308 jmc 1.2 RETURN
309 jmc 1.1 END
310 jmc 1.5
311     #ifndef REAL4_IS_SLOW
312     C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
313    
314     CBOP
315     C !ROUTINE: DIAGSTATS_TO_RL
316     C !INTERFACE:
317     SUBROUTINE DIAGSTATS_TO_RL(
318     I inpFldRL, inpFrcRL, inpFldRS, inpFrcRS,
319     O outFldRL, outFrcRL,
320     I arrType, useFract, sizF,
321     I sizI1,sizI2,sizJ1,sizJ2,sizK,sizTx,sizTy,
322     I iRun,jRun,kIn,biIn,bjIn,
323     I myThid )
324    
325     C !DESCRIPTION:
326     C Do a local copy with conversion to RL type array
327    
328     C !USES:
329     IMPLICIT NONE
330    
331     #include "EEPARAMS.h"
332     #include "SIZE.h"
333    
334     C !INPUT/OUTPUT PARAMETERS:
335     C == Routine Arguments ==
336     C inpFldRL :: input field array to convert (arrType=0,1)
337     C inpFrcRL :: input fraction array to convert (arrType=0,2)
338     C inpFldRS :: input field array to convert (arrType=2,3)
339     C inpFrcRS :: input fraction array to convert (arrType=1,3)
340     C outFldRL :: output field array
341     C outFrcRL :: output fraction array
342     C arrType :: select which array & fraction (RL/RS) to process:
343     C 0: both RL ; 1: fldRL & frcRS ; 2: fldRS,frcRL ; 3: both RS
344     C useFract :: if True, process fraction-weight
345     C sizF :: size of inpFrc array: 3rd dimension
346     C sizI1,sizI2 :: size of inpFld array: 1rst index range (min,max)
347     C sizJ1,sizJ2 :: size of inpFld array: 2nd index range (min,max)
348     C sizK :: size of inpFld array: 3rd dimension
349     C sizTx,sizTy :: size of inpFld array: tile dimensions
350     C iRun,jRun :: range of 1rst & 2nd index
351     C kIn :: level index of inpFld array to process
352     C biIn,bjIn :: tile indices of inpFld array to process
353     C myThid :: my Thread Id number
354     INTEGER sizI1,sizI2,sizJ1,sizJ2
355     INTEGER sizF,sizK,sizTx,sizTy
356     INTEGER iRun, jRun, kIn, biIn, bjIn
357     _RL inpFldRL(sizI1:sizI2,sizJ1:sizJ2,sizK,sizTx,sizTy)
358     _RL inpFrcRL(sizI1:sizI2,sizJ1:sizJ2,sizF,sizTx,sizTy)
359     _RS inpFldRS(sizI1:sizI2,sizJ1:sizJ2,sizK,sizTx,sizTy)
360     _RS inpFrcRS(sizI1:sizI2,sizJ1:sizJ2,sizF,sizTx,sizTy)
361     _RL outFldRL(1:iRun,1:jRun)
362     _RL outFrcRL(1:iRun,1:jRun)
363     INTEGER arrType
364     LOGICAL useFract
365     INTEGER myThid
366     CEOP
367    
368     C !LOCAL VARIABLES:
369     C i,j :: loop indices
370     INTEGER i, j, kFr
371    
372     IF ( arrType.LE.1 ) THEN
373     DO j=1,jRun
374     DO i=1,iRun
375     outFldRL(i,j) = inpFldRL(i,j,kIn,biIn,bjIn)
376     ENDDO
377     ENDDO
378     ELSE
379     DO j=1,jRun
380     DO i=1,iRun
381     outFldRL(i,j) = inpFldRS(i,j,kIn,biIn,bjIn)
382     ENDDO
383     ENDDO
384     ENDIF
385    
386     IF ( useFract ) THEN
387     kFr = MIN(kIn,sizF)
388     IF ( arrType.EQ.0 .OR. arrType.EQ.2 ) THEN
389     DO j=1,jRun
390     DO i=1,iRun
391     outFrcRL(i,j) = inpFrcRL(i,j,kFr,biIn,bjIn)
392     ENDDO
393     ENDDO
394     ELSE
395     DO j=1,jRun
396     DO i=1,iRun
397     outFrcRL(i,j) = inpFrcRS(i,j,kFr,biIn,bjIn)
398     ENDDO
399     ENDDO
400     ENDIF
401     ENDIF
402    
403     RETURN
404     END
405     #endif /* ndef REAL4_IS_SLOW */

  ViewVC Help
Powered by ViewVC 1.1.22