/[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.18 by molod, Wed Jul 13 17:47:21 2005 UTC revision 1.60 by jmc, Sun Jan 13 22:46:38 2013 UTC
# Line 8  CBOP 0 Line 8  CBOP 0
8  C     !ROUTINE: DIAGNOSTICS_OUT  C     !ROUTINE: DIAGNOSTICS_OUT
9    
10  C     !INTERFACE:  C     !INTERFACE:
11        SUBROUTINE  DIAGNOSTICS_OUT(        SUBROUTINE DIAGNOSTICS_OUT(
12       I     listId,       I     listId,
      I     myIter,  
13       I     myTime,       I     myTime,
14         I     myIter,
15       I     myThid )       I     myThid )
16    
17  C     !DESCRIPTION:  C     !DESCRIPTION:
# Line 26  C     !USES: Line 26  C     !USES:
26  #include "DIAGNOSTICS_SIZE.h"  #include "DIAGNOSTICS_SIZE.h"
27  #include "DIAGNOSTICS.h"  #include "DIAGNOSTICS.h"
28    
29  #ifdef ALLOW_FIZHI        INTEGER NrMax
30  #include "fizhi_SIZE.h"        PARAMETER( NrMax = numLevels )
 #else  
       INTEGER Nrphys  
       PARAMETER (Nrphys=0)  
 #endif  
   
31    
32  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
33  C     listId  :: Diagnostics list number being written  C     listId  :: Diagnostics list number being written
# Line 43  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)
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)
 C     mate  :: counter mate Id number (in available diagnostics list)  
55  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
56  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
57        INTEGER i, j, k  C     mate  :: counter mate Id number (in available diagnostics list)
58    C     mDbl  :: processing mate Id number (in case processing requires 2 diags)
59    C     mVec  :: vector mate Id number
60    C     ppFld :: post-processed diag or not (=0): =1 stored in qtmp1 ; =2 in qtmp2
61    C   isComputed :: previous post-processed diag (still available in qtmp)
62    C     nLevOutp :: number of levels to write in output file
63    C
64    C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
65    C     qtmp1 :: temporary array; used to store a copy of diag. output field.
66    C     qtmp2 :: temporary array; used to store a copy of a 2nd diag. field.
67    C-  Note: local common block no longer needed.
68    c     COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
69          _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
70          _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
71    
72          INTEGER i, j, k, lm
73        INTEGER bi, bj        INTEGER bi, bj
74        INTEGER md, ndId, ip, im        INTEGER md, ndId, nn, ip, im
75        INTEGER mate, mVec        INTEGER mate, mDbl, mVec
76        CHARACTER*8 parms1        INTEGER ppFld, isComputed
77        CHARACTER*3 mate_index        CHARACTER*10 gcode
78        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)        _RL undefRL
79        _RL qtmpsrf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        INTEGER nLevOutp, kLev
       _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)  
       _RL undef, getcon  
       EXTERNAL getcon  
       INTEGER ILNBLNK  
       EXTERNAL ILNBLNK  
       INTEGER ilen  
       integer nlevsout,nplevs  
       parameter(nplevs = 16)  
       _RL plevs1(nplevs)  
       data plevs1/ 1000.0 _d 2, 925.0 _d 2, 850.0 _d 2, 700.0 _d 2,  
      .              600.0 _d 2, 500.0 _d 2, 400.0 _d 2, 300.0 _d 2,  
      .              250.0 _d 2, 200.0 _d 2, 150.0 _d 2, 100.0 _d 2,  
      .               70.0 _d 2,  50.0 _d 2,  30.0 _d 2,  20.0 _d 2/  
       _RL plevs2(nplevs)  
       data plevs2/ 1000.0 _d 2, 950.0 _d 2, 900.0 _d 2, 850.0 _d 2,  
      .              800.0 _d 2, 750.0 _d 2, 700.0 _d 2, 600.0 _d 2,  
      .              500.0 _d 2, 400.0 _d 2, 300.0 _d 2, 250.0 _d 2,  
      .              200.0 _d 2, 150.0 _d 2, 100.0 _d 2,  50.0 _d 2/  
       _RL qprs(sNx,sNy,nplevs)  
       _RL qinp(sNx,sNy,Nr+Nrphys)  
       _RL pkz(sNx,sNy,Nr+Nrphys)  
       _RL pksrf(sNx,sNy)  
       _RL p  
       INTEGER jpoint1,ipoint1  
       INTEGER jpoint2,ipoint2  
       _RL kappa  
       logical foundp  
       data foundp /.false./  
