      PROGRAM CROSPEC                                                                                                               
C     PROGRAM CROSPEC (X,       Y,       SPEC,       INPUT,      OUTPUT,)       H2                                                  
C    1           TAPE1=X, TAPE2=Y, TAPE3=SPEC, TAPE5=INPUT,TAPE6=OUTPUT)                                                            
C     ---------------------------------------------------------------           H2                                                  
C                                                                               H2                                                  
C     NOV 02/2004 - S.KHARIN (MAKE Y OPTIONAL. ADD AR1 AND MODIFIED DANIELL     H2                                                  
C                             SPECTRAL ESTIMATORS. APPROXIMATE THE PERIODOGRAM  H2                                                  
C                             AT THE ZERO FREQUENCY BY THE VALUE AT THE CLOSEST H2                                                  
C                             FREQUENCY. ADJUST SPECTRA AND D.F. FOR TAPERING.) H2                                                  
C     JUL 05/2001 - S.KHARIN (ALLOW "EPS" TOLERANCE IN THE COHERENCE                                                                
C                             CALCULATIONS)                                                                                         
C     MAR 06/2000 - S.KHARIN                                                                                                        
C                                                                               H2                                                  
CCROSPEC - CALCUALTES VARIANCE, PHASE AND COHERENCE SPECTRA.            2  1 C  H1                                                  
C                                                                               H3                                                  
CAUTHOR  - S.KHARIN                                                             H3                                                  
C                                                                               H3                                                  
CPURPOSE - CALCULATES VARIANCE, PHASE AND COHERENCE SPECTRA OF TIME SERIES.     H3                                                  
C                                                                               H3                                                  
CINPUT FILES...                                                                 H3                                                  
C                                                                               H3                                                  
C     X = 1ST SET OF TIME SERIES IN CCRN FORMAT (KIND=4HTIME).                  H3                                                  
C     Y = (OPTIONAL) 2ND SET OF TIME SERIES IN CCRN FORMAT (KIND=4HTIME).       H3                                                  
C         THIS FILE CAN BE OMITTED, OR THE UNDERSCORE '_' CAN BE USED INSTEAD.  H3                                                  
C                                                                               H3                                                  
COUTPUT FILE...                                                                 H3                                                  
C                                                                               H3                                                  
C  SPEC = OUTPUT FILE WITH SPECTRA ESTIMATES FOR EACH TIME SERIES IN INPUT      H3                                                  
C         FILES X (AND Y). PHASE AND COHERENCE SPECTRA ARE CALCULATED ONLY      H3                                                  
C         IF THE SECOND FILE Y IS SPECIFIED.                                    H3                                                  
C                                                                               H3                                                  
C         VARIANCE SPECTRA ARE NORMALIZED SUCH THAT THEIR INTEGRAL OVER         H3                                                  
C         *ALL* (POSITIVE AND NEGATIVE) FREQUENCIES (-0.5...0.5) IS EQUAL       H3                                                  
C         (APPROXIMATELY) TO THE VARIANCE OF THE CORRESPONDING TIME SERIES.     H3                                                  
C         FOR EXAMPLE, THE POWER SPECTRUM OF A WHITE NOISE OF THE UNIT          H3                                                  
C         VARIANCE IS APPROXIMATELY EQUAL TO ONE AT ALL FREQUENCIES.            H3                                                  
C                                                                               H3                                                  
C         VAR(X)=2*SUM(SPEX_I*BW), I=1,...N, WHERE N IS NUMBER OF RESOLVED      H3                                                  
C         POSITIVE FREQUENCIES AND BW IS THE BANDWIDTH BETWEEN THE FREQUENCIES. H3                                                  
C         NORMALLY, N=NTIME/2 AND BW=1/NTIME, WHERE NTIME IS THE LENGTH         H3                                                  
C         OF TIME SERIES.                                                       H3                                                  
C                                                                               H3                                                  
C         THE RECORD ORDER AND FIELD NAMES IN FILE SPEC ARE THE FOLLOWING:      H3                                                  
C           OMEG = FREQUENCIES.                                                 H3                                                  
C           SPEX = VARIANCE SPECTRUM OF X.                                      H3                                                  
C           SPX1 = LOWER CONFIDENCE BOUND OF SPEX.                              H3                                                  
C           SPX2 = UPPER CONFIDENCE BOUND OF SPEX.                              H3                                                  
C           SPRX = VARIANCE SPECTRUM OF AR1 PROCESS FITTED TO X.                H3                                                  
C           SRX1 = UPPER CONFIDENCE BOUND OF SPRX.                              H3                                                  
C           SRX2 = LOWER CONFIDENCE BOUND OF SPRX.                              H3                                                  
C         AND OPTIONALLY, IF Y IS PRESENT:                                      H3                                                  
C           SPEY = VARIANCE SPECTRUM OF Y.                                      H3                                                  
C           SPY1 = LOWER CONFIDENCE BOUND OF SPEY.                              H3                                                  
C           SPY2 = UPPER CONFIDENCE BOUND OF SPEY.                              H3                                                  
C           SPRY = VARIANCE SPECTRUM OF AR1 PROCESS FITTED TO Y.                H3                                                  
C           SRY1 = UPPER CONFIDENCE BOUND OF SPRY.                              H3                                                  
C           SRY2 = LOWER CONFIDENCE BOUND OF SPRY.                              H3                                                  
C           COHE = SQUARED COHERENCE SPECTRUM.                                  H3                                                  
C           COH0 = CRIT. VALUE FOR STATISTICALLY SIGNIFICANT COHERENCE.         H3                                                  
C           COH1 = LOWER CONFIDENCE BOUND OF COHE.                              H3                                                  
C           COH2 = UPPER CONFIDENCE BOUND OF COHE.                              H3                                                  
C           PHAS = PHASE LAG SPECTRUM (IN DEGREES)                              H3                                                  
C           PHA1 = LOWER CONFIDENCE BOUND OF PHAS.                              H3                                                  
C           PHA2 = UPPER CONFIDENCE BOUND OF PHAS.                              H3                                                  
C           ... THE SAME AS ABOVE BUT FOR THE SECOND PAIR OF TIME SERIES ETC.   H3                                                  
C                                                                               H3                                                  
CINPUT PARAMETERS...                                                                                                                
C                                                                                                                                   
C     NTYPE  = TYPE OF THE SPECTRAL ESTIMATOR:                                  H5                                                  
C            = 0 CALCULATE VARIANCE SPECTRUM OF AR1 PROCESS ONLY                H5                                                  
C            = 1 CHUNCK ESTIMATOR,                                              H5                                                  
C            = 2 DANIELL ESTIMATOR (MOVING AVERAGE FILTER),                     H5                                                  
C            =-2 MODIFIED DANIELL ESTIMATOR (HALF WEIGHT TO THE END VALUES),    H5                                                  
C            = 3 BARTLETT ESTIMATOR,                                            H5                                                  
C            = 4 PARZEN ESTIMATOR.                                              H5                                                  
C                                                                               H5                                                  
C     LFLT   - THIS PARAMETER SPECIFIES THE SMOOTHING FACTOR OF THE RAW         H5                                                  
C              PERIODOGRAM. EQUIVALENT DEGREES OF FREEDOM (EDF) AND             H5                                                  
C              EQUIVALENT BANDWIDTH (EBW) OF A SPECTRAL ESTIMATOR ARE           H5                                                  
C                   EDF = 2*LFLT,                                               H5                                                  
C                   EBW = LFLT/NTIME.                                           H5                                                  
C            = 1, THE RESULT IS THE RAW PERIODOGRAM (NO SMOOTHING).             H5                                                  
C                                                                               H5                                                  
C            FOR CHUNK ESTIMATOR (NTYPE=1):                                     H5                                                  
C            = NUMBER OF CHUNKS.                                                H5                                                  
C                                                                               H5                                                  
C            FOR DANIELL ESTIMATOR (NTYPE=2):                                   H5                                                  
C            = LENGTH OF THE MOVING AVERAGE WINDOW TO SMOOTH THE RAW            H5                                                  
C              PERIODOGRAM. LFLT MUST BE ODD IN THIS CASE. EVEN VALUES ARE      H5                                                  
C              INCREASED BY ONE.                                                H5                                                  
C                                                                               H5                                                  
C            FOR BARTLETT OR PARZEN ESTIMATORS (NTYPE=3,4):                     H5                                                  
C              LFLT IS USED TO CALCULATE THE CUTOFF POINT OF LAG WINDOWS.       H5                                                  
C              BARTLETT CUTOFF POINT M=1.5 *NTIME/LFLT.                         H5                                                  
C              PARZEN   CUTOFF POINT M=1.86*NTIME/LFLT.                         H5                                                  
C                                                                               H5                                                  
C           AN EXCERPT FROM THE BOOK "STATISTICAL ANALYSIS IN CLIMATE RESEARCH" H5                                                  
C           BY H. VON STORCH AND F. ZWIERS:                                     H5                                                  
C              "THE CHUNK ESTIMATOR IS GENERALLY SUITABLE FOR PROBLEMS IN       H5                                                  
C              WHICH THERE IS A NATURAL CHUNK LENGTH. ... (OTHERWISE)           H5                                                  
C              WE HAVE A SLIGHT PREFERENCE FOR THE DANIELL AND PARZEN           H5                                                  
C              ESTIMATOR OVER THE BARTLETT ESTIMATOR, FOR WHICH VARIANCE        H5                                                  
C              LEAKAGE THROUGH SIDE LOBES IS MORE AN ISSUE."                    H5                                                  
C                                                                               H5                                                  
C     CONF   = CONFIDENCE INTERVAL (IN %) FOR VARIANCE, COHERENCE AND PHASE     H5                                                  
C              SPECTRA ESTIMATES. TYPICAL VALUES ARE 90. OR 95.                 H5                                                  
C                                                                               H5                                                  
C     LTREND = 0 REMOVE TIME MEAN,                                              H5                                                  
C            = 1 REMOVE LINEAR TREND.                                           H5                                                  
C                                                                               H5                                                  
C     TAPER  = 0...0.5 FRACTION OF TIME SERIES TO BE TAPERED AT EACH END.       H5                                                  
C              A SPLIT COSINE TAPER IS APPLIED TO {TAPER * LENGTH(X)} POINTS    H5                                                  
C              AT EACH END OF THE TIME SERIES. TYPICAL VALUES ARE 0.05...0.1.   H5                                                  
C                                                                               H5                                                  
C   NOTE -     THE NUMBER OF FREQUENCIES AND THE LENGTH OF OUTPUT SPECTRAL      H5                                                  
C              ESTIMATES IS EQUAL TO LCH=NTIME/(2*LFLT) FOR CHUNK ESTIMATOR AND H5                                                  
C              NTIME/2 FOR ALL OTHER ESTIMATORS.                                H5                                                  
C                                                                               H5                                                  
CEXAMPLE OF INPUT CARDS.                                                        H5                                                  
C                                                                               H5                                                  
C*CROSPEC          1        20       95.    1       0.1                         H5                                                  
C(CHUNK ESTIMATOR. TIME SERIES WILL BE DIVIDED IN 20 CHUNKS.                    H5                                                  
C CONFIDENCE INTERVAL IS 95%. REMOVE LINEAR TREND. TAPER 10% FROM EACH END.)    H5                                                  
C                                                                               H5                                                  
C*CROSPEC          2        11       95.    1       0.1                         H5                                                  
C(THE SAME AS ABOVE BUT FOR DANIELL ESTIMATOR WITH 11-POINT MOVING AVG. FILTER) H5                                                  
C--------------------------------------------------------------------------                                                         
C                                                                                                                                   
C      IMPLICIT REAL (A-H,O-Z),                                                                                                     
C     +INTEGER (I-N)                                                                                                                
      IMPLICIT NONE                                                                                                                 
      LOGICAL OK                                                                                                                    
