/[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.17 by molod, Tue Jul 12 20:25:45 2005 UTC revision 1.34 by jmc, Tue Nov 13 19:43:44 2007 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 45  CEOP Line 40  CEOP
40    
41  C     !LOCAL VARIABLES:  C     !LOCAL VARIABLES:
42  C     i,j,k :: loop indices  C     i,j,k :: loop indices
43    C     lm    :: loop index (averageCycle)
44  C     md    :: field number in the list "listId".  C     md    :: field number in the list "listId".
45  C     ndId  :: diagnostics  Id number (in available diagnostics list)  C     ndId  :: diagnostics  Id number (in available diagnostics list)
46  C     mate  :: counter mate Id number (in available diagnostics list)  C     mate  :: counter mate Id number (in available diagnostics list)
47  C     ip    :: diagnostics  pointer to storage array  C     ip    :: diagnostics  pointer to storage array
48  C     im    :: counter-mate pointer to storage array  C     im    :: counter-mate pointer to storage array
49        INTEGER i, j, k  C
50    C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
51    C     qtmp1 :: thread-shared temporary array (needs to be in common block):
52    C              to write a diagnostic field to file, copy it first from (big)
53    C              diagnostic storage qdiag into it.
54          COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
55          _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
56    
57          INTEGER i, j, k, lm
58        INTEGER bi, bj        INTEGER bi, bj
59        INTEGER md, ndId, ip, im        INTEGER md, ndId, ip, im
60        INTEGER mate, mVec        INTEGER mate, mVec
61        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)  
62        _RL undef, getcon        _RL undef, getcon
63          _RL tmpLev
64        EXTERNAL getcon        EXTERNAL getcon
65        INTEGER ILNBLNK        INTEGER ILNBLNK
66        EXTERNAL ILNBLNK        EXTERNAL ILNBLNK
67        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./  
68    
69        INTEGER ioUnit        INTEGER ioUnit
70        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
71        CHARACTER*(MAX_LEN_MBUF) suff        CHARACTER*(MAX_LEN_MBUF) suff
72        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
73          INTEGER prec, nRec
74    #ifdef ALLOW_MDSIO
75        LOGICAL glf        LOGICAL glf
76    #endif
77  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
78        INTEGER ii        INTEGER ii
79        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
       CHARACTER*(5) ctmp  
80        INTEGER CW_DIMS, NLEN        INTEGER CW_DIMS, NLEN
81        PARAMETER ( CW_DIMS = 10 )        PARAMETER ( CW_DIMS = 10 )
82        PARAMETER ( NLEN    = 80 )        PARAMETER ( NLEN    = 80 )
# Line 103  C     im    :: counter-mate pointer to s Line 84  C     im    :: counter-mate pointer to s
84        CHARACTER*(NLEN) dn(CW_DIMS)        CHARACTER*(NLEN) dn(CW_DIMS)
85        CHARACTER*(NLEN) d_cw_name        CHARACTER*(NLEN) d_cw_name
86        CHARACTER*(NLEN) dn_blnk        CHARACTER*(NLEN) dn_blnk
87        _RS ztmp(Nr+Nrphys)  #ifdef DIAG_MNC_COORD_NEEDSWORK
88          CHARACTER*(5) ctmp
89          _RS ztmp(NrMax)
90    #endif
91  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
92    
93  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
94    
95        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
96        undef = getcon('UNDEF')        undef = getcon('UNDEF')
       kappa = getcon('KAPPA')  
       glf = globalFiles  
97        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
98        ilen = ILNBLNK(fnames(listId))        ilen = ILNBLNK(fnames(listId))
99        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)
100    
 C Initialize the qtmp1 array to all undefs -- need this for unfilled levels  
       DO bj = myByLo(myThid), myByHi(myThid)  
         DO bi = myBxLo(myThid), myBxHi(myThid)  
           DO k = 1,nlevels(listId)  
             DO j = 1-OLy,sNy+OLy  
               DO i = 1-OLx,sNx+OLx  
                 qtmp1(i,j,k,bi,bj) = undef  
               ENDDO  
             ENDDO  
           ENDDO  
         ENDDO  
       ENDDO  
   
