/[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.19 by molod, Mon Jul 18 18:35:19 2005 UTC revision 1.53 by jmc, Tue Jun 14 00:18:37 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 26  C     !USES: Line 26  C     !USES:
26  #include "DIAGNOSTICS_SIZE.h"  #include "DIAGNOSTICS_SIZE.h"
27  #include "DIAGNOSTICS.h"  #include "DIAGNOSTICS.h"
28    
29  #ifdef ALLOW_FIZHI        INTEGER NrMax
30  #include "fizhi_SIZE.h"        PARAMETER( NrMax = numLevels )
 #else  
       INTEGER Nrphys  
       PARAMETER (Nrphys=0)  
 #endif  
   
31    
32  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
33  C     listId  :: Diagnostics list number being written  C     listId  :: Diagnostics list number being written
# Line 43  C     myThid  :: my Thread Id number Line 38  C     myThid  :: my Thread Id number
38        INTEGER listId, myIter, myThid        INTEGER listId, myIter, myThid
39  CEOP  CEOP
40    
41    C     !FUNCTIONS:
42          INTEGER ILNBLNK
43          EXTERNAL ILNBLNK
44    #ifdef ALLOW_FIZHI
45          _RL   getcon
46          EXTERNAL getcon
47    #endif
48    
49  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
50  C     i,j,k :: loop indices  C     i,j,k :: loop indices
51    C     bi,bj :: tile indices
52    C     lm    :: loop index (averageCycle)
53  C     md    :: field number in the list "listId".  C     md    :: field number in the list "listId".
54  C     ndId  :: diagnostics  Id number (in available diagnostics list)  C     ndId  :: diagnostics  Id number (in available diagnostics list)
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        INTEGER i, j, k  C     nLevOutp :: number of levels to write in output file
59    C
60    C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
61    C     qtmp1 :: temporary array; used to store a copy of diag. output field.
62    C     qtmp2 :: temporary array; used to store a copy of a 2nd diag. field.
63    C-  Note: local common block no longer needed.
64    c     COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
65          _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
69        INTEGER bi, bj        INTEGER bi, bj
70        INTEGER md, ndId, ip, im        INTEGER md, ndId, ip, im
71        INTEGER mate, mVec        INTEGER mate, mVec
72        CHARACTER*8 parms1        CHARACTER*10 gcode
73        CHARACTER*3 mate_index        _RL undefRL
74        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)        INTEGER nFilled
75        _RL qtmpsrf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)        INTEGER nLevOutp, kLev
       _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)  
       _RL undef, getcon  
       EXTERNAL getcon  
       INTEGER ILNBLNK  
       EXTERNAL ILNBLNK  
       INTEGER ilen  
       integer nlevsout,nplevs  
       parameter(nplevs = 16)  
       _RL plevs1(nplevs)  
       data plevs1/ 1000.0 _d 2, 925.0 _d 2, 850.0 _d 2, 700.0 _d 2,  
      .              600.0 _d 2, 500.0 _d 2, 400.0 _d 2, 300.0 _d 2,  
      .              250.0 _d 2, 200.0 _d 2, 150.0 _d 2, 100.0 _d 2,  
      .               70.0 _d 2,  50.0 _d 2,  30.0 _d 2,  20.0 _d 2/  
       _RL plevs2(nplevs)  
       data plevs2/ 1000.0 _d 2, 950.0 _d 2, 900.0 _d 2, 850.0 _d 2,  
      .              800.0 _d 2, 750.0 _d 2, 700.0 _d 2, 600.0 _d 2,  
      .              500.0 _d 2, 400.0 _d 2, 300.0 _d 2, 250.0 _d 2,  
      .              200.0 _d 2, 150.0 _d 2, 100.0 _d 2,  50.0 _d 2/  
       _RL qprs(sNx,sNy,nplevs)  
       _RL qinp(sNx,sNy,Nr+Nrphys)  
       _RL pkz(sNx,sNy,Nr+Nrphys)  
       _RL pksrf(sNx,sNy)  
       _RL p  
       INTEGER jpoint1,ipoint1  
       INTEGER jpoint2,ipoint2  
       _RL kappa  
       logical foundp  
       data foundp /.false./  
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
83          _RL     timeRec(2)
84          _RL     tmpLoc
85    #ifdef ALLOW_MDSIO
86        LOGICAL glf        LOGICAL glf
87    #endif
88  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
89        INTEGER ii        INTEGER ll, llMx, jj, jjMx
90        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
91        CHARACTER*(5) ctmp        LOGICAL useMissingValue
92        INTEGER CW_DIMS, NLEN        REAL*8 misValLoc
       PARAMETER ( CW_DIMS = 10 )  
       PARAMETER ( NLEN    = 80 )  
       INTEGER dim(CW_DIMS), ib(CW_DIMS), ie(CW_DIMS)  
       CHARACTER*(NLEN) dn(CW_DIMS)  
       CHARACTER*(NLEN) d_cw_name  
       CHARACTER*(NLEN) dn_blnk  
       _RS ztmp(Nr+Nrphys)  