C                                                                                                                                   
C     * NTMAX  = MAX LENGTH OF TIME SERIES.                                                                                         
C     * NFMAX = MAX NUMBER OF FREQUENCIES.                                                                                          
C                                                                                                                                   
      INTEGER NTMAX,NFMAX                                                                                                           
      PARAMETER (NTMAX=18528,NFMAX=NTMAX/2)                                                                                         
C                                                                                                                                   
C     ARRAYS FOR INPUT FIELDS                                                                                                       
C                                                                                                                                   
      REAL X(NTMAX),Y(NTMAX)                                                                                                        
C                                                                                                                                   
C     CROSS    = LOGICAL FLAG FOR THE PRESENCE OF THE SECOND INPUT FILE.                                                            
C                                                                                                                                   
      LOGICAL CROSS                                                                                                                 
C                                                                                                                                   
C     ARRAYS FOR OUTPUT FIELDS                                                                                                      
C                                                                                                                                   
C     FREQ     = FREQUENCIES AT WHICH POWER SPECTRA ARE ESTIMATED                                                                   
C     VARIX    = VARIANCE SPECTRUM ESTIMATE FOR X.                                                                                  
C     VARIX1,2 = LOWER AND UPPER CONFIDENCE BOUNDS OF VARIX.                                                                        
C     VARRX    = VARIANCE SPECTRUM ESTIMATE OF AR1 PROCESS FITTED TO X                                                              
C     VARRX1,2 = LOWER AND UPPER CONFIDENCE BOUNDS OF VARRX.                                                                        
C     VARIY    = VARIANCE SPECTRUM ESTIMATE FOR Y.                                                                                  
C     VARIY1,2 = LOWER AND UPPER CONFIDENCE BOUNDS OF VARIY.                                                                        
C     VARRY    = VARIANCE SPECTRUM ESTIMATE OF AR1 PROCESS FITTED TO Y                                                              
C     VARRY1,2 = LOWER AND UPPER CONFIDENCE BOUNDS OF VARRY.                                                                        
C     COHER    = COHERENCE SPECTRUM ESTIMATE.                                                                                       
C     COHER1,2 = LOWER AND UPPER CONFIDENCE BOUNDS OF COHER                                                                         
C                WHEN THE NULL HYPOTHESIS COHER=0 IS REJECTED.                                                                      
C     COHER0   = CRITICAL VALUE FOR STATISTICALLY SIGNIFICANT COHERENCE.                                                            
C     PHASE    = PHASE SPECTRUM ESTIMATE.                                                                                           
C     PHASE1,2 = LOWER AND UPPER CONFIDENCE BOUNDS OF PHASE.                                                                        
C                                                                                                                                   
      REAL FREQ(NFMAX),                                                                                                             
     1     VARIX(NFMAX),VARIX1(NFMAX),VARIX2(NFMAX),                                                                                
     2     VARIY(NFMAX),VARIY1(NFMAX),VARIY2(NFMAX),                                                                                
     3     VARRX(NFMAX),VARRX1(NFMAX),VARRX2(NFMAX),                                                                                
     4     VARRY(NFMAX),VARRY1(NFMAX),VARRY2(NFMAX),                                                                                
     5     COHER(NFMAX),COHER1(NFMAX),COHER2(NFMAX),COHER0(NFMAX),                                                                  
     6     PHASE(NFMAX),PHASE1(NFMAX),PHASE2(NFMAX)                                                                                 
C                                                                                                                                   
C     INTERNAL ARRAYS                                                                                                               
C                                                                                                                                   
      REAL AX(NFMAX),BX(NFMAX),AY(NFMAX),BY(NFMAX),FILT(-NFMAX:NFMAX),                                                              
     1     PERXX(0:NFMAX),PERYY(0:NFMAX),PERXY(0:NFMAX),PERYX(0:NFMAX)                                                              
C                                                                                                                                   
      CHARACTER*1 SGNF                                                                                                              
C                                                                                                                                   
      INTEGER IBUF,IDAT,MAXX                                                                                                        
      COMMON/ICOM/IBUF(8),IDAT(18528)                                                                                               
C                                                                                                                                   
C     UA/UST  - ARRAYS FOR KEEPING TRACK OF UNITS ASSIGNED/STATUS.                                                                  
C                                                                                                                                   
      INTEGER IUA(100),IUST(100)                                                                                                    
      COMMON /JCLPNTU/ IUA,IUST                                                                                                     
C                                                                                                                                   
      REAL EPS                                                                                                                      
      REAL PI,PI2,VARX,VARX1,VARY,VARY1,VARZX,VARZY,FR,XI,YI,XAVG,YAVG,                                                             
     1     RHOX,RHOX2,RHOY,RHOY2,FTAPER,U2,U4,ALX,BLX,ALY,BLY,SN2,                                                                  
     2     ARG,COSA,SINA,FILTS,EDF,EBW,COH1,COH2,PERXS,PERYS,CVFD,CVTD,                                                             
     3     FD,TD,STND,P,P1,P2,VARL,VARU,CVCHI2D,CVSTND,TAPER,                                                                       
     4     QUAD,COHE,CONF                                                                                                           
      INTEGER NTIME,NTIME1,I,J,N,N0,N1,NFR,LFLT,LCH,NTAPER,NR,I1,I2,                                                                
     1     ID,JJ,MCUT,NTYPE,LTREND,NFF,NFOUT                                                                                        
      INTEGER NC4TO8                                                                                                                
      DATA MAXX/18528/,EPS/1.E-12/                                                                                                  
C--------------------------------------------------------------------------                                                         
      PI=4.E0*ATAN(1.E0)                                                                                                            
      PI2=2.E0*PI                                                                                                                   