101  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
102        IF (useMNC .AND. diag_mnc) THEN        IF (useMNC .AND. diag_mnc) THEN
103          DO i = 1,MAX_LEN_FNAM          DO i = 1,MAX_LEN_FNAM
# Line 143  C       Update the record dimension by w Line 112  C       Update the record dimension by w
112          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)
113          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)
114          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)
115            CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',myIter,myThid)
116    
117    C       NOTE: at some point it would be a good idea to add a time_bounds
118    C       variable that has dimension (2,T) and clearly denotes the
119    C       beginning and ending times for each diagnostics period
120    
121          dn(1)(1:NLEN) = dn_blnk(1:NLEN)          dn(1)(1:NLEN) = dn_blnk(1:NLEN)
122          WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)          WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)
# Line 150  C       Update the record dimension by w Line 124  C       Update the record dimension by w
124          ib(1)  = 1          ib(1)  = 1
125          ie(1)  = nlevels(listId)          ie(1)  = nlevels(listId)
126    
127          CALL MNC_CW_ADD_GNAME('diag_levels', 1,          CALL MNC_CW_ADD_GNAME('diag_levels', 1,
128       &       dim, dn, ib, ie, myThid)       &       dim, dn, ib, ie, myThid)
129          CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',          CALL MNC_CW_ADD_VNAME('diag_levels', 'diag_levels',
130       &       0,0, myThid)       &       0,0, myThid)
131          CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',          CALL MNC_CW_ADD_VATTR_TEXT('diag_levels','description',
132       &       'Idicies of vertical levels within the source arrays',       &       'Idicies of vertical levels within the source arrays',
133       &       myThid)       &       myThid)
134            
135          CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,          CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
136       &       'diag_levels', levs(1,listId), myThid)       &       'diag_levels', levs(1,listId), myThid)
137    
# Line 171  C       grid.  As we start using diagnos Line 145  C       grid.  As we start using diagnos
145  C       levels, land levels, etc. the different vertical coordinate  C       levels, land levels, etc. the different vertical coordinate
146  C       dimensions will have to be taken into account.  C       dimensions will have to be taken into account.
147    
148    C       20051021 JMC & EH3 : We need to extend this so that a few
149    C       variables each defined on different grids do not have the same
150    C       vertical dimension names so we should be using a pattern such
151    C       as: Z[uml]td000000 where the 't' is the type as specified by
152    C       gdiag(10)
153    
154  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx
155          ctmp(1:5) = 'mul  '          ctmp(1:5) = 'mul  '
156          DO i = 1,3          DO i = 1,3
# Line 182  C       Now define:  Zmdxxxxxx, Zudxxxxx Line 162  C       Now define:  Zmdxxxxxx, Zudxxxxx
162  C         The following three ztmp() loops should eventually be modified  C         The following three ztmp() loops should eventually be modified
163  C         to reflect the fractional nature of levs(j,l) -- they should  C         to reflect the fractional nature of levs(j,l) -- they should
164  C         do something like:  C         do something like:
165  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))  C            ztmp(j) = rC(INT(FLOOR(levs(j,l))))
166  C                      + ( rC(INT(FLOOR(levs(j,l))))  C                      + ( rC(INT(FLOOR(levs(j,l))))
167  C                          + rC(INT(CEIL(levs(j,l)))) )  C                          + rC(INT(CEIL(levs(j,l)))) )
168  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )  C                        / ( levs(j,l) - FLOOR(levs(j,l)) )
169  C         for averaged levels.  C         for averaged levels.
# Line 218  C         for averaged levels. Line 198  C         for averaged levels.
198        ENDIF        ENDIF
199  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
200    
201    C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
202    
203        DO md = 1,nfields(listId)        DO md = 1,nfields(listId)
204          ndId = jdiag(md,listId)          ndId = jdiag(md,listId)
205          parms1 = gdiag(ndId)(1:8)          parms1 = gdiag(ndId)(1:8)
206            mate = 0
207            mVec = 0
208            IF ( parms1(5:5).EQ.'C' ) THEN
209    C-      Check for Mate of a Counter Diagnostic
210               READ(parms1,'(5X,I3)') mate
211            ELSEIF ( parms1(1:1).EQ.'U' .OR. parms1(1:1).EQ.'V' ) THEN
212    C-      Check for Mate of a Vector Diagnostic
213               READ(parms1,'(5X,I3)') mVec
214            ENDIF
215          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
216  C--     Start processing 1 Fld :  C--     Start processing 1 Fld :
217             DO lm=1,averageCycle(listId)
218    
219            ip = ABS(idiag(md,listId))            ip = ABS(idiag(md,listId)) + kdiag(ndId)*(lm-1)
220            im = mdiag(md,listId)            im = mdiag(md,listId)
221              IF (mate.GT.0) im = im + kdiag(mate)*(lm-1)
222              IF (mVec.GT.0) im = im + kdiag(mVec)*(lm-1)
223    
224            IF ( ndiag(ip,1,1).EQ.0 ) THEN            IF ( ndiag(ip,1,1).EQ.0 ) THEN
225  C-        Empty diagnostics case :  C-        Empty diagnostics case :
226    
# Line 239  C-        Empty diagnostics case : Line 234  C-        Empty diagnostics case :
234       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)       &       ' (#',md,' ) in outp.Stream: ',fnames(listId)
235              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
236       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
237              WRITE(msgBuf,'(A,I2,A)')              IF ( averageCycle(listId).GT.1 ) THEN
238       &       '- WARNING -   has not been filled (ndiag=',               WRITE(msgBuf,'(A,2(I2,A))')
239       &       ndiag(ip,1,1), ' )'       &        '- WARNING -   has not been filled (ndiag(lm=',lm,')=',
240         &                                            ndiag(ip,1,1), ' )'
241                ELSE
242                 WRITE(msgBuf,'(A,2(I2,A))')
243         &        '- WARNING -   has not been filled (ndiag=',
244         &                                            ndiag(ip,1,1), ' )'
245                ENDIF
246              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,              CALL PRINT_MESSAGE( msgBuf, errorMessageUnit,
247       &                          SQUEEZE_RIGHT, myThid)       &                          SQUEEZE_RIGHT, myThid)
248              WRITE(msgBuf,'(A)')              WRITE(msgBuf,'(A)')
# Line 264  C-        Empty diagnostics case : Line 265  C-        Empty diagnostics case :
265            ELSE            ELSE
266  C-        diagnostics is not empty :  C-        diagnostics is not empty :
267    
268              IF ( myThid.EQ.1 ) WRITE(ioUnit,'(A,I3,3A,I8,2A)')              IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN
269                  WRITE(ioUnit,'(A,I3,3A,I8,2A)')
270       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
271       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
272                  IF ( mate.GT.0 ) THEN
273              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)')  
274       &         '       use Counter Mate for  ', cdiag(ndId),       &         '       use Counter Mate for  ', cdiag(ndId),
275       &         '     Diagnostic # ',mate, '  ', cdiag(mate)       &         '     Diagnostic # ',mate, '  ', cdiag(mate)
276                  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  
277                  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
278                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I3,3A)')
279       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
280       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
281       &             ' exists '       &             ' exists '
282                  ELSE                  ELSE
283                   IF ( myThid.EQ.1 ) WRITE(ioUnit,'(3A,I3,3A)')                   WRITE(ioUnit,'(3A,I3,3A)')
284       &             '           Vector  Mate for  ', cdiag(ndId),       &             '           Vector  Mate for  ', cdiag(ndId),
285       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),       &             '     Diagnostic # ',mVec, '  ', cdiag(mVec),
286       &             ' not enabled'       &             ' not enabled'
# Line 299  C             -------------------------- Line 288  C             --------------------------
288                ENDIF                ENDIF
289              ENDIF              ENDIF
290    
291              DO bj = myByLo(myThid), myByHi(myThid)              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
292               DO bi = myBxLo(myThid), myBxHi(myThid)  C-       get all the levels (for vertical interpolation)
293                DO k = 1,nlevels(listId)                DO bj = myByLo(myThid), myByHi(myThid)
294                  CALL GETDIAG(                 DO bi = myBxLo(myThid), myBxHi(myThid)
295       I                       levs(k,listId),undef,                  DO k = 1,kdiag(ndId)
296       O                       qtmp1(1-OLx,1-OLy,k,bi,bj),                    tmpLev = k
297       I                       ndId,mate,ip,im,bi,bj,myThid)                    CALL GETDIAG(
298         I                         tmpLev,undef,
299         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
300         I                         ndId,mate,ip,im,bi,bj,myThid)
301                    ENDDO
302                   ENDDO
303                ENDDO                ENDDO
304               ENDDO              ELSE
305              ENDDO  C-       get only selected levels:
306                  DO bj = myByLo(myThid), myByHi(myThid)
307                   DO bi = myBxLo(myThid), myBxHi(myThid)
308                    DO k = 1,nlevels(listId)
309                      CALL GETDIAG(
310         I                         levs(k,listId),undef,
311         O                         qtmp1(1-OLx,1-OLy,k,bi,bj),
312         I                         ndId,mate,ip,im,bi,bj,myThid)
313                    ENDDO
314                   ENDDO
315                  ENDDO
316                ENDIF
317    
318  C-        end of empty diag / not empty block  C-        end of empty diag / not empty block
319            ENDIF            ENDIF
320    
           nlevsout = nlevels(listId)  
   