93  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
94    
95  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
96    
97    C---  set file properties
98        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
99        undef = getcon('UNDEF')        undefRL = UNSET_RL
100        kappa = getcon('KAPPA')  #ifdef ALLOW_FIZHI
101        glf = globalFiles        IF ( useFIZHI ) undefRL = getcon('UNDEF')
102    #endif
103        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
104        ilen = ILNBLNK(fnames(listId))        iLen = ILNBLNK(fnames(listId))
105        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:iLen),'.',suff(1:10)
106    C-    for now, if integrate vertically, output field has just 1 level:
107          nLevOutp = nlevels(listId)
108          IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
109    
110    C--   Set time information:
111          IF ( freq(listId).LT.0. ) THEN
112    C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
113            nTimRec = 1
114            timeRec(1) = myTime
115          ELSE
116    C-    Time-average: store the 2 edges of the time-averaging interval.
117    C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
118    C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
119            nTimRec = 2
120    
121    C-    end of time-averaging interval:
122            timeRec(2) = myTime
123    
124    C-    begining of time-averaging interval:
125    c       timeRec(1) = myTime - freq(listId)
126    C     a) find the time of the previous multiple of output freq:
127            timeRec(1) = myTime-deltaTClock*0.5 _d 0
128            timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
129            i = INT( timeRec(1) )
130            IF ( timeRec(1).LT.0. ) THEN
131              tmpLoc = FLOAT(i)
132              IF ( timeRec(1).NE.tmpLoc ) i = i - 1
133            ENDIF
134            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
135    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
136            timeRec(1) = MAX( timeRec(1), startTime )
137    
138    C     b) round off to nearest multiple of time-step:
139            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
140            i = NINT( timeRec(1) )
141    C     if just half way, NINT will return the next time-step: correct this
142            tmpLoc = FLOAT(i) - 0.5 _d 0
143            IF ( timeRec(1).EQ.tmpLoc ) i = i - 1
144            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
145    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
146          ENDIF
147    C--   Convert time to iteration number (debug)
148    c     DO i=1,nTimRec
149    c       timeRec(i) = timeRec(i)/deltaTClock
150    c     ENDDO
151    
152  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
153    C-- this is a trick to reverse the order of the loops on md (= field)
154    C   and lm (= averagePeriod): binary output: lm loop inside md loop ;
155    C                                 mnc ouput: md loop inside lm loop.
156        IF (useMNC .AND. diag_mnc) THEN        IF (useMNC .AND. diag_mnc) THEN
157          DO i = 1,MAX_LEN_FNAM          jjMx = averageCycle(listId)
158            diag_mnc_bn(i:i) = ' '          llMx = 1
159          ENDDO        ELSE
160          DO i = 1,NLEN          jjMx = 1
161            dn_blnk(i:i) = ' '          llMx = averageCycle(listId)
         ENDDO  
         WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)  
   
 C       Update the record dimension by writing the iteration number  
         CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)  
         CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',myTime,myThid)  
         CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)  
   
         dn(1)(1:NLEN) = dn_blnk(1:NLEN)  
         WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)  
         dim(1) = nlevels(listId)  
         ib(1)  = 1  
         ie(1)  = nlevels(listId)  
   
         CALL MNC_CW_ADD_GNAME('diag_levels', 1,  
      &       dim, dn, ib, ie, myThid)  
         CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',  
      &       0,0, myThid)  
         CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',  
      &       'Idicies of vertical levels within the source arrays',  
      &       myThid)  
           
         CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,  
      &       'diag_levels', levs(1,listId), myThid)  
   
         CALL MNC_CW_DEL_VNAME('diag_levels', myThid)  
         CALL MNC_CW_DEL_GNAME('diag_levels', myThid)  
   
 #ifdef DIAG_MNC_COORD_NEEDSWORK  
 C       This part has been placed in an #ifdef because, as its currently  
 C       written, it will only work with variables defined on a dynamics  
 C       grid.  As we start using diagnostics for physics grids, ice  
 C       levels, land levels, etc. the different vertical coordinate  
 C       dimensions will have to be taken into account.  
   
 C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx  
         ctmp(1:5) = 'mul  '  
         DO i = 1,3  
           dn(1)(1:NLEN) = dn_blnk(1:NLEN)  
           WRITE(dn(1),'(3a,i6.6)') 'Z',ctmp(i:i),'d',nlevels(listId)  
           CALL MNC_CW_ADD_GNAME(dn(1), 1, dim, dn, ib, ie, myThid)  
           CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)  
   
 C         The following three ztmp() loops should eventually be modified  
 C         to reflect the fractional nature of levs(j,l) -- they should  
 C         do something like:  
 C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))  
 C                      + ( rC(INT(FLOOR(levs(j,l))))  
 C                          + rC(INT(CEIL(levs(j,l)))) )  
 C                        / ( levs(j,l) - FLOOR(levs(j,l)) )  
 C         for averaged levels.  
           IF (i .EQ. 1) THEN  
             DO j = 1,nlevels(listId)  
               ztmp(j) = rC(NINT(levs(j,listId)))  
             ENDDO  
             CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',  
      &           'Dimensional coordinate value at the mid point',  
      &           myThid)  
           ELSEIF (i .EQ. 2) THEN  
             DO j = 1,nlevels(listId)  
               ztmp(j) = rF(NINT(levs(j,listId)) + 1)  
             ENDDO  
             CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',  
      &           'Dimensional coordinate value at the upper point',  
      &           myThid)  
           ELSEIF (i .EQ. 3) THEN  
             DO j = 1,nlevels(listId)  
               ztmp(j) = rF(NINT(levs(j,listId)))  
             ENDDO  
             CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',  
      &           'Dimensional coordinate value at the lower point',  
      &           myThid)  
           ENDIF  
           CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)  
           CALL MNC_CW_DEL_VNAME(dn(1), myThid)  
           CALL MNC_CW_DEL_GNAME(dn(1), myThid)  
         ENDDO  
 #endif /*  DIAG_MNC_COORD_NEEDSWORK  */  
   
