/[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.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 38  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)  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    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  C
57  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
58  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.
59  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.
60  C              diagnostic storage qdiag into it.  C-  Note: local common block no longer needed.
61        COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1  c     COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
62        _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)
63          _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
64    
65        INTEGER i, j, k, lm        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        _RL undef, getcon        CHARACTER*10 gcode
71        _RL tmpLev        _RL undefRL
72        EXTERNAL getcon        INTEGER nLevOutp, kLev
       INTEGER ILNBLNK  
       EXTERNAL ILNBLNK  
       INTEGER ilen  
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        INTEGER prec, nRec, nTimRec
80          _RL     timeRec(2)
81          _RL     tmpLoc
82  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
83        LOGICAL glf        LOGICAL glf
84  #endif  #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(NrMax)  
 #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)
       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  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
167    
168        DO md = 1,nfields(listId)         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          mate = 0          mate = 0
173          mVec = 0          mVec = 0
174          IF ( parms1(5:5).EQ.'C' ) THEN          mDbl = 0
175            ppFld = 0
176            IF ( gcode(5:5).EQ.'C' ) THEN
177  C-      Check for Mate of a Counter Diagnostic  C-      Check for Mate of a Counter Diagnostic
178             READ(parms1,'(5X,I3)') mate             mate = hdiag(ndId)
179          ELSEIF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN          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  C-      Check for Mate of a Vector Diagnostic
191             READ(parms1,'(5X,I3)') mVec             mVec = hdiag(ndId)
192          ENDIF          ENDIF
193          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
194  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
          DO lm=1,averageCycle(listId)  
