[gradsusr] Moist Potential Vorticity

[PHILBERT LUHUNGA]

Hi I am writing a script to compute moist potential vorticity using entropy temperature Theta_s, and virtual potential temperature Theta_v. Please help me to check the script if it is coded right. When I drive it, it gives me results but am not sure the way I have coded if it is right. could anyone cross check it *our script to compute MPV700hPa.
*###################Day1################################
*'open $DATA_DIR/WRF_9KM_Nonhy_Topo${RunDay}.ctl'
'open 201112200000_arw_wrfout_d01.ctl'
 
'q file''q dims''pi=3.14159'
'dtr=pi/180'
'a=6.37122*pow(10,6)'
'omega=7.2921*pow(10,-5)'
'g=9.81'
'R=287'
'Cp=1004.5'
'L=2.453*pow(10,6)'
*'L=2.453e6'
'Rv=461'
* tk    = Temp in Kelvin
* mixr  = water vapour Mixing ratio
* u     = U-wind in m/s
* v     = V-wind in m/s
* define ug=UGRDprs   U-wind in m/s in pressure level
* define vg=VGRDprs  V-wind in m/s in pressure level
*'define ua=U-ug'
*'define va=V-vg'*PV @700hPa
'set lev 800'
'define f=2*omega*sin(lat*dtr)'
'define p1=lev*100'
'define tk1=TMPprs'
'define u1=UGRDprs'
'define v1=VGRDprs'
'define rh1=rhprs'
'define thetap1=tk1*pow((100000/p1),0.286)'
'define es=6.1173*100*exp((L/Rv)*(1/273.16 - 1/tk1))'
'define ep1=rh1*es/100'
'WVmixr1=0.622*ep1/(p1)'
'define T=tk1'
'define VirT1=T*(1+0.61*WVmixr1)'
'thetav1=VirT1*pow((100000/p1),0.2857))'
'thetas1=thetap1*exp(5.87*WVmixr1)''set lev 600'
'define p2=lev*100'
'define tk2=TMPprs'
'define u2=UGRDprs'
'define v2=VGRDprs'
'define rh2=rhprs'
'define thetap2=tk2*pow((100000/p2),0.286)'
'define es=6.1173*100*exp((L/Rv)*(1/273.16 - 1/tk2))'
'define ep2=rh2*es/100'
'WVmixr2=0.622*ep2/(p2)'
'define T=tk2'
'define VirT2=T*(1+0.61*WVmixr2)'
'thetav2=VirT2*pow((100000/p2),0.2857))'
'thetas2=thetap2*exp(5.87*WVmixr2)'
********************PV at 700hPa level
'set lev 700'
'dp=p2-p1'
'dudp=(u2-u1)/dp'
'dvdp=(v2-v1)/dp'
'define p0=lev*100'
'define tk0=TMPprs'
'define u0=UGRDprs'
'define v0=VGRDprs'
'vort0=hcurl(u0,v0)'
'define Totvort=vort0+f''define alpha0=R*tk0/p0'
'define rh0=rhprs'
'theta0=tk0*pow((100000/p0),0.2857))'
'define es=6.1173*100*exp((L/Rv)*(1/273.16 - 1/tk0))'
'de
fine ep0=rh0*es/100'
'WVmixr0=0.622*ep0/(p0)'
'define T=tk0'
'define VirT0=T*(1+0.61*WVmixr0)'
'thetav0=VirT0*pow((100000/p0),0.2857))'
'thetas0=theta0*exp(5.87*WVmixr0)'
'dy=cdiff(lat,y)*dtr*a'
'dx=cdiff(lon,x)*dtr*a*cos(lat*dtr)'
'define dthetas0dx=cdiff(thetas0,x)/dx'
'define dthetas0dy=cdiff(thetas0,y)/dy'
'define dthetas0dp=(thetas2-thetas1)/dp'
'define Totdthedtas=dthetas0dx+dthetas0dy+dthetas0dp''define dthetav0dx=cdiff(thetav0,x)/dx'
'define dthetav0dy=cdiff(thetav0,y)/dy'
'define dthetav0dp=(thetav2-thetav1)/dp'
'define Totdthedtav=dthetav0dx+dthetav0dy+dthetav0dp''term11= -g*Totvort*dthetav0dp*pow(10,6)'
'term22=-g*dudp*dthetav0dy*pow(10,6)'
'term33=g*dVdp*dthetav0dx*pow(10,6)'
'PV7002=term1+term2+term3''term1= -g*Totvort*dthetas0dp*pow(10,6)''term2=-g*dudp*dthetas0dy*pow(10,6)'
'term3=g*dVdp*dthetas0dx*pow(10,6)'
'PV700=term1+term2+term3'*'term1=Totvort*dthetas0dp*e6)'
*'term2=dudp*dthetas0dy*e6'
*'term3=-dVdp*dthetas0dx*e6'
'PV700=term1+term2+term3''set t 7'
'set gxout shaded'
'set mpdset hires'
'set display color white'
'set vpage 1.0 8.5 2.0 8.5'
'c'
*'d mp'd
'd theta_s'
*'exec colors'
*'run colorset.gs'
*'set ylint 12'
*'set xlint 12'
'cbarhov2.gs'
'draw title 700hPa-Moist Potential Vorticity ''printim 700_mv.png'
'set grads off'
*'quit'*'define MV=TotVort*ept(tmpprs,rhprs.prs)'
  AS EVER
LUHUNGA

--------------------
Philbert Modest Luhunga
University of Pretoria 
Department of Geography,Geoinformatics and Meteorology
Private Bag X20 Hatfield 0028 South Africa
Tel +27 (0) 12 420 5164
Fax +27 (0) 12 420 6385
Mobile:+ 27826228060
Email address: philuhunga at yahoo.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20141002/51d6cc01/attachment-0001.html