/[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.61 by jmc, Wed Feb 6 21:25:26 2013 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 41  CEOP Line 41  CEOP
41  C     !FUNCTIONS:  C     !FUNCTIONS:
42        INTEGER ILNBLNK        INTEGER ILNBLNK
43        EXTERNAL ILNBLNK        EXTERNAL ILNBLNK
 #ifdef ALLOW_FIZHI  
       _RL   getcon  
       EXTERNAL getcon  
 #endif  
44    
45  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
46  C     i,j,k :: loop indices  C     i,j,k :: loop indices
# Line 52  C     bi,bj :: tile indices Line 48  C     bi,bj :: tile indices
48  C     lm    :: loop index (averageCycle)  C     lm    :: loop index (averageCycle)
49  C     md    :: field number in the list "listId".  C     md    :: field number in the list "listId".
50  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)  
51  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
52  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
53    C     mate  :: counter mate Id number (in available diagnostics list)
54    C     mDbl  :: processing mate Id number (in case processing requires 2 diags)
55    C     mVec  :: vector mate Id number
56    C     ppFld :: post-processed diag or not (=0): =1 stored in qtmp1 ; =2 in qtmp2
57    C   isComputed :: previous post-processed diag (still available in qtmp)
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          INTEGER ppFld, isComputed
73        CHARACTER*10 gcode        CHARACTER*10 gcode
74        _RL undef        _RL undefRL
75        _RL tmpLev        INTEGER nLevOutp, kLev
       INTEGER ilen  
       INTEGER nLevOutp  
76    
77          INTEGER iLen
78        INTEGER ioUnit        INTEGER ioUnit
79        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
80        CHARACTER*(MAX_LEN_MBUF) suff        CHARACTER*(MAX_LEN_MBUF) suff
81        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
82        INTEGER prec, nRec, nTimRec        INTEGER prec, nRec, nTimRec
83        _RL     timeRec(2)        _RL     timeRec(2)
84          _RL     tmpLoc
85  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
86        LOGICAL glf        LOGICAL glf
87  #endif  #endif
88  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
       INTEGER ll, llMx, jj, jjMx  
       INTEGER ii, klev  
89        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)  
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 = misValFlt(listId)
97  #ifdef ALLOW_FIZHI  
 c     IF ( useFIZHI ) undef = getcon('UNDEF')  
       undef = getcon('UNDEF')  
 #endif  
