/[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.34 by jmc, Tue Nov 13 19:43:44 2007 UTC revision 1.61 by jmc, Wed Feb 6 21:25:26 2013 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,
      I     myIter,  
13       I     myTime,       I     myTime,
14         I     myIter,
15       I     myThid )       I     myThid )
16    
17  C     !DESCRIPTION:  C     !DESCRIPTION:
# Line 38  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    
45  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
46  C     i,j,k :: loop indices  C     i,j,k :: loop indices
47    C     bi,bj :: tile indices
48  C     lm    :: loop index (averageCycle)  C     lm    :: loop index (averageCycle)
49  C     md    :: field number in the list "listId".  C     md    :: field number in the list "listId".
50  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)  
51  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
52  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
53    C     mate  :: counter mate Id number (in available diagnostics list)
54    C     mDbl  :: processing mate Id number (in case processing requires 2 diags)
55    C     mVec  :: vector mate Id number
56    C     ppFld :: post-processed diag or not (=0): =1 stored in qtmp1 ; =2 in qtmp2
57    C   isComputed :: previous post-processed diag (still available in qtmp)
58    C     nLevOutp :: number of levels to write in output file
59  C  C
60  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
61  C     qtmp1 :: thread-shared temporary array (needs to be in common block):  C     qtmp1 :: temporary array; used to store a copy of diag. output field.
62  C              to write a diagnostic field to file, copy it first from (big)  C     qtmp2 :: temporary array; used to store a copy of a 2nd diag. field.
63  C              diagnostic storage qdiag into it.  C-  Note: local common block no longer needed.
64        COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1  c     COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
65        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
66          _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
67    
68        INTEGER i, j, k, lm        INTEGER i, j, k, lm
69        INTEGER bi, bj        INTEGER bi, bj
70        INTEGER md, ndId, ip, im        INTEGER md, ndId, nn, ip, im
71        INTEGER mate, mVec        INTEGER mate, mDbl, mVec
72        CHARACTER*8 parms1        INTEGER ppFld, isComputed
73        _RL undef, getcon        CHARACTER*10 gcode
74        _RL tmpLev        _RL undefRL
75        EXTERNAL getcon        INTEGER nLevOutp, kLev
       INTEGER ILNBLNK  
       EXTERNAL ILNBLNK  
       INTEGER ilen  
76    
77          INTEGER iLen
78        INTEGER ioUnit        INTEGER ioUnit
79        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
80        CHARACTER*(MAX_LEN_MBUF) suff        CHARACTER*(MAX_LEN_MBUF) suff
81        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
82        INTEGER prec, nRec        INTEGER prec, nRec, nTimRec
83          _RL     timeRec(2)
84          _RL     tmpLoc
85  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
86        LOGICAL glf        LOGICAL glf
87  #endif  #endif
88  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
       INTEGER ii  
89        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(NrMax)  
 #endif  
90  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
91    
92  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
93    
94    C---  set file properties
95        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
96        undef = getcon('UNDEF')        undefRL = misValFlt(listId)
97    
98        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
99        ilen = ILNBLNK(fnames(listId))        iLen = ILNBLNK(fnames(listId))
100        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
101    C-    for now, if integrate vertically, output field has just 1 level:
102          nLevOutp = nlevels(listId)
103          IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
104    
105    C--   Set time information:
106          IF ( freq(listId).LT.0. ) THEN
107    C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
108            nTimRec = 1
109            timeRec(1) = myTime
110          ELSE
111    C-    Time-average: store the 2 edges of the time-averaging interval.
112    C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
113    C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
114            nTimRec = 2
115    
116    C-    end of time-averaging interval:
117            timeRec(2) = myTime
118    
119    C-    begining of time-averaging interval:
120    c       timeRec(1) = myTime - freq(listId)
121    C     a) find the time of the previous multiple of output freq:
122            timeRec(1) = myTime-deltaTClock*0.5 _d 0
123            timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
124            i = INT( timeRec(1) )
125            IF ( timeRec(1).LT.0. ) THEN
126              tmpLoc = FLOAT(i)
127              IF ( timeRec(1).NE.tmpLoc ) i = i - 1
128            ENDIF
129            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
130    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
131            timeRec(1) = MAX( timeRec(1), startTime )
132    
133    C     b) round off to nearest multiple of time-step:
134            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
135            i = NINT( timeRec(1) )
136    C     if just half way, NINT will return the next time-step: correct this
137            tmpLoc = FLOAT(i) - 0.5 _d 0
138            IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
139            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
140    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
141          ENDIF
142    C--   Convert time to iteration number (debug)
143    c     DO i=1,nTimRec
144    c       timeRec(i) = timeRec(i)/deltaTClock
145    c     ENDDO
146    
147  #ifdef ALLOW_MNC  C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
148        IF (useMNC .AND. diag_mnc) THEN        DO lm=1,averageCycle(listId)
         DO i = 1,MAX_LEN_FNAM  
           diag_mnc_bn(i:i) = ' '  
         ENDDO  
         DO i = 1,NLEN  
           dn_blnk(i:i) = ' '  
         ENDDO  
         WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)  
   
 C       Update the record dimension by writing the iteration number  
         CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)  
         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  */  
