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

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

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


Revision 1.56 - (hide annotations) (download)
Thu Jun 23 15:29:01 2011 UTC (12 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.55: +8 -10 lines
Extend the use of "missing_value" setting from data.diagnostics
for the case of counter-diagnostics (where fraction-weight is zero).
This applies both to MNC and MDS output file.

1 jmc 1.56 C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.55 2011/06/22 19:09:40 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: DIAGNOSTICS_OUT
9    
10     C !INTERFACE:
11 jmc 1.47 SUBROUTINE DIAGNOSTICS_OUT(
12 jmc 1.15 I listId,
13 jmc 1.50 I myTime,
14 jmc 1.1 I myIter,
15     I myThid )
16    
17     C !DESCRIPTION:
18     C Write output for diagnostics fields.
19 jmc 1.15
20 jmc 1.1 C !USES:
21 jmc 1.3 IMPLICIT NONE
22 jmc 1.1 #include "SIZE.h"
23     #include "EEPARAMS.h"
24     #include "PARAMS.h"
25 edhill 1.7 #include "GRID.h"
26 jmc 1.3 #include "DIAGNOSTICS_SIZE.h"
27     #include "DIAGNOSTICS.h"
28 jmc 1.1
29 jmc 1.30 INTEGER NrMax
30 jmc 1.33 PARAMETER( NrMax = numLevels )
31 jmc 1.1
32     C !INPUT PARAMETERS:
33 jmc 1.15 C listId :: Diagnostics list number being written
34 jmc 1.3 C myIter :: current iteration number
35 jmc 1.15 C myTime :: current time of simulation (s)
36 jmc 1.3 C myThid :: my Thread Id number
37 edhill 1.14 _RL myTime
38 jmc 1.15 INTEGER listId, myIter, myThid
39 jmc 1.1 CEOP
40    
41 jmc 1.40 C !FUNCTIONS:
42     INTEGER ILNBLNK
43     EXTERNAL ILNBLNK
44     #ifdef ALLOW_FIZHI
45     _RL getcon
46     EXTERNAL getcon
47     #endif
48    
49 jmc 1.3 C !LOCAL VARIABLES:
50 jmc 1.15 C i,j,k :: loop indices
51 jmc 1.47 C bi,bj :: tile indices
52 jmc 1.29 C lm :: loop index (averageCycle)
53 jmc 1.15 C md :: field number in the list "listId".
54     C ndId :: diagnostics Id number (in available diagnostics list)
55     C mate :: counter mate Id number (in available diagnostics list)
56     C ip :: diagnostics pointer to storage array
57     C im :: counter-mate pointer to storage array
58 jmc 1.47 C nLevOutp :: number of levels to write in output file
59 jmc 1.32 C
60     C-- COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
61 jmc 1.52 C qtmp1 :: temporary array; used to store a copy of diag. output field.
62     C qtmp2 :: temporary array; used to store a copy of a 2nd diag. field.
63     C- Note: local common block no longer needed.
64     c COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
65 jmc 1.32 _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
66 jmc 1.52 _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
67 jmc 1.32
68 jmc 1.43 INTEGER i, j, k, lm
69 jmc 1.15 INTEGER bi, bj
70 jmc 1.54 INTEGER md, ndId, nn, ip, im
71     INTEGER mate, mDbl, mVec
72 jmc 1.35 CHARACTER*10 gcode
73 jmc 1.52 _RL undefRL
74     INTEGER nLevOutp, kLev
75    
76 jmc 1.50 INTEGER iLen
77 jmc 1.6 INTEGER ioUnit
78 jmc 1.11 CHARACTER*(MAX_LEN_FNAM) fn
79 jmc 1.1 CHARACTER*(MAX_LEN_MBUF) suff
80 jmc 1.3 CHARACTER*(MAX_LEN_MBUF) msgBuf
81 jmc 1.44 INTEGER prec, nRec, nTimRec
82     _RL timeRec(2)
83 jmc 1.52 _RL tmpLoc
84 jmc 1.29 #ifdef ALLOW_MDSIO
85 jmc 1.3 LOGICAL glf
86 jmc 1.29 #endif
87 jmc 1.1 #ifdef ALLOW_MNC
88     CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
89 jmc 1.56 LOGICAL missingValFillsMask
90 jmc 1.1 #endif /* ALLOW_MNC */
91    
92 jmc 1.3 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
93    
94 jmc 1.44 C--- set file properties
95 jmc 1.6 ioUnit= standardMessageUnit
96 jmc 1.52 undefRL = UNSET_RL
97 jmc 1.40 #ifdef ALLOW_FIZHI
98 jmc 1.52 IF ( useFIZHI ) undefRL = getcon('UNDEF')
99 jmc 1.40 #endif
100 jmc 1.56 IF ( misvalFlt(listId).NE.UNSET_RL ) undefRL = misvalFlt(listId)
101    
102 jmc 1.1 WRITE(suff,'(I10.10)') myIter
103 jmc 1.50 iLen = ILNBLNK(fnames(listId))
104     WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
105 jmc 1.47 C- for now, if integrate vertically, output field has just 1 level:
106     nLevOutp = nlevels(listId)
107     IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
108 jmc 1.1
109 jmc 1.44 C-- Set time information:
110     IF ( freq(listId).LT.0. ) THEN
111     C- Snap-shot: store a unique time (which is consistent with State-Var timing)
112     nTimRec = 1
113     timeRec(1) = myTime
114     ELSE
115     C- Time-average: store the 2 edges of the time-averaging interval.
116     C this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
117     C tendencies) timing. For State-Var, this is shifted by + halt time-step.
118     nTimRec = 2
119    
120     C- end of time-averaging interval:
121     timeRec(2) = myTime
122    
123     C- begining of time-averaging interval:
124     c timeRec(1) = myTime - freq(listId)
125     C a) find the time of the previous multiple of output freq:
126     timeRec(1) = myTime-deltaTClock*0.5 _d 0
127     timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
128     i = INT( timeRec(1) )
129 jmc 1.46 IF ( timeRec(1).LT.0. ) THEN
130 jmc 1.52 tmpLoc = FLOAT(i)
131     IF ( timeRec(1).NE.tmpLoc ) i = i - 1
132 jmc 1.46 ENDIF
133 jmc 1.44 timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
134     c if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
135     timeRec(1) = MAX( timeRec(1), startTime )
136    
137     C b) round off to nearest multiple of time-step:
138     timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
139     i = NINT( timeRec(1) )
140     C if just half way, NINT will return the next time-step: correct this
141 jmc 1.52 tmpLoc = FLOAT(i) - 0.5 _d 0
142     IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
143 jmc 1.44 timeRec(1) = baseTime + deltaTClock*FLOAT(i)
144     c if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
145     ENDIF
146 jmc 1.46 C-- Convert time to iteration number (debug)
147     c DO i=1,nTimRec
148     c timeRec(i) = timeRec(i)/deltaTClock
149     c ENDDO
150 jmc 1.44
151 jmc 1.55 C-- Place the loop on lm (= averagePeriod) outside the loop on md (= field):
152     DO lm=1,averageCycle(listId)
153    
154 jmc 1.1 #ifdef ALLOW_MNC
155 jmc 1.41 IF (useMNC .AND. diag_mnc) THEN
156 jmc 1.50 CALL DIAGNOSTICS_MNC_SET(
157 jmc 1.55 I nLevOutp, listId, lm,
158 jmc 1.56 O diag_mnc_bn, missingValFillsMask,
159     I undefRL, myTime, myIter, myThid )
160 jmc 1.41 ENDIF
161 jmc 1.1 #endif /* ALLOW_MNC */
162    
163 jmc 1.29 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
164    
165 jmc 1.41 DO md = 1,nfields(listId)
166 jmc 1.15 ndId = jdiag(md,listId)
167 jmc 1.35 gcode = gdiag(ndId)(1:10)
168 jmc 1.29 mate = 0
169     mVec = 0
170 jmc 1.54 mDbl = 0
171 jmc 1.35 IF ( gcode(5:5).EQ.'C' ) THEN
172 jmc 1.29 C- Check for Mate of a Counter Diagnostic
173 jmc 1.35 mate = hdiag(ndId)
174 jmc 1.54 ELSEIF ( gcode(5:5).EQ.'P' ) THEN
175     C- Also load the mate (if stored) for Post-Processing
176     nn = ndId
177     DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
178     nn = hdiag(nn)
179     ENDDO
180     IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
181 jmc 1.35 ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
182 jmc 1.29 C- Check for Mate of a Vector Diagnostic
183 jmc 1.36 mVec = hdiag(ndId)
184 jmc 1.29 ENDIF
185 jmc 1.35 IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
186 jmc 1.3 C-- Start processing 1 Fld :
187    
188 jmc 1.29 ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
189 jmc 1.15 im = mdiag(md,listId)
190 jmc 1.29 IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
191 jmc 1.54 IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
192 jmc 1.29 IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
193    
194 jmc 1.54 IF ( ndiag(ip,1,1).EQ.0 ) THEN
195 jmc 1.3 C- Empty diagnostics case :
196    
197     _BEGIN_MASTER( myThid )
198     WRITE(msgBuf,'(A,I10)')
199     & '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
200 jmc 1.15 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
201 jmc 1.3 & SQUEEZE_RIGHT, myThid)
202 jmc 1.35 WRITE(msgBuf,'(A,I6,3A,I4,2A)')
203 jmc 1.15 & '- WARNING - diag.#',ndId, ' : ',flds(md,listId),
204     & ' (#',md,' ) in outp.Stream: ',fnames(listId)
205     CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
206 jmc 1.3 & SQUEEZE_RIGHT, myThid)
207 jmc 1.29 IF ( averageCycle(listId).GT.1 ) THEN
208 jmc 1.35 WRITE(msgBuf,'(A,2(I3,A))')
209 jmc 1.29 & '- WARNING - has not been filled (ndiag(lm=',lm,')=',
210     & ndiag(ip,1,1), ' )'
211     ELSE
212 jmc 1.35 WRITE(msgBuf,'(A,2(I3,A))')
213 jmc 1.29 & '- WARNING - has not been filled (ndiag=',
214     & ndiag(ip,1,1), ' )'
215     ENDIF
216 jmc 1.15 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
217 jmc 1.3 & SQUEEZE_RIGHT, myThid)
218     WRITE(msgBuf,'(A)')
219     & 'WARNING DIAGNOSTICS_OUT => write ZEROS instead'
220 jmc 1.15 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
221 jmc 1.3 & SQUEEZE_RIGHT, myThid)
222     _END_MASTER( myThid )
223     DO bj = myByLo(myThid), myByHi(myThid)
224     DO bi = myBxLo(myThid), myBxHi(myThid)
225 jmc 1.47 DO k = 1,nLevOutp
226 jmc 1.3 DO j = 1-OLy,sNy+OLy
227     DO i = 1-OLx,sNx+OLx
228     qtmp1(i,j,k,bi,bj) = 0. _d 0
229     ENDDO
230     ENDDO
231     ENDDO
232     ENDDO
233     ENDDO
234    
235 jmc 1.54 ELSE
236 jmc 1.3 C- diagnostics is not empty :
237    
238 jmc 1.49 IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
239 jmc 1.54 IF ( gcode(5:5).EQ.'P' ) THEN
240     WRITE(ioUnit,'(A,I6,7A,I8,2A)')
241     & ' Post-Processing Diag # ', ndId, ' ', cdiag(ndId),
242     & ' Parms: ',gdiag(ndId)
243     IF ( mDbl.EQ.0 ) THEN
244     WRITE(ioUnit,'(2(3A,I6,A,I8))') ' from diag: ',
245     & cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
246     ELSE
247     WRITE(ioUnit,'(2(3A,I6,A,I8))') ' from diag: ',
248     & cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
249     & ' and diag: ',
250     & cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
251     ENDIF
252     ELSE
253     WRITE(ioUnit,'(A,I6,3A,I8,2A)')
254 jmc 1.15 & ' Computing Diagnostic # ', ndId, ' ', cdiag(ndId),
255     & ' Counter:',ndiag(ip,1,1),' Parms: ',gdiag(ndId)
256 jmc 1.54 ENDIF
257 jmc 1.29 IF ( mate.GT.0 ) THEN
258 jmc 1.35 WRITE(ioUnit,'(3A,I6,2A)')
259 jmc 1.15 & ' use Counter Mate for ', cdiag(ndId),
260     & ' Diagnostic # ',mate, ' ', cdiag(mate)
261 jmc 1.29 ELSEIF ( mVec.GT.0 ) THEN
262 jmc 1.15 IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN
263 jmc 1.35 WRITE(ioUnit,'(3A,I6,3A)')
264 jmc 1.15 & ' Vector Mate for ', cdiag(ndId),
265     & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
266     & ' exists '
267 jmc 1.3 ELSE
268 jmc 1.35 WRITE(ioUnit,'(3A,I6,3A)')
269 jmc 1.15 & ' Vector Mate for ', cdiag(ndId),
270     & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
271     & ' not enabled'
272 jmc 1.3 ENDIF
273     ENDIF
274 jmc 1.6 ENDIF
275 jmc 1.3
276 jmc 1.52 IF ( fflags(listId)(2:2).EQ.' ' ) THEN
277     C- get only selected levels:
278 jmc 1.30 DO bj = myByLo(myThid), myByHi(myThid)
279     DO bi = myBxLo(myThid), myBxHi(myThid)
280 jmc 1.52 DO k = 1,nlevels(listId)
281     kLev = NINT(levs(k,listId))
282     CALL DIAGNOSTICS_GET_DIAG(
283     I kLev, undefRL,
284 jmc 1.30 O qtmp1(1-OLx,1-OLy,k,bi,bj),
285 jmc 1.54 I ndId, mate, ip, im, bi, bj, myThid )
286 jmc 1.30 ENDDO
287     ENDDO
288     ENDDO
289 jmc 1.54 IF ( mDbl.GT.0 ) THEN
290     DO bj = myByLo(myThid), myByHi(myThid)
291     DO bi = myBxLo(myThid), myBxHi(myThid)
292     DO k = 1,nlevels(listId)
293     kLev = NINT(levs(k,listId))
294     CALL DIAGNOSTICS_GET_DIAG(
295     I kLev, undefRL,
296     O qtmp2(1-OLx,1-OLy,k,bi,bj),
297     I mDbl, 0, im, 0, bi, bj, myThid )
298     ENDDO
299     ENDDO
300     ENDDO
301     ENDIF
302 jmc 1.30 ELSE
303 jmc 1.52 C- get all the levels (for vertical post-processing)
304 jmc 1.30 DO bj = myByLo(myThid), myByHi(myThid)
305     DO bi = myBxLo(myThid), myBxHi(myThid)
306 jmc 1.52 CALL DIAGNOSTICS_GET_DIAG(
307     I 0, undefRL,
308     O qtmp1(1-OLx,1-OLy,1,bi,bj),
309 jmc 1.54 I ndId, mate, ip, im, bi, bj, myThid )
310 jmc 1.30 ENDDO
311 jmc 1.3 ENDDO
312 jmc 1.54 IF ( mDbl.GT.0 ) THEN
313     DO bj = myByLo(myThid), myByHi(myThid)
314     DO bi = myBxLo(myThid), myBxHi(myThid)
315     DO k = 1,nlevels(listId)
316     CALL DIAGNOSTICS_GET_DIAG(
317     I 0, undefRL,
318     O qtmp2(1-OLx,1-OLy,k,bi,bj),
319     I mDbl, 0, im, 0, bi, bj, myThid )
320     ENDDO
321     ENDDO
322     ENDDO
323     ENDIF
324 jmc 1.30 ENDIF
325 jmc 1.1
326 molod 1.17 C-----------------------------------------------------------------------
327 jmc 1.47 C-- Apply specific post-processing (e.g., interpolate) before output
328 molod 1.17 C-----------------------------------------------------------------------
329 jmc 1.47 IF ( fflags(listId)(2:2).EQ.'P' ) THEN
330     C- Do vertical interpolation:
331     IF ( fluidIsAir ) THEN
332 jmc 1.29 C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
333 jmc 1.47 CALL DIAGNOSTICS_INTERP_VERT(
334     I listId, md, ndId, ip, im, lm,
335 jmc 1.52 U qtmp1, qtmp2,
336     I undefRL, myTime, myIter, myThid )
337 jmc 1.47 ELSE
338     WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
339     & 'INTERP_VERT not allowed in this config'
340     CALL PRINT_ERROR( msgBuf , myThid )
341     STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
342     ENDIF
343     ENDIF
344     IF ( fflags(listId)(2:2).EQ.'I' ) THEN
345     C- Integrate vertically: for now, output field has just 1 level:
346     CALL DIAGNOSTICS_SUM_LEVELS(
347     I listId, md, ndId, ip, im, lm,
348     U qtmp1,
349 jmc 1.52 I undefRL, myTime, myIter, myThid )
350 jmc 1.47 ENDIF
351 jmc 1.54 IF ( gcode(5:5).EQ.'P' ) THEN
352     C- Do Post-Processing:
353     IF ( flds(md,listId).EQ.'PhiVEL '
354     c & .OR. flds(md,listId).EQ.'PsiVEL '
355     & ) THEN
356     CALL DIAGNOSTICS_CALC_PHIVEL(
357     I listId, md, ndId, ip, im, lm,
358     U qtmp1, qtmp2,
359     I myTime, myIter, myThid )
360     ELSE
361     WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
362     & 'unknown Processing method for diag="',cdiag(ndId),'"'
363     CALL PRINT_ERROR( msgBuf , myThid )
364     STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
365     ENDIF
366     ENDIF
367 jmc 1.47
368     C-- End of empty diag / not-empty diag block
369 jmc 1.29 ENDIF
370 molod 1.17
371 jmc 1.47 C-- Ready to write field "md", element "lm" in averageCycle(listId)
372 jmc 1.31
373     C- write to binary file, using MDSIO pkg:
374 jmc 1.34 IF ( diag_mdsio ) THEN
375 jmc 1.55 c nRec = lm + (md-1)*averageCycle(listId)
376     nRec = md + (lm-1)*nfields(listId)
377 jmc 1.30 C default precision for output files
378     prec = writeBinaryPrec
379     C fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
380     IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
381     IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
382 jmc 1.34 C a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R
383     CALL WRITE_REC_LEV_RL(
384     I fn, prec,
385 jmc 1.47 I NrMax, 1, nLevOutp,
386 jmc 1.34 I qtmp1, -nRec, myIter, myThid )
387 jmc 1.3 ENDIF
388 jmc 1.1
389     #ifdef ALLOW_MNC
390 jmc 1.3 IF (useMNC .AND. diag_mnc) THEN
391 jmc 1.50 CALL DIAGNOSTICS_MNC_OUT(
392 jmc 1.56 I NrMax, nLevOutp, listId, ndId, mate,
393 jmc 1.50 I diag_mnc_bn,
394 jmc 1.56 I missingValFillsMask, undefRL,
395 jmc 1.50 I qtmp1,
396     I myTime, myIter, myThid )
397 jmc 1.3 ENDIF
398 jmc 1.1 #endif /* ALLOW_MNC */
399    
400 jmc 1.15 C-- end of Processing Fld # md
401 jmc 1.3 ENDIF
402 jmc 1.41 ENDDO
403    
404 jmc 1.55 C-- end loop on lm counter (= averagePeriod)
405 jmc 1.3 ENDDO
406 jmc 1.1
407 jmc 1.31 #ifdef ALLOW_MDSIO
408     IF (diag_mdsio) THEN
409 jmc 1.48 C- Note: temporary: since it is a pain to add more arguments to
410 jmc 1.31 C all MDSIO S/R, uses instead this specific S/R to write only
411     C meta files but with more informations in it.
412 jmc 1.34 glf = globalFiles
413 jmc 1.55 nRec = averageCycle(listId)*nfields(listId)
414 jmc 1.31 CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
415 jmc 1.47 & 0, 0, nLevOutp, ' ',
416 jmc 1.44 & nfields(listId), flds(1,listId), nTimRec, timeRec,
417 jmc 1.31 & nRec, myIter, myThid)
418     ENDIF
419     #endif /* ALLOW_MDSIO */
420    
421 jmc 1.15 RETURN
422 jmc 1.3 END
423 jmc 1.15
424 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

  ViewVC Help
Powered by ViewVC 1.1.22