/[MITgcm]/MITgcm/pkg/diagnostics/diagnostics_out.F
ViewVC logotype

Diff of /MITgcm/pkg/diagnostics/diagnostics_out.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.16 by edhill, Thu Jul 7 15:32:35 2005 UTC revision 1.49 by jmc, Mon Jun 6 15:42:58 2011 UTC
# Line 8  CBOP 0 Line 8  CBOP 0
8  C     !ROUTINE: DIAGNOSTICS_OUT  C     !ROUTINE: DIAGNOSTICS_OUT
9    
10  C     !INTERFACE:  C     !INTERFACE:
11        SUBROUTINE  DIAGNOSTICS_OUT(        SUBROUTINE DIAGNOSTICS_OUT(
12       I     listId,       I     listId,
13       I     myIter,       I     myIter,
14       I     myTime,       I     myTime,
# Line 26  C     !USES: Line 26  C     !USES:
26  #include "DIAGNOSTICS_SIZE.h"  #include "DIAGNOSTICS_SIZE.h"
27  #include "DIAGNOSTICS.h"  #include "DIAGNOSTICS.h"
28    
29  #ifdef ALLOW_FIZHI        INTEGER NrMax
30  #include "fizhi_SIZE.h"        PARAMETER( NrMax = numLevels )
 #else  
       INTEGER Nrphys  
       PARAMETER (Nrphys=0)  
 #endif  
   
31    
32  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
33  C     listId  :: Diagnostics list number being written  C     listId  :: Diagnostics list number being written
# Line 43  C     myThid  :: my Thread Id number Line 38  C     myThid  :: my Thread Id number
38        INTEGER listId, myIter, myThid        INTEGER listId, myIter, myThid
39  CEOP  CEOP
40    
41    C     !FUNCTIONS:
42          INTEGER ILNBLNK
43          EXTERNAL ILNBLNK
44    #ifdef ALLOW_FIZHI
45          _RL   getcon
46          EXTERNAL getcon
47    #endif
48    
49  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
50  C     i,j,k :: loop indices  C     i,j,k :: loop indices
51    C     bi,bj :: tile indices
52    C     lm    :: loop index (averageCycle)
53  C     md    :: field number in the list "listId".  C     md    :: field number in the list "listId".
54  C     ndId  :: diagnostics  Id number (in available diagnostics list)  C     ndId  :: diagnostics  Id number (in available diagnostics list)
55  C     mate  :: counter mate Id number (in available diagnostics list)  C     mate  :: counter mate Id number (in available diagnostics list)
56  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
57  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
58        INTEGER i, j, k  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 :: thread-shared temporary array (needs to be in common block):
62    C              to write a diagnostic field to file, copy it first from (big)
63    C              diagnostic storage qdiag into it.
64          COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
65          _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
66    
67          INTEGER i, j, k, lm
68        INTEGER bi, bj        INTEGER bi, bj
69        INTEGER md, ndId, ip, im        INTEGER md, ndId, ip, im
70        INTEGER mate, mVec        INTEGER mate, mVec
71        CHARACTER*8 parms1        CHARACTER*10 gcode
72        CHARACTER*3 mate_index        _RL undef
73        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)        _RL tmpLev
       _RL undef, getcon  
       EXTERNAL getcon  
       INTEGER ILNBLNK  
       EXTERNAL ILNBLNK  
74        INTEGER ilen        INTEGER ilen
75          INTEGER nLevOutp
76    
77        INTEGER ioUnit        INTEGER ioUnit
78        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
79        CHARACTER*(MAX_LEN_MBUF) suff        CHARACTER*(MAX_LEN_MBUF) suff
80        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
81          INTEGER prec, nRec, nTimRec
82          _RL     timeRec(2)
83    #ifdef ALLOW_MDSIO
84        LOGICAL glf        LOGICAL glf
85    #endif
86  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
87        INTEGER ii        INTEGER ll, llMx, jj, jjMx
88          INTEGER ii, klev
89        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
       CHARACTER*(5) ctmp  
90        INTEGER CW_DIMS, NLEN        INTEGER CW_DIMS, NLEN
91        PARAMETER ( CW_DIMS = 10 )        PARAMETER ( CW_DIMS = 10 )
92        PARAMETER ( NLEN    = 80 )        PARAMETER ( NLEN    = 80 )
# Line 79  C     im    :: counter-mate pointer to s Line 94  C     im    :: counter-mate pointer to s
94        CHARACTER*(NLEN) dn(CW_DIMS)        CHARACTER*(NLEN) dn(CW_DIMS)
95        CHARACTER*(NLEN) d_cw_name        CHARACTER*(NLEN) d_cw_name
96        CHARACTER*(NLEN) dn_blnk        CHARACTER*(NLEN) dn_blnk
97        _RS ztmp(Nr+Nrphys)  #ifdef DIAG_MNC_COORD_NEEDSWORK
98          CHARACTER*(5) ctmp
99          _RS ztmp(NrMax)
100    #endif
101          LOGICAL useMissingValue, useMisValForThisDiag
102          REAL*8 misvalLoc
103          REAL*8 misval_r8(2)
104          REAL*4 misval_r4(2)
105          INTEGER misvalIntLoc, misval_int(2)
106  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
107    
108  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
109    
110    C---  set file properties
111        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
112          undef = UNSET_RL
113    #ifdef ALLOW_FIZHI
114    c     IF ( useFIZHI ) undef = getcon('UNDEF')
115        undef = getcon('UNDEF')        undef = getcon('UNDEF')
116        glf = globalFiles  #endif
117        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
118        ilen = ILNBLNK(fnames(listId))        ilen = ILNBLNK(fnames(listId))
119        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)
120    C-    for now, if integrate vertically, output field has just 1 level:
121          nLevOutp = nlevels(listId)
122          IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
123    
124    C--   Set time information:
125          IF ( freq(listId).LT.0. ) THEN
126    C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
127            nTimRec = 1
128            timeRec(1) = myTime
129          ELSE
130    C-    Time-average: store the 2 edges of the time-averaging interval.
131    C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
132    C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
133            nTimRec = 2
134    
135    C-    end of time-averaging interval:
136            timeRec(2) = myTime
137    
138    C-    begining of time-averaging interval:
139    c       timeRec(1) = myTime - freq(listId)
140    C     a) find the time of the previous multiple of output freq:
141            timeRec(1) = myTime-deltaTClock*0.5 _d 0
142            timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
143            i = INT( timeRec(1) )
144            IF ( timeRec(1).LT.0. ) THEN
145              tmpLev = FLOAT(i)
146              IF ( timeRec(1).NE.tmpLev ) i = i - 1
147            ENDIF
148            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
149    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
150            timeRec(1) = MAX( timeRec(1), startTime )
151    
152    C     b) round off to nearest multiple of time-step:
153            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
154            i = NINT( timeRec(1) )
155    C     if just half way, NINT will return the next time-step: correct this
156            tmpLev = FLOAT(i) - 0.5 _d 0
157            IF ( timeRec(1).EQ.tmpLev ) i = i - 1
158            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
159    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
160          ENDIF
161    C--   Convert time to iteration number (debug)
162    c     DO i=1,nTimRec
163    c       timeRec(i) = timeRec(i)/deltaTClock
164    c     ENDDO
165    
166  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
167    C-- this is a trick to reverse the order of the loops on md (= field)
168    C   and lm (= averagePeriod): binary output: lm loop inside md loop ;
169    C                                 mnc ouput: md loop inside lm loop.
170        IF (useMNC .AND. diag_mnc) THEN        IF (useMNC .AND. diag_mnc) THEN
171          DO i = 1,MAX_LEN_FNAM          jjMx = averageCycle(listId)
172            diag_mnc_bn(i:i) = ' '          llMx = 1
173          ENDDO        ELSE
174          DO i = 1,NLEN          jjMx = 1
175            dn_blnk(i:i) = ' '          llMx = averageCycle(listId)
176          ENDDO        ENDIF
177          WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)        DO jj=1,jjMx
178    
179           IF (useMNC .AND. diag_mnc) THEN
180    C     Handle missing value attribute (land points)
181             useMissingValue = .FALSE.
182    #ifdef DIAGNOSTICS_MISSING_VALUE
183             useMissingValue = .TRUE.
184    #endif /* DIAGNOSTICS_MISSING_VALUE */
185             IF ( misvalFlt(listId) .NE. UNSET_RL ) THEN
186              misvalLoc = misvalFlt(listId)
187             ELSE
188              misvalLoc = undef
189             ENDIF
190    C     Defaults to UNSET_I
191             misvalIntLoc = misvalInt(listId)
192             DO ii=1,2
193    C         misval_r4(ii)  = UNSET_FLOAT4
194    C         misval_r8(ii)  = UNSET_FLOAT8
195              misval_r4(ii)  = misvalLoc
196              misval_r8(ii)  = misvalLoc
197              misval_int(ii) = UNSET_I
198             ENDDO
199             DO i = 1,MAX_LEN_FNAM
200               diag_mnc_bn(i:i) = ' '
201             ENDDO
202             DO i = 1,NLEN
203               dn_blnk(i:i) = ' '
204             ENDDO
205             WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)
206    
207  C       Update the record dimension by writing the iteration number  C       Update the record dimension by writing the iteration number
208          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)           klev = myIter + jj - jjMx
209          CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',myTime,myThid)           tmpLev = myTime + deltaTClock*(jj -jjMx)
210          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)           CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)
211             CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',tmpLev,myThid)
212          dn(1)(1:NLEN) = dn_blnk(1:NLEN)           CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)
213          WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)           CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',klev,myThid)
214          dim(1) = nlevels(listId)  
215          ib(1)  = 1  C       NOTE: at some point it would be a good idea to add a time_bounds
216          ie(1)  = nlevels(listId)  C       variable that has dimension (2,T) and clearly denotes the
217    C       beginning and ending times for each diagnostics period
218          CALL MNC_CW_ADD_GNAME('diag_levels', 1,  
219       &       dim, dn, ib, ie, myThid)           dn(1)(1:NLEN) = dn_blnk(1:NLEN)
220          CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',           WRITE(dn(1),'(a,i6.6)') 'Zmd', nLevOutp
221       &       0,0, myThid)           dim(1) = nLevOutp
222          CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',           ib(1)  = 1
223       &       'Idicies of vertical levels within the source arrays',           ie(1)  = nLevOutp
224       &       myThid)  
225                     CALL MNC_CW_ADD_GNAME('diag_levels', 1,
226          CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,       &        dim, dn, ib, ie, myThid)
227       &       'diag_levels', levs(1,listId), myThid)           CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',
228         &        0,0, myThid)
229             CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',
230         &        'Idicies of vertical levels within the source arrays',
231         &        myThid)
232    C     suppress the missing value attribute (iflag = 0)
233             IF (useMissingValue)
234         &       CALL MNC_CW_VATTR_MISSING('diag_levels', 0,
235         I       misval_r8, misval_r4, misval_int,
236         I       myThid )
237    
238          CALL MNC_CW_DEL_VNAME('diag_levels', myThid)           CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
239          CALL MNC_CW_DEL_GNAME('diag_levels', myThid)       &        'diag_levels', levs(1,listId), myThid)
240    
241             CALL MNC_CW_DEL_VNAME('diag_levels', myThid)
242             CALL MNC_CW_DEL_GNAME('diag_levels', myThid)
243    
244  #ifdef DIAG_MNC_COORD_NEEDSWORK  #ifdef DIAG_MNC_COORD_NEEDSWORK
245  C       This part has been placed in an #ifdef because, as its currently  C       This part has been placed in an #ifdef because, as its currently
# Line 133  C       grid.  As we start using diagnos Line 248  C       grid.  As we start using diagnos
248  C       levels, land levels, etc. the different vertical coordinate  C       levels, land levels, etc. the different vertical coordinate
249  C       dimensions will have to be taken into account.  C       dimensions will have to be taken into account.
250    
251    C       20051021 JMC & EH3 : We need to extend this so that a few
252    C       variables each defined on different grids do not have the same
253    C       vertical dimension names so we should be using a pattern such
254    C       as: Z[uml]td000000 where the 't' is the type as specified by
255    C       gdiag(10)
256    
257  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx
258          ctmp(1:5) = 'mul  '           ctmp(1:5) = 'mul  '
259          DO i = 1,3           DO i = 1,3
260            dn(1)(1:NLEN) = dn_blnk(1:NLEN)             dn(1)(1:NLEN) = dn_blnk(1:NLEN)
261            WRITE(dn(1),'(3a,i6.6)') 'Z',ctmp(i:i),'d',nlevels(listId)             WRITE(dn(1),'(3a,i6.6)') 'Z',ctmp(i:i),'d',nlevels(listId)
262            CALL MNC_CW_ADD_GNAME(dn(1), 1, dim, dn, ib, ie, myThid)             CALL MNC_CW_ADD_GNAME(dn(1), 1, dim, dn, ib, ie, myThid)
263            CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)             CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)
264    
265  C         The following three ztmp() loops should eventually be modified  C         The following three ztmp() loops should eventually be modified
266  C         to reflect the fractional nature of levs(j,l) -- they should  C         to reflect the fractional nature of levs(j,l) -- they should
267  C         do something like:  C         do something like:
268  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))
269  C                      + ( rC(INT(FLOOR(levs(j,l))))  C                      + ( rC(INT(FLOOR(levs(j,l))))
270  C                          + rC(INT(CEIL(levs(j,l)))) )  C                          + rC(INT(CEIL(levs(j,l)))) )
271  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )
272  C         for averaged levels.  C         for averaged levels.
273            IF (i .EQ. 1) THEN             IF (i .EQ. 1) THEN
274              DO j = 1,nlevels(listId)               DO j = 1,nlevels(listId)
275                ztmp(j) = rC(NINT(levs(j,listId)))                 ztmp(j) = rC(NINT(levs(j,listId)))
276              ENDDO               ENDDO
277              CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',               CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
278       &           'Dimensional coordinate value at the mid point',       &            'Dimensional coordinate value at the mid point',
279       &           myThid)       &            myThid)
280            ELSEIF (i .EQ. 2) THEN             ELSEIF (i .EQ. 2) THEN
281              DO j = 1,nlevels(listId)               DO j = 1,nlevels(listId)
282                ztmp(j) = rF(NINT(levs(j,listId)) + 1)                 ztmp(j) = rF(NINT(levs(j,listId)) + 1)
283              ENDDO               ENDDO
284              CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',               CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
285       &           'Dimensional coordinate value at the upper point',       &            'Dimensional coordinate value at the upper point',
286       &           myThid)       &            myThid)
287            ELSEIF (i .EQ. 3) THEN             ELSEIF (i .EQ. 3) THEN
288              DO j = 1,nlevels(listId)               DO j = 1,nlevels(listId)
289                ztmp(j) = rF(NINT(levs(j,listId)))                 ztmp(j) = rF(NINT(levs(j,listId)))
290              ENDDO               ENDDO
291              CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',               CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
292       &           'Dimensional coordinate value at the lower point',       &            'Dimensional coordinate value at the lower point',
293       &           myThid)       &            myThid)
294            ENDIF             ENDIF
295            CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)  C     suppress the missing value attribute (iflag = 0)
296            CALL MNC_CW_DEL_VNAME(dn(1), myThid)             IF (useMissingValue)
297            CALL MNC_CW_DEL_GNAME(dn(1), myThid)       &          CALL MNC_CW_VATTR_MISSING(dn(1), 0,
298          ENDDO       I          misval_r8, misval_r4, misval_int,
299         I          myThid )
300               CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)
301               CALL MNC_CW_DEL_VNAME(dn(1), myThid)
302               CALL MNC_CW_DEL_GNAME(dn(1), myThid)
303             ENDDO
304  #endif /*  DIAG_MNC_COORD_NEEDSWORK  */  #endif /*  DIAG_MNC_COORD_NEEDSWORK  */
305    
306        ENDIF         ENDIF
307  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
308    
309        DO md = 1,nfields(listId)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
310    
311           DO md = 1,nfields(listId)
312          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
313          parms1 = gdiag(ndId)(1:8)          gcode = gdiag(ndId)(1:10)
314          IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN          mate = 0
315            mVec = 0
316            IF ( gcode(5:5).EQ.'C' ) THEN
317    C-      Check for Mate of a Counter Diagnostic
318               mate = hdiag(ndId)
319            ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
320    C-      Check for Mate of a Vector Diagnostic
321               mVec = hdiag(ndId)
322            ENDIF
323            IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
324  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
325    #ifdef ALLOW_MNC
326             DO ll=1,llMx
327              lm = jj+ll-1
328    #else
329             DO lm=1,averageCycle(listId)
330    #endif
331    
332            ip = ABS(idiag(md,listId))            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
333            im = mdiag(md,listId)            im = mdiag(md,listId)
334              IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
335              IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
336    
337            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ndiag(ip,1,1).EQ.0 ) THEN
338  C-        Empty diagnostics case :  C-        Empty diagnostics case :
339    
# Line 196  C-        Empty diagnostics case : Line 342  C-        Empty diagnostics case :
342       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
343              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
344       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
345              WRITE(msgBuf,'(A,I4,3A,I3,2A)')              WRITE(msgBuf,'(A,I6,3A,I4,2A)')
346       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
347       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
348              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
349       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
350              WRITE(msgBuf,'(A,I2,A)')              IF ( averageCycle(listId).GT.1 ) THEN
351       &       '- WARNING -   has not been filled (ndiag=',               WRITE(msgBuf,'(A,2(I3,A))')
352       &       ndiag(ip,1,1), ' )'       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
353         &                                            ndiag(ip,1,1), ' )'
354                ELSE
355                 WRITE(msgBuf,'(A,2(I3,A))')
356         &        '- WARNING -   has not been filled (ndiag=',
357         &                                            ndiag(ip,1,1), ' )'
358                ENDIF
359              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
360       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
361              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
# Line 213  C-        Empty diagnostics case : Line 365  C-        Empty diagnostics case :
365              _END_MASTER( myThid )              _END_MASTER( myThid )
366              DO bj = myByLo(myThid), myByHi(myThid)              DO bj = myByLo(myThid), myByHi(myThid)
367                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
368                  DO k = 1,nlevels(listId)                  DO k = 1,nLevOutp
369                    DO j = 1-OLy,sNy+OLy                    DO j = 1-OLy,sNy+OLy
370                      DO i = 1-OLx,sNx+OLx                      DO i = 1-OLx,sNx+OLx
371                        qtmp1(i,j,k,bi,bj) = 0. _d 0                        qtmp1(i,j,k,bi,bj) = 0. _d 0
# Line 226  C-        Empty diagnostics case : Line 378  C-        Empty diagnostics case :
378            ELSE            ELSE
379  C-        diagnostics is not empty :  C-        diagnostics is not empty :
380    
381              IF ( myThid.EQ.1 ) WRITE(ioUnit,'(A,I3,3A,I8,2A)')              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
382                  WRITE(ioUnit,'(A,I6,3A,I8,2A)')
383       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
384       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
385                  IF ( mate.GT.0 ) THEN
386              IF ( parms1(5:5).EQ.'C' ) THEN                 WRITE(ioUnit,'(3A,I6,2A)')
 C             Check for Mate of a Counter Diagnostic  
 C             --------------------------------------  
               mate_index = parms1(6:8)  
               READ (mate_index,'(I3)') mate  
               IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,2A)')  
