/[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.44 by jmc, Thu Jan 7 01:09:40 2010 UTC revision 1.57 by jmc, Mon Jun 27 22:27:23 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 48  C     !FUNCTIONS: Line 48  C     !FUNCTIONS:
48    
49  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
50  C     i,j,k :: loop indices  C     i,j,k :: loop indices
51    C     bi,bj :: tile indices
52  C     lm    :: loop index (averageCycle)  C     lm    :: loop index (averageCycle)
53  C     md    :: field number in the list "listId".  C     md    :: field number in the list "listId".
54  C     ndId  :: diagnostics  Id number (in available diagnostics list)  C     ndId  :: diagnostics  Id number (in available diagnostics list)
55  C     mate  :: counter mate Id number (in available diagnostics list)  C     mate  :: counter mate Id number (in available diagnostics list)
56  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
57  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
58    C     nLevOutp :: number of levels to write in output file
59  C  C
60  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
61  C     qtmp1 :: thread-shared temporary array (needs to be in common block):  C     qtmp1 :: temporary array; used to store a copy of diag. output field.
62  C              to write a diagnostic field to file, copy it first from (big)  C     qtmp2 :: temporary array; used to store a copy of a 2nd diag. field.
63  C              diagnostic storage qdiag into it.  C-  Note: local common block no longer needed.
64        COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1  c     COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
65        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
66          _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
67    
68        INTEGER i, j, k, lm        INTEGER i, j, k, lm
69        INTEGER bi, bj        INTEGER bi, bj
70        INTEGER md, ndId, ip, im        INTEGER md, ndId, nn, ip, im
71        INTEGER mate, mVec        INTEGER mate, mDbl, mVec
72        CHARACTER*10 gcode        CHARACTER*10 gcode
73        _RL undef        _RL undefRL
74        _RL tmpLev        INTEGER nLevOutp, kLev
       INTEGER ilen  
75    
76          INTEGER iLen
77        INTEGER ioUnit        INTEGER ioUnit
78        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
79        CHARACTER*(MAX_LEN_MBUF) suff        CHARACTER*(MAX_LEN_MBUF) suff
80        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
81        INTEGER prec, nRec, 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
       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)  
89  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
90    
91  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
92    
93  C---  set file properties  C---  set file properties
94        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
95        undef = UNSET_RL        undefRL = UNSET_RL
96  #ifdef ALLOW_FIZHI  #ifdef ALLOW_FIZHI
97  c     IF ( useFIZHI ) undef = getcon('UNDEF')        IF ( useFIZHI ) undefRL = getcon('UNDEF')
       undef = getcon('UNDEF')  
