      SUBROUTINE CCAERR (NSTEP,NEOFX,NEOFY,CCA,CCA0,CCA1,CCA2)
C
C     WRITTEN BY E.ZORITA, MPI FOR METEOROLOGY, HAMBURG
C     MODIFIED BY SLAVA KHARIN, MAR 09/99.
C
C     THIS SUBROUTINE CALCULATES THE UNBIASED ESTIMATION OF
C     CANONICAL CORRELATIONS AND THEIR 95% ERROR INTERVALS, BASED       
C     ON THE FORMULAS GIVEN IN W. GLYNN AND R. MUIRHEAD (1978)          
C     J. MULTIVARIATE ANALYSIS 8, 468-478:                              
C                                                                       
C     X(I)=Z(I)-1/(2NL R(I)){N1+N2-2+R(I)**2+S(1-R(I)**2)  X            
C                                                                       
C     SUM {R(J)**2/(R(I)**2-R(J)**2)} + O(N**(-2))                      
C     I=/J                                                              
C                                                                       
C     WHERE X(I)=ATANH (NEWR(I))   AND    Z(I)=ATANH(R(I))              
C     N1 AND N2 ARE THE DIMENSIONS OF THE TWO FIELDS, R(I) ARE THE      
C     UNCORRECTED CANONICAL CORELATIONS (ABSOLUTE VALUES) AND NL IS THE  
C     TIME SERIES LENGTH.                                               
C                                                                       
C     THE 95% ERROR INTERVAL FOR THE VARIABLES X(I) ARE GIVEN BY:       
C                                                                       
C     ERR(I)=+- 2/SQRT(NL)            
C
C     INPUT PARAMETERS:
C     * NSTEP           TIME SERIES LENGTH,
C     * NEOFX,NEOFY     DIMENSIONS OF VECTORS,
C     * CCA(N)          CANONICAL CORRELATIONS. 
C     OUTPUT PARAMETERS:
C     * CCA0(N)         UNBIASED CANONICAL CORRELATIONS,
C     * CCA1(N)         UPPER BOUND OF THE 95% CONFIDENCE INTERVAL,
C     * CCA2(N)         LOWER BOUND OF THE 95% CONFIDENCE INTERVAL.
C
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      REAL CCA(1),CCA0(1),CCA1(1),CCA2(1)
C
      NEOF=MAX(NEOFX,NEOFY)
      DO I=1,NEOF
         S=0.E0
         DO J=1,NEOF
            IF (J.NE.I) S=S+CCA(J)**2/(CCA(I)**2-CCA(J)**2)
         ENDDO
         S=2.E0*S*(1.E0-CCA(I)**2)
         S=S+CCA(I)**2+NEOFX+NEOFY-2
         S=S/(FLOAT(2*NSTEP)*CCA(I))
         CCA0(I)=0.5E0*LOG((1.E0+CCA(I))/(1.E0-CCA(I)))-S
         CCA1(I)=TANH(CCA0(I)+2.E0/SQRT(FLOAT(NSTEP)))
         CCA2(I)=TANH(CCA0(I)-2.E0/SQRT(FLOAT(NSTEP)))
         CCA0(I)=TANH(CCA0(I))
      ENDDO
      RETURN
      END