162        ENDIF        ENDIF
163          DO jj=1,jjMx
164    
165           IF (useMNC .AND. diag_mnc) THEN
166             misValLoc = undefRL
167             IF ( misvalFlt(listId).NE.UNSET_RL )
168         &        misValLoc = misvalFlt(listId)
169             CALL DIAGNOSTICS_MNC_SET(
170         I                    nLevOutp, listId, jj,
171         O                    diag_mnc_bn, useMissingValue,
172         I                    misValLoc, myTime, myIter, myThid )
173           ENDIF
174  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
175    
176        DO md = 1,nfields(listId)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
177    
178           DO md = 1,nfields(listId)
179          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
180          parms1 = gdiag(ndId)(1:8)          gcode = gdiag(ndId)(1:10)
181          IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN          mate = 0
182            mVec = 0
183            IF ( gcode(5:5).EQ.'C' ) THEN
184    C-      Check for Mate of a Counter Diagnostic
185               mate = hdiag(ndId)
186            ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
187    C-      Check for Mate of a Vector Diagnostic
188               mVec = hdiag(ndId)
189            ENDIF
190            IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
191  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
192    #ifdef ALLOW_MNC
193             DO ll=1,llMx
194              lm = jj+ll-1
195    #else
196             DO lm=1,averageCycle(listId)
197    #endif
198    
199            ip = ABS(idiag(md,listId))            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
200            im = mdiag(md,listId)            im = mdiag(md,listId)
201            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
202              IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
203    
204              nFilled = ndiag(ip,1,1)
205              IF ( flds(md,listId).EQ.'PhiVEL  ' ) THEN
206                  CALL DIAGNOSTICS_CALC_PHIVEL(
207         I                         listId, md, ndId, ip, im, lm,
208         U                         nFilled, qtmp1, qtmp2,
209         I                         myTime, myIter, myThid )
210              ENDIF
211    
212    c         IF ( ndiag(ip,1,1).EQ.0 ) THEN
213              IF ( nFilled.EQ.0 ) THEN
214  C-        Empty diagnostics case :  C-        Empty diagnostics case :
215    
216              _BEGIN_MASTER( myThid )              _BEGIN_MASTER( myThid )
# Line 221  C-        Empty diagnostics case : Line 218  C-        Empty diagnostics case :
218       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
219              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
220       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
221              WRITE(msgBuf,'(A,I4,3A,I3,2A)')              WRITE(msgBuf,'(A,I6,3A,I4,2A)')
222       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
223       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
224              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
225       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
226              WRITE(msgBuf,'(A,I2,A)')              IF ( averageCycle(listId).GT.1 ) THEN
227       &       '- WARNING -   has not been filled (ndiag=',               WRITE(msgBuf,'(A,2(I3,A))')
228       &       ndiag(ip,1,1), ' )'       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
229         &                                            ndiag(ip,1,1), ' )'
230                ELSE
231                 WRITE(msgBuf,'(A,2(I3,A))')
232         &        '- WARNING -   has not been filled (ndiag=',
233         &                                            ndiag(ip,1,1), ' )'
234                ENDIF
235              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
236       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
237              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
# Line 238  C-        Empty diagnostics case : Line 241  C-        Empty diagnostics case :
241              _END_MASTER( myThid )              _END_MASTER( myThid )
242              DO bj = myByLo(myThid), myByHi(myThid)              DO bj = myByLo(myThid), myByHi(myThid)
243                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
244                  DO k = 1,nlevels(listId)                  DO k = 1,nLevOutp
245                    DO j = 1-OLy,sNy+OLy                    DO j = 1-OLy,sNy+OLy
246                      DO i = 1-OLx,sNx+OLx                      DO i = 1-OLx,sNx+OLx
247                        qtmp1(i,j,k,bi,bj) = 0. _d 0                        qtmp1(i,j,k,bi,bj) = 0. _d 0
# Line 248  C-        Empty diagnostics case : Line 251  C-        Empty diagnostics case :
251                ENDDO                ENDDO
252              ENDDO              ENDDO
253    
254            ELSE  c         ELSE
255              ELSEIF ( ndiag(ip,1,1).NE.0 ) THEN
256  C-        diagnostics is not empty :  C-        diagnostics is not empty :
257    
258              IF ( myThid.EQ.1 ) WRITE(ioUnit,'(A,I3,3A,I8,2A)')              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
259                  WRITE(ioUnit,'(A,I6,3A,I8,2A)')
260       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
261       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
262                  IF ( mate.GT.0 ) THEN
263              IF ( parms1(5:5).EQ.'C' ) THEN                 WRITE(ioUnit,'(3A,I6,2A)')
 C             Check for Mate of a Counter Diagnostic  
 C             --------------------------------------  
               mate_index = parms1(6:8)  
               READ (mate_index,'(I3)') mate  
               IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,2A)')  
