#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 5 ]
then
  echo " "
  echo "            *****     CANCOR     *****      "
  echo " "
  echo "     The program performs canonical correlation"
  echo " "
  echo "Syntax: cancor <file1> <file2> <file3> <file4> <file5> <int>"
  echo " "
  echo "      <file1> - input file 1"
  echo "      <file2> - input file 2"
  echo "      <file3> - computed EOFs for field 1"
  echo "      <file4> - computed EOFs for field 2"
  echo "      <file5> - output file"
  echo "       <int>  - nr. of considered EOFs"
  echo "  "
  exit 0
else
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 echo "file4="$4
 echo "file5="$5
 echo "Nr. of considered EOFs="$6
 cp $1 fort.1
 cp $2 fort.2
fi
cat > cancor.f <<EOF
      PROGRAM Pcancor
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER (I-N)

      INTEGER FGRENZ,GGRENZ,NEOFF,NEOFG

      PRINT *,' '
      PRINT *,'     *****    CANCOR     *****'
      PRINT *,'   '

C--  READING THE NUMBER OF EOF FROM STANDARD IN --

C     PRINT *,' ENTER NUMBER OF REGARDED EOF FOR COMPUTATION OF THE '
C     PRINT *,' CANONICAL CORRELATION FOR THE FIRST FIELD (FIELD F). '
C     PRINT *,' DEFAULT IS THE NUMBER OF THE FIRST EOF WHICH EXPLAINS '
C     PRINT *,' ONLY 1% OF THE TOTAL VARIANCE (NEOFF=0). '
C     PRINT *,'   '
      READ (5,*,END=111) NEOFF
C     PRINT *,' ENTER NUMBER OF REGARDED EOF FOR COMPUTATION OF THE '
C     PRINT *,' CANONICAL CORRELATION FOR THE SECOND FIELD (FIELD G). '
C     PRINT *,' DEFAULT IS THE NUMBER OF THE FIRST EOF WHICH EXPLAINS '
C     PRINT *,' ONLY 1% OF THE TOTAL VARIANCE (NEOFG=0). '
C     PRINT *,'   '
      READ (5,*,END=112) NEOFG
      GOTO 113
  111 NEOFF=0
  112 NEOFG=0   

C-- COMPUTATION OF THE FIELD DIMENSIONS AND IF NECESSARY COMPUTATIN OF
C   EIGENVALUES, EOFS AND PRICIPAL COMPONENTS FOR FIELD F AND G. --

  113 NCF=3
      NCG=4

C-- FOR FIELD F --

      NCHANF=1
      CALL CCAPLES (NCHANF,NDIMF,NSETF,IEOF)
      CALL CCAPLES (NCF,NDIM,NSET,IEOF)
C-- the eof-file is empty (IEOF=0) ----
      IF (IEOF.EQ.0) THEN
        PRINT *,'  THE SUBROUTINE CCAMKEOF IS CALLED'
        NMAXF=MAX0(NDIMF,NSETF)
        CALL CCAMKEOF (NDIMF,NSETF,NMAXF,NCHANF) 
      ENDIF
C-- compute the number of the last eigenvalue with more than 1%of the
C-----------  total variance  ----- 
      CALL EOFNUMB (NCF,FGRENZ,NDIMF)
      IF (NEOFF.EQ.0) NEOFF=FGRENZ

C-- THE SAME FOR FIELD G --

      IEOF=0
      NCHANG=2
      CALL CCAPLES (NCHANG,NDIMG,NSETG,IEOF)
      CALL CCAPLES (NCG,NDIM,NSET,IEOF)
      IF (IEOF.EQ.0) THEN
        PRINT *,'  THE SUBROUTINE CCAMKEOF IS CALLED'
        NMAXG=MAX0(NDIMG,NSETG)
        CALL CCAMKEOF (NDIMG,NSETG,NMAXG,NCHANG)     
      ENDIF

      CALL EOFNUMB (NCG,GGRENZ,NDIMG)
      IF (NEOFG.EQ.0) NEOFG=GGRENZ
      MAXGR=MAX0(FGRENZ,GGRENZ)
      
      PRINT *,NEOFF,' EIGENVALUES OF THE FIRST FIELD (F) AND '
      PRINT *,NEOFG,' OF THE SECOND FIELD (G) WILL BE REGARDED. '

C-- COMPUTATION OF THE CANONICAL CORRELATION 

      IF (NSETF.NE.NSETG) THEN
        PRINT *,'*** ERROR *** LENGTH OF TIMESERIES NOT EQUAL:',
     +        NSETF,'.NE.',NSETG
        STOP 'ERROR'
      ENDIF 

      CALL CANCOR (NDIMF,NDIMG,NSETF,NEOFF,NEOFG)

      STOP 
      END
C******************************************************************
      SUBROUTINE CCAPLES (NCHAN,NDIM,NSET,IEOF)
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C
      INTEGER I1,I2,I3


      REWIND NCHAN
      READ (NCHAN,END=10) ITIME,IVAR,ILEV,ISIZE
      REWIND NCHAN

      IEOF=0
      IPC =0
      IEIG=0
      NSET=0
      NDIM=0
    1 READ (NCHAN,END=100) I1,I2,I3,I4
      NDIM=MAX0(NDIM,I4)	
      READ (NCHAN,END=101) 
      IF (I1.EQ.-1.AND.(I2.EQ.503.OR.I2.EQ.504)) THEN
        IEIG=IEIG+1
      ELSE IF (I2.EQ.503.AND.I1.GT.-1) THEN
        IEOF=IEOF+1
        NEOF=I4
      ELSE IF (I2.EQ.504.AND.I1.GT.-1) THEN
        IPC =IPC +1
        NPC=I4
      ENDIF
      NSET=NSET+1
      GOTO1
      
  101 PRINT *,' MORE HEADER THEN DATA '     
  100 PRINT 1000,NCHAN,NDIM,NSET
      IF (IEIG.NE.0 .OR. IPC.NE.0 .OR. IEOF.NE.0) THEN
        PRINT 2000,NCHAN,IEIG,IEOF,IPC
        NSET=NPC
        NDIM=NEOF
      ENDIF
      REWIND NCHAN
      RETURN    
   10 PRINT *,'INPUT FILE ON TAPE',NCHAN,' IS EMPTY'
      IF (NCHAN.EQ.3.OR.NCHAN.EQ.4)RETURN
      STOP 'MISTAKE'
 1000 FORMAT (1X,' TAPE=',I2,',  NDIM=',I5,',  NSET=',I5)
 2000 FORMAT (1X,' THE FILE ON TAPE ',I2,' CONTAINS:',/ 
     +       ,I8,' DATARECORDS WITH EIGENVALUES '/
     +       ,I8,' DATARECORDS WITH EOFS '/
     +       ,I8,' DATARECORDS WITH PRINCIPAL COMPONENTS ',/)
      END

C************************************************************************
      SUBROUTINE EOFNUMB (NCHAN,NEOF,NDIM)
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER (I-N)

      REAL EXVAR,TESTVAR,ONE,EVSUM,EXSUM
      REAL, DIMENSION(:), ALLOCATABLE :: EIGEN
      INTEGER NEOF,I1,I2,I3,I4

      ALLOCATE (EIGEN(NDIM))

C      IF (error .NE. 0) THEN
C     &    STOP 'Not enough storage for data; stopping...'
C      ENDIF
      
C      POINTER (NPEOF,EIGEN)
C      ALLOCATE (NPEOF,NDIM)

      DATA ONE /1./

      REWIND NCHAN
    1 READ (NCHAN,END=10) I1,I2,I3,I4
      READ (NCHAN) (EIGEN(I),I=1,I4)
      IF (I1.NE.-1) GOTO 1
      
      EXVAR=0
      DO 100 I=1,I4
        EXVAR=EXVAR+EIGEN(I)
  100 CONTINUE

      PRINT *, 'sum of all ', I4,' eigenvalues: ',EXVAR
      IENDE=10
      IF (I4.LT.IENDE) IENDE=I4
      PRINT *, 'the first ',IENDE,' eigenvalues and explained ',
     +         "variances"
      DO 51 J=1,IENDE
  51    PRINT '(I10,2G12.4)', J, EIGEN(J), EIGEN(J)*100./EXVAR
      EVSUM=0.0
      TESTVAR=EXVAR*0.01	
      DO 200 I=1,I4
        IF (EIGEN(I).LT.TESTVAR) THEN
          NEOF=I-1
          EXSUM=EVSUM*100./EXVAR
          PRINT 1000,NEOF,TESTVAR,EXVAR,NEOF,EXSUM
          RETURN
        ELSE
          EVSUM=EVSUM+EIGEN(I)
        ENDIF
  200 CONTINUE
 
      NEOF=I4
      PRINT *,' ALL EIGENVALUES HAVE MORE THEN 1% OF THE EXPLAINED ',
     +        'VARIANCE'
      RETURN

