/[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.24 by molod, Fri Nov 18 23:46:39 2005 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))
95        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)
96    
 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  
   
97  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
98        IF (useMNC .AND. diag_mnc) THEN        IF (useMNC .AND. diag_mnc) THEN
99          DO i = 1,MAX_LEN_FNAM          DO i = 1,MAX_LEN_FNAM
# Line 143  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,1,1,'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 316  C-        end of empty diag / not empty Line 286  C-        end of empty diag / not empty
286            nlevsout = nlevels(listId)            nlevsout = nlevels(listId)
287    
288  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
289  C             Check to see if we need to interpolate before output  C         Check to see if we need to interpolate before output
290  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
291  C  (we are still inside field exist if sequence and field do loop)           IF ( fflags(listId)(2:2).EQ.'P' ) THEN
292  C  C-        Do vertical interpolation:
293              CALL DIAGNOSTICS_INTERP_VERT(
294           if(fflags(listId)(2:2).eq.'P') then       I                     listId, md, ndId, ip, im,
295         U                     nlevsout,
296  c If nonlinear free surf is active, need averaged pressures       U                     qtmp1,
297  #ifdef NONLIN_FRSURF       I                     undef,
298            if(select_rStar.GT.0)then       I                     myTime, myIter, myThid )
299             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  
300    
301  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
302  C         Prepare for mdsio optionality  C         Prepare for mdsio optionality
303            IF (diag_mdsio) THEN            IF (diag_mdsio) THEN
304              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  
305  C             Force it to be 32-bit precision  C             Force it to be 32-bit precision
306                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,'RL',                CALL MDSWRITEFIELD_NEW(fn,precFloat32,glf,.FALSE.,
307       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)
308              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN
309  C             Force it to be 64-bit precision  C             Force it to be 64-bit precision
310                CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,'RL',                CALL MDSWRITEFIELD_NEW(fn,precFloat64,glf,.FALSE.,
311       &             Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)       &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)
312                ELSE
313    C             This is the old default behavior
314                  CALL MDSWRITEFIELD_NEW(fn,writeBinaryPrec,glf,.FALSE.,
315         &             'RL',Nr+Nrphys,nlevsout,qtmp1,md,myIter,myThid)
316              ENDIF              ENDIF
317            ENDIF            ENDIF
318  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
# Line 495  C           XY dimensions Line 364  C           XY dimensions
364              ENDIF              ENDIF
365                            
366  C           Z is special since it varies  C           Z is special since it varies
367              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevsout
368              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
369       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN
370                WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevsout
371              ENDIF              ENDIF
372              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
373       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN
374                WRITE(dn(3),'(a,i6.6)') 'Zld', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zld', nlevsout
375              ENDIF              ENDIF
376              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
377       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
378                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevsout
379              ENDIF              ENDIF
380              dim(3) = Nr+Nrphys              dim(3) = Nr+Nrphys
381              ib(3)  = 1              ib(3)  = 1
382              ie(3)  = nlevels(listId)              ie(3)  = nlevsout
383    
384  C           Time dimension  C           Time dimension
385              dn(4)(1:1) = 'T'              dn(4)(1:1) = 'T'
# Line 527  C           Time dimension Line 396  C           Time dimension
396              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',              CALL MNC_CW_ADD_VATTR_TEXT(cdiag(ndId),'units',
397       &             udiag(ndId),myThid)       &             udiag(ndId),myThid)
398    
399              IF ((fflags(listId)(1:1) .EQ. ' ')              IF ( ( (writeBinaryPrec .EQ. precFloat32)
400         &           .AND. (fflags(listId)(1:1) .NE. 'D')
401         &           .AND. (fflags(listId)(1:1) .NE. 'R') )
402       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN       &           .OR. (fflags(listId)(1:1) .EQ. 'R')) THEN
403                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('R',diag_mnc_bn,0,0,
404       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
405              ELSEIF (fflags(listId)(1:1) .EQ. 'D') THEN              ELSEIF ( (writeBinaryPrec .EQ. precFloat64)
406         &             .OR. (fflags(listId)(1:1) .EQ. 'D') ) THEN
407                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,                CALL MNC_CW_RL_W('D',diag_mnc_bn,0,0,
408       &             cdiag(ndId), qtmp1, myThid)       &             cdiag(ndId), qtmp1, myThid)
409              ENDIF              ENDIF
Legend:
Removed from v.1.17  
changed lines
  Added in v.1.24
  ViewVC Help
Powered by ViewVC 1.1.22