/[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.62 by gforget, Wed Jan 11 00:22:48 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,
      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,jLen
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          IF ( diag_mdsio.AND.(diagMdsDir.NE.' ') ) THEN
102            jLen = ILNBLNK(diagMdsDir)
103            WRITE( fn, '(5A)' ) diagMdsDir(1:jLen),'/',
104         &    fnames(listId)(1:iLen),'.',suff(1:10)
105          ENDIF
106    C-    for now, if integrate vertically, output field has just 1 level:
107          nLevOutp = nlevels(listId)
108          IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
109    
110    C--   Set time information:
111          IF ( freq(listId).LT.0. ) THEN
112    C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
113            nTimRec = 1
114            timeRec(1) = myTime
115          ELSE
116    C-    Time-average: store the 2 edges of the time-averaging interval.
117    C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
118    C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
119            nTimRec = 2
120    
121    C-    end of time-averaging interval:
122            timeRec(2) = myTime
123    
124    C-    begining of time-averaging interval:
125    c       timeRec(1) = myTime - freq(listId)
126    C     a) find the time of the previous multiple of output freq:
127            timeRec(1) = myTime-deltaTClock*0.5 _d 0
128            timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
129            i = INT( timeRec(1) )
130            IF ( timeRec(1).LT.0. ) THEN
131              tmpLoc = FLOAT(i)
132              IF ( timeRec(1).NE.tmpLoc ) i = i - 1
133            ENDIF
134            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
135    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
136            timeRec(1) = MAX( timeRec(1), startTime )
137    
138    C     b) round off to nearest multiple of time-step:
139            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
140            i = NINT( timeRec(1) )
141    C     if just half way, NINT will return the next time-step: correct this
142            tmpLoc = FLOAT(i) - 0.5 _d 0
143            IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
144            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
145    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
146          ENDIF
147    C--   Convert time to iteration number (debug)
148    c     DO i=1,nTimRec
149    c       timeRec(i) = timeRec(i)/deltaTClock
150    c     ENDDO
151    
152  #ifdef ALLOW_MNC  C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
153        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  */  
154    
155        ENDIF  #ifdef ALLOW_MNC
156           IF (useMNC .AND. diag_mnc) THEN
157             CALL DIAGNOSTICS_MNC_SET(
158         I                    nLevOutp, listId, lm,
159         O                    diag_mnc_bn,
160         I                    undefRL, myTime, myIter, myThid )
161           ENDIF
162  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
163    
164  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
165    
166        DO md = 1,nfields(listId)         isComputed = 0
167           DO md = 1,nfields(listId)
168          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
169          parms1 = gdiag(ndId)(1:8)          gcode = gdiag(ndId)(1:10)
170          mate = 0          mate = 0
171          mVec = 0          mVec = 0
172          IF ( parms1(5:5).EQ.'C' ) THEN          mDbl = 0
173            ppFld = 0
174            IF ( gcode(5:5).EQ.'C' ) THEN
175  C-      Check for Mate of a Counter Diagnostic  C-      Check for Mate of a Counter Diagnostic
176             READ(parms1,'(5X,I3)') mate             mate = hdiag(ndId)
177          ELSEIF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN          ELSEIF ( gcode(5:5).EQ.'P' ) THEN
178               ppFld = 1
179               IF ( gdiag(hdiag(ndId))(5:5).EQ.'P' ) ppFld = 2
180    C-      Also load the mate (if stored) for Post-Processing
181               nn = ndId
182               DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
183                 nn = hdiag(nn)
184               ENDDO
185               IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
186    c          write(0,*) ppFld,' ndId=', ndId, nn, mDbl, isComputed
187            ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
188  C-      Check for Mate of a Vector Diagnostic  C-      Check for Mate of a Vector Diagnostic
189             READ(parms1,'(5X,I3)') mVec             mVec = hdiag(ndId)
190          ENDIF          ENDIF
191          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
192  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
          DO lm=1,averageCycle(listId)  