C                                                                                                                                   
C     * ASSIGN FILES TO FORTRAN UNITS                                                                                               
C                                                                                                                                   
      NFF=5                                                                                                                         
      CALL JCLPNT(NFF,1,2,3,5,6)                                                                                                    
      NFF=NFF-2                                                                                                                     
      REWIND 1                                                                                                                      
      REWIND 2                                                                                                                      
      CROSS=.FALSE.                                                                                                                 
      IF (NFF.EQ.3) THEN                                                                                                            
        CROSS=.TRUE.                                                                                                                
        REWIND 3                                                                                                                    
        WRITE(6,'(A)')                                                                                                              
     1       ' *** CROSS-SPECTRAL ANALYSIS OF TWO TIME SERIES. ***'                                                                 
        NFOUT=3                                                                                                                     
      ELSE                                                                                                                          
        WRITE(6,'(A)')                                                                                                              
     1       ' *** SPECTRAL ANALYSIS OF ONE TIME SERIES. ***'                                                                       
        IF (IUST(2).EQ.0) THEN                                                                                                      
          NFOUT=3                                                                                                                   
        ELSE                                                                                                                        
          NFOUT=2                                                                                                                   
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * READ INPUT PARAMETERS FROM CARD                                                                                             
C                                                                                                                                   
      READ(5,5010,END=900) NTYPE,LFLT,CONF,LTREND,TAPER                         H4                                                  
C                                                                                                                                   
C     * PRINT INPUT PARAMETERS                                                                                                      
C                                                                                                                                   
      WRITE (6,6000) NTYPE,LFLT,CONF,LTREND,TAPER                                                                                   
C                                                                                                                                   
C     * CHECK INPUT PARAMETERS                                                                                                      
C                                                                                                                                   
      IF (NTYPE.EQ.0) THEN                                                                                                          
        WRITE(6,'(A)')' AR1 ESTIMATOR.'                                                                                             
      ELSE IF (NTYPE.EQ.1) THEN                                                                                                     
        WRITE(6,'(A)')' CHUNK ESTIMATOR.'                                                                                           
      ELSE IF (NTYPE.EQ.2) THEN                                                                                                     
        WRITE(6,'(A)')' DANIELL ESTIMATOR.'                                                                                         
      ELSE IF (NTYPE.EQ.-2) THEN                                                                                                    
        WRITE(6,'(A)')' MODIFIED DANIELL ESTIMATOR.'                                                                                
      ELSE IF (NTYPE.EQ.3) THEN                                                                                                     
        WRITE(6,'(A)')' BARTLETT ESTIMATOR.'                                                                                        
      ELSE IF (NTYPE.EQ.4) THEN                                                                                                     
        WRITE(6,'(A)')' PARZEN ESTIMATOR.'                                                                                          
      ELSE                                                                                                                          
        WRITE(6,'(A)')                                                                                                              
     1       ' *** ERROR: INVALID NTYPE.'                                                                                           
        CALL                                       XIT('CROSPEC',-1)                                                                
      ENDIF                                                                                                                         
      IF (ABS(NTYPE).EQ.2 .AND. MOD(LFLT,2).EQ.0) THEN                                                                              
        LFLT=LFLT+1                                                                                                                 
        WRITE(6,'(A/A,I5)')                                                                                                         
     1       ' *** WARNING: LFLT MUST BE ODD FOR DANIELL SPEC.',                                                                    
     2       '     LFLT IS INCREASED BY ONE. NEW LFLT=',LFLT                                                                        
      ENDIF                                                                                                                         
      IF (NTYPE.EQ.-2.AND.LFLT.EQ.1) THEN                                                                                           
        WRITE(6,'(A)')                                                                                                              
     1       ' *** ERROR: LFLT MUST BE > 1 FOR',                                                                                    
     2       ' MODIFIED DANIELL ESTIMATOR'                                                                                          
        CALL                                       XIT('CROSPEC',-2)                                                                
      ENDIF                                                                                                                         
      IF (TAPER.LT.0.E0 .OR. TAPER.GT.0.5E0) THEN                                                                                   
        WRITE(6,'(A/A/A)')                                                                                                          
     1       ' *** ERROR: INVALID TAPER VALUE.',                                                                                    
     2       '     TAPER VALUE MUST BE IN 0...0.5.',                                                                                
     3       '     TYPICAL VALUES ARE 0.05, OR 0.1.'                                                                                
        CALL                                       XIT('CROSPEC',-3)                                                                
      ENDIF                                                                                                                         
      IF (CONF.LT.50.E0 .OR. CONF.GE.100.E0) THEN                                                                                   
        WRITE(6,'(A/A)')                                                                                                            
     1       ' *** ERROR: INVALID CONFIDENCE INTERVAL.',                                                                            
     2       '     TYPICAL VALUES ARE 90., 95., OR 99.'                                                                             
        CALL                                       XIT('CROSPEC',-4)                                                                
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * DO FOR ALL RECORDS IN INPUT FILES                                                                                           
C                                                                                                                                   
      NR=0                                                                                                                          
 100  CONTINUE                                                                                                                      
C                                                                                                                                   
C     * READ NEXT RECORD FROM FILE X                                                                                                
C                                                                                                                                   
      CALL GETFLD2(1,X,NC4TO8("TIME"),-1,-1,-1,IBUF,MAXX,OK)                                                                        
      IF(.NOT.OK) THEN                                                                                                              
        IF(NR.EQ.0) THEN                                                                                                            
          CALL                                     XIT('CROSPEC',-5)                                                                
        ENDIF                                                                                                                       
        WRITE(6,'(I5,A)') NR,' RECORDS PROCCESSED.'                                                                                 
        CALL                                       XIT('CROSPEC',0)                                                                 
      ENDIF                                                                                                                         
      NR=NR+1                                                                                                                       
      IF (NR.EQ.1) THEN                                                                                                             
        CALL PRTLAB(IBUF)                                                                                                           
        NTIME=IBUF(5)                                                                                                               
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * ABORT IF TIME SERIES LENGTH HAS CHANGED                                                                                     
C                                                                                                                                   
      IF (NR.GT.1.AND.NTIME.NE.IBUF(5)) CALL       XIT('CROSPEC',-6)                                                                
                                                                                                                                    
C                                                                                                                                   
C     * READ NEXT RECORD FROM FILE Y                                                                                                
C                                                                                                                                   
      IF (CROSS) THEN                                                                                                               
        CALL GETFLD2(2,Y,NC4TO8("TIME"),-1,-1,-1,IBUF,MAXX,OK)                                                                      
        IF(.NOT.OK) THEN                                                                                                            
          WRITE(6,'(2A)')                                                                                                           
     1         ' *** ERROR:',                                                                                                       
     2         ' NUMBER OF RECORDS IN FILE Y IS LESS THAN IN FILE X.'                                                               
          CALL                                     XIT('CROSPEC',-7)                                                                
        ENDIF                                                                                                                       
        IF (NR.EQ.1) CALL PRTLAB(IBUF)                                                                                              
C                                                                                                                                   
C       * MAKE SURE THAT THE TIME SERIES ARE OF THE SAME LENGTH                                                                     
C                                                                                                                                   
        IF (NTIME.NE.IBUF(5)) THEN                                                                                                  
          WRITE(6,'(A)')                                                                                                            
     1         ' *** ERROR: TIME SERIES LENGTHS ARE DIFFERENT.'                                                                     
          CALL                                     XIT('CROSPEC',-8)                                                                
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      IF (NR.EQ.1)                                                                                                                  
     1     WRITE(6,'(A,I5)')' LENGTH OF TIME SERIES NTIME=',NTIME                                                                   
C                                                                                                                                   
C     * SUBTRACT TIME MEAN FROM X                                                                                                   
C                                                                                                                                   
      XAVG=0.E0                                                                                                                     
      DO I=1,NTIME                                                                                                                  
        XAVG=XAVG+X(I)                                                                                                              
      ENDDO                                                                                                                         
      XAVG=XAVG/FLOAT(NTIME)                                                                                                        
      DO I=1,NTIME                                                                                                                  
        X(I)=X(I)-XAVG                                                                                                              
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * SUBTRACT LINEAR TREND FROM X                                                                                                
C                                                                                                                                   
      IF (LTREND.EQ.1) THEN                                                                                                         
        XI=0.E0                                                                                                                     
        DO I=1,NTIME                                                                                                                
          XI=XI+X(I)*FLOAT(I)                                                                                                       
        ENDDO                                                                                                                       
        SN2=FLOAT(NTIME)*(FLOAT(NTIME)*FLOAT(NTIME)-1.E0)/12.E0                                                                     
        ALX=XI/SN2                                                                                                                  
        BLX=-ALX*(FLOAT(NTIME)+1.E0)*0.5E0                                                                                          
        DO I=1,NTIME                                                                                                                
          X(I)=X(I)-ALX*FLOAT(I)-BLX                                                                                                
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * SUBTRACT TIME MEAN FROM Y                                                                                                   
C                                                                                                                                   
      IF (CROSS) THEN                                                                                                               
        YAVG=0.E0                                                                                                                   
        DO I=1,NTIME                                                                                                                
          YAVG=YAVG+Y(I)                                                                                                            
        ENDDO                                                                                                                       
        YAVG=YAVG/FLOAT(NTIME)                                                                                                      
        DO I=1,NTIME                                                                                                                
          Y(I)=Y(I)-YAVG                                                                                                            
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * REMOVE LINEAR TREND FROM Y                                                                                                
C                                                                                                                                   
        IF (LTREND.EQ.1) THEN                                                                                                       
          YI=0.E0                                                                                                                   
          DO I=1,NTIME                                                                                                              
            YI=YI+Y(I)*FLOAT(I)                                                                                                     
          ENDDO                                                                                                                     
          ALY=YI/SN2                                                                                                                
          BLY=-ALY*(FLOAT(NTIME)+1.E0)*0.5E0                                                                                        
          DO I=1,NTIME                                                                                                              
            Y(I)=Y(I)-ALY*FLOAT(I)-BLY                                                                                              
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     CHUNK LENGTH                                                                                                                  
C                                                                                                                                   
      IF (NTYPE.EQ.0.OR.NTYPE.EQ.1) THEN                                                                                            
        LCH=NTIME/LFLT                                                                                                              
        NTIME1=LCH*LFLT                                                                                                             
      ELSE                                                                                                                          
        LCH=NTIME                                                                                                                   
        NTIME1=NTIME                                                                                                                
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * COMPUTE FREQUENCIES                                                                                                         
C                                                                                                                                   
      NFR=LCH/2                                                                                                                     
      DO I=1,NFR                                                                                                                    
        FREQ(I)=FLOAT(I)/FLOAT(LCH)                                                                                                 
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * CALCULATE VARIANCE SPECTRA OF FITTED AR1 PROCESS                                                                            
C                                                                                                                                   
      IF (NTYPE.EQ.0.OR.NTYPE.EQ.1) THEN                                                                                            
