爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

搜索
查看: 4756|回复: 2

[分享资料] 【求助】同样的计算公式,为什么存储后画图和直接画图形状不一样?

[复制链接]
发表于 2021-12-8 23:28:45 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
存储数据脚本:
'reinit'
'open L:/DZQDL/2014datafiles/fnl_2014.ctl'


'set t 16'
'set lev 750'
'define tmpx=TMPprs'
'define rhx=RHprs'
'define prsx=lev'
'define esx=(6.112*exp(17.67*(tmpx-273.15)/(tmpx-29.65)))'
'define qx=rhx*(0.62197*esx/(prsx-esx))/100.'
'define ex=prsx*qx/(0.62197+qx)+1e-10'
'define tlclx=55.0+2840.0/(3.5*log(tmpx)-log(ex)-4.805)'
'define thetax=tmpx*pow((1000/prsx),(0.2854*(1.0-0.28*qx)))'
'define eqtx=thetax*exp(((3376./tlclx)-2.54)*qx*(1.0+0.81*qx))'
'define vx=VGRDprs'
'define ux=UGRDprs'

'set lev 650'
'define tmpy=TMPprs'
'define rhy=RHprs'
'define prsy=lev'
'define esy=(6.112*exp(17.67*(tmpy-273.15)/(tmpy-29.65)))'
'define qy=rhy*(0.62197*esy/(prsy-esy))/100.'
'define ey=prsy*qy/(0.62197+qy)+1e-10'
'define tlcly=55.0+2840.0/(3.5*log(tmpy)-log(ey)-4.805)'
'define thetay=tmpy*pow((1000/prsy),(0.2854*(1.0-0.28*qy)))'
'define eqty=thetay*exp(((3376./tlcly)-2.54)*qy*(1.0+0.81*qy))'
'define vy=VGRDprs'
'define uy=UGRDprs'

'set lev 700'
'define tmpz=TMPprs'
'define rhz=RHprs'
'define prsz=lev'
'define esz=(6.112*exp(17.67*(tmpz-273.15)/(tmpz-29.65)))'
'define qz=rhz*(0.62197*esz/(prsz-esz))/100.'
'define ez=prsz*qz/(0.62197+qz)+1e-10'
'define tlclz=55.0+2840.0/(3.5*log(tmpz)-log(ez)-4.805)'
'define thetaz=tmpz*pow((1000/prsz),(0.2854*(1.0-0.28*qz)))'
'define eqtz=thetaz*exp(((3376./tlclz)-2.54)*qz*(1.0+0.81*qz))'
'define dtx=cdiff(eqtz,x)'
'define dty=cdiff(eqtz,y)'                              
'define dx=cdiff(lon,x)*3.1416/180*cos(lat*3.1416/180)*6.371e6'
'define dy=cdiff(lat,y)*3.1416/180*6.371e6'

'define eqtxy=eqtx-eqty'
'define vxy=vx-vy'
'define uxy=ux-uy'

'define dp=10000'
'define mp1=-9.8*absvprs*eqtxy/dp'
'define mp2=9.8*((vxy/dp)*(dtx/dx)-(uxy/dp)*(dty/dy))'
'define mpv=mp1+mp2'
'define mpv2=mp2*1E6'
'define mpv1=mp1*1E6'

'set fwrite L:\DZQDL\rain\mpv20146620.grd'
'set gxout fwrite'
'set lev 700'
'set t 1'
'set lon 0 360.5'
'set lat -89.5 89.5'
'd mpv1'
'disable fwrite'
;
直接画图:
'reinit'
'open L:/DZQDL/2014datafiles/fnl_2014.ctl'


;i=1
;while(i<=31)
'set t 16'
'set gxout contour'
*'set display white color'
'set mpdset cnworld'
'set lon 90 110'
'set lat 15 35'
'set xlopts -10'
'set ylopts -10'
'set grid off'
'set xlint 10'
'set ylint 5'
'set lev 750'
'set ccolor 1'
'set cthick 4'
'define tmpx=TMPprs'
'define rhx=RHprs'
'define prsx=lev'
'define esx=(6.112*exp(17.67*(tmpx-273.15)/(tmpx-29.65)))'
'define qx=rhx*(0.62197*esx/(prsx-esx))/100.'
'define ex=prsx*qx/(0.62197+qx)+1e-10'
'define tlclx=55.0+2840.0/(3.5*log(tmpx)-log(ex)-4.805)'
'define thetax=tmpx*pow((1000/prsx),(0.2854*(1.0-0.28*qx)))'
'define eqtx=thetax*exp(((3376./tlclx)-2.54)*qx*(1.0+0.81*qx))'
'define vx=VGRDprs'
'define ux=UGRDprs'

