GRADS计算MPV的GS

zhi726

一直在网络上找GRADS计算MPV的程序，但一直没找到，最后自己写了一个，图形感觉差不多，但还没经过严格的检验，现在放在这里，大家觉得有问题请指正。
由于不知道怎么一次完成lev的循环，所以程序分为两步：
第一步计算不同层次的相当位温eqt，然后fwrite；
第二步计算不同层次的湿位涡，然后fwrite；
如果要显示还要再写一个GS文件

eqt文件：
'reinit'
'open E:\20140901\fnl\fnl.ctl'
'set mpdset cnworld'
'set grid off'
'set grads off'
**************************eqt caculate start********************************
'set gxout fwrite'
'set fwrite E:\20140901\fnl\mpv\eqt.dat'
'set lon 100 120'
'set lat 20 40'
i=1
while(i<=9)
'set t 'i''
j=1
while(j<=21)
'set z 'j''
'define prs=lev'
'define es=(6.112*exp(17.67*(TMPprs-273.15)/(TMPprs-29.65)))'
'define q=RHprs*(0.62197*es/(prs-es))/100.'
'define e=(prs/100.)*q/(0.62197+q)+1e-10'
'define tlcl=55.0+2840.0/(3.5*log(TMPprs)-log(e)-4.805)'
'define theta=TMPprs*pow((1000./prs),(0.2854*(1.0-0.28*q)))'
'define eqt=theta*exp(((3376./tlcl)-2.54)*q*(1.0+0.81*q))'
'set undef -9999'
'd eqt'
'undefine prs'
'undefine es'
'undefine q'
'undefine e'
'undefine tlcl'
'undefine theta'
'undefine eqt'
j=j+1
endwhile
i=i+1
endwhile
**************************eqt caculate end*********************************
'disable fwrite'
'reinit'
'reinit'

mpv生成文件：
'reinit'
'open E:\20140901\fnl\fnl.ctl'
'open E:\20140901\fnl\mpv\eqt.ctl'
'set mpdset cnworld'
'set grid off'
'set grads off'
**************************mpv caculate start********************************
'set gxout fwrite'
'set fwrite E:\20140901\fnl\mpv\mpv.dat'
'set lon 100 120'
'set lat 20 40'
'set z 1'
'set t 1'
i=1
while(i<=9)
'set t 'i''
j=2
while(j<=20)
'set z 'j''
'define vor=hcurl(UGRDprs,VGRDprs)'
'define f=2*7.272e-5*sin(lat*3.14159/180)'
'define g=9.8'
'define dthep=eqt.2(z-1)-eqt.2(z+1)'
'define dthex=cdiff(eqt.2,x)'
'define dthey=cdiff(eqt.2,y)'
'define dup=UGRDprs.1(z-1)-UGRDprs.1(z+1)'
'define dvp=VGRDprs.1(z-1)-VGRDprs.1(z+1)'
'define er=6.371e6'
'define dx=cdiff(lon,x)*3.14159/180*er*cos(lat*3.14159/180)'
'define dy=cdiff(lat,y)*3.1416/180*er'
'define dp=(lev(z-1)-lev(z+1))*100'
'define mpv1=-1.*g*vor*dthep/dp-g*f*dthep/dp'
'define mpv2=g*((dvp/dp)*(dthex/dx)-(dup/dp)*(dthey/dy))'
'define mpv=(mpv1+mpv2)'
'set undef -9999'
'd mpv1'
j=j+1
endwhile
j=2
while(j<=20)
'set z 'j''
'define vor=hcurl(UGRDprs,VGRDprs)'
'define f=2*7.272*sin(lat*3.14159/180)'
'define g=9.8'
'define dthep=eqt.2(z-1)-eqt.2(z+1)'
'define dthex=cdiff(eqt.2,x)'
'define dthey=cdiff(eqt.2,y)'
'define dup=UGRDprs.1(z-1)-UGRDprs.1(z+1)'
'define dvp=VGRDprs.1(z-1)-VGRDprs.1(z+1)'
'define er=6.371e6'
'define dx=cdiff(lon,x)*3.14159/180*er*cos(lat*3.14159/180)'
'define dy=cdiff(lat,y)*3.1416/180*er'
'define dp=(lev(z-1)-lev(z+1))*100'
'define mpv1=-1.*g*(vor+f)*dthep/dp'
'define mpv2=g*((dvp/dp)*(dthex/dx)-(dup/dp)*(dthey/dy))'
'define mpv=(mpv1+mpv2)'
'set undef -9999'
'd mpv2'
j=j+1
endwhile
j=2
while(j<=20)
'set z 'j''
'define vor=hcurl(UGRDprs,VGRDprs)'
'define f=2*7.272*sin(lat*3.14159/180)'
'define g=9.8'
'define dthep=eqt.2(z-1)-eqt.2(z+1)'
'define dp=(lev(z-1)-lev(z+1))*100'
'define thep=dthep/dp'
'set undef -9999'
'd thep'
j=j+1
endwhile
**************************mpv caculate end*********************************
i=i+1
endwhile
'disable fwrite'
'reinit'
'reinit'

传说中的谁

论坛原来也有的，可以对比一下
感谢楼主分享

river

感谢分享，有人分享过可以拿来试试看一不一样

zhi726

，早知道有就不用花那么多时间去查公式了

nzy718d

谢谢lz

秀才哥哥

请问楼主，在你计算MPV的那一部分，d了三个变量，一个是mpv1，一个是mpv2，为何两次mpv1的计算不一样？还有最后一个thep，是啥？

Kazoo

感谢分享，正要用

sunflower

感谢楼主分享，学习一下哈

69846421

谢谢分享，赞赞赞赞！

69846421

谢谢分享，赞赞赞赞！