C+++ CAUSE AN ABORT ++

   10 PRINT *,' FILE ON TAPE ',NCHAN,' CONTAINS NO EIGENVALUES '
      PRINT *,' *** ERROR IN ROUTINE EOFNUMB ***'
      PRINT *,' *CANCOR*: FATAL END '
      CALL MINONE(ONE,XNEGO)
      ONE=SQRT(XNEGO)
      STOP
 1000 FORMAT (' THE ',I5,' VALUE IS THE LAST EIGENVALUE WITH MORE',
     +        ' THEN 1% (',E12.4,')',/,' OF THE EXPLAINED VARIANCE (',
     +        E12.4,'). THE',I5,' EIGENVALUES ',/'EXPLAIN ',F5.2,
     +       ' PERCENT OF THE TOTAL VARIANCE',/)
      END
C*****************************************************************
      SUBROUTINE MINONE(XIN,XOUT)
        XOUT=-XIN
      RETURN
      END
C*****************************************************************
      SUBROUTINE CCAMKEOF (IE,NTS,NMAX,IN1)  
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C
      PARAMETER(NTIME=60,UNDEF=0.9E+10,UNDEFE=0.9*UNDEF)
      LOGICAL TRANS
      INTEGER IN1,OUT1,OUT2,OUT3,OUT4,NIN1
C
C IF YOU USE THIS PROGRAM NOT ON A CRAY-2 DELETE THE
C FOLLOWING THREE LINES
      PARAMETER(IE77=10000,NTS77=15000)
      INTEGER ID(4), IID(4)
      COMMON / CHANNEL / NIN1,OUT1,OUT2,OUT3,OUT4 
C     INTEGER  IDEFX(IE), IDEFT(NTS)
C      REAL X(NMAX),CM(IE),XS(NMAX,NTIME)
C      POINTER (NPX,X),(NPCM,CM),(NPXS,XS),(NPIDEFX,IDEFX),
C     &        (NPIDEFT,IDEFT)
C      ALLOCATE (NPX,NMAX)
C      ALLOCATE (NPCM,IE)
C      ALLOCATE(NPXS,NMAX*NTIME)
C      ALLOCATE(NPIDEFX,IE)
C      ALLOCATE(NPIDEFT,NTS)

      REAL, DIMENSION(:), ALLOCATABLE :: X
      REAL, DIMENSION(:), ALLOCATABLE :: CM
      REAL, DIMENSION(:,:), ALLOCATABLE :: XS
      INTEGER, DIMENSION(:), ALLOCATABLE :: IDEFX
      INTEGER, DIMENSION(:), ALLOCATABLE :: IDEFT

      ALLOCATE (X(NMAX),CM(IE),XS(NMAX,NTIME),IDEFX(IE),IDEFT(NTS))

C AND ACTIVATE THE FOLLOWING LINES
C      PARAMETER(IE77=500,NTS77=3000)
C      INTEGER ID(4), IID(4), IDEFX(IE77), IDEFT(NTS77)
C      REAL X(NTS77),XS(NTS77,NTIME),CM(IE77)
C      COMMON / CHANNEL / NIN1,OUT1,OUT2,OUT3,OUT4 
C NOW YOU HAVE DONE ALL CHANGES FOR STANDARD FORTRAN.
C
C
      DO 80 IN=1,NMAX
        X(IN)=0.0
   80 CONTINUE
      DO 81 IN=1,NMAX
        DO 81 IS=1,NTIME 
          XS(IN,IS)=0.0
   81 CONTINUE
      DO 100 I=1,IE
        IDEFX(I)=0
        CM(I)=0.0
  100 CONTINUE
      DO 110 IT=1,NTS
         IDEFT(IT)=0
  110 CONTINUE
C
      IM = 0
      TRANS=.FALSE.
      NIN1=IN1
      OUT1=IN1+2
      OUT2=IN1+29
      OUT3=IN1+39
      OUT4=IN1+49
C
C   check input file IN1. comparism of the field size and timeseries length.
C   check of the number of missing values.
C
      REWIND IN1
      READ(IN1,END=400) IID
      IE4 = IID(4)
      
      PRINT *, 'first record:  ', IID
      REWIND IN1
   11 READ(IN1,END=200) ID
        READ(IN1) (X(J),J=1,ID(4))
        IF(ID(2).NE.IID(2)) GOTO 401
        IF(ID(3).NE.IID(3)) GOTO 402
        IF(ID(4).NE.IID(4)) GOTO 403
        IM = IM + 1
        DO 12 J = 1,ID(4)
          IF(X(J).GT.UNDEFE) GOTO 12
          CM(J) = CM(J) + X(J)
          IDEFT(IM)=IDEFT(IM)+1
          IDEFX(J)=IDEFX(J)+1
   12   CONTINUE
        GOTO 11
C
  200 PRINT *, 'last record:   ', ID(1)
      PRINT *, 'number of records found: ',IM
      PRINT *, 'vector length: ', IE4
      ITEST2=IM-(IM*30/100)
      IXX=0
      DO 13 J = 1,IE4
        IF (IDEFX(J).LT.ITEST2) THEN
         CM(J)=UNDEF
         IDEFX(J)=0
        ELSE
         CM(J) = CM(J)/IDEFX(J)
         IXX=IXX+1
         IDEFX(J)=1 
        ENDIF
   13 CONTINUE
      ITEST1=IXX-(IXX*30/100)
      ITT=0

      ITRANS=1
      DO 150 IT=1,IM
        IF (IDEFT(IT).GT.ITEST1) THEN
C
C*** THIS LINE WAS INCLUDED BECAUSE IT_S NOT CLEAR IF ISOLATED DEFAULT
C*** VALUES IN THE DATA SET CAUSE TROUBLE AT THE ESTIMATION OF EOFS IF THE
C*** TRANSFORMATION IS PERFORMED. 
C
C          IF (IDEFT(IT).NE.IXX) ITRANS=0
C***
C
          ITT=ITT+1
          IDEFT(IT)=1
        ELSE
          IDEFT(IT)=0
        ENDIF
  150 CONTINUE
      IGAP=0
      DO 151 J=1,IE4
        IF (IDEFX(J).EQ.0) THEN
          IGAP=IGAP+1
        ELSE
          CM(J-IGAP)=CM(J)
        ENDIF
  151 CONTINUE

      PRINT *,' number of gridpoints with less than ',IM-ITEST2,
     &        ' default values: ',IXX
      PRINT *,' number of timesteps with less than ',IXX-ITEST1,
     &        ' default values: ',ITT
C
      IF(IXX.GT.ITT) TRANS=.TRUE.
      IF(ITRANS.EQ.0) TRANS=.FALSE.
C
C  if number of samples less than spatial degrees of freedom,
C  solve transposed problem. To do so, an auxilary file is generated
C  in EXTRA format with transposed data matrix, on tape10
C 
      IMIT=1
      REWIND IN1
      REWIND 10
      REWIND 11
      DO 14 I = 1,IM
        IF (IDEFT(I).EQ.0)THEN
          READ(IN1)
          READ(IN1)
        ELSE
          READ(IN1) ID
          READ(IN1) (X(J),J=1,ID(4))
          IGAP=0
          DO 202 J=1,ID(4)
            IF (IDEFX(J).EQ.0)THEN
             IGAP=IGAP+1
            ELSE
             X(J-IGAP)=X(J)
            ENDIF
  202     CONTINUE
          WRITE(11) ID
          WRITE(11) (X(J),J=1,IXX)
          IF(.NOT.TRANS) THEN
           DO 201 J=1,IXX
             IF(X(J).LT.UNDEFE) X(J)=X(J)-CM(J)
  201      CONTINUE
           WRITE(10) ID
           WRITE(10) (X(J),J=1,IXX)
          ENDIF
        ENDIF
   14 CONTINUE   

   15 IF(.NOT.TRANS) THEN

        IF (ITRANS.EQ.0) THEN
         PRINT *,'Default values in data-matrix-->no transformation'
        ELSE
         PRINT *,'Data-matrix is NOT transposed: IX=',IXX,
     &           ' < IT=',ITT
        ENDIF
        LS = IXX
        NS = ITT
        IF(LS.GT.IE77) GOTO 404

      ELSE

