计算画图和提取结果出图不一致

紫雨微澜

已解决！~
感谢各位提点~
计算直接出图和计算提出结果出图为什么结果差这么多？哪位大侠指教一下首先，是计算假相当位温直接出图，图和GS如下。
'reinit'
'sdfopen I:\daily\air.1979.nc'
'sdfopen I:\daily\rhum.1979.nc'
'set display color white'
'c'
'set xlopts 1 8 0.25'
'set ylopts 1 8 0.25'
'set clopts 1 8 0.25'
'set lon 70 140'
'set lat 15 55'
'set t 257'
'set lev 850'
'define t0=air.1'
'define rh=rhum.2'
'define prs=lev'
'define es=(6.112*exp(17.67*(t0-273.15)/(t0-29.65)))'
'define q=rh*(0.62197*es/(prs-es))/100.'
'define e=prs*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(t0)-log(e)-4.805)'
'define theta=t0*pow((1000/prs),(0.2854*(1.0-0.28*q)))'
'define eqt=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
'define thse=eqt-273.15'

'set ccolor 1'
'set cthick 16'
'set grid off'
'set grads off'
'd thse'
'printim I:\daily\ts.png'



其次，是计算假相当位温后先提取出来成grd,然后再画图，图和GS如下。
'reinit'

'open I:\daily\thse.ctl'
'set display color white'
'c'
'set xlopts 1 8 0.25'
'set ylopts 1 8 0.25'
'set clopts 1 8 0.25'

'set lon 70 140'
'set lat 15 55'
'set t 257'
'set lev 850'

'set ccolor 1'
'set cthick 16'
'set grid off'
'set grads off'

'd thse'
'printim I:\daily\thse-grd.png'



fwrtie假相当位温:
'reinit'
'sdfopen I:\daily\air.1979.nc'
'sdfopen I:\daily\rhum.1979.nc'
'set gxout fwrite'
'set fwrite I:\daily\thse.grd'

zz=1
  while (zz<=8)
  if zz=1   
'set lev 1000'
endif
     if zz=2   
'set lev 925'
endif
     if zz=3   
'set lev 850'  
endif
    if zz=4   
'set lev 700'
endif
     if zz=5   
'set lev 600'
endif
     if zz=6   
'set lev 500'
endif
  if zz=7   
'set lev 400'
endif
  if zz=8   
'set lev 300'
endif

    tt=1
  while (tt<=365)
'set t 'tt''

'define t0=air.1'
'define rh=rhum.2'
'define prs=lev'
'define es=(6.112*exp(17.67*(t0-273.15)/(t0-29.65)))'
'define q=rh*(0.62197*es/(prs-es))/100.'
'define e=prs*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(t0)-log(e)-4.805)'
'define theta=t0*pow((1000/prs),(0.2854*(1.0-0.28*q)))'
'define eqt=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
'define thse=eqt-273.15'
'd thse'
tt=tt+1
endwhile
zz=zz+1
endwhile

'disable fwrite'
'reinit'

@暴雨@版主  

river

提取资料你的循环都反了吧，外循环应该是时间

紫雨微澜

river 发表于 2013-11-20 17:25

                               
登录/注册后可看大图

提取资料你的循环都反了吧，外循环应该是时间

不是外循环高度或时间都可以吗？
我改下试试

紫雨微澜

river 发表于 2013-11-20 17:25

                               
登录/注册后可看大图

提取资料你的循环都反了吧，外循环应该是时间

晕，和之前上面2个图都不一样

lqouc

那把你的改过的gs什么的都贴出来看看吧。
话说你在循环高度的时候有必要用那么多判断么？直接set z不就好了？这个和set lev没区别的吧。

紫雨微澜

lqouc 发表于 2013-11-20 17:45

                               
登录/注册后可看大图

那把你的改过的gs什么的都贴出来看看吧。
话说你在循环高度的时候有必要用那么多判断么？直接set z不就好了 ...

因为之前是用Z的，但是计算时候pre=lev,用z结果不对，所以就用了判断
'reinit'
'sdfopen I:\daily\air.1979.nc'
'sdfopen I:\daily\rhum.1979.nc'
'set gxout fwrite'
'set fwrite I:\daily\thse2.grd'

    tt=1
  while (tt<=365)
'set t 'tt''

zz=1
  while (zz<=8)
  if zz=1   
'set lev 1000'
endif
     if zz=2   
'set lev 925'
endif
     if zz=3   
'set lev 850'  
endif
    if zz=4   
'set lev 700'
endif
     if zz=5   
'set lev 600'
endif
     if zz=6   
'set lev 500'
endif
  if zz=7   
'set lev 400'
endif
  if zz=8   
'set lev 300'
endif
  

*'set lon 70 140'
*'set lat 15 55'

'define t0=air.1'
'define rh=rhum.2'
'define prs=lev'
'define es=(6.112*exp(17.67*(t0-273.15)/(t0-29.65)))'
'define q=rh*(0.62197*es/(prs-es))/100.'
'define e=prs*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(t0)-log(e)-4.805)'
'define theta=t0*pow((1000/prs),(0.2854*(1.0-0.28*q)))'
'define eqt=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
'define thse=eqt-273.15'
'd thse'

zz=zz+1
endwhile
tt=tt+1
endwhile

'disable fwrite'
'reinit'

river

紫雨微澜 发表于 2013-11-20 17:27

                               
登录/注册后可看大图

不是外循环高度或时间都可以吗？
我改下试试

这个是grads的规定，资料排序必须那样排好吧

lqouc

嗯嗯，看起来好像真的是没啥问题啊。

紫雨微澜

lqouc 发表于 2013-11-20 17:59

                               
登录/注册后可看大图

嗯嗯，看起来好像真的是没啥问题啊。

难道是fwrite的对精度有影响？但也不会差这么多吧

紫雨微澜

river 发表于 2013-11-20 17:58

                               
登录/注册后可看大图

这个是grads的规定，资料排序必须那样排好吧

学习了，从来没想过这个问题