下载文档原格式
以下技巧总结都是笔者从学习实践过程中总结出来的，基本的问题。
不求全面，希望对读者学习有用，如果有问题，敬请留言指正，以促进交流学习！1、软件综述：grads软件是一款绘图软件除了绘制图形，还可以提取数据，主要应用是在大气科学中，当然只要是数据处理成grads能够读取的数据文件就可以进行相关绘图。
软件版本问题，软件本身不是很大，我接触到1.8、1.9、2.0版本的，1.8版本的安装很多情况还要修改环境变量、1.9版本的不识别‘sdfopen’命令，最稳定的版本是2.0版本，所以笔者推荐学习者安装2.0版本，选择默认安装路径就可以。
2、文件类型简述：grads处理的是网格数据，可以处理的数据类型有：grd、grib、nc（海洋常用的数据），cdf（雷达卫星数据），其中nc、cdf数据都是自带描述文件，不需要ctl，grib数据要通过命令生成ctl、index数据才可以调用，常用的是grd数据，需要ctl。
3、数据文件转换：grads软件识别的数据是二进制无格式数据，文件类型是‘binary’，写入和生成时是不需要格式的如read(20) sst(i,j,iz,it),20为文件号，通常是十进制数据与grd数据间转换，这里给一个grd转换成txt数据的fortran程序：parameter(nx=56,ny=41,nz=1,nt=360)dimension sst(nx,ny,nz,nt)real sstopen(15,file='sst.grd',form='binary') ！固定的用form=‘binary’就是二进制数据open(16,file='sst.txt') ！新建txt文件do it=1,ntdo iz=1,nzread(15) ((sst(i,j,iz,it),i=1,nx),j=1,ny) ！read后只有文件号，数据是无格式的enddoenddodo it=1,ntdo iz=1,nzwrite(16,*) ((sst(i,j,iz,it),i=1,nx),j=1,ny) ！输出时是txt文件可直接看的数据，有格式输出，有*enddoenddoclose(15)close(16)end写程序时：注意格点数要与数据对应，如：上程序对应的数据是经度90~200，纬度-20~60，时间：1971.01~2000.12共360个月的海面温度数据，数据格点精度2*2 ，nx=(200-90)/2+1,ny=(60-(-20))/2+1,nt=360，nz=1，大气的数据要根据数据的层次确定几层。
GrADS绘图软件实用手册2002年1月目录第一章 GrADS绘图软件概述1.GrADS绘图软件简介2.Internet上的GrADS资源2.1GrADS在Internet上的主页2.2 windows环境下GrADS资源3.GrADS绘图软件的安装（windows环境）3.1在windows环境下安装GrADS软件包3.2X server 的安装第二章 GrADS绘图模板1.GrADS示例演示1.1 启动GrADS1.2 退出GrADS1.3 示例演示GrADS命令的使用2.GrADS绘图模板3.GrADS模板的高级应用3.1GrADS描述语言3.2GrADS高级模板的应用第三章 GrADS数据格式1.格点数据描述文件1.1 数据描述文件各项解释1.2 生成model.le.dat和model.le.ctl文件的程序代码片段2.站点数据的格式附录1.如何精确控制图形输出的尺寸—Landscape纸型2.台站资料的显示3.Linux环境下的安装第一章 GrADS绘图软件概述1GrADS绘图软件简介The Grid Analysis and Display System(GrADS) 是一套应用广泛、使用方便的科学数据绘图软件包。
其主要特点：z GrADS属于自由软件，可以从Internet上免费获得。
z可运行于各种Windows 和Unix工作平台。
z GrADS可用于4D数据的分析。
既经度、纬度、层（气压层、高度层等）和时间/xyzt 4维。
数据可以是格点化的数据或离散点数据。
GrADS特别适用于气象类数据的分析。
但也完全可以用于更广泛类型的数据分析。
z GrADS有多种显示方式：等值线、流线、矢量图、风矢量图、站点填图、折线图、直方图等多种两维图形。
z可处理多种数据格式的数据。
GRIB、NetCDF、HDF-SDS等通用数据格式和系统自定义的一种二进制数据格式。
z采用命令行输入的方式交互式地显示图形。
Grads站点资料的多层次绘图看了很多关于如何使用grads和fortran转换站点资料的帖子，大多数都集中在地面或者单层的，对于高空的多个层次很少有提到，而且大多数有些模糊，不好理解，我把自己做的一些整理传上来，供大家参考参考。
欢迎大家指出理解和程序的错误。
一，首先是数据的说明选用micaps高空观测资料850hpa,700hpa,600 hpa,500hpa，200hpa，100 hpa 六个层次作为例子。
为了方便，时间层次就选择1个时间。
micaps高空观测资料的说明如下：diamond 2 10年07月28日20时200百帕高空观测10 07 28 20 200 419因此定义mt=419,nz=4数据每行的排列分别是：站号，经度，纬度，高度……..具体详见micaps 的使用手册，这里不再赘述。
根据这些定义变量。
二，fortran读取站点资料注：1.根据micaps的高空资料定义相应的变量，由于850-100这6个层次的站点是相同的，站号，经纬度等只需定义成一维的，如lat(mt)，而lev(nz)表示有6个层次，也为一维；高度场等又有mt个站的资料又有nz个层次，因此定义为2维的。
2.由于micaps资料分别在不同的文件内，文件名重复，我修改为100_2820,200_2820,500_2820…….然后使用trim函数顺序读取资料。
3.为了更好的说明地面和高空资料的读取，将地形高度high(mt)设置为地面的量，nlev定义为7，即高空六层，地面一层4.以输出地面一个变量，高空三个变量为例。
先进行站点的循环，再进行高空层次的循环。
5.前面我们定义了lev(nz)，现在给lev(nz)赋值，使用select case.,分别赋值为850,700,600,500,200,100。
在使用是发现如果不进行赋值，后续工作也无法进行，会出现错误。
程序如下：program stationparameter (mt=419,nz=6) ! mt 总站数，nz总层次character stid(mt)*8integer flag,nlevreal lat(mt),lon(mt),tim,high(mt),number(mt),lev(nz) realheight(mt,nz),temp(mt,nz),ttd(mt,nz),ud(mt,nz),us(mt,nz),u(mt,nz),v(m t,nz),td(mt,nz)！顺序读取资料character stationfile(6)*8data stationfile/'850','700','600','500','200','100'/do k=1,nzopen(k,file=''//trim(stationfile(k))//'_2820.000')read(k,*)read(k,*)！跳过micaps开头两行的说明文字do i=1,mtread(k,*)stid(i),lon(i),lat(i),high(i),number(i),height(i,k),temp(i,k ),ttd(i,k),ud(i,k),us(i,k)end doend do !结束k的循环tim=0.0flag=1nlev=7 ！高空+地面open(40,file='20.grd',form='binary')do i=1,mt ！站点循环write(40)stid(i),lat(i),lon(i),tim,nlev,flag,high(i)do k=1,nz ！高空循环！对lev(nz)进行赋值select case(k)case(1)lev(k)=850.0case(2)lev(k)=700.0case(3)lev(k)=600.0case(4)lev(k)=500.0case(5)lev(k)=200.0case defaultlev(k)=100.0end select！在lev(k)后将要输出的所有高空变量一次性输出write(40)lev(k),height(i,k),temp(i,k),ttd(i,k)end doend donlev=0write(40)stid(mt-1),lat(mt-1),lon(mt-1),tim,nlev,flag close(40)end三，ctl及map文件的生成注：看到相关帖子说在站点资料的ctl文件中也要设置zdef，我试了下，设置和不设置都是可以的Ctl文件如下：dset f:\20.grddtype stationstnmap f:\20.mapundef 9999.0title 20zdef 6 levels 850 700 600 500 200 100tdef 1 linear 28jul2010 1DYvars 4high 0 99 Terrain Heighth 6 99 heightt 6 99 temperaturettd 6 99 ttdendvarsCtl文件编写完成后，就可以在grads中生成map文件了，相关命令如下：！ stnmap -i f:/20.ctl注意，斜杠一定是“/”,方向不能反第三步结束后，站点资料的读取就结束了，可以通过grads打开Ctl文件看到相关站点的值了。
利用Grads画站点图（contour、shaded、grfill等）1.利用Fortran程序将数据输出为grd格式。
参考程序如下：parameter(num=160) （站点数）character*8 sta(num) (站名名数组，可任意）dimension xlon(num),ylat(num),rc(num) (经、纬度数组及其上对应的数值) open(30,file='cor.grd',form='unformatted') (工作站)open(30,file='cor.grd',form='binary') (微机)do 100 lev=1,20 (共输出20个时间上的观测或结果)tim=0.0ilev=1nflag=1do 10 i=1,numsta(i)=char(i)write(30) sta(i),ylat(i),xlon(i),tim,ilev,nflagwrite(30) rc(i)10 continueilev=0write(30)sta(num),ylat(num),xlon(num),tim,ilev,nflag (每个时次的结束) 100continueclose(30)2.ctl文件(创建与步骤1中输出的站点数据配对的station.ctl文件)：dset cor.grdformat sequentialdtype stationstnmap cor.mapundef -9.99e33tdef 20 linear jun1958 1movars 1r 0 99 correlationendvars3.为插值函数准备格点数据grid.grd及对应的grid.ctldset grid.grdundef -99.0xdef 71 linear 70 1ydef 41 linear 15 1zdef 1 linear 1 1tdef 20 linear jun1958 1movars 1g 0 99 grid data prepared for oacres functionendvars注意：两个ctl（station.ct l和grid.ctl）文件中的时间要严格一致。
![[转]grads中站点处理和oacresmaskout函数应用心得](https://uimg.taocdn.com/43f4826100f69e3143323968011ca300a7c3f655.webp)
[转]grads中站点处理和oacresmaskout函数应⽤⼼得SJ 发表于: 2011-1-05 13:22 来源:昨⼉刚刚帮同事画了两张降⽔图（不同时间段的），很是费劲，中间碰到若⼲问题，但好在最后基本达到了满意的效果，下⾯把这次画图中遇到的问题及解决办法贴出来，供⼤家参考或指正：1，数据处理的问题由于micaps的24⼩时降⽔资料中，并不是每个时次的⽂件⾥的站点个数都是相同的，就更不⽤提两个⽂件⾥有完全相同的站点了，这样就给各个站点降⽔的累加处理造成了极⼤的困难（如果⼀开始就有⼈或专门的机构来做这件事，将⼤⼤提⾼科研⼈员的⼯作效率），但本帖不打算过多描述数据处理的情形，因此这个问题就到此为⽌。
2，站点数据中打⽹格和maskout的联系降⽔的站点资料处理完毕并转成grads能识别的⼆进制格式后，就需要给站点数据打⽹格了，这个⽹格是⽤于站点资料差值使⽤的（如果⽤grads，就要⽤格点资料），此时就涉及到要打⼀个什么分辨率的⽹格了，在这次画图过程中，我⼀共打了两个全球⽹格，⼀个是0.1X0.1（细⽹格）的另⼀个是0.5X0.5（粗⽹格）的。
然后按部就班的使⽤了如下的语句：'set mpdset china' * 设定⼀个地图，⽤于只画中国范围内的物理量'define a=oacres(g,tt.2,XXXX)' * 插值函数，把变量tt插值到⽹格g上，其中XXXX代表影响半径'define a1=maskout(a,a-8)' * maskout（n1,n2）⾥第⼀个变量n1是我们要画的物理量，第⼆个变量n2的意思是只画当其值⼤于零时的物理量n1，* 当然n1和n2是⼀⼀对应的，⽽这个例⼦中的意思是，只画物理量a的值⼤于4的部分'define aa=smth9(a1)' * 平滑函数'run province-basemap china aa' * 只画中国范围内物理量下⾯通过格点的分辨率-2种，和oacres中的影响半径做组合实验：（1）粗⽹格+XXXX（50,30,20,10,5,1）见图⼀（2）细⽹格+XXXX（20,10,5,1）见图⼆（3）细⽹格+XXXX（50,30,20,10）见图三（4）粗⽹格+XXXX（50,30,20）50是grads中默认最⼤的数值见图四从这四张图中，可以看到⽐较⼤的差别，这也是我画图时思考的⼀个顺序。
由于很多同志是做业务的，接触的都是站点资料，这里将详细的介绍如何利用站点资料，画出某变量的等值线图。
第一步：站点资料（txt格式，即十进制格式）转成二进制格式GrADS画等值线图是有要求的，即只能画格点上的二进制数据。
首先，将观测资料按如下格式放在txt文件中：第一列为某站点的纬度，第二列为该站点的经度，第三列为该站点某月平均某一个指数，也就是变量，它可以是降水、气温等。
然后，将该数据转换为二进制，使用以下fortran程序：c----------------------------------------------------------c This program is used to change the format of ascii datac into the form of binary ,which is supported by the GrADSc----------------------------------------------------------program mainreal, dimension(160) :: lat, lon, indopen(1,file='d:\drought\index.txt',status='old')do i=1,160read(1,*) lat(i), lon(i),ind(i)enddoclose(1)call stntogrd(ind)endsubroutine stntogrd(x)character*8 stid(160)do i=1,160stid(i)=char(i)enddoopen(3,file='e:\drought\index.grd',form='binary')tim=0.0nlev=1nflag=1do i=1,160write(3) stid(i),lat(i),lon(i),tim,nlev,nflag,x(i)enddonlev=0write(3) stid(i-1),lat(i-1),lon(i-1),tim,nlev,nflagclose(3)returnend于是得到了二进制格式的变量index.grd,接着我们需要给该二进制数据写一个描述性文件如下（新建一个写字板，取名为"index.ctl"）：dset d:/drought/index.grddtype stationstnmap d:/drought/index.mapundef -999.0title drought indextdef 1 linear Jul1951 1movars 1ind 0 99 drought indexendvars到此，二进制数据文件的描述文件写完了，然后，在GrADS命令窗口输入命令“！stnmap”如：ga_> !stnmap系统会提示如：“Enter stn ctl filename:” 输入: d:/drought/index.ctl 系统会生成"index.map"文件。
GrADS绘图软件使用手3第三章GrADS数据格式每一组GrADS数据应至少包括两组数据文件，数据描述文件—ASCII 码和数据文件—二进制，数据的真正存放地。
数据文件中只是用户数据的有序排放，而关于数据种类、排放次序等是单独放在一个文件中的称—数据描述文件。
而象GRIB和NETCDF等通用数据格式，以上两者是存于同一个文件的—或称为自定义/自解释格式数据。
但考虑到GrADS传统，对这类自定义格式数据仍将生成相应的数据描述文件。
上一章中我们已使用过了这样的一组数据。
以此为例，介绍用户如何按GrADS的格式，将自己的数据生成相应的数据文件和数据描述文件。
1.格点数据描述文件model.le.ctl文件清单：以某开始的行为注解行。
1.1数据描述文件各项解释1.DSET数据文件名定义与此数据描述文件相对应的数据文件名。
若两者位于同一目录，前面的路经可以省略或以“^”开始，代表两者位于同一目录。
若不在同一目录下，应给出路经参数。
如：c:/pcgrad/ample/model.le.dat56注意路经的给法与DOS不同，而与UNI某环境一致，便于移植！或c:\\pcgrad\\ample\\model.le.dat,两种都行。
2.TITLE数据文件说明文字串。
3.UNDEFvaule定义缺测值。
一般给一很大的正/负值，表示，当取值超过这一正值/低于定义的负值，认为该值无效。
（GrADS采用跳过或用周围有效点的值处理。
）4.OPTIONS这里定义了与二进制存储有关的选项，二进制存储的一大特点是可移植性差，因此通过keyword项来增加可移植性。
若keyword省略，则OPTIONS也可省略。
可取：equential:顺序无格式方式。
yrev:Y维与YDEF定义相反方式存放。
zrev:Z维与ZDEF定义相反方式存放。
big_endian:如数据是在un,gi,hpcray机器上生成的，而目前不在此类机器上使用。
用grads画站点资料的等值线分类：专业绘图 | 时间：2010-01-07 10:06 | 阅读：24人/次 | 发布者：laiwf一先把站点文件转化为二进制文件假定站点资料为station lat lon var….…. ….…..1 单个时间变量，单个要素bintd.f90program bintdparameter (st=3276) ! st means nums of stationcharacter stid(st)*8integer nflag,nlevreal lat(st),lon(st),std(st),timopen(30,file='02070108.txt')open(40,file='td.dat',form='binary')do i=1,stread(30,'(a5,2f8.2,f7.1)')stid(i),lat(i),lon(i),std(i)write(*,*)stid(i),lat(i),lon(i),std(i)end doclose(30)tim=0.0nflag=1nlev=1do i=1,stwrite(40)stid(i),lat(i),lon(i),tim,nlev,nflag,std(i)end donlev=0write(40)stid(3276),lat(3276),lon(3276),tim,nlev,nflag close(40)end2 多个时间变量，单个要素binstat.f90假定站点资料为time station lat lon var……. …………….program binstatcharacter stid*8integer nflag,nlevreal lat,lon,std,timopen(30,file='st020120.txt')open(40,file='bist020120.dat',form='binary')iflag=0i=010 read(30,'(i3,1x,a5,2f8.2,f7.1)',end=90)id,stid,lat,lon,varif(iflag.eq.0)theniflag=1ims=imids=idend ifif(ids.ne.id)thennlev=0write(40)stid,lat,lon,tim,nlev,nflagi=i+1write(*,*)iend ifims=imids=idtim=0.0nflag=1nlev=1write(40)stid,lat,lon,tim,nlev,nflagwrite(40)vargo to 1090 continuenlev=0write(40)stid,lat,lon,tim,nlev,nflagstopend二创建站点文件的描述文件1单个时间变量单个要素 td.ctldset td.datdtype stationstnmap td.mapundef 999.9title daily tdtdef 1 linear Jan2002 1movars 1td 0 99 tdpoint dataendvars2 多个时间变量单个要素 stat.ctldset bist020120.datdtype stationstnmap svar020120.mapundef 999.9title daily svartdef 31 linear Jan2002 1movars 1svar 0 99 svars dataendvars三stnmap –i xxx.ctl四生成相应的格点文件 wgrid.f90program wgrid! to get gridsreal gtd(50,50)! redefine the grids numbers when necessaryopen(10,file='/Models/czh/wgrid.dat',form='binary') do i=1,50do j=1,50gtd(i,j)=0end doend dowrite(10)((gtd(i,j),i=1,50),j=1,50)end五生成格点文件的描述文件grid.ctldset wgrid.dattitle gtdundef -9.99E33xdef 50 linear 115.07 0.18ydef 50 linear 39.29 0.18zdef 1 linear 1 1tdef 1 linear jan2002 1movars 1gtd 0 99 grid tdendvars六生成最后的gs文件stat.gs1 显示一张图'reinit''open gtd.ctl''open td.ctl''set mpdset cnworld''set lon 102.2 125.8''set lat 31.4 46.4''set t 1''define tds=oacres(gtd,std.2)''set cint 0.5''set gxout shaded''d tds''set gxout contour''d tds''printim td.gif white';2 同时显示多个时刻的图'reinit''open gtd.ctl''open 200201.ctl''set mpdset cnworld''set lon 115.17 117.672''set lat 39.29 41.216'tm=1while(tm<=31)'set t '%tm'define statvar=oacres(gtd,svar.2)' 'set cint 0.5''set gxout shaded''d statvar''set gxout contour''d statvar''printim statvar'% tm % '.gif white' 'd statvar''c'tm=tm+1endwhile。
Fortran+GrADS站点作图详解mofangbao@气象家园帖子导读：以前论坛也有过相关的帖子，但是到目前为止，仍然不断的有朋友看了那些帖子之后还是不知道如何下手，于是我打算为纯粹的新手写一份说明，这份文档会非常详细，因此适合刚接触GrADS的朋友来参考。
帖子介绍了fortran文件中读写数据的基本知识，grads的站点文件的存放格式，如何通过fortran来写出这个格式的文件，如何建立站点文件的CTL文件，如何为站点数据生成映射文件，站点映射文件的作用，如何通过站点CTL文件来预览站点数据，如何用fortran生成grads支持的格点文件，如何给格点配置CTL文件，如何编写GS文件来进行插值，如何绘制等值线以及阴影图，如何屏蔽区域外的图形等内容。
如果你对以上问题存在疑惑，希望这个帖子能够对你有所帮助，因此，看完这份文档之后，您应该可以做到这些：（会的就跳过吧，我也是犹豫了很久要不要写的，算是新手礼包的一部分吧）开始之前，先了解站点数据的基本内容一般情况下，我们看到的站点数据是这样的：50353 126.39 51.43 100 8950632 121.55 48.46 100 5450527 119.45 49.13 100 2950434 121.41 50.29 100 4650557 125.14 49.10 100 8350745 123.55 47.23 100 8050756 126.58 47.26 100 7450788 131.59 47.14 100 3650873 130.17 46.49 100 3150978 130.57 45.17 100 59第一列为站号，第二、三列为经纬度，第四、五列为其他数据。
当然，也可能是把经纬度单独存放，然后数据单独存放，比如气候中心免费开放的160站的月平均温度降水值。
如果是micaps资料可以查看micaps用户手册。
这个例子中使用的就是micaps的第三类数据格式，示例文件可以从这里下载：（资料下载）除了了解这些，还必须了解该数据的缺测值用什么来表示，这个在后面的数据描述文件中要用到。
1、在fortran中读入这些资料在使用fortran之前，请确保你已经安装了fortran编译器，如果你是xp系统，建议安装CFV6，如果是win7系统，那么可以安装microsoft fortran powerstation4.0（以下简称4.0），这两个软件在气象资料站均有下载，安装完成后，请在你磁盘的某个文件夹下面新建一个文件夹用于本次作图，如果使用的是CVF6，该文件夹的路径中不要包含中文名。
新建文件夹完成后，从开始菜单打开fortran编程窗口，然后点击file->new，新建一个free format的自由格式fortran文件,4.0的直接建立一个TextFile即可：CVF中新建fortran文件4.0中新建TextFile即可注意：CVF中先选路径和文件名，然后会建立该文件，4.0中直接建立一个临时文件，当你按保存按钮时才会提示你选择保存的路径，如果你之前用4.0的fortran 编译过其他文件，那么请先关闭工作区再新建，从file->close workspace可以关闭。
当然也可以直接建立一个文本文件，然后重命名为sta2grd.f90，记住，要修改后缀为.f90，表示该文件为自由格式的fortran源文件。
现在就可以开始着手读取数据了。
（1）先把fortran的开头和结尾写了，并且用implicit none声明不能使用未定义的变量，同时添加一个注释，自由格式中可以用一个英文的!来表示该行该符号往后均为注释内容，不参与编译：Program sta2grdImplicit none!这里是程序的变量声明!变量声明结束!程序开始!程序结束End（2）定义几个变量Fortran中常用的变量类型有整型：integer，单精度实型real，双精度实型double，字符型character，关于fortran变量的更多内容请参考彭国伦的fortran95程序设计。
通过观察我们要读取的文件：diamond 3 2011年05月降水量值0 0 0 0 0 0 1 0 0 1 100050353 126.39 51.43 100 8950632 121.55 48.46 100 54…第一行和第二行对我们转换没有什么作用，因此我们只需要定义必须要用到的变量，站号stid，类型为字符型，长度为8，其余的经度lon、纬度lat，一个无意义的变量（列中的100）var和降水数据rain，均为实型，原因在下面的GrADS中的存放方式有说明。
Character*8 stidReal lon,lat,var, rain（3）打开数据文件并读入数据Fortran中使用open语句来实现打开文件，同时指定一个打开的文件号和文件路径，以及其他一些指定参数，下面的语句即可打开这次示例的文件：Open(1,file=’raindata.txt’,status=’old’)这里的1是一个文件号，是读写时需要用到的，file后的是文件名，直接写文件名表示在和程序同一个文件夹，status表示当前文件的状态，old表示该文件已经存在，这样就能避免通过write语句误写入数据覆盖了数据文件。
（4）开始读入数据到变量中在fortran中使用read语句来读入数据，每读取一行文件的指针就会自动向下移动一行（关于文件指针指向这一行你可以想象成该行被选中，要从该行读取数据了）比如read(1,*)var，就表示从打开的1号文件中按默认格式读取数据到var中。
我们先跳过两行“没用”的数据：Read(1,*)Read(1,*)这样后面不接任何变量就表示把这两行跳过去了，现在文件指针在第三行的开头了，也就是我们需要的数据了，来读一行：Read(1,*)stid,lon,lat,var,rainFortran 会自动把这五个数值分别给后面的五个变量了，现在文件的指针已经在第四行的开头了，我们需要继续读取数据，当然不能够又写一行，应该继续使用上面这一行来读取，可以通过goto语句实现，goto语句的含义是让程序跳转到指定的标号执行，这里就应该是去刚才的那个语句执行，所以在刚才的Read(1,*)stid,lon,lat,var,rain 前面加上一个标号10，然后在下面写上一句goto 10，如下所示：10 Read(1,*)stid,lon,lat,var,rainGoto 10这样，程序就会不断的往下读取了，但是，文件的行数是有限的，这样做一旦到文件结束，fortran就会报出运行时错误，通常是end-of-file…的错误，这时候可以通过END标签来避免这种情况的出现。
End标签能够让read语句读取到文件末尾的时候自动跳转到指定的标号去执行，再次修改这个读取的语句：10 Read(1,*,end=100)stid,lon,lat,var,rainGoto 10100 contine这样，当文件读取结束时，就会到100 continue来执行了，continue的意思就是接着按顺序往下执行，本身没有什么含义。
我们现在只是读取数据，没有对数据进行任何输出和处理，读取完毕后应该及时关闭这个文件，使用close方法即可。
在100 continue的下面加上close(1)，1同样是open里的那个文件号。
这段程序就完成了这个示例文件的读取，每读取一次，上面几个变量的值就会被刷新一次，如果需要验证读取是否正确，可以在标号为10的read语句和goto 10之前加上一个输出到屏幕的命令，然后用pause命令暂停执行，以便检查：10 Read(1,*,end=100)stid,lon,lat,var,rainPrint*,stid,lon,lat,var,rainpauseGoto 10现在，这段程序应该是这样子的：Program sta2grdImplicit none!这里是程序的变量声明Character*8 stidReal lon,lat,var, rain!变量声明结束!程序开始Open(1,file=’raindata.txt’,status=’old’)read(1,*)read(1,*)10 Read(1,*,end=100)stid,lon,lat,var,rainPrint*,stid,lon,lat,var,rainpauseGoto 10100 continueClose(1)!程序结束End这时候就可以编译运行一下了，4.0中点击窗口中左上角的build按钮（shift+f8），弹出的对话框中选择是，然后下面的消息窗口会提示:sta2grd.exe - 0 error(s), 0 warning(s)，表示编译通过了，可以开始执行了，选择build菜单的Execute sta2grd.exe，就可以开始执行了；CVF中直接点击红色的感叹号左侧有两个朝下箭头的按钮编译，无错后点击红色感叹号可以执行，每行输出后都会暂停，等按下回车后会继续读取下一行并暂停，经过对比，发现读取成功，然后删除用于测试的print*和pause语句：4.0中的编译2、按照GrADS的规定将数据进行输出接下来所增加的语句都在10号语句和goto 10语句之间，因为这时候几个变量的值是刚刚从文件中读出来，所以要及时把他们写出去。
在输出为GrADS可用格式之前，必须先了解GrADS中的站点数据是如何存放的，只有知道GrADS中规定的存放格式，才能够明白输出时该如何写。
GrADS中的站点数据也是二进制格式，并且按照时间来分组，（本例子中仅有一个时次，多时次的相信你看完后也会写出来了），每个时次下面包括了N个站该时次的数据，每个站的信息由一个表示该站的文件头信息以及数据值组成，每个时次结束后写入一个特殊的标记表示这个时次结束，文件头信息结构如下：struct reportheader {char id[8]; /* 站号 */float lat; /* 纬度*/float lon; /* 经度 */float t; /* 相对于该时次的相对单位值 */int nlev; /* 总层次 */int flag; /* 表示有无高空资料 */};上面的定义是C语言中的定义，转换到fortran中就是我们开始定义的几个变量，但是还没有定义t(避免和总时次混淆，定义为tim)、nlev和flag，所以先去上面的定义变量的地方补上这上个变量，把real那一行改为：real lon,lat,var, rain,tim ，在real行下面增加一行：integer nlev,flag。
