/[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.29 - (show annotations) (download)
Mon Jun 5 18:17:23 2006 UTC (17 years, 11 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint58l_post, mitgcm_mapl_00, checkpoint58r_post, checkpoint58n_post, checkpoint58h_post, checkpoint58q_post, checkpoint58j_post, checkpoint58i_post, checkpoint58g_post, checkpoint58o_post, checkpoint58k_post, checkpoint58s_post, checkpoint58p_post, checkpoint58m_post
Changes since 1.28: +71 -44 lines
Implement periodic averaging diagnostics (e.g., mean seasonal cycle,
 mean diurnal cycle)

1 C $Header: /u/gcmpack/MITgcm/pkg/diagnostics/diagnostics_out.F,v 1.28 2006/02/07 15:52:02 edhill 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 myIter,
14 I myTime,
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 #ifdef ALLOW_FIZHI
30 #include "fizhi_SIZE.h"
31 #else
32 INTEGER Nrphys
33 PARAMETER (Nrphys=0)
34 #endif
35
36
37 C !INPUT PARAMETERS:
38 C listId :: Diagnostics list number being written
39 C myIter :: current iteration number
40 C myTime :: current time of simulation (s)
41 C myThid :: my Thread Id number
42 _RL myTime
43 INTEGER listId, myIter, myThid
44 CEOP
45
46 C !LOCAL VARIABLES:
47 C i,j,k :: loop 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 mate :: counter mate Id number (in available diagnostics list)
52 C ip :: diagnostics pointer to storage array
53 C im :: counter-mate pointer to storage array
54 INTEGER i, j, k, lm
55 INTEGER bi, bj
56 INTEGER md, ndId, ip, im
57 INTEGER mate, mVec
58 CHARACTER*8 parms1
59 _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)
60 _RL undef, getcon
61 EXTERNAL getcon
62 INTEGER ILNBLNK
63 EXTERNAL ILNBLNK
64 INTEGER ilen
65 INTEGER nlevsout
66
67 INTEGER ioUnit
68 CHARACTER*(MAX_LEN_FNAM) fn
69 CHARACTER*(MAX_LEN_MBUF) suff
70 CHARACTER*(MAX_LEN_MBUF) msgBuf
71 #ifdef ALLOW_MDSIO
72 LOGICAL glf
73 INTEGER nRec
74 #endif
75 #ifdef ALLOW_MNC
76 INTEGER ii
77 CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
78 INTEGER CW_DIMS, NLEN
79 PARAMETER ( CW_DIMS = 10 )
80 PARAMETER ( NLEN = 80 )
81 INTEGER dim(CW_DIMS), ib(CW_DIMS), ie(CW_DIMS)
82 CHARACTER*(NLEN) dn(CW_DIMS)
83 CHARACTER*(NLEN) d_cw_name
84 CHARACTER*(NLEN) dn_blnk
85 #ifdef DIAG_MNC_COORD_NEEDSWORK
86 CHARACTER*(5) ctmp
87 _RS ztmp(Nr+Nrphys)
88 #endif
89 #endif /* ALLOW_MNC */
90
91 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
92
93 ioUnit= standardMessageUnit
94 undef = getcon('UNDEF')
95 WRITE(suff,'(I10.10)') myIter
96 ilen = ILNBLNK(fnames(listId))
97 WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)
98
99 #ifdef ALLOW_MNC
100 IF (useMNC .AND. diag_mnc) THEN
101 DO i = 1,MAX_LEN_FNAM
102 diag_mnc_bn(i:i) = ' '
103 ENDDO
104 DO i = 1,NLEN
105 dn_blnk(i:i) = ' '
106 ENDDO
107 WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)
108
109 C Update the record dimension by writing the iteration number
110 CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)
111 CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',myTime,myThid)
112 CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)
113 CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',myIter,myThid)
114
115 C NOTE: at some point it would be a good idea to add a time_bounds
116 C variable that has dimension (2,T) and clearly denotes the
117 C beginning and ending times for each diagnostics period
118
119 dn(1)(1:NLEN) = dn_blnk(1:NLEN)
120 WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)
121 dim(1) = nlevels(listId)
122 ib(1) = 1
123 ie(1) = nlevels(listId)
124
125 CALL MNC_CW_ADD_GNAME('diag_levels', 1,
126 & dim, dn, ib, ie, myThid)
127 CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',
128 & 0,0, myThid)
129 CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',
130 & 'Idicies of vertical levels within the source arrays',
131 & myThid)
132
133 CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
134 & 'diag_levels', levs(1,listId), myThid)
135
136 CALL MNC_CW_DEL_VNAME('diag_levels', myThid)
137 CALL MNC_CW_DEL_GNAME('diag_levels', myThid)
138
139 #ifdef DIAG_MNC_COORD_NEEDSWORK
140 C This part has been placed in an #ifdef because, as its currently
141 C written, it will only work with variables defined on a dynamics
142 C grid. As we start using diagnostics for physics grids, ice
143 C levels, land levels, etc. the different vertical coordinate
144 C dimensions will have to be taken into account.
145
146 C 20051021 JMC & EH3 : We need to extend this so that a few
147 C variables each defined on different grids do not have the same
148 C vertical dimension names so we should be using a pattern such
149 C as: Z[uml]td000000 where the 't' is the type as specified by
150 C gdiag(10)
151
152 C Now define: Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx
153 ctmp(1:5) = 'mul '
154 DO i = 1,3
155 dn(1)(1:NLEN) = dn_blnk(1:NLEN)
156 WRITE(dn(1),'(3a,i6.6)') 'Z',ctmp(i:i),'d',nlevels(listId)
157 CALL MNC_CW_ADD_GNAME(dn(1), 1, dim, dn, ib, ie, myThid)
158 CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)
159
160 C The following three ztmp() loops should eventually be modified
161 C to reflect the fractional nature of levs(j,l) -- they should
162 C do something like:
163 C ztmp(j) = rC(INT(FLOOR(levs(j,l))))
164 C + ( rC(INT(FLOOR(levs(j,l))))
165 C + rC(INT(CEIL(levs(j,l)))) )
166 C / ( levs(j,l) - FLOOR(levs(j,l)) )
167 C for averaged levels.
168 IF (i .EQ. 1) THEN
169 DO j = 1,nlevels(listId)
170 ztmp(j) = rC(NINT(levs(j,listId)))
171 ENDDO
172 CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
173 & 'Dimensional coordinate value at the mid point',
174 & myThid)
175 ELSEIF (i .EQ. 2) THEN
176 DO j = 1,nlevels(listId)
177 ztmp(j) = rF(NINT(levs(j,listId)) + 1)
178 ENDDO
179 CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
180 & 'Dimensional coordinate value at the upper point',
181 & myThid)
182 ELSEIF (i .EQ. 3) THEN
183 DO j = 1,nlevels(listId)
184 ztmp(j) = rF(NINT(levs(j,listId)))
185 ENDDO
186 CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
187 & 'Dimensional coordinate value at the lower point',
188 & myThid)
189 ENDIF
190 CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)
191 CALL MNC_CW_DEL_VNAME(dn(1), myThid)
192 CALL MNC_CW_DEL_GNAME(dn(1), myThid)
193 ENDDO
194 #endif /* DIAG_MNC_COORD_NEEDSWORK */
195
196 ENDIF
197 #endif /* ALLOW_MNC */
198
199 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
200
201 DO md = 1,nfields(listId)
202 ndId = jdiag(md,listId)
203 parms1 = gdiag(ndId)(1:8)
204 mate = 0
205 mVec = 0
206 IF ( parms1(5:5).EQ.'C' ) THEN
207 C- Check for Mate of a Counter Diagnostic
208 READ(parms1,'(5X,I3)') mate
209 ELSEIF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN
210 C- Check for Mate of a Vector Diagnostic
211 READ(parms1,'(5X,I3)') mVec
212 ENDIF
213 IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN
214 C-- Start processing 1 Fld :
215 DO lm=1,averageCycle(listId)
216
217 ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
218 im = mdiag(md,listId)
219 IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
220 IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
221
222 IF ( ndiag(ip,1,1).EQ.0 ) THEN
223 C- Empty diagnostics case :
224
225 _BEGIN_MASTER( myThid )
226 WRITE(msgBuf,'(A,I10)')
227 & '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
228 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
229 & SQUEEZE_RIGHT, myThid)
230 WRITE(msgBuf,'(A,I4,3A,I3,2A)')
231 & '- WARNING - diag.#',ndId, ' : ',flds(md,listId),
232 & ' (#',md,' ) in outp.Stream: ',fnames(listId)
233 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
234 & SQUEEZE_RIGHT, myThid)
235 IF ( averageCycle(listId).GT.1 ) THEN
236 WRITE(msgBuf,'(A,2(I2,A))')
237 & '- WARNING - has not been filled (ndiag(lm=',lm,')=',
238 & ndiag(ip,1,1), ' )'
239 ELSE
240 WRITE(msgBuf,'(A,2(I2,A))')
241 & '- WARNING - has not been filled (ndiag=',
242 & ndiag(ip,1,1), ' )'
243 ENDIF
244 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
245 & SQUEEZE_RIGHT, myThid)
246 WRITE(msgBuf,'(A)')
247 & 'WARNING DIAGNOSTICS_OUT => write ZEROS instead'
248 CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
249 & SQUEEZE_RIGHT, myThid)
250 _END_MASTER( myThid )
251 DO bj = myByLo(myThid), myByHi(myThid)
252 DO bi = myBxLo(myThid), myBxHi(myThid)
253 DO k = 1,nlevels(listId)
254 DO j = 1-OLy,sNy+OLy
255 DO i = 1-OLx,sNx+OLx
256 qtmp1(i,j,k,bi,bj) = 0. _d 0
257 ENDDO
258 ENDDO
259 ENDDO
260 ENDDO
261 ENDDO
262
263 ELSE
264 C- diagnostics is not empty :
265
266 IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN
267 WRITE(ioUnit,'(A,I3,3A,I8,2A)')
268 & ' Computing Diagnostic # ', ndId, ' ', cdiag(ndId),
269 & ' Counter:',ndiag(ip,1,1),' Parms: ',gdiag(ndId)
270 IF ( mate.GT.0 ) THEN
271 WRITE(ioUnit,'(3A,I3,2A)')
272 & ' use Counter Mate for ', cdiag(ndId),
273 & ' Diagnostic # ',mate, ' ', cdiag(mate)
274 ELSEIF ( mVec.GT.0 ) THEN
275 IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN
276 WRITE(ioUnit,'(3A,I3,3A)')
277 & ' Vector Mate for ', cdiag(ndId),
278 & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
279 & ' exists '
280 ELSE
281 WRITE(ioUnit,'(3A,I3,3A)')
282 & ' Vector Mate for ', cdiag(ndId),
283 & ' Diagnostic # ',mVec, ' ', cdiag(mVec),
284 & ' not enabled'
285 ENDIF
286 ENDIF
287 ENDIF
288
289 DO bj = myByLo(myThid), myByHi(myThid)
290 DO bi = myBxLo(myThid), myBxHi(myThid)
291 DO k = 1,nlevels(listId)
292 CALL GETDIAG(
293 I levs(k,listId),undef,
294 O qtmp1(1-OLx,1-OLy,k,bi,bj),
295 I ndId,mate,ip,im,bi,bj,myThid)
296 ENDDO
297 ENDDO
298 ENDDO
299
300 C- end of empty diag / not empty block
301 ENDIF
302
303 nlevsout = nlevels(listId)
304
305 C-----------------------------------------------------------------------
306 C Check to see if we need to interpolate before output
307 C-----------------------------------------------------------------------
308 IF ( fflags(listId)(2:2).EQ.'P' ) THEN
309 C- Do vertical interpolation:
310 c IF ( fluidIsAir ) THEN
311 C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
312 C find some problems with 5-levels AIM => use it only with FIZHI
313 IF ( useFIZHI ) THEN
314 CALL DIAGNOSTICS_INTERP_VERT(
315 I listId, md, ndId, ip, im,
316 U nlevsout,
317 U qtmp1,
318 I undef,
319 I myTime, myIter, myThid )
320 ELSE
321 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
322 & 'INTERP_VERT not safe in this config'
323 CALL PRINT_ERROR( msgBuf , myThid )
324 WRITE(msgBuf,'(2A,I3,2A)') 'DIAGNOSTICS_OUT: ',
325 & ' for list l=', listId, ', filename: ', fnames(listId)
326 CALL PRINT_ERROR( msgBuf , myThid )
327 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
328 ENDIF
329 ENDIF
330
331 #ifdef ALLOW_MDSIO
332 C Prepare for mdsio optionality
333 IF (diag_mdsio) THEN
334 glf = globalFiles
335 nRec = lm + (md-1)*averageCycle(listId)
336 IF (fflags(listId)(1:1) .EQ. 'R') THEN
337 C Force it to be 32-bit precision
338 CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,.FALSE.,
339 & 'RL',Nr+Nrphys,nlevsout,qtmp1,nRec,myIter,myThid)
340 ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN
341 C Force it to be 64-bit precision
342 CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,.FALSE.,
343 & 'RL',Nr+Nrphys,nlevsout,qtmp1,nRec,myIter,myThid)
344 ELSE
345 C This is the old default behavior
346 CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,.FALSE.,
347 & 'RL',Nr+Nrphys,nlevsout,qtmp1,nRec,myIter,myThid)
348 ENDIF
349 ENDIF
350 #endif /* ALLOW_MDSIO */
351
352 #ifdef ALLOW_MNC
353 IF (useMNC .AND. diag_mnc) THEN
354
355 _BEGIN_MASTER( myThid )
356
357 DO ii = 1,CW_DIMS
358 d_cw_name(1:NLEN) = dn_blnk(1:NLEN)
359 dn(ii)(1:NLEN) = dn_blnk(1:NLEN)
360 ENDDO
361
362 C Note that the "d_cw_name" variable is a hack that hides a
363 C subtlety within MNC. Basically, each MNC-wrapped file is
364 C caching its own concept of what each "grid name" (that is, a
365 C dimension group name) means. So one cannot re-use the same
366 C "grid" name for different collections of dimensions within a
367 C given file. By appending the "ndId" values to each name, we
368 C guarantee uniqueness within each MNC-produced file.
369 WRITE(d_cw_name,'(a,i6.6)') 'd_cw_',ndId
370
371 C XY dimensions
372 dim(1) = sNx + 2*OLx
373 dim(2) = sNy + 2*OLy
374 ib(1) = OLx + 1
375 ib(2) = OLy + 1
376 IF (gdiag(ndId)(2:2) .EQ. 'M') THEN
377 dn(1)(1:2) = 'X'
378 ie(1) = OLx + sNx
379 dn(2)(1:2) = 'Y'
380 ie(2) = OLy + sNy
381 ELSEIF (gdiag(ndId)(2:2) .EQ. 'U') THEN
382 dn(1)(1:3) = 'Xp1'
383 ie(1) = OLx + sNx + 1
384 dn(2)(1:2) = 'Y'
385 ie(2) = OLy + sNy
386 ELSEIF (gdiag(ndId)(2:2) .EQ. 'V') THEN
387 dn(1)(1:2) = 'X'
388 ie(1) = OLx + sNx
389 dn(2)(1:3) = 'Yp1'
390 ie(2) = OLy + sNy + 1
391 ELSEIF (gdiag(ndId)(2:2) .EQ. 'Z') THEN
392 dn(1)(1:3) = 'Xp1'
393 ie(1) = OLx + sNx + 1
394 dn(2)(1:3) = 'Yp1'
395 ie(2) = OLy + sNy + 1
396 ENDIF
397
398 C Z is special since it varies
399 WRITE(dn(3),'(a,i6.6)') 'Zd', nlevsout
400 IF ( (gdiag(ndId)(10:10) .EQ. 'R')
401 & .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN
402 WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevsout
403 ENDIF
404 IF ( (gdiag(ndId)(10:10) .EQ. 'R')
405 & .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN
406 WRITE(dn(3),'(a,i6.6)') 'Zld', nlevsout
407 ENDIF
408 IF ( (gdiag(ndId)(10:10) .EQ. 'R')
409 & .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
410 WRITE(dn(3),'(a,i6.6)') 'Zud', nlevsout
411 ENDIF
412 dim(3) = Nr+Nrphys
413 ib(3) = 1
414 ie(3) = nlevsout
415
416 C Time dimension
417 dn(4)(1:1) = 'T'
418 dim(4) = -1
419 ib(4) = 1
420 ie(4) = 1
421
422 CALL MNC_CW_ADD_GNAME(d_cw_name, 4,
423 & dim, dn, ib, ie, myThid)
424 CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,
425 & 4,5, myThid)
426 CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',
427 & tdiag(ndId),myThid)
428 CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',
429 & udiag(ndId),myThid)
430
431 C Per the observations of Baylor, this has been commented out
432 C until we have code that can write missing_value attributes
433 C in a way thats compatible with most of the more popular
434 C netCDF tools including ferret. Using all-zeros completely
435 C breaks ferret.
436
437 C CALL MNC_CW_ADD_VATTR_DBL(cdiag(ndId),'missing_value',
438 C & 0.0 _d 0,myThid)
439
440 IF ( ( (writeBinaryPrec .EQ. precFloat32)
441 & .AND. (fflags(listId)(1:1) .NE. 'D')
442 & .AND. (fflags(listId)(1:1) .NE. 'R') )
443 & .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN
444 CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,
445 & cdiag(ndId), qtmp1, myThid)
446 ELSEIF ( (writeBinaryPrec .EQ. precFloat64)
447 & .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN
448 CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
449 & cdiag(ndId), qtmp1, myThid)
450 ENDIF
451
452 CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)
453 CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)
454
455 _END_MASTER( myThid )
456
457 ENDIF
458 #endif /* ALLOW_MNC */
459
460 ENDDO
461 C-- end of Processing Fld # md
462 ENDIF
463 ENDDO
464
465 RETURN
466 END
467
468 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

  ViewVC Help
Powered by ViewVC 1.1.22