/[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.57 - (hide annotations) (download)
Mon Jun 27 22:27:23 2011 UTC (12 years, 10 months ago) by jmc
Branch: MAIN
Changes since 1.56: +4 -7 lines
use run-time parameter "useMissingValue" to fill land-point (i.e., where mask=0)
 with MissingValue (MNC output file only).
This replaces CPP-option DIAGNOSTICS_MISSING_VALUE.

1 jmc 1.57 C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.56 2011/06/23 15:29:01 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     #endif /* ALLOW_MNC */
90    
91 jmc 1.3 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
92    
93 jmc 1.44 C--- set file properties
94 jmc 1.6 ioUnit= standardMessageUnit
95 jmc 1.52 undefRL = UNSET_RL
96 jmc 1.40 #ifdef ALLOW_FIZHI
97 jmc 1.52 IF ( useFIZHI ) undefRL = getcon('UNDEF')
98 jmc 1.40 #endif
99 jmc 1.56 IF ( misvalFlt(listId).NE.UNSET_RL ) undefRL = misvalFlt(listId)
100    
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.50 CALL DIAGNOSTICS_MNC_SET(
156 jmc 1.55 I nLevOutp, listId, lm,
157 jmc 1.57 O diag_mnc_bn,
158 jmc 1.56 I undefRL, myTime, myIter, myThid )
159 jmc 1.41 ENDIF
160 jmc 1.1 #endif /* ALLOW_MNC */
161    
162 jmc 1.29 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
163    
164 jmc 1.41 DO md = 1,nfields(listId)
165 jmc 1.15 ndId = jdiag(md,listId)
166 jmc 1.35 gcode = gdiag(ndId)(1:10)
167 jmc 1.29 mate = 0
168     mVec = 0
169 jmc 1.54 mDbl = 0
170 jmc 1.35 IF ( gcode(5:5).EQ.'C' ) THEN
171 jmc 1.29 C- Check for Mate of a Counter Diagnostic
172 jmc 1.35 mate = hdiag(ndId)
173 jmc 1.54 ELSEIF ( gcode(5:5).EQ.'P' ) THEN
174     C- Also load the mate (if stored) for Post-Processing
175     nn = ndId
176     DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
177     nn = hdiag(nn)
178     ENDDO
179     IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
180 jmc 1.35 ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
181 jmc 1.29 C- Check for Mate of a Vector Diagnostic
182 jmc 1.36 mVec = hdiag(ndId)
183 jmc 1.29 ENDIF
184 jmc 1.35 IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
185 jmc 1.3 C-- Start processing 1 Fld :
186    
187 jmc 1.29 ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
188 jmc 1.15 im = mdiag(md,listId)
189 jmc 1.29 IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
190 jmc 1.54 IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
191 jmc 1.29 IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
192    
193 jmc 1.54 IF ( ndiag(ip,1,1).EQ.0 ) THEN
194 jmc 1.3 C- Empty diagnostics case :
195    
196     _BEGIN_MASTER( myThid )
197     WRITE(msgBuf,'(A,I10)')
198     & '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
199 jmc 1.15 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
200 jmc 1.3 & SQUEEZE_RIGHT, myThid)
201 jmc 1.35 WRITE(msgBuf,'(A,I6,3A,I4,2A)')
202 jmc 1.15 & '- WARNING - diag.#',ndId, ' : ',flds(md,listId),
203     & ' (#',md,' ) in outp.Stream: ',fnames(listId)
204     CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
205 jmc 1.3 & SQUEEZE_RIGHT, myThid)
206 jmc 1.29 IF ( averageCycle(listId).GT.1 ) THEN
207 jmc 1.35 WRITE(msgBuf,'(A,2(I3,A))')
208 jmc 1.29 & '- WARNING - has not been filled (ndiag(lm=',lm,')=',
209     & ndiag(ip,1,1), ' )'
210     ELSE
211 jmc 1.35 WRITE(msgBuf,'(A,2(I3,A))')
212 jmc 1.29 & '- WARNING - has not been filled (ndiag=',
213     & ndiag(ip,1,1), ' )'
214     ENDIF
215 jmc 1.15 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
216 jmc 1.3 & SQUEEZE_RIGHT, myThid)
217     WRITE(msgBuf,'(A)')
218     & 'WARNING DIAGNOSTICS_OUT => write ZEROS instead'
219 jmc 1.15 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
220 jmc 1.3 & SQUEEZE_RIGHT, myThid)
221     _END_MASTER( myThid )
222     DO bj = myByLo(myThid), myByHi(myThid)
223     DO bi = myBxLo(myThid), myBxHi(myThid)
224 jmc 1.47 DO k = 1,nLevOutp
225 jmc 1.3 DO j = 1-OLy,sNy+OLy
226     DO i = 1-OLx,sNx+OLx
227     qtmp1(i,j,k,bi,bj) = 0. _d 0
228     ENDDO
229     ENDDO
230     ENDDO
231     ENDDO
232     ENDDO
233    
234 jmc 1.54 ELSE
235 jmc 1.3 C- diagnostics is not empty :
236    
237 jmc 1.49 IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
238 jmc 1.54 IF ( gcode(5:5).EQ.'P' ) THEN
239     WRITE(ioUnit,'(A,I6,7A,I8,2A)')
240     & ' Post-Processing Diag # ', ndId, ' ', cdiag(ndId),
241     & ' Parms: ',gdiag(ndId)
242     IF ( mDbl.EQ.0 ) THEN
243     WRITE(ioUnit,'(2(3A,I6,A,I8))') ' from diag: ',
244     & cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
245     ELSE
246     WRITE(ioUnit,'(2(3A,I6,A,I8))') ' from diag: ',
247     & cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
248     & ' and diag: ',
249     & cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
250     ENDIF
251     ELSE
252     WRITE(ioUnit,'(A,I6,3A,I8,2A)')
253 jmc 1.15 & ' Computing Diagnostic # ', ndId, ' ', cdiag(ndId),
254     & ' Counter:',ndiag(ip,1,1),' Parms: ',gdiag(ndId)
255 jmc 1.54 ENDIF
256 jmc 1.29 IF ( mate.GT.0 ) THEN
257 jmc 1.35 WRITE(ioUnit,'(3A,I6,2A)')
258 jmc 1.15 & ' use Counter Mate for ', cdiag(ndId),
259     & ' Diagnostic # ',mate, ' ', cdiag(mate)
260 jmc 1.29 ELSEIF ( mVec.GT.0 ) THEN
261 jmc 1.15 IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN
262 jmc 1.35 WRITE(ioUnit,'(3A,I6,3A)')
263 jmc 1.15 & ' Vector Mate for ', cdiag(ndId),
264     & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
265     & ' exists '
266 jmc 1.3 ELSE
267 jmc 1.35 WRITE(ioUnit,'(3A,I6,3A)')
268 jmc 1.15 & ' Vector Mate for ', cdiag(ndId),
269     & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
270     & ' not enabled'
271 jmc 1.3 ENDIF
272     ENDIF
273 jmc 1.6 ENDIF
274 jmc 1.3
275 jmc 1.52 IF ( fflags(listId)(2:2).EQ.' ' ) THEN
276     C- get only selected levels:
277 jmc 1.30 DO bj = myByLo(myThid), myByHi(myThid)
278     DO bi = myBxLo(myThid), myBxHi(myThid)
279 jmc 1.52 DO k = 1,nlevels(listId)
280     kLev = NINT(levs(k,listId))
281     CALL DIAGNOSTICS_GET_DIAG(
282     I kLev, undefRL,
283 jmc 1.30 O qtmp1(1-OLx,1-OLy,k,bi,bj),
284 jmc 1.54 I ndId, mate, ip, im, bi, bj, myThid )
285 jmc 1.30 ENDDO
286     ENDDO
287     ENDDO
288 jmc 1.54 IF ( mDbl.GT.0 ) THEN
289     DO bj = myByLo(myThid), myByHi(myThid)
290     DO bi = myBxLo(myThid), myBxHi(myThid)
291     DO k = 1,nlevels(listId)
292     kLev = NINT(levs(k,listId))
293     CALL DIAGNOSTICS_GET_DIAG(
294     I kLev, undefRL,
295     O qtmp2(1-OLx,1-OLy,k,bi,bj),
296     I mDbl, 0, im, 0, bi, bj, myThid )
297     ENDDO
298     ENDDO
299     ENDDO
300     ENDIF
301 jmc 1.30 ELSE
302 jmc 1.52 C- get all the levels (for vertical post-processing)
303 jmc 1.30 DO bj = myByLo(myThid), myByHi(myThid)
304     DO bi = myBxLo(myThid), myBxHi(myThid)
305 jmc 1.52 CALL DIAGNOSTICS_GET_DIAG(
306     I 0, undefRL,
307     O qtmp1(1-OLx,1-OLy,1,bi,bj),
308 jmc 1.54 I ndId, mate, ip, im, bi, bj, myThid )
309 jmc 1.30 ENDDO
310 jmc 1.3 ENDDO
311 jmc 1.54 IF ( mDbl.GT.0 ) THEN
312     DO bj = myByLo(myThid), myByHi(myThid)
313     DO bi = myBxLo(myThid), myBxHi(myThid)
314     DO k = 1,nlevels(listId)
315     CALL DIAGNOSTICS_GET_DIAG(
316     I 0, 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 jmc 1.30 ENDIF
324 jmc 1.1
325 molod 1.17 C-----------------------------------------------------------------------
326 jmc 1.47 C-- Apply specific post-processing (e.g., interpolate) before output
327 molod 1.17 C-----------------------------------------------------------------------
328 jmc 1.47 IF ( fflags(listId)(2:2).EQ.'P' ) THEN
329     C- Do vertical interpolation:
330     IF ( fluidIsAir ) THEN
331 jmc 1.29 C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
332 jmc 1.47 CALL DIAGNOSTICS_INTERP_VERT(
333     I listId, md, ndId, ip, im, lm,
334 jmc 1.52 U qtmp1, qtmp2,
335     I undefRL, myTime, myIter, myThid )
336 jmc 1.47 ELSE
337     WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
338     & 'INTERP_VERT not allowed in this config'
339     CALL PRINT_ERROR( msgBuf , myThid )
340     STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
341     ENDIF
342     ENDIF
343     IF ( fflags(listId)(2:2).EQ.'I' ) THEN
344     C- Integrate vertically: for now, output field has just 1 level:
345     CALL DIAGNOSTICS_SUM_LEVELS(
346     I listId, md, ndId, ip, im, lm,
347     U qtmp1,
348 jmc 1.52 I undefRL, myTime, myIter, myThid )
349 jmc 1.47 ENDIF
350 jmc 1.54 IF ( gcode(5:5).EQ.'P' ) THEN
351     C- Do Post-Processing:
352     IF ( flds(md,listId).EQ.'PhiVEL '
353     c & .OR. flds(md,listId).EQ.'PsiVEL '
354     & ) THEN
355     CALL DIAGNOSTICS_CALC_PHIVEL(
356     I listId, md, ndId, ip, im, lm,
357     U qtmp1, qtmp2,
358     I myTime, myIter, myThid )
359     ELSE
360     WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
361     & 'unknown Processing method for diag="',cdiag(ndId),'"'
362     CALL PRINT_ERROR( msgBuf , myThid )
363     STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
364     ENDIF
365     ENDIF
366 jmc 1.47
367     C-- End of empty diag / not-empty diag block
368 jmc 1.29 ENDIF
369 molod 1.17
370 jmc 1.47 C-- Ready to write field "md", element "lm" in averageCycle(listId)
371 jmc 1.31
372     C- write to binary file, using MDSIO pkg:
373 jmc 1.34 IF ( diag_mdsio ) THEN
374 jmc 1.55 c nRec = lm + (md-1)*averageCycle(listId)
375     nRec = md + (lm-1)*nfields(listId)
376 jmc 1.30 C default precision for output files
377     prec = writeBinaryPrec
378     C fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
379     IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
380     IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
381 jmc 1.34 C a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R
382     CALL WRITE_REC_LEV_RL(
383     I fn, prec,
384 jmc 1.47 I NrMax, 1, nLevOutp,
385 jmc 1.34 I qtmp1, -nRec, myIter, myThid )
386 jmc 1.3 ENDIF
387 jmc 1.1
388     #ifdef ALLOW_MNC
389 jmc 1.3 IF (useMNC .AND. diag_mnc) THEN
390 jmc 1.50 CALL DIAGNOSTICS_MNC_OUT(
391 jmc 1.56 I NrMax, nLevOutp, listId, ndId, mate,
392 jmc 1.57 I diag_mnc_bn, qtmp1,
393     I undefRL, myTime, myIter, myThid )
394 jmc 1.3 ENDIF
395 jmc 1.1 #endif /* ALLOW_MNC */
396    
397 jmc 1.15 C-- end of Processing Fld # md
398 jmc 1.3 ENDIF
399 jmc 1.41 ENDDO
400    
401 jmc 1.55 C-- end loop on lm counter (= averagePeriod)
402 jmc 1.3 ENDDO
403 jmc 1.1
404 jmc 1.31 #ifdef ALLOW_MDSIO
405     IF (diag_mdsio) THEN
406 jmc 1.48 C- Note: temporary: since it is a pain to add more arguments to
407 jmc 1.31 C all MDSIO S/R, uses instead this specific S/R to write only
408     C meta files but with more informations in it.
409 jmc 1.34 glf = globalFiles
410 jmc 1.55 nRec = averageCycle(listId)*nfields(listId)
411 jmc 1.31 CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
412 jmc 1.47 & 0, 0, nLevOutp, ' ',
413 jmc 1.44 & nfields(listId), flds(1,listId), nTimRec, timeRec,
414 jmc 1.31 & nRec, myIter, myThid)
415     ENDIF
416     #endif /* ALLOW_MDSIO */
417    
418 jmc 1.15 RETURN
419 jmc 1.3 END
420 jmc 1.15
421 jmc 1.1 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

  ViewVC Help
Powered by ViewVC 1.1.22