/[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.44 by jmc, Thu Jan 7 01:09:40 2010 UTC
# 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     lm    :: loop index (averageCycle)
52  C     md    :: field number in the list "listId".  C     md    :: field number in the list "listId".
53  C     ndId  :: diagnostics  Id number (in available diagnostics list)  C     ndId  :: diagnostics  Id number (in available diagnostics list)
54  C     mate  :: counter mate Id number (in available diagnostics list)  C     mate  :: counter mate Id number (in available diagnostics list)
55  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
56  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
57        INTEGER i, j, k  C
58    C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
59    C     qtmp1 :: thread-shared temporary array (needs to be in common block):
60    C              to write a diagnostic field to file, copy it first from (big)
61    C              diagnostic storage qdiag into it.
62          COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
63          _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
64    
65          INTEGER i, j, k, lm
66        INTEGER bi, bj        INTEGER bi, bj
67        INTEGER md, ndId, ip, im        INTEGER md, ndId, ip, im
68        INTEGER mate, mVec        INTEGER mate, mVec
69        CHARACTER*8 parms1        CHARACTER*10 gcode
70        CHARACTER*3 mate_index        _RL undef
71        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)        _RL tmpLev
       _RL qtmpsrf(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)  
       _RL qtmp2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)  
       _RL undef, getcon  
       EXTERNAL getcon  
       INTEGER ILNBLNK  
       EXTERNAL ILNBLNK  
72        INTEGER ilen        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./  
73    
74        INTEGER ioUnit        INTEGER ioUnit
75        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
76        CHARACTER*(MAX_LEN_MBUF) suff        CHARACTER*(MAX_LEN_MBUF) suff
77        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
78          INTEGER prec, nRec, nTimRec
79          _RL     timeRec(2)
80    #ifdef ALLOW_MDSIO
81        LOGICAL glf        LOGICAL glf
82    #endif
83  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
84        INTEGER ii        INTEGER ll, llMx, jj, jjMx
85          INTEGER ii, klev
86        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
       CHARACTER*(5) ctmp  
87        INTEGER CW_DIMS, NLEN        INTEGER CW_DIMS, NLEN
88        PARAMETER ( CW_DIMS = 10 )        PARAMETER ( CW_DIMS = 10 )
89        PARAMETER ( NLEN    = 80 )        PARAMETER ( NLEN    = 80 )
# Line 103  C     im    :: counter-mate pointer to s Line 91  C     im    :: counter-mate pointer to s
91        CHARACTER*(NLEN) dn(CW_DIMS)        CHARACTER*(NLEN) dn(CW_DIMS)
92        CHARACTER*(NLEN) d_cw_name        CHARACTER*(NLEN) d_cw_name
93        CHARACTER*(NLEN) dn_blnk        CHARACTER*(NLEN) dn_blnk
94        _RS ztmp(Nr+Nrphys)  #ifdef DIAG_MNC_COORD_NEEDSWORK
95          CHARACTER*(5) ctmp
96          _RS ztmp(NrMax)
97    #endif
98          LOGICAL useMissingValue, useMisValForThisDiag
99          REAL*8 misvalLoc
100          REAL*8 misval_r8(2)
101          REAL*4 misval_r4(2)
102          INTEGER misvalIntLoc, misval_int(2)
103  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
104    
105  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
106    
107    C---  set file properties
108        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
109          undef = UNSET_RL
110    #ifdef ALLOW_FIZHI
111    c     IF ( useFIZHI ) undef = getcon('UNDEF')
112        undef = getcon('UNDEF')        undef = getcon('UNDEF')
113        kappa = getcon('KAPPA')  #endif
       glf = globalFiles  