C                                                                                                                                   
C       * IF THE ESTIMATOR IS CHUNK, REMOVE LFLT-CHUNK TIME MEANS WHEN                                                              
C       * CALCULATING VARIANCES AND LAG-1 COVARIANCES                                                                               
C                                                                                                                                   
        VARX=0.E0                                                                                                                   
        VARX1=0.E0                                                                                                                  
        DO N=1,LFLT                                                                                                                 
          N0=(N-1)*LCH                                                                                                              
          XAVG=0.E0                                                                                                                 
          DO I=1,LCH                                                                                                                
            XAVG=XAVG+X(N0+I)                                                                                                       
          ENDDO                                                                                                                     
          XAVG=XAVG/FLOAT(LCH)                                                                                                      
          DO I=1,LCH                                                                                                                
            VARX=VARX+(X(N0+I)-XAVG)**2                                                                                             
          ENDDO                                                                                                                     
          DO I=1,LCH-1                                                                                                              
            VARX1=VARX1+(X(N0+I)-XAVG)*(X(N0+I+1)-XAVG)                                                                             
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
        VARX=VARX/FLOAT(NTIME1)                                                                                                     
        VARX1=VARX1/FLOAT(NTIME1)                                                                                                   
        IF(VARX.GT.0.E0)THEN                                                                                                        
          RHOX=VARX1/VARX                                                                                                           
          RHOX2=RHOX*RHOX                                                                                                           
          VARZX=VARX*(1.E0-RHOX2)                                                                                                   
          DO I=1,NFR                                                                                                                
            VARRX(I)=VARZX/(1.E0-2.E0*RHOX*COS(PI2*FREQ(I))+RHOX2)                                                                  
          ENDDO                                                                                                                     
        ELSE                                                                                                                        
          DO I=1,NFR                                                                                                                
            VARRX(I)=0.E0                                                                                                           
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
        IF (CROSS) THEN                                                                                                             
          VARY=0.E0                                                                                                                 
          VARY1=0.E0                                                                                                                
          DO N=1,LFLT                                                                                                               
            N0=(N-1)*LCH                                                                                                            
            YAVG=0.E0                                                                                                               
            DO I=1,LCH                                                                                                              
              YAVG=YAVG+Y(N0+I)                                                                                                     
            ENDDO                                                                                                                   
            YAVG=YAVG/FLOAT(LCH)                                                                                                    
            DO I=1,LCH                                                                                                              
              VARY=VARY+(Y(N0+I)-YAVG)**2                                                                                           
            ENDDO                                                                                                                   
            DO I=1,LCH-1                                                                                                            
              VARY1=VARY1+(Y(N0+I)-YAVG)*(Y(N0+I+1)-YAVG)                                                                           
            ENDDO                                                                                                                   
          ENDDO                                                                                                                     
          VARY=VARY/FLOAT(NTIME1)                                                                                                   
          VARY1=VARY1/FLOAT(NTIME1)                                                                                                 
          IF(VARY.GT.0.E0)THEN                                                                                                      
            RHOY=VARY1/VARY                                                                                                         
            RHOY2=RHOY*RHOY                                                                                                         
            VARZY=VARY*(1.E0-RHOY2)                                                                                                 
            DO I=1,NFR                                                                                                              
              VARRY(I)=VARZY/(1.E0-2.E0*RHOY*COS(PI2*FREQ(I))+RHOY2)                                                                
            ENDDO                                                                                                                   
          ELSE                                                                                                                      
            DO I=1,NFR                                                                                                              
              VARRY(I)=0.E0                                                                                                         
            ENDDO                                                                                                                   
          ENDIF                                                                                                                     
        ENDIF                                                                                                                       
      ELSE                                                                                                                          
        VARX=0.E0                                                                                                                   
        DO I=1,NTIME1                                                                                                               
          VARX=VARX+X(I)*X(I)                                                                                                       
        ENDDO                                                                                                                       
        VARX=VARX/FLOAT(NTIME1)                                                                                                     
        IF(VARX.GT.0.E0)THEN                                                                                                        
          RHOX=0.E0                                                                                                                 
          DO I=1,NTIME1-1                                                                                                           
            RHOX=RHOX+X(I)*X(I+1)                                                                                                   
          ENDDO                                                                                                                     
          RHOX=RHOX/VARX/FLOAT(NTIME1)                                                                                              
          RHOX2=RHOX*RHOX                                                                                                           
          VARZX=VARX*(1.E0-RHOX2)                                                                                                   
          DO I=1,NFR                                                                                                                
            VARRX(I)=VARZX/(1.E0-2.E0*RHOX*COS(PI2*FREQ(I))+RHOX2)                                                                  
          ENDDO                                                                                                                     
        ELSE                                                                                                                        
          DO I=1,NFR                                                                                                                
            VARRX(I)=0.E0                                                                                                           
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
        IF (CROSS) THEN                                                                                                             
          VARY=0.E0                                                                                                                 
          DO I=1,NTIME1                                                                                                             
            VARY=VARY+Y(I)*Y(I)                                                                                                     
          ENDDO                                                                                                                     
          VARY=VARY/FLOAT(NTIME1)                                                                                                   
          IF(VARX.GT.0.E0)THEN                                                                                                      
            RHOY=0.E0                                                                                                               
            DO I=1,NTIME1-1                                                                                                         
              RHOY=RHOY+Y(I)*Y(I+1)                                                                                                 
            ENDDO                                                                                                                   
            RHOY=RHOY/VARY/FLOAT(NTIME1)                                                                                            
            RHOY2=RHOY*RHOY                                                                                                         
            VARZY=VARY*(1.E0-RHOY2)                                                                                                 
            DO I=1,NFR                                                                                                              
              VARRY(I)=VARZY/(1.E0-2.E0*RHOY*COS(PI2*FREQ(I))+RHOY2)                                                                
            ENDDO                                                                                                                   
          ELSE                                                                                                                      
            DO I=1,NFR                                                                                                              
              VARRY(I)=0.E0                                                                                                         
            ENDDO                                                                                                                   
          ENDIF                                                                                                                     
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
      IF (NTYPE.EQ.0) GOTO 200                                                                                                      
C                                                                                                                                   
C     * TAPER THE DATA AND CALCULATE FOURIER COEFFICIENTS OF X AND Y                                                                
C     * (SPLIT-COSINE TAPER)                                                                                                        
C                                                                                                                                   
      NTAPER=LCH*TAPER                                                                                                              
      IF (NTAPER.NE.0) FTAPER=1.E0/FLOAT(NTAPER)                                                                                    
      IF(NR.EQ.1)WRITE(6,'(A,I5,A)')' TAPER', NTAPER,                                                                               
     1     ' POINTS AT EACH END.'                                                                                                   
      IF (NTYPE.EQ.1) THEN                                                                                                          