C
        LS = ITT
        NS = IXX
        IF(LS.GT.IE77 .OR. NS.GT.NTS77) GOTO 404
        PRINT *, 'data matrix is transposed'
        REWIND 10
        DO 20 K = 0,IXX-1,NTIME
          REWIND 11
          KE = MIN0(NTIME,IXX-K)
          IF(K.EQ.0) THEN
            DO 21 I = 1,ITT
              READ(11) 
   21         READ(11) (XS(I,J),J=1,KE)
          ELSE
            DO 22 I = 1,ITT
              READ(11) 
   22         READ(11) (X(J),J=1,K), (XS(I,J),J=1,KE)
          ENDIF
          DO 23 J = 1,KE
            DO 300 I=1,ITT
              IF (XS(I,J).LT.UNDEFE) XS(I,J)=XS(I,J)-CM(K+J)
  300       CONTINUE
            WRITE(10) K+J, IID(2), IID(3), IM
   23       WRITE(10) (XS(I,J),I=1,IM)
   20   CONTINUE

      ENDIF 

C-- Computation of the eigenvalues, eofs and pcs

      CALL EVEOFPC (LS,IE,NS,NTS,CM,IDEFX,IDEFT,XS,IM,TRANS,NMAX,ITT,X)

      RETURN
C
C  fata ends
  400 PRINT *, 'the input file is empty '
      STOP 'MKEOF'
  401 PRINT *, 'variable code changes IVAR1(',IID(2),').NE.IVAR(',
     +          ID(2),')'
      STOP 'MKEOF'
  402 PRINT *, 'level code changes ILEV1(',IID(3),').NE.ILEV(',
     +          ID(3),')'
      STOP 'MKEOF'
  403 PRINT *, 'variable field length ILEN1(',IID(4),').NE.ILEN(',
     +          ID(4),')'
      STOP 'MKEOF'
  404 PRINT *, 'too many samples: ', NS,'>',NTS77
      PRINT *, 'or: vector is too long: ', LS,'>',IE77
      STOP 'MKEOF'
      END
C*********************************************************************
      SUBROUTINE EVEOFPC (LS,IE,NS,NTS,CM,IDEFX,IDEFT,XS,IM,
     +                    TRANS,NMAX,ITT,X)
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER (I-N)


      PARAMETER(NTIME=60,UNDEF=0.9E+10,UNDEFE=0.9*UNDEF)
      LOGICAL TRANS
      INTEGER IN1,OUT1,OUT2,OUT3,OUT4
C
C IF YOU USE THIS PROGRAM NOT ON A CRAY-2 DELETE THE
C FOLLOWING FOUR LINES
      PARAMETER(IE77=10000,NTS77=15000)
      INTEGER ID(4), IID(4), IDEFX(IE), IDEFT(NTS)

      REAL X(NMAX),CM(IE),XS(NMAX,NTIME)

      COMMON / CHANNEL / IN1,OUT1,OUT2,OUT3,OUT4

C      REAL X(NMAX),E(IE),CV(LS,LS),EV(LS,LS),CM(IE),PC(IE),
C     &     XS(NMAX,NTIME) 
C      POINTER (NPE,E),(NPCV,CV),(NPEV,EV),(NPPC,PC)
C      ALLOCATE (NPE,IE)
C      ALLOCATE (NPEV,LS*LS)
C      ALLOCATE (NPCV,LS*LS)
C      ALLOCATE (NPPC,IE)

      REAL, DIMENSION(:), ALLOCATABLE :: E
      REAL, DIMENSION(:), ALLOCATABLE :: PC
      REAL, DIMENSION(:,:), ALLOCATABLE :: CV
      REAL, DIMENSION(:,:), ALLOCATABLE :: EV

      ALLOCATE (E(IE),CV(LS,LS),EV(LS,LS),PC(IE))


C AND ACTIVATE THE FOLLOWING LINES
C      PARAMETER(IE77=500,NTS77=3000)
C      INTEGER ID(4), IID(4), IDEFX(IE77), IDEFT(NTS77)
C      REAL X(NTS77),E(IE77),CV(IE77,IE77),EV(IE77,IE77),
C     &     XS(NTS77,NTIME),CM(IE77),PC(IE77)
C      COMMON / CHANNEL / IN1,OUT1,OUT2,OUT3,OUT4 
C NOW YOU HAVE DONE ALL CHANGES FOR STANDARD FORTRAN.
C

      DO 88 I=1,IE
        E(I)=0.0
        PC(I)=0.0
   88 CONTINUE

      DO 100 I=1,LS
        DO 100 J=1,LS
          CV(I,J)=0.0
          EV(I,J)=0.0
  100 CONTINUE
C
C   estimation of the covarianzmatrix CV(I,J)
C
      REWIND 10
      DO 31 I = 1,NS
        READ(10) ID
        READ(10) (X(J),J=1,LS)
        DO 31 J=1,LS
          DO 31 K=1,LS
            IF (X(J).LT.UNDEFE .AND. X(K).LT.UNDEFE) THEN
             CV(J,K) = CV(J,K) + X(J)*X(K)
            ENDIF
   31   CONTINUE
C
      DO 36 J=1,LS
        DO 36 K = 1,LS
  36      CV(J,K) = CV(J,K)/ITT

Cc
Cc  derivation of the eigenvalues and -vectors with a nag-lib routine
      IFAIL=0
C      CALL F02ABF(CV,IE,LS,E,EV,IE,PC,IFAIL)
C
C on the cray-2 use the optimized cray eispack routines in PPPCOM
C
      CALL PPPCOM (LS,LS,CV,E,PC,EV)
C
C  eigenvalues
      REWIND OUT1
      REWIND OUT2
      WRITE(OUT2) -1, 503, 9999, LS
      WRITE(OUT2) (E(J),J=LS,1,-1)
      WRITE(OUT1) -1, 503, 9999, LS
      WRITE(OUT1) (E(J),J=LS,1,-1)
C
C  eigenvectors, possibly back transformation
C
      REWIND OUT3 
C-- with transformation ----
      IF(TRANS) THEN
        DO 73 K = 0,LS-1,NTIME
          REWIND 10
          KE = MIN0(NTIME,LS-K)
          DO 75 I = 1,NS
            READ(10) ID
            READ(10) (X(L),L=1,LS)
            DO 74 J = 1,KE
   74         XS(I,J) = 0.
            DO 75 J = 1,KE
              DO 76 L = 1,LS
                IF (X(L).LT.UNDEFE) THEN 
                 XS(I,J) = XS(I,J) + X(L)*EV(L,LS+1-K-J)              
                ENDIF
   76         CONTINUE
   75       CONTINUE
          DO 69 J = 1,KE
            XN = 0.
            DO 68 I = 1,NS
              XN = XN + XS(I,J)**2
   68       CONTINUE
            XN=1./SQRT(XN)
            IGAP=0
            EOFM=0
            DO 67 I=1,IE
              IF (IDEFX(I).EQ.0) THEN
                IGAP=IGAP+1
                PC(I)=UNDEF
              ELSE
                PC(I)=XS(I-IGAP,J)*XN
                EOFM=PC(I)+EOFM
              ENDIF
   67       CONTINUE
            IF (EOFM.LT.0) THEN 
              DO 66 I=1,IE
C--- always positiv mean value of all eofs ----
                IF (PC(I).LT.UNDEFE) PC(I)=PC(I)*(-1.)
   66         CONTINUE
            ENDIF
            WRITE(OUT3) K+J, 503, 9999, IE
            WRITE(OUT3) (PC(I),I=1,IE)
            WRITE(OUT1) K+J, 503, 9999, IE
            WRITE(OUT1) (PC(I),I=1,IE)
   69     CONTINUE
   73   CONTINUE

      ELSE  

C-- NO back transformation ---

        DO 70 J = LS,1,-1
            IGAP = 0
            EOFM=0
            DO 71 I=1,IE
              IF (IDEFX(I).EQ.0) THEN
                IGAP=IGAP+1
                PC(I)=UNDEF
              ELSE
                PC(I)=EV(I-IGAP,J)
                EOFM=EOFM+PC(I)
              ENDIF
   71       CONTINUE