264       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
265       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
266                  ELSEIF ( mVec.GT.0 ) THEN
             ELSE  
               mate = 0  
   
 C             Check for Mate of a Vector Diagnostic  
 C             -------------------------------------  
               IF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN  
                 mate_index = parms1(6:8)  
                 READ (mate_index,'(I3)') mVec  
267                  IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN                  IF ( im.GT.0 .AND. ndiag(MAX(1,im),1,1).GT.0 ) THEN
268                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
269       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
270       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
271       &             ' exists '       &             ' exists '
272                  ELSE                  ELSE
273                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
274       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
275       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
276       &             ' not enabled'       &             ' not enabled'
# Line 286  C             -------------------------- Line 278  C             --------------------------
278                ENDIF                ENDIF
279              ENDIF              ENDIF
280    
281              DO bj = myByLo(myThid), myByHi(myThid)              IF ( fflags(listId)(2:2).EQ.' ' ) THEN
282               DO bi = myBxLo(myThid), myBxHi(myThid)  C-       get only selected levels:
283                DO k = 1,nlevels(listId)                DO bj = myByLo(myThid), myByHi(myThid)
284                  CALL GETDIAG(                 DO bi = myBxLo(myThid), myBxHi(myThid)
285       I                       levs(k,listId),undef,                  DO k = 1,nlevels(listId)
286       O                       qtmp1(1-OLx,1-OLy,k,bi,bj),                    kLev = NINT(levs(k,listId))
287       I                       ndId,mate,ip,im,bi,bj,myThid)                    CALL DIAGNOSTICS_GET_DIAG(
288         I                         kLev, undefRL,
289         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
290         I                         ndId,mate,ip,im,bi,bj,myThid)
291                    ENDDO
292                   ENDDO
293                ENDDO                ENDDO
294               ENDDO              ELSE
295              ENDDO  C-       get all the levels (for vertical post-processing)
296                  DO bj = myByLo(myThid), myByHi(myThid)
297  C-        end of empty diag / not empty block                 DO bi = myBxLo(myThid), myBxHi(myThid)
298            ENDIF                    CALL DIAGNOSTICS_GET_DIAG(
299         I                         0, undefRL,
300            nlevsout = nlevels(listId)       O                         qtmp1(1-OLx,1-OLy,1,bi,bj),
301         I                         ndId,mate,ip,im,bi,bj,myThid)
302                   ENDDO
303                  ENDDO
304                ENDIF
305    
306  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
307  C             Check to see if we need to interpolate before output  C--     Apply specific post-processing (e.g., interpolate) before output
308  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
309  C  (we are still inside field exist if sequence and field do loop)              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
310  C  C-          Do vertical interpolation:
311                 IF ( fluidIsAir ) THEN
312           if(fflags(listId)(2:2).eq.'P') then  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
313                  CALL DIAGNOSTICS_INTERP_VERT(
314  c If nonlinear free surf is active, need averaged pressures       I                         listId, md, ndId, ip, im, lm,
315  #ifdef NONLIN_FRSURF       U                         qtmp1, qtmp2,
316            if(select_rStar.GT.0)then       I                         undefRL, myTime, myIter, myThid )
317             call diagnostics_get_pointers('RSURF   ',ipoint1,jpoint1,               ELSE
318       .                                                           myThid)                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
319             call diagnostics_get_pointers('PRESSURE',ipoint2,jpoint2,       &           'INTERP_VERT not allowed in this config'
320       .                                                           myThid)                 CALL PRINT_ERROR( msgBuf , myThid )
321  C if fizhi is being  used, may need to get physics grid pressures                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
322  #ifdef ALLOW_FIZHI               ENDIF
323             if(gdiag(ndId)(10:10) .EQ. 'L')then              ENDIF
324             call diagnostics_get_pointers('FIZPRES ',ipoint2,jpoint2,              IF ( fflags(listId)(2:2).EQ.'I' ) THEN
325       .                                                           myThid)  C-          Integrate vertically: for now, output field has just 1 level:
326             endif                CALL DIAGNOSTICS_SUM_LEVELS(
327  #endif       I                         listId, md, ndId, ip, im, lm,
328             if( jpoint1.ne.0 .and. jpoint2.ne.0) foundp = .true.       U                         qtmp1,
329         I                         undefRL, myTime, myIter, myThid )
330             if(.not. foundp) then              ENDIF
             WRITE(msgBuf,'(3A)') 'DIAGNOSTICS_OUT: ',  
      .    ' Have asked for pressure interpolation but have not ',  
      .    ' Activated surface and 3D pressure diagnostic, ',  
      .    ' RSURF and PRESSURE'  
             CALL PRINT_ERROR( msgBuf , myThid )  
             STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'  
            endif  
   
            DO bj = myByLo(myThid), myByHi(myThid)  
             DO bi = myBxLo(myThid), myBxHi(myThid)  
              call getdiag(1.,undef,qtmpsrf(1-OLx,1-OLy,bi,bj),  
      .                       jpoint1,0,ipoint1,0,bi,bj,myThid)  
             ENDDO  
            ENDDO  
            DO bj = myByLo(myThid), myByHi(myThid)  
             DO bi = myBxLo(myThid), myBxHi(myThid)  
              DO k = 1,nlevels(listId)  
               call getdiag(levs(k,listId),undef,  
      .          qtmp2(1-OLx,1-OLy,k,bi,bj),jpoint2,0,ipoint2,0,  
      .          bi,bj,myThid)  
              ENDDO  
             ENDDO  
            ENDDO  
           endif  
 #else  
 C If nonlinear free surf is off, get pressures from rC and rF arrays  
           DO bj = myByLo(myThid), myByHi(myThid)  
            DO bi = myBxLo(myThid), myBxHi(myThid)  
             DO j = 1-OLy,sNy+OLy  
              DO i = 1-OLx,sNx+OLx  
               qtmpsrf(i,j,bi,bj) = rF(1)  
              ENDDO  
             ENDDO  
             DO j = 1-OLy,sNy+OLy  
              DO i = 1-OLx,sNx+OLx  
               DO k = 1,nlevels(listId)  
                qtmp2(i,j,k,bi,bj) = rC(k)  
               ENDDO  
              ENDDO  
             ENDDO  
            ENDDO  
           ENDDO  
 #endif  
 C Load p to the kappa into a temporary array  
           nlevsout = nplevs  
           DO bj = myByLo(myThid), myByHi(myThid)  
            DO bi = myBxLo(myThid), myBxHi(myThid)  
             DO j = 1,sNy  
              DO i = 1,sNx  
               pksrf(i,j) = qtmpsrf(i,j,bi,bj) ** kappa  
               DO k = 1,nlevels(listId)  
                if(gdiag(ndId)(10:10).eq.'R') then  
                 if(hFacC(i,j,nlevels(listId)-k+1,bi,bj).ne.0.) then  
                  qinp(i,j,k) =  qtmp1(i,j,nlevels(listId)-k+1,bi,bj)  
                 else  
                  qinp(i,j,k) =  undef  
                 endif  
                 pkz(i,j,k) = qtmp2(i,j,nlevels(listId)-k+1,bi,bj)**kappa  
                elseif(gdiag(ndId)(10:10).eq.'L') then  
                 qinp(i,j,k) =  qtmp1(i,j,k,bi,bj)  
                 pkz(i,j,k) = qtmp2(i,j,k,bi,bj)**kappa  
                endif  
               ENDDO  
              ENDDO  
             ENDDO  
   
             DO k = 1,nplevs  
              if(fflags(listId)(3:3).eq.'1') then  
               p = plevs1(k)  
              elseif(fflags(listId)(3:3).eq.'2')then  
               p = plevs2(k)  
              endif  
              call prestopres(qprs(1,1,k),qinp,pkz,pksrf,0.,p,sNx,sNy,  
      .                                                 nlevels(listId) )  
             ENDDO  
