- 积分
- 15777
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-5-15
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
想请问一下。我要用barnes做721暴雨中尺度滤波,我一边编一边检验发现到了求带通滤波那部分时,一开始引用子程序是求得r是没问题的,但是一进if语句,r的值就变了,全变成0了
最后程序运行时,出现了这个
麻烦各位大能帮帮看一下
程序在下面
program main
implicit none
integer,parameter::x=10,y=10
real,parameter::C1=30,C2=3000,G=0.3
real,external::distance
integer::i,j,k,l
real::r
real::hgt(x,y),F1(x,y),F2(x,y),F(x,y),D1(x,y),D2(x,y), &
E1(x,y),E2(x,y),F11(x,y),F21(x,y),F12(x,y),F22(x,y),t1,t2,w1,w2
!读高度场
open(unit=10,file='f:/hgt500.grd',form='unformatted',access='direct',recl=1,status='old')
k=0
do j=1,y
do i=1,x
k=k+1
read(unit=10,rec=k) hgt(i,j)
enddo
enddo
close(10)
! write(*,*) hgt(1:10,1:10)
!求带通滤波
do j=1,y
do i=1,x
w1=0
t1=0
w2=0
t2=0
do k=1,y
do l=1,x
r=distance(30.0+(j-1)*1.0,30.0+(k-1)*1.0,(100+i-1)*1.0,(100+l-1)*1.0)
if(r<=sqrt(C1))then
w1=w1+exp(-r**2/(4.0*C1))
t1=t1+hgt(k,l)*exp(-r**2/(4.0*C1))
print*,r
endif
if(r<=sqrt(C2))then
w2=w2+exp(-r**2/(4.0*C2))
t2=t2+hgt(k,l)*exp(-r**2/(4.0*C2))
endif
enddo
enddo
f1(i,j)=t1/w1
f2(i,j)=t2/w2
enddo
enddo
print *,'ok!!'
f=f1-f2
!输出
open(unit=10,file='f:/bl.grd',form='binary')
do j=1,y
do i=1,x
write(10) F(i,j)
enddo
enddo
close(10)
end program main
!求距离子程序
function distance(lat1,lat2,lon1,lon2)
implicit none
real,parameter::R=6371.004,pi=3.141
real::lat1,lat2,lon1,lon2
real::distance
distance=R*acos(cos(lat1*pi/180)*cos(lat2*pi/180)*cos((lon1-lon2)*pi/180)+sin(lat1*pi/180)*sin(lat2*pi/180))
return
end function distance
|
|