C--- always positiv mean value of all eofs ----
            IF (EOFM.LT.0) THEN 
              DO 72 I=1,IE
                IF (PC(I).LT.UNDEFE) PC(I)=PC(I)*(-1.)
   72         CONTINUE
            ENDIF
          WRITE(OUT3) LS+1-J, 503, 9999, IE
          WRITE(OUT3) (PC(K),K=1,IE)
          WRITE(OUT1) LS+1-J, 503, 9999, IE
          WRITE(OUT1) (PC(K),K=1,IE)
   70 CONTINUE
      ENDIF
C
C   the principal components
C
      REWIND IN1
      REWIND 10

      DO 40 I = 1,IM
       READ(IN1) ID
       READ(IN1) (X(JJ),JJ=1,ID(4))
       REWIND OUT3
       J=0

       DO 41 JK=0,LS-1,NTIME
         NTIMAX=MIN0(NTIME,LS-JK)
         DO 45 KK=1,NTIMAX
           READ (OUT3) IID
           READ (OUT3) (XS(IA,KK),IA=1,IID(4))             
   45    CONTINUE

       IF (IDEFT(I).EQ.0) THEN
        DO 700 JJ=1,NTIMAX
          J=J+1
          PC(J)=UNDEF
  700   CONTINUE

       ELSE

        DO 43 JJ = 1,NTIMAX
          J=J+1
          PC(J) = 0.
          IGAP=0
          DO 42 K = 1,ID(4)
            IF (IDEFX(K).EQ.0) THEN
              IGAP=IGAP+1
            ELSE
              IF (X(K).LT.UNDEFE .AND. XS(K,JJ).LT.UNDEFE .AND. 
     +            CM(K-IGAP).LT.UNDEFE) THEN
                PC(J) = PC(J) +(X(K)-CM(K-IGAP))*XS(K,JJ)
              ENDIF 
            ENDIF
   42     CONTINUE
   43   CONTINUE

       ENDIF

   41  CONTINUE

        WRITE(10) ID(1), 504, 9999, LS
        WRITE(10) (PC(JJ),JJ=1,LS)
   40 CONTINUE
      CALL PCTP (LS,IM,OUT1,OUT4,X,XS,NTIME,NMAX)
      RETURN
C
C  fata ends
  400 PRINT *, 'input file 1 is empty '
      STOP 'MKEOF'
  401 PRINT *, 'variable code changes IVAR1(',IID(2),').NE.IVAR(',
     +          ID(2),')'
      STOP 'MKEOF'
  402 PRINT *, 'level code changes ILEV1(',IID(3),').NE.ILEV(',
     +          ID(3),')'
      STOP 'MKEOF'
  403 PRINT *, 'variable field length ILEN1(',IID(4),').NE.ILEN(',
     +          ID(4),')'
      STOP 'MKEOF'
  404 PRINT *, 'too many samples: ', NS,'>',NTS77
      PRINT *, 'or: vector is too long: ', LS,'>',IE77
      STOP 'MKEOF'
      END
C*****************************************************************
      SUBROUTINE PPPCOM (NM,N,RM,DVECT,EVECT,ZMATR)
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C-------------------------------------------------------------------
C     COMPUTE EIGENVALUES AND EIGENVECTORS (IN DECREASING EIGENVALUE
C     ORDER) OF A REAL SYMMETRIC (POSITIVE DEFINITE?) MATRIX.
C     ROUTINES "TRED2" AND "IMTQL2" OF THE EISPACK-PACKAGE ARE CALLED.
C
C --- PARAMETERS
C
C --- NM....... FIRST DIMENSION OF MATRICES
C --- N........ DIMENSION OF PROBLEM
C --- RM....... GIVEN MATRIX
C --- DVECT(J). EIGENVALUE NO. J
C --- EVECT.... AUXILIARY ARRAY
C --- ZMATR(I,J) CONTAINS EIGENVECTOR NO. J (I=1,...,N)
C-------------------------------------------------------------------
      DIMENSION RM(NM,N),DVECT(N),EVECT(N),ZMATR(NM,N)
      CALL TRED2(NM,N,RM,DVECT,EVECT,ZMATR)
      CALL IMTQL2(NM,N,DVECT,EVECT,ZMATR,IERR)
C
      IF(IERR.NE.0) PRINT 100, IERR
C
C --- CHECK EIGENVALUES AND EIGENVECTORS
C
      DO 3 K=1,N
      DV = DVECT(K)
      EMAX = 0.
      DO 2 I=1,N
      XZ=0.0
      DO 1 J=1,N
    1 XZ   = XZ+RM(I,J)*ZMATR(J,K)
    2 EMAX = AMAX1(EMAX , ABS(XZ - DV*ZMATR(I,K)))
    3 IF(EMAX.GE.(DV*1.E-8).AND.DV.GT.1.E-4) PRINT 200,N+1-K,EMAX,DV
C
  100 FORMAT(' PPPCOM, ERROR TYPE...',I4)
  200 FORMAT(' PPPCOM, NO. OF EIGENVALUE',I4,', MAX. ERROR'
     &        ,E10.2,' EIGENVALUE ',E14.6)
      RETURN
      END
C*****************************************************************
	SUBROUTINE TRED2 (NM,N,A,D,E,Z)
        IMPLICIT REAL (A-H,O-Z)
        IMPLICIT INTEGER (I-N)
        INTEGER I,J,K,L,N,II,NM,JP1
        REAL A(NM,N),D(N),E(N),Z(NM,N)
        REAL F,G,H,HH,SCALE
        REAL SQRT,ABS,SIGN

	DO 100 I=1,N
	  DO 100 J=1,I
	    Z(I,J)=A(I,J)
  100   CONTINUE

        IF (N.EQ.1) GO TO 320
C********** FOR I=N STEP -1 UNTIL 2 DO -- *****************
	DO 300 II=2,N
	  I=N+2-II
	  L=I-1
	  H=0.0
	  SCALE=0.0
	  IF(L.LT.2)GO TO 130
C******* SCALE ROW (ALGOL TOL THEN NOT NEEDED) ************
      	  DO 120 K=1,L
  120	    SCALE=SCALE+ABS(Z(I,K))

	  IF (SCALE.NE.0.0) GOTO 140
  130	  E(I)=Z(I,L)
       	  GOTO 290

  140	  DO 150 K=1,L
	    Z(I,K)=Z(I,K)/SCALE
	    H=H+Z(I,K)*Z(I,K)
  150	  CONTINUE

	  F=Z(I,L)
	  G=-SIGN(SQRT(H),F)
	  E(I)=SCALE*G
      	  H=H-F*G
	  Z(I,L)=F-G
          F=0.0

	DO 240 J=1,L
	  Z(J,I)=Z(I,J)/H
	  G=0.0
C*********  FORM ELEMENT OF A*U *******************
	  DO 180 K=1,J
  180	  G=G+Z(J,K)*Z(I,K)

	  JP1=J+1
	  IF(L.LT.JP1)GOTO 220

	  DO 200 K=JP1,L
  200	  G=G+Z(K,J)*Z(I,K)
C ********  FORM ELEMENT OF P **********************
  220  	  E(J)=G/H
	  F=F+E(J)*Z(I,J)
  240   CONTINUE

	HH=F/(H+H)
C ******** FORM REDUCED A **************************
	DO 260 J=1,L
	  F=Z(I,J)
	  G=E(J) - HH*F
	  E(J)=G

	  DO 260 K=1,J
	    Z(J,K)=Z(J,K)-F*E(K)-G*Z(I,K)
  260	CONTINUE

  290	D(I)=H
  300 CONTINUE

  320 D(1)=0.0
      E(1)=0.0
C ********* ACCUMULATION OF TRANSFORMATION MATRICES ******
      DO 500 I=1,N
	L=I-1
	IF (D(I).EQ.0.0) GOTO 380
	  
	DO 360 J=1,L
	  G=0.0

	  DO 340 K=1,L
  340	    G=G+Z(I,K)*Z(K,J)

	  DO 360 K=1,L
 	    Z(K,J)=Z(K,J)-G*Z(K,I)
  360	CONTINUE

  380	D(I)=Z(I,I)
	Z(I,I)=1.0
	IF (L.LT.1) GOTO 500

	DO 400 J=1,L
	  Z(I,J)=0.0
	  Z(J,I)=0.0
  400	CONTINUE

  500 CONTINUE

      RETURN
      END
