/[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.55 by jmc, Wed Jun 22 19:09:40 2011 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    #ifdef ALLOW_FIZHI
45          _RL   getcon
46          EXTERNAL getcon
47    #endif
48    
49  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
50  C     i,j,k :: loop indices  C     i,j,k :: loop indices
51    C     bi,bj :: tile indices
52  C     lm    :: loop index (averageCycle)  C     lm    :: loop index (averageCycle)
53  C     md    :: field number in the list "listId".  C     md    :: field number in the list "listId".
54  C     ndId  :: diagnostics  Id number (in available diagnostics list)  C     ndId  :: diagnostics  Id number (in available diagnostics list)
55  C     mate  :: counter mate Id number (in available diagnostics list)  C     mate  :: counter mate Id number (in available diagnostics list)
56  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
57  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
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        CHARACTER*10 gcode
73        _RL undef, getcon        _RL undefRL
74        _RL tmpLev        INTEGER nLevOutp, kLev
       EXTERNAL getcon  
       INTEGER ILNBLNK  
       EXTERNAL ILNBLNK  
       INTEGER ilen  
75    
76          INTEGER iLen
77        INTEGER ioUnit        INTEGER ioUnit
78        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
79        CHARACTER*(MAX_LEN_MBUF) suff        CHARACTER*(MAX_LEN_MBUF) suff
80        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
81        INTEGER prec, nRec        INTEGER prec, nRec, nTimRec
82          _RL     timeRec(2)
83          _RL     tmpLoc
84  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
85        LOGICAL glf        LOGICAL glf
86  #endif  #endif
87  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
       INTEGER ii  
88        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
89        INTEGER CW_DIMS, NLEN        LOGICAL useMissingValue
90        PARAMETER ( CW_DIMS = 10 )        REAL*8 misValLoc
       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  
91  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
92    
93  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
94    
95    C---  set file properties
96        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
97        undef = getcon('UNDEF')        undefRL = UNSET_RL
98    #ifdef ALLOW_FIZHI
99          IF ( useFIZHI ) undefRL = getcon('UNDEF')
100    #endif
101        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
102        ilen = ILNBLNK(fnames(listId))        iLen = ILNBLNK(fnames(listId))
103        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
104    C-    for now, if integrate vertically, output field has just 1 level:
105          nLevOutp = nlevels(listId)
106          IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
107    
108    C--   Set time information:
109          IF ( freq(listId).LT.0. ) THEN
110    C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
111            nTimRec = 1
112            timeRec(1) = myTime
113          ELSE
114    C-    Time-average: store the 2 edges of the time-averaging interval.
115    C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
116    C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
117            nTimRec = 2
118    
119    C-    end of time-averaging interval:
120            timeRec(2) = myTime
121    
122    C-    begining of time-averaging interval:
123    c       timeRec(1) = myTime - freq(listId)
124    C     a) find the time of the previous multiple of output freq:
125            timeRec(1) = myTime-deltaTClock*0.5 _d 0
126            timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
127            i = INT( timeRec(1) )
128            IF ( timeRec(1).LT.0. ) THEN
129              tmpLoc = FLOAT(i)
130              IF ( timeRec(1).NE.tmpLoc ) i = i - 1
131            ENDIF
132            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
133    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
134            timeRec(1) = MAX( timeRec(1), startTime )
135    
136    C     b) round off to nearest multiple of time-step:
137            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
138            i = NINT( timeRec(1) )
139    C     if just half way, NINT will return the next time-step: correct this
140            tmpLoc = FLOAT(i) - 0.5 _d 0
141            IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
142            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
143    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
144          ENDIF
145    C--   Convert time to iteration number (debug)
146    c     DO i=1,nTimRec
147    c       timeRec(i) = timeRec(i)/deltaTClock
148    c     ENDDO
149    
150  #ifdef ALLOW_MNC  C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
151        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  */  
152    
153        ENDIF  #ifdef ALLOW_MNC
154           IF (useMNC .AND. diag_mnc) THEN
155             misValLoc = undefRL
156             IF ( misvalFlt(listId).NE.UNSET_RL )
157         &        misValLoc = misvalFlt(listId)
158             CALL DIAGNOSTICS_MNC_SET(
159         I                    nLevOutp, listId, lm,
160         O                    diag_mnc_bn, useMissingValue,
161         I                    misValLoc, myTime, myIter, myThid )
162           ENDIF
163  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
164    
165  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
166    
167        DO md = 1,nfields(listId)         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            IF ( gcode(5:5).EQ.'C' ) THEN
174  C-      Check for Mate of a Counter Diagnostic  C-      Check for Mate of a Counter Diagnostic
175             READ(parms1,'(5X,I3)') mate             mate = hdiag(ndId)
176          ELSEIF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN          ELSEIF ( gcode(5:5).EQ.'P' ) THEN
177    C-      Also load the mate (if stored) for Post-Processing
178               nn = ndId
179               DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
180                 nn = hdiag(nn)
181               ENDDO
182               IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
183            ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
184  C-      Check for Mate of a Vector Diagnostic  C-      Check for Mate of a Vector Diagnostic
185             READ(parms1,'(5X,I3)') mVec             mVec = hdiag(ndId)
186          ENDIF          ENDIF
187          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
188  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
          DO lm=1,averageCycle(listId)  