98  #endif  #endif
99          IF ( misvalFlt(listId).NE.UNSET_RL ) undefRL = misvalFlt(listId)
100    
101        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
102        ilen = ILNBLNK(fnames(listId))        iLen = ILNBLNK(fnames(listId))
103        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
104    C-    for now, if integrate vertically, output field has just 1 level:
105          nLevOutp = nlevels(listId)
106          IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
107    
108  C--   Set time information:  C--   Set time information:
109        IF ( freq(listId).LT.0. ) THEN        IF ( freq(listId).LT.0. ) THEN
# Line 135  C     a) find the time of the previous m Line 125  C     a) find the time of the previous m
125          timeRec(1) = myTime-deltaTClock*0.5 _d 0          timeRec(1) = myTime-deltaTClock*0.5 _d 0
126          timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)          timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
127          i = INT( timeRec(1) )          i = INT( timeRec(1) )
128            IF ( timeRec(1).LT.0. ) THEN
129              tmpLoc = FLOAT(i)
130              IF ( timeRec(1).NE.tmpLoc ) i = i - 1
131            ENDIF
132          timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)          timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
133  c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock  c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
134          timeRec(1) = MAX( timeRec(1), startTime )          timeRec(1) = MAX( timeRec(1), startTime )
# Line 143  C     b) round off to nearest multiple o Line 137  C     b) round off to nearest multiple o
137          timeRec(1) = (timeRec(1)-baseTime)/deltaTClock          timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
138          i = NINT( timeRec(1) )          i = NINT( timeRec(1) )
139  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
140          tmpLev = FLOAT(i) - 0.5 _d 0          tmpLoc = FLOAT(i) - 0.5 _d 0
141          IF ( timeRec(1).EQ.tmpLev ) i = i - 1          IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
142          timeRec(1) = baseTime + deltaTClock*FLOAT(i)          timeRec(1) = baseTime + deltaTClock*FLOAT(i)
143  c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock  c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
144        ENDIF        ENDIF
145    C--   Convert time to iteration number (debug)
146    c     DO i=1,nTimRec
147    c       timeRec(i) = timeRec(i)/deltaTClock
148    c     ENDDO
149    
150  #ifdef ALLOW_MNC  C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
151  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  
152    
153    #ifdef ALLOW_MNC
154         IF (useMNC .AND. diag_mnc) THEN         IF (useMNC .AND. diag_mnc) THEN
155  C     Handle missing value attribute (land points)           CALL DIAGNOSTICS_MNC_SET(
156           useMissingValue = .FALSE.       I                    nLevOutp, listId, lm,
157  #ifdef DIAGNOSTICS_MISSING_VALUE       O                    diag_mnc_bn,
158           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', 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  */  
   
159         ENDIF         ENDIF
160  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
161    
# Line 299  C---+----1----+----2----+----3----+----4 Line 166  C---+----1----+----2----+----3----+----4
166          gcode = gdiag(ndId)(1:10)          gcode = gdiag(ndId)(1:10)
167          mate = 0          mate = 0
168          mVec = 0          mVec = 0
169            mDbl = 0
170          IF ( gcode(5:5).EQ.'C' ) THEN          IF ( gcode(5:5).EQ.'C' ) THEN
171  C-      Check for Mate of a Counter Diagnostic  C-      Check for Mate of a Counter Diagnostic
172             mate = hdiag(ndId)             mate = hdiag(ndId)
173            ELSEIF ( gcode(5:5).EQ.'P' ) THEN
174    C-      Also load the mate (if stored) for Post-Processing
175               nn = ndId
176               DO WHILE ( gdiag(nn)(5:5).EQ.'P' )
177                 nn = hdiag(nn)
178               ENDDO
179               IF ( mdiag(md,listId).NE.0 ) mDbl = hdiag(nn)
180          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
181  C-      Check for Mate of a Vector Diagnostic  C-      Check for Mate of a Vector Diagnostic
182             mVec = hdiag(ndId)             mVec = hdiag(ndId)
183          ENDIF          ENDIF
184          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
185  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  
186    
187            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
188            im = mdiag(md,listId)            im = mdiag(md,listId)
189            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
190              IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
191            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
192    
193            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ndiag(ip,1,1).EQ.0 ) THEN
# Line 351  C-        Empty diagnostics case : Line 221  C-        Empty diagnostics case :
221              _END_MASTER( myThid )              _END_MASTER( myThid )
222              DO bj = myByLo(myThid), myByHi(myThid)              DO bj = myByLo(myThid), myByHi(myThid)
223                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
224                  DO k = 1,nlevels(listId)                  DO k = 1,nLevOutp
225                    DO j = 1-OLy,sNy+OLy                    DO j = 1-OLy,sNy+OLy
226                      DO i = 1-OLx,sNx+OLx                      DO i = 1-OLx,sNx+OLx
227                        qtmp1(i,j,k,bi,bj) = 0. _d 0                        qtmp1(i,j,k,bi,bj) = 0. _d 0
# Line 364  C-        Empty diagnostics case : Line 234  C-        Empty diagnostics case :
234            ELSE            ELSE
235  C-        diagnostics is not empty :  C-        diagnostics is not empty :
236    
237              IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
238                WRITE(ioUnit,'(A,I6,3A,I8,2A)')                IF ( gcode(5:5).EQ.'P' ) THEN
239                   WRITE(ioUnit,'(A,I6,7A,I8,2A)')
240         &         ' Post-Processing Diag # ', ndId, '  ', cdiag(ndId),
241         &         '   Parms: ',gdiag(ndId)
242                   IF ( mDbl.EQ.0 ) THEN
243                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
244         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
245                   ELSE
246                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
247         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
248         &          ' and diag: ',
249         &            cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
250                   ENDIF
251                  ELSE
252                   WRITE(ioUnit,'(A,I6,3A,I8,2A)')
253       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
254       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
255                  ENDIF
256                IF ( mate.GT.0 ) THEN                IF ( mate.GT.0 ) THEN
257                 WRITE(ioUnit,'(3A,I6,2A)')                 WRITE(ioUnit,'(3A,I6,2A)')
258       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
# Line 387  C-        diagnostics is not empty : Line 272  C-        diagnostics is not empty :
272                ENDIF                ENDIF
273              ENDIF              ENDIF
274    
275              IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
276  C-       get all the levels (for vertical interpolation)  C-       get only selected levels:
277                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
278                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
279                  DO k = 1,kdiag(ndId)                  DO k = 1,nlevels(listId)
280                    tmpLev = k                    kLev = NINT(levs(k,listId))
281                    CALL GETDIAG(                    CALL DIAGNOSTICS_GET_DIAG(
282       I                         tmpLev,undef,       I                         kLev, undefRL,
283       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
284       I                         ndId,mate,ip,im,bi,bj,myThid)       I                         ndId, mate, ip, im, bi, bj, myThid )
285                  ENDDO                  ENDDO
286                 ENDDO                 ENDDO
287                ENDDO                ENDDO
288                  IF ( mDbl.GT.0 ) THEN
289                   DO bj = myByLo(myThid), myByHi(myThid)
290                    DO bi = myBxLo(myThid), myBxHi(myThid)
291                     DO k = 1,nlevels(listId)
292                      kLev = NINT(levs(k,listId))
293                      CALL DIAGNOSTICS_GET_DIAG(
294         I                         kLev, undefRL,
295         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
296         I                         mDbl, 0, im, 0, bi, bj, myThid )
297                     ENDDO
298                    ENDDO
299                   ENDDO
300                  ENDIF
301              ELSE              ELSE
302  C-       get only selected levels:  C-       get all the levels (for vertical post-processing)
303                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
304                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
305                  DO k = 1,nlevels(listId)                    CALL DIAGNOSTICS_GET_DIAG(
306                    CALL GETDIAG(       I                         0, undefRL,
307       I                         levs(k,listId),undef,       O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
308       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  
309                 ENDDO                 ENDDO
310                ENDDO                ENDDO
311                  IF ( mDbl.GT.0 ) THEN
312                   DO bj = myByLo(myThid), myByHi(myThid)
313                    DO bi = myBxLo(myThid), myBxHi(myThid)
314                     DO k = 1,nlevels(listId)
315                      CALL DIAGNOSTICS_GET_DIAG(
316         I                         0, 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              ENDIF              ENDIF
324    
 C-        end of empty diag / not empty block  
           ENDIF  
   
325  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
326  C         Check to see if we need to interpolate before output  C--     Apply specific post-processing (e.g., interpolate) before output
327  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
328            IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
329  C-        Do vertical interpolation:  C-          Do vertical interpolation:
330             IF ( fluidIsAir ) THEN               IF ( fluidIsAir ) THEN
331  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);
332              CALL DIAGNOSTICS_INTERP_VERT(                CALL DIAGNOSTICS_INTERP_VERT(
333       I                     listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
334       U                     qtmp1,       U                         qtmp1, qtmp2,
335       I                     undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
336             ELSE               ELSE
337               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
338       &         'INTERP_VERT not allowed in this config'       &           'INTERP_VERT not allowed in this config'
339               CALL PRINT_ERROR( msgBuf , myThid )                 CALL PRINT_ERROR( msgBuf , myThid )
340               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
341             ENDIF               ENDIF
342                ENDIF
343                IF ( fflags(listId)(2:2).EQ.'I' ) THEN
344    C-          Integrate vertically: for now, output field has just 1 level:
345                  CALL DIAGNOSTICS_SUM_LEVELS(
346         I                         listId, md, ndId, ip, im, lm,
347         U                         qtmp1,
348         I                         undefRL, myTime, myIter, myThid )
349                ENDIF
350                IF ( gcode(5:5).EQ.'P' ) THEN
351    C-          Do Post-Processing:
352                 IF ( flds(md,listId).EQ.'PhiVEL  '
353    c    &       .OR. flds(md,listId).EQ.'PsiVEL  '
354         &          ) THEN
355                  CALL DIAGNOSTICS_CALC_PHIVEL(
356         I                         listId, md, ndId, ip, im, lm,
357         U                         qtmp1, qtmp2,
358         I                         myTime, myIter, myThid )
359                 ELSE
360                   WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
361         &           'unknown Processing method for diag="',cdiag(ndId),'"'
362                   CALL PRINT_ERROR( msgBuf , myThid )
363                   STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
364                 ENDIF
365                ENDIF
366    
367    C--     End of empty diag / not-empty diag block
368            ENDIF            ENDIF
369    
370  C--    Ready to write field "md", element "lm" in averageCycle(listId)  C--     Ready to write field "md", element "lm" in averageCycle(listId)
371    
372  C-        write to binary file, using MDSIO pkg:  C-        write to binary file, using MDSIO pkg:
373            IF ( diag_mdsio ) THEN            IF ( diag_mdsio ) THEN
374              nRec = lm + (md-1)*averageCycle(listId)  c           nRec = lm + (md-1)*averageCycle(listId)
375                nRec = md + (lm-1)*nfields(listId)
376  C           default precision for output files  C           default precision for output files
377              prec = writeBinaryPrec              prec = writeBinaryPrec
378  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 449  C           fFlag(1)=R(or D): force it t Line 381  C           fFlag(1)=R(or D): force it t
381  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
382              CALL WRITE_REC_LEV_RL(              CALL WRITE_REC_LEV_RL(
383       I                            fn, prec,       I                            fn, prec,
384       I                            NrMax, 1, nlevels(listId),       I                            NrMax, 1, nLevOutp,
385       I                            qtmp1, -nRec, myIter, myThid )       I                            qtmp1, -nRec, myIter, myThid )
386            ENDIF            ENDIF
387    
388  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
389            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
390                CALL DIAGNOSTICS_MNC_OUT(
391              _BEGIN_MASTER( myThid )       I                       NrMax, nLevOutp, listId, ndId, mate,
392         I                       diag_mnc_bn, qtmp1,
393              DO ii = 1,CW_DIMS       I                       undefRL, myTime, myIter, myThid )
               d_cw_name(1:NLEN) = dn_blnk(1:NLEN)  
               dn(ii)(1:NLEN) = dn_blnk(1:NLEN)  
             ENDDO  
   
 C           Note that the "d_cw_name" variable is a hack that hides a  
 C           subtlety within MNC.  Basically, each MNC-wrapped file is  
 C           caching its own concept of what each "grid name" (that is, a  
 C           dimension group name) means.  So one cannot re-use the same  
 C           "grid" name for different collections of dimensions within a  
 C           given file.  By appending the "ndId" values to each name, we  
 C           guarantee uniqueness within each MNC-produced file.  
             WRITE(d_cw_name,'(a,i6.6)') 'd_cw_',ndId  
   
 C           XY dimensions  
             dim(1)       = sNx + 2*OLx  
             dim(2)       = sNy + 2*OLy  
             ib(1)        = OLx + 1  
             ib(2)        = OLy + 1  
             IF (gdiag(ndId)(2:2) .EQ. 'M') THEN  
               dn(1)(1:2) = 'X'  
               ie(1)      = OLx + sNx  
               dn(2)(1:2) = 'Y'  
               ie(2)      = OLy + sNy  
             ELSEIF (gdiag(ndId)(2:2) .EQ. 'U') THEN  
               dn(1)(1:3) = 'Xp1'  
               ie(1)      = OLx + sNx + 1  
               dn(2)(1:2) = 'Y'  
               ie(2)      = OLy + sNy  
             ELSEIF (gdiag(ndId)(2:2) .EQ. 'V') THEN  
               dn(1)(1:2) = 'X'  
               ie(1)      = OLx + sNx  
               dn(2)(1:3) = 'Yp1'  
               ie(2)      = OLy + sNy + 1  
             ELSEIF (gdiag(ndId)(2:2) .EQ. 'Z') THEN  
               dn(1)(1:3) = 'Xp1'  
               ie(1)      = OLx + sNx + 1  
               dn(2)(1:3) = 'Yp1'  
               ie(2)      = OLy + sNy + 1  
             ENDIF  
   
 C           Z is special since it varies  
             WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevels(listId)  
             ENDIF  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zld', nlevels(listId)  
             ENDIF  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)  
             ENDIF  
             dim(3) = NrMax  
             ib(3)  = 1  
             ie(3)  = nlevels(listId)  
   
 C           Time dimension  
             dn(4)(1:1) = 'T'  
             dim(4) = -1  
             ib(4)  = 1  
             ie(4)  = 1  
   
             CALL MNC_CW_ADD_GNAME(d_cw_name, 4,  
      &             dim, dn, ib, ie, myThid)  
             CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,  
      &             4,5, myThid)  
             CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',  
      &             tdiag(ndId),myThid)  
             CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',  
      &             udiag(ndId),myThid)  
   
 C     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 ( 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 )  
   