C*****************************************************************
      SUBROUTINE IMTQL2 (NM,N,D,E,Z,IERR)

      INTEGER I,J,K,L,M,N,II,NM,MML,IERR
      REAL D(N),E(N),Z(NM,N)
      REAL B,C,F,G,P,R,S,MACHEP
      REAL SQRT,ABS,SIGN

C ****** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING
C        THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.

C      MACHEP=PRECIS
C OLD MACHEP      MACHEP=0.00000000000001
      MACHEP=0.00000000000001
      IERR=0
      IF (N.EQ.1) GOTO 1001

      DO 100 I=2,N
  100 E(I-1)=E(I)

      E(N)=0.0
      DO 240 L=1,N
        J=0
C **** LOOK FOR SMALL SUB-DIAGONAL ELEMENT ******************
  105   DO 110 M=L,N
	  IF(M.EQ.N) GOTO 120
          IF(ABS(E(M)).LE.MACHEP*(ABS(D(M))+ABS(D(M+1)))) GOTO 120
  110   CONTINUE

  120   P=D(L)
        IF(M.EQ.L) GOTO 240
	IF(J.EQ.30) GOTO 1000
	J=J+1
C ******* FORM SHIFT *****************************************
	G=(D(L+1)-P)/(2.0+E(L))
	R=SQRT(G*G+1.0)
	G=D(M)-P+E(L)/(G+SIGN(R,G))
	S=1.0
	C=1.0
	P=0.0
	MML=M-L
C ********** FOR I=M-1 STEP -1 UNTIL L DO -- *******************
  	DO 200 II=1,MML
	  I=M-II
	  F=S*E(I)
	  B=C*E(I)
	  IF(ABS(F).LT.ABS(G))GOTO 150
	  C=G/F
	  R=SQRT(C*C+1.0)
	  E(I+1)=F*R
	  S=1.0/R
	  C=C*S
	  GOTO 160
  150	  S=F/G
	  R=SQRT(S*S+1.0)
	  E(I+1)=G*R
	  C=1.0/R
	  S=S*C
  160	  G=D(I+1)-P
	  R=(D(I)-G)*S+2.0*C*B
	  P=S*R
	  D(I+1)=G+P
	  G=C*R-B
C ********  FORM VECTOR  ****************************
	  DO 180 K=1,N
	    F=Z(K,I+1)
	    Z(K,I+1)=S*Z(K,I)+C*F
	    Z(K,I)=C*Z(K,I)-S*F
  180	  CONTINUE

  200   CONTINUE
	D(L)=D(L)-P
	E(L)=G
	E(M)=0.0
	GOTO 105
  240 CONTINUE

C ******** ORDER EIGENVALUES AND EIGENVECTORS **********
      DO 300 II=2,N
	I=II-1
	K=I
	P=D(I)

	DO 260 J=II,N
	  IF(D(J).GE.P) GOTO 260
	  K=J
	  P=D(J)
  260   CONTINUE

	IF (K.EQ.I) GOTO 300
	D(K)=D(I)
	D(I)=P

	DO 280 J=1,N
	  P=Z(J,I)
	  Z(J,I)=Z(J,K)
	  Z(J,K)=P
  280   CONTINUE
  300 CONTINUE
      GOTO 1001
C ****** SET ERROR NO CONVERGENCE AFTER 30 ITERATIONS *******
 1000 IERR=L
 1001 RETURN
      END
C*****************************************************************
      SUBROUTINE PCTP (IE,NTS,OUT1,OUT4,X,XS,NTIME,NMAX)
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C
C   TransPosition of data matrix given in EXTRA format 
C   data have to be homogenous
      PARAMETER(IE77=10000,NTS77=15000)
      REAL X(NMAX),XS(NMAX,NTIME)
      INTEGER ID(4),IID(4),OUT1,OUT4
C
      IM=NTS
      LS=IE
      REWIND OUT4

      DO 100 I=0,LS,NTIME
        II=I+NTIME
        IF (II.GT.LS) II=LS
        NTMIN=MIN0(NTIME,LS-I)
        REWIND 10
        DO 110 K = 1,IM
          READ(10) ID
          READ(10) (X(J),J=1,II)
          DO 111 L=1,NTMIN
            XS(K,L) = X(II-NTMIN+L)
  111     CONTINUE
  110   CONTINUE
        DO 112 L=1,NTMIN
          WRITE(OUT1) I+L,ID(2),ID(3),IM
          WRITE(OUT1) (XS(K,L),K=1,IM) 
          WRITE(OUT4) I+L,ID(2),ID(3),IM
          WRITE(OUT4) (XS(K,L),K=1,IM) 
  112   CONTINUE
  100 CONTINUE

      RETURN
      END
C#######################################################################
      SUBROUTINE CANCOR (NDF,NDG,NSER,NEOFF,NEOFG)
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER (I-N)

      PARAMETER (UNDEF=0.9E+10,UNDEFE=0.99*UNDEF)

      INTEGER IFAIL,IT,IXA,IXB,ICAN,IEOF,OUT5,OUT6
      REAL    FAKTOR,TESTW,EWERT,FGSUM(2)

C      REAL    TS1(NSER),TS2(NSER),CMAT(NEOFF,NEOFG),EVAL(NEOFF),
C     +        PMAT(NEOFF,NEOFF),FHILF(NDG),EVECA(NEOFF,NEOFF),
C     +        CTMAT(NEOFG,NEOFF),EVECB(NEOFG,NEOFF),
C     +        PC(NEOFG,NSER),TSC(NSER,NEOFG),EPSI(2,NEOFF),
C     +        EOFVAL(NEOFG),MAP(NDG,NEOFF),TFAK(NDG,NEOFF),
C     +        FAKTOR,TESTW,EWERT,FGSUM(2),FSUM(NDG),
C     +        VARI(NDG,NEOFG)
C      POINTER (NPTS1,TS1),(NPTS2,TS2),(NPCMAT,CMAT),(NPEVAL,EVAL),
C     +    (NPPMAT,PMAT),(NPFHILF,FHILF),(NPEVECA,EVECA),(NPCTMAT,CTMAT)
C     +    ,(NPEVECB,EVECB),(NPPC,PC),(NPTSC,TSC),(NPEPSI,EPSI),
C     +    (NPEOFVAL,EOFVAL),(NPMAP,MAP),(NPTFAK,TFAK),(NPFSUM,FSUM),
C     +    (NPVARI,VARI)
C      ALLOCATE (NPTS1,NSER)
C      ALLOCATE (NPTS2,NSER)
C      ALLOCATE (NPCMAT,NEOFF*NEOFG)
C      ALLOCATE (NPEVAL,NEOFF)
C      ALLOCATE (NPPMAT,NEOFF*NEOFF)
C      ALLOCATE (NPFHILF,NDG)
C      ALLOCATE (NPEVECA,NEOFF*NEOFF)
C      ALLOCATE (NPCTMAT,NEOFG*NEOFF)
C      ALLOCATE (NPEVECB,NEOFF*NEOFG)
C      ALLOCATE (NPPC,NEOFG*NSER)
C      ALLOCATE (NPTSC,NSER*NEOFG)
C      ALLOCATE (NPEPSI,2*NEOFF)
C      ALLOCATE (NPEOFVAL,NEOFG)
C      ALLOCATE (NPMAP,NDG*NEOFF)
C      ALLOCATE (NPTFAK,NDG*NEOFF)
C      ALLOCATE (NPFSUM,NDG)
C      ALLOCATE (NPVARI,NDG*NEOFG)

      REAL, DIMENSION(:), ALLOCATABLE :: TS1
      REAL, DIMENSION(:), ALLOCATABLE :: TS2
      REAL, DIMENSION(:), ALLOCATABLE :: EVAL
      REAL, DIMENSION(:), ALLOCATABLE :: FHILF
      REAL, DIMENSION(:), ALLOCATABLE :: EOFVAL
      REAL, DIMENSION(:), ALLOCATABLE :: FSUM

      REAL, DIMENSION(:,:), ALLOCATABLE :: CMAT
      REAL, DIMENSION(:,:), ALLOCATABLE :: PMAT
      REAL, DIMENSION(:,:), ALLOCATABLE :: EVECA
      REAL, DIMENSION(:,:), ALLOCATABLE :: CTMAT
      REAL, DIMENSION(:,:), ALLOCATABLE :: EVECB
      REAL, DIMENSION(:,:), ALLOCATABLE :: PC
      REAL, DIMENSION(:,:), ALLOCATABLE :: TSC
      REAL, DIMENSION(:,:), ALLOCATABLE :: EPSI
      REAL, DIMENSION(:,:), ALLOCATABLE :: MAP
      REAL, DIMENSION(:,:), ALLOCATABLE :: TFAK
      REAL, DIMENSION(:,:), ALLOCATABLE :: VARI

      ALLOCATE (TS1(NSER),TS2(NSER),EVAL(NEOFF),FHILF(NDG))
      ALLOCATE (EOFVAL(NEOFG),FSUM(NDG))
      ALLOCATE (CMAT(NEOFF,NEOFG),PMAT(NEOFF,NEOFF))
      ALLOCATE (EVECA(NEOFF,NEOFG),CTMAT(NEOFG,NEOFF))
      ALLOCATE (EVECB(NEOFG,NEOFG),PC(NEOFG,NSER))
      ALLOCATE (TSC(NSER,NEOFG),EPSI(2,NEOFF))
      ALLOCATE (MAP(NDG,NEOFF),TFAK(NDG,NEOFF))
      ALLOCATE (VARI(NDG,NEOFG))

