program pca_vaganov
c
c
c     reduce vaganov dataset (60 records) to its dominant pcs...
c
      parameter (ipcmax=40)
      parameter (maxproxy = 500)
      parameter (maxstat = 500)
      parameter (mmax = 500)
      parameter (max0 = mmax+20)
      parameter (nlarge=8000)
      parameter (nmax1 = 700)
      parameter (nmax0=690)
      parameter (nmax2=250)
      parameter (ntrnmax = 80)
      parameter (pi=3.14159265359d0)
      parameter (half=0.5d0)
      parameter (outofbounds=-999.d0)
c
      real *8 aprox(nmax0,mmax),wgt(mmax),
     $        weight(mmax),weight0(max0)
      real *8 weightprx(maxproxy)
      real *8 aint(nmax0,max0)
      real *8 aver0(nmax0,max0)
      real *8 sd(max0),datave(max0)
      real *8 proxave(maxproxy),sdprox(maxproxy),
     $          sdprox2(maxproxy)
      real *8 mean(max0),annmean(max0)
      real *8 S(mmax),S0(mmax),partsum(mmax+1),sum(mmax+1)
      real *8 SS(ipcmax)
      integer index0(max0)
      integer istart(maxproxy),ikeep(maxproxy)
      integer igood(max0)
      real *8 eof(mmax,ipcmax),tpc(nmax0,ipcmax)
      real *8 rpc(ipcmax)
      complex *16 AA(nmax1,mmax),UU(nmax1,nmax1),VV(mmax,mmax)
      real *8 B0(ntrnmax),x0(ipcmax),
     $        work0(ipcmax)
c
      real *8 beta(maxproxy,ipcmax)
c
      real *8 zero,dt,one,small,adum
      real *8 sumtot
      character name*53,name1*40,name2*40
      character flnm*60(maxproxy)
      character *20 nome0
      character *20 nome(maxproxy)
      character *24 nome2(maxproxy)
      character *48 flnm2
c
      zero = 0.d0
      one = 1.d0
      small = 1e-8
c
c     dt    = time interval spacing in units of years
c     nscan = length of time series in units of "dt"
c     iabv  = number of time series in analysis (ie, "spatial" array)
c
c     define limits for entire joint proxy/instrumental dataset
c
c     "nlow" will define reference year #1
c
      IP=0
      IP2 = 0
      nlow = 1300
      nhigh = 1980
      iproxmax = nhigh
      ntot = nhigh-nlow+1
c
c     define limits for raw data and training interval
c
      itrainmin0 = 1902
      itrainmin = itrainmin0
      itrainmax0 = nhigh
      itrainmax = itrainmax0
      ntrain = itrainmax0-itrainmin0+1
      iyearmax = 1993
      irawmin = 1902
      iearlyinst = 1820
      irawmax = nhigh
      nraw = irawmax-irawmin+1
c      
      do i=1,maxproxy
         istart(i)=99999
      end do
c
      open (unit=1,
     $   file='noamer.inf',
     $   status='old')
c
555   format (f5.2,1a20)
      i=0
1010  i=i+1
         if (i.gt.maxproxy) goto 415
         read (1,555,end=415) weightprx(i),nome(i)
         if (weightprx(i).lt.zero) then
            nome0=nome(i)
            i=i-1
            goto 1010
         else
           if (weightprx(i).lt.small) then
               weightprx(i)=small
           else
              if (iunif.eq.1) weightprx(i)=one
           endif
         endif
         if (weightprx(i).ge.zero) then
c
         flnm(i) = 
     $      nome(i)
         write (6,*) flnm(i)
         endif
      goto 1010
415   continue
      nproxy = i-1
      close (unit=1)
