/[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.41 by jmc, Sun Jan 25 20:22:57 2009 UTC revision 1.49 by jmc, Mon Jun 6 15:42:58 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,
13       I     myIter,       I     myIter,
14       I     myTime,       I     myTime,
# Line 48  C     !FUNCTIONS: Line 48  C     !FUNCTIONS:
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)  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    C     nLevOutp :: number of levels to write in output file
59  C  C
60  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)  C--   COMMON /LOCAL_DIAGNOSTICS_OUT/ local common block (for multi-threaded)
61  C     qtmp1 :: thread-shared temporary array (needs to be in common block):  C     qtmp1 :: thread-shared temporary array (needs to be in common block):
# Line 62  C              diagnostic storage qdiag Line 64  C              diagnostic storage qdiag
64        COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1        COMMON /LOCAL_DIAGNOSTICS_OUT/ qtmp1
65        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)        _RL qtmp1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,NrMax,nSx,nSy)
66    
67        INTEGER i, j, k, lm, klev        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
# Line 70  C              diagnostic storage qdiag Line 72  C              diagnostic storage qdiag
72        _RL undef        _RL undef
73        _RL tmpLev        _RL tmpLev
74        INTEGER ilen        INTEGER ilen
75          INTEGER nLevOutp
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        INTEGER prec, nRec, nTimRec
82          _RL     timeRec(2)
83  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
84        LOGICAL glf        LOGICAL glf
85  #endif  #endif
86  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
87        INTEGER ll, llMx, jj, jjMx        INTEGER ll, llMx, jj, jjMx
88        INTEGER ii        INTEGER ii, klev
89        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn        CHARACTER*(MAX_LEN_FNAM) diag_mnc_bn
90        INTEGER CW_DIMS, NLEN        INTEGER CW_DIMS, NLEN
91        PARAMETER ( CW_DIMS = 10 )        PARAMETER ( CW_DIMS = 10 )
# Line 103  C              diagnostic storage qdiag Line 107  C              diagnostic storage qdiag
107    
108  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
109    
110    C---  set file properties
111        ioUnit= standardMessageUnit        ioUnit= standardMessageUnit
112        undef = UNSET_RL        undef = UNSET_RL
113  #ifdef ALLOW_FIZHI  #ifdef ALLOW_FIZHI
# Line 112  c     IF ( useFIZHI ) undef = getcon('UN Line 117  c     IF ( useFIZHI ) undef = getcon('UN
117        WRITE(suff,'(I10.10)') myIter        WRITE(suff,'(I10.10)') myIter
118        ilen = ILNBLNK(fnames(listId))        ilen = ILNBLNK(fnames(listId))
119        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)        WRITE( fn, '(A,A,A)' ) fnames(listId)(1:ilen),'.',suff(1:10)
120    C-    for now, if integrate vertically, output field has just 1 level:
121          nLevOutp = nlevels(listId)
122          IF ( fflags(listId)(2:2).EQ.'I' ) nLevOutp = 1
123    
124    C--   Set time information:
125          IF ( freq(listId).LT.0. ) THEN
126    C-    Snap-shot: store a unique time (which is consistent with State-Var timing)
127            nTimRec = 1
128            timeRec(1) = myTime
129          ELSE
130    C-    Time-average: store the 2 edges of the time-averaging interval.
131    C      this time is consitent with intermediate Var (i.e., non-state, e.g, flux,
132    C      tendencies) timing. For State-Var, this is shifted by + halt time-step.
133            nTimRec = 2
134    
135    C-    end of time-averaging interval:
136            timeRec(2) = myTime
137    
138    C-    begining of time-averaging interval:
139    c       timeRec(1) = myTime - freq(listId)
140    C     a) find the time of the previous multiple of output freq:
141            timeRec(1) = myTime-deltaTClock*0.5 _d 0
142            timeRec(1) = (timeRec(1)-phase(listId))/freq(listId)
143            i = INT( timeRec(1) )
144            IF ( timeRec(1).LT.0. ) THEN
145              tmpLev = FLOAT(i)
146              IF ( timeRec(1).NE.tmpLev ) i = i - 1
147            ENDIF
148            timeRec(1) = phase(listId) + freq(listId)*FLOAT(i)
149    c       if ( listId.eq.2 ) write(0,*) 'f',i,timeRec(1)/deltaTClock
150            timeRec(1) = MAX( timeRec(1), startTime )
151    
152    C     b) round off to nearest multiple of time-step:
153            timeRec(1) = (timeRec(1)-baseTime)/deltaTClock
154            i = NINT( timeRec(1) )
155    C     if just half way, NINT will return the next time-step: correct this
156            tmpLev = FLOAT(i) - 0.5 _d 0
157            IF ( timeRec(1).EQ.tmpLev ) i = i - 1
158            timeRec(1) = baseTime + deltaTClock*FLOAT(i)
159    c       if ( listId.eq.2 ) write(0,*) i,timeRec(1)/deltaTClock
160          ENDIF
161    C--   Convert time to iteration number (debug)
162    c     DO i=1,nTimRec
163    c       timeRec(i) = timeRec(i)/deltaTClock
164    c     ENDDO
165    
166  #ifdef ALLOW_MNC  #ifdef ALLOW_MNC
167  C-- this is a trick to reverse the order of the loops on md (= field)  C-- this is a trick to reverse the order of the loops on md (= field)
# Line 167  C       variable that has dimension (2,T Line 217  C       variable that has dimension (2,T
217  C       beginning and ending times for each diagnostics period  C       beginning and ending times for each diagnostics period
218    
219           dn(1)(1:NLEN) = dn_blnk(1:NLEN)           dn(1)(1:NLEN) = dn_blnk(1:NLEN)
220           WRITE(dn(1),'(a,i6.6)') 'Zmd', nlevels(listId)           WRITE(dn(1),'(a,i6.6)') 'Zmd', nLevOutp
221           dim(1) = nlevels(listId)           dim(1) = nLevOutp
222           ib(1)  = 1           ib(1)  = 1
223           ie(1)  = nlevels(listId)           ie(1)  = nLevOutp
224    
225           CALL MNC_CW_ADD_GNAME('diag_levels', 1,           CALL MNC_CW_ADD_GNAME('diag_levels', 1,
226       &        dim, dn, ib, ie, myThid)       &        dim, dn, ib, ie, myThid)
# Line 315  C-        Empty diagnostics case : Line 365  C-        Empty diagnostics case :
365              _END_MASTER( myThid )              _END_MASTER( myThid )
366              DO bj = myByLo(myThid), myByHi(myThid)              DO bj = myByLo(myThid), myByHi(myThid)
367                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
368                  DO k = 1,nlevels(listId)                  DO k = 1,nLevOutp
369                    DO j = 1-OLy,sNy+OLy                    DO j = 1-OLy,sNy+OLy
370                      DO i = 1-OLx,sNx+OLx                      DO i = 1-OLx,sNx+OLx
371                        qtmp1(i,j,k,bi,bj) = 0. _d 0                        qtmp1(i,j,k,bi,bj) = 0. _d 0
# Line 328  C-        Empty diagnostics case : Line 378  C-        Empty diagnostics case :
378            ELSE            ELSE
379  C-        diagnostics is not empty :  C-        diagnostics is not empty :
380    
381              IF ( debugLevel.GE.debLevA .AND. myThid.EQ.1 ) THEN              IF ( debugLevel.GE.debLevB .AND. myThid.EQ.1 ) THEN
382                WRITE(ioUnit,'(A,I6,3A,I8,2A)')                WRITE(ioUnit,'(A,I6,3A,I8,2A)')
383       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),       &         ' Computing Diagnostic # ', ndId, '  ', cdiag(ndId),
384       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)       &         '     Counter:',ndiag(ip,1,1),'   Parms: ',gdiag(ndId)
# Line 351  C-        diagnostics is not empty : Line 401  C-        diagnostics is not empty :
401                ENDIF                ENDIF
402              ENDIF              ENDIF
403    
404              IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).NE.' ' ) THEN
405  C-       get all the levels (for vertical interpolation)  C-       get all the levels (for vertical post-processing)
406                DO bj = myByLo(myThid), myByHi(myThid)                DO bj = myByLo(myThid), myByHi(myThid)
407                 DO bi = myBxLo(myThid), myBxHi(myThid)                 DO bi = myBxLo(myThid), myBxHi(myThid)
408                  DO k = 1,kdiag(ndId)                  DO k = 1,kdiag(ndId)
# Line 378  C-       get only selected levels: Line 428  C-       get only selected levels:
428                ENDDO                ENDDO
429              ENDIF              ENDIF
430    
 C-        end of empty diag / not empty block  
           ENDIF  
   
