      PROGRAM CANCOR                                                                                                                
C     PROGRAM CANCOR(X,       Y,       XEOF,       YEOF,        CCA,            H2                                                  
C    1                                             INPUT,       OUTPUT, )       H2                                                  
C    2        TAPE11=X,TAPE12=Y,TAPE13=XEOF,TAPE14=YEOF, TAPE15=CCA,                                                                
C    3                                       TAPE5=INPUT, TAPE6=OUTPUT)                                                             
C     -----------------------------------------------------------------         H2                                                  
C                                                                               H2                                                  
C     DEC 05/05 - SLAVA KHARIN (CORRECT SMALL BUG AFFECTING CAPX/CCPX USE.)     H2                                                  
C     JUL 03/03 - A.IRWIN,F.MAJAESS (REVISED FOR F90 AND LINUX PLATFORMS)                                                           
C     MAR 27/03 - SLAVA KHARIN (FIX A MISSING VALUE BUG. CORRECT NOTATION.)                                                         
C     MAR 09/99 - SLAVA KHARIN.                                                                                                     
C                                                                               H2                                                  
CCANCOR  - PERFORMS CANONICAL CORRELATION ANALYSIS.                     4  1 C  H1                                                  
C                                                                               H3                                                  
CAUTHOR  - SLAVA KHARIN                                                         H3                                                  
C                                                                               H3                                                  
CPURPOSE - PERFORMS CANONICAL CORRELATION ANALYSIS OF TWO FIELDS.               H3                                                  
C                                                                               H3                                                  
CINPUT FILES...                                                                 H3                                                  
C                                                                               H3                                                  
C     X    = INPUT FILE WITH FIELDS OF THE 1ST VARIABLE.                        H3                                                  
C     Y    = INPUT FILE WITH FIELDS OF THE 2ND VARIABLE.                        H3                                                  
C     XEOF = INPUT FILE WITH EOFS AND PRINCIPAL COMPONENTS OF THE 1ST VARIABLE  H3                                                  
C            EXACTLY THEY COME FROM THE MKEOF PROGRAM.                          H3                                                  
C            THE EOFS MUST BE CALCULATED FOR THE COVARIANCE MATRIX (ICOV=0 IN   H3                                                  
C            MKEOF.)  EOF NORMALIZTION DOES NOT MATTER.                         H3                                                  
C     YEOF = THE SAME AS XEOF BUT FOR THE 2ND VARIABLE.                         H3                                                  
C                                                                               H3                                                  
COUTPUT FILE...                                                                 H3                                                  
C                                                                               H3                                                  
C     CCA = OUTPUT FILE WITH CANONICAL PATTERNS AND TIME SERIES IN CCRN FORMAT. H3                                                  
C       1ST RECORD CONTAINS CANONICAL CORRELATIONS  (VARIABLE NAME = 4HCCOR)    H3                                                  
C       2ND RECORD CONTAINS 1ST CANONICAL COEFFICIENTS OF X  (NAME = 4HCAPX)    H3                                                  
C       3RD RECORD CONTAINS 1ST CANONICAL TIME SERIES OF X   (NAME = 4HCCTX)    H3                                                  
C                               CCTX=DOT PRODUCT <CAPX,X>                       H3                                                  
C       4TH RECORD CONTAINS 1ST CAN. EXPANSION PATTERN OF X  (NAME = 4HCCPX)    H3                                                  
C                               X=SUM(CCPX*CCTX)                                H3                                                  
C       5TH RECORD CONTAINS 1ST CANONICAL COEFFICIENTS OF Y  (NAME = 4HCAPY)    H3                                                  
C       6TH RECORD CONTAINS 1ST CANONICAL TIME SERIES OF Y   (NAME = 4HCCTY)    H3                                                  
C                               CCTY=DOT PRODUCT <CAPY,X>                       H3                                                  
C       7TH RECORD CONTAINS 1ST CAN. EXPANISON PATTERN OF Y  (NAME = 4HCCPY)    H3                                                  
C                               Y=SUM(CCPY*CCTY)                                H3                                                  
C       ...                 2ND ...                                             H3                                                  
C                                                                               H3                                                  
C   NOTE - IF YOU USE AREA WEIGHTING OF INPUT FIELDS BEFORE RUNNING             H3                                                  
C          MKEOF+CANCOR, DO NOT FORGET TO REMOVE THIS WEIGHTING FROM PATTERNS   H3                                                  
C          (CCPX, CCPY, CAPX AND CAPY) AFTERWARDS.                              H3                                                  
C                                                                                                                                   
CINPUT PARAMETERS:                                                                                                                  
C                                                                               H5                                                  
C     NEOFX   =    THE NUMBER OF EOFS FOR X.                                    H5                                                  
C             = 0  USE ALL EOFS FROM FILE EOFX.                                 H5                                                  
C     NEOFY   =    THE NUMBER OF EOFS FOR Y.                                    H5                                                  
C             = 0  USE ALL EOFS FROM FILE EOFY.                                 H5                                                  
C     SPVALX  =    REAL NUMBER USED TO INDICATE MISSING VALUES IN THE DATA      H5                                                  
C                  FIELDS OF X.                                                 H5                                                  
C     SPVALY  =    REAL NUMBER USED TO INDICATE MISSING VALUES IN THE DATA      H5                                                  
C                  FIELDS OF Y.                                                 H5                                                  
C                                                                               H5                                                  
CEXAMPLES OF INPUT CARDS...                                                     H5                                                  
C                                                                               H5                                                  
C*CANCOR      5    3    1.E+38    1.E+38                                        H5                                                  
C                                                                               H5                                                  
C (PERFORM CCA FOR THE FIRST 5 EOFS/PCS OF THE FIRST VARIABLE                   H5                                                  
C  AND FOR THE FIRST 3 EOFS/PCS OF THE SECOND VARIABLE.)                        H5                                                  
C--------------------------------------------------------------------------                                                         
C                                                                                                                                   
      IMPLICIT REAL (A-H,O-Z),                                                                                                      
     +INTEGER (I-N)                                                                                                                 
      LOGICAL OK                                                                                                                    
