利用WRF的grads输出结果画气压垂直剖面的地形遮挡的问题

莎朗嘿哟YY

我现在是想要利用WRF的grads输出结果画气压垂直剖面的地形遮挡，然后看了http://bbs.06climate.com/forum.php?mod=viewthread&tid=32926这个帖子之后，自己画了一个，结果图出来是这样的，然后grads产生了warning.这个gs里面的子程序中有一些不能理解，我直接套用的那个帖子的子程序。还请各位大神帮忙指点一下。谢谢！！！！

gs如下：
'reinit'
'open 1900HGT.ctl'
'enable print 1900.gmf'
'set grads off '
'set grid off'
'set z 1 19'
lon1=104.0
lon2=108.0
lat1=30.0
lat2=30.0
'define pressure=lev'
'define vor=hcurl(U,V)*100000'
*'set lev 1000 200'
'set lat 'lat1''
'set lon 'lon1' 'lon2''
'set font 1'
nlev=19
'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 clevs -50 -40 -30 -20 -10 0 10 20 30 40 50'
'd smth9(smth9(smth9(vor)))'
'run axis.gs -type B -position o -start 104 -end 108 -interval 0.5 -lsize 0.15 -lthick 0.8 -lfont 4 -suffix `3.'
'run axis.gs -type L -position o -start 1000 -end 200 -interval -100 -lsize 0.15 -lthick 0.8 -lfont 4 '
'cbarn 0.9 1 10.0 4.5'
ln=lon1
lt=lat1
lt1=lat1
ln1=lon1
'set z 1'
rc=terrain(nlev,lt,ln1,lt1,lon1,lon2,lat1,lat2)
'print'
'disable print'
;
**************************************************************

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 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

  new_hgt = (1000-pres)/100.0*8.75/8
  if(new_hgt <= 0)
   return 1
  else
   return new_hgt+1  
  endif
endif
return

小虹

把gs内容贴出来嘛，不然还要花钱下载

莎朗嘿哟YY

小虹 发表于 2015-3-6 15:21
把gs内容贴出来嘛，不然还要花钱下载

已经贴出来了，麻烦大神们帮忙看一下！谢谢

90后

画阴影时设了虚页面（'set parea 3 9.75 1 7.75' ），画完阴影后叠加地形不知道还有没有用，位置和尺寸貌似没对上。

river

可以去掉'set parea 3 9.75 1 7.75'  先试一下

莎朗嘿哟YY

river 发表于 2015-3-7 08:49
可以去掉'set parea 3 9.75 1 7.75'  先试一下

把这个去掉之后还是这样的错误。我这个贴子里的连接中那个楼主有这样的注释

1. function terrain(nlev,ln,ln1,lt1,lon1,lon2,lat1,lat2)
   
2. 
   
3. * this terrain adjustment is specified for the parea used here
   
4. * when scaling the terrain height if parea is different
   
5. * use ymax-ymin in parea(eg. 7.75-1=6.75)

复制代码

且那个楼主还说，注意修改note后面的部分；

1. **NOTE: 1000.0 is pressure(z=1),100.0 is Y_direction interveral(hPa)
   
2. ***** 9.75 should be 6.75=7.75-1 ,but 6.75 was too low
   
3. ***** 9 was Y-direction tab number 9=(1000.0 - 100)/100.0
   
4. *****   qiuxx   2015/1/19
   
5. 
   
6.   new_hgt = (1000-pres)/100.0*8.75/8
   
7.   if(new_hgt <= 0)
   
8.    return 1
   
9.   else
   
10.    return new_hgt+1  
    
11.   endif
    
12. endif
    
13. return

复制代码

这一块我不是很明白。这个new_hgt的公式是和页面设置有关系吗？

莎朗嘿哟YY

river 发表于 2015-3-7 08:49
可以去掉'set parea 3 9.75 1 7.75'  先试一下

我已经找到错误了，低级错误。子函数中的一个变量写错了，。。。

ogions

因为Y方向parea设置从1到7.75，所以理论上应该是6.75，但是实际效果来看，可能6.75偏小，需要稍调大一点比较好，这个值需要试验才好，根据地形高度，大致可以判断接近哪个气压层。

莎朗嘿哟YY

ogions 发表于 2015-3-7 22:23
因为Y方向parea设置从1到7.75，所以理论上应该是6.75，但是实际效果来看，可能6.75偏小，需要稍调大一点比 ...

嗯，好的，谢谢！

wolfwang

不能做个固定地形的剖面图？