/[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.62 - (show annotations) (download)
Wed Jan 11 00:22:48 2017 UTC (7 years, 3 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint66e, checkpoint66d, checkpoint66c
Changes since 1.61: +7 -2 lines
- add run-time variable diagMdsDir to specify a directory where diagnostics will be written when using mds
- note: cannot be used with either pkg/mnc or mdsioLocalDir that alredy use subirectories

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

  ViewVC Help
Powered by ViewVC 1.1.22