- 积分
- 26
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2014-5-1
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
我-想要利用WRF的grads输出结果画气压垂直剖面的地形遮挡,然后看了http://bbs.06climate.com/forum.php?mod=viewthread&tid=32926这个帖子之后,自己画了一个,结果图出来是这样的,还请各位大神帮忙指点一下。谢谢!!!!
'reinit'
'open /20170802_0xxxd1.ctl'
it = 30
'set font 4'
while(it<= 50)
'set grads off '
'set grid off'
'set z 1 16'
lon1=118
lon2=126
lat1=40.3
lat2=40.3
'set lon 118 126'
'set lat 40.3'
*'set lev 1000 200'
'set t 'it' '
'define pres=lev'
*'define tk=tc+273.15'
'define es=(6.112*exp(17.67*(tk-273.15)/(tk-29.65)))'
'define q=rh*(0.62197*es/(pres-0.378*es))/100.'
'define e=pres*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(tk)-log(e)-4.805)'
'define the=tk*pow((1000/pres),(0.2854*(1.0-0.28*q)))'
'define eqt=the*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
*nlev=16
'set xlint 0.5'
'set xlpos -20'
'set ylpos -20'
'set xlopts 1 6 0.15'
'set ylopts 1 6 0.15'
'set clopts 1 6 0.15'
'set gxout shaded'
'set csmooth on'
'set parea 3 9.75 1 7.75'
*'set cmin 5'
*'set cint 10'
'set rgb 201 166 242 143'
'set rgb 202 61 186 61'
'set rgb 203 97 184 255'
'set rgb 204 0 0 255'
'set rgb 205 250 0 250'
'set rgb 206 128 0 64'
'set clevs 5 15 25 35 45 55'
'set ccols 0 201 202 203 204 205 206'
'd dbz'
'set gxout contour'
'set cint 4'
'set clskip 2'
'set line 1 1 4'
'd eqt'
'cbarn 0.9 1 10.0 4.5'
'set gxout contour'
'set cthick 20'
'set ccolor 1'
'set line 1 1 30'
'set arrscl 0.5 40'
'd skip(v,8,1);w '
'run axis.gs -type B -position o -start 118 -end 126 -interval 2 -lsize 0.15 -lthick 0.8 -suffix `3.'
*'run axis.gs -type b -position i -start 118 -end 126 -suffix `3.'
'run axis.gs -type L -position o -start 1000 -end 200 -interval -100 -lsize 0.15 -lthick 0.8 '
*'run axis.gs -type b -position i -start 118 -end 126 -suffix `3.'
*'run axis.gs -type L -position o -start 20 -end 55 -suffix `3.'
ln=lon1
lt=lat1
lt1=lat1
ln1=lon1
'set z 1'
rc=terrain1(nlev,lt,ln1,lt1,lon1,lon2,lat1,lat2)
*'gxprint test.ps white'
'printim /sya/g2/jinwei/run/plot/2017080200__theta_v_w_'it'-f1(d1).gif white x1024 y768'
'c'
it=it+1
endwhile
'quit'
**************************************************************
function terrain(nlev,ln,ln1,lt1,lon1,lon2,lat1,lat2)
* this terrain adjustment is specified for the parea used here
* when scaling the terrain height if parea is different
* use ymax-ymin in parea(eg. 7.75-1=6.75)
while (ln<lon2)
ln=ln+0.005
lt=lat1 + (lat2-lat1)*(ln-lon1) / (lon2-lon1)
* calculates the lat/lon for which 2D terrain is extracted
'q w2xy 'ln1' 'lt1
x1=subwrd(result,3)
'set lon 'ln1
'set lat 'lt1
'set z 1'
'd hgt'
res =sublin(result,2)
mt=subwrd(res,4)
*check if the actual terrain height is higher than the domain height and
*if it is set it to that height so that terrain doesn't go outside the cross-section
* if (mt>height); mt=height;endif
* t1=mt*9.75/height+1
t1=hgt2pres(nlev,mt,ln1,lt1)
'q w2xy 'ln' 'lt
x2=subwrd(result,3)
'set lon 'ln
'set lat 'lt
'set z 1'
'd hgt'
res =sublin(result,2)
mt=subwrd(res,4)
* if (mt>height); mt=height;endif
* t2=mt*9.75/height+1
t2=hgt2pres(nlev,mt,ln,lt)
'set line 1 1 5'
'draw polyf 'x1' 1 'x1' 't1' 'x2' 't2' 'x2' 1'
lt1=lt
ln1=ln
endwhile
return
function terrain1(nlev,lt,ln1,lt1,lon1,lon2,lat1,lat2)
* this terrain adjustment is specified for the parea used here
* when scaling the terrain height if parea is different
* use ymax-ymin in parea(eg. 7.75-1=6.75)
while (lt<lat2)
lt=lt+0.005
ln=lon1 + (lon2-lon1)*(lt-lat1) / (lat2-lat1)
* calculates the lat/lon for which 2D terrain is extracted
'q w2xy 'ln1' 'lt1
x1=subwrd(result,3)
'set lon 'ln1
'set lat 'lt1
'set z 1'
'd hgt'
res =sublin(result,2)
mt=subwrd(res,4)
*check if the actual terrain height is higher than the domain height and
*if it is set it to that height so that terrain doesn't go outside the cross-section
* if (mt>height); mt=height;endif
* t1=mt*9.75/height+1
t1=hgt2pres(nlev,mt,ln1,lt1)
'q w2xy 'ln' 'lt
x2=subwrd(result,3)
'set lon 'ln
'set lat 'lt
'set z 1'
'd hgt'
res =sublin(result,2)
mt=subwrd(res,4)
* if (mt>height); mt=height;endif
* t2=mt*9.75/height+1
t2=hgt2pres(nlev,mt,ln,lt)
'set line 1 1 5'
'draw polyf 'x1' 1 'x1' 't1' 'x2' 't2' 'x2' 1'
lt1=lt
ln1=ln
endwhile
**************************************
return
function hgt2pres(nlev,mt,lon1,lat1)
'set z 1'
'set lon 'lon1
'set lat 'lat1
'd height*1000.0'
res =sublin(result,2)
h=subwrd(res,4)
if(mt <= h)
return 1
else
nl=1
h=0
while(h <= mt | h < -999999.0)
'set z 'nl
'set lon 'lon1
'set lat 'lat1
'd height*1000.0'
res =sublin(result,2)
h=subwrd(res,4)
'd lev'
res =sublin(result,2)
pres=subwrd(res,4)
nl=nl+1
endwhile
h2=h
pres2=pres
l2=nl-1
l1=nl
'set z 'l1
'set lon 'lon1
'set lat 'lat1
'd height*1000.0'
res =sublin(result,2)
h1=subwrd(res,4)
'd lev'
res =sublin(result,2)
pres1=subwrd(res,4)
tl=(pres1-pres2)/(h1-h2)
pres=(pres2) - tl*(h2-mt)
l3=nl-2
'set z 'l3
'd lev'
res =sublin(result,2)
pres3=subwrd(res,4)
h3=(pres3-pres2)/tl+h2
**NOTE: 1000.0 is pressure(z=1),100.0 is Y_direction interveral(hPa)
***** 9.75 should be 6.75=7.75-1 ,but 6.75 was too low
***** 9 was Y-direction tab number 9=(1000.0 - 100)/100.0
***** qiuxx 2015/1/19
new_hgt = (1000-pres)/100.0*8.75/8
if(new_hgt <= 0)
return 1
else
return new_hgt+1
endif
endif
return
|
-
|