C                                                                                                                                   
C     * LENMAX = MAX LENGTH OF A DATA RECORD.                                                                                       
C     * NTMAX  = MAX NUMBER OF TIME STEPS.                                                                                          
C     * NMAX   = MAX DIMENSION OF "CANONICAL CORRELATION" MATRIX.                                                                   
C              THE NUMBER OF EOFS USED FOR BOTH VARIABLES MUST BE                                                                   
C              SMALLER THAN OR EQUALS TO NMAX.                                                                                      
C                                                                                                                                   
      PARAMETER (NTMAX=18528,LENMAX=65341,NMAX=100)                                                                                 
      COMMON/ICOM/IBUF(8),IDAT(65341)                                                                                               
      COMMON/JCOM/JBUF(8),JDAT(18528)                                                                                               
C                                                                                                                                   
C     * ARRAYS FOR INPUT FIELDS                                                                                                     
C                                                                                                                                   
      REAL  X(LENMAX),Y(LENMAX),EVALX(NMAX),EVALY(NMAX)                                                                             
     1     ,PCTX(NMAX),PCTY(NMAX),PCTCX(NMAX),PCTCY(NMAX)                                                                           
     2     ,TAVGX(LENMAX),TAVGY(LENMAX),VARX(NMAX),VARY(NMAX)                                                                       
     3     ,EVECX(LENMAX,NMAX),EVECY(LENMAX,NMAX)                                                                                   
     4     ,PCX(NTMAX,NMAX),PCY(NTMAX,NMAX)                                                                                         
      EQUIVALENCE (X(1),Y(1))                                                                                                       
C                                                                                                                                   
C     * ARRAYS FOR OUTPUT FIELDS                                                                                                    
C                                                                                                                                   
      REAL  CCA(NMAX),CCA0(NMAX),CCA1(NMAX),CCA2(NMAX)                                                                              
     1     ,CCPX(LENMAX),CCPY(LENMAX),CCTX(NTMAX),CCTY(NTMAX)                                                                       
     2     ,CAPX(LENMAX),CAPY(LENMAX)                                                                                               
C                                                                                                                                   
C     * SAVE SPACE BY MAKING EQUIVALENCE                                                                                            
C                                                                                                                                   
      EQUIVALENCE (CCPX(1),CCPY(1)),(CAPX(1),CAPY(1))                                                                               
     1     ,(PCX(1,1),CCTX(1)),(PCY(1,1),CCTY(1))                                                                                   
C                                                                                                                                   
C     INTERNAL ARRAYS                                                                                                               
C                                                                                                                                   
      REAL  COV(NMAX,NMAX),SNGX(NMAX,NMAX),SNGY(NMAX,NMAX)                                                                          
      INTEGER NORDX(NMAX),NORDY(NMAX)                                                                                               
C                                                                                                                                   
      DATA MAXX/65341/, MAXT/18528/, SMALL/1.E-12/,BIG/1.E+39/                                                                      
C--------------------------------------------------------------------------                                                         
C                                                                                                                                   
C     * ASSIGN FILES TO FORTRAN UNITS                                                                                               
C                                                                                                                                   
      NFF=7                                                                                                                         
      CALL JCLPNT(NFF,11,12,13,14,15,5,6)                                                                                           
      REWIND 11                                                                                                                     
      REWIND 12                                                                                                                     
      REWIND 13                                                                                                                     
      REWIND 14                                                                                                                     
      REWIND 15                                                                                                                     
C                                                                                                                                   
C     * READ INPUT PARAMETERS FROM CARD                                                                                             
C                                                                                                                                   
      READ(5,5010,END=900) NEOFX, NEOFY, SPVALX, SPVALY                         H4                                                  
C                                                                                                                                   
C     * WRITE OUT THE INPUT PARAMETERS TO STANDARD OUTPUT                                                                           
C                                                                                                                                   
      WRITE(6,6010) NEOFX,NEOFY,SPVALX,SPVALY                                                                                       
C                                                                                                                                   
C     * READ EVS, EOFS AND PCS:                                                                                                     
C                                                                                                                                   
C     * EVS                                                                                                                         
C                                                                                                                                   
      CALL GETFLD2(13,EVALX,NC4TO8("TIME"),-1,NC4TO8("EVAL"),-1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF(.NOT.OK)THEN                                                                                                               
         CALL                                      XIT('CANCOR',-1)                                                                 
      ENDIF                                                                                                                         
      CALL PRTLAB(IBUF)                                                                                                             
      NEOFMX=IBUF(5)                                                                                                                
      IF (NEOFX.EQ.0) THEN                                                                                                          
         WRITE(6,'(A)') ' NEOFX=0. USE ALL INPUT EOFS.'                                                                             
         NEOFX=NEOFMX                                                                                                               
      ENDIF                                                                                                                         
      IF (NEOFX .GT. NMAX) THEN                                                                                                     
         WRITE(6,'(A,I5,A)')                                                                                                        
     1        ' *** ERROR: NUMBER OF EOFS NEOFX=',NEOFX,                                                                            
     2        ' IS TOO LARGE.'                                                                                                      
         CALL                                      XIT('CANCOR',-2)                                                                 
      ENDIF                                                                                                                         
      IF (NEOFMX .LT. NEOFX) THEN                                                                                                   
         WRITE(6,6020) NEOFMX                                                                                                       
         NEOFX=NEOFMX                                                                                                               
      ENDIF                                                                                                                         
                                                                                                                                    
      CALL GETFLD2(14,EVALY,NC4TO8("TIME"),-1,NC4TO8("EVAL"),-1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF(.NOT.OK)THEN                                                                                                               
         CALL                                      XIT('CANCOR',-3)                                                                 
      ENDIF                                                                                                                         
      CALL PRTLAB(IBUF)                                                                                                             
      NEOFMY=IBUF(5)                                                                                                                
      IF (NEOFY.EQ.0)THEN                                                                                                           
         WRITE(6,'(A)') ' NEOFY=0. USE ALL INPUT EOFS.'                                                                             
         NEOFY=NEOFMY                                                                                                               
      ENDIF                                                                                                                         
      IF (NEOFY .GT. NMAX) THEN                                                                                                     
         WRITE(6,'(A,I5,A)')                                                                                                        
     1        ' *** ERROR: NUMBER OF EOFS NEOFY=',NEOFY,                                                                            
     2        ' IS TOO LARGE.'                                                                                                      
         CALL                                      XIT('CANCOR',-4)                                                                 
      ENDIF                                                                                                                         
      IF (NEOFMY .LT. NEOFY) THEN                                                                                                   
         WRITE(6,6025) NEOFMY                                                                                                       
         NEOFY=NEOFMY                                                                                                               
      ENDIF                                                                                                                         
                                                                                                                                    
      NCCA=MIN(NEOFX,NEOFY)                                                                                                         
C                                                                                                                                   
C     * PCT (THE PERCENTAGES OF EXPLAINED VARIANCE)                                                                                 
C                                                                                                                                   
      CALL GETFLD2(13,PCTX,NC4TO8("TIME"),-1,NC4TO8(" PCT"),-1,                                                                     
     +                                            IBUF,MAXX,OK)                                                                     
      IF(.NOT.OK)THEN                                                                                                               
         CALL                                      XIT('CANCOR',-5)                                                                 
      ENDIF                                                                                                                         
      CALL PRTLAB(IBUF)                                                                                                             
                                                                                                                                    
      CALL GETFLD2(14,PCTY,NC4TO8("TIME"),-1,NC4TO8(" PCT"),-1,                                                                     
     +                                            IBUF,MAXX,OK)                                                                     
      IF(.NOT.OK)THEN                                                                                                               
         CALL                                      XIT('CANCOR',-6)                                                                 
      ENDIF                                                                                                                         
      CALL PRTLAB(IBUF)                                                                                                             
C                                                                                                                                   
C     * PCTC (THE CUMULATIVE PERCENTAGES)                                                                                           
C                                                                                                                                   
      CALL GETFLD2(13,PCTCX,NC4TO8("TIME"),-1,NC4TO8("PCTC"),-1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF(.NOT.OK)THEN                                                                                                               
         CALL                                      XIT('CANCOR',-7)                                                                 
      ENDIF                                                                                                                         
      CALL PRTLAB(IBUF)                                                                                                             
                                                                                                                                    
      CALL GETFLD2(14,PCTCY,NC4TO8("TIME"),-1,NC4TO8("PCTC"),-1,                                                                    
     +                                             IBUF,MAXX,OK)                                                                    
      IF(.NOT.OK)THEN                                                                                                               
         CALL                                      XIT('CANCOR',-8)                                                                 
      ENDIF                                                                                                                         
      CALL PRTLAB(IBUF)                                                                                                             
C                                                                                                                                   
C     * TIME MEAN                                                                                                                   
C                                                                                                                                   
      CALL GETFLD2(13,TAVGX,-1,-1,NC4TO8("TAVG"),-1,IBUF,MAXX,OK)                                                                   
      IF(.NOT.OK)THEN                                                                                                               
         CALL                                      XIT('CANCOR',-9)                                                                 
      ENDIF                                                                                                                         
      KINDX=IBUF(1)                                                                                                                 
      LEVX=IBUF(4)                                                                                                                  
      NLGX=IBUF(5)                                                                                                                  
      NLATX=IBUF(6)                                                                                                                 
      KHEMX=IBUF(7)                                                                                                                 
      NWDSX=NLGX*NLATX                                                                                                              
      IF (KINDX.EQ.NC4TO8("SPEC") .OR.                                                                                              
     +    KINDX.EQ.NC4TO8("FOUR")     ) NWDSX=NWDSX*2                                                                               
      CALL PRTLAB(IBUF)                                                                                                             
                                                                                                                                    
      CALL GETFLD2(14,TAVGY,-1,-1,NC4TO8("TAVG"),-1,IBUF,MAXX,OK)                                                                   
      IF(.NOT.OK)THEN                                                                                                               
         CALL                                      XIT('CANCOR',-10)                                                                
      ENDIF                                                                                                                         
      KINDY=IBUF(1)                                                                                                                 
      LEVY=IBUF(4)                                                                                                                  
      NLGY=IBUF(5)                                                                                                                  
      NLATY=IBUF(6)                                                                                                                 
      KHEMY=IBUF(7)                                                                                                                 
      NWDSY=NLGY*NLATY                                                                                                              
      IF (KINDY.EQ.NC4TO8("SPEC") .OR.                                                                                              
     +    KINDY.EQ.NC4TO8("FOUR")     ) NWDSY=NWDSY*2                                                                               
      CALL PRTLAB(IBUF)                                                                                                             
C                                                                                                                                   
C     * EOFS & PCS                                                                                                                  
C                                                                                                                                   
      DO NE=1,NEOFX                                                                                                                 
         CALL GETFLD2(13,EVECX(1,NE),-1,-1,NC4TO8("EVEC"),-1,                                                                       
     +                                               IBUF,MAXX,OK)                                                                  
         IF(.NOT.OK)THEN                                                                                                            
            CALL                                   XIT('CANCOR',-11)                                                                
         ENDIF                                                                                                                      
         IF (IBUF(5) .NE. NLGX .OR. IBUF(6) .NE. NLATX)                                                                             
     1        CALL                                 XIT('CANCOR',-12)                                                                
         IF(NE .EQ.1) CALL PRTLAB(IBUF)                                                                                             
                                                                                                                                    
         CALL GETFLD2(13,PCX(1,NE),NC4TO8("TIME"),-1,NC4TO8("  PC"),                                                                
     +                                              -1,JBUF,MAXT,OK)                                                                
         IF(.NOT.OK)THEN                                                                                                            
            CALL                                   XIT('CANCOR',-13)                                                                
         ENDIF                                                                                                                      
         IF(NE .EQ.1) CALL PRTLAB(JBUF)                                                                                             
         NSTEP=JBUF(5)                                                                                                              
      ENDDO                                                                                                                         
                                                                                                                                    
      DO NE=1,NEOFY                                                                                                                 
         CALL GETFLD2(14,EVECY(1,NE),-1,-1,NC4TO8("EVEC"),-1,                                                                       
     +                                          IBUF,MAXX,OK)                                                                       
         IF(.NOT.OK)THEN                                                                                                            
            CALL                                   XIT('CANCOR',-14)                                                                
         ENDIF                                                                                                                      
         IF (IBUF(5) .NE. NLGY .OR. IBUF(6) .NE. NLATY)                                                                             
     1        CALL                                 XIT('CANCOR',-15)                                                                
         IF(NE .EQ.1) CALL PRTLAB(IBUF)                                                                                             
                                                                                                                                    
         CALL GETFLD2(14,PCY(1,NE),NC4TO8("TIME"),-1,NC4TO8("  PC"),                                                                
     +                                              -1,JBUF,MAXT,OK)                                                                
         IF(.NOT.OK)THEN                                                                                                            
            CALL                                   XIT('CANCOR',-16)                                                                
         ENDIF                                                                                                                      
         IF(NE .EQ.1) CALL PRTLAB(JBUF)                                                                                             
         IF (JBUF(5) .NE. NSTEP)                                                                                                    
     1        CALL                                 XIT('CANCOR',-17)                                                                
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * CALCULATE TOTAL VARIANCES FROM EVAL(1) AND PCT(1)                                                                           
C                                                                                                                                   
      VARTOTX=EVALX(1)/PCTX(1)                                                                                                      
      VARTOTY=EVALY(1)/PCTY(1)                                                                                                      
C                                                                                                                                   
C     * PRINT EVAL AND EXPLAINED VARIANCES                                                                                          
C                                                                                                                                   
      WRITE(6,6030) VARTOTX,NEOFX                                                                                                   
      DO NE=1,NEOFX                                                                                                                 
         WRITE(6,6035) NE,EVALX(NE),100.E0*PCTX(NE),100.E0*PCTCX(NE)                                                                
      ENDDO                                                                                                                         
      WRITE(6,6040) VARTOTY,NEOFY                                                                                                   
      DO NE=1,NEOFY                                                                                                                 
         WRITE(6,6045) NE,EVALY(NE),100.E0*PCTY(NE),100.E0*PCTCY(NE)                                                                
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * NORMALIZE EOFS AND PCS PROPERLY                                                                                             
C     * (THE EOFS MUST BE NORMALIZED BY THE SQRT OF EIGENVALUES AND                                                                 
C     *  THE PCS MUST BE STANDARDIZED BY ONE).                                                                                      
C                                                                                                                                   
      DO NE=1,NEOFX                                                                                                                 
         SUM=0.E0                                                                                                                   
         DO NT=1,NSTEP                                                                                                              
            SUM=SUM+PCX(NT,NE)**2                                                                                                   
         ENDDO                                                                                                                      
         SUM=SQRT(SUM/FLOAT(NSTEP))                                                                                                 
         DO IJ=1,NWDSX                                                                                                              
            IF (ABS(EVECX(IJ,NE)-SPVALX).GT.SMALL) THEN                                                                             
               EVECX(IJ,NE)=EVECX(IJ,NE)*SUM                                                                                        
            ENDIF                                                                                                                   
         ENDDO                                                                                                                      
         DO NT=1,NSTEP                                                                                                              
            PCX(NT,NE)=PCX(NT,NE)/SUM                                                                                               
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
                                                                                                                                    
      DO NE=1,NEOFY                                                                                                                 
         SUM=0.E0                                                                                                                   
         DO NT=1,NSTEP                                                                                                              
            SUM=SUM+PCY(NT,NE)**2                                                                                                   
         ENDDO                                                                                                                      
         SUM=SQRT(SUM/FLOAT(NSTEP))                                                                                                 
         DO IJ=1,NWDSY                                                                                                              
            IF (ABS(EVECY(IJ,NE)-SPVALY).GT.SMALL) THEN                                                                             
               EVECY(IJ,NE)=EVECY(IJ,NE)*SUM                                                                                        
            ENDIF                                                                                                                   
         ENDDO                                                                                                                      
         DO NT=1,NSTEP                                                                                                              
            PCY(NT,NE)=PCY(NT,NE)/SUM                                                                                               
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * CALCULATE CROSS-CORRELATION MATRIX <PCX,PCY>                                                                                
C                                                                                                                                   
      DO NE2=1,NEOFY                                                                                                                
         DO NE1=1,NEOFX                                                                                                             
            COV(NE1,NE2)=0.E0                                                                                                       
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
      DO NE2=1,NEOFY                                                                                                                
         DO NE1=1,NEOFX                                                                                                             
            DO NT=1,NSTEP                                                                                                           
               COV(NE1,NE2)=COV(NE1,NE2)+PCX(NT,NE1)*PCY(NT,NE2)                                                                    
            ENDDO                                                                                                                   
            COV(NE1,NE2)=COV(NE1,NE2)/FLOAT(NSTEP)                                                                                  
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * CALCULATE "SINGULAR" MATRIX FOR X                                                                                           
C                                                                                                                                   
      DO NE2=1,NEOFX                                                                                                                
         DO NE1=1,NEOFX                                                                                                             
            SNGX(NE1,NE2)=0.E0                                                                                                      
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
      DO NE2=1,NEOFX                                                                                                                
         DO NE1=1,NEOFX                                                                                                             
            DO NE3=1,NEOFY                                                                                                          
               SNGX(NE1,NE2)=SNGX(NE1,NE2)+COV(NE1,NE3)*COV(NE2,NE3)                                                                
            ENDDO                                                                                                                   
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * SOLVE EIGENVALUE PROBLEM                                                                                                    
C                                                                                                                                   
      CALL TRED2(SNGX,NEOFX,NMAX,CCA,PCTX)                                                                                          
      CALL TQLI1(CCA,PCTX,NEOFX,NMAX,SNGX)                                                                                          
C                                                                                                                                   
C     * ORDER THE EIGENVALUES FROM MAX TO MIN                                                                                       
C                                                                                                                                   
      ELAST=BIG                                                                                                                     
      DO NE=1,NEOFX                                                                                                                 
         EMAX=-BIG                                                                                                                  
         DO NE1=1,NEOFX                                                                                                             
            IF (CCA(NE1).GT.EMAX .AND. CCA(NE1).LT.ELAST) THEN                                                                      
               EMAX=CCA(NE1)                                                                                                        
               NORDX(NE)=NE1                                                                                                        
            ENDIF                                                                                                                   
         ENDDO                                                                                                                      
         ELAST=EMAX                                                                                                                 
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * CALCULATE "SINGULAR" MATRIX FOR Y                                                                                           
C                                                                                                                                   
      DO NE2=1,NEOFY                                                                                                                
         DO NE1=1,NEOFY                                                                                                             
            SNGY(NE1,NE2)=0.E0                                                                                                      
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
      DO NE2=1,NEOFY                                                                                                                
         DO NE1=1,NEOFY                                                                                                             
            DO NE3=1,NEOFX                                                                                                          
               SNGY(NE1,NE2)=SNGY(NE1,NE2)+COV(NE3,NE1)*COV(NE3,NE2)                                                                
            ENDDO                                                                                                                   
         ENDDO                                                                                                                      
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * SOLVE EIGENVALUE PROBLEM                                                                                                    
C                                                                                                                                   
      CALL TRED2(SNGY,NEOFY,NMAX,CCA,PCTY)                                                                                          
      CALL TQLI1(CCA,PCTY,NEOFY,NMAX,SNGY)                                                                                          
C                                                                                                                                   
C     * ORDER THE EIGENVALUES FROM MAX TO MIN                                                                                       
C                                                                                                                                   
      ELAST=BIG                                                                                                                     
      DO NE=1,NEOFY                                                                                                                 
         EMAX=-BIG                                                                                                                  
         DO NE1=1,NEOFY                                                                                                             
            IF (CCA(NE1).GT.EMAX .AND. CCA(NE1).LT.ELAST) THEN                                                                      
               EMAX=CCA(NE1)                                                                                                        
               NORDY(NE)=NE1                                                                                                        
            ENDIF                                                                                                                   
         ENDDO                                                                                                                      
         ELAST=EMAX                                                                                                                 
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * SAVE CANONICAL CORRELATIONS (NAME=4HCCOR)                                                                                   
C                                                                                                                                   
      DO NE=1,NCCA                                                                                                                  
         CCA(NORDY(NE))=SQRT(CCA(NORDY(NE)))                                                                                        
         CCA0(NE)=CCA(NORDY(NE))                                                                                                    
      ENDDO                                                                                                                         
      CALL SETLAB(IBUF,NC4TO8("TIME"),1,NC4TO8("CCOR"),0,NCCA,1,0,1)                                                                
      CALL PUTFLD2(15,CCA0,IBUF,MAXX)                                                                                               
                                                                                                                                    
C                                                                                                                                   
C     * TRANSFORM CANONICAL COEFFICIENTS FROM EOF SPACE TO PHYSICAL SPACE                                                           
C                                                                                                                                   
      DO NE=1,NCCA                                                                                                                  
C                                                                                                                                   
C       * CALCULATE CANONICAL COEFFICIENTS OF X                                                                                     
C                                                                                                                                   
         DO IJ=1,NWDSX                                                                                                              
            CAPX(IJ)=0.E0                                                                                                           
            DO NE1=1,NEOFX                                                                                                          
               IF (ABS(EVECX(IJ,NE1)-SPVALX).GT.SMALL) THEN                                                                         
                  CAPX(IJ)=CAPX(IJ)+SNGX(NE1,NORDX(NE))/EVALX(NE1)                                                                  
     1                 *EVECX(IJ,NE1)                                                                                               
               ELSE                                                                                                                 
                  CAPX(IJ)=SPVALX                                                                                                   
               ENDIF                                                                                                                
            ENDDO                                                                                                                   
         ENDDO                                                                                                                      
C                                                                                                                                   
C        * CALCULATE CANONICAL CANONICAL TIME SERIES = <CAPX,X>                                                                     
C                                                                                                                                   
         REWIND 11                                                                                                                  
         DO NT=1,NSTEP                                                                                                              
            CALL GETFLD2(11,X,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                             
            IF(.NOT.OK)THEN                                                                                                         
               CALL                                XIT('CANCOR',-18)                                                                
            ENDIF                                                                                                                   
            IF (IBUF(5) .NE. NLGX .OR. IBUF(6) .NE. NLATX)                                                                          
     1           CALL                              XIT('CANCOR',-19)                                                                
            CCTX(NT)=0.E0                                                                                                           
            DO IJ=1,NWDSX                                                                                                           
               IF ( ABS(CAPX(IJ)-SPVALX).GT.SMALL .AND.                                                                             
     1              ABS(   X(IJ)-SPVALX).GT.SMALL)                                                                                  
     2              CCTX(NT)=CCTX(NT)+CAPX(IJ)*(X(IJ)-TAVGX(IJ))                                                                    
            ENDDO                                                                                                                   
         ENDDO                                                                                                                      
C                                                                                                                                   
C        * CALCULATE CANONICAL EXPANSION PATTERN OF X                                                                               
C                                                                                                                                   
         DO IJ=1,NWDSX                                                                                                              
            CCPX(IJ)=0.E0                                                                                                           
            DO NE1=1,NEOFX                                                                                                          
               IF (ABS(EVECX(IJ,NE1)-SPVALX).GT.SMALL) THEN                                                                         
                  CCPX(IJ)=CCPX(IJ)+SNGX(NE1,NORDX(NE))*EVECX(IJ,NE1)                                                               
               ELSE                                                                                                                 
                  CCPX(IJ)=SPVALX                                                                                                   
               ENDIF                                                                                                                
            ENDDO                                                                                                                   
         ENDDO                                                                                                                      
C                                                                                                                                   
C        * CALCULATE VARIANCE EXPLAINED BY THE EXPANSION PATTERN OF X                                                               
C                                                                                                                                   
         VARX(NE)=0.E0                                                                                                              
         DO IJ=1,NWDSX                                                                                                              
            IF (ABS(CCPX(IJ)-SPVALX).GT.SMALL)                                                                                      
     1           VARX(NE)=VARX(NE)+CCPX(IJ)**2                                                                                      
         ENDDO                                                                                                                      
C                                                                                                                                   
C        * SAVE CANONICAL COEFFICIENTS (NAME=4HCAPX)                                                                                
C                                                                                                                                   
         CALL SETLAB(IBUF,KINDX,NE,NC4TO8("CAPX"),LEVX,NLGX,                                                                        
     +                                        NLATX,KHEMX,1)                                                                        
         CALL PUTFLD2(15,CAPX,IBUF,MAXX)                                                                                            
C                                                                                                                                   
C        * SAVE CANONICAL TIME SERIES (NAME=4HCCTX)                                                                                 
C                                                                                                                                   
         CALL SETLAB(JBUF,NC4TO8("TIME"),NE,NC4TO8("CCTX"),LEVX,                                                                    
     +                                          NSTEP,1,KHEMX,1)                                                                    
         CALL PUTFLD2(15,CCTX,JBUF,MAXT)                                                                                            
C                                                                                                                                   
C        * SAVE CANONICAL EXPANSION PATTERN (NAME=4HCCPX)                                                                           
C                                                                                                                                   
         CALL SETLAB(IBUF,KINDX,NE,NC4TO8("CCPX"),LEVX,NLGX,                                                                        
     +                                        NLATX,KHEMX,1)                                                                        
         CALL PUTFLD2(15,CCPX,IBUF,MAXX)                                                                                            
C                                                                                                                                   
C        * CALCULATE CANONICAL COEFFICIENTS OF Y                                                                                    
C                                                                                                                                   
         DO IJ=1,NWDSY                                                                                                              
            CAPY(IJ)=0.E0                                                                                                           
            DO NE1=1,NEOFY                                                                                                          
               IF (ABS(EVECY(IJ,NE1)-SPVALY).GT.SMALL) THEN                                                                         
                  CAPY(IJ)=CAPY(IJ)+SNGY(NE1,NORDY(NE))/EVALY(NE1)                                                                  
     1                 *EVECY(IJ,NE1)                                                                                               
               ELSE                                                                                                                 
                  CAPY(IJ)=SPVALY                                                                                                   
               ENDIF                                                                                                                
            ENDDO                                                                                                                   
         ENDDO                                                                                                                      
C                                                                                                                                   
C        * CALCULATE CANONICAL TIME SERIES OF Y=<Y,CAPY>                                                                            
C                                                                                                                                   
         REWIND 12                                                                                                                  
         DO NT=1,NSTEP                                                                                                              
            CALL GETFLD2(12,Y,-1,-1,-1,-1,IBUF,MAXX,OK)                                                                             
            IF(.NOT.OK)THEN                                                                                                         
               CALL                                XIT('CANCOR',-20)                                                                
            ENDIF                                                                                                                   
            IF (IBUF(5) .NE. NLGY .OR. IBUF(6) .NE. NLATY)                                                                          
     1           CALL                              XIT('CANCOR',-21)                                                                
            CCTY(NT)=0.E0                                                                                                           
            DO IJ=1,NWDSY                                                                                                           
               IF ( ABS(CAPY(IJ)-SPVALY).GT.SMALL .AND.                                                                             
     1              ABS(   Y(IJ)-SPVALY).GT.SMALL)                                                                                  
     2              CCTY(NT)=CCTY(NT)+CAPY(IJ)*(Y(IJ)-TAVGY(IJ))                                                                    
            ENDDO                                                                                                                   
         ENDDO                                                                                                                      
C                                                                                                                                   
C       * CALCULATE CANONICAL EXPANSION PATTERN OF Y                                                                                
C                                                                                                                                   
         DO IJ=1,NWDSY                                                                                                              
            CCPY(IJ)=0.E0                                                                                                           
            DO NE1=1,NEOFY                                                                                                          
               IF (ABS(EVECY(IJ,NE1)-SPVALY).GT.SMALL) THEN                                                                         
                  CCPY(IJ)=CCPY(IJ)+SNGY(NE1,NORDY(NE))*EVECY(IJ,NE1)                                                               
               ELSE                                                                                                                 
                  CCPY(IJ)=SPVALY                                                                                                   
               ENDIF                                                                                                                
            ENDDO                                                                                                                   
         ENDDO                                                                                                                      
C                                                                                                                                   
C        * CALCULATE VARIANCE EXPLAINED BY THE EXPANSION PATTERN OF Y                                                               
C                                                                                                                                   
         VARY(NE)=0.E0                                                                                                              
         DO IJ=1,NWDSY                                                                                                              
            IF (ABS(CCPY(IJ)-SPVALY).GT.SMALL)                                                                                      
     1           VARY(NE)=VARY(NE)+CCPY(IJ)**2                                                                                      
         ENDDO                                                                                                                      
C                                                                                                                                   
C        * CHANGE THE SIGN OF TIME SERIES TO INSURE POSITIVE CORRELATION                                                            
C                                                                                                                                   
         SUM=0.E0                                                                                                                   
         DO NT=1,NSTEP                                                                                                              
            SUM=SUM+CCTX(NT)*CCTY(NT)                                                                                               
         ENDDO                                                                                                                      
         SGN=SIGN(1.E0,SUM)                                                                                                         
         IF (SGN .LT. 0.E0) THEN                                                                                                    
            DO IJ=1,NWDSY                                                                                                           
               IF ( ABS(CCPY(IJ)-SPVALY).GT.SMALL) THEN                                                                             
                  CCPY(IJ)=-CCPY(IJ)                                                                                                
                  CAPY(IJ)=-CAPY(IJ)                                                                                                
               ENDIF                                                                                                                
            ENDDO                                                                                                                   
         ENDIF                                                                                                                      
         DO NT=1,NSTEP                                                                                                              
            CCTY(NT)=SGN*CCTY(NT)                                                                                                   
         ENDDO                                                                                                                      
C                                                                                                                                   
C        * SAVE CANONICAL COEFFICIENTS (NAME=4HCAPY)                                                                                
C                                                                                                                                   
         CALL SETLAB(IBUF,KINDY,NE,NC4TO8("CAPY"),LEVY,NLGY,                                                                        
     +                                        NLATY,KHEMY,1)                                                                        
         CALL PUTFLD2(15,CAPY,IBUF,MAXX)                                                                                            
C                                                                                                                                   
C        * SAVE CANONICAL TIME SERIES (NAME=4HCCTY)                                                                                 
C                                                                                                                                   
         CALL SETLAB(JBUF,NC4TO8("TIME"),NE,NC4TO8("CCTY"),                                                                         
     +                                 LEVY,NSTEP,1,KHEMY,1)                                                                        
         CALL PUTFLD2(15,CCTY,JBUF,MAXT)                                                                                            
C                                                                                                                                   
C        * SAVE CANONICAL EXPANSION PATTERN (NAME=4HCCPY)                                                                           
C                                                                                                                                   
         CALL SETLAB(IBUF,KINDY,NE,NC4TO8("CCPY"),LEVY,                                                                             
     +                              NLGY,NLATY,KHEMY,1)                                                                             
         CALL PUTFLD2(15,CCPY,IBUF,MAXX)                                                                                            
      ENDDO                                                                                                                         
C                                                                                                                                   
C     * CALCULATE UNBIASED CANONICAL CORRELATIONS AND 95% ERROR INTERVAL                                                            
C                                                                                                                                   
      CALL CCAERR(NSTEP,NEOFX,NEOFY,CCA,CCA0,CCA1,CCA2)                                                                             
C                                                                                                                                   
C     * PRINT CANONICAL CORRELATIONS                                                                                                
C                                                                                                                                   
      WRITE (6,6050) NCCA                                                                                                           
      WRITE (6,6055)                                                                                                                
      DO NE=1,NCCA                                                                                                                  
         WRITE (6,6060) NE, CCA(NORDY(NE))                                                                                          
     1        ,VARX(NE)/VARTOTX*100.E0,VARY(NE)/VARTOTY*100.E0                                                                      
     2        ,CCA2(NORDY(NE)),CCA0(NORDY(NE)),CCA1(NORDY(NE))                                                                      
      ENDDO                                                                                                                         
      CALL                                         XIT('CANCOR',0)                                                                  
 900  CALL                                         XIT('CANCOR',-22)                                                                
C-------------------------------------------------------------------------                                                          
 5010 FORMAT (10X,2I5,2E10.0)                                                   H4                                                  
 6010 FORMAT (/' INPUT PARAMETERS ARE:'/                                                                                            
     1     ' NEOFX=',I5,' NEOFY=',I5,                                                                                               
     2     ' SPVALX=',1P1E12.5,' SPVALY=',1P1E12.5/)                                                                                
 6020 FORMAT (                                                                                                                      
     1     /' THE NUMBER OF EOFS IN FILE XEOF IS SMALLER THAN NEOFX.'                                                               
     2     /' NEOFX IS SET TO ',I5/)                                                                                                
 6025 FORMAT (                                                                                                                      
     1     /' THE NUMBER OF EOFS IN FILE YEOF IS SMALLER THAN NEOFY.'                                                               
     2     /' NEOFY IS SET TO ',I5/)                                                                                                
 6030 FORMAT (/' TOTAL VARIANCE OF X: ',1P1E12.5                                                                                    
     1     /' THE FIRST ',I5,' EIGENVALUES OF X:'                                                                                   
     2     /'    NUMBER     EVAL      PCT(%)    CUM.PCT(%)')                                                                        
 6035 FORMAT (5X,I5,1P1E12.4,0P2F10.3)                                                                                              
 6040 FORMAT (/' TOTAL VARIANCE OF Y: ',1P1E12.5                                                                                    
     1     /' THE FIRST ',I5,' EIGENVALUES OF Y:'                                                                                   
     2     /'    NUMBER     EVAL      PCT(%)    CUM.PCT(%)')                                                                        
 6045 FORMAT (5X,I5,1P1E12.4,0P2F10.3)                                                                                              
 6050 FORMAT (                                                                                                                      
     1     /' FIRST ', I5, ' CANONICAL CORRELATIONS'                                                                                
     2     /' AND VARIANCES (%) EXPLAINED BY CANONICAL PATTERNS:'                                                                   
     3     /' NOTE: CANONICAL CORRELATION PATTERNS ARE NOT ORTHOGONAL'                                                              
     4     /'       AND THE EXPLAINED VARIANCES ARE NOT ADDITIVE.')                                                                 
 6055 FORMAT (                                                                                                                      
     1     /'   NR.    CAN.COR.    VAR.X(%)    VAR.Y(%) '                                                                           
     2     ,'      -95%  CAN.COR.UNBIASED  +95%')                                                                                   
 6060 FORMAT (1X,I4,0P6F12.4)                                                                                                       
      END                                                                                                                           
