/[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.63 by jmc, Fri Mar 24 23:34:13 2017 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, myTime, myIter, myThid )
      I     myIter,  
      I     myTime,  
      I     myThid )  
13    
14  C     !DESCRIPTION:  C     !DESCRIPTION:
15  C     Write output for diagnostics fields.  C     Write output for diagnostics fields.
# Line 26  C     !USES: Line 23  C     !USES:
23  #include "DIAGNOSTICS_SIZE.h"  #include "DIAGNOSTICS_SIZE.h"
24  #include "DIAGNOSTICS.h"  #include "DIAGNOSTICS.h"
25    
26  #ifdef ALLOW_FIZHI        INTEGER NrMax
27  #include "fizhi_SIZE.h"        PARAMETER( NrMax = numLevels )
 #else  
       INTEGER Nrphys  
       PARAMETER (Nrphys=0)  
 #endif  
   
28    
29  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
30  C     listId  :: Diagnostics list number being written  C     listId  :: Diagnostics list number being written
# Line 43  C     myThid  :: my Thread Id number Line 35  C     myThid  :: my Thread Id number
35        INTEGER listId, myIter, myThid        INTEGER listId, myIter, myThid
36  CEOP  CEOP
37    
38    C     !FUNCTIONS:
39          INTEGER ILNBLNK
40          EXTERNAL ILNBLNK
41    
42  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
43  C     i,j,k :: loop indices  C     i,j,k :: loop indices
44    C     bi,bj :: tile indices
45    C     lm    :: loop index (averageCycle)
46  C     md    :: field number in the list "listId".  C     md    :: field number in the list "listId".
47  C     ndId  :: diagnostics  Id number (in available diagnostics list)  C     ndId  :: diagnostics  Id number (in available diagnostics list)
 C     mate  :: counter mate Id number (in available diagnostics list)  
48  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
49  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
50        INTEGER i, j, k  C     mate  :: counter mate Id number (in available diagnostics list)
51    C     mDbl  :: processing mate Id number (in case processing requires 2 diags)
52    C     mVec  :: vector mate Id number
53    C     ppFld :: post-processed diag or not (=0): =1 stored in qtmp1 ; =2 in qtmp2
54    C   isComputed :: previous post-processed diag (still available in qtmp)
55    C     nLevOutp :: number of levels to write in output file
56    C
57    C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
58    C     qtmp1 :: temporary array; used to store a copy of diag. output field.
59    C     qtmp2 :: temporary array; used to store a copy of a 2nd diag. field.
60    C-  Note: local common block no longer needed.
61    c     COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
62          _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
63          _RL qtmp2(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, nn, ip, im
68        INTEGER mate, mVec        INTEGER mate, mDbl, mVec
69        CHARACTER*8 parms1        INTEGER ppFld, isComputed
70        CHARACTER*3 mate_index        CHARACTER*10 gcode
71        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)        _RL undefRL
72        _RL undef, getcon        INTEGER nLevOutp, kLev
       EXTERNAL getcon  
       INTEGER ILNBLNK  
       EXTERNAL ILNBLNK  
       INTEGER ilen  
       INTEGER nlevsout  
73    
74          INTEGER iLen,jLen
75        INTEGER ioUnit        INTEGER ioUnit
76        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
77        CHARACTER*(MAX_LEN_MBUF) suff        CHARACTER*(10) suff
78        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
79          INTEGER prec, nRec, nTimRec
80          _RL     timeRec(2)
81          _RL     tmpLoc
82    #ifdef ALLOW_MDSIO
83        LOGICAL glf        LOGICAL glf
84    #endif
85  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
       INTEGER ii  
86        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
       INTEGER CW_DIMS, NLEN  
       PARAMETER ( CW_DIMS = 10 )  
       PARAMETER ( NLEN    = 80 )  
       INTEGER dim(CW_DIMS), ib(CW_DIMS), ie(CW_DIMS)  
       CHARACTER*(NLEN) dn(CW_DIMS)  
       CHARACTER*(NLEN) d_cw_name  
       CHARACTER*(NLEN) dn_blnk  
 #ifdef DIAG_MNC_COORD_NEEDSWORK  
       CHARACTER*(5) ctmp  
       _RS ztmp(Nr+Nrphys)  
 #endif  