149    
150        ENDIF  #ifdef ALLOW_MNC
151           IF (useMNC .AND. diag_mnc) THEN
152             CALL DIAGNOSTICS_MNC_SET(
153         I                    nLevOutp, listId, lm,
154         O                    diag_mnc_bn,
155         I                    undefRL, myTime, myIter, myThid )
156           ENDIF
157  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
158    
159  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
160    
161        DO md = 1,nfields(listId)         isComputed = 0
162           DO md = 1,nfields(listId)
163          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
164          parms1 = gdiag(ndId)(1:8)          gcode = gdiag(ndId)(1:10)
165          mate = 0          mate = 0
166          mVec = 0          mVec = 0
167          IF ( parms1(5:5).EQ.'C' ) THEN          mDbl = 0
168            ppFld = 0
169            IF ( gcode(5:5).EQ.'C' ) THEN
170  C-      Check for Mate of a Counter Diagnostic  C-      Check for Mate of a Counter Diagnostic
171             READ(parms1,'(5X,I3)') mate             mate = hdiag(ndId)
172          ELSEIF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN          ELSEIF ( gcode(5:5).EQ.'P' ) THEN
173               ppFld = 1
174               IF ( gdiag(hdiag(ndId))(5:5).EQ.'P' ) ppFld = 2
175    C-      Also load the mate (if stored) for Post-Processing
176               nn = ndId
177               DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
178                 nn = hdiag(nn)
179               ENDDO
180               IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
181    c          write(0,*) ppFld,' ndId=', ndId, nn, mDbl, isComputed
182            ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
183  C-      Check for Mate of a Vector Diagnostic  C-      Check for Mate of a Vector Diagnostic
184             READ(parms1,'(5X,I3)') mVec             mVec = hdiag(ndId)
185          ENDIF          ENDIF
186          IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN          IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
187  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
          DO lm=1,averageCycle(listId)  
