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