394            ENDIF            ENDIF
395  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
396    
 C--      end loop on lm (or ll if ALLOW_MNC) counter  
          ENDDO  
397  C--     end of Processing Fld # md  C--     end of Processing Fld # md
398          ENDIF          ENDIF
399         ENDDO         ENDDO
400    
401  #ifdef ALLOW_MNC  C--   end loop on lm counter (= averagePeriod)
 C--   end loop on jj counter  
402        ENDDO        ENDDO
 #endif  
403    
404  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
405        IF (diag_mdsio) THEN        IF (diag_mdsio) THEN
406  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
407  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
408  C     meta files but with more informations in it.  C     meta files but with more informations in it.
409              glf = globalFiles              glf = globalFiles
410              nRec = nfields(listId)*averageCycle(listId)              nRec = averageCycle(listId)*nfields(listId)
411              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
412       &              0, 0, nlevels(listId), ' ',       &              0, 0, nLevOutp, ' ',
413       &              nfields(listId), flds(1,listId), nTimRec, timeRec,       &              nfields(listId), flds(1,listId), nTimRec, timeRec,
414       &              nRec, myIter, myThid)       &              nRec, myIter, myThid)
415        ENDIF        ENDIF
Legend:
Removed from v.1.44  
changed lines
  Added in v.1.57
  ViewVC Help
Powered by ViewVC 1.1.22