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

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

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


Revision 1.56 - (show 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 C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.55 2011/06/22 19:09:40 jmc Exp $
2 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 SUBROUTINE DIAGNOSTICS_OUT(
12 I listId,
13 I myTime,
14 I myIter,
15 I myThid )
16
17 C !DESCRIPTION:
18 C Write output for diagnostics fields.
19
20 C !USES:
21 IMPLICIT NONE
22 #include "SIZE.h"
23 #include "EEPARAMS.h"
24 #include "PARAMS.h"
25 #include "GRID.h"
26 #include "DIAGNOSTICS_SIZE.h"
27 #include "DIAGNOSTICS.h"
28
29 INTEGER NrMax
30 PARAMETER( NrMax = numLevels )
31
32 C !INPUT PARAMETERS:
33 C listId :: Diagnostics list number being written
34 C myIter :: current iteration number
35 C myTime :: current time of simulation (s)
36 C myThid :: my Thread Id number
37 _RL myTime
38 INTEGER listId, myIter, myThid
39 CEOP
40
41 C !FUNCTIONS:
42 INTEGER ILNBLNK
43 EXTERNAL ILNBLNK
44 #ifdef ALLOW_FIZHI
45 _RL getcon
46 EXTERNAL getcon
47 #endif
48
49 C !LOCAL VARIABLES:
50 C i,j,k :: loop indices
51 C bi,bj :: tile indices
52 C lm :: loop index (averageCycle)
53 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 C nLevOutp :: number of levels to write in output file
59 C
60 C-- COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
61 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 _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
66 _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
67
68 INTEGER i, j, k, lm
69 INTEGER bi, bj
70 INTEGER md, ndId, nn, ip, im
71 INTEGER mate, mDbl, mVec
72 CHARACTER*10 gcode
73 _RL undefRL
74 INTEGER nLevOutp, kLev
75
76 INTEGER iLen
77 INTEGER ioUnit
78 CHARACTER*(MAX_LEN_FNAM) fn
79 CHARACTER*(MAX_LEN_MBUF) suff
80 CHARACTER*(MAX_LEN_MBUF) msgBuf
81 INTEGER prec, nRec, nTimRec
82 _RL timeRec(2)
83 _RL tmpLoc
84 #ifdef ALLOW_MDSIO
85 LOGICAL glf
86 #endif
87 #ifdef ALLOW_MNC
88 CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
89 LOGICAL missingValFillsMask
90 #endif /* ALLOW_MNC */
91
92 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
93
94 C--- set file properties
95 ioUnit= standardMessageUnit
96 undefRL = UNSET_RL
97 #ifdef ALLOW_FIZHI
98 IF ( useFIZHI ) undefRL = getcon('UNDEF')
99 #endif
100 IF ( misvalFlt(listId).NE.UNSET_RL ) undefRL = misvalFlt(listId)
101
102 WRITE(suff,'(I10.10)') myIter
103 iLen = ILNBLNK(fnames(listId))
104 WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
105 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
109 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 IF ( timeRec(1).LT.0. ) THEN
130 tmpLoc = FLOAT(i)
131 IF ( timeRec(1).NE.tmpLoc ) i = i - 1
132 ENDIF
133 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 tmpLoc = FLOAT(i) - 0.5 _d 0
142 IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
143 timeRec(1) = baseTime + deltaTClock*FLOAT(i)
144 c if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
145 ENDIF
146 C-- Convert time to iteration number (debug)
147 c DO i=1,nTimRec
148 c timeRec(i) = timeRec(i)/deltaTClock
149 c ENDDO
150
151 C-- Place the loop on lm (= averagePeriod) outside the loop on md (= field):
152 DO lm=1,averageCycle(listId)
153
154 #ifdef ALLOW_MNC
155 IF (useMNC .AND. diag_mnc) THEN
156 CALL DIAGNOSTICS_MNC_SET(
157 I nLevOutp, listId, lm,
158 O diag_mnc_bn, missingValFillsMask,
159 I undefRL, myTime, myIter, myThid )
160 ENDIF
161 #endif /* ALLOW_MNC */
162
163 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
164
165 DO md = 1,nfields(listId)
166 ndId = jdiag(md,listId)
167 gcode = gdiag(ndId)(1:10)
168 mate = 0
169 mVec = 0
170 mDbl = 0
171 IF ( gcode(5:5).EQ.'C' ) THEN
172 C- Check for Mate of a Counter Diagnostic
173 mate = hdiag(ndId)
174 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 ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
182 C- Check for Mate of a Vector Diagnostic
183 mVec = hdiag(ndId)
184 ENDIF
185 IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
186 C-- Start processing 1 Fld :
187
188 ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
189 im = mdiag(md,listId)
190 IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
191 IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
192 IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
193
194 IF ( ndiag(ip,1,1).EQ.0 ) THEN
195 C- Empty diagnostics case :
196
197 _BEGIN_MASTER( myThid )
198 WRITE(msgBuf,'(A,I10)')
199 & '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
200 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
201 & SQUEEZE_RIGHT, myThid)
202 WRITE(msgBuf,'(A,I6,3A,I4,2A)')
203 & '- WARNING - diag.#',ndId, ' : ',flds(md,listId),
204 & ' (#',md,' ) in outp.Stream: ',fnames(listId)
205 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
206 & SQUEEZE_RIGHT, myThid)
207 IF ( averageCycle(listId).GT.1 ) THEN
208 WRITE(msgBuf,'(A,2(I3,A))')
209 & '- WARNING - has not been filled (ndiag(lm=',lm,')=',
210 & ndiag(ip,1,1), ' )'
211 ELSE
212 WRITE(msgBuf,'(A,2(I3,A))')
213 & '- WARNING - has not been filled (ndiag=',
214 & ndiag(ip,1,1), ' )'
215 ENDIF
216 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
217 & SQUEEZE_RIGHT, myThid)
218 WRITE(msgBuf,'(A)')
219 & 'WARNING DIAGNOSTICS_OUT => write ZEROS instead'
220 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
221 & SQUEEZE_RIGHT, myThid)
222 _END_MASTER( myThid )
223 DO bj = myByLo(myThid), myByHi(myThid)
224 DO bi = myBxLo(myThid), myBxHi(myThid)
225 DO k = 1,nLevOutp
226 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 ELSE
236 C- diagnostics is not empty :
237
238 IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
239 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 & ' Computing Diagnostic # ', ndId, ' ', cdiag(ndId),
255 & ' Counter:',ndiag(ip,1,1),' Parms: ',gdiag(ndId)
256 ENDIF
257 IF ( mate.GT.0 ) THEN
258 WRITE(ioUnit,'(3A,I6,2A)')
259 & ' use Counter Mate for ', cdiag(ndId),
260 & ' Diagnostic # ',mate, ' ', cdiag(mate)
261 ELSEIF ( mVec.GT.0 ) THEN
262 IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN
263 WRITE(ioUnit,'(3A,I6,3A)')
264 & ' Vector Mate for ', cdiag(ndId),
265 & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
266 & ' exists '
267 ELSE
268 WRITE(ioUnit,'(3A,I6,3A)')
269 & ' Vector Mate for ', cdiag(ndId),
270 & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
271 & ' not enabled'
272 ENDIF
273 ENDIF
274 ENDIF
275
276 IF ( fflags(listId)(2:2).EQ.' ' ) THEN
277 C- get only selected levels:
278 DO bj = myByLo(myThid), myByHi(myThid)
279 DO bi = myBxLo(myThid), myBxHi(myThid)
280 DO k = 1,nlevels(listId)
281 kLev = NINT(levs(k,listId))
282 CALL DIAGNOSTICS_GET_DIAG(
283 I kLev, undefRL,
284 O qtmp1(1-OLx,1-OLy,k,bi,bj),
285 I ndId, mate, ip, im, bi, bj, myThid )
286 ENDDO
287 ENDDO
288 ENDDO
289 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 ELSE
303 C- get all the levels (for vertical post-processing)
304 DO bj = myByLo(myThid), myByHi(myThid)
305 DO bi = myBxLo(myThid), myBxHi(myThid)
306 CALL DIAGNOSTICS_GET_DIAG(
307 I 0, undefRL,
308 O qtmp1(1-OLx,1-OLy,1,bi,bj),
309 I ndId, mate, ip, im, bi, bj, myThid )
310 ENDDO
311 ENDDO
312 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 ENDIF
325
326 C-----------------------------------------------------------------------
327 C-- Apply specific post-processing (e.g., interpolate) before output
328 C-----------------------------------------------------------------------
329 IF ( fflags(listId)(2:2).EQ.'P' ) THEN
330 C- Do vertical interpolation:
331 IF ( fluidIsAir ) THEN
332 C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
333 CALL DIAGNOSTICS_INTERP_VERT(
334 I listId, md, ndId, ip, im, lm,
335 U qtmp1, qtmp2,
336 I undefRL, myTime, myIter, myThid )
337 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 I undefRL, myTime, myIter, myThid )
350 ENDIF
351 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
368 C-- End of empty diag / not-empty diag block
369 ENDIF
370
371 C-- Ready to write field "md", element "lm" in averageCycle(listId)
372
373 C- write to binary file, using MDSIO pkg:
374 IF ( diag_mdsio ) THEN
375 c nRec = lm + (md-1)*averageCycle(listId)
376 nRec = md + (lm-1)*nfields(listId)
377 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 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 I NrMax, 1, nLevOutp,
386 I qtmp1, -nRec, myIter, myThid )
387 ENDIF
388
389 #ifdef ALLOW_MNC
390 IF (useMNC .AND. diag_mnc) THEN
391 CALL DIAGNOSTICS_MNC_OUT(
392 I NrMax, nLevOutp, listId, ndId, mate,
393 I diag_mnc_bn,
394 I missingValFillsMask, undefRL,
395 I qtmp1,
396 I myTime, myIter, myThid )
397 ENDIF
398 #endif /* ALLOW_MNC */
399
400 C-- end of Processing Fld # md
401 ENDIF
402 ENDDO
403
404 C-- end loop on lm counter (= averagePeriod)
405 ENDDO
406
407 #ifdef ALLOW_MDSIO
408 IF (diag_mdsio) THEN
409 C- Note: temporary: since it is a pain to add more arguments to
410 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 glf = globalFiles
413 nRec = averageCycle(listId)*nfields(listId)
414 CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
415 & 0, 0, nLevOutp, ' ',
416 & nfields(listId), flds(1,listId), nTimRec, timeRec,
417 & nRec, myIter, myThid)
418 ENDIF
419 #endif /* ALLOW_MDSIO */
420
421 RETURN
422 END
423
424 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

  ViewVC Help
Powered by ViewVC 1.1.22