c     
      do j=1,nproxy
         ii = j
         open (unit=9,file=flnm(j),status='old')
         iflag = 0
         do i=1,nmax0
            aprox(i,ii)=-99999999.0
         end do
         do i=1,nlarge
            read (9,*,end=88) ayear,adat
            iyr = int(ayear+0.5)
            if (iflag.eq.0) then
                istart(j)=iyr
                iflag = 1
            endif
            if ((iyr.ge.nlow).and.(iyr.le.nhigh)) then
               iyear = iyr-nlow+1
               aprox(iyear,ii)=adat
            endif
         end do
88       close (unit=9)
c
c        extrapolate last data point if series ends before end of
c        training interval
c
         iyear = nhigh-nlow+2
737      iyear = iyear-1
         if (aprox(iyear,ii).lt.-99999990.0) goto 737
         do iyr = iyear+1,nhigh-nlow+1
            aprox(iyr,ii)=aprox(iyear,ii)
         end do
c       
      end do
c
c
      write (6,*) 'read in: '
      write (6,*) nproxy, ' time series'
c
c
c
9999  continue
c
c
      do j=1,nproxy
c 
       proxave(j)=zero
       do i=1,ntrain
          iyear = itrainmin+i-nlow
          proxave(j)=proxave(j)+aprox(iyear,j)
       end do
       proxave(j)=proxave(j)/float(ntrain)
c
c      remove 1902-19xx mean from training data
c
       do i=nlow,nhigh
          iyear = i-nlow+1
          aprox(iyear,j)=aprox(iyear,j)-proxave(j)
       end do
c
c      standardize proxy data
c
       amean = 0.0
       sdprox(j)=0.0
       do i=1,ntrain
         iyear = itrainmin+i-nlow
         sdprox(j)=sdprox(j)+aprox(iyear,j)**2
       end do
       sdprox(j)=sqrt(sdprox(j)/float(ntrain))
       do i=nlow,nhigh
          iyear = i-nlow+1
          aprox(iyear,j)=aprox(iyear,j)/sdprox(j)
       end do
      end do
c
c     now restandardized by detrended normalized variance
c
      do j=1,nproxy
c 
       sum1 = zero
       sum2 = zero
       sum5 = zero
       sum6 = zero
       do i=1,ntrain
          iyear = itrainmin+i-nlow
          sum1 = sum1 + float(i)
          sum2 = sum2 + float(i**2)
          sum5 = sum5 + aprox(iyear,j)
          sum6 = sum6 + aprox(iyear,j)*float(i)
       end do
c
       slope = (ntrain*sum6-sum1*sum5)
     $       /(ntrain*sum2-sum1**2)
       ainter = (sum2*sum5-sum1*sum6)/
     $      (ntrain*sum2-sum1**2)
c
c
       amean = 0.0
       sdprox2(j)=0.0
       do i=1,ntrain
         iyear = itrainmin+i-nlow
         standprx=aprox(iyear,j)-slope*float(i)
     $      - ainter
         sdprox2(j)=sdprox2(j)+standprx**2
       end do
       sdprox2(j)=sqrt(sdprox2(j)/float(ntrain))
       do i=nlow,nhigh
          iyear = i-nlow+1
          aprox(iyear,j)=aprox(iyear,j)/sdprox2(j)
       end do
       sdprox(j)=sdprox(j)*sdprox2(j)
      end do
c
      write (6,*) 'year proxies must date back to:'
      read (5,*) iproxmin
c
c
c     select proxy records to be used
c
      if (iset.eq.0) then
        mkeep = 0
        do i=1,nproxy
         if (istart(i).le.iproxmin) then
            mkeep = mkeep + 1
            ikeep(mkeep)=iabv+i
         endif
        end do
        write (6,*) 'total of ',mkeep,
     $      ' training records back to ',iproxmin
      endif
c 
c
      ipc = 15
      im1 = iproxmin
      im2 = iproxmax
      N0 = iproxmax-iproxmin+1
      jj = 0
      do j=1,nproxy
         if (istart(j).le.iproxmin) then
           jj = jj+1
           i = 0
           do i2=im1,im2
              i = i+1
              iyear = i2-nlow+1
              adum = aprox(iyear,j)
              AA(i,jj)=dcmplx(adum,zero)
           end do
         end if
      end do
      M0 = jj
      NU0 = M0
      NV0 = M0
