/[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.25 by edhill, Fri Nov 25 06:12:31 2005 UTC revision 1.45 by jmc, Thu Jan 7 23:48:01 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. ) i = i - 1
139            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
140    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
141            timeRec(1) = MAX( timeRec(1), startTime )
142    
143    C     b) round off to nearest multiple of time-step:
144            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
145            i = NINT( timeRec(1) )
146    C     if just half way, NINT will return the next time-step: correct this
147            tmpLev = FLOAT(i) - 0.5 _d 0
148            IF ( timeRec(1).EQ.tmpLev ) i = i - 1
149            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
150    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
151          ENDIF
152    
153  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
154    C-- this is a trick to reverse the order of the loops on md (= field)
155    C   and lm (= averagePeriod): binary output: lm loop inside md loop ;
156    C                                 mnc ouput: md loop inside lm loop.
157        IF (useMNC .AND. diag_mnc) THEN        IF (useMNC .AND. diag_mnc) THEN
158          DO i = 1,MAX_LEN_FNAM          jjMx = averageCycle(listId)
159            diag_mnc_bn(i:i) = ' '          llMx = 1
160          ENDDO        ELSE
161          DO i = 1,NLEN          jjMx = 1
162            dn_blnk(i:i) = ' '          llMx = averageCycle(listId)
163          ENDDO        ENDIF
164          WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)        DO jj=1,jjMx
165    
166           IF (useMNC .AND. diag_mnc) THEN
167    C     Handle missing value attribute (land points)
168             useMissingValue = .FALSE.
169    #ifdef DIAGNOSTICS_MISSING_VALUE
170             useMissingValue = .TRUE.
171    #endif /* DIAGNOSTICS_MISSING_VALUE */
172             IF ( misvalFlt(listId) .NE. UNSET_RL ) THEN
173              misvalLoc = misvalFlt(listId)
174             ELSE
175              misvalLoc = undef
176             ENDIF
177    C     Defaults to UNSET_I
178             misvalIntLoc = misvalInt(listId)
179             DO ii=1,2
180    C         misval_r4(ii)  = UNSET_FLOAT4
181    C         misval_r8(ii)  = UNSET_FLOAT8
182              misval_r4(ii)  = misvalLoc
183              misval_r8(ii)  = misvalLoc
184              misval_int(ii) = UNSET_I
185             ENDDO
186             DO i = 1,MAX_LEN_FNAM
187               diag_mnc_bn(i:i) = ' '
188             ENDDO
189             DO i = 1,NLEN
190               dn_blnk(i:i) = ' '
191             ENDDO
192             WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)
193    
194  C       Update the record dimension by writing the iteration number  C       Update the record dimension by writing the iteration number
195          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)           klev = myIter + jj - jjMx
196          CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',myTime,myThid)           tmpLev = myTime + deltaTClock*(jj -jjMx)
197          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)           CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)
198          CALL MNC_CW_I_W_S('I',diag_mnc_bn,1,1,'iter',myIter,myThid)           CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',tmpLev,myThid)
199             CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)
200             CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',klev,myThid)
201    
202  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
203  C       variable that has dimension (2,T) and clearly denotes the  C       variable that has dimension (2,T) and clearly denotes the
204  C       beginning and ending times for each diagnostics period  C       beginning and ending times for each diagnostics period
205    
206          dn(1)(1:NLEN) = dn_blnk(1:NLEN)           dn(1)(1:NLEN) = dn_blnk(1:NLEN)
207          WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)           WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)
208          dim(1) = nlevels(listId)           dim(1) = nlevels(listId)
209          ib(1)  = 1           ib(1)  = 1
210          ie(1)  = nlevels(listId)           ie(1)  = nlevels(listId)
211    
212          CALL MNC_CW_ADD_GNAME('diag_levels', 1,           CALL MNC_CW_ADD_GNAME('diag_levels', 1,
213       &       dim, dn, ib, ie, myThid)       &        dim, dn, ib, ie, myThid)
214          CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',           CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',
215       &       0,0, myThid)       &        0,0, myThid)
216          CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',           CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',
217       &       'Idicies of vertical levels within the source arrays',       &        'Idicies of vertical levels within the source arrays',
218       &       myThid)       &        myThid)
219            C     suppress the missing value attribute (iflag = 0)
220          CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,           IF (useMissingValue)
221       &       'diag_levels', levs(1,listId), myThid)       &       CALL MNC_CW_VATTR_MISSING('diag_levels', 0,
222         I       misval_r8, misval_r4, misval_int,
223         I       myThid )
224    
225          CALL MNC_CW_DEL_VNAME('diag_levels', myThid)           CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
226          CALL MNC_CW_DEL_GNAME('diag_levels', myThid)       &        'diag_levels', levs(1,listId), myThid)
227    
228             CALL MNC_CW_DEL_VNAME('diag_levels', myThid)
229             CALL MNC_CW_DEL_GNAME('diag_levels', myThid)
230    
231  #ifdef DIAG_MNC_COORD_NEEDSWORK  #ifdef DIAG_MNC_COORD_NEEDSWORK
232  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 242  C       as: Z[uml]td000000 where the 't'
242  C       gdiag(10)  C       gdiag(10)
243    
244  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx
245          ctmp(1:5) = 'mul  '           ctmp(1:5) = 'mul  '
246          DO i = 1,3           DO i = 1,3
247            dn(1)(1:NLEN) = dn_blnk(1:NLEN)             dn(1)(1:NLEN) = dn_blnk(1:NLEN)
248            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)
249            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)
250            CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)             CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)
251    
252  C         The following three ztmp() loops should eventually be modified  C         The following three ztmp() loops should eventually be modified
253  C         to reflect the fractional nature of levs(j,l) -- they should  C         to reflect the fractional nature of levs(j,l) -- they should
254  C         do something like:  C         do something like:
255  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))
256  C                      + ( rC(INT(FLOOR(levs(j,l))))  C                      + ( rC(INT(FLOOR(levs(j,l))))
257  C                          + rC(INT(CEIL(levs(j,l)))) )  C                          + rC(INT(CEIL(levs(j,l)))) )
258  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )
259  C         for averaged levels.  C         for averaged levels.
260            IF (i .EQ. 1) THEN             IF (i .EQ. 1) THEN
261              DO j = 1,nlevels(listId)               DO j = 1,nlevels(listId)
262                ztmp(j) = rC(NINT(levs(j,listId)))                 ztmp(j) = rC(NINT(levs(j,listId)))
263              ENDDO               ENDDO
264              CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',               CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
265       &           'Dimensional coordinate value at the mid point',       &            'Dimensional coordinate value at the mid point',
266       &           myThid)       &            myThid)
267            ELSEIF (i .EQ. 2) THEN             ELSEIF (i .EQ. 2) THEN
268              DO j = 1,nlevels(listId)               DO j = 1,nlevels(listId)
269                ztmp(j) = rF(NINT(levs(j,listId)) + 1)                 ztmp(j) = rF(NINT(levs(j,listId)) + 1)
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 upper point',       &            'Dimensional coordinate value at the upper point',
273       &           myThid)       &            myThid)
274            ELSEIF (i .EQ. 3) THEN             ELSEIF (i .EQ. 3) THEN
275              DO j = 1,nlevels(listId)               DO j = 1,nlevels(listId)
276                ztmp(j) = rF(NINT(levs(j,listId)))                 ztmp(j) = rF(NINT(levs(j,listId)))
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 lower point',       &            'Dimensional coordinate value at the lower point',
280       &           myThid)       &            myThid)
281            ENDIF             ENDIF
282            CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)  C     suppress the missing value attribute (iflag = 0)
283            CALL MNC_CW_DEL_VNAME(dn(1), myThid)             IF (useMissingValue)
284            CALL MNC_CW_DEL_GNAME(dn(1), myThid)       &          CALL MNC_CW_VATTR_MISSING(dn(1), 0,
285          ENDDO       I          misval_r8, misval_r4, misval_int,
286         I          myThid )
287               CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)
288               CALL MNC_CW_DEL_VNAME(dn(1), myThid)
289               CALL MNC_CW_DEL_GNAME(dn(1), myThid)
290             ENDDO
291  #endif /*  DIAG_MNC_COORD_NEEDSWORK  */  #endif /*  DIAG_MNC_COORD_NEEDSWORK  */
292    
293        ENDIF         ENDIF
294  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
295    
296        DO md = 1,nfields(listId)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
297    
298           DO md = 1,nfields(listId)
299          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
300          parms1 = gdiag(ndId)(1:8)          gcode = gdiag(ndId)(1:10)
301          IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN          mate = 0
302            mVec = 0
303            IF ( gcode(5:5).EQ.'C' ) THEN
304    C-      Check for Mate of a Counter Diagnostic
305               mate = hdiag(ndId)
306            ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
307    C-      Check for Mate of a Vector Diagnostic
308               mVec = hdiag(ndId)
309            ENDIF
310            IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
311  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
312    #ifdef ALLOW_MNC
313             DO ll=1,llMx
314              lm = jj+ll-1
315    #else
316             DO lm=1,averageCycle(listId)
317    #endif
318    
319            ip = ABS(idiag(md,listId))            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
320            im = mdiag(md,listId)            im = mdiag(md,listId)
321              IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
322              IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
323    
324            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ndiag(ip,1,1).EQ.0 ) THEN
325  C-        Empty diagnostics case :  C-        Empty diagnostics case :
326    
# Line 210  C-        Empty diagnostics case : Line 329  C-        Empty diagnostics case :
329       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
330              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
331       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
332              WRITE(msgBuf,'(A,I4,3A,I3,2A)')              WRITE(msgBuf,'(A,I6,3A,I4,2A)')
333       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
334       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
335              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
336       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
337              WRITE(msgBuf,'(A,I2,A)')              IF ( averageCycle(listId).GT.1 ) THEN
338       &       '- WARNING -   has not been filled (ndiag=',               WRITE(msgBuf,'(A,2(I3,A))')
339       &       ndiag(ip,1,1), ' )'       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
340         &                                            ndiag(ip,1,1), ' )'
341                ELSE
342                 WRITE(msgBuf,'(A,2(I3,A))')
343         &        '- WARNING -   has not been filled (ndiag=',
344         &                                            ndiag(ip,1,1), ' )'
345                ENDIF
346              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
347       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
348              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
# Line 240  C-        Empty diagnostics case : Line 365  C-        Empty diagnostics case :
365            ELSE            ELSE
366  C-        diagnostics is not empty :  C-        diagnostics is not empty :
367    
368              IF ( myThid.EQ.1 ) WRITE(ioUnit,'(A,I3,3A,I8,2A)')              IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN
369                  WRITE(ioUnit,'(A,I6,3A,I8,2A)')
370       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
371       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
372                  IF ( mate.GT.0 ) THEN
373              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)')  
374       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
375       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
376                  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  
377                  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
378                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
379       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
380       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
381       &             ' exists '       &             ' exists '
382                  ELSE                  ELSE
383                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
384       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
385       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
386       &             ' not enabled'       &             ' not enabled'
# Line 275  C             -------------------------- Line 388  C             --------------------------
388                ENDIF                ENDIF
389              ENDIF              ENDIF
390    
391              DO bj = myByLo(myThid), myByHi(myThid)              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
392               DO bi = myBxLo(myThid), myBxHi(myThid)  C-       get all the levels (for vertical interpolation)
393                DO k = 1,nlevels(listId)                DO bj = myByLo(myThid), myByHi(myThid)
394                  CALL GETDIAG(                 DO bi = myBxLo(myThid), myBxHi(myThid)
395       I                       levs(k,listId),undef,                  DO k = 1,kdiag(ndId)
396       O                       qtmp1(1-OLx,1-OLy,k,bi,bj),                    tmpLev = k
397       I                       ndId,mate,ip,im,bi,bj,myThid)                    CALL GETDIAG(
398         I                         tmpLev,undef,
399         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
400         I                         ndId,mate,ip,im,bi,bj,myThid)
401                    ENDDO
402                   ENDDO
403                ENDDO                ENDDO
404               ENDDO              ELSE
405              ENDDO  C-       get only selected levels:
406                  DO bj = myByLo(myThid), myByHi(myThid)
407                   DO bi = myBxLo(myThid), myBxHi(myThid)
408                    DO k = 1,nlevels(listId)
409                      CALL GETDIAG(
410         I                         levs(k,listId),undef,
411         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
412         I                         ndId,mate,ip,im,bi,bj,myThid)
413                    ENDDO
414                   ENDDO
415                  ENDDO
416                ENDIF
417    
418  C-        end of empty diag / not empty block  C-        end of empty diag / not empty block
419            ENDIF            ENDIF
420    
           nlevsout = nlevels(listId)  
   