C                                                                                                                                   
C       * TAPER EACH CHUNK                                                                                                          
C                                                                                                                                   
        DO N=1,LFLT                                                                                                                 
          N0=(N-1)*LCH                                                                                                              
          N1=(N-1)*NFR                                                                                                              
          DO I=1,NTAPER                                                                                                             
            X(N0+I)=X(N0+I)*0.5E0*                                                                                                  
     1           (1.E0-COS(PI*(FLOAT(I)-0.5E0)*FTAPER))                                                                             
          ENDDO                                                                                                                     
          DO I=LCH-NTAPER+1,LCH                                                                                                     
            X(N0+I)=X(N0+I)*0.5E0*                                                                                                  
     1           (1.E0-COS(PI*(FLOAT(LCH-I)+0.5E0)*FTAPER))                                                                         
          ENDDO                                                                                                                     
          IF (CROSS) THEN                                                                                                           
            DO I=1,NTAPER                                                                                                           
              Y(N0+I)=Y(N0+I)*0.5E0*                                                                                                
     1             (1.E0-COS(PI*(FLOAT(I)-0.5E0)*FTAPER))                                                                           
            ENDDO                                                                                                                   
            DO I=LCH-NTAPER+1,LCH                                                                                                   
              Y(N0+I)=Y(N0+I)*0.5E0*                                                                                                
     1             (1.E0-COS(PI*(FLOAT(LCH-I)+0.5E0)*FTAPER))                                                                       
            ENDDO                                                                                                                   
          ENDIF                                                                                                                     
C                                                                                                                                   
C         * (SLOW) FOURIER TRANSFORM                                                                                                
C                                                                                                                                   
          DO I=1,NFR                                                                                                                
            AX(N1+I)=0.E0                                                                                                           
            BX(N1+I)=0.E0                                                                                                           
            FR=PI2*FREQ(I)                                                                                                          
            DO J=1,LCH                                                                                                              
              ARG=FLOAT(J)*FR                                                                                                       
              COSA=COS(ARG)                                                                                                         
              SINA=SIN(ARG)                                                                                                         
              AX(N1+I)=AX(N1+I)+X(N0+J)*COSA                                                                                        
              BX(N1+I)=BX(N1+I)+X(N0+J)*SINA                                                                                        
            ENDDO                                                                                                                   
          ENDDO                                                                                                                     
          IF (CROSS) THEN                                                                                                           
            DO I=1,NFR                                                                                                              
              AY(N1+I)=0.E0                                                                                                         
              BY(N1+I)=0.E0                                                                                                         
              FR=PI2*FREQ(I)                                                                                                        
              DO J=1,LCH                                                                                                            
                ARG=FLOAT(J)*FR                                                                                                     
                COSA=COS(ARG)                                                                                                       
                SINA=SIN(ARG)                                                                                                       
                AY(N1+I)=AY(N1+I)+Y(N0+J)*COSA                                                                                      
                BY(N1+I)=BY(N1+I)+Y(N0+J)*SINA                                                                                      
              ENDDO                                                                                                                 
            ENDDO                                                                                                                   
          ENDIF                                                                                                                     
        ENDDO                                                                                                                       
      ELSE                                                                                                                          
C                                                                                                                                   
C       * TAPER TIME SERIES ENDS                                                                                                    
C                                                                                                                                   
        DO I=1,NTAPER                                                                                                               
          X(I)=X(I)*0.5E0*(1.E0-COS(PI*(FLOAT(I)-0.5E0)*FTAPER))                                                                    
        ENDDO                                                                                                                       
        DO I=LCH-NTAPER+1,LCH                                                                                                       
          X(I)=X(I)*0.5E0*(1.E0-COS(PI*(FLOAT(LCH-I)+0.5E0)*FTAPER))                                                                
        ENDDO                                                                                                                       
        IF (CROSS) THEN                                                                                                             
          DO I=1,NTAPER                                                                                                             
            Y(I)=Y(I)*0.5E0*(1.E0-COS(PI*(FLOAT(I)-0.5E0)*FTAPER))                                                                  
          ENDDO                                                                                                                     
          DO I=LCH-NTAPER+1,LCH                                                                                                     
            Y(I)=Y(I)*0.5E0*(1.E0-COS(PI*(FLOAT(LCH-I)+0.5E0)*FTAPER))                                                              
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * (SLOW) FOURIER TRANSFORM                                                                                                  
C                                                                                                                                   
        DO I=1,NFR                                                                                                                  
          AX(I)=0.E0                                                                                                                
          BX(I)=0.E0                                                                                                                
          FR=PI2*FREQ(I)                                                                                                            
          DO J=1,LCH                                                                                                                
            ARG=FLOAT(J)*FR                                                                                                         
            COSA=COS(ARG)                                                                                                           
            SINA=SIN(ARG)                                                                                                           
            AX(I)=AX(I)+X(J)*COSA                                                                                                   
            BX(I)=BX(I)+X(J)*SINA                                                                                                   
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
        IF (CROSS) THEN                                                                                                             
          DO I=1,NFR                                                                                                                
            AY(I)=0.E0                                                                                                              
            BY(I)=0.E0                                                                                                              
            FR=PI2*FREQ(I)                                                                                                          
            DO J=1,LCH                                                                                                              
              ARG=FLOAT(J)*FR                                                                                                       
              COSA=COS(ARG)                                                                                                         
              SINA=SIN(ARG)                                                                                                         
              AY(I)=AY(I)+Y(J)*COSA                                                                                                 
              BY(I)=BY(I)+Y(J)*SINA                                                                                                 
            ENDDO                                                                                                                   
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * CALCULATE FILTER KERNELS                                                                                                    
C                                                                                                                                   
      IF (NTYPE.EQ.1.OR.NTYPE.EQ.2) THEN                                                                                            
C                                                                                                                                   
C       * KERNEL OF DANIELL ESTIMATOR (MOVING AVERAGE) AND OF                                                                       
C       * CHUNK ESTIMATOR                                                                                                           
C                                                                                                                                   
        DO I=-NTIME/2,NTIME/2                                                                                                       
          FILT(I)=0.E0                                                                                                              
        ENDDO                                                                                                                       
        FILT(0)=1.E0/FLOAT(LFLT)                                                                                                    
        DO I=-(LFLT-1)/2,LFLT/2                                                                                                     
          FILT(I)=FILT(0)                                                                                                           
        ENDDO                                                                                                                       
      ELSE IF(NTYPE.EQ.-2) THEN                                                                                                     
C                                                                                                                                   
C       * KERNEL OF MODIFIED DANIELL ESTIMATOR WITH HALF WEIGHT                                                                     
C       * TO THE END VALUES (AS IN S-PLUS AND R PACKAGES).                                                                          
C                                                                                                                                   
        DO I=-NTIME/2,NTIME/2                                                                                                       
          FILT(I)=0.E0                                                                                                              
        ENDDO                                                                                                                       
        FILT(0)=1.E0/FLOAT(LFLT-1)                                                                                                  
        DO I=-(LFLT-1)/2+1,LFLT/2-1                                                                                                 
          FILT(I)=FILT(0)                                                                                                           
        ENDDO                                                                                                                       
        FILT(-(LFLT-1)/2)=FILT(0)*0.5E0                                                                                             
        FILT(LFLT/2)=FILT(0)*0.5E0                                                                                                  
      ELSE IF (NTYPE.EQ.3) THEN                                                                                                     
C                                                                                                                                   
C       * KERNEL OF BARTLETT ESTIMATOR                                                                                              
C                                                                                                                                   
        MCUT=1.5E0*FLOAT(NTIME)/FLOAT(LFLT)                                                                                         
        FILT(0)=FLOAT(MCUT)/FLOAT(NTIME)                                                                                            
        DO I=1,NTIME/2                                                                                                              
          FILT(I)=FILT(0)*(SIN(PI*FREQ(I)*FLOAT(MCUT))                                                                              
     1         /(PI*FREQ(I)*FLOAT(MCUT)))**2                                                                                        
          FILT(-I)=FILT(I)                                                                                                          
        ENDDO                                                                                                                       
      ELSE IF (NTYPE.EQ.4) THEN                                                                                                     
C                                                                                                                                   
C       * COMPUTE KERNEL FOR PARZEN ESTIMATOR                                                                                       
C                                                                                                                                   
        MCUT=1.86E0*FLOAT(NTIME)/FLOAT(LFLT)                                                                                        
        FILT(0)=0.75E0*FLOAT(MCUT)/FLOAT(NTIME)                                                                                     
        DO I=1,NTIME/2                                                                                                              
          FILT(I)=FILT(0)*(SIN(PI*FREQ(I)*FLOAT(MCUT)/2.E0)                                                                         
     1         /(PI*FREQ(I)*FLOAT(MCUT)/2.E0))**4                                                                                   
          FILT(-I)=FILT(I)                                                                                                          
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     RE-NORMALIZE KERNEL (MUST SUM TO 1)                                                                                           
C                                                                                                                                   
      FILTS=0.E0                                                                                                                    
      DO I=-NTIME/2,NTIME/2                                                                                                         
        FILTS=FILTS+FILT(I)                                                                                                         
      ENDDO                                                                                                                         
      DO I=-NTIME/2,NTIME/2                                                                                                         
        FILT(I)=FILT(I)/FILTS                                                                                                       
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * CALCULATE EQUIVALENT DEGREES OF FREEDOM                                                                                     
C     * (AS IN 'DF.KERNEL' FROM R-PROJECT)                                                                                          
C                                                                                                                                   
      EDF=0.E0                                                                                                                      
      DO I=-NTIME/2,NTIME/2                                                                                                         
        EDF=EDF+FILT(I)**2                                                                                                          
      ENDDO                                                                                                                         
      EDF=2.E0/EDF                                                                                                                  
C                                                                                                                                   
C     * ADJUST EDF FOR TAPERING (AS IN 'SPEC.TAPER' FROM R-PROJECT)                                                                 
C                                                                                                                                   
      U2=(1.E0-(5.E0/8.E0)*TAPER*2.E0)                                                                                              
      U4=(1.E0-(93.E0/128.E0)*TAPER*2.E0)                                                                                           
      EDF=EDF/(U4/U2**2)                                                                                                            
C                                                                                                                                   
C     * CALCULATE EQUIVALENT BADNWIDTH                                                                                              
C     * (AS IN VON STORCH AND ZWIERS)                                                                                               
C                                                                                                                                   
      EBW=EDF/2.E0/FLOAT(NTIME1)                                                                                                    
C                                                                                                                                   
C     * CALCULATE EQUIVALENT BADNWIDTH                                                                                              
C     * (AS IN 'BANDWIDTH.KERNEL' FROM R-PROJECT)                                                                                   
C                                                                                                                                   
C      EBW=0.E0                                                                                                                     
C      DO I=-NTIME/2,NTIME/2                                                                                                        
C         EBW=EBW+(1.E0/12.E0+FLOAT(I)**2)*FILT(I)                                                                                  
C      ENDDO                                                                                                                        
C      EBW=SQRT(EBW)/FLOAT(NTIME1)                                                                                                  
C                                                                                                                                   
C     * CALCULATE CRITICAL VALUES                                                                                                   
C                                                                                                                                   
      P=CONF/100.E0                                                                                                                 
      P1=(1.E0+P)/2.E0                                                                                                              
      P2=(1.E0-P)/2.E0                                                                                                              
      VARL=EDF/CVCHI2D(P1,EDF)                                                                                                      
      VARU=EDF/CVCHI2D(P2,EDF)                                                                                                      
      STND=CVSTND(P1)                                                                                                               
      IF (LFLT.GT.1) THEN                                                                                                           
        FD=CVFD(P,2.E0,EDF-2.E0)                                                                                                    
        TD=CVTD(P1,EDF-2.E0)                                                                                                        
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * COMPUTE PERIODOGRAMS                                                                                                        
C                                                                                                                                   
      IF (NTYPE.EQ.1) THEN                                                                                                          
C                                                                                                                                   
C       * FOR CHUNK ESTIMATOR                                                                                                       
C                                                                                                                                   
        PERXS=0.E0                                                                                                                  
        DO N=1,LFLT                                                                                                                 
          DO I=(N-1)*NFR+1,N*NFR                                                                                                    
            PERXX(I)=(AX(I)*AX(I)+BX(I)*BX(I))/FLOAT(LCH)                                                                           
            PERXS=PERXS+2.E0*PERXX(I)                                                                                               
          ENDDO                                                                                                                     
          IF(MOD(LCH,2).EQ.0)PERXS=PERXS-PERXX(N*NFR)                                                                               
        ENDDO                                                                                                                       
        PERXS=PERXS/FLOAT(NTIME1)                                                                                                   
        IF (CROSS) THEN                                                                                                             
          PERYS=0.0                                                                                                                 
          DO N=1,LFLT                                                                                                               
            DO I=(N-1)*NFR+1,N*NFR                                                                                                  
              PERYY(I)=(AY(I)*AY(I)+BY(I)*BY(I))/FLOAT(LCH)                                                                         
              PERYS=PERYS+2.E0*PERYY(I)                                                                                             
              IF(MOD(LCH,2).EQ.0)PERYS=PERYS-PERYY(N*NFR)                                                                           
              PERXY(I)=(AX(I)*AY(I)+BX(I)*BY(I))/FLOAT(LCH)                                                                         
              PERYX(I)=(AY(I)*BX(I)-AX(I)*BY(I))/FLOAT(LCH)                                                                         
            ENDDO                                                                                                                   
          ENDDO                                                                                                                     
          PERYS=PERYS/FLOAT(NTIME1)                                                                                                 
        ENDIF                                                                                                                       
      ELSE                                                                                                                          
C                                                                                                                                   
C       * FOR ALL OTHER ESTIMATORS                                                                                                  
C                                                                                                                                   
        PERXS=0.E0                                                                                                                  
        DO I=1,NFR                                                                                                                  
          PERXX(I)=(AX(I)*AX(I)+BX(I)*BX(I))/FLOAT(LCH)                                                                             
          PERXS=PERXS+2.E0*PERXX(I)                                                                                                 
        ENDDO                                                                                                                       
        PERXX(0)=PERXX(1)                                                                                                           
        IF(MOD(LCH,2).EQ.0)PERXS=PERXS-PERXX(NFR)                                                                                   
        PERXS=PERXS/FLOAT(NTIME)                                                                                                    
        IF (CROSS) THEN                                                                                                             
          PERYS=0.E0                                                                                                                
          DO I=1,NFR                                                                                                                
            PERYY(I)=(AY(I)*AY(I)+BY(I)*BY(I))/FLOAT(LCH)                                                                           
            PERYS=PERYS+2.E0*PERYY(I)                                                                                               
            PERXY(I)=(AX(I)*AY(I)+BX(I)*BY(I))/FLOAT(LCH)                                                                           
            PERYX(I)=(AY(I)*BX(I)-AX(I)*BY(I))/FLOAT(LCH)                                                                           
          ENDDO                                                                                                                     
          PERYY(0)=PERYY(1)                                                                                                         
          PERXY(0)=PERXY(1)                                                                                                         
          PERYX(0)=PERYX(1)                                                                                                         
          IF(MOD(LCH,2).EQ.0)PERYS=PERYS-PERYY(NFR)                                                                                 
          PERYS=PERYS/FLOAT(NTIME)                                                                                                  
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * SMOOTH THE PERIODOGRAMS                                                                                                     
C                                                                                                                                   
      IF (NTYPE.EQ.1) THEN                                                                                                          
C                                                                                                                                   
C       * CHUNK ESTIMATOR                                                                                                           
C                                                                                                                                   
        DO I=1,NFR                                                                                                                  
          VARIX(I)=0.E0                                                                                                             
          DO N=1,LFLT                                                                                                               
            VARIX(I)=VARIX(I)+PERXX((N-1)*NFR+I)                                                                                    
          ENDDO                                                                                                                     
          VARIX(I)=VARIX(I)/FLOAT(LFLT)                                                                                             
        ENDDO                                                                                                                       
        IF (CROSS) THEN                                                                                                             
          DO I=1,NFR                                                                                                                
            VARIY(I)=0.E0                                                                                                           
            COHER(I)=0.E0                                                                                                           
            PHASE(I)=0.E0                                                                                                           
            DO N=1,LFLT                                                                                                             
              VARIY(I)=VARIY(I)+PERYY((N-1)*NFR+I)                                                                                  
              COHER(I)=COHER(I)+PERXY((N-1)*NFR+I)                                                                                  
              PHASE(I)=PHASE(I)+PERYX((N-1)*NFR+I)                                                                                  
            ENDDO                                                                                                                   
            VARIY(I)=VARIY(I)/FLOAT(LFLT)                                                                                           
            COHER(I)=COHER(I)/FLOAT(LFLT)                                                                                           
            PHASE(I)=PHASE(I)/FLOAT(LFLT)                                                                                           
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
      ELSE                                                                                                                          