195    
196            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
197            im = mdiag(md,listId)            im = mdiag(md,listId)
198            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)            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)            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
201    
202            IF ( ndiag(ip,1,1).EQ.0 ) THEN            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              IF ( averageCycle(listId).GT.1 ) THEN              IF ( averageCycle(listId).GT.1 ) THEN
227               WRITE(msgBuf,'(A,2(I2,A))')               WRITE(msgBuf,'(A,2(I3,A))')
228       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
229       &                                            ndiag(ip,1,1), ' )'       &                                            ndiag(ip,1,1), ' )'
230              ELSE              ELSE
231               WRITE(msgBuf,'(A,2(I2,A))')               WRITE(msgBuf,'(A,2(I3,A))')
232       &        '- WARNING -   has not been filled (ndiag=',       &        '- WARNING -   has not been filled (ndiag=',
233       &                                            ndiag(ip,1,1), ' )'       &                                            ndiag(ip,1,1), ' )'
234              ENDIF              ENDIF
# Line 252  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 264  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 ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
259                WRITE(ioUnit,'(A,I3,3A,I8,2A)')                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 ( mate.GT.0 ) THEN                IF ( mate.GT.0 ) THEN
278                 WRITE(ioUnit,'(3A,I3,2A)')                 WRITE(ioUnit,'(3A,I6,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                ELSEIF ( mVec.GT.0 ) THEN
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                   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                   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 288  C-        diagnostics is not empty : Line 293  C-        diagnostics is not empty :
293                ENDIF                ENDIF
294              ENDIF              ENDIF
295    
296              IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
297  C-       get all the levels (for vertical interpolation)  C-       get only selected levels:
298                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
299                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
300                  DO k = 1,kdiag(ndId)                  DO k = 1,nlevels(listId)
301                    tmpLev = k                    kLev = NINT(levs(k,listId))
302                    CALL GETDIAG(                    CALL DIAGNOSTICS_GET_DIAG(
303       I                         tmpLev,undef,       I                         kLev, undefRL,
304       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
305       I                         ndId,mate,ip,im,bi,bj,myThid)       I                         ndId, mate, ip, im, bi, bj, myThid )
306                  ENDDO                  ENDDO
307                 ENDDO                 ENDDO
308                ENDDO                ENDDO
309                  IF ( mDbl.GT.0 ) THEN
310                   DO bj = myByLo(myThid), myByHi(myThid)
311                    DO bi = myBxLo(myThid), myBxHi(myThid)
312                     DO k = 1,nlevels(listId)
313                      kLev = NINT(levs(k,listId))
314                      CALL DIAGNOSTICS_GET_DIAG(
315         I                         kLev, undefRL,
316         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
317         I                         mDbl, 0, im, 0, bi, bj, myThid )
318                     ENDDO
319                    ENDDO
320                   ENDDO
321                  ENDIF
322              ELSE              ELSE
323  C-       get only selected levels:  C-       get all the levels (for vertical post-processing)
324                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
325                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
326                  DO k = 1,nlevels(listId)                    CALL DIAGNOSTICS_GET_DIAG(
327                    CALL GETDIAG(       I                         0, undefRL,
328       I                         levs(k,listId),undef,       O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
329       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  
330                 ENDDO                 ENDDO
331                ENDDO                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
343    
 C-        end of empty diag / not empty block  
           ENDIF  
   
344  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
345  C         Check to see if we need to interpolate before output  C--     Apply specific post-processing (e.g., interpolate) before output
346  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
347            IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
348  C-        Do vertical interpolation:  C-          Do vertical interpolation:
349             IF ( fluidIsAir ) THEN               IF ( fluidIsAir ) THEN
350  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);
351              CALL DIAGNOSTICS_INTERP_VERT(                CALL DIAGNOSTICS_INTERP_VERT(
352       I                     listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
353       U                     qtmp1,       U                         qtmp1, qtmp2,
354       I                     undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
355             ELSE               ELSE
356               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
357       &         'INTERP_VERT not allowed in this config'       &           'INTERP_VERT not allowed in this config'
358               CALL PRINT_ERROR( msgBuf , myThid )                 CALL PRINT_ERROR( msgBuf , myThid )
359               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
360             ENDIF               ENDIF
361                ENDIF
362                IF ( fflags(listId)(2:2).EQ.'I' ) THEN
363    C-          Integrate vertically: for now, output field has just 1 level:
364                  CALL DIAGNOSTICS_SUM_LEVELS(
365         I                         listId, md, ndId, ip, im, lm,
366         U                         qtmp1,
367         I                         undefRL, myTime, myIter, myThid )
368                ENDIF
369                IF ( ppFld.GE.1 ) THEN
370    C-          Do Post-Processing:
371                 IF ( flds(md,listId).EQ.'PhiVEL  '
372         &       .OR. flds(md,listId).EQ.'PsiVEL  '
373         &          ) THEN
374                  CALL DIAGNOSTICS_CALC_PHIVEL(
375         I                         listId, md, ndId, ip, im, lm,
376         I                         NrMax,
377         U                         qtmp1, qtmp2,
378         I                         myTime, myIter, myThid )
379                  isComputed = ndId
380                 ELSE
381                   WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
382         &           'unknown Processing method for diag="',cdiag(ndId),'"'
383                   CALL PRINT_ERROR( msgBuf , myThid )
384                   STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
385                 ENDIF
386                ENDIF
387    
388    C--     End of empty diag / not-empty diag block
389            ENDIF            ENDIF
390    
391  C--    Ready to write field "md", element "lm" in averageCycle(listId)  C--     Ready to write field "md", element "lm" in averageCycle(listId)
392    
393  C-        write to binary file, using MDSIO pkg:  C-        write to binary file, using MDSIO pkg:
394            IF ( diag_mdsio ) THEN            IF ( diag_mdsio ) THEN
395              nRec = lm + (md-1)*averageCycle(listId)  c          nRec = lm + (md-1)*averageCycle(listId)
396  C           default precision for output files             nRec = md + (lm-1)*nfields(listId)
397              prec = writeBinaryPrec  C         default precision for output files
398  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision             prec = writeBinaryPrec
399              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32  C         fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
400              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64             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  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(              CALL WRITE_REC_LEV_RL(
405       I                            fn, prec,       I                            fn, prec,
406       I                            NrMax, 1, nlevels(listId),       I                            NrMax, 1, nLevOutp,
407       I                            qtmp1, -nRec, myIter, myThid )       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            ENDIF
415    
416  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
417            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
418               IF ( ppFld.LE.1 ) THEN
419              _BEGIN_MASTER( myThid )              CALL DIAGNOSTICS_MNC_OUT(
420         I                       NrMax, nLevOutp, listId, ndId, mate,
421              DO ii = 1,CW_DIMS       I                       diag_mnc_bn, qtmp1,
422                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)       I                       undefRL, myTime, myIter, myThid )
423                dn(ii)(1:NLEN) = dn_blnk(1:NLEN)             ELSE
424              ENDDO              CALL DIAGNOSTICS_MNC_OUT(
425         I                       NrMax, nLevOutp, listId, ndId, mate,
426  C           Note that the "d_cw_name" variable is a hack that hides a       I                       diag_mnc_bn, qtmp2,
427  C           subtlety within MNC.  Basically, each MNC-wrapped file is       I                       undefRL, myTime, myIter, myThid )
428  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 )  
   
429            ENDIF            ENDIF
430  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
431    
          ENDDO  
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  #ifdef ALLOW_MDSIO
440        IF (diag_mdsio) THEN        IF (diag_mdsio) THEN
441  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
442  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
443  C     meta files but with more informations in it.  C     meta files but with more informations in it.
444              glf = globalFiles              glf = globalFiles
445              nRec = nfields(listId)*averageCycle(listId)              nRec = averageCycle(listId)*nfields(listId)
446              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
447       &              0, 0, nlevels(listId), ' ',       &              0, 0, nLevOutp, ' ',
448       &              nfields(listId), flds(1,listId), 1, myTime,       &              nfields(listId), flds(1,listId),
449         &              nTimRec, timeRec, undefRL,
450       &              nRec, myIter, myThid)       &              nRec, myIter, myThid)
451        ENDIF        ENDIF
452  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
Legend:
Removed from v.1.34  
changed lines
  Added in v.1.63
  ViewVC Help
Powered by ViewVC 1.1.22