请教WRF的grads输出结果画气压垂直剖面的地形遮挡

lnyk_jw

我-想要利用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

lnyk_jw

地形出不来.

小其其格

地形绘制函数可以参照这个http://bbs.06climate.com/forum.php?mod=viewthread&tid=22506

出图效果：

                               
登录/注册后可看大图

lnyk_jw

非常谢谢你给予的提示，

术舟

好东西，不用金币