      PROGRAM TIMAVG                                                                                                                
C     PROGRAM TIMAVG (X1,       AVG,       DEV,       VAR,                      C2                                                  
C    1                X2,       ...,       X12,                OUTPUT,  )       C2                                                  
C    2         TAPE11=X1, TAPE2=AVG, TAPE3=DEV, TAPE4=VAR,                      C2                                                  
C    3         TAPE12=X1,       ...,TAPE22=X12,           TAPE6=OUTPUT)         C2                                                  
C     ----------------------------------------------------------------          C2                                                  
C                                                                               C2                                                  
C     NOV 24/09 - S.KHARIN (ADD UP TO 12 OPTIONAL ADDITIONAL INPUT FILES)       C2                                                  
C     DEC 12/06 - F.MAJAESS (BASE PACKING DENSITY ON OVERALL SETTING INSTEAD                                                        
C                            OF JUST THAT IN THE FIRST RECORD)                                                                      
C     JUN 24/05 - S.KHARIN (FORCE AN ABORT ON "TIME" RECORD KIND)                                                                   
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)                                                           
C     NOV 17/00 - S.KHARIN (FIX THE NEGATIVE VAR BUG. ABORT IF INPUT FILE                                                           
C                           CONTAINS DIFFERENT VARIABLE NAMES.)                                                                     
C     APR 11/00 - S.KHARIN (SIMPLIFY CODE. RELAX RESTRICTION ON LEVEL NUMBER.                                                       
C                           OPTIONALLY CALCULATE DEVIATIONS AND VARIANCE.)                                                          
C     NOV 06/95 - F.MAJAESS (REVISE DEFAULT PACKING DENSITY VALUE)                                                                  
C     JAN 29/92 - E. CHAN  (CONVERT HOLLERITH LITERALS TO ASCII, REPLACE                                                            
C                           BUFFERED I/O, REPLACE CALL ASSIGN WITH OPEN,                                                            
C                           AND REPLACE CALL RELEASE WITH A CALL SYSTEM                                                             
C                           TO REMOVE THE FILE)                                                                                     
C     MAY 24/88 - F.MAJAESS (INCREASE THE DIMENSION TO 75000; T30-16 LEVELS)                                                        
C     MAR 25/88 - F.MAJAESS (AVOID USING A SCRATCH DISK FILES IF POSSIBLE                                                           
C                            AND DO ONE LEVEL COMPUTATION AT A TIME)                                                                
C     MAR 11/88 - F.MAJAESS (COSMETIC MODIFICATIONS)                                                                                
C     MAR 10/88 - F.MAJAESS (CORRECT HANDLING OF AN EMPTY FILE "SERA")                                                              
C     FEB 20/85 - B.DUGAS. (RESTRUCTURE TO USE SCRATCH DISK FILE)                                                                   
C     MAY 07/80 - J.D.HENDERSON                                                                                                     
C                                                                               C2                                                  
CTIMAVG  - COMPUTES MULTI-LEVEL TIME AVERAGE, DEVIATIONS AND VARIANCE  12  3    C1                                                  
C                                                                               C3                                                  
CAUTHOR  - J.D.HENDERSON                                                        C3                                                  
C                                                                               C3                                                  
CPURPOSE - COMPUTES THE 2-D TIME AVERAGE FOR EACH LEVEL OF A SERIES OF          C3                                                  
C          (MULTI-LEVEL) SETS OF DATA HAVING THE SAME SIZE AND TYPE,            C3                                                  
C          AND OPTIONALLY, VARIANCE AND DEVIATIONS FROM THE TIME MEAN.          C3                                                  
C          NOTE - MAXIMUM NUMBER OF LEVELS MAXL(=$L$).                          C3                                                  
C                                                                               C3                                                  
CINPUT FILE...                                                                  C3                                                  
C                                                                               C3                                                  
C        X1 = SERIES OF MULTI-LEVEL SETS, (DATA MAY BE REAL OR COMPLEX)         C3                                                  
C             ALL RECORDS MUST HAVE THE SAME KIND, NAME AND DIMENSIONS.         C3                                                  
C   X2...X12= (OPTIONAL) ADDITIONAL TIME SERIES OF MULTI-LEVEL SETS.            C3                                                  
C             TIME MEAN OF ALL INPUT FILES IS COMPUTED.                         C3                                                  
C             THEY MUST BE SPECIFIED AFTER OUTPUT VAR FILE ON THE COMMAND LINE. C3                                                  
C                                                                               C3                                                  
COUTPUT FILE...                                                                 C3                                                  
C                                                                               C3                                                  
C       AVG = MULTI-LEVEL SET WHICH IS THE MEAN OF ALL THE SETS IN FILE X.      C3                                                  
C       DEV = (OPTIONAL) DEVIATIONS OF X FROM THE TIME MEAN AVG.                C3                                                  
C       VAR = (OPTIONAL) THE VARIANCE OF ALL THE SETS IN X.                     C3                                                  
C                                                                               C3                                                  
CEXAMPLES:                                                                      C3                                                  
C     timavg x tx         # CALCULATE MULTI-LEVEL TIME AVERAGE.                 C3                                                  
C     timavg x tx xp      # CALCULATE MULTI-LEVEL TIME AVERAGE AND DEVIATIONS.  C3                                                  
C     timavg x tx xp txp2 # CALCULATE MULTI-LEVEL TIME AVERAGE, DEVIATIONS AND  C3                                                  
C                         # VARIANCE .                                          C3                                                  
C     timavg x tx _ txp2  # CALCULATE MULTI-LEVEL TIME AVERAGE AND VARIANCE.    C3                                                  
C                                                                               C3                                                  
C     timavg x1 tx _ _ x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12                      C3                                                  
C     (CALCULATE MULTI-LEVEL TIME AVERAGE OF ALL INPUT FILES x1,...,x12         C3                                                  
C-----------------------------------------------------------------------------                                                      
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMMON/BLANCK/AVG(6534100),VAR(6534100),                                                                                      
     +              X(65341),GAVG(100)                                                                                              
                                                                                                                                    
      LOGICAL OK,SPEC,TIMDEV,TIMVAR                                                                                                 
      INTEGER LEV(100),JBUF(8)                                                                                                      
      COMMON/ICOM/IBUF(8),IDAT(65341)                                                                                               
C                                                                                                                                   
C     * IUA/IUST  - ARRAY KEEPING TRACK OF UNITS ASSIGNED/STATUS.                                                                   
C                                                                                                                                   
      INTEGER IUA(100),IUST(100)                                                                                                    
      COMMON /JCLPNTU/ IUA,IUST                                                                                                     
                                                                                                                                    
      DATA MAXX,MAXG,MAXL / 65341,6534100,100/                                                                                      
      DATA MAXRSZ /65341/                                                                                                           
C---------------------------------------------------------------------                                                              
      NF=16                                                                                                                         
      CALL JCLPNT(NF,11,2,3,4,12,13,14,15,16,17,18,19,20,21,22,6)                                                                   
      IF (NF.LT.3) CALL                            XIT('TIMAVG',-1)                                                                 
      REWIND 11                                                                                                                     
      REWIND 2                                                                                                                      
      TIMDEV=.FALSE.                                                                                                                
      IF (IUST(3).EQ.1) THEN                                                                                                        
        WRITE(6,6005)                                                                                                               
        TIMDEV=.TRUE.                                                                                                               
        REWIND 3                                                                                                                    
      ENDIF                                                                                                                         
      TIMVAR=.FALSE.                                                                                                                
      IF (IUST(4).EQ.1) THEN                                                                                                        
        WRITE(6,6006)                                                                                                               
        TIMVAR=.TRUE.                                                                                                               
        REWIND 4                                                                                                                    
      ENDIF                                                                                                                         
      NINP=1                                                                                                                        
      DO I=2,12                                                                                                                     
        IF (IUST(I+10).EQ.1) THEN                                                                                                   
          NINP=NINP+1                                                                                                               
        ENDIF                                                                                                                       
      ENDDO                                                                                                                         
      WRITE(6,*)'NUMBER OF INPUT FILES=',NINP                                                                                       
                                                                                                                                    
C     * FIND THE NUMBER OF LEVELS, TYPE OF DATA, ETC...                                                                             
                                                                                                                                    
      CALL FILEV(LEV,NLEV,IBUF,11)                                                                                                  
      IF (NLEV.EQ.0) THEN                                                                                                           
        WRITE(6,6010)                                                                                                               
        CALL                                       XIT('TIMAVG',-2)                                                                 
      ENDIF                                                                                                                         
      IF (NLEV.GT.MAXL) CALL                       XIT('TIMAVG',-3)                                                                 
      NWDS=IBUF(5)*IBUF(6)                                                                                                          
      KIND=IBUF(1)                                                                                                                  
      IF (KIND.EQ.NC4TO8("TIME")) CALL             XIT('TIMAVG',-4)                                                                 
      NAME=IBUF(3)                                                                                                                  
      NPACK=MIN(2,IBUF(8))                                                                                                          
      SPEC=(KIND.EQ.NC4TO8("FOUR").OR.KIND.EQ.NC4TO8("SPEC"))                                                                       
      IF (SPEC) NWDS=NWDS*2                                                                                                         
      IF (NWDS.GT.MAXRSZ.OR.NWDS*NLEV.GT.MAXG)                                                                                      
     +     CALL                                    XIT('TIMAVG',-5)                                                                 
      WRITE(6,6020) NAME,NLEV,(LEV(L),L=1,NLEV)                                                                                     
      DO I=1,8                                                                                                                      
        JBUF(I)=IBUF(I)                                                                                                             
      ENDDO                                                                                                                         
      CALL PRTLAB (IBUF)                                                                                                            
C                                                                                                                                   
C     * INITIALIZE AVG AND VAR                                                                                                      
C                                                                                                                                   
      DO I=1,NLEV*NWDS                                                                                                              
        AVG(I)=0.E0                                                                                                                 
        VAR(I)=0.E0                                                                                                                 
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * PROCESS ALL INPUT FILES.                                                                                                    
C                                                                                                                                   
      NSETS=0                                                                                                                       
      DO NF=1,NINP                                                                                                                  
        NSETSNF=0                                                                                                                   
C                                                                                                                                   
C       * TIMESTEP AND LEVEL LOOPS.                                                                                                 
C                                                                                                                                   
 200    CONTINUE                                                                                                                    
        DO L=1,NLEV                                                                                                                 
          IW=(L-1)*NWDS                                                                                                             
C                                                                                                                                   
C         * GET THE NEXT FIELD FROM FILE 11 AND CHECK THE LABEL.                                                                    
C                                                                                                                                   
          CALL GETFLD2(NF+10,X,-1,-1,-1,LEV(L),IBUF,MAXX,OK)                                                                        
          IF (.NOT.OK) THEN                                                                                                         
            IF (L.EQ.1) GO TO 300                                                                                                   
C                                                                                                                                   
C           * (MULTI-LEVEL) SET IN FILE NF IS NOT COMPLETE. ABORT.                                                                  
C                                                                                                                                   
            CALL PRTLAB (IBUF)                                                                                                      
            WRITE(6,6030) NAME,LEV(L)                                                                                               
            CALL                                   XIT('TIMAVG',-6)                                                                 
          ENDIF                                                                                                                     
C                                                                                                                                   
C         * MAKE SURE THAT ALL RECORDS HAVE THE SAME NAME, KIND AND DIMENSIONS                                                      
C                                                                                                                                   
          CALL CMPLBL(0,IBUF,0,JBUF,OK)                                                                                             
          IF (.NOT.OK .OR. IBUF(3).NE.JBUF(3)) THEN                                                                                 
            CALL PRTLAB (IBUF)                                                                                                      
            CALL PRTLAB (JBUF)                                                                                                      
            CALL                                   XIT('TIMAVG',-7)                                                                 
          ENDIF                                                                                                                     
          NPACK=MIN(NPACK,IBUF(8))                                                                                                  
C                                                                                                                                   
C         * CALCULATE THE TIME AVERAGE BY ACCUMULATING THE CURRENT FIELD.                                                           
C                                                                                                                                   
          DO I=1,NWDS                                                                                                               
            AVG(IW+I)=AVG(IW+I)+X(I)                                                                                                
          ENDDO                                                                                                                     
C                                                                                                                                   
C         * ACCUMULATE SQUARED FIELDS                                                                                               
C                                                                                                                                   
          IF (TIMVAR) THEN                                                                                                          
            IF(NSETS.EQ.0)THEN                                                                                                      
C                                                                                                                                   
C             * CALCULATE THE AVERAGED VALUE AT EACH LEVEL FOR THE FIRST SET.                                                       
C             * WE SUBSTRACT THIS VALUE FROM ALL OTHER FIELDS TO INCREASE                                                           
C             * NUMERICAL ACCURACY OF VARIANCE COMPUTATIONS.                                                                        
C                                                                                                                                   
              GAVG(L)=0.E0                                                                                                          
              DO I=1,NWDS                                                                                                           
                GAVG(L)=GAVG(L)+X(I)                                                                                                
              ENDDO                                                                                                                 
              GAVG(L)=GAVG(L)/FLOAT(NWDS)                                                                                           
            ENDIF                                                                                                                   
            DO I=1,NWDS                                                                                                             
              VAR(IW+I)=VAR(IW+I)+(X(I)-GAVG(L))**2                                                                                 
            ENDDO                                                                                                                   
          ENDIF                                                                                                                     
        ENDDO                   ! L=1,NLEV                                                                                          
        NSETS=NSETS+1                                                                                                               
        NSETSNF=NSETSNF+1                                                                                                           
        GO TO 200                                                                                                                   
 300    CONTINUE                                                                                                                    
        WRITE(6,6040) NSETSNF,NAME,NF                                                                                               
        IF (NSETSNF.EQ.0) THEN                                                                                                      
          WRITE(6,6050)                                                                                                             
          CALL                                     XIT('TIMAVG',-8)                                                                 
        ENDIF                                                                                                                       
      ENDDO                     ! NF=11,NINP+10                                                                                     
C                                                                                                                                   
C     * CALCULATE THE AVERAGE FOR EACH LEVEL.                                                                                       
C                                                                                                                                   
      FNI=1.E0/FLOAT(NSETS)                                                                                                         
      DO L=1,NLEV                                                                                                                   
        IW=(L-1)*NWDS                                                                                                               
        DO I=1,NWDS                                                                                                                 
          AVG(IW+I)=AVG(IW+I)*FNI                                                                                                   
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * PUT VARIANCE TO FILE 4 (PACKED AT HIGHEST 2:1).                                                                           
C       * (WE SAVE 'VAR' BEFORE 'AVG' TO AVOID ACCURACY LOSS WHEN NPACK>1)                                                          
C                                                                                                                                   
        IBUF(2)=NSETS                                                                                                               
        IBUF(4)=LEV(L)                                                                                                              
        IBUF(8)=NPACK                                                                                                               
        IF (TIMVAR) THEN                                                                                                            
          DO I=1,NWDS                                                                                                               
            VAR(IW+I)=MAX(0.E0,VAR(IW+I)*FNI-                                                                                       
     &           (AVG(IW+I)-GAVG(L))**2)                                                                                            
          ENDDO                                                                                                                     
          CALL PUTFLD2(4,VAR(IW+1),IBUF,MAXX)                                                                                       
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * PUT THE TIME AVERAGE TO FILE 2 (PACKED AT HIGHEST 2:1).                                                                   
C                                                                                                                                   
        CALL PUTFLD2(2,AVG(IW+1),IBUF,MAXX)                                                                                         
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * CALCULATE TIME DEVIATIONS, IF NECESSARY                                                                                     
C                                                                                                                                   
      IF (.NOT.TIMDEV) GO TO 900                                                                                                    
C                                                                                                                                   
C     * INPUT FILES, TIMESTEP AND LEVEL LOOPS.                                                                                      
C                                                                                                                                   
      DO NF=1,NINP                                                                                                                  
        REWIND NF+10                                                                                                                
 400    CONTINUE                                                                                                                    
        DO L=1,NLEV                                                                                                                 
          IW=(L-1)*NWDS                                                                                                             
C                                                                                                                                   
C         * GET THE NEXT FIELD FROM FILE NF                                                                                         
C                                                                                                                                   
          CALL GETFLD2(NF+10,X,KIND,-1,NAME,LEV(L),IBUF,MAXX,OK)                                                                    
          IF (.NOT.OK) GO TO 450                                                                                                    
C                                                                                                                                   
C         * CALCULATE THE TIME DEVIATIONS OF THE CURRENT FIELD.                                                                     
C         * WRITE THIS DEVIATION TO FILE 3.                                                                                         
C                                                                                                                                   
          DO I=1,NWDS                                                                                                               
            X(I)=X(I)-AVG(IW+I)                                                                                                     
          ENDDO                                                                                                                     
          CALL PUTFLD2(3,X,IBUF,MAXX)                                                                                               
        ENDDO                   ! L=1,NLEV                                                                                          
        GOTO 400                                                                                                                    
 450    CONTINUE                                                                                                                    
      ENDDO                     ! NF=1,NINP                                                                                         
C                                                                                                                                   
C     * NORMAL EXIT.                                                                                                                
C                                                                                                                                   
 900  CONTINUE                                                                                                                      
      WRITE(6,6060) NSETS,NAME                                                                                                      
      CALL PRTLAB (IBUF)                                                                                                            
      CALL                                         XIT('TIMAVG',0)                                                                  
                                                                                                                                    
C---------------------------------------------------------------------                                                              
 6005 FORMAT('0 CALCULATE ALSO DEVIATIONS FROM THE MEAN.')                                                                          
 6006 FORMAT('0 CALCULATE ALSO THE VARIANCE.')                                                                                      
 6010 FORMAT('0..TIMAVG INPUT FILE IS EMPTY')                                                                                       
 6020 FORMAT('0NAME =',A4/'0NLEVS =',I5/                                                                                            
     1       '0LEVELS = ',15I6/100(10X,15I6/))                                                                                      
 6030 FORMAT('0..TIMAVG INPUT ERROR - NAME,L=',2X,A4,I5)                                                                            
 6040 FORMAT('0TIMAVG PROCESSED ',I10,' SETS OF ',A4,' IN FILE ',I5)                                                                
 6050 FORMAT('0..TIMAVG TIME SERIES NAMES INCORRECT OR EMPTY FILE ',I5)                                                             
 6060 FORMAT('0TIMAVG PROCESSED ',I10,' SETS OF ',A4,' IN TOTAL')                                                                   
      END                                                                                                                           
