/[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.60 - (show annotations) (download)
Sun Jan 13 22:46:38 2013 UTC (11 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint64c
Changes since 1.59: +3 -2 lines
- add missing value argument to S/R MDS_WR_METAFILES argument list

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

  ViewVC Help
Powered by ViewVC 1.1.22