C--- PARAMETERCONTROL 

      IF (NEOFG.LT.NEOFF) GOTO 777
      IF (NDG.LT.NDF) GOTO 888

C--- SETTING OF TAPENAMES AND VARIABELS

      KANF=1
      KANG=2
      IN1 =3
      IN2 =4
      OUT5=21
      OUT6=21
      KO1F=60
      KO2F=70
      KO3F=80
      KO2G=71
      KO3G=81

C---- COMPUTATION OF CMAT(I,J)=<PCF(IT,I),PCG(IT,J)> ------------------
C----------------------------------------------------------------------

      IANZ=0
      IEV=0

C---- READ PCF FIELDS --------

    1 REWIND (IN2)
    2 READ (IN2,END=3) I1,I2,I3,I4
        IF(I2.EQ.504 .AND. I1.LE.NEOFG .AND. I1.GT.-1 .AND. 
     +     IEV.NE.0) THEN 
            IANZ=IANZ+1
            FHILF(I1)=1/SQRT(FHILF(I1))
            READ (IN2) (TSC(I,I1),I=1,I4)
            DO 90 I=1,I4
              IF (TSC(I,I1).LT.UNDEFE) THEN
                TSC(I,I1)=TSC(I,I1)*FHILF(I1)
              ENDIF 
   90       CONTINUE
	ELSE IF (I1.EQ.-1 .AND. IEV.EQ.0.AND.(I2.EQ.503.OR.I2.EQ.504))
     +           THEN
            IEV=1
            READ (IN2) (FHILF(I),I=1,NEOFG)
            GOTO 1
 	ELSE
            READ (IN2)
	ENDIF
      GOTO 2

    3 IF (IEV.EQ.0) THEN
        PRINT *,'** FOUND NO EIGENVALUE ON TAPE ',IN2,' **'
        STOP 'ERROR'
      ENDIF
      IF (IANZ.NE.NEOFG) THEN
        PRINT *,'** FOUND NOT ENOUGH PRINCIPAL COMPONENTS ON TAPE ',
     +           IN2,' **'
        STOP 'ERROR'
      ENDIF
      PRINT *,' READED ',NEOFG,' PRINCIPAL COMPONENTS FROM UNIT ',IN2

C---- read pcg fields ---

      IANZ=0
      IEV=0
    4 REWIND (IN1)
    5 READ (IN1,END=6) I1,I2,I3,I4
        IF(I2.EQ.504 .AND. I1.LE.NEOFF .AND. I1.GT.-1 .AND. 
     +     IEV.NE.0) THEN 
            IANZ=IANZ+1
            EOFVAL(I1)=1/SQRT(EOFVAL(I1))
            READ (IN1) (PC(I1,I),I=1,I4)
            DO 95 IT=1,NSER
              IF (PC(I1,IT).LT.UNDEFE) THEN
                PC(I1,IT)=PC(I1,IT)*EOFVAL(I1)
                TS1(IT)=PC(I1,IT)
              ELSE
                TS1(IT)=UNDEF
              ENDIF 
   95       CONTINUE
            DO 96 L=1,NEOFG
              DO 97 IT=1,NSER
                TS2(IT)=TSC(IT,L) 
   97         CONTINUE
              CMAT(I1,L)=RCCV(TS1,TS2,NSER,0)
   96       CONTINUE
	ELSE IF (I1.EQ.-1 .AND. IEV.EQ.0.AND.(I2.EQ.503.OR.I2.EQ.504))
     +           THEN
            IEV=1
            READ (IN1) (EOFVAL(I),I=1,NEOFG)
            GOTO 4
 	ELSE
            READ (IN1) 
	ENDIF
      GOTO 5

    6 IF (IEV.EQ.0) THEN
        PRINT *,'** FOUND NO EIGENVALUE ON TAPE ',IN1,' **'
        STOP 'ERROR'
      ENDIF
      IF (IANZ.NE.NEOFF) THEN
        PRINT *,'** FOUND NOT ENOUGH PRINCIPAL COMPONENTS ON TAPE ',
     +           IN1,' **'
        STOP 'ERROR'
      ENDIF
      PRINT *,' READED ',NEOFF,' PRINCIPAL COMPONENTS FROM UNIT ',IN1

C-- calculate the transposed matrix cmat''

      DO 300 I=1,NEOFG
      DO 300 J=1,NEOFF
        IF(CMAT(I,J).GT.UNDEFE) GOTO 666
        CTMAT(I,J)=CMAT(J,I)
  300 CONTINUE

C---- CALCULATION OF PMAT=CMAT*CMAT'' --------------------------------------
C--------------------------------------------------------------------------
C  CRAY SUBROUTINE    ---- MATRIX MULTIPLICATION -----

      CALL MXM (CMAT,NEOFF,CTMAT,NEOFG,PMAT,NEOFF)

C---- BERECHNUNG VON LAMBDA**2 UND EVECA(IXA,ICAN) ------------------
C---- CALCULATE EIGENVALUES AND EIGENVECTORS --

      CALL PPPCOM (NEOFF,NEOFF,PMAT,EVAL,FHILF,EVECA)

C--------------------------------------

      NEOFF1=NEOFF+1
      DO 200 IXA=1,NEOFF
        DO 201 ICAN=1,NEOFF
          FHILF(ICAN)=EVECA(IXA,NEOFF1-ICAN)
  201   CONTINUE
        DO 202 ICAN=1,NEOFF
          EVECA(IXA,ICAN)=FHILF(ICAN)
  202   CONTINUE
  200 CONTINUE
      DO 209  ICAN=1,NEOFF
  209 FHILF(ICAN)=EVAL(ICAN)

      CALL TROS (FHILF,NEOFF)

C-- create an info-printout ---

      NEMAX =MIN0(NEOFF,10)
      ES = 0
      DO 50 J = 1,NEOFF
  50    ES = ES + FHILF(J)
      PRINT *, NEOFF,' correlations found '
      PRINT *, 'the first ',NEOFF,' canonical correlation coefficents'
      DO 51 J = 1,NEMAX 
  51    PRINT '(I10,2G12.4)', J, FHILF(J)

	PRINT *,'the correlations between the  canonical time series'
	DO 52 J=1,NEOFF
	SS=SQRT(FHILF(J))
 52	PRINT '(I10,2G12.4)', J, SS  


C--- write correlation on tapes ---

        KST=-1
        KEV=520
        IDEF=9999
	REWIND (KO1F)
	REWIND (OUT6)
        WRITE (85) KST,KEV,IDEF,NEOFF
        WRITE (OUT6) KST,KEV,IDEF,NEOFF
        WRITE (OUT6) (FHILF(ICAN),ICAN=1,NEOFF)
	REWIND (OUT5)
        WRITE (OUT5) KST,KEV,IDEF,NEOFF
        WRITE (OUT5) (FHILF(ICAN),ICAN=1,NEOFF)
        WRITE (KO1F) KST,KEV,IDEF,NEOFF
        WRITE (KO1F) (FHILF(ICAN),ICAN=1,NEOFF)

