多变量单时次观测数据grads绘图求助

liuzf406

各位朋友, 我有清风版主给的代码和数据,可以绘出图来, 但我用自己的数据给图失败. 我的数据格式如下:
Precipitation isotope data  July 2012         
station Lon Lat Alt d18O d2H Dex Annprecip Anntemp Dd18O
50353   117.10  39.06  3.00  -7.70  -49.97  11.64  496.20  14.48  -1.06
50632   121.48  31.22  4.00  -8.05  -50.81  13.61  1111.00  15.79  -4.10
50527   113.32  23.13  7.00  -5.90  -41.01  6.20  1723.50  22.18  -2.96
50434   110.21  20.02  15.00  -5.95  -39.55  8.10  1655.39  24.92  -3.75
50557   119.17  26.05  16.00  -6.63  -42.90  9.76  1223.63  19.98  -2.11
..........

第一步读数据没问题, 第二步按grads格式输出也没问题, 第三步生成站点数据的ctl文件也没问题, 第四步可以看到数据.
但在第五步出错了, 我想应该是插值生成格点时出错了, 我就直接用了清风大人给的代码:
Program gen_grid
Implicit none
Integer,parameter::x=33,y=17
Integer i,j
Real grid
grid=1.0
Open(1,file='grid.grd',status='replace',form='binary')
Do i=1,y
  Do j=1,x
    Write(1)grid
  Enddo
Enddo
Close(1)
End
但这里看不出如何输出Alt、 d18O、 d2H 、Dex 、Annprecip 、Anntemp 、Dd18O等变量的格点啊。 结果就出现如下错误：

mofangbao

问题就出在你以为没问题的地方

做个霸气的木头

你grid.ctl文件出错了，贴上来看看，貌似少了什么~

liuzf406

谢谢楼上两位啊, 我再检查下下,把东西贴上来求救!

liuzf406

贴出来我做的过程,请大家帮我诊断下怎么错了.
第1步,数据格式如下:
Precipitation isotope data  July 2012                                                                        
station        Lon        Lat        Alt        d18O        d2H        Dex        Annp        Annt        Dd18O
50353   117.10         39.06         3.00         -7.70         -49.97         11.64         496.20         14.48         -1.06
50632   121.48         31.22         4.00         -8.05         -50.81         13.61         1111.00         15.79         -4.10
.........................................

数据有些保密要求, 所以不能全部显示!

谢谢

liuzf406

第2步,用fortran读这些数据, 代码如下(似乎没有错,但我的数据中有缺失值,  设为-32766.00,对于其,代码里面没有体现出来,这里有必要处理吗?如何处理? ):
program sta2grd
implicit none
Character*8 stid
Real lon,lat,alt,d18o,d2h,dex,annp,annt,dd18o,tim
Integer nlev,flag
tim=0.0
nlev=1
flag=1
Open(1,file='isotope.txt',status='old')
Open(2,file='sta.grd',status='replace',form='binary')
Read(1,*)
Read(1,*)
10 Read(1,*,end=50)stid,lon,lat,alt,d18o,d2h,dex,annp,annt,dd18o
Write(2)stid,lat,lon,tim,nlev,flag,alt,d18o,d2h,dex,annp,annt,dd18o
Goto 10
50 continue
close(1)
nlev=0
Write(2)stid,lat,lon,tim,nlev,flag
Close(2)
End

liuzf406

第3步按照 GrADS的规定将数据进行输出, ctl 文件如下:
DSET C:/isotope_grads/sta.grd
DTYPE station
STNMAP C:/isotope_grads/sta.map
UNDEF -32766.00
TITLE Isotope Data
TDEF 1 linear Jul2012 1mo
VARS 7
alt   0 99 elevation of the station
d18o  0 99 stable oxygen isotope value
d2h   0 99 stable hydrogen isotope value
dex   0 99 deuterium of excess
annp  0 99 annual precepitation amount
annt  0 99 mean annula temperature
dd18o 0 99 difference in d18o between summer and winter
ENDVARS
可以输出, 见下图:

liuzf406

第4步, 查看各变量,显示各变量, 似乎没问题, 见下图:

如d18o, dex等变更均可以显示, 见下图:

liuzf406

第5步,生成格点的背景场，为插值做准备, 这一步出现了问题, 由于有7个变量,我就加了个k, 似乎也没有错误, 见下面代码:
Program gen_grid
Implicit none
Integer i,j,k
Integer,parameter::x=33,y=17
Real grid
grid=1.0
Open(1,file='grid.grd',status='replace',form='binary')
Do k=1,7
Do i=1,y
  Do j=1,x
    Write(1)grid
  Enddo
Enddo
Enddo
Close(1)
End

liuzf406

第6步,为生成的格点场配置CTL文件, 如下:
DSET C:/isotope_grads/grid.grd
TITLE Gridded Data
UNDEF -32766.00  
XDEF 33 LINEAR 70.0 2.5
YDEF 17 LINEAR 15.0 2.5
ZDEF 1 LEVELS 1000
TDEF 1 LINEAR Jul2012 1mo
VARS 7
alt   0 99 elevation of the station
d18o  0 99 stable oxygen isotope value
d2h   0 99 stable hydrogen isotope value
dex   0 99 deuterium of excess
annp  0 99 annual precepitation amount
annt  0 99 mean annula temperature
dd18o 0 99 difference in d18o between summer and winter
ENDVARS
可以在grads中打开此ctl文件,且可以显示每个变量(如alt, d18o,d2h..., 但值格点值均为1),似乎也没有错.见下图: