/[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.32 by jmc, Fri Dec 29 23:57:15 2006 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          INTEGER NrMax
30  #ifdef ALLOW_FIZHI  #ifdef ALLOW_FIZHI
31  #include "fizhi_SIZE.h"  #include "fizhi_SIZE.h"
32          PARAMETER( NrMax = Nr+Nrphys )
33  #else  #else
34        INTEGER Nrphys        PARAMETER( NrMax = Nr )
       PARAMETER (Nrphys=0)  
35  #endif  #endif
36    
   
37  C     !INPUT PARAMETERS:  C     !INPUT PARAMETERS:
38  C     listId  :: Diagnostics list number being written  C     listId  :: Diagnostics list number being written
39  C     myIter  :: current iteration number  C     myIter  :: current iteration number
# Line 45  CEOP Line 45  CEOP
45    
46  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
47  C     i,j,k :: loop indices  C     i,j,k :: loop indices
48    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)
51  C     mate  :: counter mate Id number (in available diagnostics list)  C     mate  :: counter mate Id number (in available diagnostics list)
52  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
53  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
54        INTEGER i, j, k  C
55    C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
56    C     qtmp1 :: thread-shared temporary array (needs to be in common block):
57    C              to write a diagnostic field to file, copy it first from (big)
58    C              diagnostic storage qdiag into it.
59          COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
60          _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
61    
62          INTEGER i, j, k, lm
63        INTEGER bi, bj        INTEGER bi, bj
64        INTEGER md, ndId, ip, im        INTEGER md, ndId, ip, im
65        INTEGER mate, mVec        INTEGER mate, mVec
66        CHARACTER*8 parms1        CHARACTER*8 parms1
       CHARACTER*3 mate_index  
       _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)  
       _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)  
67        _RL undef, getcon        _RL undef, getcon
68          _RL tmpLev
69        EXTERNAL getcon        EXTERNAL getcon
70        INTEGER ILNBLNK        INTEGER ILNBLNK
71        EXTERNAL 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    #ifdef ALLOW_MDSIO
79        LOGICAL glf        LOGICAL glf
80          INTEGER nRec
81          INTEGER prec
82    #endif
83  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
84        INTEGER ii        INTEGER ii
85        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
       CHARACTER*(5) ctmp  