80    
81          INTEGER iLen
82        INTEGER ioUnit        INTEGER ioUnit
83        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
84        CHARACTER*(MAX_LEN_MBUF) suff        CHARACTER*(MAX_LEN_MBUF) suff
85        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
86          INTEGER prec, nRec, nTimRec
87          _RL     timeRec(2)
88          _RL     tmpLoc
89    #ifdef ALLOW_MDSIO
90        LOGICAL glf        LOGICAL glf
91    #endif
92  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
       INTEGER ii  
93        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
       CHARACTER*(5) ctmp  
       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  
       _RS ztmp(Nr+Nrphys)  
94  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
95    
96  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
97    
98    C---  set file properties
99        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
100        undef = getcon('UNDEF')        undefRL = UNSET_RL
101        kappa = getcon('KAPPA')  #ifdef ALLOW_FIZHI
102        glf = globalFiles        IF ( useFIZHI ) undefRL = getcon('UNDEF')
103    #endif
104          IF ( misvalFlt(listId).NE.UNSET_RL ) undefRL = misvalFlt(listId)
105    
106        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
107        ilen = ILNBLNK(fnames(listId))        iLen = ILNBLNK(fnames(listId))
108        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
109    C-    for now, if integrate vertically, output field has just 1 level:
110          nLevOutp = nlevels(listId)
111          IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
112    
113    C--   Set time information:
114          IF ( freq(listId).LT.0. ) THEN
115    C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
116            nTimRec = 1
117            timeRec(1) = myTime
118          ELSE
119    C-    Time-average: store the 2 edges of the time-averaging interval.
120    C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
121    C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
122            nTimRec = 2
123    
124    C-    end of time-averaging interval:
125            timeRec(2) = myTime
126    
127    C-    begining of time-averaging interval:
128    c       timeRec(1) = myTime - freq(listId)
129    C     a) find the time of the previous multiple of output freq:
130            timeRec(1) = myTime-deltaTClock*0.5 _d 0
131            timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
132            i = INT( timeRec(1) )
133            IF ( timeRec(1).LT.0. ) THEN
134              tmpLoc = FLOAT(i)
135              IF ( timeRec(1).NE.tmpLoc ) i = i - 1
136            ENDIF
137            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
138    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
139            timeRec(1) = MAX( timeRec(1), startTime )
140    
141    C     b) round off to nearest multiple of time-step:
142            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
143            i = NINT( timeRec(1) )
144    C     if just half way, NINT will return the next time-step: correct this
145            tmpLoc = FLOAT(i) - 0.5 _d 0
146            IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
147            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
148    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
149          ENDIF
150    C--   Convert time to iteration number (debug)
151    c     DO i=1,nTimRec
152    c       timeRec(i) = timeRec(i)/deltaTClock
153    c     ENDDO
154    
155  #ifdef ALLOW_MNC  C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
156        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)  
   
         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       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  */  
