/[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.27 by edhill, Mon Feb 6 21:20:23 2006 UTC
# Line 57  C     im    :: counter-mate pointer to s Line 57  C     im    :: counter-mate pointer to s
57        CHARACTER*8 parms1        CHARACTER*8 parms1
58        CHARACTER*3 mate_index        CHARACTER*3 mate_index
59        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+Nrphys,nSx,nSy)        _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)  
60        _RL undef, getcon        _RL undef, getcon
61        EXTERNAL getcon        EXTERNAL getcon
62        INTEGER ILNBLNK        INTEGER ILNBLNK
63        EXTERNAL ILNBLNK        EXTERNAL ILNBLNK
64        INTEGER ilen        INTEGER ilen
65        integer nlevsout,nplevs        INTEGER nlevsout
       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./  
66    
67        INTEGER ioUnit        INTEGER ioUnit
68        CHARACTER*(MAX_LEN_FNAM) fn        CHARACTER*(MAX_LEN_FNAM) fn
# Line 95  C     im    :: counter-mate pointer to s Line 72  C     im    :: counter-mate pointer to s
72  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
73        INTEGER ii        INTEGER ii
74        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
       CHARACTER*(5) ctmp  
75        INTEGER CW_DIMS, NLEN        INTEGER CW_DIMS, NLEN
76        PARAMETER ( CW_DIMS = 10 )        PARAMETER ( CW_DIMS = 10 )
77        PARAMETER ( NLEN    = 80 )        PARAMETER ( NLEN    = 80 )
# Line 103  C     im    :: counter-mate pointer to s Line 79  C     im    :: counter-mate pointer to s
79        CHARACTER*(NLEN) dn(CW_DIMS)        CHARACTER*(NLEN) dn(CW_DIMS)
80        CHARACTER*(NLEN) d_cw_name        CHARACTER*(NLEN) d_cw_name
81        CHARACTER*(NLEN) dn_blnk        CHARACTER*(NLEN) dn_blnk
82    #ifdef DIAG_MNC_COORD_NEEDSWORK
83          CHARACTER*(5) ctmp
84        _RS ztmp(Nr+Nrphys)        _RS ztmp(Nr+Nrphys)
85    #endif
86  #endif /*  ALLOW_MNC  */  #endif /*  ALLOW_MNC  */
87    
88  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
89    
90        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
91        undef = getcon('UNDEF')        undef = getcon('UNDEF')
       kappa = getcon('KAPPA')  
92        glf = globalFiles        glf = globalFiles
93        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
94        ilen = ILNBLNK(fnames(listId))        ilen = ILNBLNK(fnames(listId))
# Line 130  C       Update the record dimension by w Line 108  C       Update the record dimension by w
108          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)          CALL MNC_CW_SET_UDIM(diag_mnc_bn, -1, myThid)
109          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)
110          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)          CALL MNC_CW_SET_UDIM(diag_mnc_bn, 0, myThid)
111            CALL MNC_CW_I_W_S('I',diag_mnc_bn,0,0,'iter',myIter,myThid)
112    
113    C       NOTE: at some point it would be a good idea to add a time_bounds
114    C       variable that has dimension (2,T) and clearly denotes the
115    C       beginning and ending times for each diagnostics period
116    
117          dn(1)(1:NLEN) = dn_blnk(1:NLEN)          dn(1)(1:NLEN) = dn_blnk(1:NLEN)
118          WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)          WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)
# Line 158  C       grid.  As we start using diagnos Line 141  C       grid.  As we start using diagnos
141  C       levels, land levels, etc. the different vertical coordinate  C       levels, land levels, etc. the different vertical coordinate
142  C       dimensions will have to be taken into account.  C       dimensions will have to be taken into account.
143    
144    C       20051021 JMC & EH3 : We need to extend this so that a few
145    C       variables each defined on different grids do not have the same
146    C       vertical dimension names so we should be using a pattern such
147    C       as: Z[uml]td000000 where the 't' is the type as specified by
148    C       gdiag(10)
149    
150  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx  C       Now define:  Zmdxxxxxx, Zudxxxxxx, Zldxxxxxx
151          ctmp(1:5) = 'mul  '          ctmp(1:5) = 'mul  '
152          DO i = 1,3          DO i = 1,3
# Line 303  C-        end of empty diag / not empty Line 292  C-        end of empty diag / not empty
292            nlevsout = nlevels(listId)            nlevsout = nlevels(listId)
293    
294  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
295  C             Check to see if we need to interpolate before output  C         Check to see if we need to 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              CALL DIAGNOSTICS_INTERP_VERT(
300           if(fflags(listId)(2:2).eq.'P') then       I                     listId, md, ndId, ip, im,
301         U                     nlevsout,
302  c If nonlinear free surf is active, need averaged pressures       U                     qtmp1,
303  #ifdef NONLIN_FRSURF       I                     undef,
304            if(select_rStar.GT.0)then       I                     myTime, myIter, myThid )
305             call diagnostics_get_pointers('RSURF   ',ipoint1,jpoint1,           ENDIF
      .                                                           myThid)  
            call diagnostics_get_pointers('PRESSURE',ipoint2,jpoint2,  
      .                                                           myThid)  
 C if fizhi is being  used, may need to get physics grid pressures  
 #ifdef ALLOW_FIZHI  
            if(gdiag(ndId)(10:10) .EQ. 'L')then  
            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  
   
          endif  