'set lev 650'
'define tmpy=TMPprs'
'define rhy=RHprs'
'define prsy=lev'
'define esy=(6.112*exp(17.67*(tmpy-273.15)/(tmpy-29.65)))'
'define qy=rhy*(0.62197*esy/(prsy-esy))/100.'
'define ey=prsy*qy/(0.62197+qy)+1e-10'
'define tlcly=55.0+2840.0/(3.5*log(tmpy)-log(ey)-4.805)'
'define thetay=tmpy*pow((1000/prsy),(0.2854*(1.0-0.28*qy)))'
'define eqty=thetay*exp(((3376./tlcly)-2.54)*qy*(1.0+0.81*qy))'
'define vy=VGRDprs'
'define uy=UGRDprs'

'set lev 700'
'define tmpz=TMPprs'
'define rhz=RHprs'
'define prsz=lev'
'define esz=(6.112*exp(17.67*(tmpz-273.15)/(tmpz-29.65)))'
'define qz=rhz*(0.62197*esz/(prsz-esz))/100.'
'define ez=prsz*qz/(0.62197+qz)+1e-10'
'define tlclz=55.0+2840.0/(3.5*log(tmpz)-log(ez)-4.805)'
'define thetaz=tmpz*pow((1000/prsz),(0.2854*(1.0-0.28*qz)))'
'define eqtz=thetaz*exp(((3376./tlclz)-2.54)*qz*(1.0+0.81*qz))'
'define dtx=cdiff(eqtz,x)'
'define dty=cdiff(eqtz,y)'
'define dx=cdiff(lon,x)*3.1416/180*cos(lat*3.1416/180)*6.371e6'
'define dy=cdiff(lat,y)*3.1416/180*6.371e6'

'define eqtxy=eqtx-eqty'
'define vxy=vx-vy'
'define uxy=ux-uy'

'define dp=10000'
'define mp1=-9.8*absvprs*eqtxy/dp'
'define mp2=9.8*((vxy/dp)*(dtx/dx)-(uxy/dp)*(dty/dy))'
'define mpv=mp1+mp2'
'define mpv2=mp2*1E6'
'define mpv1=mp1*1E6'

'set grads off'
'set lon 96 107'
'set lat 20 31'

'set mpdset cnworld'
'set map 9 1 1'
'draw map'
'set gxout shaded'
*'run  L:\DZQDL\rainbow1.gs'
'd mpv1'
'run cbarn_interp.gs'

'set gxout contour'
'set ccolor 2'
'd mpv2'

'q dim'
it=sublin(result,5)
it=subwrd(it,6)
'run axis.gs -type b -position i -label on -interval 3 -sinterval 2 -lsize 0.2  -voffset -0.2'
'run axis.gs -type l -position i -label on -interval 3 -sinterval 2 -lsize 0.2'
'draw title mpv1 700hpa t 'it''
'printim L:/DZQDL/2014datafiles/mpv12h700/mpv12_700hpa2_'it'.jpg white'
;'c'
;i=i+1
;endwhile
'disable gxprint'
*

存储后画MPV1

存储后画MPV1

直接画MPV1(阴影)

直接画MPV1(阴影)

左图为存储后画的MPV1,右图为直接画的MPV1(阴影)
可以看出两个纹路走向是不一样的,不知道是什么原因

密码修改失败请联系微信:mofangbao
发表于 2021-12-9 10:18:15 | 显示全部楼层
把两个数据都读出来,算下插值,确认下存的对不对
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

 楼主| 发表于 2021-12-9 22:47:46 | 显示全部楼层
croton 发表于 2021-12-9 10:18
把两个数据都读出来,算下插值,确认下存的对不对

检查了一下就是在存储这步上有问题,但是不知道怎么处理
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表