188    
189            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
190            im = mdiag(md,listId)            im = mdiag(md,listId)
191            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
192              IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
193            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
194    
195            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ppFld.EQ.2 .AND. isComputed.EQ.hdiag(ndId) ) THEN
196    C-        Post-Processed diag from an other Post-Processed diag -and-
197    C         both of them have just been calculated and are still stored in qtmp:
198    C         => skip computation and just write qtmp2
199                IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
200                   WRITE(ioUnit,'(A,I6,3A,I6)')
201         &         '  get Post-Proc. Diag # ', ndId, '  ', cdiag(ndId),
202         &         ' from previous computation of Diag # ', isComputed
203                ENDIF
204                isComputed = 0
205              ELSEIF ( ndiag(ip,1,1).EQ.0 ) THEN
206  C-        Empty diagnostics case :  C-        Empty diagnostics case :
207                isComputed = 0
208    
209              _BEGIN_MASTER( myThid )              _BEGIN_MASTER( myThid )
210              WRITE(msgBuf,'(A,I10)')              WRITE(msgBuf,'(A,I10)')
211       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
212              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
213       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
214              WRITE(msgBuf,'(A,I4,3A,I3,2A)')              WRITE(msgBuf,'(A,I6,3A,I4,2A)')
215       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
216       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
217              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
218       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
219              IF ( averageCycle(listId).GT.1 ) THEN              IF ( averageCycle(listId).GT.1 ) THEN
220               WRITE(msgBuf,'(A,2(I2,A))')               WRITE(msgBuf,'(A,2(I3,A))')
221       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
222       &                                            ndiag(ip,1,1), ' )'       &                                            ndiag(ip,1,1), ' )'
223              ELSE              ELSE
224               WRITE(msgBuf,'(A,2(I2,A))')               WRITE(msgBuf,'(A,2(I3,A))')
225       &        '- WARNING -   has not been filled (ndiag=',       &        '- WARNING -   has not been filled (ndiag=',
226       &                                            ndiag(ip,1,1), ' )'       &                                            ndiag(ip,1,1), ' )'
227              ENDIF              ENDIF
# Line 252  C-        Empty diagnostics case : Line 234  C-        Empty diagnostics case :
234              _END_MASTER( myThid )              _END_MASTER( myThid )
235              DO bj = myByLo(myThid), myByHi(myThid)              DO bj = myByLo(myThid), myByHi(myThid)
236                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
237                  DO k = 1,nlevels(listId)                  DO k = 1,nLevOutp
238                    DO j = 1-OLy,sNy+OLy                    DO j = 1-OLy,sNy+OLy
239                      DO i = 1-OLx,sNx+OLx                      DO i = 1-OLx,sNx+OLx
240                        qtmp1(i,j,k,bi,bj) = 0. _d 0                        qtmp1(i,j,k,bi,bj) = 0. _d 0
# Line 264  C-        Empty diagnostics case : Line 246  C-        Empty diagnostics case :
246    
247            ELSE            ELSE
248  C-        diagnostics is not empty :  C-        diagnostics is not empty :
249                isComputed = 0
250    
251              IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
252                WRITE(ioUnit,'(A,I3,3A,I8,2A)')                IF ( ppFld.GE.1 ) THEN
253                   WRITE(ioUnit,'(A,I6,7A,I8,2A)')
254         &         ' Post-Processing Diag # ', ndId, '  ', cdiag(ndId),
255         &         '   Parms: ',gdiag(ndId)
256                   IF ( mDbl.EQ.0 ) THEN
257                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
258         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
259                   ELSE
260                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
261         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
262         &          ' and diag: ',
263         &            cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
264                   ENDIF
265                  ELSE
266                   WRITE(ioUnit,'(A,I6,3A,I8,2A)')
267       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
268       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
269                  ENDIF
270                IF ( mate.GT.0 ) THEN                IF ( mate.GT.0 ) THEN
271                 WRITE(ioUnit,'(3A,I3,2A)')                 WRITE(ioUnit,'(3A,I6,2A)')
272       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
273       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
274                ELSEIF ( mVec.GT.0 ) THEN                ELSEIF ( mVec.GT.0 ) THEN
275                  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
276                   WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
277       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
278       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
279       &             ' exists '       &             ' exists '
280                  ELSE                  ELSE
281                   WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
282       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
283       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
284       &             ' not enabled'       &             ' not enabled'
# Line 288  C-        diagnostics is not empty : Line 286  C-        diagnostics is not empty :
286                ENDIF                ENDIF
287              ENDIF              ENDIF
288    
289              IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
290  C-       get all the levels (for vertical interpolation)  C-       get only selected levels:
291                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
292                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
293                  DO k = 1,kdiag(ndId)                  DO k = 1,nlevels(listId)
294                    tmpLev = k                    kLev = NINT(levs(k,listId))
295                    CALL GETDIAG(                    CALL DIAGNOSTICS_GET_DIAG(
296       I                         tmpLev,undef,       I                         kLev, undefRL,
297       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
298       I                         ndId,mate,ip,im,bi,bj,myThid)       I                         ndId, mate, ip, im, bi, bj, myThid )
299                  ENDDO                  ENDDO
300                 ENDDO                 ENDDO
301                ENDDO                ENDDO
302                  IF ( mDbl.GT.0 ) THEN
303                   DO bj = myByLo(myThid), myByHi(myThid)
304                    DO bi = myBxLo(myThid), myBxHi(myThid)
305                     DO k = 1,nlevels(listId)
306                      kLev = NINT(levs(k,listId))
307                      CALL DIAGNOSTICS_GET_DIAG(
308         I                         kLev, undefRL,
309         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
310         I                         mDbl, 0, im, 0, bi, bj, myThid )
311                     ENDDO
312                    ENDDO
313                   ENDDO
314                  ENDIF
315              ELSE              ELSE
316  C-       get only selected levels:  C-       get all the levels (for vertical post-processing)
317                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
318                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
319                  DO k = 1,nlevels(listId)                    CALL DIAGNOSTICS_GET_DIAG(
320                    CALL GETDIAG(       I                         0, undefRL,
321       I                         levs(k,listId),undef,       O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
322       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),       I                         ndId, mate, ip, im, bi, bj, myThid )
      I                         ndId,mate,ip,im,bi,bj,myThid)  
                 ENDDO  
