| 
 
	积分2210贡献 精华在线时间 小时注册时间2011-9-15最后登录1970-1-1 
 | 
 
 
 楼主|
发表于 2013-11-1 18:07:56
|
显示全部楼层 
| 本帖最后由 lqouc 于 2013-11-1 18:28 编辑 
 附上程序。
 ****相关*****文件1为单纯的一列数据,文件2为空间点,每个点都有同样长时间序列的数据*********
 parameter(nt=34,nx=144,ny=73)
 !nt::time series, nx:X方向格点数,ny:Y方向格点数
 real f(nt),ff(nx,ny,nt)    !分别存放两个原始数据
 real Y(nt)
 open(1,file='sjxs.prewinter.grd',form='binary')    !!!!!!!!!
 do i=1,nt
 read(1) f(i)
 print*,f(i)
 enddo
 close(1)          !读入文件1
 open(2,file='chi-summer200.grd',
 $form='binary') !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 do k=1,nt
 do j=1,ny
 do i=1,nx
 read(2) ff(i,j,k)
 ! print*,ff(i,j,k)
 enddo
 enddo
 enddo
 close(2)          !读入文件2
 open(5,file='r.grd',form='binary') !结果输出
 open(6,file='t.grd',form='binary')
 
 do 50 jj=1,ny
 do 40 ii=1,nx
 do t=1,nt
 Y(t)=ff(ii,jj,t)
 enddo
 call correlation(f,Y,nt,r,tt)
 write(5) r
 write(6) tt
 40    continue
 50    continue
 close(5)
 close(6)
 end
 ccc计算相关系数  CCCCCCCC
 subroutine correlation(X,Y,nt,RHO,TTT)
 real X(nt),Y(nt)
 real*8 SX,SY,SXX,SYY,SXY,FN
 SX=0.0
 SY=0.0
 SXX=0.0
 SYY=0.0
 SXY=0.0
 FN=0.0
 do 20 i=1,nt
 if(Y(i).ne.-9.99E+33)then
 FN=FN+1
 SX=SX+X(i)
 SY=SY+Y(i)
 SXX=SXX+X(i)*X(i)
 SYY=SYY+Y(i)*Y(i)
 SXY=SXY+X(i)*Y(i)
 endif
 20    continue
 if(FN.ne.0)then
 RHO=(SXY-SX*SY/FN)/sqrt((SXX-SX*SX/FN)
 &  *(SYY-SY*SY/FN))
 else
 RHO=-9.99E+33
 endif
 if(RHO.ne.-9.99E+33)then
 TTT=RHO*sqrt(30.0)/sqrt(1-RHO*RHO)
 else
 TTT=-9.99E+33
 endif
 end
 
 
 | 
 |