87  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
88    
89  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
90    
91    C---  set file properties
92        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
93        undef = getcon('UNDEF')        undefRL = misValFlt(listId)
       glf = globalFiles  
       WRITE(suff,'(I10.10)') myIter  
       ilen = ILNBLNK(fnames(listId))  
       WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)  
94    
95  #ifdef ALLOW_MNC        IF ( rwSuffixType.EQ.0 ) THEN
96        IF (useMNC .AND. diag_mnc) THEN          WRITE(suff,'(I10.10)') myIter
97          DO i = 1,MAX_LEN_FNAM        ELSE
98            diag_mnc_bn(i:i) = ' '          CALL RW_GET_SUFFIX( suff, myTime, myIter, myThid )
99          ENDDO        ENDIF
100          DO i = 1,NLEN        iLen = ILNBLNK(fnames(listId))
101            dn_blnk(i:i) = ' '        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff
102          ENDDO        IF ( diag_mdsio.AND.(diagMdsDir.NE.' ') ) THEN
103          WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)          jLen = ILNBLNK(diagMdsDir)
104            WRITE( fn, '(5A)' ) diagMdsDir(1:jLen),'/',
105  C       Update the record dimension by writing the iteration number       &                       fnames(listId)(1:iLen),'.',suff
106          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)        ENDIF
         CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',myTime,myThid)  
         CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)  
         CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',myIter,myThid)  
   
 C       NOTE: at some point it would be a good idea to add a time_bounds  
 C       variable that has dimension (2,T) and clearly denotes the  
 C       beginning and ending times for each diagnostics period  
   
         dn(1)(1:NLEN) = dn_blnk(1:NLEN)  
         WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)  
         dim(1) = nlevels(listId)  
         ib(1)  = 1  
         ie(1)  = nlevels(listId)  
   
         CALL MNC_CW_ADD_GNAME('diag_levels', 1,  
      &       dim, dn, ib, ie, myThid)  
         CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',  
      &       0,0, myThid)  
         CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',  
      &       'Idicies of vertical levels within the source arrays',  
      &       myThid)  
           
         CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,  
      &       'diag_levels', levs(1,listId), myThid)  
   
         CALL MNC_CW_DEL_VNAME('diag_levels', myThid)  
         CALL MNC_CW_DEL_GNAME('diag_levels', myThid)  
   
 #ifdef DIAG_MNC_COORD_NEEDSWORK  
 C       This part has been placed in an #ifdef because, as its currently  
 C       written, it will only work with variables defined on a dynamics  
 C       grid.  As we start using diagnostics for physics grids, ice  
 C       levels, land levels, etc. the different vertical coordinate  
 C       dimensions will have to be taken into account.  
   
 C       20051021 JMC & EH3 : We need to extend this so that a few  
 C       variables each defined on different grids do not have the same  
 C       vertical dimension names so we should be using a pattern such  
 C       as: Z[uml]td000000 where the 't' is the type as specified by  
 C       gdiag(10)  
   
 C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx  
         ctmp(1:5) = 'mul  '  
         DO i = 1,3  
           dn(1)(1:NLEN) = dn_blnk(1:NLEN)  
           WRITE(dn(1),'(3a,i6.6)') 'Z',ctmp(i:i),'d',nlevels(listId)  
           CALL MNC_CW_ADD_GNAME(dn(1), 1, dim, dn, ib, ie, myThid)  
           CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)  
   
 C         The following three ztmp() loops should eventually be modified  
 C         to reflect the fractional nature of levs(j,l) -- they should  
 C         do something like:  
 C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))  
 C                      + ( rC(INT(FLOOR(levs(j,l))))  
 C                          + rC(INT(CEIL(levs(j,l)))) )  
 C                        / ( levs(j,l) - FLOOR(levs(j,l)) )  
 C         for averaged levels.  
           IF (i .EQ. 1) THEN  
             DO j = 1,nlevels(listId)  
               ztmp(j) = rC(NINT(levs(j,listId)))  
             ENDDO  
             CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',  
      &           'Dimensional coordinate value at the mid point',  
      &           myThid)  
           ELSEIF (i .EQ. 2) THEN  
             DO j = 1,nlevels(listId)  
               ztmp(j) = rF(NINT(levs(j,listId)) + 1)  
             ENDDO  
             CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',  
      &           'Dimensional coordinate value at the upper point',  
      &           myThid)  
           ELSEIF (i .EQ. 3) THEN  
             DO j = 1,nlevels(listId)  
               ztmp(j) = rF(NINT(levs(j,listId)))  
             ENDDO  
             CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',  
      &           'Dimensional coordinate value at the lower point',  
      &           myThid)  
           ENDIF  
           CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)  
           CALL MNC_CW_DEL_VNAME(dn(1), myThid)  
           CALL MNC_CW_DEL_GNAME(dn(1), myThid)  
         ENDDO  
 #endif /*  DIAG_MNC_COORD_NEEDSWORK  */  
