/[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.64 - (show annotations) (download)
Sun Jul 23 00:24:18 2017 UTC (6 years, 9 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint66o, checkpoint66n, checkpoint66m, checkpoint66l, checkpoint66k, checkpoint66j, HEAD
Changes since 1.63: +2 -2 lines
allows for negative "jdiag" (interpret |jdiag| instead)

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

  ViewVC Help
Powered by ViewVC 1.1.22