306    
307  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
308  C         Prepare for mdsio optionality  C         Prepare for mdsio optionality
309            IF (diag_mdsio) THEN            IF (diag_mdsio) THEN
310              IF (fflags(listId)(1:1) .EQ. ' ') THEN              IF (fflags(listId)(1:1) .EQ. 'R') THEN
 C             This is the old default behavior  
               CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,'RL',  
      &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)  
             ELSEIF (fflags(listId)(1:1) .EQ. 'R') THEN  
311  C             Force it to be 32-bit precision  C             Force it to be 32-bit precision
312                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,'RL',                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,.FALSE.,
313       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)
314              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN
315  C             Force it to be 64-bit precision  C             Force it to be 64-bit precision
316                CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,'RL',                CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,.FALSE.,
317       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)
318                ELSE
319    C             This is the old default behavior
320                  CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,.FALSE.,
321         &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)
322              ENDIF              ENDIF
323            ENDIF            ENDIF
324  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
# Line 482  C           XY dimensions Line 370  C           XY dimensions
370              ENDIF              ENDIF
371                            
372  C           Z is special since it varies  C           Z is special since it varies
373              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevsout
374              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
375       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN
376                WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevsout
377              ENDIF              ENDIF
378              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
379       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN
380                WRITE(dn(3),'(a,i6.6)') 'Zld', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zld', nlevsout
381              ENDIF              ENDIF
382              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
383       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
384                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevsout
385              ENDIF              ENDIF
386              dim(3) = Nr+Nrphys              dim(3) = Nr+Nrphys
387              ib(3)  = 1              ib(3)  = 1
388              ie(3)  = nlevels(listId)              ie(3)  = nlevsout
389    
390  C           Time dimension  C           Time dimension
391              dn(4)(1:1) = 'T'              dn(4)(1:1) = 'T'
# Line 513  C           Time dimension Line 401  C           Time dimension
401       &             tdiag(ndId),myThid)       &             tdiag(ndId),myThid)
402              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',
403       &             udiag(ndId),myThid)       &             udiag(ndId),myThid)
404                CALL MNC_CW_ADD_VATTR_DBL(cdiag(ndId),'missing_value',
405         &             0.0 _d 0,myThid)
406    
407              IF ((fflags(listId)(1:1) .EQ. ' ')              IF ( ( (writeBinaryPrec .EQ. precFloat32)
408         &           .AND. (fflags(listId)(1:1) .NE. 'D')
409         &           .AND. (fflags(listId)(1:1) .NE. 'R') )
410       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN
411                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,
412       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
413              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              ELSEIF ( (writeBinaryPrec .EQ. precFloat64)
414         &             .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN
415                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
416       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
417              ENDIF              ENDIF
Legend:
Removed from v.1.19  
changed lines
  Added in v.1.27
  ViewVC Help
Powered by ViewVC 1.1.22