107    
108    C-    for now, if integrate vertically, output field has just 1 level:
109          nLevOutp = nlevels(listId)
110          IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
111    
112    C--   Set time information:
113          IF ( freq(listId).LT.0. ) THEN
114    C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
115            nTimRec = 1
116            timeRec(1) = myTime
117          ELSE
118    C-    Time-average: store the 2 edges of the time-averaging interval.
119    C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
120    C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
121            nTimRec = 2
122    
123    C-    end of time-averaging interval:
124            timeRec(2) = myTime
125    
126    C-    begining of time-averaging interval:
127    c       timeRec(1) = myTime - freq(listId)
128    C     a) find the time of the previous multiple of output freq:
129            timeRec(1) = myTime-deltaTClock*0.5 _d 0
130            timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
131            i = INT( timeRec(1) )
132            IF ( timeRec(1).LT.0. ) THEN
133              tmpLoc = FLOAT(i)
134              IF ( timeRec(1).NE.tmpLoc ) i = i - 1
135            ENDIF
136            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
137    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
138            timeRec(1) = MAX( timeRec(1), startTime )
139    
140    C     b) round off to nearest multiple of time-step:
141            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
142            i = NINT( timeRec(1) )
143    C     if just half way, NINT will return the next time-step: correct this
144            tmpLoc = FLOAT(i) - 0.5 _d 0
145            IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
146            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
147    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
148        ENDIF        ENDIF
149    C--   Convert time to iteration number (debug)
150    c     DO i=1,nTimRec
151    c       timeRec(i) = timeRec(i)/deltaTClock
152    c     ENDDO
153    
154    C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
155          DO lm=1,averageCycle(listId)
156    
157    #ifdef ALLOW_MNC
158           IF (useMNC .AND. diag_mnc) THEN
159             CALL DIAGNOSTICS_MNC_SET(
160         I                    nLevOutp, listId, lm,
161         O                    diag_mnc_bn,
162         I                    undefRL, myTime, myIter, myThid )
163           ENDIF
164  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
165    
166        DO md = 1,nfields(listId)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
167    
168           isComputed = 0
169           DO md = 1,nfields(listId)
170          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
171          parms1 = gdiag(ndId)(1:8)          gcode = gdiag(ndId)(1:10)
172          IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN          mate = 0
173            mVec = 0
174            mDbl = 0
175            ppFld = 0
176            IF ( gcode(5:5).EQ.'C' ) THEN
177    C-      Check for Mate of a Counter Diagnostic
178               mate = hdiag(ndId)
179            ELSEIF ( gcode(5:5).EQ.'P' ) THEN
180               ppFld = 1
181               IF ( gdiag(hdiag(ndId))(5:5).EQ.'P' ) ppFld = 2
182    C-      Also load the mate (if stored) for Post-Processing
183               nn = ndId
184               DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
185                 nn = hdiag(nn)
186               ENDDO
187               IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
188    c          write(0,*) ppFld,' ndId=', ndId, nn, mDbl, isComputed
189            ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
190    C-      Check for Mate of a Vector Diagnostic
191               mVec = hdiag(ndId)
192            ENDIF
193            IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
194  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
195    
196            ip = ABS(idiag(md,listId))            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
197            im = mdiag(md,listId)            im = mdiag(md,listId)
198            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
199              IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
200              IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
201    
202              IF ( ppFld.EQ.2 .AND. isComputed.EQ.hdiag(ndId) ) THEN
203    C-        Post-Processed diag from an other Post-Processed diag -and-
204    C         both of them have just been calculated and are still stored in qtmp:
205    C         => skip computation and just write qtmp2
206                IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
207                   WRITE(ioUnit,'(A,I6,3A,I6)')
208         &         '  get Post-Proc. Diag # ', ndId, '  ', cdiag(ndId),
209         &         ' from previous computation of Diag # ', isComputed
210                ENDIF
211                isComputed = 0
212              ELSEIF ( ndiag(ip,1,1).EQ.0 ) THEN
213  C-        Empty diagnostics case :  C-        Empty diagnostics case :
214                isComputed = 0
215    
216              _BEGIN_MASTER( myThid )              _BEGIN_MASTER( myThid )
217              WRITE(msgBuf,'(A,I10)')              WRITE(msgBuf,'(A,I10)')
218       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
219              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
220       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
221              WRITE(msgBuf,'(A,I4,3A,I3,2A)')              WRITE(msgBuf,'(A,I6,3A,I4,2A)')
222       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
223       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
224              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
225       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
226              WRITE(msgBuf,'(A,I2,A)')              IF ( averageCycle(listId).GT.1 ) THEN
227       &       '- WARNING -   has not been filled (ndiag=',               WRITE(msgBuf,'(A,2(I3,A))')
228       &       ndiag(ip,1,1), ' )'       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
229         &                                            ndiag(ip,1,1), ' )'
230                ELSE
231                 WRITE(msgBuf,'(A,2(I3,A))')
232         &        '- WARNING -   has not been filled (ndiag=',
233         &                                            ndiag(ip,1,1), ' )'
234                ENDIF
235              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
236       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
237              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
# Line 227  C-        Empty diagnostics case : Line 241  C-        Empty diagnostics case :
241              _END_MASTER( myThid )              _END_MASTER( myThid )
242              DO bj = myByLo(myThid), myByHi(myThid)              DO bj = myByLo(myThid), myByHi(myThid)
243                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
244                  DO k = 1,nlevels(listId)                  DO k = 1,nLevOutp
245                    DO j = 1-OLy,sNy+OLy                    DO j = 1-OLy,sNy+OLy
246                      DO i = 1-OLx,sNx+OLx                      DO i = 1-OLx,sNx+OLx
247                        qtmp1(i,j,k,bi,bj) = 0. _d 0                        qtmp1(i,j,k,bi,bj) = 0. _d 0
# Line 239  C-        Empty diagnostics case : Line 253  C-        Empty diagnostics case :
253    
254            ELSE            ELSE
255  C-        diagnostics is not empty :  C-        diagnostics is not empty :
256                isComputed = 0
257    
258              IF ( myThid.EQ.1 ) WRITE(ioUnit,'(A,I3,3A,I8,2A)')              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
259                  IF ( ppFld.GE.1 ) THEN
260                   WRITE(ioUnit,'(A,I6,7A,I8,2A)')
261         &         ' Post-Processing Diag # ', ndId, '  ', cdiag(ndId),
262         &         '   Parms: ',gdiag(ndId)
263                   IF ( mDbl.EQ.0 ) THEN
264                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
265         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
266                   ELSE
267                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
268         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
269         &          ' and diag: ',
270         &            cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
271                   ENDIF
272                  ELSE
273                   WRITE(ioUnit,'(A,I6,3A,I8,2A)')
274       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
275       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
276                  ENDIF
277              IF ( parms1(5:5).EQ.'C' ) THEN                IF ( mate.GT.0 ) THEN
278  C             Check for Mate of a Counter Diagnostic                 WRITE(ioUnit,'(3A,I6,2A)')
 C             --------------------------------------  
               mate_index = parms1(6:8)  
               READ (mate_index,'(I3)') mate  
               IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,2A)')  