157    
158        ENDIF  #ifdef ALLOW_MNC
159           IF (useMNC .AND. diag_mnc) THEN
160             CALL DIAGNOSTICS_MNC_SET(
161         I                    nLevOutp, listId, lm,
162         O                    diag_mnc_bn,
163         I                    undefRL, myTime, myIter, myThid )
164           ENDIF
165  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
166    
167        DO md = 1,nfields(listId)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
168    
169           isComputed = 0
170           DO md = 1,nfields(listId)
171          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
172          parms1 = gdiag(ndId)(1:8)          gcode = gdiag(ndId)(1:10)
173          IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN          mate = 0
174            mVec = 0
175            mDbl = 0
176            ppFld = 0
177            IF ( gcode(5:5).EQ.'C' ) THEN
178    C-      Check for Mate of a Counter Diagnostic
179               mate = hdiag(ndId)
180            ELSEIF ( gcode(5:5).EQ.'P' ) THEN
181               ppFld = 1
182               IF ( gdiag(hdiag(ndId))(5:5).EQ.'P' ) ppFld = 2
183    C-      Also load the mate (if stored) for Post-Processing
184               nn = ndId
185               DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
186                 nn = hdiag(nn)
187               ENDDO
188               IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
189    c          write(0,*) ppFld,' ndId=', ndId, nn, mDbl, isComputed
190            ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
191    C-      Check for Mate of a Vector Diagnostic
192               mVec = hdiag(ndId)
193            ENDIF
194            IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
195  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
196    
197            ip = ABS(idiag(md,listId))            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
198            im = mdiag(md,listId)            im = mdiag(md,listId)
199            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
200              IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
201              IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
202    
203              IF ( ppFld.EQ.2 .AND. isComputed.EQ.hdiag(ndId) ) THEN
204    C-        Post-Processed diag from an other Post-Processed diag -and-
205    C         both of them have just been calculated and are still stored in qtmp:
206    C         => skip computation and just write qtmp2
207                IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
208                   WRITE(ioUnit,'(A,I6,3A,I6)')
209         &         '  get Post-Proc. Diag # ', ndId, '  ', cdiag(ndId),
210         &         ' from previous computation of Diag # ', isComputed
211                ENDIF
212                isComputed = 0
213              ELSEIF ( ndiag(ip,1,1).EQ.0 ) THEN
214  C-        Empty diagnostics case :  C-        Empty diagnostics case :
215                isComputed = 0
216    
217              _BEGIN_MASTER( myThid )              _BEGIN_MASTER( myThid )
218              WRITE(msgBuf,'(A,I10)')              WRITE(msgBuf,'(A,I10)')
219       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
220              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
221       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
222              WRITE(msgBuf,'(A,I4,3A,I3,2A)')              WRITE(msgBuf,'(A,I6,3A,I4,2A)')
223       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
224       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
225              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
226       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
227              WRITE(msgBuf,'(A,I2,A)')              IF ( averageCycle(listId).GT.1 ) THEN
228       &       '- WARNING -   has not been filled (ndiag=',               WRITE(msgBuf,'(A,2(I3,A))')
229       &       ndiag(ip,1,1), ' )'       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
230         &                                            ndiag(ip,1,1), ' )'
231                ELSE
232                 WRITE(msgBuf,'(A,2(I3,A))')
233         &        '- WARNING -   has not been filled (ndiag=',
234         &                                            ndiag(ip,1,1), ' )'
235                ENDIF
236              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
237       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
238              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
# Line 238  C-        Empty diagnostics case : Line 242  C-        Empty diagnostics case :
242              _END_MASTER( myThid )              _END_MASTER( myThid )
243              DO bj = myByLo(myThid), myByHi(myThid)              DO bj = myByLo(myThid), myByHi(myThid)
244                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
245                  DO k = 1,nlevels(listId)                  DO k = 1,nLevOutp
246                    DO j = 1-OLy,sNy+OLy                    DO j = 1-OLy,sNy+OLy
247                      DO i = 1-OLx,sNx+OLx                      DO i = 1-OLx,sNx+OLx
248                        qtmp1(i,j,k,bi,bj) = 0. _d 0                        qtmp1(i,j,k,bi,bj) = 0. _d 0
# Line 250  C-        Empty diagnostics case : Line 254  C-        Empty diagnostics case :
254    
255            ELSE            ELSE
256  C-        diagnostics is not empty :  C-        diagnostics is not empty :
257                isComputed = 0
258    
259              IF ( myThid.EQ.1 ) WRITE(ioUnit,'(A,I3,3A,I8,2A)')              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
260                  IF ( ppFld.GE.1 ) THEN
261                   WRITE(ioUnit,'(A,I6,7A,I8,2A)')
262         &         ' Post-Processing Diag # ', ndId, '  ', cdiag(ndId),
263         &         '   Parms: ',gdiag(ndId)
264                   IF ( mDbl.EQ.0 ) THEN
265                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
266         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
267                   ELSE
268                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
269         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
270         &          ' and diag: ',
271         &            cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
272                   ENDIF
273                  ELSE
274                   WRITE(ioUnit,'(A,I6,3A,I8,2A)')
275       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
276       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
277                  ENDIF
278              IF ( parms1(5:5).EQ.'C' ) THEN                IF ( mate.GT.0 ) THEN
279  C             Check for Mate of a Counter Diagnostic                 WRITE(ioUnit,'(3A,I6,2A)')
 C             --------------------------------------  
               mate_index = parms1(6:8)  
               READ (mate_index,'(I3)') mate  
               IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,2A)')  
280       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
281       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
282                  ELSEIF ( mVec.GT.0 ) THEN
             ELSE  
               mate = 0  
   
 C             Check for Mate of a Vector Diagnostic  
 C             -------------------------------------  
               IF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN  
                 mate_index = parms1(6:8)  
                 READ (mate_index,'(I3)') mVec  