387       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
388       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
389                  ELSEIF ( mVec.GT.0 ) THEN
             ELSE  
               mate = 0  
   
 C             Check for Mate of a Vector Diagnostic  
 C             -------------------------------------  
               IF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN  
                 mate_index = parms1(6:8)  
                 READ (mate_index,'(I3)') mVec  
390                  IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN                  IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN
391                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
392       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
393       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
394       &             ' exists '       &             ' exists '
395                  ELSE                  ELSE
396                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
397       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
398       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
399       &             ' not enabled'       &             ' not enabled'
# Line 261  C             -------------------------- Line 401  C             --------------------------
401                ENDIF                ENDIF
402              ENDIF              ENDIF
403    
404              DO bj = myByLo(myThid), myByHi(myThid)              IF ( fflags(listId)(2:2).NE.' ' ) THEN
405               DO bi = myBxLo(myThid), myBxHi(myThid)  C-       get all the levels (for vertical post-processing)
406                DO k = 1,nlevels(listId)                DO bj = myByLo(myThid), myByHi(myThid)
407                  CALL GETDIAG(                 DO bi = myBxLo(myThid), myBxHi(myThid)
408       I                       levs(k,listId),undef,                  DO k = 1,kdiag(ndId)
409       O                       qtmp1(1-OLx,1-OLy,k,bi,bj),                    tmpLev = k
410       I                       ndId,mate,ip,im,bi,bj,myThid)                    CALL GETDIAG(
411         I                         tmpLev,undef,
412         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
413         I                         ndId,mate,ip,im,bi,bj,myThid)
414                    ENDDO
415                   ENDDO
416                ENDDO                ENDDO
417               ENDDO              ELSE
418              ENDDO  C-       get only selected levels:
419                  DO bj = myByLo(myThid), myByHi(myThid)
420                   DO bi = myBxLo(myThid), myBxHi(myThid)
421                    DO k = 1,nlevels(listId)
422                      CALL GETDIAG(
423         I                         levs(k,listId),undef,
424         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
425         I                         ndId,mate,ip,im,bi,bj,myThid)
426                    ENDDO
427                   ENDDO
428                  ENDDO
429                ENDIF
430    
431    C-----------------------------------------------------------------------
432    C--     Apply specific post-processing (e.g., interpolate) before output
433    C-----------------------------------------------------------------------
434                IF ( fflags(listId)(2:2).EQ.'P' ) THEN
435    C-          Do vertical interpolation:
436                 IF ( fluidIsAir ) THEN
437    C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
438                  CALL DIAGNOSTICS_INTERP_VERT(
439         I                         listId, md, ndId, ip, im, lm,
440         U                         qtmp1,
441         I                         undef, myTime, myIter, myThid )
442                 ELSE
443                   WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
444         &           'INTERP_VERT not allowed in this config'
445                   CALL PRINT_ERROR( msgBuf , myThid )
446                   STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
447                 ENDIF
448                ENDIF
449                IF ( fflags(listId)(2:2).EQ.'I' ) THEN
450    C-          Integrate vertically: for now, output field has just 1 level:
451                  CALL DIAGNOSTICS_SUM_LEVELS(
452         I                         listId, md, ndId, ip, im, lm,
453         U                         qtmp1,
454         I                         undef, myTime, myIter, myThid )
455                ENDIF
456    
457  C-        end of empty diag / not empty block  C--     End of empty diag / not-empty diag block
458            ENDIF            ENDIF
459    
460  #ifdef ALLOW_MDSIO  C--     Ready to write field "md", element "lm" in averageCycle(listId)
461  C         Prepare for mdsio optionality  
462            IF (diag_mdsio) THEN  C-        write to binary file, using MDSIO pkg:
463              IF (fflags(listId)(1:1) .EQ. ' ') THEN            IF ( diag_mdsio ) THEN
464  C             This is the old default behavior              nRec = lm + (md-1)*averageCycle(listId)
465                CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,'RL',  C           default precision for output files
466       &             Nr+Nrphys,nlevels(listId),qtmp1,md,myIter,myThid)              prec = writeBinaryPrec
467              ELSEIF (fflags(listId)(1:1) .EQ. 'R') THEN  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
468  C             Force it to be 32-bit precision              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
469                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,'RL',              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
470       &             Nr+Nrphys,nlevels(listId),qtmp1,md,myIter,myThid)  C         a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R
471              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              CALL WRITE_REC_LEV_RL(
472  C             Force it to be 64-bit precision       I                            fn, prec,
473                CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,'RL',       I                            NrMax, 1, nLevOutp,
474       &             Nr+Nrphys,nlevels(listId),qtmp1,md,myIter,myThid)       I                            qtmp1, -nRec, myIter, myThid )
             ENDIF  