279       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
280       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
281                  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  
282                  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
283                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
284       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
285       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
286       &             ' exists '       &             ' exists '
287                  ELSE                  ELSE
288                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
289       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
290       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
291       &             ' not enabled'       &             ' not enabled'
# Line 275  C             -------------------------- Line 293  C             --------------------------
293                ENDIF                ENDIF
294              ENDIF              ENDIF
295    
296              DO bj = myByLo(myThid), myByHi(myThid)              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
297               DO bi = myBxLo(myThid), myBxHi(myThid)  C-       get only selected levels:
298                DO k = 1,nlevels(listId)                DO bj = myByLo(myThid), myByHi(myThid)
299                  CALL GETDIAG(                 DO bi = myBxLo(myThid), myBxHi(myThid)
300       I                       levs(k,listId),undef,                  DO k = 1,nlevels(listId)
301       O                       qtmp1(1-OLx,1-OLy,k,bi,bj),                    kLev = NINT(levs(k,listId))
302       I                       ndId,mate,ip,im,bi,bj,myThid)                    CALL DIAGNOSTICS_GET_DIAG(
303         I                         kLev, undefRL,
304         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
305         I                         ndId, mate, ip, im, bi, bj, myThid )
306                    ENDDO
307                   ENDDO
308                ENDDO                ENDDO
309               ENDDO                IF ( mDbl.GT.0 ) THEN
310              ENDDO                 DO bj = myByLo(myThid), myByHi(myThid)
311                    DO bi = myBxLo(myThid), myBxHi(myThid)
312  C-        end of empty diag / not empty block                   DO k = 1,nlevels(listId)
313            ENDIF                    kLev = NINT(levs(k,listId))
314                      CALL DIAGNOSTICS_GET_DIAG(
315            nlevsout = nlevels(listId)       I                         kLev, undefRL,
316         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
317  C-----------------------------------------------------------------------       I                         mDbl, 0, im, 0, bi, bj, myThid )
318  C         Check to see if we need to interpolate before output                   ENDDO
319  C-----------------------------------------------------------------------                  ENDDO
320           IF ( fflags(listId)(2:2).EQ.'P' ) THEN                 ENDDO
321  C-        Do vertical interpolation:                ENDIF
           CALL DIAGNOSTICS_INTERP_VERT(  
      I                     listId, md, ndId, ip, im,  
      U                     nlevsout,  
      U                     qtmp1,  
      I                     undef,  
      I                     myTime, myIter, myThid )  
          ENDIF  
   
 #ifdef ALLOW_MDSIO  
 C         Prepare for mdsio optionality  
           IF (diag_mdsio) THEN  
             IF (fflags(listId)(1:1) .EQ. 'R') THEN  
 C             Force it to be 32-bit precision  
               CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,.FALSE.,  
      &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)  
             ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN  
 C             Force it to be 64-bit precision  
               CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,.FALSE.,  
      &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)  
