/[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.49 by jmc, Mon Jun 6 15:42:58 2011 UTC revision 1.56 by jmc, Thu Jun 23 15:29:01 2011 UTC
# Line 10  C     !ROUTINE: DIAGNOSTICS_OUT Line 10  C     !ROUTINE: DIAGNOSTICS_OUT
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 58  C     im    :: counter-mate pointer to s Line 58  C     im    :: counter-mate pointer to s
58  C     nLevOutp :: number of levels to write in output file  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*10 gcode        CHARACTER*10 gcode
73        _RL undef        _RL undefRL
74        _RL tmpLev        INTEGER nLevOutp, kLev
       INTEGER ilen  
       INTEGER nLevOutp  
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, nTimRec        INTEGER prec, nRec, nTimRec
82        _RL     timeRec(2)        _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 ll, llMx, jj, jjMx  
       INTEGER ii, klev  
88        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
89        INTEGER CW_DIMS, NLEN        LOGICAL missingValFillsMask
       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)  
90  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
91    
92  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
93    
94  C---  set file properties  C---  set file properties
95        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
96        undef = UNSET_RL        undefRL = UNSET_RL
97  #ifdef ALLOW_FIZHI  #ifdef ALLOW_FIZHI
98  c     IF ( useFIZHI ) undef = getcon('UNDEF')        IF ( useFIZHI ) undefRL = getcon('UNDEF')
       undef = getcon('UNDEF')  
99  #endif  #endif
100          IF ( misvalFlt(listId).NE.UNSET_RL ) undefRL = misvalFlt(listId)
101    
102        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
103        ilen = ILNBLNK(fnames(listId))        iLen = ILNBLNK(fnames(listId))
104        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
105  C-    for now, if integrate vertically, output field has just 1 level:  C-    for now, if integrate vertically, output field has just 1 level:
106        nLevOutp = nlevels(listId)        nLevOutp = nlevels(listId)
107        IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1        IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
# Line 142  C     a) find the time of the previous m Line 127  C     a) find the time of the previous m
127          timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)          timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
128          i = INT( timeRec(1) )          i = INT( timeRec(1) )
129          IF ( timeRec(1).LT.0. ) THEN          IF ( timeRec(1).LT.0. ) THEN
130            tmpLev = FLOAT(i)            tmpLoc = FLOAT(i)
131            IF ( timeRec(1).NE.tmpLev ) i = i - 1            IF ( timeRec(1).NE.tmpLoc ) i = i - 1
132          ENDIF          ENDIF
133          timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)          timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
134  c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock  c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
# Line 153  C     b) round off to nearest multiple o Line 138  C     b) round off to nearest multiple o
138          timeRec(1) = (timeRec(1)-baseTime)/deltaTClock          timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
139          i = NINT( timeRec(1) )          i = NINT( timeRec(1) )
140  C     if just half way, NINT will return the next time-step: correct this  C     if just half way, NINT will return the next time-step: correct this
141          tmpLev = FLOAT(i) - 0.5 _d 0          tmpLoc = FLOAT(i) - 0.5 _d 0
142          IF ( timeRec(1).EQ.tmpLev ) i = i - 1          IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
143          timeRec(1) = baseTime + deltaTClock*FLOAT(i)          timeRec(1) = baseTime + deltaTClock*FLOAT(i)
144  c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock  c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
145        ENDIF        ENDIF
# Line 163  c     DO i=1,nTimRec Line 148  c     DO i=1,nTimRec
148  c       timeRec(i) = timeRec(i)/deltaTClock  c       timeRec(i) = timeRec(i)/deltaTClock
149  c     ENDDO  c     ENDDO
150    
151  #ifdef ALLOW_MNC  C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
152  C-- this is a trick to reverse the order of the loops on md (= field)        DO lm=1,averageCycle(listId)
 C   and lm (= averagePeriod): binary output: lm loop inside md loop ;  
 C                                 mnc ouput: md loop inside lm loop.  
       IF (useMNC .AND. diag_mnc) THEN  
         jjMx = averageCycle(listId)  
         llMx = 1  
       ELSE  
         jjMx = 1  
         llMx = averageCycle(listId)  
       ENDIF  
       DO jj=1,jjMx  