283                  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
284                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
285       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
286       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
287       &             ' exists '       &             ' exists '
288                  ELSE                  ELSE
289                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
290       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
291       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
292       &             ' not enabled'       &             ' not enabled'
# Line 286  C             -------------------------- Line 294  C             --------------------------
294                ENDIF                ENDIF
295              ENDIF              ENDIF
296    
297              DO bj = myByLo(myThid), myByHi(myThid)              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
298               DO bi = myBxLo(myThid), myBxHi(myThid)  C-       get only selected levels:
299                DO k = 1,nlevels(listId)                DO bj = myByLo(myThid), myByHi(myThid)
300                  CALL GETDIAG(                 DO bi = myBxLo(myThid), myBxHi(myThid)
301       I                       levs(k,listId),undef,                  DO k = 1,nlevels(listId)
302       O                       qtmp1(1-OLx,1-OLy,k,bi,bj),                    kLev = NINT(levs(k,listId))
303       I                       ndId,mate,ip,im,bi,bj,myThid)                    CALL DIAGNOSTICS_GET_DIAG(
304         I                         kLev, undefRL,
305         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
306         I                         ndId, mate, ip, im, bi, bj, myThid )
307                    ENDDO
308                   ENDDO
309                ENDDO                ENDDO
310               ENDDO                IF ( mDbl.GT.0 ) THEN
311              ENDDO                 DO bj = myByLo(myThid), myByHi(myThid)
312                    DO bi = myBxLo(myThid), myBxHi(myThid)
313  C-        end of empty diag / not empty block                   DO k = 1,nlevels(listId)
314            ENDIF                    kLev = NINT(levs(k,listId))
315                      CALL DIAGNOSTICS_GET_DIAG(
316            nlevsout = nlevels(listId)       I                         kLev, undefRL,
317         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
318         I                         mDbl, 0, im, 0, bi, bj, myThid )
319                     ENDDO
320                    ENDDO
321                   ENDDO
322                  ENDIF
323                ELSE
324    C-       get all the levels (for vertical post-processing)
325                  DO bj = myByLo(myThid), myByHi(myThid)
326                   DO bi = myBxLo(myThid), myBxHi(myThid)
327                      CALL DIAGNOSTICS_GET_DIAG(
328         I                         0, undefRL,
329         O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
330         I                         ndId, mate, ip, im, bi, bj, myThid )
331                   ENDDO
332                  ENDDO
333                  IF ( mDbl.GT.0 ) THEN
334                   DO bj = myByLo(myThid), myByHi(myThid)
335                    DO bi = myBxLo(myThid), myBxHi(myThid)
336                      CALL DIAGNOSTICS_GET_DIAG(
337         I                         0, undefRL,
338         O                         qtmp2(1-OLx,1-OLy,1,bi,bj),
339         I                         mDbl, 0, im, 0, bi, bj, myThid )
340                    ENDDO
341                   ENDDO
342                  ENDIF
343                ENDIF
344    
345  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
346  C             Check to see if we need to interpolate before output  C--     Apply specific post-processing (e.g., interpolate) before output
347  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
348  C  (we are still inside field exist if sequence and field do loop)              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
349  C  C-          Do vertical interpolation:
350                 IF ( fluidIsAir ) THEN
351           if(fflags(listId)(2:2).eq.'P') then  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
352                  CALL DIAGNOSTICS_INTERP_VERT(
353  c If nonlinear free surf is active, need averaged pressures       I                         listId, md, ndId, ip, im, lm,
354  #ifdef NONLIN_FRSURF       U                         qtmp1, qtmp2,
355            if(select_rStar.GT.0)then       I                         undefRL, myTime, myIter, myThid )
356             call diagnostics_get_pointers('RSURF   ',ipoint1,jpoint1,               ELSE
357       .                                                           myThid)                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
358             call diagnostics_get_pointers('PRESSURE',ipoint2,jpoint2,       &           'INTERP_VERT not allowed in this config'
359       .                                                           myThid)                 CALL PRINT_ERROR( msgBuf , myThid )
360  C if fizhi is being  used, may need to get physics grid pressures                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
361  #ifdef ALLOW_FIZHI               ENDIF
362             if(gdiag(ndId)(10:10) .EQ. 'L')then              ENDIF
363             call diagnostics_get_pointers('FIZPRES ',ipoint2,jpoint2,              IF ( fflags(listId)(2:2).EQ.'I' ) THEN
364       .                                                           myThid)  C-          Integrate vertically: for now, output field has just 1 level:
365             endif                CALL DIAGNOSTICS_SUM_LEVELS(
366  #endif       I                         listId, md, ndId, ip, im, lm,
367             if( jpoint1.ne.0 .and. jpoint2.ne.0) foundp = .true.       U                         qtmp1,
368         I                         undefRL, myTime, myIter, myThid )
369             if(.not. foundp) then              ENDIF
370              WRITE(msgBuf,'(3A)') 'DIAGNOSTICS_OUT: ',              IF ( ppFld.GE.1 ) THEN
371       .    ' Have asked for pressure interpolation but have not ',  C-          Do Post-Processing:
372       .    ' Activated surface and 3D pressure diagnostic, ',               IF ( flds(md,listId).EQ.'PhiVEL  '
373       .    ' RSURF and PRESSURE'       &       .OR. flds(md,listId).EQ.'PsiVEL  '
374              CALL PRINT_ERROR( msgBuf , myThid )       &          ) THEN
375              STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'                CALL DIAGNOSTICS_CALC_PHIVEL(
376             endif       I                         listId, md, ndId, ip, im, lm,
377         I                         NrMax,
378             DO bj = myByLo(myThid), myByHi(myThid)       U                         qtmp1, qtmp2,
379              DO bi = myBxLo(myThid), myBxHi(myThid)       I                         myTime, myIter, myThid )
380               call getdiag(1,undef,qtmpsrf(1-OLx,1-OLy,bi,bj),                isComputed = ndId
381       .                       jpoint1,0,ipoint1,0,bi,bj,myThid)               ELSE
382              ENDDO                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
383             ENDDO       &           'unknown Processing method for diag="',cdiag(ndId),'"'
384             DO bj = myByLo(myThid), myByHi(myThid)                 CALL PRINT_ERROR( msgBuf , myThid )
385              DO bi = myBxLo(myThid), myBxHi(myThid)                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
386               DO k = 1,nlevels(listId)               ENDIF
387                call getdiag(levs(k,listId),undef,              ENDIF
      .          qtmp2(1-OLx,1-OLy,k,bi,bj),jpoint2,0,ipoint2,0,  
      .          bi,bj,myThid)  
              ENDDO  
             ENDDO  
            ENDDO  
           endif  
 #else  
 C If nonlinear free surf is off, get pressures from rC and rF arrays  
           DO bj = myByLo(myThid), myByHi(myThid)  
            DO bi = myBxLo(myThid), myBxHi(myThid)  
             DO j = 1-OLy,sNy+OLy  
              DO i = 1-OLx,sNx+OLx  
               qtmpsrf(i,j,bi,bj) = rF(1)  
              ENDDO  
             ENDDO  
             DO j = 1-OLy,sNy+OLy  
              DO i = 1-OLx,sNx+OLx  
               DO k = 1,nlevels(listId)  
                qtmp2(i,j,k,bi,bj) = rC(k)  
               ENDDO  
              ENDDO  
             ENDDO  
            ENDDO  
           ENDDO  
 #endif  
 C Load p to the kappa into a temporary array  
           nlevsout = nplevs  
           DO bj = myByLo(myThid), myByHi(myThid)  
            DO bi = myBxLo(myThid), myBxHi(myThid)  
             DO j = 1,sNy  
              DO i = 1,sNx  
               pksrf(i,j) = qtmpsrf(i,j,bi,bj) ** kappa  
               DO k = 1,nlevels(listId)  
                if(gdiag(ndId)(10:10).eq.'R') then  
                 if(hFacC(i,j,nlevels(listId)-k+1,bi,bj).ne.0.) then  
                  qinp(i,j,k) =  qtmp1(i,j,nlevels(listId)-k+1,bi,bj)  
                 else  
                  qinp(i,j,k) =  undef  
                 endif  
                 pkz(i,j,k) = qtmp2(i,j,nlevels(listId)-k+1,bi,bj)**kappa  
                elseif(gdiag(ndId)(10:10).eq.'L') then  
                 qinp(i,j,k) =  qtmp1(i,j,k,bi,bj)  
                 pkz(i,j,k) = qtmp2(i,j,k,bi,bj)**kappa  
                endif  
               ENDDO  
              ENDDO  
             ENDDO  
   
             DO k = 1,nplevs  
              if(fflags(listId)(3:3).eq.'1') then  
               p = plevs1(k)  
              elseif(fflags(listId)(3:3).eq.'2')then  
               p = plevs2(k)  
              endif  
              call prestopres(qprs(1,1,k),qinp,pkz,pksrf,0.,p,sNx,sNy,  
      .                                                 nlevels(listId) )  
             ENDDO  