323                 ENDDO                 ENDDO
324                ENDDO                ENDDO
325                  IF ( mDbl.GT.0 ) THEN
326                   DO bj = myByLo(myThid), myByHi(myThid)
327                    DO bi = myBxLo(myThid), myBxHi(myThid)
328                      CALL DIAGNOSTICS_GET_DIAG(
329         I                         0, undefRL,
330         O                         qtmp2(1-OLx,1-OLy,1,bi,bj),
331         I                         mDbl, 0, im, 0, bi, bj, myThid )
332                    ENDDO
333                   ENDDO
334                  ENDIF
335              ENDIF              ENDIF
336    
 C-        end of empty diag / not empty block  
           ENDIF  
   
337  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
338  C         Check to see if we need to interpolate before output  C--     Apply specific post-processing (e.g., interpolate) before output
339  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
340            IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
341  C-        Do vertical interpolation:  C-          Do vertical interpolation:
342             IF ( fluidIsAir ) THEN               IF ( fluidIsAir ) THEN
343  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
344              CALL DIAGNOSTICS_INTERP_VERT(                CALL DIAGNOSTICS_INTERP_VERT(
345       I                     listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
346       U                     qtmp1,       U                         qtmp1, qtmp2,
347       I                     undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
348             ELSE               ELSE
349               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
350       &         'INTERP_VERT not allowed in this config'       &           'INTERP_VERT not allowed in this config'
351               CALL PRINT_ERROR( msgBuf , myThid )                 CALL PRINT_ERROR( msgBuf , myThid )
352               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
353             ENDIF               ENDIF
354                ENDIF
355                IF ( fflags(listId)(2:2).EQ.'I' ) THEN
356    C-          Integrate vertically: for now, output field has just 1 level:
357                  CALL DIAGNOSTICS_SUM_LEVELS(
358         I                         listId, md, ndId, ip, im, lm,
359         U                         qtmp1,
360         I                         undefRL, myTime, myIter, myThid )
361                ENDIF
362                IF ( ppFld.GE.1 ) THEN
363    C-          Do Post-Processing:
364                 IF ( flds(md,listId).EQ.'PhiVEL  '
365         &       .OR. flds(md,listId).EQ.'PsiVEL  '
366         &          ) THEN
367                  CALL DIAGNOSTICS_CALC_PHIVEL(
368         I                         listId, md, ndId, ip, im, lm,
369         I                         NrMax,
370         U                         qtmp1, qtmp2,
371         I                         myTime, myIter, myThid )
372                  isComputed = ndId
373                 ELSE
374                   WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
375         &           'unknown Processing method for diag="',cdiag(ndId),'"'
376                   CALL PRINT_ERROR( msgBuf , myThid )
377                   STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
378                 ENDIF
379                ENDIF
380    
381    C--     End of empty diag / not-empty diag block
382            ENDIF            ENDIF
383    
384  C--    Ready to write field "md", element "lm" in averageCycle(listId)  C--     Ready to write field "md", element "lm" in averageCycle(listId)
385    
386  C-        write to binary file, using MDSIO pkg:  C-        write to binary file, using MDSIO pkg:
387            IF ( diag_mdsio ) THEN            IF ( diag_mdsio ) THEN
388              nRec = lm + (md-1)*averageCycle(listId)  c          nRec = lm + (md-1)*averageCycle(listId)
389  C           default precision for output files             nRec = md + (lm-1)*nfields(listId)
390              prec = writeBinaryPrec  C         default precision for output files
391  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision             prec = writeBinaryPrec
392              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32  C         fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
393              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64             IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
394               IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
395  C         a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R  C         a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R
396               IF ( ppFld.LE.1 ) THEN
397              CALL WRITE_REC_LEV_RL(              CALL WRITE_REC_LEV_RL(
398       I                            fn, prec,       I                            fn, prec,
399       I                            NrMax, 1, nlevels(listId),       I                            NrMax, 1, nLevOutp,
400       I                            qtmp1, -nRec, myIter, myThid )       I                            qtmp1, -nRec, myIter, myThid )
401               ELSE
402                CALL WRITE_REC_LEV_RL(
403         I                            fn, prec,
404         I                            NrMax, 1, nLevOutp,
405         I                            qtmp2, -nRec, myIter, myThid )
406               ENDIF
407            ENDIF            ENDIF
408    
409  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
410            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
411               IF ( ppFld.LE.1 ) THEN
412              _BEGIN_MASTER( myThid )              CALL DIAGNOSTICS_MNC_OUT(
413         I                       NrMax, nLevOutp, listId, ndId, mate,
414              DO ii = 1,CW_DIMS       I                       diag_mnc_bn, qtmp1,
415                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)       I                       undefRL, myTime, myIter, myThid )
416                dn(ii)(1:NLEN) = dn_blnk(1:NLEN)             ELSE
417              ENDDO              CALL DIAGNOSTICS_MNC_OUT(
418         I                       NrMax, nLevOutp, listId, ndId, mate,
419  C           Note that the "d_cw_name" variable is a hack that hides a       I                       diag_mnc_bn, qtmp2,
420  C           subtlety within MNC.  Basically, each MNC-wrapped file is       I                       undefRL, myTime, myIter, myThid )
421  C           caching its own concept of what each "grid name" (that is, a             ENDIF
 C           dimension group name) means.  So one cannot re-use the same  
 C           "grid" name for different collections of dimensions within a  
 C           given file.  By appending the "ndId" values to each name, we  
 C           guarantee uniqueness within each MNC-produced file.  
             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', nlevels(listId)  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevels(listId)  
             ENDIF  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zld', nlevels(listId)  
             ENDIF  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)  
             ENDIF  
             dim(3) = NrMax  
             ib(3)  = 1  
             ie(3)  = nlevels(listId)  
   
 C           Time dimension  
             dn(4)(1:1) = 'T'  
             dim(4) = -1  
             ib(4)  = 1  
             ie(4)  = 1  
   
             CALL MNC_CW_ADD_GNAME(d_cw_name, 4,  
      &             dim, dn, ib, ie, myThid)  
             CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,  
      &             4,5, myThid)  
             CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',  
      &             tdiag(ndId),myThid)  
             CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',  
      &             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)  
             ENDIF  
   
             CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)  
             CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)  
   
             _END_MASTER( myThid )  
   
