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