c$lager c program svdsv154 C THIS PROGRAM USES SVD FOR ANALYSING relationship between C a scalar FIELD series and a vector FIELD series C M:LENTH OF TIME SERIES C N1:GRID-POINTS NUMBER of scalar field f; C N2:2 times of GRID-POINTS NUMBER of vector field g; C KS=-1:SELF; KS=0:DEPATURE; KS=1:STANDaRD DEPATURE.But only used c in ks=o or 1. c kv=1:standard g for vector V;kv=2:standard g for scalar u,v. C KVl:NUMBER OF S.VALUES WILL BE OUTPUT C KVT:NUMBER OF S.VECTORS AND TIME SERIES WILL BE OUTPUT c ff,gg is undefine value of f,g. c evf,evg=s.v. for f,g; c tcf,tcg=time coefficents for evf,evg. c r=correlation between tcf and tcg C ER(KVl,1)=EIGENVALUE lambda. C ER(KVl,2)=ACCUMULATE LAMbDA C ER(KVl,3)=THE SUM OF COMPONENTS VECTORS PROJECTED ONTO EIGENVACTOR. C ER(KVl,4)=ACCUMULATE ER(KVl,3) c *************************************************************** c * This program used SVD between SSTA and SLP anomaly * c * Region: 20n-60n,dLon*dLat=5*4,Nx*Ny=18*12=216. * c * Time:1955.1.-1998.1, 1yr.,M=44. * c * Climatic standard depature. * c *************************************************************** PARAMETER(n1=275,n2=297,min12=275,m=44,KS=-1,KVl=10,KVT=10,mt=100) PARAMETER(ff=9.99E+33,gg=-999.0,nx=25,ny=11,nnx=33,nny=9) DIMENSION f(n1,m),g(n2,m) dimension nf(n1),ng(n2) dimension fg(n1,n2) c **************************** dimension a(min12,min12),s(min12,min12),er(kvl,4) c **************************** dimension avf(n1),df(n1) dimension evf(n1,kvt),tcf(m,kvt) c **************************** dimension avg(n2),dg(n2) dimension evg(n2,kvt),tcg(m,kvt) c **************************** dimension pf(kvt),pg(kvt),r(kvt) c **************************** dimension ersqrt(kvl,4),ermo(100,kvl) dimension data1(nx,ny),data2(nnx,nny) CCCC INPUT DATA f(n1,m),g(n2,m) ccccccccccccccccccc c open(11,file='d:\ptsvd2\grd\sjp1.grd',form='unformatted', c !access='direct',recl=nx*ny*4) c do 132 it=1,m c read(11,rec=it)((data1(i,j),i=1,nx),j=1,ny) c do 132 jj=1,ny c do 132 ii=1,nx c kkkk=nx*(jj-1)+ii c f(kkkk,it)=data1(ii,jj) c132 continue c close(11) open (16,file='e:\svd1\mtls.dat') ix=133 s=1.0 am=0.0 do 1000 kk=1,mt do 7 i=1,n1 do 7 it=1,m call rands2(ix,yfl,s,am,v) 7 f(i,it)=v CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC open(11,file='e:\svd1\npslp1.grd',form='unformatted', !access='direct',recl=nnx*nny) do 133 it=1,m read(11,rec=it)((data2(i,j),i=1,nnx),j=1,nny) do 133 jj=1,nny do 133 ii=1,nnx kkkk=nnx*(jj-1)+ii g(kkkk,it)=data2(ii,jj) 133 continue close(11) CCCCCCCCCCCCCCCC The end of INPUT DATA CCCCCCCCCCCCCCCCCCC call test1(n1,n2,m,ff,gg,f,g,nf,ng) call transf(n1,n2,m,f,g,nf,ng,avf,avg,df,dg,ks) call forma(n1,n2,min12,m,f,g,fg,a) c call cms(n1,m,f,sumf) c call cms(n2,m,g,sumg) call jcb(min12,a,s,0.0001) c when n1=min12,s is s.vectors of f; c when n2=min12,s is s.vectors of g. if (min12.eq.n1) then CALL ARRANG(KVl,kvt,N1,A,ER,S,evf) else CALL ARRANG(KVl,kvt,N2,A,ER,S,evg) endif c CALL test2(kvt,n1,n2,min12,fg,s,a,evf,evg) c call ctc(kvt,n1,m,f,evf,tcf) c call ctc(kvt,n2,m,g,evg,tcg) c call outer(kvl,er) c CALL test3(n1,n2,ff,gg,nf,ng,evf,evg,kvt) c call outvt(kvt,n1,n2,m,evf,evg,tcf,tcg) c call cpr(m,kvt,sumf,sumg,tcf,tcg,pf,pg,r) do 110 j=1,kvl 110 ermo(kk,j)=er(j,3) write(16,179)(ermo(kk,j),j=1,kvl) 179 format(10f7.3) cccccccccccccccccccccccccc 1000 continue c********************************************* c call test4(kvl,min12,A,ERsqrt) c call OUTER(KVL,ERsqrt) c********************************************* write (16,*) '------------' do j=1,kvl c=0 do k1=mt,1,-1 do k2=k1,mt IF(ermo(K2,j).LT.ermo(K2+1,j)) THEN C=ermo(K2+1,j) ermo(K2+1,j)=ermo(K2,j) ermo(K2,j)=C endif enddo enddo enddo do k=1,6 c 输出降序排列后最大的6个 write(16,179)(ermo(k,j),j=1,kvl) enddo close (16) stop END c********************************************* subroutine rands2(ix,yfl,s,am,v) a=0.0 do 1 i=1,12 call randu2(ix,yfl) 1 a=a+yfl v=(a-6.0)*s+am return end subroutine randu2(ix,yfl) if(ix.eq.0)ix=67107 ix=125*ix ix=ix-ix/2796203*2796203 yfl=float(ix) yfl=yfl/2796203.0 return end c********************************************* SUBROUTINE Test1(n1,n2,m,ff,gg,f,g,nf,ng) DIMENSION F(N1,M),g(n2,m) DIMENSION nF(N1),ng(n2) do i=1,n1 nf(i)=0.0 enddo do i=1,n2 ng(i)=0.0 enddo do i=1,n1 do k=1,m if(f(i,k).eq.ff)then f(i,k)=0.0 nf(i)=nf(i)+1 endif enddo enddo do i=1,n2 do k=1,m if(g(i,k).eq.gg)then g(i,k)=0.0 ng(i)=ng(i)+1 endif enddo enddo return end SUBROUTINE TRANSF(n1,n2,m,f,g,nf,ng,avf,avg,df,dg,ks) C THIS SUBROUTINE PROVIDES INITIAL F and g BY KS and kv. DIMENSION F(N1,M),g(n2,m),AVF(N1),avg(n2) DIMENSION DF(N1),dg(n2) DIMENSION nF(N1),ng(n2) if(ks.eq.-1)then goto 30 endif do i=1,n1 avf(i)=0.0 enddo do i=1,n2 avg(i)=0.0 enddo if(ks.eq.0)then goto 5 endif do i=1,n1 df(i)=0.0 enddo do i=1,n2 dg(i)=0.0 enddo cccccccccccccccccc 5 continue DO 141 I=1,N1 if(nf(i).ne.0) goto 141 do 12 j=1,m 12 AVF(I)=AVF(I)+F(I,J) AVF(I)=AVF(I)/M DO 14 J=1,M 14 F(I,J)=F(I,J)-AVF(I) 141 CONTINUE DO 181 I=1,N2 if(ng(i).ne.0) goto 181 DO 16 J=1,M 16 AVg(I)=AVg(I)+g(I,J) AVg(I)=AVg(I)/M DO 18 J=1,M 18 g(I,J)=g(I,J)-AVg(I) 181 CONTINUE IF(KS.EQ.0) THEN RETURN ELSE DO 241 I=1,N1 if(nf(i).ne.0) goto 241 DO 22 J=1,M 22 DF(I)=DF(I)+F(I,J)*F(I,J) DF(I)=SQRT(DF(I)/M) DO 24 J=1,M 24 F(I,J)=F(I,J)/DF(I) 241 CONTINUE DO 281 I=1,N2 if(ng(i).ne.0) goto 281 DO 26 J=1,M 26 Dg(I)=Dg(I)+g(I,J)*g(I,J) Dg(I)=SQRT(Dg(I)/M) DO 28 J=1,M 28 g(I,J)=g(I,J)/Dg(I) 281 CONTINUE ENDIF 30 CONTINUE RETURN END SUBROUTINE FORMA(N1,n2,min12,M,f,g,fg,a) C THIS SUBROUTINE FORMS fg and a bY f and g DIMENSION F(N1,M),g(n2,m) dimension fg(n1,n2) dimension a(min12,min12) DO I=1,n1 DO J=1,n2 fg(I,J)=0.0 DO jt=1,m fg(I,J)=fg(I,J)+F(i,jt)*g(j,jt) enddo enddo enddo if(n1.eq.min12)then do i=1,n1 do j=i,n1 a(i,j)=0.0 do k=1,n2 a(i,j)=a(i,j)+fg(i,k)*fg(j,k) enddo a(j,i)=a(i,j) enddo enddo else do i=1,n2 do j=i,n2 a(i,j)=0.0 do k=1,n1 a(i,j)=a(i,j)+fg(k,i)*fg(k,j) enddo a(j,i)=a(i,j) enddo enddo endif RETURN END subroutine cms(n,m,f,sum) c This subrutine computs total variance of F or G. dimension f(n,m) sum=0.0 do i=1,n do j=1,m sum=sum+f(i,j)*f(i,j) enddo enddo return end SUBROUTINE JCB(N,A,S,EPS) C THIS SUBROUTINE COMPUTS EIGENVALUES AND EIGENVECTORS OF A DIMENSION A(N,N),S(N,N) DO 30 I=1,N DO 30 J=1,I IF(I-J) 20,10,20 10 S(I,J)=1. GO TO 30 20 S(I,J)=0. S(J,I)=0. 30 CONTINUE G=0. DO 40 I=2,N I1=I-1 DO 40 J=1,I1 40 G=G+2.*A(I,J)*A(I,J) S1=SQRT(G) S2=EPS/FLOAT(N)*S1 S3=S1 L=0 50 S3=S3/FLOAT(N) 60 DO 130 IQ=2,N IQ1=IQ-1 DO 130 IP=1,IQ1 IF(ABS(A(IP,IQ)).LT.S3) GOTO 130 L=1 V1=A(IP,IP) V2=A(IP,IQ) V3=A(IQ,IQ) U=0.5*(V1-V3) IF(U.EQ.0.0) G=1. IF(ABS(U).GE.1E-10) G=-SIGN(1.,U)*V2/SQRT(V2*V2+U*U) ST=G/SQRT(2.*(1.+SQRT(1.-G*G))) CT=SQRT(1.-ST*ST) DO 110 I=1,N G=A(I,IP)*CT-A(I,IQ)*ST A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT A(I,IP)=G G=S(I,IP)*CT-S(I,IQ)*ST S(I,IQ)=S(I,IP)*ST+S(I,IQ)*CT 110 S(I,IP)=G DO 120 I=1,N A(IP,I)=A(I,IP) 120 A(IQ,I)=A(I,IQ) G=2.*V2*ST*CT A(IP,IP)=V1*CT*CT+V3*ST*ST-G A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST) A(IQ,IP)=A(IP,IQ) 130 CONTINUE IF(L-1) 150,140,150 140 L=0 GO TO 60 150 IF(S3.GT.S2) GOTO 50 RETURN end SUBROUTINE ARRANG(KV,kvt,N,A,ER,S,ev) C THIS SUBROUTINE PROVIDES A SERIES OF EIGENVALUES C FROM MAX TO MIN,and a series of standard S.vectors_ev. DIMENSION A(N,N),ER(KV,4),S(N,N),ev(n,kvt) TR=0.0 DO 200 I=1,N 200 TR=TR+A(I,I) N1=N-1 DO 210 K1=N1,1,-1 DO 210 K2=K1,N1 IF(a(K2,k2).LT.a(K2+1,k2+1)) THEN C=a(K2+1,k2+1) a(K2+1,k2+1)=a(K2,k2) a(K2,k2)=C DO 205 I=1,N C=S(I,K2+1) S(I,K2+1)=S(I,K2) S(I,K2)=C 205 CONTINUE ENDIF 210 CONTINUE c************* standard S.vectors *************** do k=1,kvt c=0.0 do i=1,n c=c+s(i,k)*s(i,k) enddo c=sqrt(c) do i=1,n ev(i,k)=s(i,k)/c enddo enddo c**************************** ER(1,1)=a(1,1) ER(1,2)=a(1,1) DO 220 I=2,KV er(i,1)=a(i,i) ER(I,2)=ER(I-1,2)+ER(I,1) 220 CONTINUE DO 230 I=1,KV ER(I,3)=ER(I,1)/TR ER(I,4)=ER(I,2)/TR 230 CONTINUE c WRITE(9,250) TR c WRITE(6,250) TR 250 FORMAT(/5X,'TOTAL SQUARE ERROR=',F20.5) RETURN END subroutine test2(kvt,n1,n2,min12,fg,s,a,evf,evg) c This subroutine computs long standard s.vectors c by short standard s.vectors. dimension fg(n1,n1),s(min12,min12),a(min12,min12) dimension evf(n1,kvt),evg(n2,kvt) if(n1.eq.min12)then do k=1,kvt do i=1,n1 enddo do i=1,n2 evg(i,k)=0.0 do j=1,n1 evg(i,k)=evg(i,k)+fg(j,i)*evf(j,k)/sqrt(a(k,k)) enddo enddo enddo else do k=1,kvt do i=1,n2 evg(i,k)=s(i,k) enddo do i=1,n1 evf(i,k)=0.0 do j=1,n2 evf(i,k)=evf(i,k)+fg(i,j)*evg(j,k)/sqrt(a(k,k)) enddo enddo enddo endif c write(6,*)((evf(i,k),k=1,kvt),i=1,n1) c write(6,*)((evg(i,k),k=1,kvt),i=1,n2) return end SUBROUTINE ctc(KVt,N,m,z,ev,tc) C THIS SUBROUTINE computs time coefficents of standard s.vector_ev. dimension z(n,m),ev(n,kvt) dimension tc(m,kvt) do k=1,kvt do j=1,m tc(j,k)=0.0 do i=1,n tc(j,k)=tc(j,k)+z(i,j)*ev(i,k) enddo enddo enddo return end subroutine cpr(m,kvt,sumf,sumg,tcf,tcg,pf,pg,r) dimension tcf(m,kvt),tcg(m,kvt) dimension pf(kvt),pg(kvt),r(kvt) do k=1,kvt pf(k)=0.0 pg(k)=0.0 r(k)=0.0 do it=1,m pf(k)=pf(k)+tcf(it,k)*tcf(it,k) pg(k)=pg(k)+tcg(it,k)*tcg(it,k) r(k)=r(k)+tcf(it,k)*tcg(it,k) enddo c write(6,*)pf(k),pg(k),r(k) r(k)=r(k)/sqrt(pf(k))/sqrt(pg(k)) pf(k)=pf(k)/sumf pg(k)=pg(k)/sumg c write(6,*)pf(k),pg(k),r(k) enddo return end SUBROUTINE OUTER(KVL,ER) C THIS SUBROUTINE PRINTS ARRAY ER C ER(KV,1) FOR SEQUENCE OF EIGENVALUE FROM BIG TO SMALL C ER(KV,2) FOR EIGENVALUE FROM BIG TO SMALL C ER(KV,3) FOR SMALL LO=(LAMDA/TOTAL VARIANCE) C ER(KV,4) FOR BIG LO=SUM OF SMALL LO) DIMENSION ER(KVL,4) WRITE(6,510) 510 FORMAT(/20X,'EIGENVALUE AND ANALYSIS ERROR1') WRITE(6,520) 520 FORMAT(10X,1HH,8X,5HLAMDA,10X,6HSLAMDA,11X,2HPH,12X,3HSPH) WRITE(6,530) (IS,(ER(IS,J),J=1,4),IS=1,KVL) 530 FORMAT(1X,I10,4F15.5) WRITE(6,540) 540 FORMAT(//) RETURN END SUBROUTINE test3(N1,n2,ff,gg,nf,ng,evf,evg,kvt) c this subroutine sent undefine value ff,gg to evf,evg. dimension nf(n1),ng(n2),evf(n1,kvt),evg(n2,kvt) do i=1,n1 if(nf(i).ne.0)then do k=1,kvt evf(i,k)=ff enddo endif enddo do i=1,n2 if(ng(i).ne.0)then do k=1,kvt evg(i,k)=gg enddo endif enddo return end SUBROUTINE OUTVT(KVT,N1,n2,M,evf,evg,tcf,tcg) C THIS SUBROUTINE PRINTS STANDARD EIGENVECTORS C AND ITS TIME-COEFFICENT SERIES DIMENSION evf(n1,kvt),evg(n2,kvt),tcf(m,kvt),tcg(m,kvt) WRITE(6,560) 560 FORMAT(20X,'Left STANDARD EIGENVECTORS') WRITE(6,570) (IS,IS=1,KVT) 570 FORMAT(13X,15I12) DO 550 I=1,N1 WRITE(6,580) I,(evf(I,JS),JS=1,KVT) 580 FORMAT(10X,I3,15F12.3,/) 550 continue WRITE(6,720) 720 FORMAT(//) WRITE(6,610) 610 FORMAT(20X,'TIME-COEFFICENT SERIES OF L. S. E.') WRITE(6,620) (IS,IS=1,KVT) 620 FORMAT(13X,15I12) DO 600 J=1,M WRITE(6,630) J,(tcF(j,IS),IS=1,KVT) 630 FORMAT(10X,I3,15F12.3) 600 continue write(6,640) 640 format(//) WRITE(6,1560) 1560 FORMAT(20X,'Right STANDARD EIGENVECTORS') WRITE(6,1570) (IS,IS=1,KVT) 1570 FORMAT(13X,15I12) DO 1550 I=1,N2 WRITE(6,1580) I,(evg(I,JS),JS=1,KVT) 1580 FORMAT(10X,I3,15F12.3,/) 1550 continue WRITE(6,1720) 1720 FORMAT(//) WRITE(6,1610) 1610 FORMAT(20X,'TIME-COEFFICENT SERIES OF R. S. E.') WRITE(6,1620) (IS,IS=1,KVT) 1620 FORMAT(13X,15I12) DO 1600 J=1,M WRITE(6,1630) J,(tcg(j,IS),IS=1,KVT) 1630 FORMAT(10X,I3,15F12.3) 1600 cONTINUE write(6,1800) 1800 format(//) RETURN END subroutine outspr(kvt,sumf,sumg,pf,pg,r) dimension pf(kvt),pg(kvt),r(kvt) write(6,1810) 1810 format(20x,'total square of departual or standard depatual') write(6,1820)sumf 1820 format(20x,'total square of F =',15f12.3) write(6,1830)sumg 1830 format(20x,'total square of G =',15f12.3) write(6,1840) 1840 format(//) write(6,1850) 1850 format(20x,'p and r') write(6,1860)(is,is=1,kvt) 1860 format(10x,15i12) write(6,1870)(pf(is),is=1,kvt) 1870 format(9x,3hp.L,15f12.3) write(6,1880)(pg(is),is=1,kvt) 1880 format(9x,3hp.R,15f12.3) 1890 format(10x,1hr,1x,15f12.3) write(6,1900) 1900 format(//) return end cccccccccccccccccccccccccccccccccccccccccccccccc subroutine test4(kvL,n,A,ER) c This subroutine changes lambda to sqrt(lambda) and computs er. DIMENSION A(N,N),ER(KVL,4) TR=0.0 DO 200 I=1,N 200 TR=TR+A(I,I) tr=sqrt(tr) c********************* do i=1,n if (a(i,i).LT.0.)then a(i,i)=0. else a(i,i)=sqrt(a(i,i)) endif enddo c********************* ER(1,1)=a(1,1) ER(1,2)=a(1,1) DO 220 I=2,KVL er(i,1)=a(i,i) ER(I,2)=ER(I-1,2)+ER(I,1) 220 CONTINUE DO 230 I=1,KVL ER(I,3)=ER(I,1)/TR ER(I,4)=ER(I,2)/TR 230 CONTINUE WRITE(9,250) TR WRITE(6,250) TR 250 FORMAT(/5X,'TOTAL model of FG=',F20.5) RETURN END