c
c
      call CSVD(AA,nmax1,mmax,N0,M0,IP0,NU0,NV0,S,UU,VV)
c
c     collect eigenvalues, fractional variances explained..
c
      sumtot = zero
      do i=1,M0
           partsum(i)=zero
           sumtot = sumtot + S(i)**2
           do j=i,M0
             partsum(i)=partsum(i)+S(j)**2
           end do
      end do
c
      open (unit=7,file='eigen.out',status='unknown')
      do i=1,M0
         eigen = S(i)**2/sumtot
         write (7,434) i,
     $     eigen,one-real(partsum(i+1)/sumtot)
      end do
434   format (i4,2f9.4)
c
c
      do i=1,ipc
        id1 = i/10
        id2 = i-id1*10
        name1 = 'eof'
     $     //char(48+id1)//char(48+id2)//'.out'
        open (unit=30+i,file=name1,status='unknown')
        do j=1,M0
           eof(j,i)=real(VV(j,i))
           write (30+i,*) j,real(VV(j,i))
        end do
        close (unit=30+i)
      end do
c
      do i=1,ipc
        id1 = i/10
        id2 = i-id1*10
        name2 = 'pc'
     $     //char(48+id1)//char(48+id2)//'.out'
        open (unit=30+i,file=name2,status='unknown')
        do j=1,N0
           jj = iproxmin+j-1
           write (30+i,*) jj,real(UU(j,i))
        end do
        close (unit=30+i)
      end do
c
c
      stop
      end
c
c
      SUBROUTINE CSVD (A,MMAX,NMAX,M,N,IP,NU,NV,S,U,V)
C
C   Singular value decomposition of an  M by N  complex matrix  A,
C   where M .GT. N .  The singular values are stored in the vector
C   S. The first NU columns of the M by M unitary matrix U and the
C   first NV columns of the N by N unitary matrix V  that minimize
C   Det(A-USV*) are also computed.
C
C
C         P.A. Businger and G.H. Golub, "Singular Value Decomposition
C         of a Complex Matrix," Communications of the ACM, vol. 12,
C         pp. 564-565, October 1969.
C
C   This algorithm is reprinted by permission, Association for
C   Computing Machinery; copyright 1969.
C
      COMPLEX *16 A(MMAX,NMAX),U(MMAX,MMAX),V(NMAX,NMAX),Q,R
      REAL *8 S(NMAX),B(100000),C(100000),T(100000),
     $   zero,one,ETA,TOL,
     $        EPS
      ETA = 1.2d-7
      TOL = 2.4d-32
      zero = 0.d0
      one = 1.d0
      NP=N+IP
      N1=N+1
C   Householder reduction
      C(1)=zero
      K=1
10    K1=K+1
C   Elimination of A(I,K) , I=K+1,...,M
      Z=zero
      DO 20 I=K,M
20      Z=Z+REAL(A(I,K))**2+dIMAG(A(I,K))**2
      B(K)=zero
      IF (Z .LE. TOL)  GO TO 70
      Z=SQRT(Z)
      B(K)=Z
      W=CdABS(A(K,K))
      Q=dcmplx(one,zero)
      IF (W .NE. zero)  Q=A(K,K)/W
      A(K,K)=Q*(Z+W)
      IF (K .EQ. NP)  GO TO 70
      DO 50 J=K1,NP
        Q=dcmplx(zero,zero)
        DO 30 I=K,M
30        Q=Q+dCONJG(A(I,K))*A(I,J)
        Q=Q/(Z*(Z+W))
        DO 40 I=K,M
40        A(I,J)=A(I,J)-Q*A(I,K)
50      CONTINUE
C   Phase transformation
      Q=-dCONJG(A(K,K))/CdABS(A(K,K))
      DO 60 J=K1,NP