153    
154    #ifdef ALLOW_MNC
155         IF (useMNC .AND. diag_mnc) THEN         IF (useMNC .AND. diag_mnc) THEN
156  C     Handle missing value attribute (land points)           CALL DIAGNOSTICS_MNC_SET(
157           useMissingValue = .FALSE.       I                    nLevOutp, listId, lm,
158  #ifdef DIAGNOSTICS_MISSING_VALUE       O                    diag_mnc_bn, missingValFillsMask,
159           useMissingValue = .TRUE.       I                    undefRL, myTime, myIter, myThid )
 #endif /* DIAGNOSTICS_MISSING_VALUE */  
          IF ( misvalFlt(listId) .NE. UNSET_RL ) THEN  
           misvalLoc = misvalFlt(listId)  
          ELSE  
           misvalLoc = undef  
          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  
          klev = myIter + jj - jjMx  
          tmpLev = myTime + deltaTClock*(jj -jjMx)  
          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)  
          CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',tmpLev,myThid)  
          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)  
          CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',klev,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', nLevOutp  
          dim(1) = nLevOutp  
          ib(1)  = 1  
          ie(1)  = nLevOutp  
   
          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  */  
   
160         ENDIF         ENDIF
161  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
162    
# Line 313  C---+----1----+----2----+----3----+----4 Line 167  C---+----1----+----2----+----3----+----4
167          gcode = gdiag(ndId)(1:10)          gcode = gdiag(ndId)(1:10)
168          mate = 0          mate = 0
169          mVec = 0          mVec = 0
170            mDbl = 0
171          IF ( gcode(5:5).EQ.'C' ) THEN          IF ( gcode(5:5).EQ.'C' ) THEN
172  C-      Check for Mate of a Counter Diagnostic  C-      Check for Mate of a Counter Diagnostic
173             mate = hdiag(ndId)             mate = hdiag(ndId)
174            ELSEIF ( gcode(5:5).EQ.'P' ) THEN
175    C-      Also load the mate (if stored) for Post-Processing
176               nn = ndId
177               DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
178                 nn = hdiag(nn)
179               ENDDO
180               IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
181          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
182  C-      Check for Mate of a Vector Diagnostic  C-      Check for Mate of a Vector Diagnostic
183             mVec = hdiag(ndId)             mVec = hdiag(ndId)
184          ENDIF          ENDIF
185          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
186  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
 #ifdef ALLOW_MNC  
          DO ll=1,llMx  
           lm = jj+ll-1  
 #else  
          DO lm=1,averageCycle(listId)  
 #endif  