421  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
422  C         Check to see if we need to interpolate before output  C         Check to see if we need to interpolate before output
423  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
424           IF ( fflags(listId)(2:2).EQ.'P' ) THEN            IF ( fflags(listId)(2:2).EQ.'P' ) THEN
425  C-        Do vertical interpolation:  C-        Do vertical interpolation:
426            CALL DIAGNOSTICS_INTERP_VERT(             IF ( fluidIsAir ) THEN
427       I                     listId, md, ndId, ip, im,  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
428       U                     nlevsout,              CALL DIAGNOSTICS_INTERP_VERT(
429         I                     listId, md, ndId, ip, im, lm,
430       U                     qtmp1,       U                     qtmp1,
431       I                     undef,       I                     undef, myTime, myIter, myThid )
432       I                     myTime, myIter, myThid )             ELSE
433           ENDIF               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
434         &         'INTERP_VERT not allowed in this config'
435                 CALL PRINT_ERROR( msgBuf , myThid )
436                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
437               ENDIF
438              ENDIF
439    
440  #ifdef ALLOW_MDSIO  C--    Ready to write field "md", element "lm" in averageCycle(listId)
441  C         Prepare for mdsio optionality  
442            IF (diag_mdsio) THEN  C-        write to binary file, using MDSIO pkg:
443              IF (fflags(listId)(1:1) .EQ. 'R') THEN            IF ( diag_mdsio ) THEN
444  C             Force it to be 32-bit precision              nRec = lm + (md-1)*averageCycle(listId)
445                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,.FALSE.,  C           default precision for output files
446       &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)              prec = writeBinaryPrec
447              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
448  C             Force it to be 64-bit precision              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
449                CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,.FALSE.,              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
450       &             '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
451              ELSE              CALL WRITE_REC_LEV_RL(
452  C             This is the old default behavior       I                            fn, prec,
453                CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,.FALSE.,       I                            NrMax, 1, nlevels(listId),
454       &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       I                            qtmp1, -nRec, myIter, myThid )
             ENDIF  