322              ELSE              ELSE
323  C             This is the old default behavior  C-       get all the levels (for vertical post-processing)
324                CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,.FALSE.,                DO bj = myByLo(myThid), myByHi(myThid)
325       &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
326                      CALL DIAGNOSTICS_GET_DIAG(
327         I                         0, undefRL,
328         O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
329         I                         ndId, mate, ip, im, bi, bj, myThid )
330                   ENDDO
331                  ENDDO
332                  IF ( mDbl.GT.0 ) THEN
333                   DO bj = myByLo(myThid), myByHi(myThid)
334                    DO bi = myBxLo(myThid), myBxHi(myThid)
335                      CALL DIAGNOSTICS_GET_DIAG(
336         I                         0, undefRL,
337         O                         qtmp2(1-OLx,1-OLy,1,bi,bj),
338         I                         mDbl, 0, im, 0, bi, bj, myThid )
339                    ENDDO
340                   ENDDO
341                  ENDIF
342              ENDIF              ENDIF
           ENDIF  
 #endif /*  ALLOW_MDSIO  */  
343    
344  #ifdef ALLOW_MNC  C-----------------------------------------------------------------------
345            IF (useMNC .AND. diag_mnc) THEN  C--     Apply specific post-processing (e.g., interpolate) before output
346    C-----------------------------------------------------------------------
347              _BEGIN_MASTER( myThid )              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
348    C-          Do vertical interpolation:
349              DO ii = 1,CW_DIMS               IF ( fluidIsAir ) THEN
350                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
351                dn(ii)(1:NLEN) = dn_blnk(1:NLEN)                CALL DIAGNOSTICS_INTERP_VERT(
352              ENDDO       I                         listId, md, ndId, ip, im, lm,
353         U                         qtmp1, qtmp2,
354  C           Note that the "d_cw_name" variable is a hack that hides a       I                         undefRL, myTime, myIter, myThid )
355  C           subtlety within MNC.  Basically, each MNC-wrapped file is               ELSE
356  C           caching its own concept of what each "grid name" (that is, a                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
357  C           dimension group name) means.  So one cannot re-use the same       &           'INTERP_VERT not allowed in this config'
358  C           "grid" name for different collections of dimensions within a                 CALL PRINT_ERROR( msgBuf , myThid )
359  C           given file.  By appending the "ndId" values to each name, we                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
360  C           guarantee uniqueness within each MNC-produced file.               ENDIF
             WRITE(d_cw_name,'(a,i6.6)') 'd_cw_',ndId  
   
 C           XY dimensions  
             dim(1)       = sNx + 2*OLx  
             dim(2)       = sNy + 2*OLy  
             ib(1)        = OLx + 1  
             ib(2)        = OLy + 1  
             IF (gdiag(ndId)(2:2) .EQ. 'M') THEN  
               dn(1)(1:2) = 'X'  
               ie(1)      = OLx + sNx  
               dn(2)(1:2) = 'Y'  
               ie(2)      = OLy + sNy  
             ELSEIF (gdiag(ndId)(2:2) .EQ. 'U') THEN  
               dn(1)(1:3) = 'Xp1'  
               ie(1)      = OLx + sNx + 1  
               dn(2)(1:2) = 'Y'  
               ie(2)      = OLy + sNy  
             ELSEIF (gdiag(ndId)(2:2) .EQ. 'V') THEN  
               dn(1)(1:2) = 'X'  
               ie(1)      = OLx + sNx  
               dn(2)(1:3) = 'Yp1'  
               ie(2)      = OLy + sNy + 1  
             ELSEIF (gdiag(ndId)(2:2) .EQ. 'Z') THEN  
               dn(1)(1:3) = 'Xp1'  
               ie(1)      = OLx + sNx + 1  
               dn(2)(1:3) = 'Yp1'  
               ie(2)      = OLy + sNy + 1  
             ENDIF  
               
 C           Z is special since it varies  
             WRITE(dn(3),'(a,i6.6)') 'Zd', nlevsout  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevsout  
             ENDIF  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zld', nlevsout  