98        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
99        ilen = ILNBLNK(fnames(listId))        iLen = ILNBLNK(fnames(listId))
100        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
101  C-    for now, if integrate vertically, output field has just 1 level:  C-    for now, if integrate vertically, output field has just 1 level:
102        nLevOutp = nlevels(listId)        nLevOutp = nlevels(listId)
103        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 123  C     a) find the time of the previous m
123          timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)          timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
124          i = INT( timeRec(1) )          i = INT( timeRec(1) )
125          IF ( timeRec(1).LT.0. ) THEN          IF ( timeRec(1).LT.0. ) THEN
126            tmpLev = FLOAT(i)            tmpLoc = FLOAT(i)
127            IF ( timeRec(1).NE.tmpLev ) i = i - 1            IF ( timeRec(1).NE.tmpLoc ) i = i - 1
128          ENDIF          ENDIF
129          timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)          timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
130  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 134  C     b) round off to nearest multiple o
134          timeRec(1) = (timeRec(1)-baseTime)/deltaTClock          timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
135          i = NINT( timeRec(1) )          i = NINT( timeRec(1) )
136  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
137          tmpLev = FLOAT(i) - 0.5 _d 0          tmpLoc = FLOAT(i) - 0.5 _d 0
138          IF ( timeRec(1).EQ.tmpLev ) i = i - 1          IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
139          timeRec(1) = baseTime + deltaTClock*FLOAT(i)          timeRec(1) = baseTime + deltaTClock*FLOAT(i)
140  c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock  c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
141        ENDIF        ENDIF
# Line 163  c     DO i=1,nTimRec Line 144  c     DO i=1,nTimRec
144  c       timeRec(i) = timeRec(i)/deltaTClock  c       timeRec(i) = timeRec(i)/deltaTClock
145  c     ENDDO  c     ENDDO
146    
147  #ifdef ALLOW_MNC  C--   Place the loop on lm (= averagePeriod) outside the loop on md (= field):
148  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  
149    
150    #ifdef ALLOW_MNC
151         IF (useMNC .AND. diag_mnc) THEN         IF (useMNC .AND. diag_mnc) THEN
152  C     Handle missing value attribute (land points)           CALL DIAGNOSTICS_MNC_SET(
153           useMissingValue = .FALSE.       I                    nLevOutp, listId, lm,
154  #ifdef DIAGNOSTICS_MISSING_VALUE       O                    diag_mnc_bn,
155           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  */  
   
156         ENDIF         ENDIF
157  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
158    
159  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
160    
161           isComputed = 0
162         DO md = 1,nfields(listId)         DO md = 1,nfields(listId)
163          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
164          gcode = gdiag(ndId)(1:10)          gcode = gdiag(ndId)(1:10)
165          mate = 0          mate = 0
166          mVec = 0          mVec = 0
167            mDbl = 0
168            ppFld = 0
169          IF ( gcode(5:5).EQ.'C' ) THEN          IF ( gcode(5:5).EQ.'C' ) THEN
170  C-      Check for Mate of a Counter Diagnostic  C-      Check for Mate of a Counter Diagnostic
171             mate = hdiag(ndId)             mate = hdiag(ndId)
172            ELSEIF ( gcode(5:5).EQ.'P' ) THEN
173               ppFld = 1
174               IF ( gdiag(hdiag(ndId))(5:5).EQ.'P' ) ppFld = 2
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    c          write(0,*) ppFld,' ndId=', ndId, nn, mDbl, isComputed
182          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
183  C-      Check for Mate of a Vector Diagnostic  C-      Check for Mate of a Vector Diagnostic
184             mVec = hdiag(ndId)             mVec = hdiag(ndId)
185          ENDIF          ENDIF
186          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
187  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  
188    
189            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
190            im = mdiag(md,listId)            im = mdiag(md,listId)
191            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
192              IF (mDbl.GT.0) im = im + kdiag(mDbl)*(lm-1)
193            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)            IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
194    
195            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ppFld.EQ.2 .AND. isComputed.EQ.hdiag(ndId) ) THEN
196    C-        Post-Processed diag from an other Post-Processed diag -and-
197    C         both of them have just been calculated and are still stored in qtmp:
198    C         => skip computation and just write qtmp2
199                IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
200                   WRITE(ioUnit,'(A,I6,3A,I6)')
201         &         '  get Post-Proc. Diag # ', ndId, '  ', cdiag(ndId),
202         &         ' from previous computation of Diag # ', isComputed
203                ENDIF
204                isComputed = 0
205              ELSEIF ( ndiag(ip,1,1).EQ.0 ) THEN
206  C-        Empty diagnostics case :  C-        Empty diagnostics case :
207                isComputed = 0
208    
209              _BEGIN_MASTER( myThid )              _BEGIN_MASTER( myThid )
210              WRITE(msgBuf,'(A,I10)')              WRITE(msgBuf,'(A,I10)')
# Line 377  C-        Empty diagnostics case : Line 246  C-        Empty diagnostics case :
246    
247            ELSE            ELSE
248  C-        diagnostics is not empty :  C-        diagnostics is not empty :
249                isComputed = 0
250    
251              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
252                WRITE(ioUnit,'(A,I6,3A,I8,2A)')                IF ( ppFld.GE.1 ) THEN
253                   WRITE(ioUnit,'(A,I6,7A,I8,2A)')
254         &         ' Post-Processing Diag # ', ndId, '  ', cdiag(ndId),
255         &         '   Parms: ',gdiag(ndId)
256                   IF ( mDbl.EQ.0 ) THEN
257                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
258         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1)
259                   ELSE
260                    WRITE(ioUnit,'(2(3A,I6,A,I8))') '   from diag: ',
261         &            cdiag(nn), ' (#', nn, ') Cnt=', ndiag(ip,1,1),
262         &          ' and diag: ',
263         &            cdiag(mDbl),' (#',mDbl,') Cnt=',ndiag(im,1,1)
264                   ENDIF
265                  ELSE
266                   WRITE(ioUnit,'(A,I6,3A,I8,2A)')
267       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
268       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
269                  ENDIF
270                IF ( mate.GT.0 ) THEN                IF ( mate.GT.0 ) THEN
271                 WRITE(ioUnit,'(3A,I6,2A)')                 WRITE(ioUnit,'(3A,I6,2A)')
272       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
# Line 401  C-        diagnostics is not empty : Line 286  C-        diagnostics is not empty :
286                ENDIF                ENDIF
287              ENDIF              ENDIF
288    
289              IF ( fflags(listId)(2:2).NE.' ' ) THEN              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
290  C-       get all the levels (for vertical post-processing)  C-       get only selected levels:
291                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
292                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
293                  DO k = 1,kdiag(ndId)                  DO k = 1,nlevels(listId)
294                    tmpLev = k                    kLev = NINT(levs(k,listId))
295                    CALL GETDIAG(                    CALL DIAGNOSTICS_GET_DIAG(
296       I                         tmpLev,undef,       I                         kLev, undefRL,
297       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),       O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
298       I                         ndId,mate,ip,im,bi,bj,myThid)       I                         ndId, mate, ip, im, bi, bj, myThid )
299                  ENDDO                  ENDDO
300                 ENDDO                 ENDDO
301                ENDDO                ENDDO
302                  IF ( mDbl.GT.0 ) THEN
303                   DO bj = myByLo(myThid), myByHi(myThid)
304                    DO bi = myBxLo(myThid), myBxHi(myThid)
305                     DO k = 1,nlevels(listId)
306                      kLev = NINT(levs(k,listId))
307                      CALL DIAGNOSTICS_GET_DIAG(
308         I                         kLev, undefRL,
309         O                         qtmp2(1-OLx,1-OLy,k,bi,bj),
310         I                         mDbl, 0, im, 0, bi, bj, myThid )
311                     ENDDO
312                    ENDDO
313                   ENDDO
314                  ENDIF
315              ELSE              ELSE
316  C-       get only selected levels:  C-       get all the levels (for vertical post-processing)
317                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
318                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
319                  DO k = 1,nlevels(listId)                    CALL DIAGNOSTICS_GET_DIAG(
320                    CALL GETDIAG(       I                         0, undefRL,
321       I                         levs(k,listId),undef,       O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
322       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  
323                 ENDDO                 ENDDO
324                ENDDO                ENDDO
325                  IF ( mDbl.GT.0 ) THEN
326                   DO bj = myByLo(myThid), myByHi(myThid)
327                    DO bi = myBxLo(myThid), myBxHi(myThid)
328                      CALL DIAGNOSTICS_GET_DIAG(
329         I                         0, undefRL,
330         O                         qtmp2(1-OLx,1-OLy,1,bi,bj),
331         I                         mDbl, 0, im, 0, bi, bj, myThid )
332                    ENDDO
333                   ENDDO
334                  ENDIF
335              ENDIF              ENDIF
336    
337  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
# Line 437  C-          Do vertical interpolation: Line 343  C-          Do vertical interpolation:
343  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);
344                CALL DIAGNOSTICS_INTERP_VERT(                CALL DIAGNOSTICS_INTERP_VERT(
345       I                         listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
346       U                         qtmp1,       U                         qtmp1, qtmp2,
347       I                         undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
348               ELSE               ELSE
349                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
350       &           'INTERP_VERT not allowed in this config'       &           'INTERP_VERT not allowed in this config'
# Line 451  C-          Integrate vertically: for no Line 357  C-          Integrate vertically: for no
357                CALL DIAGNOSTICS_SUM_LEVELS(                CALL DIAGNOSTICS_SUM_LEVELS(
358       I                         listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
359       U                         qtmp1,       U                         qtmp1,
360       I                         undef, myTime, myIter, myThid )       I                         undefRL, myTime, myIter, myThid )
361                ENDIF
362                IF ( ppFld.GE.1 ) THEN
363    C-          Do Post-Processing:
364                 IF ( flds(md,listId).EQ.'PhiVEL  '
365         &       .OR. flds(md,listId).EQ.'PsiVEL  '
366         &          ) THEN
367                  CALL DIAGNOSTICS_CALC_PHIVEL(
368         I                         listId, md, ndId, ip, im, lm,
369         I                         NrMax,
370         U                         qtmp1, qtmp2,
371         I                         myTime, myIter, myThid )
372                  isComputed = ndId
373                 ELSE
374                   WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
375         &           'unknown Processing method for diag="',cdiag(ndId),'"'
376                   CALL PRINT_ERROR( msgBuf , myThid )
377                   STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
378                 ENDIF
379              ENDIF              ENDIF
380    
381  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 385  C--     Ready to write field "md", eleme
385    
386  C-        write to binary file, using MDSIO pkg:  C-        write to binary file, using MDSIO pkg:
387            IF ( diag_mdsio ) THEN            IF ( diag_mdsio ) THEN
388              nRec = lm + (md-1)*averageCycle(listId)  c          nRec = lm + (md-1)*averageCycle(listId)
389  C           default precision for output files             nRec = md + (lm-1)*nfields(listId)
390              prec = writeBinaryPrec  C         default precision for output files
391  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision             prec = writeBinaryPrec
392              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32  C         fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
393              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64             IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
394               IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
395  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
396               IF ( ppFld.LE.1 ) THEN
397              CALL WRITE_REC_LEV_RL(              CALL WRITE_REC_LEV_RL(
398       I                            fn, prec,       I                            fn, prec,
399       I                            NrMax, 1, nLevOutp,       I                            NrMax, 1, nLevOutp,
400       I                            qtmp1, -nRec, myIter, myThid )       I                            qtmp1, -nRec, myIter, myThid )
401               ELSE
402                CALL WRITE_REC_LEV_RL(
403         I                            fn, prec,
404         I                            NrMax, 1, nLevOutp,
405         I                            qtmp2, -nRec, myIter, myThid )
406               ENDIF
407            ENDIF            ENDIF
408    
409  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
410            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
411               IF ( ppFld.LE.1 ) THEN
412              _BEGIN_MASTER( myThid )              CALL DIAGNOSTICS_MNC_OUT(
413         I                       NrMax, nLevOutp, listId, ndId, mate,
414              DO ii = 1,CW_DIMS       I                       diag_mnc_bn, qtmp1,
415                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)       I                       undefRL, myTime, myIter, myThid )
416                dn(ii)(1:NLEN) = dn_blnk(1:NLEN)             ELSE
417              ENDDO              CALL DIAGNOSTICS_MNC_OUT(
418         I                       NrMax, nLevOutp, listId, ndId, mate,
419  C           Note that the "d_cw_name" variable is a hack that hides a       I                       diag_mnc_bn, qtmp2,
420  C           subtlety within MNC.  Basically, each MNC-wrapped file is       I                       undefRL, myTime, myIter, myThid )
421  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', 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 )  
   
422            ENDIF            ENDIF
423  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
424    
 C--      end loop on lm (or ll if ALLOW_MNC) counter  
          ENDDO  
425  C--     end of Processing Fld # md  C--     end of Processing Fld # md
426          ENDIF          ENDIF
427         ENDDO         ENDDO
428    
429  #ifdef ALLOW_MNC  C--   end loop on lm counter (= averagePeriod)
 C--   end loop on jj counter  
430        ENDDO        ENDDO
 #endif  
431    
432  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
433        IF (diag_mdsio) THEN        IF (diag_mdsio) THEN
# Line 625  C-    Note: temporary: since it is a pai Line 435  C-    Note: temporary: since it is a pai
435  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
436  C     meta files but with more informations in it.  C     meta files but with more informations in it.
437              glf = globalFiles              glf = globalFiles
438              nRec = nfields(listId)*averageCycle(listId)              nRec = averageCycle(listId)*nfields(listId)
439              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
440       &              0, 0, nLevOutp, ' ',       &              0, 0, nLevOutp, ' ',
441       &              nfields(listId), flds(1,listId), nTimRec, timeRec,       &              nfields(listId), flds(1,listId),
442         &              nTimRec, timeRec, undefRL,
443       &              nRec, myIter, myThid)       &              nRec, myIter, myThid)
444        ENDIF        ENDIF
445  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
Legend:
Removed from v.1.49  
changed lines
  Added in v.1.61
  ViewVC Help
Powered by ViewVC 1.1.22