C                                                                                                                                   
C       * APPLY FILTER TO PERIODOGRAMS                                                                                              
C                                                                                                                                   
        DO I=1,NFR                                                                                                                  
          VARIX(I)=0.E0                                                                                                             
          DO JJ=-NFR,NFR                                                                                                            
            J=ABS(I+JJ)                                                                                                             
            IF (J .GT. NFR) J=LCH-J                                                                                                 
            VARIX(I)=VARIX(I)+PERXX(J)*FILT(JJ)                                                                                     
          ENDDO                                                                                                                     
        ENDDO                                                                                                                       
        IF (CROSS) THEN                                                                                                             
          DO I=1,NFR                                                                                                                
            VARIY(I)=0.E0                                                                                                           
            COHER(I)=0.E0                                                                                                           
            PHASE(I)=0.E0                                                                                                           
            DO JJ=-NFR,NFR                                                                                                          
              J=ABS(I+JJ)                                                                                                           
              IF (J .GT. NFR) J=LCH-J                                                                                               
              VARIY(I)=VARIY(I)+PERYY(J)*FILT(JJ)                                                                                   
              COHER(I)=COHER(I)+PERXY(J)*FILT(JJ)                                                                                   
              PHASE(I)=PHASE(I)+PERYX(J)*FILT(JJ)                                                                                   
            ENDDO                                                                                                                   
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * ADJUST VARIANCE SPECTRA FOR TAPERING                                                                                        
C                                                                                                                                   
      DO I=1,NFR                                                                                                                    
        VARIX(I)=VARIX(I)/U2                                                                                                        
      ENDDO                                                                                                                         
      IF (CROSS) THEN                                                                                                               
        DO I=1,NFR                                                                                                                  
          VARIY(I)=VARIY(I)/U2                                                                                                      
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * COMPLETE CALCULATIONS OF COHERENCE AND PHASE SPECTRA                                                                        
C                                                                                                                                   
      IF (CROSS) THEN                                                                                                               
        DO I=1,NFR                                                                                                                  
          COHE=COHER(I)                                                                                                             
          QUAD=PHASE(I)                                                                                                             
          PHASE(I)=180.E0/PI*ATAN2(-QUAD,COHE)                                                                                      
          COHER(I)=(COHE*COHE+QUAD*QUAD)/(VARIX(I)*VARIY(I))                                                                        
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * PRINT SOME DIAGNOSTICS                                                                                                      
C                                                                                                                                   
      IF (NR.EQ.1) THEN                                                                                                             
        IF (NTYPE.EQ.1) THEN                                                                                                        
          WRITE(6,6010) LFLT,LCH                                                                                                    
        ELSE IF (ABS(NTYPE).EQ.2) THEN                                                                                              
          WRITE(6,6012) LFLT                                                                                                        
        ELSE IF (NTYPE.EQ.3) THEN                                                                                                   
          WRITE(6,6014) LFLT,MCUT                                                                                                   
        ELSE IF (NTYPE.EQ.4) THEN                                                                                                   
          WRITE(6,6016) LFLT,MCUT                                                                                                   
        ENDIF                                                                                                                       
        WRITE(6,6018) EDF, EBW                                                                                                      
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * CALCULATE CONFIDENCE INTERVALS FOR VARIANCE SPECTRA                                                                         
C                                                                                                                                   
      IF (NR.EQ.1) WRITE (6,6020) CONF,VARL,VARU                                                                                    
      DO N=1,NFR                                                                                                                    
        VARIX1(N)=VARIX(N)*VARL                                                                                                     
        VARIX2(N)=VARIX(N)*VARU                                                                                                     
        VARRX1(N)=VARRX(N)/VARU                                                                                                     
        VARRX2(N)=VARRX(N)/VARL                                                                                                     
      ENDDO                                                                                                                         
      IF (CROSS) THEN                                                                                                               
        DO N=1,NFR                                                                                                                  
          VARIY1(N)=VARIY(N)*VARL                                                                                                   
          VARIY2(N)=VARIY(N)*VARU                                                                                                   
          VARRY1(N)=VARRY(N)/VARU                                                                                                   
          VARRY2(N)=VARRY(N)/VARL                                                                                                   
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * SIGNIFICANCE TEST FOR COHERENCE.                                                                                          
C       * H0: COHERENCE =0 VS. HA: COHERENCE > 0                                                                                    
C                                                                                                                                   
        IF (LFLT.GT.1) THEN                                                                                                         
          COHER0(1)=2.E0*FD/(2.E0*FD+EDF-2.E0)                                                                                      
        ELSE                                                                                                                        
          COHER0(1)=1.E0                                                                                                            
        ENDIF                                                                                                                       
        DO N=2,NFR                                                                                                                  
          COHER0(N)=COHER0(1)                                                                                                       
        ENDDO                                                                                                                       
        IF (NR.EQ.1) WRITE (6,6030) 100.E0-CONF,COHER0(1)                                                                           
C                                                                                                                                   
C       * CALCULATE CONFIDENCE INTERVALS FOR THE COHERENCE SPECTRUM                                                                 
C       * (MAKES SENSE ONLY IF THE NULL HYPOTHESIS COH=0 IS REJECTED)                                                               
C                                                                                                                                   
        DO N=1,NFR                                                                                                                  
          COH1=MIN(1.E0-EPS,MAX(0.E0,SQRT(COHER(N))-STND/SQRT(EDF)))                                                                
          COH2=MIN(1.E0-EPS,SQRT(COHER(N))+STND/SQRT(EDF))                                                                          
          COHER1(N)=(TANH(0.5E0*LOG((1.E0+COH1)/(1.E0-COH1))))**2                                                                   
          COHER2(N)=(TANH(0.5E0*LOG((1.E0+COH2)/(1.E0-COH2))))**2                                                                   
        ENDDO                                                                                                                       
C                                                                                                                                   
C       * CALCULATE CONFIDENCE INTERVALS FOR THE PHASE SPECTRUM                                                                     
C                                                                                                                                   
        IF (LFLT.GT.1) THEN                                                                                                         
          DO N=1,NFR                                                                                                                
            ARG=(180.E0/PI)*                                                                                                        
     1           ASIN(MIN(1.E0,((1.E0/COHER(N))-1.E0)*(TD/(EDF-2.E0))))                                                             
            PHASE1(N)=PHASE(N)-ARG                                                                                                  
            PHASE2(N)=PHASE(N)+ARG                                                                                                  
          ENDDO                                                                                                                     
        ELSE                                                                                                                        
          DO N=1,NFR                                                                                                                
            PHASE1(N)=PHASE(N)-180.E0                                                                                               
            PHASE2(N)=PHASE(N)+180.E0                                                                                               
          ENDDO                                                                                                                     
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
                                                                                                                                    
      WRITE(6,6035) NR,VARX,PERXS                                                                                                   
      IF (CROSS) WRITE(6,6036) VARY,PERYS                                                                                           
C                                                                                                                                   
C     * PRINT ESTIMATES (ONE PER BANDWIDTH)                                                                                         
C                                                                                                                                   
      WRITE(6,'(A)')                                                                                                                
     1     '  PRINT OUT EVERY LFLT-TH VALUE:'                                                                                       
      IF (CROSS) THEN                                                                                                               
        WRITE(6,'(A,F6.2,A)')                                                                                                       
     1       '  (COHERENCE SIGNIFICANT AT THE',100.E0-CONF,                                                                         
     2       '% LEVEL IS MARKED BY AN ASTERISK)'                                                                                    
        WRITE(6,6040) (P2*100.E0,P1*100.E0,I=1,2)                                                                                   
      ELSE                                                                                                                          
        WRITE(6,6045) P2*100.E0,P1*100.E0                                                                                           
      ENDIF                                                                                                                         
      IF (NTYPE.EQ.1) THEN                                                                                                          
        I1=1                                                                                                                        
        I2=NFR                                                                                                                      
        ID=1                                                                                                                        
      ELSE                                                                                                                          
        I1=LFLT                                                                                                                     
        I2=NFR                                                                                                                      
        ID=LFLT                                                                                                                     
      ENDIF                                                                                                                         
      IF (CROSS) THEN                                                                                                               
        DO I=I1,I2,ID                                                                                                               
          IF (COHER(I).GE.COHER0(1)+EPS) THEN                                                                                       
            SGNF='*'                                                                                                                
          ELSE                                                                                                                      
            SGNF=' '                                                                                                                
          ENDIF                                                                                                                     
          WRITE(6,6050) I,FREQ(I),1.E0/FREQ(I),                                                                                     
     1         VARIX1(I),VARIX(I),VARIX2(I),                                                                                        
     2         VARIY1(I),VARIY(I),VARIY2(I),                                                                                        
     3         COHER(I),SGNF,PHASE(I)                                                                                               
        ENDDO                                                                                                                       
      ELSE                                                                                                                          
        DO I=I1,I2,ID                                                                                                               
          WRITE(6,6050) I,FREQ(I),1.E0/FREQ(I),                                                                                     
     1         VARIX1(I),VARIX(I),VARIX2(I)                                                                                         
        ENDDO                                                                                                                       
      ENDIF                                                                                                                         
                                                                                                                                    
 200  CONTINUE                                                                                                                      
      WRITE(6,6060) RHOX                                                                                                            
      IF(CROSS)WRITE(6,6060) RHOY                                                                                                   
C                                                                                                                                   
C     * SAVE RESULTS                                                                                                                
C                                                                                                                                   
      IBUF(2)=NR                                                                                                                    
      IBUF(5)=NFR                                                                                                                   
      IBUF(8)=1                                                                                                                     
C                                                                                                                                   
C     * FREQUENCIES                                                                                                                 
C                                                                                                                                   
      IBUF(3)=NC4TO8("OMEG")                                                                                                        
      CALL PUTFLD2(NFOUT,FREQ,IBUF,MAXX)                                                                                            
C                                                                                                                                   
C     * VARIANCE SPECTRUM OF X                                                                                                      
C                                                                                                                                   
      IF (NTYPE.NE.0) THEN                                                                                                          
        IBUF(3)=NC4TO8("SPEX")                                                                                                      
        CALL PUTFLD2(NFOUT,VARIX,IBUF,MAXX)                                                                                         
C                                                                                                                                   
C       * LOWER CONFIDENCE BOUND OF X SPECTRUM                                                                                      
C                                                                                                                                   
        IBUF(3)=NC4TO8("SPX1")                                                                                                      
        CALL PUTFLD2(NFOUT,VARIX1,IBUF,MAXX)                                                                                        
