- 积分
- 5600
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-10-15
- 最后登录
- 1970-1-1
|
GrADS
系统平台: |
|
问题截图: |
- |
问题概况: |
运行出来的EOF1的时间系数很大,大概在50左右,但是EOF2,EOF3的时间系数大概都在正负零点几的范围 |
我看过提问的智慧: |
看过 |
自己思考时长(天): |
1 |
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
本帖最后由 yh1223 于 2014-5-28 13:18 编辑
运行出来的EOF1的时间系数很大,大概在50左右,但是EOF2,EOF3的时间系数大概都在正负零点几的范围,这个结果正常吗?还是说我的程序或者数据有问题?EOF程序:
PROGRAM EOF
C THIS PROGRAM USES EOF FOR ANALYSING TIME SERIES
C OF METEOROLOGICAL FIELD
C M:LENTH OF TIME SERIES !!!!!!!!!! m:时间序列长度
C N:NUMBER OF GRID-POINTS !!!!!!!!!! n:格点数
C KS=-1:SELF; KS=0:DEPATURE; KS=1:STANDERDLIZED DEPATURE
C KV:NUMBER OF EIGENVALUES WILL BE OUTPUT
C KVT:NUMBER OF EIGENVECTORS AND TIME SERIES WILL BE OUTPUT
C MNH=MIN(M,N)
C EGVT=EIGENVACTORS, ECOF=TIME COEFFICIENTS FOR EGVT.
C ER(KV,1)=LAMDA,LAMDA EIGENVALUE
C ER(KV,2)=ACCUMULATE LAMDA
C ER(KV,3)=THE SUM OF COMPONENTS VECTORS PROJECTED ONTO
c EIGENVACTOR.
C ER(KV,4)=ACCUMULATE ER(KV,3)
C
PARAMETER(M=35,N=144*13,MNH=35,KS=-1,KV=3,KVT=3,pi=3.1415926)
C
DIMENSION F(N,M),A(MNH,MNH),S(MNH,MNH),ER(MNH,4),
* DF(N),V(MNH),AVF(N),EGVT(N,KVT),ECOF(M,KVT),s1(3),av(3)
*
dimension h(144,13,35)
open(10,file='D:\lw\suh\data\hgt10.197811-201303.grd',
&form='binary')
open(20,file='D:\lw\suh\eof\egvt-hgt10.grd',
cform='binary')
open(30,file='D:\lw\suh\eof\t-hgt10.txt')
open(16,file='D:\lw\suh\eof\eof-hgt10.txt')
do it=1,35
do j=1,13
do i=1,144
read(10)h(i,j,it)
enddo
enddo
enddo
print*,'ok'
do it=1,35
do j=1,13
do i=1,144
F(i+(j-1)*144,it)=h(i,j,it)
enddo
enddo
enddo
cccccccccccccccccc读数据
!三维数据根据需要转化为二维
CALL TRANSF(N,M,F,AVF,DF,KS)
write(*,*)'ok program 1'
CALL FORMA(N,M,MNH,F,A)
write(*,*)'ok program 2'
CALL JCB(MNH,A,S,0.00001)
write(*,*)'ok program 3'
CALL ARRANG(KV,MNH,A,ER,S)
write(*,*)'ok program 4'
CALL TCOEFF(KVT,KV,N,M,MNH,S,F,V,ER)
write(*,*)'ok program 5'
CALL OUTER(KV,ER,MNH)
write(*,*)'ok program 6'
CALL OUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF)
write(*,*)'ok program 7'
c transf(m,n,f,ks)
c ------------------------------
c 处理资料为距平、标准方差或不变
c-----*----------------------------------------------------6---------7--
c Preprocessing data to provide a field by ks.
c input: m,n,f
c m: number of grid-points
c n: lenth of time series
c f(m,n): the raw spatial-temporal seires
c ks: contral parameter
c ks=-1: self, i.e., for the raw time series;
c ks=0: departure, i.e., for the anomaly time series from climatological mean;
c ks=1: normalized departure, i.e., for the normalized time series;
c output: f
c f(m,n): output field based on the control parameter ks.
c work variables: fw(n)
subroutine transf(m,n,f,ks)
dimension f(m,n)
dimension fw(n),wn(m) !work array
c
i0=0
do i=1,m
do j=1,n
fw(j)=f(i,j)
enddo
call meanvar(n,fw,af,sf,vf)
if(sf.eq.0.)then
i0=i0+1
wn(i0)=i
endif
enddo
c
if(i0.ne.0)then
write(*,*)'**** FAULT ****'
write(*,*)' The program cannot go on because
*the original field has invalid data.'
write(*,*)' There are totally ',i0,
* ' gridpionts with invalid data.'
write(*,*)' These invalid data are those whose standard variance
*equal zero.'
write(*,*)' The array WN stores the positions of those invalid
*grid-points. You must pick off those invalid data from the orignal
*field and then reinput a new field to calculate its EOFs.'
write(*,*)'**** FAULT ****'
stop
endif
c
c
if(ks.eq.-1)return
c
if(ks.eq.0)then !anomaly of f
do i=1,m
do j=1,n
fw(j)=f(i,j)
enddo
call meanvar(n,fw,af,sf,vf)
do j=1,n
f(i,j)=f(i,j)-af
enddo
enddo
return
endif
c
if(ks.eq.1)then !normalizing f
do i=1,m
do j=1,n
fw(j)=f(i,j)
enddo
call meanvar(n,fw,af,sf,vf)
do j=1,n
f(i,j)=(f(i,j)-af)/sf
enddo
enddo
endif
return
end
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
ccccccccccccc存储数据
do j=1,3
do i=1,m
av(j)=av(j)+ecof(i,j)/m
enddo
do i=1,m
s1(j)=s1(j)+(ecof(i,j)-av(j))**2/m
enddo
s1(j)=sqrt(s1(j))
do i=1,m
ecof(i,j)=ecof(i,j)/(s1(j))
enddo
enddo
do j=1,m
write(30,*) (ecof(j,i),i=1,3)
enddo
close(30)
do it=1,kvt
do j=1,n
egvt(j,it)=egvt(j,it)*s1(it)
enddo;enddo
do it=1,kvt
do j=1,n
write(20)egvt(j,it)
enddo;enddo
close(20)
write(*,*)'ok 8'
cccccccccccc
END
ccccccccccccccccccccccccc子程序
SUBROUTINE TRANSF(N,M,F,AVF,DF,KS)
C THIS SUBROUTINE PROVIDES INITIAL F BY KS
DIMENSION F(N,M),AVF(N),DF(N)
DO 5 I=1,N
AVF(I)=0.0
5 DF(I)=0.0
IF(KS) 30,10,10
10 DO 14 I=1,N
DO 12 J=1,M
12 AVF(I)=AVF(I)+F(I,J)
AVF(I)=AVF(I)/M
DO 14 J=1,M
F(I,J)=F(I,J)-AVF(I)
14 CONTINUE
IF(KS.EQ.0) THEN
RETURN
ELSE
DO 24 I=1,N
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
F(I,J)=F(I,J)/DF(I)
24 CONTINUE
ENDIF
30 CONTINUE
RETURN
END
SUBROUTINE FORMA(N,M,MNH,F,A)
C THIS SUBROUTINE FORMS A BY F
DIMENSION F(N,M),A(MNH,MNH)
IF(M-N) 40,50,50
40 DO 44 I=1,MNH
DO 44 J=I,MNH
A(I,J)=0.0
DO 42 IS=1,N
42 A(I,J)=A(I,J)+F(IS,I)*F(IS,J)
A(J,I)=A(I,J)
44 CONTINUE
RETURN
50 DO 54 I=1,MNH
DO 54 J=I,MNH
A(I,J)=0.0
DO 52 JS=1,M
52 A(I,J)=A(I,J)+F(I,JS)*F(J,JS)
A(J,I)=A(I,J)
54 CONTINUE
RETURN
END
SUBROUTINE JCB(N,A,S,EPS)
C THIS SUBROUTINE COMPUTS EIGENVALUES AND standard 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,MNH,A,ER,S)
C THIS SUBROUTINE PROVIDES A SERIES OF EIGENVALUES
C FROM MAX TO MIN
DIMENSION A(MNH,MNH),ER(MNH,4),S(MNH,MNH)
TR=0.0
DO 200 I=1,MNH
TR=TR+A(I,I)
200 ER(I,1)=A(I,I)
MNH1=MNH-1
DO 210 K1=MNH1,1,-1
DO 210 K2=K1,MNH1
IF(ER(K2,1).LT.ER(K2+1,1)) THEN
C=ER(K2+1,1)
ER(K2+1,1)=ER(K2,1)
ER(K2,1)=C
DO 205 I=1,MNH
C=S(I,K2+1)
S(I,K2+1)=S(I,K2)
S(I,K2)=C
205 CONTINUE
ENDIF
210 CONTINUE
ER(1,2)=ER(1,1)
DO 220 I=2,KV
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
WRITE(*,250) TR
250 FORMAT(/5X,'TOTAL SQUARE ERROR=',F20.5)
RETURN
END
SUBROUTINE TCOEFF(KVT,KV,N,M,MNH,S,F,V,ER)
C THIS SUBROUTINE PROVIDES STANDARD EIGENVECTORS (M.GE.N,SAVED IN S;
C M.LT.N,SAVED IN F) AND ITS TIME COEFFICENTS SERIES (M.GE.N,
C SAVED IN F; M.LT.N,SAVED IN S)
DIMENSION S(MNH,MNH),F(N,M),V(MNH),ER(MNH,4)
IF(N.LE.M) THEN
DO 390 J=1,M
DO 370 I=1,N
V(I)=F(I,J)
F(I,J)=0.
370 CONTINUE
DO 380 IS=1,KVT
DO 380 I=1,N
380 F(IS,J)=F(IS,J)+V(I)*S(I,IS)
390 CONTINUE
ELSE
DO 410 I=1,N
DO 400 J=1,M
V(J)=F(I,J)
F(I,J)=0.
400 CONTINUE
DO 410 JS=1,KVT
DO 410 J=1,M
F(I,JS)=F(I,JS)+V(J)*S(J,JS)
410 CONTINUE
DO 430 JS=1,KVT
DO 420 J=1,M
S(J,JS)=S(J,JS)*SQRT(ER(JS,1))
420 CONTINUE
DO 430 I=1,N
F(I,JS)=F(I,JS)/SQRT(ER(JS,1))
430 CONTINUE
ENDIF
RETURN
END
SUBROUTINE OUTER(KV,ER,MNH)
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(MNH,4)
WRITE(16,510)
510 FORMAT(/10X,'EIGENVALUE AND ANALYSIS ERROR')
WRITE(16,520)
520 FORMAT(10X,1HH,8X,5HLAMDA,10X,6HSLAMDA,11X,2HPH,12X,3HSPH)
WRITE(16,530) (IS,(ER(IS,J),J=1,4),IS=1,KV)
530 FORMAT(1X,I10,4F15.5)
WRITE(16,540)
540 FORMAT(//)
RETURN
END
SUBROUTINE OUTVT(KVT,N,M,MNH,S,F,EGVT,ECOF)
C THIS SUBROUTINE PRINTS STANDARD EIGENVECTORS
C AND ITS TIME-COEFFICENT SERIES
DIMENSION F(N,M),S(MNH,MNH),EGVT(N,KVT),ECOF(M,KVT)
WRITE(16,560)
560 FORMAT(10X,'STANDARD EIGENVECTORS')
WRITE(16,570) (IS,IS=1,KVT)
570 FORMAT(3X,10i7)
DO 550 I=1,N
IF(M.GE.N) THEN
WRITE(16,580) I,(S(I,JS),JS=1,KVT)
580 FORMAT(1X,I3,10F7.3,/)
DO 11 JS=1,KVT
EGVT(I,JS)=S(I,JS)
11 CONTINUE
ELSE
WRITE(16,590) I,(F(I,JS),JS=1,KVT)
590 FORMAT(1X,I5,10F7.3)
DO 12 JS=1,KVT
EGVT(I,JS)=F(I,JS)
12 CONTINUE
ENDIF
550 CONTINUE
C WRITE(16,590) I,(F(I,JS),JS=1,KVT)
! WRITE(20)((F(I,JS),i=1,n),JS=1,KVT)
WRITE(16,720)
720 FORMAT(//)
WRITE(16,610)
610 FORMAT(10X,'TIME-COEFFICENT SERIES OF S. E.')
WRITE(16,620) (IS,IS=1,KVT)
620 FORMAT(3X,5i12)
DO 600 J=1,M
IF(M.GE.N) THEN
WRITE(16,630) J,(f(is,j),is=1,kvt)
630 FORMAT(1X,I3,5F12.3)
DO 13 IS=1,KVT
ECOF(J,IS)=F(IS,J)
13 CONTINUE
ELSE
WRITE(16,640) J,(S(J,IS),IS=1,KVT)
640 FORMAT(1X,I3,10F12.3)
DO 14 IS=1,KVT
ECOF(J,IS)=S(J,IS)
14 CONTINUE
ENDIF
600 CONTINUE
C WRITE(30)((S(J,IS),j=1,m),IS=1,KVT)
RETURN
END
32
|
|