475            ENDIF            ENDIF
 #endif /*  ALLOW_MDSIO  */  
476    
477  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
478            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
# Line 318  C           XY dimensions Line 498  C           XY dimensions
498              dim(2)       = sNy + 2*OLy              dim(2)       = sNy + 2*OLy
499              ib(1)        = OLx + 1              ib(1)        = OLx + 1
500              ib(2)        = OLy + 1              ib(2)        = OLy + 1
501              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN
502                dn(1)(1:2) = 'X'                dn(1)(1:2) = 'X'
503                ie(1)      = OLx + sNx                ie(1)      = OLx + sNx
504                dn(2)(1:2) = 'Y'                dn(2)(1:2) = 'Y'
# Line 339  C           XY dimensions Line 519  C           XY dimensions
519                dn(2)(1:3) = 'Yp1'                dn(2)(1:3) = 'Yp1'
520                ie(2)      = OLy + sNy + 1                ie(2)      = OLy + sNy + 1
521              ENDIF              ENDIF
522                
523  C           Z is special since it varies  C           Z is special since it varies
524              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)              WRITE(dn(3),'(a,i6.6)') 'Zd', nLevOutp
525              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
526       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN
527                WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zmd', nLevOutp
528              ENDIF              ENDIF
529              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
530       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN
531                WRITE(dn(3),'(a,i6.6)') 'Zld', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zld', nLevOutp
532              ENDIF              ENDIF
533              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
534       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
535                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zud', nLevOutp
536              ENDIF              ENDIF
537              dim(3) = Nr+Nrphys              dim(3) = NrMax
538              ib(3)  = 1              ib(3)  = 1
539              ie(3)  = nlevels(listId)              ie(3)  = nLevOutp
540    
541  C           Time dimension  C           Time dimension
542              dn(4)(1:1) = 'T'              dn(4)(1:1) = 'T'
# Line 364  C           Time dimension Line 544  C           Time dimension
544              ib(4)  = 1              ib(4)  = 1
545              ie(4)  = 1              ie(4)  = 1
546    
547              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,
548       &             dim, dn, ib, ie, myThid)       &             dim, dn, ib, ie, myThid)
549              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,
550       &             4,5, myThid)       &             4,5, myThid)
551              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',
552       &             tdiag(ndId),myThid)       &             tdiag(ndId),myThid)
553              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',
554       &             udiag(ndId),myThid)       &             udiag(ndId),myThid)
555    
556              IF ((fflags(listId)(1:1) .EQ. ' ')  C     Missing values only for scalar diagnostics at mass points (so far)
557       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN              useMisValForThisDiag = useMissingValue
558         &           .AND.gdiag(ndId)(1:2).EQ.'SM'
559                IF ( useMisValForThisDiag ) THEN
560    C     assign missing values and set flag for adding the netCDF atttibute
561                 CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 2,
562         I            misval_r8, misval_r4, misval_int,
563         I            myThid )
564    C     and now use the missing values for masking out the land points
565    C     note: better to use 2-D mask if kdiag <> Nr or vert.integral
566                 DO bj = myByLo(myThid), myByHi(myThid)
567                  DO bi = myBxLo(myThid), myBxHi(myThid)
568                   DO k = 1,nLevOutp
569                    klev = NINT(levs(k,listId))
570                    IF ( fflags(listId)(2:2).EQ.'I' ) kLev = 1
571                    DO j = 1-OLy,sNy+OLy
572                     DO i = 1-OLx,sNx+OLx
573                      IF ( maskC(i,j,klev,bi,bj) .EQ. 0. )
574         &                 qtmp1(i,j,k,bi,bj) = misvalLoc
575                     ENDDO
576                    ENDDO
577                   ENDDO
578                  ENDDO
579                 ENDDO
580                ELSE
581    C     suppress the missing value attribute (iflag = 0)
582    C     Note: We have to call the following subroutine for each mnc that has
583    C     been created "on the fly" by mnc_cw_add_vname and will be deleted
584    C     by mnc_cw_del_vname, because all of these variables use the same
585    C     identifier so that mnc_cw_vfmv(indv) needs to be overwritten for
586    C     each of these variables
587                 CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 0,
588         I            misval_r8, misval_r4, misval_int,
589         I            myThid )
590                ENDIF
591    
592                IF (  ((writeBinaryPrec .EQ. precFloat32)
593         &            .AND. (fflags(listId)(1:1) .NE. 'D'))
594         &             .OR. (fflags(listId)(1:1) .EQ. 'R') ) THEN
595                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,
596       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
597              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              ELSEIF ( (writeBinaryPrec .EQ. precFloat64)
598         &             .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN
599                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
600       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
601              ENDIF              ENDIF
602                
603              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)
604              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)
605    
# Line 390  C           Time dimension Line 608  C           Time dimension
608            ENDIF            ENDIF
609  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
610    
611    C--      end loop on lm (or ll if ALLOW_MNC) counter
612             ENDDO
613  C--     end of Processing Fld # md  C--     end of Processing Fld # md
614          ENDIF          ENDIF
615           ENDDO
616    
617    #ifdef ALLOW_MNC
618    C--   end loop on jj counter
619        ENDDO        ENDDO
620    #endif
621    
622    #ifdef ALLOW_MDSIO
623          IF (diag_mdsio) THEN
624    C-    Note: temporary: since it is a pain to add more arguments to
625    C     all MDSIO S/R, uses instead this specific S/R to write only
626    C     meta files but with more informations in it.
627                glf = globalFiles
628                nRec = nfields(listId)*averageCycle(listId)
629                CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
630         &              0, 0, nLevOutp, ' ',
631         &              nfields(listId), flds(1,listId), nTimRec, timeRec,
632         &              nRec, myIter, myThid)
633          ENDIF
634    #endif /*  ALLOW_MDSIO  */
635    
636        RETURN        RETURN
637        END        END
Legend:
Removed from v.1.16  
changed lines
  Added in v.1.49
  ViewVC Help
Powered by ViewVC 1.1.22