60      A(K,J)=Q*A(K,J)
C   Elimination of A(K,J) , J=K+2,...,N
70    IF (K .EQ. N)  GO TO 140
      Z=zero
      DO 80 J=K1,N
80      Z=Z+REAL(A(K,J))**2+dIMAG(A(K,J))**2
      C(K1)=zero
      IF (Z .LE. TOL)  GO TO 130
      Z=SQRT(Z)
      C(K1)=Z
      W=CdABS(A(K,K1))
      Q=dcmplx(one,zero)
      IF (W .NE. zero)  Q=A(K,K1)/W
      A(K,K1)=Q*(Z+W)
      DO 110 I=K1,M
        Q=dcmplx(zero,zero)
        DO 90 J=K1,N
90        Q=Q+dCONJG(A(K,J))*A(I,J)
        Q=Q/(Z*(Z+W))
        DO 100 J=K1,N
100       A(I,J)=A(I,J)-Q*A(K,J)
110     CONTINUE
C   Phase transformation
      Q=-dCONJG(A(K,K1))/CdABS(A(K,K1))
      DO 120 I=K1,M
120     A(I,K1)=A(I,K1)*Q
130   K=K1
      GO TO 10
C   Tolerance for negligible elements
140   EPS=zero
      DO 150 K=1,N
        S(K)=B(K)
        T(K)=C(K)
150   EPS=DMAX1(EPS,S(K)+T(K))
      EPS=EPS*ETA
C   Initialization of U and V
      IF (NU .EQ. 0)  GO TO 180
      DO 170 J=1,NU
        DO 160 I=1,M
160       U(I,J)=dcmplx(zero,zero)
170     U(J,J)=dcmplx(one,zero)
180   IF (NV .EQ. 0)  GO TO 210
      DO 200 J=1,NV
        DO 190 I=1,N
190       V(I,J)=dcmplx(zero,zero)
200     V(J,J)=dcmplx(one,zero)
C   QR diagonalization
210   DO 380 KK=1,N
        K=N1-KK
C   Test for split
220     DO 230 LL=1,K
          L=K+1-LL
          IF (ABS(T(L)) .LE. EPS)  GO TO 290
          IF (ABS(S(L-1)) .LE. EPS)  GO TO 240
230       CONTINUE
C   Cancellation of B(L)
240     CS=zero
        SN=one
        L1=L-1
        DO 280 I=L,K
          F=SN*T(I)
          T(I)=CS*T(I)
          IF (ABS(F) .LE. EPS)  GO TO 290
          H=S(I)
          W=SQRT(F*F+H*H)
          S(I)=W
          CS=H/W
          SN=-F/W
          IF (NU .EQ. 0)  GO TO 260
          DO 250 J=1,N
            X=REAL(U(J,L1))
            Y=REAL(U(J,I))
            U(J,L1)=dCMPLX(X*CS+Y*SN,0.)
250         U(J,I)=dCMPLX(Y*CS-X*SN,0.)
260       IF (NP .EQ. N)  GO TO 280
          DO 270 J=N1,NP
            Q=A(L1,J)
            R=A(I,J)
            A(L1,J)=Q*CS+R*SN
270         A(I,J)=R*CS-Q*SN
280       CONTINUE
C   Test for convergence
290     W=S(K)
        IF (L .EQ. K)  GO TO 360
C   Origin shift
        X=S(L)
        Y=S(K-1)
        G=T(K-1)
        H=T(K)
        F=((Y-W)*(Y+W)+(G-H)*(G+H))/(2.d0*H*Y)
        G=SQRT(F*F+one)
        IF (F .LT. zero)  G=-G
        F=((X-W)*(X+W)+(Y/(F+G)-H)*H)/X
