CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC               
C     TS*PF.FILANC                                      C               
C                                                        C              
C     THIS IS A LANCZOS FILTER PROGRAM                   C              
C                                                        C              
C       UNIT   3 = INPUT DATA SET                        C              
C       UNIT  12 = LOW PASS DATA.....OUTPUT              C              
C       UNIT  13 = FILTER WEIGHTS....OUTPUT              C              
c	
c	dt is time step in hrs
c	n  is the number of data
c	val is the time series on input and the filtered
c		series on output
c	riww is the # of filter weights used
c	tco  is a filter design parameter close to the
c		half power point
CAKM                                                      C             
C11-11-85CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC              
      SUBROUTINE FILTER (DT,N,VAL,RIWW,TCO)                             
C                                                                       
      REAL VAL(*),LLPVAL(99000),H(99000),DT,TCO
      INTEGER N,IWW                                                     
      CHARACTER*80 H1                                                   
      CHARACTER*32 H2                                                   
C                                                                       
      D(Z)=SIN(Z)/Z                                                     
C                                                                       
      PI = 4.*ATAN(1.)                                                  
      PI2 = 2.*PI                                                       
C                                                                       
C                                                                       
                                                                        
      IWW = RIWW/DT                                                     
      if (int(iww/2)*2.eq.iww) iww = iww+1
      IWW2=(IWW-1)/2                                                    
      T=FLOAT(IWW-1)                                                    
C                                                                       
        KKK=0                                                           
C                                                                       
      DO 25 I=1,IWW2                                                    
         LLPVAL(I)=0.0                                                  
25    CONTINUE                                                          
      DO 26 I=N-IWW2+1,N                                                
          LLPVAL(I)=0.0                                                 
26    CONTINUE                                                          
C                                                                       
      OMEGA=PI2/TCO*DT                                                  
      H0=OMEGA/PI                                                       
      CON=PI2/FLOAT(IWW)                                                
C                                                                       
C     COMPUTE WINDOW WEIGHTS                                            
C                                                                       
      SUM=H0                                                            
      DO 30 I=1,IWW2                                                    
         H(I)=H0*D(FLOAT(I)*OMEGA)*D(FLOAT(I)*CON)                      
         SUM=SUM+2.0*H(I)                                               
   30 CONTINUE                                                          
C                                                                       
C     NORMALIZE EACH WEIGHT BY THE SUM OF ALL WEIGHTS                   
C                                                                       
      H0=H0/SUM                                                         
      DO 35 I=1,IWW2                                                    
         H(I)=H(I)/SUM                                                  
199	format(i4,f10.4,3f10.1)
   35 CONTINUE                                                          
C                                                                       
C                                                                       
      DO 55 I=IWW2+1,N-IWW2                                             
         K=I                                                            
         SUM=0.0                                                        
         TEMP=H0*VAL(I)                                                 
         DO 50 J=1,IWW2                                                 
            TEMP=TEMP+H(J)*(VAL(I+J)+VAL(I-J))                          
   50    CONTINUE                                                       
         LLPVAL(I)=TEMP                                                 
         TEMP=0                                                         
   55 CONTINUE                                                          
C                                                                       
      DO 56 I=1,N                                                       
         VAL(I) = LLPVAL(I)                                             
56    CONTINUE                                                          
      RETURN                                                            
      END                                                               
C                                                                       
C  LANCZOS BAND PASS FILTER                                             
C  09-24-1991 AMUE                                                      
C                                                                       
      SUBROUTINE BFILTER(DT,NDT,U,RIWW,TCO,RIWW1,TCO1)                  
      REAL DT,U(*),RIWW,TCO,riww1,tco1                                             
      REAL U1(99000),U2(99000)                                            
C                                                                       
	write(6,*)'low-pass: ',tco1,riww1
	write(6,*)'high-pass: ',tco,riww
      DO 1 I=1,NDT                                                      
         U1(I) = U(I)                                                   
         U2(I) = U(I)                                                   
1     CONTINUE                                                          
C                                                                       
C     TCO1 (DATA BELOW RETAINED, HIGH) > TCO (DATA ABOVE RETAINED, LOW) 
C                                                                       
      CALL FILTER(DT,NDT,U1,RIWW1,TCO1)                                 
      CALL FILTER(DT,NDT,U2,RIWW,TCO)                                   
      IWW2 = INT(RIWW1/DT-1)/2                                          
C                                                                       
      DO 2 I=IWW2+1,NDT-IWW2                                            
         U(I) = U2(I)-U1(I)                                             
2     CONTINUE                                                          
      DO 3 I=1,IWW2                                                     
         U(I) = 0.                                                      
3     CONTINUE                                                          
      DO 4 I=NDT-IWW2+1,NDT                                             
         U(I) = 0.                                                      
4     CONTINUE                                                          
      RETURN                                                            
      END                                                               
C                                                                       
C  LANCZOS HIGH-PASS FILTER                                             
C  09-25-1991 AMUE                                                      
C                                                                       
      SUBROUTINE HFILTER(DT,NDT,U,RIWW,TCO)                             
      REAL DT,U(*),RIWW,TCO                                             
      REAL U1(99000)                                                     
C                                                                       
      DO 1 I=1,NDT                                                      
         U1(I) = U(I)                                                   
1     CONTINUE                                                          
C                                                                       
      CALL FILTER(DT,NDT,U1,RIWW,TCO)                                   
      IWW2 = INT(RIWW/DT-1)/2                                           
C                                                                       
      DO 2 I=IWW2+1,NDT-IWW2                                            
         U(I) = U(I)-U1(I)                                              
2     CONTINUE                                                          
      DO 3 I=1,IWW2                                                     
         U(I) = 0.                                                      
3     CONTINUE                                                          
      DO 4 I=NDT-IWW2+1,NDT                                             
         U(I) = 0.                                                      
4     CONTINUE                                                          
      RETURN                                                            
      END