/[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.55 - (hide 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 jmc 1.55 C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.54 2011/06/21 18:00:48 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.50 LOGICAL useMissingValue
90     REAL*8 misValLoc
91 jmc 1.1 #endif /* ALLOW_MNC */
92    
93 jmc 1.3 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
94    
95 jmc 1.44 C--- set file properties
96 jmc 1.6 ioUnit= standardMessageUnit
97 jmc 1.52 undefRL = UNSET_RL
98 jmc 1.40 #ifdef ALLOW_FIZHI
99 jmc 1.52 IF ( useFIZHI ) undefRL = getcon('UNDEF')
100 jmc 1.40 #endif
101 jmc 1.1 WRITE(suff,'(I10.10)') myIter
102 jmc 1.50 iLen = ILNBLNK(fnames(listId))
103     WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
104 jmc 1.47 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 jmc 1.1
108 jmc 1.44 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 jmc 1.46 IF ( timeRec(1).LT.0. ) THEN
129 jmc 1.52 tmpLoc = FLOAT(i)
130     IF ( timeRec(1).NE.tmpLoc ) i = i - 1
131 jmc 1.46 ENDIF
132 jmc 1.44 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 jmc 1.52 tmpLoc = FLOAT(i) - 0.5 _d 0
141     IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
142 jmc 1.44 timeRec(1) = baseTime + deltaTClock*FLOAT(i)
143     c if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
144     ENDIF
145 jmc 1.46 C-- Convert time to iteration number (debug)
146     c DO i=1,nTimRec
147     c timeRec(i) = timeRec(i)/deltaTClock
148     c ENDDO
149 jmc 1.44
150 jmc 1.55 C-- Place the loop on lm (= averagePeriod) outside the loop on md (= field):
151     DO lm=1,averageCycle(listId)
152    
153 jmc 1.1 #ifdef ALLOW_MNC
154 jmc 1.41 IF (useMNC .AND. diag_mnc) THEN
155 jmc 1.52 misValLoc = undefRL
156 jmc 1.51 IF ( misvalFlt(listId).NE.UNSET_RL )
157     & misValLoc = misvalFlt(listId)
158 jmc 1.50 CALL DIAGNOSTICS_MNC_SET(
159 jmc 1.55 I nLevOutp, listId, lm,
160 jmc 1.51 O diag_mnc_bn, useMissingValue,
161     I misValLoc, myTime, myIter, myThid )
162 jmc 1.41 ENDIF
163 jmc 1.1 #endif /* ALLOW_MNC */
164    
165 jmc 1.29 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
166    
167 jmc 1.41 DO md = 1,nfields(listId)
168 jmc 1.15 ndId = jdiag(md,listId)
169 jmc 1.35 gcode = gdiag(ndId)(1:10)
170 jmc 1.29 mate = 0
171     mVec = 0
172 jmc 1.54 mDbl = 0
173 jmc 1.35 IF ( gcode(5:5).EQ.'C' ) THEN
174 jmc 1.29 C- Check for Mate of a Counter Diagnostic
175 jmc 1.35 mate = hdiag(ndId)
176 jmc 1.54 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 jmc 1.35 ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
184 jmc 1.29 C- Check for Mate of a Vector Diagnostic
185 jmc 1.36 mVec = hdiag(ndId)
186 jmc 1.29 ENDIF
187 jmc 1.35 IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
188 jmc 1.3 C-- Start processing 1 Fld :
189    
190 jmc 1.29 ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
191 jmc 1.15 im = mdiag(md,listId)
192 jmc 1.29 IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
193 jmc 1.54 IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
194 jmc 1.29 IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
195    
196 jmc 1.54 IF ( ndiag(ip,1,1).EQ.0 ) THEN
197 jmc 1.3 C- Empty diagnostics case :
198    
199     _BEGIN_MASTER( myThid )
200     WRITE(msgBuf,'(A,I10)')
201     & '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
202 jmc 1.15 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
203 jmc 1.3 & SQUEEZE_RIGHT, myThid)
204 jmc 1.35 WRITE(msgBuf,'(A,I6,3A,I4,2A)')
205 jmc 1.15 & '- WARNING - diag.#',ndId, ' : ',flds(md,listId),
206     & ' (#',md,' ) in outp.Stream: ',fnames(listId)
207     CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
208 jmc 1.3 & SQUEEZE_RIGHT, myThid)
209 jmc 1.29 IF ( averageCycle(listId).GT.1 ) THEN
210 jmc 1.35 WRITE(msgBuf,'(A,2(I3,A))')
211 jmc 1.29 & '- WARNING - has not been filled (ndiag(lm=',lm,')=',
212     & ndiag(ip,1,1), ' )'
213     ELSE
214 jmc 1.35 WRITE(msgBuf,'(A,2(I3,A))')
215 jmc 1.29 & '- WARNING - has not been filled (ndiag=',
216     & ndiag(ip,1,1), ' )'
217     ENDIF
218 jmc 1.15 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
219 jmc 1.3 & SQUEEZE_RIGHT, myThid)
220     WRITE(msgBuf,'(A)')
221     & 'WARNING DIAGNOSTICS_OUT => write ZEROS instead'
222 jmc 1.15 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
223 jmc 1.3 & SQUEEZE_RIGHT, myThid)
224     _END_MASTER( myThid )
225     DO bj = myByLo(myThid), myByHi(myThid)
226     DO bi = myBxLo(myThid), myBxHi(myThid)
227 jmc 1.47 DO k = 1,nLevOutp
228 jmc 1.3 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 jmc 1.54 ELSE
238 jmc 1.3 C- diagnostics is not empty :
239    
240 jmc 1.49 IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
241 jmc 1.54 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 jmc 1.15 & ' Computing Diagnostic # ', ndId, ' ', cdiag(ndId),
257     & ' Counter:',ndiag(ip,1,1),' Parms: ',gdiag(ndId)
258 jmc 1.54 ENDIF
259 jmc 1.29 IF ( mate.GT.0 ) THEN
260 jmc 1.35 WRITE(ioUnit,'(3A,I6,2A)')
261 jmc 1.15 & ' use Counter Mate for ', cdiag(ndId),
262     & ' Diagnostic # ',mate, ' ', cdiag(mate)
263 jmc 1.29 ELSEIF ( mVec.GT.0 ) THEN
264 jmc 1.15 IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN
265 jmc 1.35 WRITE(ioUnit,'(3A,I6,3A)')
266 jmc 1.15 & ' Vector Mate for ', cdiag(ndId),
267     & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
268     & ' exists '
269 jmc 1.3 ELSE
270 jmc 1.35 WRITE(ioUnit,'(3A,I6,3A)')
271 jmc 1.15 & ' Vector Mate for ', cdiag(ndId),
272     & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
273     & ' not enabled'
274 jmc 1.3 ENDIF
275     ENDIF
276 jmc 1.6 ENDIF
277 jmc 1.3
278 jmc 1.52 IF ( fflags(listId)(2:2).EQ.' ' ) THEN
279     C- get only selected levels:
280 jmc 1.30 DO bj = myByLo(myThid), myByHi(myThid)
281     DO bi = myBxLo(myThid), myBxHi(myThid)
282 jmc 1.52 DO k = 1,nlevels(listId)
283     kLev = NINT(levs(k,listId))
284     CALL DIAGNOSTICS_GET_DIAG(
285     I kLev, undefRL,
286 jmc 1.30 O qtmp1(1-OLx,1-OLy,k,bi,bj),
287 jmc 1.54 I ndId, mate, ip, im, bi, bj, myThid )
288 jmc 1.30 ENDDO
289     ENDDO
290     ENDDO
291 jmc 1.54 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 jmc 1.30 ELSE
305 jmc 1.52 C- get all the levels (for vertical post-processing)
306 jmc 1.30 DO bj = myByLo(myThid), myByHi(myThid)
307     DO bi = myBxLo(myThid), myBxHi(myThid)
308 jmc 1.52 CALL DIAGNOSTICS_GET_DIAG(
309     I 0, undefRL,
310     O qtmp1(1-OLx,1-OLy,1,bi,bj),
311 jmc 1.54 I ndId, mate, ip, im, bi, bj, myThid )
312 jmc 1.30 ENDDO
313 jmc 1.3 ENDDO
314 jmc 1.54 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 jmc 1.30 ENDIF
327 jmc 1.1
328 molod 1.17 C-----------------------------------------------------------------------
329 jmc 1.47 C-- Apply specific post-processing (e.g., interpolate) before output
330 molod 1.17 C-----------------------------------------------------------------------
331 jmc 1.47 IF ( fflags(listId)(2:2).EQ.'P' ) THEN
332     C- Do vertical interpolation:
333     IF ( fluidIsAir ) THEN
334 jmc 1.29 C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
335 jmc 1.47 CALL DIAGNOSTICS_INTERP_VERT(
336     I listId, md, ndId, ip, im, lm,
337 jmc 1.52 U qtmp1, qtmp2,
338     I undefRL, myTime, myIter, myThid )
339 jmc 1.47 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 jmc 1.52 I undefRL, myTime, myIter, myThid )
352 jmc 1.47 ENDIF
353 jmc 1.54 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 jmc 1.47
370     C-- End of empty diag / not-empty diag block
371 jmc 1.29 ENDIF
372 molod 1.17
373 jmc 1.47 C-- Ready to write field "md", element "lm" in averageCycle(listId)
374 jmc 1.31
375     C- write to binary file, using MDSIO pkg:
376 jmc 1.34 IF ( diag_mdsio ) THEN
377 jmc 1.55 c nRec = lm + (md-1)*averageCycle(listId)
378     nRec = md + (lm-1)*nfields(listId)
379 jmc 1.30 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 jmc 1.34 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 jmc 1.47 I NrMax, 1, nLevOutp,
388 jmc 1.34 I qtmp1, -nRec, myIter, myThid )
389 jmc 1.3 ENDIF
390 jmc 1.1
391     #ifdef ALLOW_MNC
392 jmc 1.3 IF (useMNC .AND. diag_mnc) THEN
393 jmc 1.50 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 jmc 1.3 ENDIF
400 jmc 1.1 #endif /* ALLOW_MNC */
401    
402 jmc 1.15 C-- end of Processing Fld # md
403 jmc 1.3 ENDIF
404 jmc 1.41 ENDDO
405    
406 jmc 1.55 C-- end loop on lm counter (= averagePeriod)
407 jmc 1.3 ENDDO
408 jmc 1.1
409 jmc 1.31 #ifdef ALLOW_MDSIO
410     IF (diag_mdsio) THEN
411 jmc 1.48 C- Note: temporary: since it is a pain to add more arguments to
412 jmc 1.31 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 jmc 1.34 glf = globalFiles
415 jmc 1.55 nRec = averageCycle(listId)*nfields(listId)
416 jmc 1.31 CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
417 jmc 1.47 & 0, 0, nLevOutp, ' ',
418 jmc 1.44 & nfields(listId), flds(1,listId), nTimRec, timeRec,
419 jmc 1.31 & nRec, myIter, myThid)
420     ENDIF
421     #endif /* ALLOW_MDSIO */
422    
423 jmc 1.15 RETURN
424 jmc 1.3 END
425 jmc 1.15
426 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

  ViewVC Help
Powered by ViewVC 1.1.22