187    
188            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
189            im = mdiag(md,listId)            im = mdiag(md,listId)
190            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
191              IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
192            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
193    
194            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ndiag(ip,1,1).EQ.0 ) THEN
# Line 379  C-        Empty diagnostics case : Line 236  C-        Empty diagnostics case :
236  C-        diagnostics is not empty :  C-        diagnostics is not empty :
237    
238              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
239                WRITE(ioUnit,'(A,I6,3A,I8,2A)')                IF ( gcode(5:5).EQ.'P' ) THEN
240                   WRITE(ioUnit,'(A,I6,7A,I8,2A)')
241         &         ' Post-Processing Diag # ', ndId, '  ', cdiag(ndId),
242         &         '   Parms: ',gdiag(ndId)
243                   IF ( mDbl.EQ.0 ) THEN
244                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
245         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
246                   ELSE
247                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
248         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
249         &          ' and diag: ',
250         &            cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
251                   ENDIF
252                  ELSE
253                   WRITE(ioUnit,'(A,I6,3A,I8,2A)')
254       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
255       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
256                  ENDIF
257                IF ( mate.GT.0 ) THEN                IF ( mate.GT.0 ) THEN
258                 WRITE(ioUnit,'(3A,I6,2A)')                 WRITE(ioUnit,'(3A,I6,2A)')
259       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
# Line 401  C-        diagnostics is not empty : Line 273  C-        diagnostics is not empty :
273                ENDIF                ENDIF
274              ENDIF              ENDIF
275    
276              IF ( fflags(listId)(2:2).NE.' ' ) THEN              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
277  C-       get all the levels (for vertical post-processing)  C-       get only selected levels:
278                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
279                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
280                  DO k = 1,kdiag(ndId)                  DO k = 1,nlevels(listId)
281                    tmpLev = k                    kLev = NINT(levs(k,listId))
282                    CALL GETDIAG(                    CALL DIAGNOSTICS_GET_DIAG(
283       I                         tmpLev,undef,       I                         kLev, undefRL,
284       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
285       I                         ndId,mate,ip,im,bi,bj,myThid)       I                         ndId, mate, ip, im, bi, bj, myThid )
286                  ENDDO                  ENDDO
287                 ENDDO                 ENDDO
288                ENDDO                ENDDO
289                  IF ( mDbl.GT.0 ) THEN
290                   DO bj = myByLo(myThid), myByHi(myThid)
291                    DO bi = myBxLo(myThid), myBxHi(myThid)
292                     DO k = 1,nlevels(listId)
293                      kLev = NINT(levs(k,listId))
294                      CALL DIAGNOSTICS_GET_DIAG(
295         I                         kLev, undefRL,
296         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
297         I                         mDbl, 0, im, 0, bi, bj, myThid )
298                     ENDDO
299                    ENDDO
300                   ENDDO
301                  ENDIF
302              ELSE              ELSE
303  C-       get only selected levels:  C-       get all the levels (for vertical post-processing)
304                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
305                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
306                  DO k = 1,nlevels(listId)                    CALL DIAGNOSTICS_GET_DIAG(
307                    CALL GETDIAG(       I                         0, undefRL,
308       I                         levs(k,listId),undef,       O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
309       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  
310                 ENDDO                 ENDDO
311                ENDDO                ENDDO
312                  IF ( mDbl.GT.0 ) THEN
313                   DO bj = myByLo(myThid), myByHi(myThid)
314                    DO bi = myBxLo(myThid), myBxHi(myThid)
315                     DO k = 1,nlevels(listId)
316                      CALL DIAGNOSTICS_GET_DIAG(
317         I                         0, undefRL,
318         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
319         I                         mDbl, 0, im, 0, bi, bj, myThid )
320                     ENDDO
321                    ENDDO
322                   ENDDO
323                  ENDIF
324              ENDIF              ENDIF
325    
326  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
# Line 437  C-          Do vertical interpolation: Line 332  C-          Do vertical interpolation:
332  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);
333                CALL DIAGNOSTICS_INTERP_VERT(                CALL DIAGNOSTICS_INTERP_VERT(
334       I                         listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
335       U                         qtmp1,       U                         qtmp1, qtmp2,
336       I                         undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
337               ELSE               ELSE
338                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
339       &           'INTERP_VERT not allowed in this config'       &           'INTERP_VERT not allowed in this config'
# Line 451  C-          Integrate vertically: for no Line 346  C-          Integrate vertically: for no
346                CALL DIAGNOSTICS_SUM_LEVELS(                CALL DIAGNOSTICS_SUM_LEVELS(
347       I                         listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
348       U                         qtmp1,       U                         qtmp1,
349       I                         undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
350                ENDIF
351                IF ( gcode(5:5).EQ.'P' ) THEN
352    C-          Do Post-Processing:
353                 IF ( flds(md,listId).EQ.'PhiVEL  '
354    c    &       .OR. flds(md,listId).EQ.'PsiVEL  '
355         &          ) THEN
356                  CALL DIAGNOSTICS_CALC_PHIVEL(
357         I                         listId, md, ndId, ip, im, lm,
358         U                         qtmp1, qtmp2,
359         I                         myTime, myIter, myThid )
360                 ELSE
361                   WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
362         &           'unknown Processing method for diag="',cdiag(ndId),'"'
363                   CALL PRINT_ERROR( msgBuf , myThid )
364                   STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
365                 ENDIF
366              ENDIF              ENDIF
367    
368  C--     End of empty diag / not-empty diag block  C--     End of empty diag / not-empty diag block
# Line 461  C--     Ready to write field "md", eleme Line 372  C--     Ready to write field "md", eleme
372    
373  C-        write to binary file, using MDSIO pkg:  C-        write to binary file, using MDSIO pkg:
374            IF ( diag_mdsio ) THEN            IF ( diag_mdsio ) THEN
375              nRec = lm + (md-1)*averageCycle(listId)  c           nRec = lm + (md-1)*averageCycle(listId)
376                nRec = md + (lm-1)*nfields(listId)
377  C           default precision for output files  C           default precision for output files
378              prec = writeBinaryPrec              prec = writeBinaryPrec
379  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 476  C         a hack not to write meta files Line 388  C         a hack not to write meta files
388    
389  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
390            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
391                CALL DIAGNOSTICS_MNC_OUT(
392              _BEGIN_MASTER( myThid )       I                       NrMax, nLevOutp, listId, ndId, mate,
393         I                       diag_mnc_bn,
394              DO ii = 1,CW_DIMS       I                       missingValFillsMask, undefRL,
395                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)       I                       qtmp1,
396                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', nLevOutp  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zmd', nLevOutp  
             ENDIF  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zld', nLevOutp  
             ENDIF  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zud', nLevOutp  
             ENDIF  
             dim(3) = NrMax  
             ib(3)  = 1  
             ie(3)  = nLevOutp  
   
 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  
 C     note: better to use 2-D mask if kdiag <> Nr or vert.integral  
              DO bj = myByLo(myThid), myByHi(myThid)  
               DO bi = myBxLo(myThid), myBxHi(myThid)  
                DO k = 1,nLevOutp  
                 klev = NINT(levs(k,listId))  
                 IF ( fflags(listId)(2:2).EQ.'I' ) kLev = 1  
                 DO j = 1-OLy,sNy+OLy  
                  DO i = 1-OLx,sNx+OLx  
                   IF ( maskC(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'))  
      &             .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 )  
   
397            ENDIF            ENDIF
398  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
399    
 C--      end loop on lm (or ll if ALLOW_MNC) counter  
          ENDDO  
400  C--     end of Processing Fld # md  C--     end of Processing Fld # md
401          ENDIF          ENDIF
402         ENDDO         ENDDO
403    
404  #ifdef ALLOW_MNC  C--   end loop on lm counter (= averagePeriod)
 C--   end loop on jj counter  
405        ENDDO        ENDDO
 #endif  
406    
407  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
408        IF (diag_mdsio) THEN        IF (diag_mdsio) THEN
# Line 625  C-    Note: temporary: since it is a pai Line 410  C-    Note: temporary: since it is a pai
410  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
411  C     meta files but with more informations in it.  C     meta files but with more informations in it.
412              glf = globalFiles              glf = globalFiles
413              nRec = nfields(listId)*averageCycle(listId)              nRec = averageCycle(listId)*nfields(listId)
414              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
415       &              0, 0, nLevOutp, ' ',       &              0, 0, nLevOutp, ' ',
416       &              nfields(listId), flds(1,listId), nTimRec, timeRec,       &              nfields(listId), flds(1,listId), nTimRec, timeRec,
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.56
  ViewVC Help
Powered by ViewVC 1.1.22