388    
389              DO j = 1,sNy  C--     End of empty diag / not-empty diag block
390               DO i = 1,sNx            ENDIF
               DO k = 1,nlevsout  
                qtmp1(i,j,k,bi,bj) =  qprs(i,j,k)  
                if(qtmp1(i,j,k,bi,bj).eq.undef) qtmp1(i,j,k,bi,bj) = 0.  
               ENDDO  
              ENDDO  
             ENDDO  
            ENDDO  
           ENDDO  
391    
392           endif  C--     Ready to write field "md", element "lm" in averageCycle(listId)
393    
394  #ifdef ALLOW_MDSIO  C-        write to binary file, using MDSIO pkg:
395  C         Prepare for mdsio optionality            IF ( diag_mdsio ) THEN
396            IF (diag_mdsio) THEN  c          nRec = lm + (md-1)*averageCycle(listId)
397              IF (fflags(listId)(1:1) .EQ. ' ') THEN             nRec = md + (lm-1)*nfields(listId)
398  C             This is the old default behavior  C         default precision for output files
399                CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,'RL',             prec = writeBinaryPrec
400       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)  C         fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
401              ELSEIF (fflags(listId)(1:1) .EQ. 'R') THEN             IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
402  C             Force it to be 32-bit precision             IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
403                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,'RL',  C         a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R
404       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)             IF ( ppFld.LE.1 ) THEN
405              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              CALL WRITE_REC_LEV_RL(
406  C             Force it to be 64-bit precision       I                            fn, prec,
407                CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,'RL',       I                            NrMax, 1, nLevOutp,
408       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       I                            qtmp1, -nRec, myIter, myThid )
409              ENDIF             ELSE
410                CALL WRITE_REC_LEV_RL(
411         I                            fn, prec,
412         I                            NrMax, 1, nLevOutp,
413         I                            qtmp2, -nRec, myIter, myThid )
414               ENDIF
415            ENDIF            ENDIF
 #endif /*  ALLOW_MDSIO  */  
