/[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.55 - (show annotations) (download)
Wed Jun 22 19:09:40 2011 UTC (12 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.54: +9 -28 lines
Change the ordering of fields and time-period in MDS output file
 when using periodic averaging: now writes one time record with the full set
 of fields for this time period, then the next time reccord (it used to be
 one field with the full set of time records then the next field)
Note: this is now similar to NetCDF (MNC) output file.

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

  ViewVC Help
Powered by ViewVC 1.1.22