114        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
115        ilen = ILNBLNK(fnames(listId))        ilen = ILNBLNK(fnames(listId))
116        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)
117    
118    C--   Set time information:
119          IF ( freq(listId).LT.0. ) THEN
120    C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
121            nTimRec = 1
122            timeRec(1) = myTime
123          ELSE
124    C-    Time-average: store the 2 edges of the time-averaging interval.
125    C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
126    C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
127            nTimRec = 2
128    
129    C-    end of time-averaging interval:
130            timeRec(2) = myTime
131    
132    C-    begining of time-averaging interval:
133    c       timeRec(1) = myTime - freq(listId)
134    C     a) find the time of the previous multiple of output freq:
135            timeRec(1) = myTime-deltaTClock*0.5 _d 0
136            timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
137            i = INT( timeRec(1) )
138            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
139    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
140            timeRec(1) = MAX( timeRec(1), startTime )
141    
142    C     b) round off to nearest multiple of time-step:
143            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
144            i = NINT( timeRec(1) )
145    C     if just half way, NINT will return the next time-step: correct this
146            tmpLev = FLOAT(i) - 0.5 _d 0
147            IF ( timeRec(1).EQ.tmpLev ) i = i - 1
148            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
149    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
150          ENDIF
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)
162          ENDDO        ENDIF
163          WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)        DO jj=1,jjMx
164    
165           IF (useMNC .AND. diag_mnc) THEN
166    C     Handle missing value attribute (land points)
167             useMissingValue = .FALSE.
168    #ifdef DIAGNOSTICS_MISSING_VALUE
169             useMissingValue = .TRUE.
170    #endif /* DIAGNOSTICS_MISSING_VALUE */
171             IF ( misvalFlt(listId) .NE. UNSET_RL ) THEN
172              misvalLoc = misvalFlt(listId)
173             ELSE
174              misvalLoc = undef
175             ENDIF
176    C     Defaults to UNSET_I
177             misvalIntLoc = misvalInt(listId)
178             DO ii=1,2
179    C         misval_r4(ii)  = UNSET_FLOAT4
180    C         misval_r8(ii)  = UNSET_FLOAT8
181              misval_r4(ii)  = misvalLoc
182              misval_r8(ii)  = misvalLoc
183              misval_int(ii) = UNSET_I
184             ENDDO
185             DO i = 1,MAX_LEN_FNAM
186               diag_mnc_bn(i:i) = ' '
187             ENDDO
188             DO i = 1,NLEN
189               dn_blnk(i:i) = ' '
190             ENDDO
191             WRITE( diag_mnc_bn, '(A)' ) fnames(listId)(1:ilen)
192    
193  C       Update the record dimension by writing the iteration number  C       Update the record dimension by writing the iteration number
194          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)           klev = myIter + jj - jjMx
195          CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',myTime,myThid)           tmpLev = myTime + deltaTClock*(jj -jjMx)
196          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)           CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)
197             CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',tmpLev,myThid)
198          dn(1)(1:NLEN) = dn_blnk(1:NLEN)           CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)
199          WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)           CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',klev,myThid)
200          dim(1) = nlevels(listId)  
201          ib(1)  = 1  C       NOTE: at some point it would be a good idea to add a time_bounds
202          ie(1)  = nlevels(listId)  C       variable that has dimension (2,T) and clearly denotes the
203    C       beginning and ending times for each diagnostics period
204          CALL MNC_CW_ADD_GNAME('diag_levels', 1,  
205       &       dim, dn, ib, ie, myThid)           dn(1)(1:NLEN) = dn_blnk(1:NLEN)
206          CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',           WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)
207       &       0,0, myThid)           dim(1) = nlevels(listId)
208          CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',           ib(1)  = 1
209       &       'Idicies of vertical levels within the source arrays',           ie(1)  = nlevels(listId)
210       &       myThid)  
211                     CALL MNC_CW_ADD_GNAME('diag_levels', 1,
212          CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,       &        dim, dn, ib, ie, myThid)
213       &       'diag_levels', levs(1,listId), myThid)           CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',
214         &        0,0, myThid)
215             CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',
216         &        'Idicies of vertical levels within the source arrays',
217         &        myThid)
218    C     suppress the missing value attribute (iflag = 0)
219             IF (useMissingValue)
220         &       CALL MNC_CW_VATTR_MISSING('diag_levels', 0,
221         I       misval_r8, misval_r4, misval_int,
222         I       myThid )
223    
224             CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
225         &        'diag_levels', levs(1,listId), myThid)
226    
227          CALL MNC_CW_DEL_VNAME('diag_levels', myThid)           CALL MNC_CW_DEL_VNAME('diag_levels', myThid)
228          CALL MNC_CW_DEL_GNAME('diag_levels', myThid)           CALL MNC_CW_DEL_GNAME('diag_levels', myThid)
229    
230  #ifdef DIAG_MNC_COORD_NEEDSWORK  #ifdef DIAG_MNC_COORD_NEEDSWORK
231  C       This part has been placed in an #ifdef because, as its currently  C       This part has been placed in an #ifdef because, as its currently
# Line 158  C       grid.  As we start using diagnos Line 234  C       grid.  As we start using diagnos
234  C       levels, land levels, etc. the different vertical coordinate  C       levels, land levels, etc. the different vertical coordinate
235  C       dimensions will have to be taken into account.  C       dimensions will have to be taken into account.
236    
237    C       20051021 JMC & EH3 : We need to extend this so that a few
238    C       variables each defined on different grids do not have the same
239    C       vertical dimension names so we should be using a pattern such
240    C       as: Z[uml]td000000 where the 't' is the type as specified by
241    C       gdiag(10)
242    
243  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx
244          ctmp(1:5) = 'mul  '           ctmp(1:5) = 'mul  '
245          DO i = 1,3           DO i = 1,3
246            dn(1)(1:NLEN) = dn_blnk(1:NLEN)             dn(1)(1:NLEN) = dn_blnk(1:NLEN)
247            WRITE(dn(1),'(3a,i6.6)') 'Z',ctmp(i:i),'d',nlevels(listId)             WRITE(dn(1),'(3a,i6.6)') 'Z',ctmp(i:i),'d',nlevels(listId)
248            CALL MNC_CW_ADD_GNAME(dn(1), 1, dim, dn, ib, ie, myThid)             CALL MNC_CW_ADD_GNAME(dn(1), 1, dim, dn, ib, ie, myThid)
249            CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)             CALL MNC_CW_ADD_VNAME(dn(1), dn(1), 0,0, myThid)
250    
251  C         The following three ztmp() loops should eventually be modified  C         The following three ztmp() loops should eventually be modified
252  C         to reflect the fractional nature of levs(j,l) -- they should  C         to reflect the fractional nature of levs(j,l) -- they should
253  C         do something like:  C         do something like:
254  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))
255  C                      + ( rC(INT(FLOOR(levs(j,l))))  C                      + ( rC(INT(FLOOR(levs(j,l))))
256  C                          + rC(INT(CEIL(levs(j,l)))) )  C                          + rC(INT(CEIL(levs(j,l)))) )
257  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )
258  C         for averaged levels.  C         for averaged levels.
259            IF (i .EQ. 1) THEN             IF (i .EQ. 1) THEN
260              DO j = 1,nlevels(listId)               DO j = 1,nlevels(listId)
261                ztmp(j) = rC(NINT(levs(j,listId)))                 ztmp(j) = rC(NINT(levs(j,listId)))
262              ENDDO               ENDDO
263              CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',               CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
264       &           'Dimensional coordinate value at the mid point',       &            'Dimensional coordinate value at the mid point',
265       &           myThid)       &            myThid)
266            ELSEIF (i .EQ. 2) THEN             ELSEIF (i .EQ. 2) THEN
267              DO j = 1,nlevels(listId)               DO j = 1,nlevels(listId)
268                ztmp(j) = rF(NINT(levs(j,listId)) + 1)                 ztmp(j) = rF(NINT(levs(j,listId)) + 1)
269              ENDDO               ENDDO
270              CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',               CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
271       &           'Dimensional coordinate value at the upper point',       &            'Dimensional coordinate value at the upper point',
272       &           myThid)       &            myThid)
273            ELSEIF (i .EQ. 3) THEN             ELSEIF (i .EQ. 3) THEN
274              DO j = 1,nlevels(listId)               DO j = 1,nlevels(listId)
275                ztmp(j) = rF(NINT(levs(j,listId)))                 ztmp(j) = rF(NINT(levs(j,listId)))
276              ENDDO               ENDDO
277              CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',               CALL MNC_CW_ADD_VATTR_TEXT(dn(1),'description',
278       &           'Dimensional coordinate value at the lower point',       &            'Dimensional coordinate value at the lower point',
279       &           myThid)       &            myThid)
280            ENDIF             ENDIF
281            CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)  C     suppress the missing value attribute (iflag = 0)
282            CALL MNC_CW_DEL_VNAME(dn(1), myThid)             IF (useMissingValue)
283            CALL MNC_CW_DEL_GNAME(dn(1), myThid)       &          CALL MNC_CW_VATTR_MISSING(dn(1), 0,
284          ENDDO       I          misval_r8, misval_r4, misval_int,
285         I          myThid )
286               CALL MNC_CW_RS_W('D',diag_mnc_bn,0,0, dn(1), ztmp, myThid)
287               CALL MNC_CW_DEL_VNAME(dn(1), myThid)
288               CALL MNC_CW_DEL_GNAME(dn(1), myThid)
289             ENDDO
290  #endif /*  DIAG_MNC_COORD_NEEDSWORK  */  #endif /*  DIAG_MNC_COORD_NEEDSWORK  */
291    
292        ENDIF         ENDIF
293  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
294    
295        DO md = 1,nfields(listId)  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
296    
297           DO md = 1,nfields(listId)
298          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
299          parms1 = gdiag(ndId)(1:8)          gcode = gdiag(ndId)(1:10)
300          IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN          mate = 0
301            mVec = 0
302            IF ( gcode(5:5).EQ.'C' ) THEN
303    C-      Check for Mate of a Counter Diagnostic
304               mate = hdiag(ndId)
305            ELSEIF ( gcode(1:1).EQ.'U' .OR. gcode(1:1).EQ.'V' ) THEN
306    C-      Check for Mate of a Vector Diagnostic
307               mVec = hdiag(ndId)
308            ENDIF
309            IF ( idiag(md,listId).NE.0 .AND. gcode(5:5).NE.'D' ) THEN
310  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
311    #ifdef ALLOW_MNC
312             DO ll=1,llMx
313              lm = jj+ll-1
314    #else
315             DO lm=1,averageCycle(listId)
316    #endif
317    
318            ip = ABS(idiag(md,listId))            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
319            im = mdiag(md,listId)            im = mdiag(md,listId)
320              IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
321              IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
322    
323            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ndiag(ip,1,1).EQ.0 ) THEN
324  C-        Empty diagnostics case :  C-        Empty diagnostics case :
325    
# Line 221  C-        Empty diagnostics case : Line 328  C-        Empty diagnostics case :
328       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter       &        '- WARNING - from DIAGNOSTICS_OUT at iter=', myIter
329              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
330       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
331              WRITE(msgBuf,'(A,I4,3A,I3,2A)')              WRITE(msgBuf,'(A,I6,3A,I4,2A)')
332       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),       &       '- WARNING -   diag.#',ndId, ' : ',flds(md,listId),
333       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
334              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
335       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
336              WRITE(msgBuf,'(A,I2,A)')              IF ( averageCycle(listId).GT.1 ) THEN
337       &       '- WARNING -   has not been filled (ndiag=',               WRITE(msgBuf,'(A,2(I3,A))')
338       &       ndiag(ip,1,1), ' )'       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
339         &                                            ndiag(ip,1,1), ' )'
340                ELSE
341                 WRITE(msgBuf,'(A,2(I3,A))')
342         &        '- WARNING -   has not been filled (ndiag=',
343         &                                            ndiag(ip,1,1), ' )'
344                ENDIF
345              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
346       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
347              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
# Line 251  C-        Empty diagnostics case : Line 364  C-        Empty diagnostics case :
364            ELSE            ELSE
365  C-        diagnostics is not empty :  C-        diagnostics is not empty :
366    
367              IF ( myThid.EQ.1 ) WRITE(ioUnit,'(A,I3,3A,I8,2A)')              IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN
368                  WRITE(ioUnit,'(A,I6,3A,I8,2A)')
369       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
370       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
371                  IF ( mate.GT.0 ) THEN
372              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)')  
373       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
374       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
375                  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  
376                  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
377                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
378       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
379       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
380       &             ' exists '       &             ' exists '
381                  ELSE                  ELSE
382                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I6,3A)')
383       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
384       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
385       &             ' not enabled'       &             ' not enabled'
# Line 286  C             -------------------------- Line 387  C             --------------------------
387                ENDIF                ENDIF
388              ENDIF              ENDIF
389    
390              DO bj = myByLo(myThid), myByHi(myThid)              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
391               DO bi = myBxLo(myThid), myBxHi(myThid)  C-       get all the levels (for vertical interpolation)
392                DO k = 1,nlevels(listId)                DO bj = myByLo(myThid), myByHi(myThid)
393                  CALL GETDIAG(                 DO bi = myBxLo(myThid), myBxHi(myThid)
394       I                       levs(k,listId),undef,                  DO k = 1,kdiag(ndId)
395       O                       qtmp1(1-OLx,1-OLy,k,bi,bj),                    tmpLev = k
396       I                       ndId,mate,ip,im,bi,bj,myThid)                    CALL GETDIAG(
397         I                         tmpLev,undef,
398         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
399         I                         ndId,mate,ip,im,bi,bj,myThid)
400                    ENDDO
401                   ENDDO
402                ENDDO                ENDDO
403               ENDDO              ELSE
404              ENDDO  C-       get only selected levels:
405                  DO bj = myByLo(myThid), myByHi(myThid)
406                   DO bi = myBxLo(myThid), myBxHi(myThid)
407                    DO k = 1,nlevels(listId)
408                      CALL GETDIAG(
409         I                         levs(k,listId),undef,
410         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
411         I                         ndId,mate,ip,im,bi,bj,myThid)
412                    ENDDO
413                   ENDDO
414                  ENDDO
415                ENDIF
416    
417  C-        end of empty diag / not empty block  C-        end of empty diag / not empty block
418            ENDIF            ENDIF
419    
           nlevsout = nlevels(listId)  
   
