- 积分
- 1594
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-12-13
- 最后登录
- 1970-1-1
|
楼主 |
发表于 2015-12-11 19:39:12
|
显示全部楼层
parameter(nt=24)
integer::nlev,nflag !!nt为24个时次,全国117个固定测站
real::lon(117),lat(117),k(117,24),var(117)
character*8 sta(117)
character*14 filename(nt)
real ksum(117),kave(117)
integer num(117)
open(1,file='filename.txt')
do i=1,nt
read(1,*) filename(i)
!print*,'Filename:',filename(i);pause
enddo
close(1)
open(2,file='k2010.grd',form='binary')
do j=1,nt
open(3,file=filename(j))
do i=1,3
read(3,*)
end do
do i=1,117
read(3,*) sta(i),lon(i),lat(i),var(i),k(i,j)
! print *,k(i,j);pause
enddo
enddo
ksum=0.0
num=0
do i=1,117
do j=1,24
if (ABS(k(i,j)-9999).LT.1D-2) then
num(i)=num(i)+1
continue
else
ksum(i)=ksum(i)+k(i,j)
endif
enddo
kave(i)=ksum(i)/real(24-num(i))
tim=0.0;nlev=1;nflag=1
write(2) sta(i),lat(i),lon(i),tim,nlev,nflag,kave(i)
enddo
nlev=0
write(2) sta(118),lat(118),lon(118),tim,nlev,nflag
!一个时次的k场输入完毕
close(3)
close(2)
end
好的,不过贴太长怕没人看啊。
if (ABS(k(i,j)-9999).LT.1D-2) then 这个本来写的
if(k(i,j)==9999.)又怕浮点类型不能这样比值,就改成绝对值<1e-2了
那个write(2) sta(118),lat(118),lon(118),tim,nlev,nflag,想表示读取站点文件结束。还没试到那里,后面调试再继续改
|
|