422            ENDIF            ENDIF
423  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
424    
          ENDDO  
425  C--     end of Processing Fld # md  C--     end of Processing Fld # md
426          ENDIF          ENDIF
427           ENDDO
428    
429    C--   end loop on lm counter (= averagePeriod)
430        ENDDO        ENDDO
431    
432  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
433        IF (diag_mdsio) THEN        IF (diag_mdsio) THEN
434  C-    Note: temporary: since it's a pain to add more arguments to  C-    Note: temporary: since it is a pain to add more arguments to
435  C     all MDSIO S/R, uses instead this specific S/R to write only  C     all MDSIO S/R, uses instead this specific S/R to write only
436  C     meta files but with more informations in it.  C     meta files but with more informations in it.
437              glf = globalFiles              glf = globalFiles
438              nRec = nfields(listId)*averageCycle(listId)              nRec = averageCycle(listId)*nfields(listId)
439              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
440       &              0, 0, nlevels(listId), ' ',       &              0, 0, nLevOutp, ' ',
441       &              nfields(listId), flds(1,listId), 1, myTime,       &              nfields(listId), flds(1,listId),
442         &              nTimRec, timeRec, undefRL,
443       &              nRec, myIter, myThid)       &              nRec, myIter, myThid)
444        ENDIF        ENDIF
445  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
Legend:
Removed from v.1.34  
changed lines
  Added in v.1.61
  ViewVC Help
Powered by ViewVC 1.1.22