331    
332              DO j = 1,sNy  C--     End of empty diag / not-empty diag block
333               DO i = 1,sNx            ENDIF
               DO k = 1,nlevsout  
                qtmp1(i,j,k,bi,bj) =  qprs(i,j,k)  
                if(qtmp1(i,j,k,bi,bj).eq.undef) qtmp1(i,j,k,bi,bj) = 0.  
               ENDDO  
              ENDDO  
             ENDDO  
            ENDDO  
           ENDDO  
334    
335           endif  C--     Ready to write field "md", element "lm" in averageCycle(listId)
336    
337  #ifdef ALLOW_MDSIO  C-        write to binary file, using MDSIO pkg:
338  C         Prepare for mdsio optionality            IF ( diag_mdsio ) THEN
339            IF (diag_mdsio) THEN              nRec = lm + (md-1)*averageCycle(listId)
340              IF (fflags(listId)(1:1) .EQ. ' ') THEN  C           default precision for output files
341  C             This is the old default behavior              prec = writeBinaryPrec
342                CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,'RL',  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
343       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
344              ELSEIF (fflags(listId)(1:1) .EQ. 'R') THEN              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
345  C             Force it to be 32-bit precision  C         a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R
346                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,'RL',              CALL WRITE_REC_LEV_RL(
347       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       I                            fn, prec,
348              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN       I                            NrMax, 1, nLevOutp,
349  C             Force it to be 64-bit precision       I                            qtmp1, -nRec, myIter, myThid )
               CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,'RL',  
      &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)  
             ENDIF  
350            ENDIF            ENDIF
 #endif /*  ALLOW_MDSIO  */  
351    
352  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
353            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
354                CALL DIAGNOSTICS_MNC_OUT(
355              _BEGIN_MASTER( myThid )       I                       NrMax, nLevOutp, listId, ndId,
356         I                       diag_mnc_bn,
357              DO ii = 1,CW_DIMS       I                       useMissingValue, misValLoc,
358                d_cw_name(1:NLEN) = dn_blnk(1:NLEN)       I                       qtmp1,
359                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', nlevels(listId)  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevels(listId)  
             ENDIF  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zld', nlevels(listId)  
             ENDIF  
             IF ( (gdiag(ndId)(10:10) .EQ. 'R')  
      &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN  
               WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)  
             ENDIF  
             dim(3) = Nr+Nrphys  
             ib(3)  = 1  
             ie(3)  = nlevels(listId)  
   
 C           Time dimension  
             dn(4)(1:1) = 'T'  
             dim(4) = -1  
             ib(4)  = 1  
             ie(4)  = 1  
   
             CALL MNC_CW_ADD_GNAME(d_cw_name, 4,  
      &             dim, dn, ib, ie, myThid)  
             CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,  
      &             4,5, myThid)  
             CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',  
      &             tdiag(ndId),myThid)  
             CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',  
      &             udiag(ndId),myThid)  
   
             IF ((fflags(listId)(1:1) .EQ. ' ')  
      &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN  
               CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,  
      &             cdiag(ndId), qtmp1, myThid)  
             ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN  
               CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,  
      &             cdiag(ndId), qtmp1, myThid)  
             ENDIF  
               
             CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)  
             CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)  
   
             _END_MASTER( myThid )  
   
360            ENDIF            ENDIF
361  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
362    
363    C--      end loop on lm (or ll if ALLOW_MNC) counter
364             ENDDO
365  C--     end of Processing Fld # md  C--     end of Processing Fld # md
366          ENDIF          ENDIF
367           ENDDO
368    
369    #ifdef ALLOW_MNC
370    C--   end loop on jj counter
371        ENDDO        ENDDO
372    #endif
373    
374    #ifdef ALLOW_MDSIO
375          IF (diag_mdsio) THEN
376    C-    Note: temporary: since it is a pain to add more arguments to
377    C     all MDSIO S/R, uses instead this specific S/R to write only
378    C     meta files but with more informations in it.
379                glf = globalFiles
380                nRec = nfields(listId)*averageCycle(listId)
381                CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
382         &              0, 0, nLevOutp, ' ',
383         &              nfields(listId), flds(1,listId), nTimRec, timeRec,
384         &              nRec, myIter, myThid)
385          ENDIF
386    #endif /*  ALLOW_MDSIO  */
387    
388        RETURN        RETURN
389        END        END
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.53
  ViewVC Help
Powered by ViewVC 1.1.22