416    
417  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
418            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
419               IF ( ppFld.LE.1 ) THEN
420              _BEGIN_MASTER( myThid )              CALL DIAGNOSTICS_MNC_OUT(
421         I                       NrMax, nLevOutp, listId, ndId, mate,
422              DO ii = 1,CW_DIMS       I                       diag_mnc_bn, qtmp1,
423                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)       I                       undefRL, myTime, myIter, myThid )
424                dn(ii)(1:NLEN) = dn_blnk(1:NLEN)             ELSE
425              ENDDO              CALL DIAGNOSTICS_MNC_OUT(
426         I                       NrMax, nLevOutp, listId, ndId, mate,
427  C           Note that the "d_cw_name" variable is a hack that hides a       I                       diag_mnc_bn, qtmp2,
428  C           subtlety within MNC.  Basically, each MNC-wrapped file is       I                       undefRL, myTime, myIter, myThid )
429  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) = Nr+Nrphys  
             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)  
   
             IF ((fflags(listId)(1:1) .EQ. ' ')  
      &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN  
               CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,  
      &             cdiag(ndId), qtmp1, myThid)  
             ELSEIF (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 )  
   
430            ENDIF            ENDIF
431  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
432    
433  C--     end of Processing Fld # md  C--     end of Processing Fld # md
434          ENDIF          ENDIF
435           ENDDO
436    
437    C--   end loop on lm counter (= averagePeriod)
438        ENDDO        ENDDO
439    
440    #ifdef ALLOW_MDSIO
441          IF (diag_mdsio) THEN
442    C-    Note: temporary: since it is a pain to add more arguments to
443    C     all MDSIO S/R, uses instead this specific S/R to write only
444    C     meta files but with more informations in it.
445                glf = globalFiles
446                nRec = averageCycle(listId)*nfields(listId)
447                CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
448         &              0, 0, nLevOutp, ' ',
449         &              nfields(listId), flds(1,listId),
450         &              nTimRec, timeRec, undefRL,
451         &              nRec, myIter, myThid)
452          ENDIF
453    #endif /*  ALLOW_MDSIO  */
454    
455        RETURN        RETURN
456        END        END
457    

Legend:
Removed from v.1.18  
changed lines
  Added in v.1.60

  ViewVC Help
Powered by ViewVC 1.1.22