C---- CALCULATION OF THE U TIMESERIES -----------------------------

      DO 400 IT=1,NSER
      DO 400 ICAN=1,NEOFF
         SUMMEA=0.0
         DO 401 IX=1,NEOFF
           IF (PC(IX,IT).LT.UNDEFE .AND. EVECA(IX,ICAN).LT.UNDEFE)THEN
             SUMMEA=SUMMEA+PC(IX,IT)*EVECA(IX,ICAN)
           ENDIF
  401    CONTINUE
         IF (SUMMEA.NE.0.0)THEN
           TSC(IT,ICAN)=SUMMEA
         ELSE
           TSC(IT,ICAN)=UNDEF
         ENDIF
  400 CONTINUE

      OPEN (UNIT=KO3F,FORM='UNFORMATTED')
      REWIND KO3F
      DO 581 ICAN=1,NEOFF
        WRITE (KO3F) ICAN,522,IDEF,NSER
        WRITE (KO3F) (TSC(IT,ICAN),IT=1,NSER)
        WRITE (OUT5) ICAN,522,IDEF,NSER
        WRITE (OUT5) (TSC(IT,ICAN),IT=1,NSER)
  581 CONTINUE

C---- CALCULATION OF THE CANONICAL MAPS OF FIELD F -----------

      DO 511 IC=1,NEOFF
        EPSI(1,IC)=0.0 
        DO 510 IX=1,NDG
          MAP(IX,IC)=0.0
          TFAK(IX,IC)=0.0
  510   CONTINUE
  511 CONTINUE

      REWIND KANF
      DO 520 IT=1,NSER
        READ (KANF) I1,I2,I3,I4
        READ (KANF) (FHILF(IX),IX=1,I4)
        DO 521 IC=1,NEOFF
          IF (TSC(IT,IC).LT.UNDEFE) THEN
            DO 522 IX=1,I4
              IF (FHILF(IX).LT.UNDEFE) THEN
                TFAK(IX,IC)=TFAK(IX,IC)+1. 
                MAP(IX,IC)=MAP(IX,IC)+FHILF(IX)*TSC(IT,IC)
              ENDIF
  522       CONTINUE
          ENDIF
  521   CONTINUE
  520 CONTINUE

      DO 530 IC=1,NEOFF
        DO 531 IX=1,NDF
          VARI(IX,IC)=0.0
          IF (IC.EQ.1) FSUM(IX)=0.0
          IF(TFAK(IX,IC).GT.(NSER*0.3))THEN
            MAP(IX,IC)=MAP(IX,IC)/TFAK(IX,IC)
          ELSE
            MAP(IX,IC)=UNDEF
          ENDIF
  531   CONTINUE
  530 CONTINUE 

C EXPLAINED VARIANCE OF EVERY CANONICAL MODE (FOR FMAPS)

      REWIND KANF
      FGSUM(1)=0.0
      DO 525 IT=1,NSER
        READ (KANF) I1,I2,I3,I4
        READ (KANF) (FHILF(IX),IX=1,I4)
        DO 528 IX=1,I4
          IF (FHILF(IX).LT.UNDEFE) THEN
            FGSUM(1)=FGSUM(1)+(FHILF(IX) * FHILF(IX))
            FSUM(IX)=FSUM(IX)+(FHILF(IX) * FHILF(IX))
          ENDIF
  528   CONTINUE 
        DO 526 IC=1,NEOFF
          IF (TSC(IT,IC).LT.UNDEFE) THEN
            DO 527 IX=1,I4
              IF (FHILF(IX).LT.UNDEFE.AND.MAP(IX,IC).LT.UNDEFE) THEN
                EPSI(1,IC)=EPSI(1,IC)+(FHILF(IX)-(TSC(IT,IC)*
     +                     MAP(IX,IC)) )**2
                VARI(IX,IC)=VARI(IX,IC)+(FHILF(IX)-(TSC(IT,IC)*
     +                     MAP(IX,IC)) )**2
              ENDIF
  527       CONTINUE
          ENDIF
  526   CONTINUE
  525 CONTINUE

      DO 540 IC=1,NEOFF
        DO 541 IX=1,I4
          IF (FSUM(IX).GT.0.0.AND.VARI(IX,IC).GT.0.0) THEN
            VARI(IX,IC)=100.*(1.-(VARI(IX,IC)/FSUM(IX)))
          ELSE
            VARI(IX,IC)=UNDEF
          ENDIF
  541   CONTINUE
  540 CONTINUE  
      OPEN (UNIT=KO2F,FORM='UNFORMATTED')
      REWIND KO2F
      DO 580 ICAN=1,NEOFF
        WRITE (KO2F) ICAN,521,IDEF,NDF
        WRITE (KO2F) (MAP(IA,ICAN),IA=1,NDF)
        WRITE (OUT5) ICAN,521,IDEF,NDF
        WRITE (OUT5) (MAP(IA,ICAN),IA=1,NDF)
        WRITE (OUT5) ICAN,525,IDEF,NDF
        WRITE (OUT5) (VARI(IA,ICAN),IA=1,NDF)
  580 CONTINUE

C---- CALCULATION OF EVECB ---------------------------------------

C  CRAY SUBROUTINE    ---------------
      CALL MXM (CTMAT,NEOFG,EVECA,NEOFF,EVECB,NEOFF)

C---------------------------------------

      DO 310 ICAN=1,NEOFF
        FAKTOR=1./SQRT(EVAL(NEOFF+1-ICAN))
        DO 311 IX=1,NEOFG
           EVECB(IX,ICAN)=EVECB(IX,ICAN)*FAKTOR
  311   CONTINUE
  310 CONTINUE

C--- READING DATA ---------------

      IANZ=0
      IEV=0
    7 REWIND (IN2)
    8 READ (IN2,END=9) I1,I2,I3,I4
        IF(I2.EQ.504 .AND. I1.LE.NEOFG .AND. I1.GT.-1) THEN 
            IANZ=IANZ+1
            READ (IN2) (PC(I1,IT),IT=1,I4)
 	ELSE IF (I1.LT.0.AND.(I2.EQ.504.OR.I2.EQ.503)) THEN
            READ (IN2) (EOFVAL(II),II=1,NEOFG)
            IEV=1
            DO 599 II=1,NEOFG
              EOFVAL(II)=1./SQRT(EOFVAL(II))
  599       CONTINUE
        ELSE
            READ (IN2) 
	ENDIF
      GOTO 8

    9 IF (IANZ.NE.NEOFG) THEN
        PRINT *,'** FOUND NOT ENOUGH PRINCIPAL COMPONENTS ON TAPE ',
     +           IN2,' **'
        STOP 'ERROR'
      ELSE IF (IEV.NE.1) THEN
        PRINT *,'** FOUND NO EIGENVALUES ON TAPE ',IN2,' **'
        STOP 'ERROR'
      ENDIF

C---- CALCULATION OF THE V TIMESERIES -----------------------------


      DO 600 IT=1,NSER
        DO 602 IX=1,NEOFG
          IF (PC(IX,IT).LT.UNDEFE) THEN
            PC(IX,IT)=PC(IX,IT)*EOFVAL(IX)
          ENDIF 
  602   CONTINUE      
      DO 600 ICAN=1,NEOFF
         SUMMEA=0.0
         DO 601 IX=1,NEOFG
           IF (PC(IX,IT).LT.UNDEFE .AND. EVECB(IX,ICAN).LT.UNDEFE)THEN
             SUMMEA=SUMMEA+(PC(IX,IT)*EVECB(IX,ICAN))
           ENDIF
  601    CONTINUE
         IF (SUMMEA.NE.0.0)THEN
           TSC(IT,ICAN)=SUMMEA
         ELSE
           TSC(IT,ICAN)=UNDEF
         ENDIF
  600 CONTINUE


	
      OPEN (UNIT=KO3G,FORM='UNFORMATTED')
      REWIND KO3G
      DO 591 ICAN=1,NEOFF
        WRITE (OUT6) ICAN,524,IDEF,NSER
        WRITE (OUT6) (TSC(IT,ICAN),IT=1,NSER)
        WRITE (KO3G) ICAN,524,IDEF,NSER
        WRITE (KO3G) (TSC(IT,ICAN),IT=1,NSER)
  591 CONTINUE