321  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
322  C             Check to see if we need to interpolate before output  C         Check to see if we need to interpolate before output
323  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
324  C  (we are still inside field exist if sequence and field do loop)            IF ( fflags(listId)(2:2).EQ.'P' ) THEN
325  C  C-        Do vertical interpolation:
326               IF ( fluidIsAir ) THEN
327           if(fflags(listId)(2:2).eq.'P') then  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
328                CALL DIAGNOSTICS_INTERP_VERT(
329  c If nonlinear free surf is active, need averaged pressures       I                     listId, md, ndId, ip, im, lm,
330  #ifdef NONLIN_FRSURF       U                     qtmp1,
331            if(select_rStar.GT.0)then       I                     undef, myTime, myIter, myThid )
332             call diagnostics_get_pointers('RSURF   ',ipoint1,jpoint1,             ELSE
333       .                                                           myThid)               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
334             call diagnostics_get_pointers('PRESSURE',ipoint2,jpoint2,       &         'INTERP_VERT not allowed in this config'
335       .                                                           myThid)               CALL PRINT_ERROR( msgBuf , myThid )
336  C if fizhi is being  used, may need to get physics grid pressures               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
337  #ifdef ALLOW_FIZHI             ENDIF
338             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  
339    
340           endif  C--    Ready to write field "md", element "lm" in averageCycle(listId)
341    
342  #ifdef ALLOW_MDSIO  C-        write to binary file, using MDSIO pkg:
343  C         Prepare for mdsio optionality            IF ( diag_mdsio ) THEN
344            IF (diag_mdsio) THEN              nRec = lm + (md-1)*averageCycle(listId)
345              IF (fflags(listId)(1:1) .EQ. ' ') THEN  C           default precision for output files
346  C             This is the old default behavior              prec = writeBinaryPrec
347                CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,'RL',  C           fFlag(1)=R(or D): force it to be 32-bit(or 64) precision
348       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)              IF ( fflags(listId)(1:1).EQ.'R' ) prec = precFloat32
349              ELSEIF (fflags(listId)(1:1) .EQ. 'R') THEN              IF ( fflags(listId)(1:1).EQ.'D' ) prec = precFloat64
350  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
351                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,'RL',              CALL WRITE_REC_LEV_RL(
352       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       I                            fn, prec,
353              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN       I                            NrMax, 1, nlevels(listId),
354  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  
355            ENDIF            ENDIF
 #endif /*  ALLOW_MDSIO  */  
356    
357  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
358            IF (useMNC .AND. diag_mnc) THEN            IF (useMNC .AND. diag_mnc) THEN
# Line 472  C           XY dimensions Line 378  C           XY dimensions
378              dim(2)       = sNy + 2*OLy              dim(2)       = sNy + 2*OLy
379              ib(1)        = OLx + 1              ib(1)        = OLx + 1
380              ib(2)        = OLy + 1              ib(2)        = OLy + 1
381              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN              IF (gdiag(ndId)(2:2) .EQ. 'M') THEN
382                dn(1)(1:2) = 'X'                dn(1)(1:2) = 'X'
383                ie(1)      = OLx + sNx                ie(1)      = OLx + sNx
384                dn(2)(1:2) = 'Y'                dn(2)(1:2) = 'Y'
# Line 493  C           XY dimensions Line 399  C           XY dimensions
399                dn(2)(1:3) = 'Yp1'                dn(2)(1:3) = 'Yp1'
400                ie(2)      = OLy + sNy + 1                ie(2)      = OLy + sNy + 1
401              ENDIF              ENDIF
402                
403  C           Z is special since it varies  C           Z is special since it varies
404              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)
405              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
# Line 508  C           Z is special since it varies Line 414  C           Z is special since it varies
414       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
415                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)
416              ENDIF              ENDIF
417              dim(3) = Nr+Nrphys              dim(3) = NrMax
418              ib(3)  = 1              ib(3)  = 1
419              ie(3)  = nlevels(listId)              ie(3)  = nlevels(listId)
420    
# Line 518  C           Time dimension Line 424  C           Time dimension
424              ib(4)  = 1              ib(4)  = 1
425              ie(4)  = 1              ie(4)  = 1
426    
427              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,              CALL MNC_CW_ADD_GNAME(d_cw_name, 4,
428       &             dim, dn, ib, ie, myThid)       &             dim, dn, ib, ie, myThid)
429              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,              CALL MNC_CW_ADD_VNAME(cdiag(ndId), d_cw_name,
430       &             4,5, myThid)       &             4,5, myThid)
431              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'description',
432       &             tdiag(ndId),myThid)       &             tdiag(ndId),myThid)
433              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',
434       &             udiag(ndId),myThid)       &             udiag(ndId),myThid)
435    
436              IF ((fflags(listId)(1:1) .EQ. ' ')  C           Per the observations of Baylor, this has been commented out
437    C           until we have code that can write missing_value attributes
438    C           in a way thats compatible with most of the more popular
439    C           netCDF tools including ferret.  Using all-zeros completely
440    C           breaks ferret.
441    
442    C           CALL MNC_CW_ADD_VATTR_DBL(cdiag(ndId),'missing_value',
443    C           &             0.0 _d 0,myThid)
444    
445                IF ( ( (writeBinaryPrec .EQ. precFloat32)
446         &           .AND. (fflags(listId)(1:1) .NE. 'D')
447         &           .AND. (fflags(listId)(1:1) .NE. 'R') )
448       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN
449                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,
450       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
451              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              ELSEIF ( (writeBinaryPrec .EQ. precFloat64)
452         &             .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN
453                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
454       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
455              ENDIF              ENDIF
456                
457              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)              CALL MNC_CW_DEL_VNAME(cdiag(ndId), myThid)
458              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)              CALL MNC_CW_DEL_GNAME(d_cw_name, myThid)
459    
# Line 544  C           Time dimension Line 462  C           Time dimension
462            ENDIF            ENDIF
463  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
464    
465             ENDDO
466  C--     end of Processing Fld # md  C--     end of Processing Fld # md
467          ENDIF          ENDIF
468        ENDDO        ENDDO
469    
470    #ifdef ALLOW_MDSIO
471          IF (diag_mdsio) THEN
472    C-    Note: temporary: since it's a pain to add more arguments to
473    C     all MDSIO S/R, uses instead this specific S/R to write only
474    C     meta files but with more informations in it.
475                glf = globalFiles
476                nRec = nfields(listId)*averageCycle(listId)
477                CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
478         &              0, 0, nlevels(listId), ' ',
479         &              nfields(listId), flds(1,listId), 1, myTime,
480         &              nRec, myIter, myThid)
481          ENDIF
482    #endif /*  ALLOW_MDSIO  */
483    
484        RETURN        RETURN
485        END        END
486    
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.34
  ViewVC Help
Powered by ViewVC 1.1.22