求助：grads画站点

Malik

最近看了站内的教程，然后找了160站数据来画图，参考站内的程序自己写了一个
program sta2grd
implicit none
integer i,j,k,nlev,flag
integer,parameter::sta=160,year=60,mon=12
real tim,lon(sta),lat(sta)
integer::stid(sta),rain(sta,year,mon),temp(sta,year,mon)
character filename1*10,filename2*10,cmon*2

open(1,file='sta.grd',status='replace',form='binary')
open(2,file='station.txt',status='old')
do i=1,sta
read(2,*)stid(i),lat(i),lon(i)
end do
close(2)

do k=1,mon
if(k<=9)then   
  write(cmon,'(i1)')k   
  cmon="0"//trim(cmon)   
else
    write(cmon,'(i2)')k  
endif
filename1="r16"//trim(cmon)//".txt"
filename2="t16"//trim(cmon)//".txt"
open(2,file=filename1,status='old')
open(3,file=filename2,status='old')
read(2,*)((rain(i,j,k),i=1,sta),j=1,year)
read(3,*)((temp(i,j,k),i=1,sta),j=1,year)
close(2)
close(3)
end do



do j=1,year
do k=1,12
do i=1,sta
tim=0.0
nlev=1
flag=1
write(1)stid(i),lat(i),lon(i),tim,nlev,flag,rain(i,j,k),temp(i,j,k)
end do
nlev=0
write(1)stid(1),lat(1),lon(1),tim,nlev,flag
end do
end do
close(1)

end

降水和温度写一起的，我估计程序最后写入的部分有问题，然后生成map文件出错
这是ctl
DSET D:\160zuotu\sta.dat
DTYPE station
STNMAP D:\160zuotu\sta.map
UNDEF -32766
TITLE 160 Monthly rain temp data
TDEF 720 linear jan1951 1mo
VARS 1
rain 0 99 rain data
temp 0 99 Surface Temperature
ENDVARS
不知道是哪里不对啊

冬日情愫+

看不出哪里错了。楼主说的错误是在哪一步出现的呢？不过提几点建议吧，建议自己读写文件的时候不要使用10以前的文件号，因为10以前的文件号fortran内部使用。另外源代码中
if(k<=9)then   
  write(cmon,'(i1)')k   
  cmon="0"//trim(cmon) 有点多余只需要一个语句就可以完成：write(cmon,'(i2.2)')k。

Malik

冬日情愫+ 发表于 2013-8-30 17:57

                               
登录/注册后可看大图

看不出哪里错了。楼主说的错误是在哪一步出现的呢？不过提几点建议吧，建议自己读写文件的时候不要使用10以 ...

这一段是复制的，字符的运算不太熟悉，我又参考了几个帖子，已经把程序和CTL都改了一下，画出来站点上没有数据，都是0，不知道怎么回事

冬日情愫+

Malik 发表于 2013-8-30 18:04

                               
登录/注册后可看大图

这一段是复制的，字符的运算不太熟悉，我又参考了几个帖子，已经把程序和CTL都改了一下，画出来站点上没有 ...

既然有数，说明应该是你的源代码有问题。integer::stid(sta),rain(sta,year,mon),temp(sta,year,mon)
你这一行的声明应该有问题吧，降水和温度应该都是实数吧，还有站点号一般都是5位，应该超出了32767了吧，建议用字符型来表示吧，省得出错。

nuister

integer::stid(sta),rain(sta,year,mon),temp(sta,year,mon)这句错了。grads里面stid必须是字符串型，而且必须是8位的。还有就是你这温度和降水咋都是整型呢？建议改成：
character*8 stid(sta)
real rain(sta,year,mon),temp(sta,year,mon)

Malik

感谢楼上两位啊！
站号改成了字符
降水和温度改成了实数型
倒数第五行改成了另一个帖子的write(1)stid(i-1),lat(i-1),lon(i-1),tim,nlev,flag
CTL第一行最后应该是.grd
CTL少了一个层次
站点上数据已经画出来了

leeol

纠结了一天的问题终于解决了，谢谢5L

269772539

谢谢楼主的分享！