C                                                                                                                                   
C       * UPPER CONFIDENCE BOUND OF X SPECTRUM                                                                                      
C                                                                                                                                   
        IBUF(3)=NC4TO8("SPX2")                                                                                                      
        CALL PUTFLD2(NFOUT,VARIX2,IBUF,MAXX)                                                                                        
C                                                                                                                                   
C       * VARIANCE SPECTRUM OF FITTED AR1 PROCESS FOR X                                                                             
C                                                                                                                                   
      ENDIF                                                                                                                         
      IBUF(3)=NC4TO8("SPRX")                                                                                                        
      CALL PUTFLD2(NFOUT,VARRX,IBUF,MAXX)                                                                                           
      IBUF(3)=NC4TO8("SRX1")                                                                                                        
      CALL PUTFLD2(NFOUT,VARRX1,IBUF,MAXX)                                                                                          
      IBUF(3)=NC4TO8("SRX2")                                                                                                        
      CALL PUTFLD2(NFOUT,VARRX2,IBUF,MAXX)                                                                                          
C                                                                                                                                   
C     * VARIANCE SPECTRUM OF Y                                                                                                      
C                                                                                                                                   
      IF (CROSS) THEN                                                                                                               
        IF (NTYPE.NE.0) THEN                                                                                                        
          IBUF(3)=NC4TO8("SPEY")                                                                                                    
          CALL PUTFLD2(NFOUT,VARIY,IBUF,MAXX)                                                                                       
C                                                                                                                                   
C         * LOWER CONFIDENCE BOUND OF Y SPECTRUM                                                                                    
C                                                                                                                                   
          IBUF(3)=NC4TO8("SPY1")                                                                                                    
          CALL PUTFLD2(NFOUT,VARIY1,IBUF,MAXX)                                                                                      
C                                                                                                                                   
C         * UPPER CONFIDENCE BOUND OF Y SPECTRUM                                                                                    
C                                                                                                                                   
          IBUF(3)=NC4TO8("SPY2")                                                                                                    
          CALL PUTFLD2(NFOUT,VARIY2,IBUF,MAXX)                                                                                      
        ENDIF                                                                                                                       
C                                                                                                                                   
C       * VARIANCE SPECTRUM OF FITTED AR1 PROCESS FOR Y                                                                             
C                                                                                                                                   
        IBUF(3)=NC4TO8("SPRY")                                                                                                      
        CALL PUTFLD2(NFOUT,VARRY,IBUF,MAXX)                                                                                         
        IBUF(3)=NC4TO8("SRY1")                                                                                                      
        CALL PUTFLD2(NFOUT,VARRY1,IBUF,MAXX)                                                                                        
        IBUF(3)=NC4TO8("SRY2")                                                                                                      
        CALL PUTFLD2(NFOUT,VARRY2,IBUF,MAXX)                                                                                        
        IF (NTYPE.NE.0) THEN                                                                                                        
C                                                                                                                                   
C         * SQUARED COHERENCE SPECTRUM                                                                                              
C                                                                                                                                   
          IBUF(3)=NC4TO8("COHE")                                                                                                    
          CALL PUTFLD2(NFOUT,COHER,IBUF,MAXX)                                                                                       
C                                                                                                                                   
C         * NON SIGNIFICANT COHERENCE                                                                                               
C                                                                                                                                   
          IBUF(3)=NC4TO8("COH0")                                                                                                    
          CALL PUTFLD2(NFOUT,COHER0,IBUF,MAXX)                                                                                      
C                                                                                                                                   
C         * LOWER CONFIDENCE BOUND OF COHERENCE                                                                                     
C                                                                                                                                   
          IBUF(3)=NC4TO8("COH1")                                                                                                    
          CALL PUTFLD2(NFOUT,COHER1,IBUF,MAXX)                                                                                      
C                                                                                                                                   
C         * UPPER CONFIDENCE BOUND OF COHERENCE                                                                                     
C                                                                                                                                   
          IBUF(3)=NC4TO8("COH2")                                                                                                    
          CALL PUTFLD2(NFOUT,COHER2,IBUF,MAXX)                                                                                      
C                                                                                                                                   
C         * PHASE LAG SPECTRUM                                                                                                      
C                                                                                                                                   
          IBUF(3)=NC4TO8("PHAS")                                                                                                    
          CALL PUTFLD2(NFOUT,PHASE,IBUF,MAXX)                                                                                       
C                                                                                                                                   
C         * LOWER CONFIDENCE BOUND OF PHASE LAG SPECTRUM                                                                            
C                                                                                                                                   
          IBUF(3)=NC4TO8("PHA1")                                                                                                    
          CALL PUTFLD2(NFOUT,PHASE1,IBUF,MAXX)                                                                                      
C                                                                                                                                   
C         * UPPER CONFIDENCE BOUND OF PHASE LAG SPECTRUM                                                                            
C                                                                                                                                   
          IBUF(3)=NC4TO8("PHA2")                                                                                                    
          CALL PUTFLD2(NFOUT,PHASE2,IBUF,MAXX)                                                                                      
        ENDIF                                                                                                                       
      ENDIF                                                                                                                         
C                                                                                                                                   
C     * REPEAT FOR THE NEXT RECORD                                                                                                  
C                                                                                                                                   
      GOTO 100                                                                                                                      
                                                                                                                                    
 900  CALL                                         XIT('CROSPEC',-9)                                                                
                                                                                                                                    
C-------------------------------------------------------------------------                                                          
 5010 FORMAT (10X,2I10,1E10.0,I5,1E10.0)                                        H4                                                  
 6000 FORMAT (                                                                                                                      
     1     /' INPUT PARAMETERS: NTYPE=',I2,' LFLT=',I5,                                                                             
     2      ' CONF. INTERVAL=',F8.2,'% LTREND=',I5,' TAPER=',F6.3)                                                                  
 6010 FORMAT (                                                                                                                      
     1     /' VARIANCE, COHER^2 AND PHASE SPECTRA (CHUNK ESTIMATOR)'                                                                
     2     /'   NUMBER OF CHUNKS: ',I10                                                                                             
     3     /'   CHUNK LENGTH    : ',I10)                                                                                            
 6012 FORMAT (                                                                                                                      
     1     /' VARIANCE, COHER^2 AND PHASE SPECTRA (DANIELL ESTIMATOR)'                                                              
     2     /'  NUMBER OF PERIODOGRAM ORDINATES USED FOR AVERAGING LFLT='                                                            
     3     ,I10)                                                                                                                    
 6014 FORMAT (                                                                                                                      
     1     /' VARIANCE, COHER^2 AND PHASE SPECTRA (BARTLETT ESTIMATOR)'                                                             
     2     /'   LFLT=',I10,' THE CUTOFF POINT M=1.5*NTIME/LFLT =',I10)                                                              
 6016 FORMAT (                                                                                                                      
     1     /' VARIANCE, COHER^2 AND PHASE SPECTRA (PARZEN ESTIMATOR)'                                                               
     2     /'   LFLT=',I10,' THE CUTOFF POINT M=1.86*NTIME/LFLT =',I10)                                                             
 6018 FORMAT ('   EQUIVALENT DEGREES OF FREEDOM : ',1PE14.6                                                                         
     1     /'   EQUIVALENT BANDWIDTH          : ',1PE14.6)                                                                          
 6020 FORMAT (                                                                                                                      
     1     /' ',F6.2,'% CONFIDENCE INTERVALS FOR VARIANCE SPECTRUM:'                                                                
     2     /'  [',1PE10.4,'*(SPECTRUM ESTIMATE) ... ',                                                                              
     3     1PE10.4,'*(SPECTRUM ESTIMATE)]')                                                                                         
 6030 FORMAT (/'   THE ',F6.3,                                                                                                      
     1     '% LEVEL CRITIAL VALUE FOR SQUARED COHERENCE = ',F6.3)                                                                   
 6035 FORMAT (                                                                                                                      
     1     /' NR =',I6,                                                                                                             
     2     /'  VAR(X)=',1PE14.6,'  MEAN PERIODOGRAM(X)=',1PE14.6)                                                                   
 6036 FORMAT ( '  VAR(Y)=',1PE14.6,'  MEAN PERIODOGRAM(Y)=',1PE14.6)                                                                
 6040 FORMAT ('    I:   FREQ    PERIOD   ',F6.2,'%   VAR(X)  ',F6.2,'%',                                                            
     1     '   ',F6.2,'%   VAR(Y)  ',F6.2,'%   COHER  PHASE')                                                                       
 6045 FORMAT ('    I:   FREQ    PERIOD   ',F6.2,'%   VAR(X)  ',F6.2,'%')                                                            
 6050 FORMAT (I5,':',1P2E9.2,1X,1P3E9.2,1X,1P3E9.2,0P1F6.2,A1,1F7.1)                                                                
 6060 FORMAT (' LAG-1 CORRELATION(X)=',E12.6)                                                                                       
 6065 FORMAT (' LAG-1 CORRELATION(Y)=',E12.6)                                                                                       
      END                                                                                                                           
