- 积分
- 2926
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2012-3-3
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 lotina 于 2013-4-17 12:53 编辑
❤----------------❤--------------------------------❤--------------------------------❤
资料是730站逐日（8天）降水数据，经过缺测处理，有数值！
根据台站转格点数据画图的步骤一步步编程画图的！但是画出来的图形数值为零！
求问原因，是不是哪里错了， 画出来的图形数值为零？ 我知道看程序很让人头疼，但还是期待大家的帮助0.0.坐等解答，万分感谢！
程序如下：
（1）站点转格点FOR 时次为8天：
! -----------------------------------------------------------
!This program is used to change the format of ascii data into
! the form of binary, which is supported by the GrADS
!------------------------------------------------------------
program STGRD
implicit none
integer,Parameter::phas=8,sn=730
integer i,j,k,td,ty,st1,st730,conter,ii
real lat(sn),lon(sn)
!integer stn(sn)
integer rgq(730,8),rgj(730,8) !------------->改为real
character(len=5) :: station ------------->改为LEN=8
integer it,nflag,nlev
real tim
character*8 stn(sn)
st1=50136;st730=59985
!--------------读入730 台站号- 经纬度------------------
open(1,file='D:\DATA\RAIN\daily\station-730-1.txt',status='old')
do i=1,730
read(1,*) stn(i),lat(i),lon(i)
end do
close(1)
print*, stn(1),lat(1),lon(1),stn(sn),lat(sn),lon(sn)
!---------------读入-位相（强年 偏西）合成后的数据-----------
open(3,file='D:\DATA\RAIN\rain.station\com\8com-gjq.txt')
read(3,*) ((rgj(i,td),i=1,730),td=1,8)
close(3)
open(4,file='D:\DATA\RAIN\rain.station\com\8com-gqq.txt')
read(4,*) ((rgq(i,td),i=1,730),td=1,8)
close(4)
print*, '-------------read over-----------'
print*, rgq
!----将文本记录转换为GRADS所支持的二进制记录-----------
open(40,file='D:\STUDY\WPSH-01\Rain\raintemp\rain.gq.grd',form='binary')
do it=1,phas
nlev=1
tim=0.0
nflag=1
do i=1,SN
write(40) stn(i),lat(i),lon(i),tim,nlev,nflag,rgq(i,it)
enddo
nlev=0
write(40) stn(sn),lat(sn),lon(sn),tim,nlev,nflag
enddo
close(40) ...........................................................................................................
open(40,file='D:\STUDY\WPSH-01\Rain\raintemp\rain.gq.txt')
do it=1,phas
nlev=1
tim=0.0
nflag=1
do i=1,SN
write(40,*) stn(i),lat(i),lon(i),tim,nlev,nflag,rgq(i,it)
enddo
enddo
close(40)
end program
（2）降水资料的描述文件：
dset D:\STUDY\WPSH-01\Rain\raintemp\rain.gq.grd
dtype station
stnmap D:\STUDY\WPSH-01\Rain\raintemp\rain.gq.map
undef 32766
title 8phases rain
tdef 8 linear 00z01jan0000 1dy
vars 1
r 0 99 rain (0.1mm)
endvars
;
！stnmap 已经成功制作了.map文件！
❤----------------❤--------------------------------❤--------------------------------❤
（3）制作网格文件FOR：
!---------------------------------------------------------
!This program is used to build a grid file of binary.
!------------------------------------------------------
program main
parameter(nx=71,ny=41,nt=8)
real g(nx,ny,nt),lat(ny),lon(nx)
open(1,file='D:\STUDY\WPSH-01\Rain\raintemp\grid.grd',form='binary')
lat(1)=15.0
lon(1)=70.0
do j=1,ny-1
lat(j+1)=lat(j)+1.0
enddo
do i=1,nx-1
lon(i+1)=lon(i)+1.0
enddo
do t=1,nt
do i=1,nx
do j=1,ny
g(i,j,t)=1
enddo
enddo
enddo
do t=1,nt
write(1)((g(i,j,t),i=1,nx),j=1,ny)
enddo
close(1)
end program
（4）网格资料的描述文件：
dset D:\STUDY\WPSH-01\Rain\raintemp\grid.grd
undef -999.0
title sample grid data
xdef 71 linear 70 1
ydef 41 linear 15 1
zdef 1 linear 1000 1
tdef 8 linear 00z01jan0000 1dy
vars 1
g 0 99 grid data
endvars;
❤----------------❤--------------------------------❤--------------------------------❤
(5) 画图.gs文件
'reinit'
'open D:\STUDY\WPSH-01\Rain\raintemp\grid.ctl'
'open D:\STUDY\WPSH-01\Rain\raintemp\rain.gq.ctl'
'enable print D:\STUDY\WPSH-01\Rain\raintemp\rain.gq.gmf'
'set mpdset hires'
i=1
while(i<=8)
'set t 'i
'set lon 70 140'
'set lat 15 55'
'define a=oacres(g,r.2,1.5)'
'define a1=maskout(a,g-0.5)'
'define aa=smth9(a1)'
'set grads off'
'set grid off'
'set gxout contour'
'set clab forced'
'set cint 0.1'
'd aa'
'draw title phases-'i
'print'
'c'
i=i+1
endwhile
'disable print'
;
（6）画出来的图，数值都是0……
❤----------------❤--------------------------------❤--------------------------------❤
重复下问题：是不是我程序哪里出错了，没有插值对，导致数据都为零了？
另外就是，处理台站降水资料一定要把缺测台站的所有值都除掉吗？做了缺测处理也不行？
|
|
窥视卡
雷达卡
发表于 2013-4-17 11:50:49
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
发表于 2013-4-17 12:05:32
发表于 2013-4-17 12:06:17