C---- CALCULATION OF THE CANONICAL MAPS OF FIELD G -----------

      FAKTOR = 1./REAL(NSER)

      DO 610 IC=1,NEOFF
        EPSI(2,IC)=0.0 
        DO 611 IX=1,NDG
          TFAK(IX,IC)=0.0
          MAP(IX,IC)=0.0
  611   CONTINUE
  610 CONTINUE

      REWIND KANG
      DO 620 IT=1,NSER
        READ (KANG) I1,I2,I3,I4
        READ (KANG) (FHILF(IX),IX=1,I4)
        DO 621 IC=1,NEOFF
          IF (TSC(IT,IC).LT.UNDEFE) THEN
            DO 622 IX=1,I4
              IF (FHILF(IX).LT.UNDEFE) THEN
                TFAK(IX,IC)=TFAK(IX,IC)+1. 
                MAP(IX,IC)=MAP(IX,IC)+FHILF(IX)*TSC(IT,IC)
              ENDIF
  622       CONTINUE
          ENDIF
  621   CONTINUE
  620 CONTINUE

      DO 630 IC=1,NEOFF
        DO 631 IX=1,NDG
          VARI(IX,IC)=0.0
          IF (IC.EQ.1) FSUM(IX)=0.0
          IF(TFAK(IX,IC).GT.(NSER*0.3))THEN
            MAP(IX,IC)=MAP(IX,IC)/TFAK(IX,IC)
          ELSE
            MAP(IX,IC)=UNDEF
          ENDIF
  631   CONTINUE
  630 CONTINUE 

C EXPLAINED VARIANCE OF EVERY CANONICAL MODE (FOR GMAPS)

      REWIND KANG
      ISUM1=0
      ISUM2=0
      FGSUM(2)=0.0
      DO 625 IT=1,NSER
        READ (KANG) I1,I2,I3,I4
        READ (KANG) (FHILF(IX),IX=1,I4)
        DO 628 IX=1,I4
          IF (FHILF(IX).LT.UNDEFE) THEN
            FGSUM(2)=FGSUM(2)+(FHILF(IX) * FHILF(IX))
            FSUM(IX)=FSUM(IX)+(FHILF(IX) * FHILF(IX))
          ENDIF
  628   CONTINUE 
        DO 626 IC=1,NEOFF
          IF (TSC(IT,IC).LT.UNDEFE) THEN
            DO 627 IX=1,I4
              IF (FHILF(IX).LT.UNDEFE.AND.MAP(IX,IC).LT.UNDEFE) THEN
                EPSI(2,IC)=EPSI(2,IC)+(FHILF(IX)-(TSC(IT,IC)*
     +                     MAP(IX,IC)) )**2
                VARI(IX,IC)=VARI(IX,IC)+(FHILF(IX)-(TSC(IT,IC)*
     +                     MAP(IX,IC)) )**2
              ENDIF
  627       CONTINUE
          ENDIF
  626   CONTINUE
  625 CONTINUE

      DO 640 IC=1,NEOFF
        DO 641 IX=1,I4
          IF (FSUM(IX).GT.0.0.AND.VARI(IX,IC).GT.0.0) THEN
            VARI(IX,IC)=100.*(1.-(VARI(IX,IC)/FSUM(IX)))
          ELSE
            VARI(IX,IC)=UNDEF
          ENDIF
  641   CONTINUE
  640 CONTINUE  

      OPEN (UNIT=KO2G,FORM='UNFORMATTED')
      REWIND KO2G
      DO 590 ICAN=1,NEOFF
        WRITE (OUT6) ICAN,523,IDEF,NDG
        WRITE (OUT6) (MAP(IB,ICAN),IB=1,NDG)
        WRITE (KO2G) ICAN,523,IDEF,NDG
        WRITE (KO2G) (MAP(IB,ICAN),IB=1,NDG)
        WRITE (OUT6) ICAN,526,IDEF,NDG
        WRITE (OUT6) (VARI(IB,ICAN),IB=1,NDG)
  590 CONTINUE

      CLOSE (KO2F)
      CLOSE (KO3F)
      CLOSE (KO2G)
      CLOSE (KO3G)

C PRINTOUT OF THE EXPLAINED VARIANCES

      PRINT*,' '
      PRINT *,'EXPLAINED VARIANCES OF THE CANONICAL MODES'
      PRINT *,' MODE  FMAPS  GMAPS '
      DO 700 IC=1,NEOFF
       EPSI(1,IC)=1.-(EPSI(1,IC)/FGSUM(1))
       EPSI(2,IC)=1.-(EPSI(2,IC)/FGSUM(2))
       PRINT 1000,IC,EPSI(1,IC)*100.,EPSI(2,IC)*100.
  700 CONTINUE
      STOP
  666 PRINT *,'--- DEFAULT VALUE IN COVARIANCEMATRIX. ---' 
      PRINT *,' A PRINCIPAL COMPONENT CONTAINS ONLY DEFAULT VALUES, '
      PRINT *,' CORECT IT !'
      STOP 'ERROR'
  777 PRINT *,'--- NUMBER OF REGARDED EOFS IS NOT CORECT ! ----'
      PRINT *,' FOR FIELD F GREATER THAN FOR FIELD G. '
      PRINT *,'    EOFF=',NEOFF,' > EOFG=',NEOFG
      STOP 'ERROR'
  888 PRINT *,'--- THE FIELD SIZE IS NOT CORECT ! ----'
      PRINT *,'  FIELD F GREATER THAN  FIELD G. '
      PRINT *,'    DIMF=',NDF,' > DIMG=',NDG
      STOP 'ERROR'
 1000 FORMAT (1X,I4,3X,F5.2,2X,F5.2)
      END
C**********************************************************
      subroutine mxm(a,nra,b,nca,c,ncb)
      implicit real (a-h,o-z)
      dimension a(nra,nca),b(nca,ncb),c(nra,ncb)
      call zero(c,nra*ncb)
      do 30 i=1,nra
         do 20 j=1,ncb
            do 10 k=1,nca
               c(i,j) = c(i,j) + a(i,k) * b(k,j)
10          continue
20       continue
30    continue
      return
      end
C**********************************************************
      SUBROUTINE ZERO (C,IE)
      REAL C(IE)
      DO 100 I=1,IE
        C(I)=0.0
  100 CONTINUE
      RETURN
      END
C*****************************************************************
      FUNCTION RCCV (TS1,TS2,NT,ITAU)
C-----------------------------------------------------------------
C     CROSS COVARIANCE FUNCTION VALUE BETWEEN TIME SERIES TS1 AND
C     TS2 OF LENGTHS NT AT DISCRETE LAG ITAU.
C     THE LAG IS IN TS2 AND ALWAYS ZERO OR POSITIVE
C-----------------------------------------------------------------
      DIMENSION TS1(NT),TS2(NT)
      RCCV = 0.
      NANZ = 0
      IF (ITAU.LT.0) THEN
        NTS=1-ITAU
        NTE=NT
      ELSE
        NTS=1
        NTE=NT-ITAU
      ENDIF
      DO 1 IT=NTS,NTE
        IF (TS1(IT).LE.0.9E9 .AND. TS2(IT+ITAU).LE.0.9E9) THEN
         NANZ = NANZ + 1
         RCCV = RCCV + TS1(IT) * TS2(IT+ITAU)
        ENDIF
    1 CONTINUE
      IF (NANZ.EQ.0)THEN
        RCCV=0.9E+10
      ELSE
        RCCV = RCCV / REAL(NANZ)
      ENDIF
      RETURN
      END
C*******************************************************************
      SUBROUTINE TROS(VECT,N)
      IMPLICIT REAL (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C-------------------------------------------------------------------
C     BRING THE COMPONENTS OF VECTOR "VECT" INTO REVERSE ORDER
C-------------------------------------------------------------------
      DIMENSION VECT(N)
      N1 = N + 1
      NH = N/2
      DO 1 I=1,NH
      P = VECT(I)
      VECT(I) = VECT(N1-I)
      VECT(N1-I) = P
    1 CONTINUE
      RETURN
      END
EOF
lf95 cancor.f -o cancor.x
cancor.x $1 $2 << M
$6
M
cp fort.3 $3
cp fort.4 $4
cp fort.21 $5
rm fort*
rm cancor.x cancor.f
echo "REGULAR END!"
exit