      PROGRAM TIMCOV                                                                                                                
C     PROGRAM TIMCOV (DEVX,       DEVY,       COV,      OUTPUT,         )       C2                                                  
C    1         TAPE11=DEVX,TAPE12=DEVY,TAPE13=COV,TAPE6=OUTPUT)                 C2                                                  
C     ---------------------------------------------------------                 C2                                                  
C                                                                               C2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)       C2                                                  
C     NOV 20/00 - S.KHARIN (MAKE DEVY OPTIONAL)                                 C2                                                  
C     APR 11/00 - S.KHARIN (SIMPLIFY CODE. RELAX RESTRICTION ON LEVEL NUMBER)                                                       
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     FEB 19/85 - B.DUGAS. (RESTRUCTURE TO USE SCRATCH DISK FILE)                                                                   
C     MAY 07/80 - J.D.HENDERSON                                                                                                     
C                                                                               C2                                                  
CTIMCOV  - COMPUTES MULTI-LEVEL COVARIANCE OF TWO DEVIATION FILES       2  1    C1                                                  
C                                                                               C3                                                  
CAUTHOR  - J.D.HENDERSON                                                        C3                                                  
C                                                                               C3                                                  
CPURPOSE - COMPUTES THE 2-D TIME COVARIANCES FOR EACH LEVEL OF TWO              C3                                                  
C          (MULTI-LEVEL) TIME SERIES DATA SETS OF DEVIATIONS, (ASSUMED          C3                                                  
C          TO BE OF THE SAME SIZE AND TYPE), FROM THEIR RESPECTIVE              C3                                                  
C          TIME MEANS.                                                          C3                                                  
C          NOTE - MAXIMUM NUMBER OF LEVELS MAXL(=$L$).                          C3                                                  
C                                                                               C3                                                  
CINPUT FILES...                                                                 C3                                                  
C                                                                               C3                                                  
C      DEVX  = FIRST  MULTI-LEVEL SERIES OF DEVIATIONS                          C3                                                  
C      DEVY  = (OPTIONAL) SECOND MULTI-LEVEL SERIES OF DEVIATIONS.              C3                                                  
C              IF NOT PRESENT, OR '_' IS SPECIFIED, ASSUME DEVX=DEVY.           C3                                                  
C                                                                               C3                                                  
COUTPUT FILE...                                                                 C3                                                  
C                                                                               C3                                                  
C        COV = MULTI-LEVEL SET OF THE (CO)VARIANCE.                             C3                                                  
C-----------------------------------------------------------------------------                                                      
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      COMMON/BLANCK/COV(1852800),X(18528),Y(18528)                                                                                  
                                                                                                                                    
      LOGICAL OK,SPEC,CV                                                                                                            
      INTEGER LEV(100)                                                                                                              
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
      COMMON/JCOM/JBUF(8),JDAT(18528)                                                                                               
                                                                                                                                    
C     * IUA/IUST  - ARRAY KEEPING TRACK OF UNITS ASSIGNED/STATUS.                                                                   
      INTEGER IUA(100),IUST(100)                                                                                                    
      COMMON /JCLPNTU/ IUA,IUST                                                                                                     
                                                                                                                                    
      DATA MAXX,MAXG,MAXL / 18528,1852800,100 /                                                                                     
      DATA MAXRSZ /18528/                                                                                                           
C---------------------------------------------------------------------                                                              
      NF=4                                                                                                                          
      CALL JCLPNT(NF,11,12,13,6)                                                                                                    
      IF (NF.LT.3) CALL                            XIT('TIMCOV',-1)                                                                 
                                                                                                                                    
      REWIND 11                                                                                                                     
      NOUT=13                                                                                                                       
      IF (NF.EQ.3) THEN                                                                                                             
         CV=.FALSE.                                                                                                                 
         IF (IUST(12).EQ.1) NOUT=12                                                                                                 
      ELSE                                                                                                                          
         CV=.TRUE.                                                                                                                  
         REWIND 12                                                                                                                  
      ENDIF                                                                                                                         
      REWIND NOUT                                                                                                                   
                                                                                                                                    
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('TIMCOV',-2)                                                                 
      ENDIF                                                                                                                         
      IF (NLEV.GT.MAXL) CALL                       XIT('TIMCOV',-3)                                                                 
      WRITE(6,6020) NLEV,(LEV(L),L=1,NLEV)                                                                                          
      NWDS=IBUF(5)*IBUF(6)                                                                                                          
      KIND=IBUF(1)                                                                                                                  
      NPACK=MIN(2,IBUF(8))                                                                                                          
      SPEC=(KIND.EQ.NC4TO8("FOUR").OR.KIND.EQ.NC4TO8("SPEC"))                                                                       
      IF (SPEC) NWDS=NWDS*2                                                                                                         
      IF (NWDS.GT.MAXRSZ) CALL                     XIT('TIMCOV',-4)                                                                 
