Fortran做EOF分析
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
!此程序为EOF程序
! 运行时要改动前面的空间、时间格点以及文件路径,ks和kvt根据自己的需要进行改动!程序中自动去除缺省值并将其写回生成数据(生成数据中缺省值为-9999.0)
!对程序中data_in到F的传递进行调整后此程序也可用于s-eof和mv-eof
!
PROGRAM EOF
IMPLICIT NONE
INTEGER,PARAMETER :: nt=12,nx=23,ny=34 ! you need change,NT为时间长度
INTEGER,PARAMETER ::M=nt,KS=0,KVT=8 !kvt为输出的模态数
! KS的设置:ks>0 计算前先将数据标准化,
! ks=0时取距平,ks<0时不进行这一步处理
INTEGER :: i,j ,MNH,N ,K,IM , m1
REAL, allocatable,dimension(:,:,:)::DA TA_IN
REAL, allocatable,dimension(:,:)::F,S,ER,A,S1,F1
CHARACTER(LEN=20) :: NOW , TRACK
REAL :: land(nx,ny), D,A VE,PT(NX,NY,kvt) ,ran1
TRACK='E:\aat\EA\' !输出的目标文件夹,默认为程序所在文件夹
call time(now)
print*, now
!!1111111111读入数据并去掉缺省值11111111111111
ALLOCA TE(DA TA_IN(NX,NY,NT))
OPEN(1,file='E:\aat\EA\aat.eof.dat',access='direct',recl=nx*ny*nt) !****修改路径
READ(1,rec=1) (((data_in(I,J,K),I=1,nx),J=1,ny),K=1,nt)
CLOSE(1)
!注意数据排列顺序
!************做纬度加权平均,中、高纬度使用,热带或小范围不必******(未验证)!do j=1,ny
!z(j)=0.+(real(j)-1.)*2.5/180.*3.1415926575 !使用时需要改动格距和起始纬度
!data_in(:,j,:)=data_in(:,j,:)*sqrt(cos(z(j)))
!enddo
land=0.0
N=NX*NY
DO I=1,nx
DO J=1,ny
DO K=1,nt
IF(abs(data_in(I,J,K))>99999.0)then !判断缺省值(注意条件)
land(I,J)=-9999.0
N=N-1
EXIT
ENDIF
ENDDO
ENDDO
ENDDO
ALLOCA TE(F(1:N,1:M))
im=0
DO I=1,nx
DO J=1,ny
IF(land(I,J)/=-9999.0)then
im=im+1
F(IM,1:m)=data_in(I,J,1:m)
ENDIF
ENDDO
ENDDO
print*, '空间点数' , nx*ny, '非缺省值空间点数:',im,N
DEALLOCA TE(DA TA_IN)
MNH=min(N,M)
ALLOCA TE( A(MNH,MNH))
ALLOCA TE(S(MNH,MNH))
ALLOCA TE(ER(mnh,6))
!222222222222222222计算过程22222222222222222222222
CALL TRANSF(N,M,F,KS) !根据KS的设置,-1时跳出,0时距平,1时标准化
print*,"**"
CALL FORMA(N,M,MNH,F,A) !求协方差矩阵A
print*,"***"
CALL JCB(MNH,A,S,0.0000001) !雅可比过关法求特征值特征向量
print*,"****" !最后这个EPS的值控制计算精度,越小精度越高
CALL ARRANG(MNH,A,ER,S) !按照特征值大小排序
print*,"*****"
DEALLOCA TE( A)
CALL TCOEFF(KVT,N,M,MNH,S,F,ER) !给出时间序列和标准化
的空间场
print*,"******"
ALLOCA TE(S1(MNH,MNH))
ALLOCA TE(F1(N,M))
!33333********数据输出**********333333333 !输出数据为标准化后的时间序列及相应的空间场
!求时间序列的标准差,时间序列除以标准差,空间乘以该标准差
IF (M>=N) THEN
DO K=1,KVT
A VE=SUM(F(K,1:M))/REAL(M)
F1(K,1:M)=F(K,1:M)-A VE
D= SQRT(SUM(F1(K,1:M)*F1(K,1:M))/REAL(M))
F(K,1:M)= F(K,1:M)/D !时间
S(K,1:N)= S(K,1:N)*D !空间
ENDDO
m1=0
DO i=1,nx
DO j=1,ny
IF(land(i,j).eq.0.0)then
m1=m1+1
PT(i,j,1:kvt)=S(1:kvt,m1)
ELSE
PT(i,j,1:kvt)=-9999.00
ENDIF
ENDDO
ENDDO
OPEN (1,file=TRIM(track)//'pt.dat',access='direct',recl=NX*NY)
DO k=1,KVT
WRITE(1,rec=k)((PT(i,j,k),i=1,nx),j=1,ny)
ENDDO
CLOSE(1)
OPEN (2,FILE=TRIM(track)//'PC.DA T',ACCESS='DIRECT',RECL=M)
DO K=1,KVT
WRITE(2,REC=K) ((F(K,1:M)))
ENDDO
CLOSE(2)
ELSE
DO K=1,KVT
A VE=SUM(S(1:M,K))/REAL(M)
S1(1:M,K)=S(1:M,K)-A VE
D= SQRT( SUM(S1(1:M,K)*S1(1:M,K))/REAL(M))
S(1:M,K)=S(1:M,K)/D !时间
F(1:N,K)=F(1:N,K)*D !空间
ENDDO
m1=0
DO i=1,nx
DO j=1,ny
IF(land(i,j).eq.0.0)then
m1=m1+1
PT(i,j,1:kvt)=F(m1,1:kvt)
ELSE
PT(i,j,1:kvt)=-9999.00
ENDIF
ENDDO
ENDDO
OPEN (1,file=TRIM(track)//'PT.dat',access='direct',recl=NX*NY)
DO k=1,KVT
WRITE(1,rec=k) ((PT(i,j,k),i=1,nx),j=1,ny)
ENDDO
CLOSE(1)
OPEN (2,FILE=TRIM(track)//'PC.DA T',ACCESS='DIRECT',RECL=KVT) DO K=1,M
WRITE(2,REC=K) ((S(K,1:KVT)))
ENDDO
CLOSE(2)
OPEN(3,FILE=TRIM(TRACK)//'PC10.TXT')
DO K=1,M
WRITE(3,'(8F16.4)') S(K,1:KVT)
ENDDO
ENDIF
call time(now)
print*, now ,'OK!'
ENDPROGRAM
!######################################!
! !
! 以下为计算过程调用的5个子程序!
! !
!######################################!
!11111111111111111!根据KS的设置进行初步处理
SUBROUTINE TRANSF(N,M,F,KS)
IMPLICIT NONE
! THIS SUBROUTINE PROVIDES INITIAL F BY KS
INTEGER ::KS, I,M,N
REAL ::F(N,M),A VF(N),DF(N)
A VF=0.0
DF=0.0
IF(KS>0 .or. KS ==0) then !根据KS的设置,-1时跳出,0时距平,1时标准化
DO I=1,N
A VF(I)=SUM(F(I,1:M)/M)
F(I,1:M)=F(I,1:M)-A VF(I)
ENDDO
IF(KS==0) RETURN
DO I=1,N
DF(I)=SUM(F(I,1:M)*F(I,1:M))
DF(I)=SQRT(DF(I)/M)
F(I,1:M)=F(I,1:M)/DF(I)
ENDDO
ENDIF
RETURN
END
!!22222222222222222222222222222222求协方差矩阵A
SUBROUTINE FORMA(N,M,MNH,F,A)
IMPLICIT NONE
! THIS SUBROUTINE FORMS A BY F
INTEGER :: I,J,M,N,MNH
REAL :: F(N,M),A(MNH,MNH)
A=0.0
IF(M<N) THEN
DO I=1,M
DO J=I,M
A(I,J)=SUM(F(1:N,I)*F(1:N,J))
A(J,I)=A(I,J)
ENDDO
ENDDO
ELSE
DO I=1,N
DO J=I,N
A(I,J)=SUM(F(I,1:M)*F(J,1:M))
A(J,I)=A(I,J)
ENDDO
ENDDO
ENDIF
RETURN
END
!!333333333333333333333333333333333333333雅可比过关法求特征值特征向量SUBROUTINE JCB(N,A,S,EPS)
IMPLICIT NONE
! THIS SUBROUTINE COMPUTS EIGENV ALUES
! AND EIGENVECTORS OF A RETUERN S
INTEGER :: I,J,K,N,L ,I1
REAL ::A(N,N),S(N,N)
REAL :: EPS,G,S1,S2,S3,V1,V2,V3,ST,CT,IP,IQ,U,IQ1
S=0.
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)
print*,"999"
S2=EPS/FLOA T(N)*S1
S3=S1
L=0
50 S3=S3/FLOA T(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)
!PRINT*,V2*V2+U*U,1.-G*G,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
!!444444444444444444444444444444444444444按照特征值大小排序SUBROUTINE ARRANG(MNH,A,ER,S)
IMPLICIT NONE
! THIS SUBROUTINE PROVIDES A SERIES OF EIGENV ALUES ! FROM MAX TO MIN
INTEGER :: MNH,K1,K2,I ,MNH1
REAL :: A(MNH,MNH),ER(mnh,6),S(MNH,MNH)
REAL :: TR,C
TR=0.0
DO I=1,MNH
TR=TR+A(I,I)
ER(I,1)=A(I,I)
ENDDO
MNH1=MNH-1
DO K1=MNH1,1,-1
DO 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 I=1,MNH
C=S(I,K2+1)
S(I,K2+1)=S(I,K2)
S(I,K2)=C
ENDDO
ENDIF
ENDDO
ENDDO
ER(1,2)=ER(1,1)
DO I=2,mnh
ER(I,2)=ER(I-1,2)+ER(I,1)
enddo
CONTINUE
er(:,5)=er(:,1)*sqrt(2/real(mnh))
er(:,6)=er(:,5)/TR*100
ER(:,3)=ER(:,1)/TR*100
ER(:,4)=ER(:,2)/TR*100
OPEN(119,file="E:\aat\EA\eigenvalues.txt")
DO i=1,mnh
WRITE(119,'(2f16.2,4f16.6)') er(i,1),er(i,2),er(i,3),er(i,4),er(i,5),er(i,6)
ENDDO
CLOSE(119)
RETURN
END
!!555555555555555555555555555求Y!给出时间序列和标准化的空间场
SUBROUTINE TCOEFF(KVT,N,M,MNH,S,F,ER)
IMPLICIT NONE
! THIS SUBROUTINE PROVIDES EIGENVECTORS (M.GE.N,SA VED IN S;
! M.L T.N,SA VED IN F) AND ITS STANDARD TIME COEFFICENTS SERIES (M.GE.N, ! SA VED IN F; M.LT.N,SA VED IN S)
INTEGER :: i,j,k,M,N,JS ,MNH,IS,KVT
REAL :: S(MNH,MNH),F(N,M),V(MNH),ER(mnh,6) REAL :: C
DO J=1,KVT
C=0.
C=SUM(S(:,J)*S(:,J))
C=SQRT(C)
S(:,J)=S(:,J)/C
ENDDO
IF(N.LE.M) THEN
DO J=1,M
V(1:N)=F(1:N,J)
F(1:N,J)=0.
DO IS=1,KVT
F(IS,J)=SUM(V(1:N)*S(1:N,IS))
ENDDO
ENDDO
ELSE
DO I=1,N
V(1:M)=F(I,1:M)
F(I,1:M)=0.
DO JS=1,KVT
F(I,JS)=SUM(V(1:M)*S(1:M,JS))
ENDDO
ENDDO
DO JS=1,KVT
S(1:M,JS)=S(1:m,JS)*SQRT(ER(JS,1))
F(1:N,JS)=F(1:N,JS)/SQRT(ER(JS,1))
ENDDO
ENDIF
RETURN
END。