361              ENDIF              ENDIF
362              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( fflags(listId)(2:2).EQ.'I' ) THEN
363       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN  C-          Integrate vertically: for now, output field has just 1 level:
364                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevsout                CALL DIAGNOSTICS_SUM_LEVELS(
365         I                         listId, md, ndId, ip, im, lm,
366         U                         qtmp1,
367         I                         undefRL, myTime, myIter, myThid )
368              ENDIF              ENDIF
369              dim(3) = Nr+Nrphys              IF ( ppFld.GE.1 ) THEN
370              ib(3)  = 1  C-          Do Post-Processing:
371              ie(3)  = nlevsout               IF ( flds(md,listId).EQ.'PhiVEL  '
372         &       .OR. flds(md,listId).EQ.'PsiVEL  '
373  C           Time dimension       &          ) THEN
374              dn(4)(1:1) = 'T'                CALL DIAGNOSTICS_CALC_PHIVEL(
375              dim(4) = -1       I                         listId, md, ndId, ip, im, lm,
376              ib(4)  = 1       I                         NrMax,
377              ie(4)  = 1       U                         qtmp1, qtmp2,
378         I                         myTime, myIter, myThid )
379              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,                isComputed = ndId
380       &             dim, dn, ib, ie, myThid)               ELSE
381              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
382       &             4,5, myThid)       &           'unknown Processing method for diag="',cdiag(ndId),'"'
383              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',                 CALL PRINT_ERROR( msgBuf , myThid )
384       &             tdiag(ndId),myThid)                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
385              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',               ENDIF
      &             udiag(ndId),myThid)  
   
 C           Per the observations of Baylor, this has been commented out  
 C           until we have code that can write missing_value attributes  
 C           in a way thats compatible with most of the more popular  
 C           netCDF tools including ferret.  Using all-zeros completely  
 C           breaks ferret.  
   
 C           CALL MNC_CW_ADD_VATTR_DBL(cdiag(ndId),'missing_value',  
 C           &             0.0 _d 0,myThid)  
   
             IF ( ( (writeBinaryPrec .EQ. precFloat32)  
      &           .AND. (fflags(listId)(1:1) .NE. 'D')  
      &           .AND. (fflags(listId)(1:1) .NE. 'R') )  
      &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN  
               CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,  
      &             cdiag(ndId), qtmp1, myThid)  
             ELSEIF ( (writeBinaryPrec .EQ. precFloat64)  
      &             .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN  
               CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,  
      &             cdiag(ndId), qtmp1, myThid)  