86        INTEGER CW_DIMS, NLEN        INTEGER CW_DIMS, NLEN
87        PARAMETER ( CW_DIMS = 10 )        PARAMETER ( CW_DIMS = 10 )
88        PARAMETER ( NLEN    = 80 )        PARAMETER ( NLEN    = 80 )
# Line 103  C     im    :: counter-mate pointer to s Line 90  C     im    :: counter-mate pointer to s
90        CHARACTER*(NLEN) dn(CW_DIMS)        CHARACTER*(NLEN) dn(CW_DIMS)
91        CHARACTER*(NLEN) d_cw_name        CHARACTER*(NLEN) d_cw_name
92        CHARACTER*(NLEN) dn_blnk        CHARACTER*(NLEN) dn_blnk
93        _RS ztmp(Nr+Nrphys)  #ifdef DIAG_MNC_COORD_NEEDSWORK
94          CHARACTER*(5) ctmp
95          _RS ztmp(NrMax)
96    #endif
97  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
98    
99  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
100    
101        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
102        undef = getcon('UNDEF')        undef = getcon('UNDEF')
       kappa = getcon('KAPPA')  
       glf = globalFiles  
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)
# Line 130  C       Update the record dimension by w Line 118  C       Update the record dimension by w
118          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)
119          CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',myTime,myThid)          CALL MNC_CW_RL_W_S('D',diag_mnc_bn,0,0,'T',myTime,myThid)
120          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)
121            CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',myIter,myThid)
122    
123    C       NOTE: at some point it would be a good idea to add a time_bounds
124    C       variable that has dimension (2,T) and clearly denotes the
125    C       beginning and ending times for each diagnostics period
126    
127          dn(1)(1:NLEN) = dn_blnk(1:NLEN)          dn(1)(1:NLEN) = dn_blnk(1:NLEN)
128          WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)          WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)
# Line 137  C       Update the record dimension by w Line 130  C       Update the record dimension by w
130          ib(1)  = 1          ib(1)  = 1
131          ie(1)  = nlevels(listId)          ie(1)  = nlevels(listId)
132    
133          CALL MNC_CW_ADD_GNAME('diag_levels', 1,          CALL MNC_CW_ADD_GNAME('diag_levels', 1,
134       &       dim, dn, ib, ie, myThid)       &       dim, dn, ib, ie, myThid)
135          CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',          CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',
136       &       0,0, myThid)       &       0,0, myThid)
137          CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',          CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',
138       &       'Idicies of vertical levels within the source arrays',       &       'Idicies of vertical levels within the source arrays',
139       &       myThid)       &       myThid)
140            
141          CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,          CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
142       &       'diag_levels', levs(1,listId), myThid)       &       'diag_levels', levs(1,listId), myThid)
143    
# Line 158  C       grid.  As we start using diagnos Line 151  C       grid.  As we start using diagnos
151  C       levels, land levels, etc. the different vertical coordinate  C       levels, land levels, etc. the different vertical coordinate
152  C       dimensions will have to be taken into account.  C       dimensions will have to be taken into account.
153    
154    C       20051021 JMC & EH3 : We need to extend this so that a few
155    C       variables each defined on different grids do not have the same
156    C       vertical dimension names so we should be using a pattern such
157    C       as: Z[uml]td000000 where the 't' is the type as specified by
158    C       gdiag(10)
159    
160  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx
161          ctmp(1:5) = 'mul  '          ctmp(1:5) = 'mul  '
162          DO i = 1,3          DO i = 1,3
# Line 169  C       Now define:  Zmdxxxxxx, Zudxxxxx Line 168  C       Now define:  Zmdxxxxxx, Zudxxxxx
168  C         The following three ztmp() loops should eventually be modified  C         The following three ztmp() loops should eventually be modified
169  C         to reflect the fractional nature of levs(j,l) -- they should  C         to reflect the fractional nature of levs(j,l) -- they should
170  C         do something like:  C         do something like:
171  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))
172  C                      + ( rC(INT(FLOOR(levs(j,l))))  C                      + ( rC(INT(FLOOR(levs(j,l))))
173  C                          + rC(INT(CEIL(levs(j,l)))) )  C                          + rC(INT(CEIL(levs(j,l)))) )
174  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )
175  C         for averaged levels.  C         for averaged levels.
# Line 205  C         for averaged levels. Line 204  C         for averaged levels.
204        ENDIF        ENDIF
205  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
206    
207    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
208    
209        DO md = 1,nfields(listId)        DO md = 1,nfields(listId)
210          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
211          parms1 = gdiag(ndId)(1:8)          parms1 = gdiag(ndId)(1:8)
212            mate = 0
213            mVec = 0
214            IF ( parms1(5:5).EQ.'C' ) THEN
215    C-      Check for Mate of a Counter Diagnostic
216               READ(parms1,'(5X,I3)') mate
217            ELSEIF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN
218    C-      Check for Mate of a Vector Diagnostic
219               READ(parms1,'(5X,I3)') mVec
220            ENDIF
221          IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN          IF ( idiag(md,listId).NE.0 .AND. parms1(5:5).NE.'D' ) THEN
222  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
223             DO lm=1,averageCycle(listId)
224    
225            ip = ABS(idiag(md,listId))            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
226            im = mdiag(md,listId)            im = mdiag(md,listId)
227              IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
228              IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
229    
230            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ndiag(ip,1,1).EQ.0 ) THEN
231  C-        Empty diagnostics case :  C-        Empty diagnostics case :
232    
# Line 226  C-        Empty diagnostics case : Line 240  C-        Empty diagnostics case :
240       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
241              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
242       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
243              WRITE(msgBuf,'(A,I2,A)')              IF ( averageCycle(listId).GT.1 ) THEN
244       &       '- WARNING -   has not been filled (ndiag=',               WRITE(msgBuf,'(A,2(I2,A))')
245       &       ndiag(ip,1,1), ' )'       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
246         &                                            ndiag(ip,1,1), ' )'
247                ELSE
248                 WRITE(msgBuf,'(A,2(I2,A))')
249         &        '- WARNING -   has not been filled (ndiag=',
250         &                                            ndiag(ip,1,1), ' )'
251                ENDIF
252              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
253       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
254              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
# Line 251  C-        Empty diagnostics case : Line 271  C-        Empty diagnostics case :
271            ELSE            ELSE
272  C-        diagnostics is not empty :  C-        diagnostics is not empty :
273    
274              IF ( myThid.EQ.1 ) WRITE(ioUnit,'(A,I3,3A,I8,2A)')              IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN
275                  WRITE(ioUnit,'(A,I3,3A,I8,2A)')
276       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
277       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
278                  IF ( mate.GT.0 ) THEN
279              IF ( parms1(5:5).EQ.'C' ) THEN                 WRITE(ioUnit,'(3A,I3,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)')  
280       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
281       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
282                  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  
283                  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
284                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I3,3A)')
285       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
286       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
287       &             ' exists '       &             ' exists '
288                  ELSE                  ELSE
289                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I3,3A)')
290       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
291       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
292       &             ' not enabled'       &             ' not enabled'
# Line 286  C             -------------------------- Line 294  C             --------------------------
294                ENDIF                ENDIF
295              ENDIF              ENDIF
296    
297              DO bj = myByLo(myThid), myByHi(myThid)              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
298               DO bi = myBxLo(myThid), myBxHi(myThid)  C-       get all the levels (for vertical interpolation)
299                DO k = 1,nlevels(listId)                DO bj = myByLo(myThid), myByHi(myThid)
300                  CALL GETDIAG(                 DO bi = myBxLo(myThid), myBxHi(myThid)
301       I                       levs(k,listId),undef,                  DO k = 1,kdiag(ndId)
302       O                       qtmp1(1-OLx,1-OLy,k,bi,bj),                    tmpLev = k
303       I                       ndId,mate,ip,im,bi,bj,myThid)                    CALL GETDIAG(
304         I                         tmpLev,undef,
305         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
306         I                         ndId,mate,ip,im,bi,bj,myThid)
307                    ENDDO
308                   ENDDO
309                ENDDO                ENDDO
310               ENDDO              ELSE
311              ENDDO  C-       get only selected levels:
312                  DO bj = myByLo(myThid), myByHi(myThid)
313                   DO bi = myBxLo(myThid), myBxHi(myThid)
314                    DO k = 1,nlevels(listId)
315                      CALL GETDIAG(
316         I                         levs(k,listId),undef,
317         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
318         I                         ndId,mate,ip,im,bi,bj,myThid)
319                    ENDDO
320                   ENDDO
321                  ENDDO
322                ENDIF
323    
324  C-        end of empty diag / not empty block  C-        end of empty diag / not empty block
325            ENDIF            ENDIF
326    
           nlevsout = nlevels(listId)  
   