C                                                                                                                                   
C     * FIND THE MAXIMAL NUMBER OF LEVELS THAT FIT IN AVG                                                                           
C                                                                                                                                   
      NLEVM=MIN(MAXG/NWDS,NLEV)                                                                                                     
      NREAD=(NLEV-1)/NLEVM+1                                                                                                        
      WRITE(6,6030) NLEVM                                                                                                           
      CALL PRTLAB (IBUF)                                                                                                            
C                                                                                                                                   
C     * PROCESS INPUT FILES SEVERAL TIMES, IF NECESSARY.                                                                            
C                                                                                                                                   
      DO N=1,NREAD                                                                                                                  
         REWIND 11                                                                                                                  
         REWIND 12                                                                                                                  
         L1=(N-1)*NLEVM+1                                                                                                           
         L2=MIN(N*NLEVM,NLEV)                                                                                                       
C                                                                                                                                   
C     * READ NLEVM-LEVEL CHUNK FOR THE FIRST TIME STEP.                                                                             
C                                                                                                                                   
         DO L=L1,L2                                                                                                                 
            IW=(L-L1)*NWDS                                                                                                          
            CALL GETFLD2(11,X,KIND,-1,-1,LEV(L),IBUF,MAXX,OK)                                                                       
            IF (.NOT.OK) THEN                                                                                                       
C                                                                                                                                   
C     * SET IN FILE 11 IS NOT COMPLETE. ABORT.                                                                                      
C                                                                                                                                   
               CALL PRTLAB (IBUF)                                                                                                   
               WRITE(6,6050) IBUF(3),LEV(L)                                                                                         
               CALL                                XIT('TIMCOV',-5)                                                                 
            ENDIF                                                                                                                   
            IF (CV) THEN                                                                                                            
               CALL GETFLD2(12,Y,KIND,-1,-1,LEV(L),JBUF,MAXX,OK)                                                                    
               IF (.NOT.OK) THEN                                                                                                    
C                                                                                                                                   
C     * SET IN FILE 12 IS NOT COMPLETE. ABORT.                                                                                      
C                                                                                                                                   
                  CALL PRTLAB (JBUF)                                                                                                
                  WRITE(6,6050) JBUF(3),LEV(L)                                                                                      
                  CALL                             XIT('TIMCOV',-6)                                                                 
               ENDIF                                                                                                                
               DO I=1,NWDS                                                                                                          
                  COV(IW+I)=X(I)*Y(I)                                                                                               
               ENDDO                                                                                                                
            ELSE                                                                                                                    
               DO I=1,NWDS                                                                                                          
                  COV(IW+I)=X(I)*X(I)                                                                                               
               ENDDO                                                                                                                
            ENDIF                                                                                                                   
         ENDDO                                                                                                                      
         NAME1=IBUF(3)                                                                                                              
         NAME2=JBUF(3)                                                                                                              
C                                                                                                                                   
C     * TIMESTEP AND LEVEL LOOPS.                                                                                                   
C                                                                                                                                   
         NSETS=1                                                                                                                    
 200     CONTINUE                                                                                                                   
         DO L=L1,L2                                                                                                                 
C                                                                                                                                   
C     * GET THE NEXT FIELD FROM FILE 11 AND CHECK THE LABEL.                                                                        
C                                                                                                                                   
            CALL GETFLD2(11,X,KIND,-1,NAME1,LEV(L),IBUF,MAXX,OK)                                                                    
            IF (.NOT.OK) THEN                                                                                                       
               IF (L.EQ.L1) GO TO 300                                                                                               
C                                                                                                                                   
C     * SET IN FILE 11 IS NOT COMPLETE. ABORT.                                                                                      
C                                                                                                                                   
               CALL PRTLAB (IBUF)                                                                                                   
               WRITE(6,6050) NAME1,LEV(L)                                                                                           
               CALL                                XIT('TIMCOV',-7)                                                                 
            ENDIF                                                                                                                   
