program pca_noamer c c c reduce north american dataset 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 = 1500) parameter (nmax0=1500) 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*40(maxproxy) character *20 nome0 character *20 nome(maxproxy) character *24 nome2(maxproxy) character *48 flnm2 C C My compiler says these are used without being defined. Give C them values that the probably get from default initialisation, just C to be sure. But with my compiler, it makes no difference. iabv=0 iunif=0 C 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 = 1000 nhigh = 1980 iproxmax = nhigh ntot = nhigh-nlow+1 c do i=1,maxproxy istart(i)=99999 end do c open (unit=1, $ file='wmc.inf', $ status='old') C C Read this in from the .inf file read(unit=1,fmt=*) itrainmin0,iproxmin print *,itrainmin0,iproxmin C 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) endif goto 1010 415 continue nproxy = i-1 close (unit=1) C C c c define limits for raw data and training interval c C itrainmin0 = 1902 ! Now read from file itrainmin = itrainmin0 itrainmax0 = nhigh itrainmax = itrainmax0 ntrain = itrainmax0-itrainmin0+1 iyearmax = 1993 C irawmin = 1902 irawmin = itrainmin0 C iearlyinst = 1820 ! Not used - remove for clarity irawmax = nhigh nraw = irawmax-irawmin+1 C C C 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:' C C Now read in. Write to verify C read (5,*) iproxmin print *,iproxmin C c c c select proxy records to be used c iset=0.0 !wmc if (iset.eq.0) then mkeep = 0 do i=1,nproxy if (istart(i).le.iproxmin) then mkeep = mkeep + 1 write (6,222) mkeep, flnm(i) ikeep(mkeep)=iabv+i endif end do write (6,*) 'total of ',mkeep, $ ' training records back to ',iproxmin endif 222 format (i5,3x,a30) 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