327  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
328  C             Check to see if we need to interpolate before output  C         Check to see if we need to interpolate before output
329  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
330  C  (we are still inside field exist if sequence and field do loop)            IF ( fflags(listId)(2:2).EQ.'P' ) THEN
331  C  C-        Do vertical interpolation:
332               IF ( fluidIsAir ) THEN
333           if(fflags(listId)(2:2).eq.'P') then  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
334                CALL DIAGNOSTICS_INTERP_VERT(
335  c If nonlinear free surf is active, need averaged pressures       I                     listId, md, ndId, ip, im, lm,
336  #ifdef NONLIN_FRSURF       U                     qtmp1,
337            if(select_rStar.GT.0)then       I                     undef, myTime, myIter, myThid )
338             call diagnostics_get_pointers('RSURF   ',ipoint1,jpoint1,             ELSE
339       .                                                           myThid)               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
340             call diagnostics_get_pointers('PRESSURE',ipoint2,jpoint2,       &         'INTERP_VERT not allowed in this config'
341       .                                                           myThid)               CALL PRINT_ERROR( msgBuf , myThid )
342  C if fizhi is being  used, may need to get physics grid pressures               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
343  #ifdef ALLOW_FIZHI             ENDIF
344             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  
345    
346           endif  C--    Ready to write field "md", element "lm" in averageCycle(listId)
347    
348  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
349  C         Prepare for mdsio optionality  C-        write to binary file, using MDSIO pkg:
350            IF (diag_mdsio) THEN            IF (diag_mdsio) THEN
351              IF (fflags(listId)(1:1) .EQ. ' ') THEN              glf = globalFiles
352  C             This is the old default behavior              nRec = lm + (md-1)*averageCycle(listId)
353                CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,'RL',  C           default precision for output files
354       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)              prec = writeBinaryPrec
355              ELSEIF (fflags(listId)(1:1) .EQ. 'R') THEN  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
356  C             Force it to be 32-bit precision              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
357                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,'RL',              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
358       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)  c           CALL MDS_WRITE_FIELD(fn,prec,glf,.FALSE.,'RL',
359              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN  c    &              NrMax,nlevels(listId),qtmp1,nRec,myIter,myThid)
360  C             Force it to be 64-bit precision  C         a hack not to write meta files now:
361                CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,'RL',              CALL MDS_WRITE_FIELD(fn,prec,glf,.FALSE.,'RL',
362       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       &              NrMax,nlevels(listId),qtmp1,-nRec,myIter,myThid)
             ENDIF  