455            ENDIF            ENDIF
 #endif /*  ALLOW_MDSIO  */  
456    
457  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
458            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
# Line 347  C           XY dimensions Line 478  C           XY dimensions
478              dim(2)       = sNy + 2*OLy              dim(2)       = sNy + 2*OLy
479              ib(1)        = OLx + 1              ib(1)        = OLx + 1
480              ib(2)        = OLy + 1              ib(2)        = OLy + 1
481              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN
482                dn(1)(1:2) = 'X'                dn(1)(1:2) = 'X'
483                ie(1)      = OLx + sNx                ie(1)      = OLx + sNx
484                dn(2)(1:2) = 'Y'                dn(2)(1:2) = 'Y'
# Line 368  C           XY dimensions Line 499  C           XY dimensions
499                dn(2)(1:3) = 'Yp1'                dn(2)(1:3) = 'Yp1'
500                ie(2)      = OLy + sNy + 1                ie(2)      = OLy + sNy + 1
501              ENDIF              ENDIF
502                
503  C           Z is special since it varies  C           Z is special since it varies
504              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevsout              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)
505              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
506       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN
507                WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevsout                WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevels(listId)
508              ENDIF              ENDIF
509              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
510       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN
511                WRITE(dn(3),'(a,i6.6)') 'Zld', nlevsout                WRITE(dn(3),'(a,i6.6)') 'Zld', nlevels(listId)
512              ENDIF              ENDIF
513              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
514       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
515                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevsout                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)
516              ENDIF              ENDIF
517              dim(3) = Nr+Nrphys              dim(3) = NrMax
518              ib(3)  = 1              ib(3)  = 1
519              ie(3)  = nlevsout              ie(3)  = nlevels(listId)
520    
521  C           Time dimension  C           Time dimension
522              dn(4)(1:1) = 'T'              dn(4)(1:1) = 'T'
# Line 393  C           Time dimension Line 524  C           Time dimension
524              ib(4)  = 1              ib(4)  = 1
525              ie(4)  = 1              ie(4)  = 1
526    
527              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,
528       &             dim, dn, ib, ie, myThid)       &             dim, dn, ib, ie, myThid)
529              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,
530       &             4,5, myThid)       &             4,5, myThid)
531              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',
532       &             tdiag(ndId),myThid)       &             tdiag(ndId),myThid)
533              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',
534       &             udiag(ndId),myThid)       &             udiag(ndId),myThid)
535    
536              IF ( ( (writeBinaryPrec .EQ. precFloat32)  C     Missing values only for scalar diagnostics at mass points (so far)
537       &           .AND. (fflags(listId)(1:1) .NE. 'D')              useMisValForThisDiag = useMissingValue
538       &           .AND. (fflags(listId)(1:1) .NE. 'R') )       &           .AND.gdiag(ndId)(1:2).EQ.'SM'
539       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN              IF ( useMisValForThisDiag ) THEN
540    C     assign missing values and set flag for adding the netCDF atttibute
541                 CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 2,
542         I            misval_r8, misval_r4, misval_int,
543         I            myThid )
544    C     and now use the missing values for masking out the land points
545                 DO bj = myByLo(myThid), myByHi(myThid)
546                  DO bi = myBxLo(myThid), myBxHi(myThid)
547                   DO k = 1,nlevels(listId)
548                    klev = NINT(levs(k,listId))
549                    DO j = 1-OLy,sNy+OLy
550                     DO i = 1-OLx,sNx+OLx
551                      IF ( maskC(i,j,klev,bi,bj) .EQ. 0. )
552         &                 qtmp1(i,j,k,bi,bj) = misvalLoc
553                     ENDDO
554                    ENDDO
555                   ENDDO
556                  ENDDO
557                 ENDDO
558                ELSE
559    C     suppress the missing value attribute (iflag = 0)
560    C     Note: We have to call the following subroutine for each mnc that has
561    C     been created "on the fly" by mnc_cw_add_vname and will be deleted
562    C     by mnc_cw_del_vname, because all of these variables use the same
563    C     identifier so that mnc_cw_vfmv(indv) needs to be overwritten for
564    C     each of these variables
565                 CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 0,
566         I            misval_r8, misval_r4, misval_int,
567         I            myThid )
568                ENDIF
569    
570                IF (  ((writeBinaryPrec .EQ. precFloat32)
571         &            .AND. (fflags(listId)(1:1) .NE. 'D'))
572         &             .OR. (fflags(listId)(1:1) .EQ. 'R') ) THEN
573                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,
574       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
575              ELSEIF ( (writeBinaryPrec .EQ. precFloat64)              ELSEIF ( (writeBinaryPrec .EQ. precFloat64)
576       &             .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN       &             .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN
577                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
578       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
579              ENDIF              ENDIF
580                
581              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)
582              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)
583    
# Line 422  C           Time dimension Line 586  C           Time dimension
586            ENDIF            ENDIF
587  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
588    
589    C--      end loop on lm (or ll if ALLOW_MNC) counter
590             ENDDO
591  C--     end of Processing Fld # md  C--     end of Processing Fld # md
592          ENDIF          ENDIF
593           ENDDO
594    
595    #ifdef ALLOW_MNC
596    C--   end loop on jj counter
597        ENDDO        ENDDO
598    #endif
599    
600    #ifdef ALLOW_MDSIO
601          IF (diag_mdsio) THEN
602    C-    Note: temporary: since it's a pain to add more arguments to
603    C     all MDSIO S/R, uses instead this specific S/R to write only
604    C     meta files but with more informations in it.
605                glf = globalFiles
606                nRec = nfields(listId)*averageCycle(listId)
607                CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
608         &              0, 0, nlevels(listId), ' ',
609         &              nfields(listId), flds(1,listId), nTimRec, timeRec,
610         &              nRec, myIter, myThid)
611          ENDIF
612    #endif /*  ALLOW_MDSIO  */
613    
614        RETURN        RETURN
615        END        END
Legend:
Removed from v.1.25  
changed lines
  Added in v.1.45
  ViewVC Help
Powered by ViewVC 1.1.22