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