/[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.39 by mlosch, Tue May 27 08:37:19 2008 UTC revision 1.58 by jmc, Fri Jul 1 18:49:35 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)
 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    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  C
64  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
65  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.
66  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.
67  C              diagnostic storage qdiag into it.  C-  Note: local common block no longer needed.
68        COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1  c     COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
69        _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)
70          _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
71    
72        INTEGER i, j, k, lm, klev        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          INTEGER ppFld, isComputed
77        CHARACTER*10 gcode        CHARACTER*10 gcode
78        _RL undef, getcon        _RL undefRL
79        _RL tmpLev        INTEGER nLevOutp, kLev
       EXTERNAL getcon  
       INTEGER ILNBLNK  
       EXTERNAL ILNBLNK  
       INTEGER ilen  
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        INTEGER prec, nRec, nTimRec
87          _RL     timeRec(2)
88          _RL     tmpLoc
89  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
90        LOGICAL glf        LOGICAL glf
91  #endif  #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
       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  
       LOGICAL useMissingValue, useMisValForThisDiag  
       REAL*8 misvalLoc  
       REAL*8 misval_r8(2)  
       REAL*4 misval_r4(2)  
       INTEGER misvalIntLoc, misval_int(2)  
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    #ifdef ALLOW_FIZHI
102          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    C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
156          DO lm=1,averageCycle(listId)
157    
158  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
159        IF (useMNC .AND. diag_mnc) THEN         IF (useMNC .AND. diag_mnc) THEN
160  C     Handle missing value attribute (land points)           CALL DIAGNOSTICS_MNC_SET(
161         useMissingValue = .FALSE.       I                    nLevOutp, listId, lm,
162  #ifdef DIAGNOSTICS_MISSING_VALUE       O                    diag_mnc_bn,
163         useMissingValue = .TRUE.       I                    undefRL, myTime, myIter, myThid )
 #endif /* DIAGNOSTICS_MISSING_VALUE */  
        IF ( misvalFlt(listId) .NE. UNSET_RL ) THEN  
         misvalLoc = misvalFlt(listId)  
        ELSE  
         misvalLoc = undef  
