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