C   QR step
        CS=one
        SN=one
        L1=L+1
        DO 350 I=L1,K
          G=T(I)
          Y=S(I)
          H=SN*G
          G=CS*G
          W=SQRT(H*H+F*F)
          T(I-1)=W
          CS=F/W
          SN=H/W
          F=X*CS+G*SN
          G=G*CS-X*SN
          H=Y*SN
          Y=Y*CS
          IF (NV .EQ. 0)  GO TO 310
          DO 300 J=1,N
            X=REAL(V(J,I-1))
            W=REAL(V(J,I))
            V(J,I-1)=dCMPLX(X*CS+W*SN,0.)
300         V(J,I)=dCMPLX(W*CS-X*SN,0.)
310       W=SQRT(H*H+F*F)
          S(I-1)=W
          CS=F/W
          SN=H/W
          F=CS*G+SN*Y
          X=CS*Y-SN*G
          IF (NU .EQ. 0)  GO TO 330
          DO 320 J=1,N
            Y=REAL(U(J,I-1))
            W=REAL(U(J,I))
            U(J,I-1)=dCMPLX(Y*CS+W*SN,0.)
320         U(J,I)=dCMPLX(W*CS-Y*SN,0.)
330       IF (N .EQ. NP)  GO TO 350
          DO 340 J=N1,NP
            Q=A(I-1,J)
            R=A(I,J)
            A(I-1,J)=Q*CS+R*SN
340         A(I,J)=R*CS-Q*SN
350       CONTINUE
        T(L)=zero
        T(K)=F
        S(K)=X
        GO TO 220
C   Convergence
360     IF (W .GE. zero)  GO TO 380
        S(K)=-W
        IF (NV .EQ. 0)  GO TO 380
        DO 370 J=1,N
370       V(J,K)=-V(J,K)
380     CONTINUE
C   Sort singular values
      DO 450 K=1,N
        G=-one
        J=K
        DO 390 I=K,N
          IF (S(I) .LE. G)  GO TO 390
          G=S(I)
          J=I
390       CONTINUE
        IF (J .EQ. K)  GO TO 450
        S(J)=S(K)
        S(K)=G
        IF (NV .EQ. 0)  GO TO 410
        DO 400 I=1,N
          Q=V(I,J)
          V(I,J)=V(I,K)
400       V(I,K)=Q
410     IF (NU .EQ. 0)  GO TO 430
        DO 420 I=1,N
          Q=U(I,J)
          U(I,J)=U(I,K)
420       U(I,K)=Q
430     IF (N .EQ. NP)  GO TO 450
        DO 440 I=N1,NP
          Q=A(J,I)
          A(J,I)=A(K,I)
440       A(K,I)=Q
450     CONTINUE
C   Back transformation
      IF (NU .EQ. 0)  GO TO 510
      DO 500 KK=1,N
        K=N1-KK
        IF (B(K) .EQ. zero)  GO TO 500
        Q=-A(K,K)/CdABS(A(K,K))
        DO 460 J=1,NU
460       U(K,J)=Q*U(K,J)
        DO 490 J=1,NU
          Q=dcmplx(zero,zero)
          DO 470 I=K,M
470         Q=Q+dCONJG(A(I,K))*U(I,J)
          Q=Q/(CdABS(A(K,K))*B(K))
          DO 480 I=K,M
480         U(I,J)=U(I,J)-Q*A(I,K)
490       CONTINUE
500     CONTINUE
510   IF (NV .EQ. 0)  GO TO 570
      IF (N .LT. 2)  GO TO 570
      DO 560 KK=2,N
        K=N1-KK
        K1=K+1
        IF (C(K1) .EQ. zero)  GO TO 560
        Q=-dCONJG(A(K,K1))/CdABS(A(K,K1))
        DO 520 J=1,NV
520       V(K1,J)=Q*V(K1,J)
        DO 550 J=1,NV
          Q=dcmplx(zero,zero)
          DO 530 I=K1,N
530         Q=Q+A(K,I)*V(I,J)
          Q=Q/(CdABS(A(K,K1))*C(K1))
          DO 540 I=K1,N
540         V(I,J)=V(I,J)-Q*dCONJG(A(K,I))
550       CONTINUE
560     CONTINUE
570   RETURN
      END