C                                                                                                                                   
C     * GET THE NEXT FIELD FROM FILE 12 AND CHECK THE LABEL.                                                                        
C                                                                                                                                   
            IF (CV) THEN                                                                                                            
               CALL GETFLD2(12,Y,KIND,-1,NAME2,LEV(L),JBUF,MAXX,OK)                                                                 
               IF (.NOT.OK) THEN                                                                                                    
                  IF (L.EQ.L1) GO TO 300                                                                                            
C                                                                                                                                   
C     * SET IN FILE 12 IS NOT COMPLETE. ABORT.                                                                                      
C                                                                                                                                   
                  CALL PRTLAB (JBUF)                                                                                                
                  WRITE(6,6050) NAME2,LEV(L)                                                                                        
                  CALL                             XIT('TIMCOV',-8)                                                                 
               ENDIF                                                                                                                
               CALL CMPLBL(0,IBUF,0,JBUF,OK)                                                                                        
               IF (.NOT.OK) THEN                                                                                                    
                  CALL PRTLAB (IBUF)                                                                                                
                  CALL PRTLAB (JBUF)                                                                                                
                  CALL                             XIT('TIMCOV',-9)                                                                 
               ENDIF                                                                                                                
C                                                                                                                                   
C        * CALCULATE THE COVARIANCE BY ACCUMULATING THE PRODUCT OF CURRENT FIELDS.                                                  
C                                                                                                                                   
               IW=(L-L1)*NWDS                                                                                                       
               DO I=1,NWDS                                                                                                          
                  COV(IW+I)=COV(IW+I)+X(I)*Y(I)                                                                                     
               ENDDO                                                                                                                
            ELSE                                                                                                                    
C                                                                                                                                   
C        * CALCULATE THE VARIANCE BY ACCUMULATING THE SQUARE OF CURRENT FIELD.                                                      
C                                                                                                                                   
               IW=(L-L1)*NWDS                                                                                                       
               DO I=1,NWDS                                                                                                          
                  COV(IW+I)=COV(IW+I)+X(I)*X(I)                                                                                     
               ENDDO                                                                                                                
            ENDIF                                                                                                                   
         ENDDO                                                                                                                      
         NSETS=NSETS+1                                                                                                              
         GO TO 200                                                                                                                  
C                                                                                                                                   
C     * CALCULATE THE AVERAGE FOR EACH LEVEL.                                                                                       
C                                                                                                                                   
 300     CONTINUE                                                                                                                   
         FNI=1.E0/FLOAT(NSETS)                                                                                                      
         DO L=L1,L2                                                                                                                 
            IW=(L-L1)*NWDS                                                                                                          
            DO I=1,NWDS                                                                                                             
               X(I)=COV(IW+I)*FNI                                                                                                   
            ENDDO                                                                                                                   
C                                                                                                                                   
C     * PUT THE TIME AVERAGE TO FILE 13 (PACKED AT HIGHEST 2:1).                                                                    
C                                                                                                                                   
            IBUF(2)=NSETS                                                                                                           
            IBUF(3)=NC4TO8(" COV")                                                                                                  
            IBUF(4)=LEV(L)                                                                                                          
            IBUF(8)=NPACK                                                                                                           
            CALL PUTFLD2(NOUT,X,IBUF,MAXX)                                                                                          
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * NORMAL EXIT.                                                                                                                
C                                                                                                                                   
      WRITE(6,6060) NSETS,NAME1,NAME2                                                                                               
      CALL PRTLAB (IBUF)                                                                                                            
      CALL                                         XIT('TIMCOV',0)                                                                  
                                                                                                                                    
C---------------------------------------------------------------------                                                              
 6010 FORMAT('0..TIMCOV INPUT FILE IS EMPTY')                                                                                       
 6020 FORMAT('0NLEVS =',I5/'0LEVELS = ',15I6/100(10X,15I6/))                                                                        
 6030 FORMAT(' ',I5,' LEVELS CAN BE KEPT IN MEMORY.')                                                                               
 6050 FORMAT('0..TIMCOV INPUT ERROR - NAME,L=',2X,A4,I5)                                                                            
 6060 FORMAT('0TIMCOV PROCESSED ',I10,' SETS OF ',A4,1X,A4)                                                                         
      END                                                                                                                           