386              ENDIF              ENDIF
               
             CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)  
             CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)  
387    
388              _END_MASTER( myThid )  C--     End of empty diag / not-empty diag block
389              ENDIF
390    
391    C--     Ready to write field "md", element "lm" in averageCycle(listId)
392    
393    C-        write to binary file, using MDSIO pkg:
394              IF ( diag_mdsio ) THEN
395    c          nRec = lm + (md-1)*averageCycle(listId)
396               nRec = md + (lm-1)*nfields(listId)
397    C         default precision for output files
398               prec = writeBinaryPrec
399    C         fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
400               IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
401               IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
402    C         a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R
403               IF ( ppFld.LE.1 ) THEN
404                CALL WRITE_REC_LEV_RL(
405         I                            fn, prec,
406         I                            NrMax, 1, nLevOutp,
407         I                            qtmp1, -nRec, myIter, myThid )
408               ELSE
409                CALL WRITE_REC_LEV_RL(
410         I                            fn, prec,
411         I                            NrMax, 1, nLevOutp,
412         I                            qtmp2, -nRec, myIter, myThid )
413               ENDIF
414              ENDIF
415    
416    #ifdef ALLOW_MNC
417              IF (useMNC .AND. diag_mnc) THEN
418               IF ( ppFld.LE.1 ) THEN
419                CALL DIAGNOSTICS_MNC_OUT(
420         I                       NrMax, nLevOutp, listId, ndId, mate,
421         I                       diag_mnc_bn, qtmp1,
422         I                       undefRL, myTime, myIter, myThid )
423               ELSE
424                CALL DIAGNOSTICS_MNC_OUT(
425         I                       NrMax, nLevOutp, listId, ndId, mate,
426         I                       diag_mnc_bn, qtmp2,
427         I                       undefRL, myTime, myIter, myThid )
428               ENDIF
429            ENDIF            ENDIF
430  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
431    
432  C--     end of Processing Fld # md  C--     end of Processing Fld # md
433          ENDIF          ENDIF
434           ENDDO
435    
436    C--   end loop on lm counter (= averagePeriod)
437        ENDDO        ENDDO
438    
439    #ifdef ALLOW_MDSIO
440          IF (diag_mdsio) THEN
441    C-    Note: temporary: since it is a pain to add more arguments to
442    C     all MDSIO S/R, uses instead this specific S/R to write only
443    C     meta files but with more informations in it.
444                glf = globalFiles
445                nRec = averageCycle(listId)*nfields(listId)
446                CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
447         &              0, 0, nLevOutp, ' ',
448         &              nfields(listId), flds(1,listId),
449         &              nTimRec, timeRec, undefRL,
450         &              nRec, myIter, myThid)
451          ENDIF
452    #endif /*  ALLOW_MDSIO  */
453    
454        RETURN        RETURN
455        END        END
456    
Legend:
Removed from v.1.28  
changed lines
  Added in v.1.63
  ViewVC Help
Powered by ViewVC 1.1.22