164         ENDIF         ENDIF
 C     Defaults to UNSET_I  
        misvalIntLoc = misvalInt(listId)  
        DO ii=1,2  
 C       misval_r4(ii)  = UNSET_FLOAT4  
 C       misval_r8(ii)  = UNSET_FLOAT8  
         misval_r4(ii)  = misvalLoc  
         misval_r8(ii)  = misvalLoc  
         misval_int(ii) = UNSET_I  
        ENDDO  
         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)  
 C     suppress the missing value attribute (iflag = 0)  
         IF (useMissingValue)  
      &       CALL MNC_CW_VATTR_MISSING('diag_levels', 0,  
      I       misval_r8, misval_r4, misval_int,  
      I       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  
 C     suppress the missing value attribute (iflag = 0)  
           IF (useMissingValue)  
      &         CALL MNC_CW_VATTR_MISSING(dn(1), 0,  
      I         misval_r8, misval_r4, misval_int,  
      I         myThid )  
           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  */  
   
       ENDIF  
165  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
166    
167  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
168    
169        DO md = 1,nfields(listId)         isComputed = 0
170           DO md = 1,nfields(listId)
171          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
172          gcode = gdiag(ndId)(1:10)          gcode = gdiag(ndId)(1:10)
173          mate = 0          mate = 0
174          mVec = 0          mVec = 0
175            mDbl = 0
176            ppFld = 0
177          IF ( gcode(5:5).EQ.'C' ) THEN          IF ( gcode(5:5).EQ.'C' ) THEN
178  C-      Check for Mate of a Counter Diagnostic  C-      Check for Mate of a Counter Diagnostic
179             mate = hdiag(ndId)             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          ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
191  C-      Check for Mate of a Vector Diagnostic  C-      Check for Mate of a Vector Diagnostic
192             mVec = hdiag(ndId)             mVec = hdiag(ndId)
193          ENDIF          ENDIF
194          IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN          IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
195  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
          DO lm=1,averageCycle(listId)  
196    
197            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
198            im = mdiag(md,listId)            im = mdiag(md,listId)
199            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)            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)            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
202    
203            IF ( ndiag(ip,1,1).EQ.0 ) THEN            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)')
# Line 286  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 298  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 ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
260                WRITE(ioUnit,'(A,I6,3A,I8,2A)')                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 ( mate.GT.0 ) THEN                IF ( mate.GT.0 ) THEN
279                 WRITE(ioUnit,'(3A,I6,2A)')                 WRITE(ioUnit,'(3A,I6,2A)')
280       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
# Line 322  C-        diagnostics is not empty : Line 294  C-        diagnostics is not empty :
294                ENDIF                ENDIF
295              ENDIF              ENDIF
296    
297              IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
298  C-       get all the levels (for vertical interpolation)  C-       get only selected levels:
299                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
300                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
301                  DO k = 1,kdiag(ndId)                  DO k = 1,nlevels(listId)
302                    tmpLev = k                    kLev = NINT(levs(k,listId))
303                    CALL GETDIAG(                    CALL DIAGNOSTICS_GET_DIAG(
304       I                         tmpLev,undef,       I                         kLev, undefRL,
305       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
306       I                         ndId,mate,ip,im,bi,bj,myThid)       I                         ndId, mate, ip, im, bi, bj, myThid )
307                  ENDDO                  ENDDO
308                 ENDDO                 ENDDO
309                ENDDO                ENDDO
310                  IF ( mDbl.GT.0 ) THEN
311                   DO bj = myByLo(myThid), myByHi(myThid)
312                    DO bi = myBxLo(myThid), myBxHi(myThid)
313                     DO k = 1,nlevels(listId)
314                      kLev = NINT(levs(k,listId))
315                      CALL DIAGNOSTICS_GET_DIAG(
316         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              ELSE
324  C-       get only selected levels:  C-       get all the levels (for vertical post-processing)
325                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
326                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
327                  DO k = 1,nlevels(listId)                    CALL DIAGNOSTICS_GET_DIAG(
328                    CALL GETDIAG(       I                         0, undefRL,
329       I                         levs(k,listId),undef,       O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
330       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  
331                 ENDDO                 ENDDO
332                ENDDO                ENDDO
333                  IF ( mDbl.GT.0 ) THEN
334                   DO bj = myByLo(myThid), myByHi(myThid)
335                    DO bi = myBxLo(myThid), myBxHi(myThid)
336                     DO k = 1,nlevels(listId)
337                      CALL DIAGNOSTICS_GET_DIAG(
338         I                         0, undefRL,
339         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
340         I                         mDbl, 0, im, 0, bi, bj, myThid )
341                     ENDDO
342                    ENDDO
343                   ENDDO
344                  ENDIF
345              ENDIF              ENDIF
346    
 C-        end of empty diag / not empty block  
           ENDIF  
   
347  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
348  C         Check to see if we need to interpolate before output  C--     Apply specific post-processing (e.g., interpolate) before output
349  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
350            IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
351  C-        Do vertical interpolation:  C-          Do vertical interpolation:
352             IF ( fluidIsAir ) THEN               IF ( fluidIsAir ) THEN
353  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);
354              CALL DIAGNOSTICS_INTERP_VERT(                CALL DIAGNOSTICS_INTERP_VERT(
355       I                     listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
356       U                     qtmp1,       U                         qtmp1, qtmp2,
357       I                     undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
358             ELSE               ELSE
359               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
360       &         'INTERP_VERT not allowed in this config'       &           'INTERP_VERT not allowed in this config'
361               CALL PRINT_ERROR( msgBuf , myThid )                 CALL PRINT_ERROR( msgBuf , myThid )
362               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
363             ENDIF               ENDIF
364                ENDIF
365                IF ( fflags(listId)(2:2).EQ.'I' ) THEN
366    C-          Integrate vertically: for now, output field has just 1 level:
367                  CALL DIAGNOSTICS_SUM_LEVELS(
368         I                         listId, md, ndId, ip, im, lm,
369         U                         qtmp1,
370         I                         undefRL, myTime, myIter, myThid )
371                ENDIF
372                IF ( ppFld.GE.1 ) THEN
373    C-          Do Post-Processing:
374                 IF ( flds(md,listId).EQ.'PhiVEL  '
375         &       .OR. flds(md,listId).EQ.'PsiVEL  '
376         &          ) THEN
377                  CALL DIAGNOSTICS_CALC_PHIVEL(
378         I                         listId, md, ndId, ip, im, lm,
379         I                         NrMax,
380         U                         qtmp1, qtmp2,
381         I                         myTime, myIter, myThid )
382                  isComputed = ndId
383                 ELSE
384                   WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
385         &           'unknown Processing method for diag="',cdiag(ndId),'"'
386                   CALL PRINT_ERROR( msgBuf , myThid )
387                   STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
388                 ENDIF
389                ENDIF
390    
391    C--     End of empty diag / not-empty diag block
392            ENDIF            ENDIF
393    
394  C--    Ready to write field "md", element "lm" in averageCycle(listId)  C--     Ready to write field "md", element "lm" in averageCycle(listId)
395    
396  C-        write to binary file, using MDSIO pkg:  C-        write to binary file, using MDSIO pkg:
397            IF ( diag_mdsio ) THEN            IF ( diag_mdsio ) THEN
398              nRec = lm + (md-1)*averageCycle(listId)  c          nRec = lm + (md-1)*averageCycle(listId)
399  C           default precision for output files             nRec = md + (lm-1)*nfields(listId)
400              prec = writeBinaryPrec  C         default precision for output files
401  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision             prec = writeBinaryPrec
402              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32  C         fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
403              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64             IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
404               IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
405  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
406               IF ( ppFld.LE.1 ) THEN
407              CALL WRITE_REC_LEV_RL(              CALL WRITE_REC_LEV_RL(
408       I                            fn, prec,       I                            fn, prec,
409       I                            NrMax, 1, nlevels(listId),       I                            NrMax, 1, nLevOutp,
410       I                            qtmp1, -nRec, myIter, myThid )       I                            qtmp1, -nRec, myIter, myThid )
411               ELSE
412                CALL WRITE_REC_LEV_RL(
413         I                            fn, prec,
414         I                            NrMax, 1, nLevOutp,
415         I                            qtmp2, -nRec, myIter, myThid )
416               ENDIF
417            ENDIF            ENDIF
418    
419  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
420            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
421               IF ( ppFld.LE.1 ) THEN
422              _BEGIN_MASTER( myThid )              CALL DIAGNOSTICS_MNC_OUT(
423         I                       NrMax, nLevOutp, listId, ndId, mate,
424              DO ii = 1,CW_DIMS       I                       diag_mnc_bn, qtmp1,
425                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)       I                       undefRL, myTime, myIter, myThid )
426                dn(ii)(1:NLEN) = dn_blnk(1:NLEN)             ELSE
427              ENDDO              CALL DIAGNOSTICS_MNC_OUT(
428         I                       NrMax, nLevOutp, listId, ndId, mate,
429  C           Note that the "d_cw_name" variable is a hack that hides a       I                       diag_mnc_bn, qtmp2,
430  C           subtlety within MNC.  Basically, each MNC-wrapped file is       I                       undefRL, myTime, myIter, myThid )
431  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     Missing values only for scalar diagnostics at mass points (so far)  
             useMisValForThisDiag = useMissingValue  
      &           .AND.gdiag(ndId)(1:2).EQ.'SM'  
             IF ( useMisValForThisDiag ) THEN  
 C     assign missing values and set flag for adding the netCDF atttibute  
              CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 2,  
      I            misval_r8, misval_r4, misval_int,  
      I            myThid )  
 C     and now use the missing values for masking out the land points  
              DO bj = myByLo(myThid), myByHi(myThid)  
               DO bi = myBxLo(myThid), myBxHi(myThid)  
                DO k = 1,nlevels(listId)  
                 klev = NINT(levs(k,listId))  
                 DO j = 1-OLy,sNy+OLy  
                  DO i = 1-OLx,sNx+OLx  
                   IF ( _hFacC(I,J,klev,bi,bj) .EQ. 0. )  
      &                 qtmp1(i,j,k,bi,bj) = misvalLoc  
                  ENDDO  
                 ENDDO  
                ENDDO  
               ENDDO  
              ENDDO  
             ELSE  
 C     suppress the missing value attribute (iflag = 0)  
 C     Note: We have to call the following subroutine for each mnc that has  
 C     been created "on the fly" by mnc_cw_add_vname and will be deleted  
 C     by mnc_cw_del_vname, because all of these variables use the same  
 C     identifier so that mnc_cw_vfmv(indv) needs to be overwritten for  
 C     each of these variables  
              CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 0,  
      I            misval_r8, misval_r4, misval_int,  
      I            myThid )  
             ENDIF  
   
             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 )  
   
432            ENDIF            ENDIF
433  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
434    
          ENDDO  
435  C--     end of Processing Fld # md  C--     end of Processing Fld # md
436          ENDIF          ENDIF
437           ENDDO
438    
439    C--   end loop on lm counter (= averagePeriod)
440        ENDDO        ENDDO
441    
442  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
443        IF (diag_mdsio) THEN        IF (diag_mdsio) THEN
444  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
445  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
446  C     meta files but with more informations in it.  C     meta files but with more informations in it.
447              glf = globalFiles              glf = globalFiles
448              nRec = nfields(listId)*averageCycle(listId)              nRec = averageCycle(listId)*nfields(listId)
449              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
450       &              0, 0, nlevels(listId), ' ',       &              0, 0, nLevOutp, ' ',
451       &              nfields(listId), flds(1,listId), 1, myTime,       &              nfields(listId), flds(1,listId), nTimRec, timeRec,
452       &              nRec, myIter, myThid)       &              nRec, myIter, myThid)
453        ENDIF        ENDIF
454  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
Legend:
Removed from v.1.39  
changed lines
  Added in v.1.58
  ViewVC Help
Powered by ViewVC 1.1.22