#! /bin/sh
#
#    Last change: Mihai Dima 2001
#
if [ $# -lt 5 ]
then
  echo " "
  echo "            *****     MKEOF     *****      "
  echo " "
  echo "       The program calculates EOFs and PCs"
  echo " "
  echo "Syntax: mkeof <file1> <file2> <file3> <file4> <file5>"
  echo " "
  echo "  <file1> - input file"
  echo "  <file2> - output file - EOFs,PCs and eigenvalues"
  echo "  <file3> - output file - eigenvalues"
  echo "  <file4> - output file - EOFs"
  echo "  <file5> - output file - PCs"
  echo "  "
  exit 0
else
 echo "file1="$1
 echo "file2="$2
 echo "file3="$3
 echo "file4="$4
 echo "file5="$5
 cp $1 fort.1
fi
cat > mkeof.f << EOF
      PROGRAM PMKEOF
      IMPLICIT REAL (A-H,O-Z)

      PRINT *,' '      
      PRINT *,'            *****     MKEOF     *****'
      PRINT *,'   '
      NCHAN=1
      CALL PARLES (NCHAN,NDIM,NSET)
      NMAX=MAX0(NDIM,NSET)
      CALL MKEOF  (NDIM,NSET,NMAX)
      PRINT *,'REGULAR END!'     
      STOP 
      END
C******************************************************************
      SUBROUTINE PARLES (NCHAN,NDIM,NSET)
      IMPLICIT REAL (A-H,O-Z)
C
      INTEGER I1,I2,I3

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

      NSET=0
      NDIM=0  
    1 READ (NCHAN,END=100) I1,I2,I3,I4
      NDIM=MAX0(NDIM,I4)
      READ (NCHAN,END=101) 
      NSET=NSET+1
      GOTO1
      
  101 PRINT *,' MORE HEADER THEN DATA '     
  100 PRINT 1000,NCHAN,NDIM,NSET
      REWIND NCHAN
      RETURN    
   10 PRINT *,'INPUT FILE ON TAPE',NCHAN,' IS EMPTY'
      STOP 'MISTAKE'
 1000 FORMAT (1X,' TAPE=',I2,',  NDIM=',I5,',  NSET=',I5)
      END

C*****************************************************************
      SUBROUTINE MKEOF (IE,NTS,NMAX)  
      IMPLICIT REAL (A-H,O-Z)

      PARAMETER(NTIME=60,UNDEF=0.9E+10,UNDEFE=0.9*UNDEF)
      LOGICAL TRANS
      INTEGER IN1,OUT1,OUT2,OUT3,OUT4
C      COMMON / CHANNEL / IN1,OUT1,OUT2,OUT3,OUT4 
C      DATA IN1,OUT1,OUT2,OUT3,OUT4 /1,2,20,21,22/

      DATA  TRANS/.FALSE./      
C
      PARAMETER(IE77=10000,NTS77=15000)
      DIMENSION ID(4), IID(4) 

C      REAL RDEFX(IE), RDEFT(NTS)
C      REAL X(NMAX),CM(IE),XS(NMAX,NTIME)
C      POINTER (NPX,X),(NPCM,CM),(NPXS,XS),(NPRDEFX,RDEFX),
C     &        (NPRDEFT,RDEFT)
C      ALLOCATE (NPX,NMAX)
C      ALLOCATE (NPCM,IE)
C      ALLOCATE(NPXS,NMAX*NTIME)
C      ALLOCATE(NPRDEFX,IE)
C      ALLOCATE(NPRDEFT,NTS)

      REAL, DIMENSION(:,:), ALLOCATABLE :: XS
      REAL, DIMENSION(:), ALLOCATABLE :: X
      REAL, DIMENSION(:), ALLOCATABLE :: CM
      REAL, DIMENSION(:), ALLOCATABLE :: RDEFX
      REAL, DIMENSION(:), ALLOCATABLE :: RDEFT

      ALLOCATE (X(NMAX),CM(IE),XS(NMAX,NTIME))
      ALLOCATE (RDEFX(IE),RDEFT(NTS))

C AND ACTIVATE THE FOLLOWING LINES
C      PARAMETER(IE77=500,NTS77=3000)
C      DIMENSION ID(4), IID(4), RDEFX(IE77), RDEFT(NTS77)
C      REAL X(NTS77),XS(NTS77,NTIME),CM(IE77)
C NOW YOU HAVE DONE ALL CHANGES FOR STANDARD FORTRAN.
C
      IN1=1
      OUT1=2
      OUT2=20
      OUT3=21
      OUT4=22

C ------------------------------------------------------
C   Get from standard input the minimum percentage <MPER>
c ------------------------------------------------------
      PRINT*,'Minimum percentage of available data (0-100):'
      READ(5,*,END=120) MPER
      GOTO 121
  120   MPER = 70
  121 IF (MPER.GT.100) MPER = 100
      IF (MPER.LT.0)   MPER = 0
      PRINT*,'Minimum percentage of available data =',MPER
C
      IM = 0
      DO 100 I=1,IE
        RDEFX(I)=0.
        CM(I)=0.0
  100 CONTINUE
      DO 110 IT=1,NTS
         RDEFT(IT)=0.
  110 CONTINUE
C
C   check input file IN1
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)
          RDEFT(IM)=RDEFT(IM)+1.
          RDEFX(J)=RDEFX(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*MPER/100.
      IXX=0
      DO 13 J = 1,IE4
        IF (RDEFX(J).LE.ITEST2) THEN
         CM(J)=UNDEF
         RDEFX(J)=0.
        ELSE
         CM(J) = CM(J)/RDEFX(J)
         IXX=IXX+1
         RDEFX(J)=1 
        ENDIF
   13 CONTINUE
      ITEST1=IXX*MPER/100.
      ITT=0

      ITRANS=1
      DO 150 IT=1,IM
        IF (RDEFT(IT).GE.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 (RDEFT(IT).NE.IXX) ITRANS=0
C***
C
          ITT=ITT+1
          RDEFT(IT)=1.
        ELSE
          RDEFT(IT)=0.
        ENDIF
  150 CONTINUE
      IGAP=0
      DO 151 J=1,IE4
        IF (RDEFX(J).EQ.0) THEN
          IGAP=IGAP+1
        ELSE
          CM(J-IGAP)=CM(J)
        ENDIF
  151 CONTINUE

      PRINT *,' number of gridpoints with less than ',ITEST2,
     &        ' default values: ',IXX
      PRINT *,' number of timesteps with less than ',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
      DO 14 I = 1,IM
        IF (RDEFT(I).LT.0.5)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 (RDEFX(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   

      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. IT=',ITT,
     +        ' < IX=',IXX
        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 

      CALL EVEOFPC (LS,IE,NS,NTS,CM,RDEFX,RDEFT,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,RDEFX,RDEFT,XS,IM,
     +                    TRANS,NMAX,ITT,X)
      IMPLICIT REAL (A-H,O-Z)


      PARAMETER(NTIME=60,UNDEF=0.9E+10,UNDEFE=0.9*UNDEF)
      LOGICAL TRANS
      INTEGER IN1,OUT1,OUT2,OUT3,OUT4

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

      PARAMETER(IE77=10000,NTS77=15000)
      DIMENSION ID(4), IID(4), RDEFX(IE), RDEFT(NTS)

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

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 :: CV
      REAL, DIMENSION(:,:), ALLOCATABLE :: EV
      REAL, DIMENSION(:), ALLOCATABLE :: PC

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

C AND ACTIVATE THE FOLLOWING LINES
C      PARAMETER(IE77=500,NTS77=3000)
C      DIMENSION ID(4), IID(4), RDEFX(IE77), RDEFT(NTS77)
C      REAL X(NTS77),E(IE77),CV(IE77,IE77),EV(IE77,IE77),
C     &     XS(NTS77,NTIME),CM(IE77),PC(IE77)
C NOW YOU HAVE DONE ALL CHANGES FOR STANDARD FORTRAN.
C
      IN1=1
      OUT1=2
      OUT2=20
      OUT3=21
      OUT4=22

      DO 100 I=1,LS
        DO 100 J=1,LS
          CV(I,J)=0.0
  100 CONTINUE
C
C   estimation of the covarianzmatrix CV(I,J)
C
C      PRINT *,'LS=',LS
      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
      IFAIL=0
      CALL PPPCOM (LS,LS,CV,E,PC,EV)
C
C  eigenvalues
      REWIND OUT1
      REWIND OUT2
C      PRINT *,'OUT1=',OUT1,' OUT2=',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)
      ES = 0
      DO 50 J = 1,LS
  50    ES = ES + E(J)
      PRINT *, 'sum of all ', LS,' eigenvalues: ',ES
      PRINT *, 'the first 10 eigenvalues and explained variances'
      DO 51 J = LS,LS-9,-1
  51    PRINT '(I10,2G12.4)', LS+1-J, E(J), E(J)*100./ES
C
C  eigenvectors, possibly back transformation
C
      REWIND OUT3 
      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 (RDEFX(I).LT.0.5) 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
                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  

        DO 70 J = LS,1,-1
            IGAP = 0
            EOFM=0
            DO 71 I=1,IE
              IF (RDEFX(I).LT.0.5) THEN
                IGAP=IGAP+1
                PC(I)=UNDEF
              ELSE
                PC(I)=EV(I-IGAP,J)
                EOFM=EOFM+PC(I)
              ENDIF
   71       CONTINUE
            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 (RDEFT(I).LT.0.5) 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 (RDEFX(K).LT.0.5) 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)
      STOP
      END
C*****************************************************************
      SUBROUTINE PPPCOM (NM,N,RM,DVECT,EVECT,ZMATR)
      IMPLICIT REAL (A-H,O-Z)
C-------------------------------------------------------------------
C     COMPUTE EIGENVALUES AND EIGENVECTORS (IN DECREASING EIGENVALUE
C     ORDER) OF A REAL SYMMETRIC (POSITIVE DEFINITE?) MATRIX.
C     ROUTINES 'CTRED2' AND 'CIMTQL2' 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, 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)
        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)

      IMPLICIT NONE
      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.

      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)
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
EOF
lf95 mkeof.f -o mkeof.x
mkeof.x $1 $2 << M
70
M
cp fort.2 $2
cp fort.20 $3
cp fort.21 $4
cp fort.22 $5
rm fort*
rm mkeof.x mkeof.f
exit