189    
190            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
191            im = mdiag(md,listId)            im = mdiag(md,listId)
192            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
193              IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
194            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
195    
196            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ndiag(ip,1,1).EQ.0 ) THEN
# Line 229  C-        Empty diagnostics case : Line 201  C-        Empty diagnostics case :
201       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
202              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
203       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
204              WRITE(msgBuf,'(A,I4,3A,I3,2A)')              WRITE(msgBuf,'(A,I6,3A,I4,2A)')
205       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
206       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
207              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
208       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
209              IF ( averageCycle(listId).GT.1 ) THEN              IF ( averageCycle(listId).GT.1 ) THEN
210               WRITE(msgBuf,'(A,2(I2,A))')               WRITE(msgBuf,'(A,2(I3,A))')
211       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
212       &                                            ndiag(ip,1,1), ' )'       &                                            ndiag(ip,1,1), ' )'
213              ELSE              ELSE
214               WRITE(msgBuf,'(A,2(I2,A))')               WRITE(msgBuf,'(A,2(I3,A))')
215       &        '- WARNING -   has not been filled (ndiag=',       &        '- WARNING -   has not been filled (ndiag=',
216       &                                            ndiag(ip,1,1), ' )'       &                                            ndiag(ip,1,1), ' )'
217              ENDIF              ENDIF
# Line 252  C-        Empty diagnostics case : Line 224  C-        Empty diagnostics case :
224              _END_MASTER( myThid )              _END_MASTER( myThid )
225              DO bj = myByLo(myThid), myByHi(myThid)              DO bj = myByLo(myThid), myByHi(myThid)
226                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
227                  DO k = 1,nlevels(listId)                  DO k = 1,nLevOutp
228                    DO j = 1-OLy,sNy+OLy                    DO j = 1-OLy,sNy+OLy
229                      DO i = 1-OLx,sNx+OLx                      DO i = 1-OLx,sNx+OLx
230                        qtmp1(i,j,k,bi,bj) = 0. _d 0                        qtmp1(i,j,k,bi,bj) = 0. _d 0
# Line 265  C-        Empty diagnostics case : Line 237  C-        Empty diagnostics case :
237            ELSE            ELSE
238  C-        diagnostics is not empty :  C-        diagnostics is not empty :
239    
240              IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
241                WRITE(ioUnit,'(A,I3,3A,I8,2A)')                IF ( gcode(5:5).EQ.'P' ) THEN
242                   WRITE(ioUnit,'(A,I6,7A,I8,2A)')
243         &         ' Post-Processing Diag # ', ndId, '  ', cdiag(ndId),
244         &         '   Parms: ',gdiag(ndId)
245                   IF ( mDbl.EQ.0 ) THEN
246                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
247         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
248                   ELSE
249                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
250         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
251         &          ' and diag: ',
252         &            cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
253                   ENDIF
254                  ELSE
255                   WRITE(ioUnit,'(A,I6,3A,I8,2A)')
256       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
257       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
258                  ENDIF
259                IF ( mate.GT.0 ) THEN                IF ( mate.GT.0 ) THEN
260                 WRITE(ioUnit,'(3A,I3,2A)')                 WRITE(ioUnit,'(3A,I6,2A)')
261       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
262       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
263                ELSEIF ( mVec.GT.0 ) THEN                ELSEIF ( mVec.GT.0 ) THEN
264                  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
265                   WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
266       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
267       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
268       &             ' exists '       &             ' exists '
269                  ELSE                  ELSE
270                   WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
271       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
272       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
273       &             ' not enabled'       &             ' not enabled'
# Line 288  C-        diagnostics is not empty : Line 275  C-        diagnostics is not empty :
275                ENDIF                ENDIF
276              ENDIF              ENDIF
277    
278              IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
279  C-       get all the levels (for vertical interpolation)  C-       get only selected levels:
280                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
281                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
282                  DO k = 1,kdiag(ndId)                  DO k = 1,nlevels(listId)
283                    tmpLev = k                    kLev = NINT(levs(k,listId))
284                    CALL GETDIAG(                    CALL DIAGNOSTICS_GET_DIAG(
285       I                         tmpLev,undef,       I                         kLev, undefRL,
286       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
287       I                         ndId,mate,ip,im,bi,bj,myThid)       I                         ndId, mate, ip, im, bi, bj, myThid )
288                  ENDDO                  ENDDO
289                 ENDDO                 ENDDO
290                ENDDO                ENDDO
291                  IF ( mDbl.GT.0 ) THEN
292                   DO bj = myByLo(myThid), myByHi(myThid)
293                    DO bi = myBxLo(myThid), myBxHi(myThid)
294                     DO k = 1,nlevels(listId)
295                      kLev = NINT(levs(k,listId))
296                      CALL DIAGNOSTICS_GET_DIAG(
297         I                         kLev, undefRL,
298         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
299         I                         mDbl, 0, im, 0, bi, bj, myThid )
300                     ENDDO
301                    ENDDO
302                   ENDDO
303                  ENDIF
304              ELSE              ELSE
305  C-       get only selected levels:  C-       get all the levels (for vertical post-processing)
306                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
307                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
308                  DO k = 1,nlevels(listId)                    CALL DIAGNOSTICS_GET_DIAG(
309                    CALL GETDIAG(       I                         0, undefRL,
310       I                         levs(k,listId),undef,       O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
311       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  
312                 ENDDO                 ENDDO
313                ENDDO                ENDDO
314                  IF ( mDbl.GT.0 ) THEN
315                   DO bj = myByLo(myThid), myByHi(myThid)
316                    DO bi = myBxLo(myThid), myBxHi(myThid)
317                     DO k = 1,nlevels(listId)
318                      CALL DIAGNOSTICS_GET_DIAG(
319         I                         0, undefRL,
320         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
321         I                         mDbl, 0, im, 0, bi, bj, myThid )
322                     ENDDO
323                    ENDDO
324                   ENDDO
325                  ENDIF
326              ENDIF              ENDIF
327    
 C-        end of empty diag / not empty block  
           ENDIF  
   
328  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
329  C         Check to see if we need to interpolate before output  C--     Apply specific post-processing (e.g., interpolate) before output
330  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
331            IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
332  C-        Do vertical interpolation:  C-          Do vertical interpolation:
333             IF ( fluidIsAir ) THEN               IF ( fluidIsAir ) THEN
334  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);
335              CALL DIAGNOSTICS_INTERP_VERT(                CALL DIAGNOSTICS_INTERP_VERT(
336       I                     listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
337       U                     qtmp1,       U                         qtmp1, qtmp2,
338       I                     undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
339             ELSE               ELSE
340               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
341       &         'INTERP_VERT not allowed in this config'       &           'INTERP_VERT not allowed in this config'
342               CALL PRINT_ERROR( msgBuf , myThid )                 CALL PRINT_ERROR( msgBuf , myThid )
343               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
344             ENDIF               ENDIF
345                ENDIF
346                IF ( fflags(listId)(2:2).EQ.'I' ) THEN
347    C-          Integrate vertically: for now, output field has just 1 level:
348                  CALL DIAGNOSTICS_SUM_LEVELS(
349         I                         listId, md, ndId, ip, im, lm,
350         U                         qtmp1,
351         I                         undefRL, myTime, myIter, myThid )
352                ENDIF
353                IF ( gcode(5:5).EQ.'P' ) THEN
354    C-          Do Post-Processing:
355                 IF ( flds(md,listId).EQ.'PhiVEL  '
356    c    &       .OR. flds(md,listId).EQ.'PsiVEL  '
357         &          ) THEN
358                  CALL DIAGNOSTICS_CALC_PHIVEL(
359         I                         listId, md, ndId, ip, im, lm,
360         U                         qtmp1, qtmp2,
361         I                         myTime, myIter, myThid )
362                 ELSE
363                   WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
364         &           'unknown Processing method for diag="',cdiag(ndId),'"'
365                   CALL PRINT_ERROR( msgBuf , myThid )
366                   STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
367                 ENDIF
368                ENDIF
369    
370    C--     End of empty diag / not-empty diag block
371            ENDIF            ENDIF
372    
373  C--    Ready to write field "md", element "lm" in averageCycle(listId)  C--     Ready to write field "md", element "lm" in averageCycle(listId)
374    
375  C-        write to binary file, using MDSIO pkg:  C-        write to binary file, using MDSIO pkg:
376            IF ( diag_mdsio ) THEN            IF ( diag_mdsio ) THEN
377              nRec = lm + (md-1)*averageCycle(listId)  c           nRec = lm + (md-1)*averageCycle(listId)
378                nRec = md + (lm-1)*nfields(listId)
379  C           default precision for output files  C           default precision for output files
380              prec = writeBinaryPrec              prec = writeBinaryPrec
381  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
# Line 350  C           fFlag(1)=R(or D): force it t Line 384  C           fFlag(1)=R(or D): force it t
384  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
385              CALL WRITE_REC_LEV_RL(              CALL WRITE_REC_LEV_RL(
386       I                            fn, prec,       I                            fn, prec,
387       I                            NrMax, 1, nlevels(listId),       I                            NrMax, 1, nLevOutp,
388       I                            qtmp1, -nRec, myIter, myThid )       I                            qtmp1, -nRec, myIter, myThid )
389            ENDIF            ENDIF
390    
391  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
392            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
393                CALL DIAGNOSTICS_MNC_OUT(
394              _BEGIN_MASTER( myThid )       I                       NrMax, nLevOutp, listId, ndId,
395         I                       diag_mnc_bn,
396              DO ii = 1,CW_DIMS       I                       useMissingValue, misValLoc,
397                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)       I                       qtmp1,
398                dn(ii)(1:NLEN) = dn_blnk(1:NLEN)       I                       myTime, myIter, myThid )
             ENDDO  
   
 C           Note that the "d_cw_name" variable is a hack that hides a  
 C           subtlety within MNC.  Basically, each MNC-wrapped file is  
 C           caching its own concept of what each "grid name" (that is, a  
 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 )  
   
399            ENDIF            ENDIF
400  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
401    
          ENDDO  
402  C--     end of Processing Fld # md  C--     end of Processing Fld # md
403          ENDIF          ENDIF
404           ENDDO
405    
406    C--   end loop on lm counter (= averagePeriod)
407        ENDDO        ENDDO
408    
409  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
410        IF (diag_mdsio) THEN        IF (diag_mdsio) THEN
411  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
412  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
413  C     meta files but with more informations in it.  C     meta files but with more informations in it.
414              glf = globalFiles              glf = globalFiles
415              nRec = nfields(listId)*averageCycle(listId)              nRec = averageCycle(listId)*nfields(listId)
416              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
417       &              0, 0, nlevels(listId), ' ',       &              0, 0, nLevOutp, ' ',
418       &              nfields(listId), flds(1,listId), 1, myTime,       &              nfields(listId), flds(1,listId), nTimRec, timeRec,
419       &              nRec, myIter, myThid)       &              nRec, myIter, myThid)
420        ENDIF        ENDIF
421  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
Legend:
Removed from v.1.34  
changed lines
  Added in v.1.55
  ViewVC Help
Powered by ViewVC 1.1.22