193    
194            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
195            im = mdiag(md,listId)            im = mdiag(md,listId)
196            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
197              IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
198            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
199    
200            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ppFld.EQ.2 .AND. isComputed.EQ.hdiag(ndId) ) THEN
201    C-        Post-Processed diag from an other Post-Processed diag -and-
202    C         both of them have just been calculated and are still stored in qtmp:
203    C         => skip computation and just write qtmp2
204                IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
205                   WRITE(ioUnit,'(A,I6,3A,I6)')
206         &         '  get Post-Proc. Diag # ', ndId, '  ', cdiag(ndId),
207         &         ' from previous computation of Diag # ', isComputed
208                ENDIF
209                isComputed = 0
210              ELSEIF ( ndiag(ip,1,1).EQ.0 ) THEN
211  C-        Empty diagnostics case :  C-        Empty diagnostics case :
212                isComputed = 0
213    
214              _BEGIN_MASTER( myThid )              _BEGIN_MASTER( myThid )
215              WRITE(msgBuf,'(A,I10)')              WRITE(msgBuf,'(A,I10)')
216       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
217              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
218       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
219              WRITE(msgBuf,'(A,I4,3A,I3,2A)')              WRITE(msgBuf,'(A,I6,3A,I4,2A)')
220       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
221       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
222              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
223       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
224              IF ( averageCycle(listId).GT.1 ) THEN              IF ( averageCycle(listId).GT.1 ) THEN
225               WRITE(msgBuf,'(A,2(I2,A))')               WRITE(msgBuf,'(A,2(I3,A))')
226       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
227       &                                            ndiag(ip,1,1), ' )'       &                                            ndiag(ip,1,1), ' )'
228              ELSE              ELSE
229               WRITE(msgBuf,'(A,2(I2,A))')               WRITE(msgBuf,'(A,2(I3,A))')
230       &        '- WARNING -   has not been filled (ndiag=',       &        '- WARNING -   has not been filled (ndiag=',
231       &                                            ndiag(ip,1,1), ' )'       &                                            ndiag(ip,1,1), ' )'
232              ENDIF              ENDIF
# Line 252  C-        Empty diagnostics case : Line 239  C-        Empty diagnostics case :
239              _END_MASTER( myThid )              _END_MASTER( myThid )
240              DO bj = myByLo(myThid), myByHi(myThid)              DO bj = myByLo(myThid), myByHi(myThid)
241                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
242                  DO k = 1,nlevels(listId)                  DO k = 1,nLevOutp
243                    DO j = 1-OLy,sNy+OLy                    DO j = 1-OLy,sNy+OLy
244                      DO i = 1-OLx,sNx+OLx                      DO i = 1-OLx,sNx+OLx
245                        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 251  C-        Empty diagnostics case :
251    
252            ELSE            ELSE
253  C-        diagnostics is not empty :  C-        diagnostics is not empty :
254                isComputed = 0
255    
256              IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
257                WRITE(ioUnit,'(A,I3,3A,I8,2A)')                IF ( ppFld.GE.1 ) THEN
258                   WRITE(ioUnit,'(A,I6,7A,I8,2A)')
259         &         ' Post-Processing Diag # ', ndId, '  ', cdiag(ndId),
260         &         '   Parms: ',gdiag(ndId)
261                   IF ( mDbl.EQ.0 ) THEN
262                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
263         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
264                   ELSE
265                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
266         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
267         &          ' and diag: ',
268         &            cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
269                   ENDIF
270                  ELSE
271                   WRITE(ioUnit,'(A,I6,3A,I8,2A)')
272       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
273       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
274                  ENDIF
275                IF ( mate.GT.0 ) THEN                IF ( mate.GT.0 ) THEN
276                 WRITE(ioUnit,'(3A,I3,2A)')                 WRITE(ioUnit,'(3A,I6,2A)')
277       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
278       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
279                ELSEIF ( mVec.GT.0 ) THEN                ELSEIF ( mVec.GT.0 ) THEN
280                  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
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       &             ' exists '       &             ' exists '
285                  ELSE                  ELSE
286                   WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
287       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
288       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
289       &             ' not enabled'       &             ' not enabled'
# Line 288  C-        diagnostics is not empty : Line 291  C-        diagnostics is not empty :
291                ENDIF                ENDIF
292              ENDIF              ENDIF
293    
294              IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
295  C-       get all the levels (for vertical interpolation)  C-       get only selected levels:
296                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
297                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
298                  DO k = 1,kdiag(ndId)                  DO k = 1,nlevels(listId)
299                    tmpLev = k                    kLev = NINT(levs(k,listId))
300                    CALL GETDIAG(                    CALL DIAGNOSTICS_GET_DIAG(
301       I                         tmpLev,undef,       I                         kLev, undefRL,
302       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
303       I                         ndId,mate,ip,im,bi,bj,myThid)       I                         ndId, mate, ip, im, bi, bj, myThid )
304                  ENDDO                  ENDDO
305                 ENDDO                 ENDDO
306                ENDDO                ENDDO
307                  IF ( mDbl.GT.0 ) THEN
308                   DO bj = myByLo(myThid), myByHi(myThid)
309                    DO bi = myBxLo(myThid), myBxHi(myThid)
310                     DO k = 1,nlevels(listId)
311                      kLev = NINT(levs(k,listId))
312                      CALL DIAGNOSTICS_GET_DIAG(
313         I                         kLev, undefRL,
314         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
315         I                         mDbl, 0, im, 0, bi, bj, myThid )
316                     ENDDO
317                    ENDDO
318                   ENDDO
319                  ENDIF
320              ELSE              ELSE
321  C-       get only selected levels:  C-       get all the levels (for vertical post-processing)
322                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
323                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
324                  DO k = 1,nlevels(listId)                    CALL DIAGNOSTICS_GET_DIAG(
325                    CALL GETDIAG(       I                         0, undefRL,
326       I                         levs(k,listId),undef,       O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
327       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  
328                 ENDDO                 ENDDO
329                ENDDO                ENDDO
330                  IF ( mDbl.GT.0 ) THEN
331                   DO bj = myByLo(myThid), myByHi(myThid)
332                    DO bi = myBxLo(myThid), myBxHi(myThid)
333                      CALL DIAGNOSTICS_GET_DIAG(
334         I                         0, undefRL,
335         O                         qtmp2(1-OLx,1-OLy,1,bi,bj),
336         I                         mDbl, 0, im, 0, bi, bj, myThid )
337                    ENDDO
338                   ENDDO
339                  ENDIF
340              ENDIF              ENDIF
341    
 C-        end of empty diag / not empty block  
           ENDIF  
   
342  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
343  C         Check to see if we need to interpolate before output  C--     Apply specific post-processing (e.g., interpolate) before output
344  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
345            IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
346  C-        Do vertical interpolation:  C-          Do vertical interpolation:
347             IF ( fluidIsAir ) THEN               IF ( fluidIsAir ) THEN
348  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);
349              CALL DIAGNOSTICS_INTERP_VERT(                CALL DIAGNOSTICS_INTERP_VERT(
350       I                     listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
351       U                     qtmp1,       U                         qtmp1, qtmp2,
352       I                     undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
353             ELSE               ELSE
354               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
355       &         'INTERP_VERT not allowed in this config'       &           'INTERP_VERT not allowed in this config'
356               CALL PRINT_ERROR( msgBuf , myThid )                 CALL PRINT_ERROR( msgBuf , myThid )
357               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
358             ENDIF               ENDIF
359                ENDIF
360                IF ( fflags(listId)(2:2).EQ.'I' ) THEN
361    C-          Integrate vertically: for now, output field has just 1 level:
362                  CALL DIAGNOSTICS_SUM_LEVELS(
363         I                         listId, md, ndId, ip, im, lm,
364         U                         qtmp1,
365         I                         undefRL, myTime, myIter, myThid )
366                ENDIF
367                IF ( ppFld.GE.1 ) THEN
368    C-          Do Post-Processing:
369                 IF ( flds(md,listId).EQ.'PhiVEL  '
370         &       .OR. flds(md,listId).EQ.'PsiVEL  '
371         &          ) THEN
372                  CALL DIAGNOSTICS_CALC_PHIVEL(
373         I                         listId, md, ndId, ip, im, lm,
374         I                         NrMax,
375         U                         qtmp1, qtmp2,
376         I                         myTime, myIter, myThid )
377                  isComputed = ndId
378                 ELSE
379                   WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
380         &           'unknown Processing method for diag="',cdiag(ndId),'"'
381                   CALL PRINT_ERROR( msgBuf , myThid )
382                   STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
383                 ENDIF
384                ENDIF
385    
386    C--     End of empty diag / not-empty diag block
387            ENDIF            ENDIF
388    
389  C--    Ready to write field "md", element "lm" in averageCycle(listId)  C--     Ready to write field "md", element "lm" in averageCycle(listId)
390    
391  C-        write to binary file, using MDSIO pkg:  C-        write to binary file, using MDSIO pkg:
392            IF ( diag_mdsio ) THEN            IF ( diag_mdsio ) THEN
393              nRec = lm + (md-1)*averageCycle(listId)  c          nRec = lm + (md-1)*averageCycle(listId)
394  C           default precision for output files             nRec = md + (lm-1)*nfields(listId)
395              prec = writeBinaryPrec  C         default precision for output files
396  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision             prec = writeBinaryPrec
397              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32  C         fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
398              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64             IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
399               IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
400  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
401               IF ( ppFld.LE.1 ) THEN
402              CALL WRITE_REC_LEV_RL(              CALL WRITE_REC_LEV_RL(
403       I                            fn, prec,       I                            fn, prec,
404       I                            NrMax, 1, nlevels(listId),       I                            NrMax, 1, nLevOutp,
405       I                            qtmp1, -nRec, myIter, myThid )       I                            qtmp1, -nRec, myIter, myThid )
406               ELSE
407                CALL WRITE_REC_LEV_RL(
408         I                            fn, prec,
409         I                            NrMax, 1, nLevOutp,
410         I                            qtmp2, -nRec, myIter, myThid )
411               ENDIF
412            ENDIF            ENDIF
413    
414  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
415            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
416               IF ( ppFld.LE.1 ) THEN
417              _BEGIN_MASTER( myThid )              CALL DIAGNOSTICS_MNC_OUT(
418         I                       NrMax, nLevOutp, listId, ndId, mate,
419              DO ii = 1,CW_DIMS       I                       diag_mnc_bn, qtmp1,
420                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)       I                       undefRL, myTime, myIter, myThid )
421                dn(ii)(1:NLEN) = dn_blnk(1:NLEN)             ELSE
422              ENDDO              CALL DIAGNOSTICS_MNC_OUT(
423         I                       NrMax, nLevOutp, listId, ndId, mate,
424  C           Note that the "d_cw_name" variable is a hack that hides a       I                       diag_mnc_bn, qtmp2,
425  C           subtlety within MNC.  Basically, each MNC-wrapped file is       I                       undefRL, myTime, myIter, myThid )
426  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 )  
   
427            ENDIF            ENDIF
428  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
429    
          ENDDO  
430  C--     end of Processing Fld # md  C--     end of Processing Fld # md
431          ENDIF          ENDIF
432           ENDDO
433    
434    C--   end loop on lm counter (= averagePeriod)
435        ENDDO        ENDDO
436    
437  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
438        IF (diag_mdsio) THEN        IF (diag_mdsio) THEN
439  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
440  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
441  C     meta files but with more informations in it.  C     meta files but with more informations in it.
442              glf = globalFiles              glf = globalFiles
443              nRec = nfields(listId)*averageCycle(listId)              nRec = averageCycle(listId)*nfields(listId)
444              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
445       &              0, 0, nlevels(listId), ' ',       &              0, 0, nLevOutp, ' ',
446       &              nfields(listId), flds(1,listId), 1, myTime,       &              nfields(listId), flds(1,listId),
447         &              nTimRec, timeRec, undefRL,
448       &              nRec, myIter, myThid)       &              nRec, myIter, myThid)
449        ENDIF        ENDIF
450  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
Legend:
Removed from v.1.34  
changed lines
  Added in v.1.62
  ViewVC Help
Powered by ViewVC 1.1.22