420  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
421  C             Check to see if we need to interpolate before output  C         Check to see if we need to interpolate before output
422  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
423  C  (we are still inside field exist if sequence and field do loop)            IF ( fflags(listId)(2:2).EQ.'P' ) THEN
424  C  C-        Do vertical interpolation:
425               IF ( fluidIsAir ) THEN
426           if(fflags(listId)(2:2).eq.'P') then  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
427                CALL DIAGNOSTICS_INTERP_VERT(
428  c If nonlinear free surf is active, need averaged pressures       I                     listId, md, ndId, ip, im, lm,
429  #ifdef NONLIN_FRSURF       U                     qtmp1,
430            if(select_rStar.GT.0)then       I                     undef, myTime, myIter, myThid )
431             call diagnostics_get_pointers('RSURF   ',ipoint1,jpoint1,             ELSE
432       .                                                           myThid)               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
433             call diagnostics_get_pointers('PRESSURE',ipoint2,jpoint2,       &         'INTERP_VERT not allowed in this config'
434       .                                                           myThid)               CALL PRINT_ERROR( msgBuf , myThid )
435  C if fizhi is being  used, may need to get physics grid pressures               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
436  #ifdef ALLOW_FIZHI             ENDIF
437             if(gdiag(ndId)(10:10) .EQ. 'L')then            ENDIF
            call diagnostics_get_pointers('FIZPRES ',ipoint2,jpoint2,  
      .                                                           myThid)  
            endif  
 #endif  
            if( jpoint1.ne.0 .and. jpoint2.ne.0) foundp = .true.  
   
            if(.not. foundp) then  
             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  
   
             DO j = 1,sNy  
              DO i = 1,sNx  
               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  
