Grads 使用笔记
求导———
define dv = cdiff(v,x)
define dx = cdiff(lon,x)*3.1416/180
define du = cdiff(u*cos(lat*3.1416/180),y)
define dy = cdiff(lat,y)*3.1416/180
display (dv/dx-du/dy)/(6.37e6*cos(lat*3.1416/180))
************************************************************************************************
矢量图
之————控制矢量图中矢量标志的位置
'set gxout vector'
'set arrlab off' ＊关掉缺省的图标
'set arrscl 0.5 0.5' ＊表示在图中0.5英寸的长度代表0.5m/s
'd u;v'
'draw string 0.95 3.6 `3"'＊画箭头，使用set font 3的特殊图形
'draw string 5.75 0.8 0.5'＊写箭头下文字
或
'set arrlab off' ＊关掉缺省的图标
'set arrscl 0.5 0.5' ＊表示在图中0.5英寸的长度代表0.5m/s
'd u;v'
'draw line 5.5 1 6 1' ＊画线
'draw line 5.95 0.95 6 1'＊画箭头
'draw line 5.95 1.05 6 1'
'draw string 5.75 0.8 0.5'＊写箭头下文字
之————规定矢量图输出的阈值，使用set cmin/cmax
但要注意的是，display时，必须是
'd u;v;mag(u,v)'
或
'd u;v;sqrt(u*u,v*v)'
而不能仅仅'd u;v'
之————控制矢量密度与颜色
'set ccolor 1'*控制颜色
'd skip(u,2);v'*控制密度
*******************************************************************************************
function
之————gs重复使用到的语句，打包成function
以function main(track)开头
以return结尾
*******************************************************************************************
函数
smth9对于资料分辨率过低引起的不平滑起作用，但如果对于模式输出等高分辨率的不平滑不会有任何效果
多幅图输出
之————多幅图输出时，边缘的经纬度为了避免重叠，必须使得第n张图的最后一个坐标不输出，
如set lon 132 144
set xlint 3
set xlopts color
*因为144是（144-132）/3的公倍数，144必定会在最后一个经度坐标上出现，这样就会与第n+1个图的第一个横坐标：132重叠，解决方法有:
1.在最初的计算区域，就是分页前，就规定所选区域为[132,144)
2.分页前，依然选取[132,144]，在分页后，每幅图的display前，将区域选为[132,144)
*******************************************************************************************
经纬度对应图上位置——
之1.
lon = '133.5 '
lat = '33.8'
'q w2xy '%subwrd(lon,1)%' '%subwrd(lat,1)
x = subwrd(result,3)
y = subwrd(result,6)
之2.
lon='xxx.x xxx.x xxx.x'
lat='xx.x xx.x xx.x'
j=1
cn=999
while(j
x = subwrd(result,3)
y = subwrd(result,6)
'draw mark 5 'x' 'y' 0.1'
j=j+1
endwhile
之3.
i=1
cnt=999
while(i<=cnt)
aa=read('.txt')
aa1=sublin(aa,2)
lat=subwrd(aa1,3)
lon=subwrd(aa1,4)
'q w2xy 'lon' 'lat
x = subwrd(result,3)
y = subwrd(result,6)
········
i = i + 1
endwhile
ret=close('.txt')
*******************************************************************************************
后处理
之——色标
cbarc 在图的右上角给出了标尺
cbarn 在图的下方给出标尺
之——设置特殊字体与符号
用set font n设置字体
font 3时，可以输出很多symbols，
具体写法是:
'set font 3' *actually，this setting can be deleted，and it's strongly recommended not to write this command so that you can avoid unwanted character change
之字符：'draw string 9.09 3.03 `52900141' *then you can see bold roman numbers 2900141
之符号：'draw string xxx xxx `3#' then you can see a arrow,remember this the symbol must come together with the number "3"
*******************************************************************************************
关于画水汽通量，水汽通量为一个向量，方法如下：
set gxout shaded
d mag(q*u,q*v)
set gxout vector
d q*u;q*v （mag(q*u,q*v)只是画通量的绝对值，所以要加上d q*u;q*v才能看出通量的方向）
关于画水汽通量散度，应该为：
hdivg（q*u,q*v）
*******************************************************************************************
grads将数据写入不同类型文件：
之——存成ASCII码：
file='output.txt'
'set gxout print'
'd var'
rc=write(file,result)
rc=close(file)
之——存成binary：
'set gxout fwrite'
'set fwrite filename'
'd var'
'disable fwrite'
********************************************************************************************
对于时间坐标的设置
（比如说画时间-经度或者纬度剖面图）
可以把时间为在写数据的时候按照y维或者x维进行读写
在ctl文件里面进行相同的设置
在画图的时候：要设置'set mpdraw off'不画地图底图
'set mproj off',地图投影方式同scaled，但此时x/y轴不表示经度/纬度，x/y轴标记也变成数值，如90°W变成270，,90°S变成-90
想要坐标轴标记按照自己的标记，则需要两个命令结合使用，'set yaxis(xaxis) start end incr'和'set ylabs lab1|lab2|...'
或者'set ylevs lab1 lab2...'和'set ylabs lab1|lab2|...'配合使用。
*********************************************************************************************
Grads批处理文件名
若想要提取从1951-2006年56年nc文件中的某些数据，一个一个处理非常麻烦，这里介绍种较为简易的方法。例如想提取6-8月的位势高度资料。
'reinit'
t5=1951
*作文件名循环
while(t5<=2006)
'set gxout fwrite'
'set fwrite D:\sichuan\hgt1\'%t5%'.dat'
'sdfopen e:\ncep1\hgt\hgt.'%t5%'.nc'
t3=t5-1950
*判断是否为闰年
if
（t3=2|t3=6|t3=10|t3=14|t3=18|t3=22|t3=26|t3=30|t3=34|t3=38|t3=42|t3=46|t3=50|t3=54)
to=153
else
to=152
endif
t4=to+91
while(to<=t4)
'set t 'to
t1=1
while(t1<=12)
'set z 't1
'set lon 80 140'
'set lat 15 55'
'd hgt'
t1=t1+1
endwhile
to=to+1
endwhile
*这里必须先观点上述运行的文件，grads最多同时可以打开20个文件左右。
'reinit'
t5=t5+1
endwhile
'reinit'
这样可以提取你想要的年数据，然后你大可运用fortran对数据进行随心所欲的处理。
***************************************************************
'set timelab on|off'
可以关闭或者显示右下角的时间标记。
****************************************************************
!!!!重要命令
'set tlsupp year|month'
控制时间坐标轴的标记
例如'set tlsupp year'就会去掉年的标记只留下月份。
*****************************************************************
将gmf文件转换成eps文件
->!gxeps -i d:/**/**.gmf -c -o d:/***/**.eps
文件路劲一定是反斜杠