********************************************************************* * * EOF Analysis Based using the Subroutines in Numerical Recipes * * For CMAP Baiu Precipitation * * ic = 0 : covarinace matrix * ic = 1 : correlation matrix * * 09/17/04 by T. Tomita, Ph.D. * ********************************************************************* c program eofbaiu c integer ic,ix,it,nrot,mode real dat parameter (ic=0,dat=-999.0) parameter (ix=10,it=26,mode=10) real in0(ix,it),in1(ix,it),cov(ix,ix) real eval(ix),evec(ix,ix),escr(it,ix) c The x-value is for array of the eigenvector. c The y-value corresponds to # of mode. real evec2(ix,mode) real dmyn1(it),dmyn2(it) real dmy(ix) c open(10,file='C:\home\tomita\DATA\Baiu\baiu.pre.dat', & access='direct',form='unformatted', & recl=ix) open(20,file='C:\home\tomita\DATA\Baiu\baiu.eof.dat') c ccc Read Data c do 10 m=1,it mm=(m-1)*2 read(10,rec=mm+1) (in0(i,m),i=1,ix) read(10,rec=mm+2) (dmy(i),i=1,ix) c read(10,rec=mm+1) (dmy(i),i=1,ix) c read(10,rec=mm+2) (in0(i,m),i=1,ix) 10 continue c do 20 m=1,it icnt=0 do 21 i=1,ix icnt=icnt+1 in1(icnt,m)=in0(i,m) 21 continue 20 continue c ccc Covariance c do 30 n1=1,ix do 31 m=1,it dmyn1(m)=in1(n1,m) 31 continue do 35 n2=1,ix do 36 m=1,it dmyn2(m)=in1(n2,m) 36 continue c if(ic.eq.0) then call covar (dmyn1,dmyn2,it,dat,dcov) else if(ic.eq.1) then call coco (dmyn1,dmyn2,it,dat,dcov) end if cov(n1,n2)=dcov c 35 continue 30 continue c ccc Eigen Value Problem & Sorting c call jacobi(cov,ix,ix,eval,evec,nrot) call eigsrt(eval,evec,ix,ix) c do 39 n=1,mode icnt=0 do 40 i=1,ix icnt=icnt+1 evec2(i,n)=evec(icnt,n) 40 continue 39 continue c ccc Time Score c do 50 n=1,ix ! ix here: number of mode do 51 m=1,it ! m: parameter for time escr(m,n)=0.0 do 52 nn=1,ix ! ix here : number of component of the evec escr(m,n)=escr(m,n)+evec(nn,n)*in1(nn,m) 52 continue 51 continue 50 continue c ccc Total Variance c total=0.0 do 60 n=1,ix total=total+eval(n) 60 continue c ccc Write Data c do 70 n=1,mode write(20,*) eval(n),eval(n)/total*100.0 write(20,2000) (evec2(i,n),i=1,ix) write(20,2000) (escr(m,n),m=1,it) 2000 format(8f10.4) 70 continue c close(10) close(20) c stop end ******************************************************************* * * COVARIANCE SUBROUTINE * * --DUMMY ARGUMENTS-- * X : TIME SERIES X (REAL) * Y : TIME SERIES Y (REAL) * N : LENGTH OF X AND Y (INTEGER) * DAT : THE VALUE OF MISSING DATA (REAL) * COV : COVARIANCE OF X AND Y (REAL) * * 01/19/2000 BY T.TOMITA * ******************************************************************* C SUBROUTINE COVAR (X,Y,N,DAT,COV) C REAL X(N),Y(N),XM,YM,CNT,DAT,COV XM=0.0 YM=0.0 CNT=0.0 DO 10 I=1,N IF(ABS(X(I)-DAT).LT.0.001.OR. & ABS(Y(I)-DAT).LT.0.001) GO TO 10 XM=XM+X(I) YM=YM+Y(I) CNT=CNT+1 10 CONTINUE IF(CNT.EQ.0.0) GO TO 100 XM=XM/CNT YM=YM/CNT COV=0.0 DO 20 I=1,N IF(ABS(X(I)-DAT).LT.0.001.OR. & ABS(Y(I)-DAT).LT.0.001) GO TO 20 COV=COV+(X(I)-XM)*(Y(I)-YM) 20 CONTINUE IF((CNT-1.0).LE.0.0) GO TO 100 COV=COV/(CNT-1.0) GO TO 99 100 COV=DAT 99 RETURN END ******************************************************************* * * CORRELATION COEFFICIENT SUBROUTINE * * --DUMMY ARGUMENTS-- * X : TIME SERIES X DATA (REAL) * Y : TIME SERIES Y DATA (REAL) * N : LENGTH OF X AND Y (INTEGER) * DAT : THE VALUE OF MISSING DATA (REAL) * R : CORRELATION COEFFICIENT (REAL) * * 1990 1.7 (MON) BY T. TOMITA * ******************************************************************* C SUBROUTINE COCO (X,Y,N,DAT,R) C REAL X(N),Y(N),XM,YM,CNT,DAT,SXX,SYY,SXY,R XM=0 YM=0 CNT=0 DO 10 I=1,N IF(ABS(X(I)-DAT).LT.0.001.OR. & ABS(Y(I)-DAT).LT.0.001) GO TO 10 XM=XM+X(I) YM=YM+Y(I) CNT=CNT+1 10 CONTINUE IF(CNT.EQ.0) GO TO 100 XM=XM/CNT YM=YM/CNT SXX=0 SYY=0 SXY=0 DO 20 I=1,N IF(ABS(X(I)-DAT).LT.0.001.OR. & ABS(Y(I)-DAT).LT.0.001) GO TO 20 SXX=SXX+(X(I)-XM)**2 SYY=SYY+(Y(I)-YM)**2 SXY=SXY+(X(I)-XM)*(Y(I)-YM) 20 CONTINUE IF((SXX*SYY).LE.0.0) GO TO 100 R=SXY/SQRT(SXX*SYY) GO TO 99 100 R=DAT 99 RETURN END ********************************************************************* SUBROUTINE jacobi(a,n,np,d,v,nrot) INTEGER n,np,nrot,NMAX REAL a(np,np),d(np),v(np,np) PARAMETER (NMAX=500) INTEGER i,ip,iq,j REAL c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX) do 12 ip=1,n do 11 iq=1,n v(ip,iq)=0. 11 continue v(ip,ip)=1. 12 continue do 13 ip=1,n b(ip)=a(ip,ip) d(ip)=b(ip) z(ip)=0. 13 continue nrot=0 do 24 i=1,50 sm=0. do 15 ip=1,n-1 do 14 iq=ip+1,n sm=sm+abs(a(ip,iq)) 14 continue 15 continue if(sm.eq.0.)return if(i.lt.4)then tresh=0.2*sm/n**2 else tresh=0. endif do 22 ip=1,n-1 do 21 iq=ip+1,n g=100.*abs(a(ip,iq)) if((i.gt.4).and.(abs(d(ip))+ *g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then a(ip,iq)=0. else if(abs(a(ip,iq)).gt.tresh)then h=d(iq)-d(ip) if(abs(h)+g.eq.abs(h))then t=a(ip,iq)/h else theta=0.5*h/a(ip,iq) t=1./(abs(theta)+sqrt(1.+theta**2)) if(theta.lt.0.)t=-t endif c=1./sqrt(1+t**2) s=t*c tau=s/(1.+c) h=t*a(ip,iq) z(ip)=z(ip)-h z(iq)=z(iq)+h d(ip)=d(ip)-h d(iq)=d(iq)+h a(ip,iq)=0. do 16 j=1,ip-1 g=a(j,ip) h=a(j,iq) a(j,ip)=g-s*(h+g*tau) a(j,iq)=h+s*(g-h*tau) 16 continue do 17 j=ip+1,iq-1 g=a(ip,j) h=a(j,iq) a(ip,j)=g-s*(h+g*tau) a(j,iq)=h+s*(g-h*tau) 17 continue do 18 j=iq+1,n g=a(ip,j) h=a(iq,j) a(ip,j)=g-s*(h+g*tau) a(iq,j)=h+s*(g-h*tau) 18 continue do 19 j=1,n g=v(j,ip) h=v(j,iq) v(j,ip)=g-s*(h+g*tau) v(j,iq)=h+s*(g-h*tau) 19 continue nrot=nrot+1 endif 21 continue 22 continue do 23 ip=1,n b(ip)=b(ip)+z(ip) d(ip)=b(ip) z(ip)=0. 23 continue 24 continue pause 'too many iterations in jacobi' return END ********************************************************************* SUBROUTINE eigsrt(d,v,n,np) INTEGER n,np REAL d(np),v(np,np) INTEGER i,j,k REAL p do 13 i=1,n-1 k=i p=d(i) do 11 j=i+1,n if(d(j).ge.p)then k=j p=d(j) endif 11 continue if(k.ne.i)then d(k)=d(i) d(i)=p do 12 j=1,n p=v(j,i) v(j,i)=v(j,k) v(j,k)=p 12 continue endif 13 continue return END