438    
439           endif  C--    Ready to write field "md", element "lm" in averageCycle(listId)
440    
441  #ifdef ALLOW_MDSIO  C-        write to binary file, using MDSIO pkg:
442  C         Prepare for mdsio optionality            IF ( diag_mdsio ) THEN
443            IF (diag_mdsio) THEN              nRec = lm + (md-1)*averageCycle(listId)
444              IF (fflags(listId)(1:1) .EQ. ' ') THEN  C           default precision for output files
445  C             This is the old default behavior              prec = writeBinaryPrec
446                CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,'RL',  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
447       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
448              ELSEIF (fflags(listId)(1:1) .EQ. 'R') THEN              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
449  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
450                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,'RL',              CALL WRITE_REC_LEV_RL(
451       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       I                            fn, prec,
452              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN       I                            NrMax, 1, nlevels(listId),
453  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  
454            ENDIF            ENDIF
 #endif /*  ALLOW_MDSIO  */  
455    
456  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
457            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
# Line 459  C           XY dimensions Line 477  C           XY dimensions
477              dim(2)       = sNy + 2*OLy              dim(2)       = sNy + 2*OLy
478              ib(1)        = OLx + 1              ib(1)        = OLx + 1
479              ib(2)        = OLy + 1              ib(2)        = OLy + 1
480              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN
481                dn(1)(1:2) = 'X'                dn(1)(1:2) = 'X'
482                ie(1)      = OLx + sNx                ie(1)      = OLx + sNx
483                dn(2)(1:2) = 'Y'                dn(2)(1:2) = 'Y'
# Line 480  C           XY dimensions Line 498  C           XY dimensions
498                dn(2)(1:3) = 'Yp1'                dn(2)(1:3) = 'Yp1'
499                ie(2)      = OLy + sNy + 1                ie(2)      = OLy + sNy + 1
500              ENDIF              ENDIF
501                
502  C           Z is special since it varies  C           Z is special since it varies
503              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)
504              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
# Line 495  C           Z is special since it varies Line 513  C           Z is special since it varies
513       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
514                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)
515              ENDIF              ENDIF
516              dim(3) = Nr+Nrphys              dim(3) = NrMax
517              ib(3)  = 1              ib(3)  = 1
518              ie(3)  = nlevels(listId)              ie(3)  = nlevels(listId)
519    
# Line 505  C           Time dimension Line 523  C           Time dimension
523              ib(4)  = 1              ib(4)  = 1
524              ie(4)  = 1              ie(4)  = 1
525    
526              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,
527       &             dim, dn, ib, ie, myThid)       &             dim, dn, ib, ie, myThid)
528              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,
529       &             4,5, myThid)       &             4,5, myThid)
530              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',
531       &             tdiag(ndId),myThid)       &             tdiag(ndId),myThid)
532              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',
533       &             udiag(ndId),myThid)       &             udiag(ndId),myThid)
534    
535              IF ((fflags(listId)(1:1) .EQ. ' ')  C     Missing values only for scalar diagnostics at mass points (so far)
536       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN              useMisValForThisDiag = useMissingValue
537         &           .AND.gdiag(ndId)(1:2).EQ.'SM'
538                IF ( useMisValForThisDiag ) THEN
539    C     assign missing values and set flag for adding the netCDF atttibute
540                 CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 2,
541         I            misval_r8, misval_r4, misval_int,
542         I            myThid )
543    C     and now use the missing values for masking out the land points
544                 DO bj = myByLo(myThid), myByHi(myThid)
545                  DO bi = myBxLo(myThid), myBxHi(myThid)
546                   DO k = 1,nlevels(listId)
547                    klev = NINT(levs(k,listId))
548                    DO j = 1-OLy,sNy+OLy
549                     DO i = 1-OLx,sNx+OLx
550                      IF ( maskC(i,j,klev,bi,bj) .EQ. 0. )
551         &                 qtmp1(i,j,k,bi,bj) = misvalLoc
552                     ENDDO
553                    ENDDO
554                   ENDDO
555                  ENDDO
556                 ENDDO
557                ELSE
558    C     suppress the missing value attribute (iflag = 0)
559    C     Note: We have to call the following subroutine for each mnc that has
560    C     been created "on the fly" by mnc_cw_add_vname and will be deleted
561    C     by mnc_cw_del_vname, because all of these variables use the same
562    C     identifier so that mnc_cw_vfmv(indv) needs to be overwritten for
563    C     each of these variables
564                 CALL MNC_CW_VATTR_MISSING(cdiag(ndId), 0,
565         I            misval_r8, misval_r4, misval_int,
566         I            myThid )
567                ENDIF
568    
569                IF (  ((writeBinaryPrec .EQ. precFloat32)
570         &            .AND. (fflags(listId)(1:1) .NE. 'D'))
571         &             .OR. (fflags(listId)(1:1) .EQ. 'R') ) THEN
572                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,
573       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
574              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              ELSEIF ( (writeBinaryPrec .EQ. precFloat64)
575         &             .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN
576                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
577       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
578              ENDIF              ENDIF
579                
580              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)
581              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)
582    
# Line 531  C           Time dimension Line 585  C           Time dimension
585            ENDIF            ENDIF
586  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
587    
588    C--      end loop on lm (or ll if ALLOW_MNC) counter
589             ENDDO
590  C--     end of Processing Fld # md  C--     end of Processing Fld # md
591          ENDIF          ENDIF
592           ENDDO
593    
594    #ifdef ALLOW_MNC
595    C--   end loop on jj counter
596        ENDDO        ENDDO
597    #endif
598    
599    #ifdef ALLOW_MDSIO
600          IF (diag_mdsio) THEN
601    C-    Note: temporary: since it's a pain to add more arguments to
602    C     all MDSIO S/R, uses instead this specific S/R to write only
603    C     meta files but with more informations in it.
604                glf = globalFiles
605                nRec = nfields(listId)*averageCycle(listId)
606                CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
607         &              0, 0, nlevels(listId), ' ',
608         &              nfields(listId), flds(1,listId), nTimRec, timeRec,
609         &              nRec, myIter, myThid)
610          ENDIF
611    #endif /*  ALLOW_MDSIO  */
612    
613        RETURN        RETURN
614        END        END
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.44
  ViewVC Help
Powered by ViewVC 1.1.22