431  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
432  C         Check to see if we need to interpolate before output  C--     Apply specific post-processing (e.g., interpolate) before output
433  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
434            IF ( fflags(listId)(2:2).EQ.'P' ) THEN              IF ( fflags(listId)(2:2).EQ.'P' ) THEN
435  C-        Do vertical interpolation:  C-          Do vertical interpolation:
436             IF ( fluidIsAir ) THEN               IF ( fluidIsAir ) THEN
437  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);  C jmc: for now, this can only work in an atmospheric set-up (fluidIsAir);
438              CALL DIAGNOSTICS_INTERP_VERT(                CALL DIAGNOSTICS_INTERP_VERT(
439       I                     listId, md, ndId, ip, im, lm,       I                         listId, md, ndId, ip, im, lm,
440       U                     qtmp1,       U                         qtmp1,
441       I                     undef, myTime, myIter, myThid )       I                         undef, myTime, myIter, myThid )
442             ELSE               ELSE
443               WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',                 WRITE(msgBuf,'(2A)') 'DIAGNOSTICS_OUT: ',
444       &         'INTERP_VERT not allowed in this config'       &           'INTERP_VERT not allowed in this config'
445               CALL PRINT_ERROR( msgBuf , myThid )                 CALL PRINT_ERROR( msgBuf , myThid )
446               STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'                 STOP 'ABNORMAL END: S/R DIAGNOSTICS_OUT'
447             ENDIF               ENDIF
448                ENDIF
449                IF ( fflags(listId)(2:2).EQ.'I' ) THEN
450    C-          Integrate vertically: for now, output field has just 1 level:
451                  CALL DIAGNOSTICS_SUM_LEVELS(
452         I                         listId, md, ndId, ip, im, lm,
453         U                         qtmp1,
454         I                         undef, myTime, myIter, myThid )
455                ENDIF
456    
457    C--     End of empty diag / not-empty diag block
458            ENDIF            ENDIF
459    
460  C--    Ready to write field "md", element "lm" in averageCycle(listId)  C--     Ready to write field "md", element "lm" in averageCycle(listId)
461    
462  C-        write to binary file, using MDSIO pkg:  C-        write to binary file, using MDSIO pkg:
463            IF ( diag_mdsio ) THEN            IF ( diag_mdsio ) THEN
# Line 413  C           fFlag(1)=R(or D): force it t Line 470  C           fFlag(1)=R(or D): force it t
470  C         a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R  C         a hack not to write meta files now: pass -nRec < 0 to MDS_WRITE S/R
471              CALL WRITE_REC_LEV_RL(              CALL WRITE_REC_LEV_RL(
472       I                            fn, prec,       I                            fn, prec,
473       I                            NrMax, 1, nlevels(listId),       I                            NrMax, 1, nLevOutp,
474       I                            qtmp1, -nRec, myIter, myThid )       I                            qtmp1, -nRec, myIter, myThid )
475            ENDIF            ENDIF
476    
# Line 464  C           XY dimensions Line 521  C           XY dimensions
521              ENDIF              ENDIF
522    
523  C           Z is special since it varies  C           Z is special since it varies
524              WRITE(dn(3),'(a,i6.6)') 'Zd', nlevels(listId)              WRITE(dn(3),'(a,i6.6)') 'Zd', nLevOutp
525              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
526       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'M') ) THEN
527                WRITE(dn(3),'(a,i6.6)') 'Zmd', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zmd', nLevOutp
528              ENDIF              ENDIF
529              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
530       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'L') ) THEN
531                WRITE(dn(3),'(a,i6.6)') 'Zld', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zld', nLevOutp
532              ENDIF              ENDIF
533              IF ( (gdiag(ndId)(10:10) .EQ. 'R')              IF ( (gdiag(ndId)(10:10) .EQ. 'R')
534       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN       &           .AND. (gdiag(ndId)(9:9) .EQ. 'U') ) THEN
535                WRITE(dn(3),'(a,i6.6)') 'Zud', nlevels(listId)                WRITE(dn(3),'(a,i6.6)') 'Zud', nLevOutp
536              ENDIF              ENDIF
537              dim(3) = NrMax              dim(3) = NrMax
538              ib(3)  = 1              ib(3)  = 1
539              ie(3)  = nlevels(listId)              ie(3)  = nLevOutp
540    
541  C           Time dimension  C           Time dimension
542              dn(4)(1:1) = 'T'              dn(4)(1:1) = 'T'
# Line 505  C     assign missing values and set flag Line 562  C     assign missing values and set flag
562       I            misval_r8, misval_r4, misval_int,       I            misval_r8, misval_r4, misval_int,
563       I            myThid )       I            myThid )
564  C     and now use the missing values for masking out the land points  C     and now use the missing values for masking out the land points
565    C     note: better to use 2-D mask if kdiag <> Nr or vert.integral
566               DO bj = myByLo(myThid), myByHi(myThid)               DO bj = myByLo(myThid), myByHi(myThid)
567                DO bi = myBxLo(myThid), myBxHi(myThid)                DO bi = myBxLo(myThid), myBxHi(myThid)
568                 DO k = 1,nlevels(listId)                 DO k = 1,nLevOutp
569                  klev = NINT(levs(k,listId))                  klev = NINT(levs(k,listId))
570                    IF ( fflags(listId)(2:2).EQ.'I' ) kLev = 1
571                  DO j = 1-OLy,sNy+OLy                  DO j = 1-OLy,sNy+OLy
572                   DO i = 1-OLx,sNx+OLx                   DO i = 1-OLx,sNx+OLx
573                    IF ( _hFacC(I,J,klev,bi,bj) .EQ. 0. )                    IF ( maskC(i,j,klev,bi,bj) .EQ. 0. )
574       &                 qtmp1(i,j,k,bi,bj) = misvalLoc       &                 qtmp1(i,j,k,bi,bj) = misvalLoc
575                   ENDDO                   ENDDO
576                  ENDDO                  ENDDO
# Line 562  C--   end loop on jj counter Line 621  C--   end loop on jj counter
621    
622  #ifdef ALLOW_MDSIO  #ifdef ALLOW_MDSIO
623        IF (diag_mdsio) THEN        IF (diag_mdsio) THEN
624  C-    Note: temporary: since it's a pain to add more arguments to  C-    Note: temporary: since it is a pain to add more arguments to
625  C     all MDSIO S/R, uses instead this specific S/R to write only  C     all MDSIO S/R, uses instead this specific S/R to write only
626  C     meta files but with more informations in it.  C     meta files but with more informations in it.
627              glf = globalFiles              glf = globalFiles
628              nRec = nfields(listId)*averageCycle(listId)              nRec = nfields(listId)*averageCycle(listId)
629              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,              CALL MDS_WR_METAFILES(fn, prec, glf, .FALSE.,
630       &              0, 0, nlevels(listId), ' ',       &              0, 0, nLevOutp, ' ',
631       &              nfields(listId), flds(1,listId), 1, myTime,       &              nfields(listId), flds(1,listId), nTimRec, timeRec,
632       &              nRec, myIter, myThid)       &              nRec, myIter, myThid)
633        ENDIF        ENDIF
634  #endif /*  ALLOW_MDSIO  */  #endif /*  ALLOW_MDSIO  */
Legend:
Removed from v.1.41  
changed lines
  Added in v.1.49
  ViewVC Help
Powered by ViewVC 1.1.22