363            ENDIF            ENDIF
364  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
365    
# Line 459  C           XY dimensions Line 387  C           XY dimensions
387              dim(2)       = sNy + 2*OLy              dim(2)       = sNy + 2*OLy
388              ib(1)        = OLx + 1              ib(1)        = OLx + 1
389              ib(2)        = OLy + 1              ib(2)        = OLy + 1
390              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN
391                dn(1)(1:2) = 'X'                dn(1)(1:2) = 'X'
392                ie(1)      = OLx + sNx                ie(1)      = OLx + sNx
393                dn(2)(1:2) = 'Y'                dn(2)(1:2) = 'Y'
# Line 480  C           XY dimensions Line 408  C           XY dimensions
408                dn(2)(1:3) = 'Yp1'                dn(2)(1:3) = 'Yp1'
409                ie(2)      = OLy + sNy + 1                ie(2)      = OLy + sNy + 1
410              ENDIF              ENDIF
411                
412  C           Z is special since it varies  C           Z is special since it varies
413              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)
414              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
# Line 495  C           Z is special since it varies Line 423  C           Z is special since it varies
423       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
424                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)
425              ENDIF              ENDIF
426              dim(3) = Nr+Nrphys              dim(3) = NrMax
427              ib(3)  = 1              ib(3)  = 1
428              ie(3)  = nlevels(listId)              ie(3)  = nlevels(listId)
429    
# Line 505  C           Time dimension Line 433  C           Time dimension
433              ib(4)  = 1              ib(4)  = 1
434              ie(4)  = 1              ie(4)  = 1
435    
436              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,
437       &             dim, dn, ib, ie, myThid)       &             dim, dn, ib, ie, myThid)
438              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,
439       &             4,5, myThid)       &             4,5, myThid)
440              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',
441       &             tdiag(ndId),myThid)       &             tdiag(ndId),myThid)
442              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',
443       &             udiag(ndId),myThid)       &             udiag(ndId),myThid)
444    
445              IF ((fflags(listId)(1:1) .EQ. ' ')  C           Per the observations of Baylor, this has been commented out
446    C           until we have code that can write missing_value attributes
447    C           in a way thats compatible with most of the more popular
448    C           netCDF tools including ferret.  Using all-zeros completely
449    C           breaks ferret.
450    
451    C           CALL MNC_CW_ADD_VATTR_DBL(cdiag(ndId),'missing_value',
452    C           &             0.0 _d 0,myThid)
453    
454                IF ( ( (writeBinaryPrec .EQ. precFloat32)
455         &           .AND. (fflags(listId)(1:1) .NE. 'D')
456         &           .AND. (fflags(listId)(1:1) .NE. 'R') )
457       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN
458                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,
459       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
460              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              ELSEIF ( (writeBinaryPrec .EQ. precFloat64)
461         &             .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN
462                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
463       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
464              ENDIF              ENDIF
465                
466              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)
467              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)
468    
# Line 531  C           Time dimension Line 471  C           Time dimension
471            ENDIF            ENDIF
472  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
473    
474             ENDDO
475  C--     end of Processing Fld # md  C--     end of Processing Fld # md
476          ENDIF          ENDIF
477        ENDDO        ENDDO
478    
479    #ifdef ALLOW_MDSIO
480          IF (diag_mdsio) THEN
481    C-    Note: temporary: since it's a pain to add more arguments to
482    C     all MDSIO S/R, uses instead this specific S/R to write only
483    C     meta files but with more informations in it.
484                nRec = nfields(listId)*averageCycle(listId)
485                CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
486         &              0, 0, nlevels(listId), ' ',
487         &              nfields(listId), flds(1,listId), 1, myTime,
488         &              nRec, myIter, myThid)
489          ENDIF
490    #endif /*  ALLOW_MDSIO  */
491    
492        RETURN        RETURN
493        END        END
494    
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.32
  ViewVC Help
Powered by ViewVC 1.1.22