[gradsusr] how to calculate Potential Vorticity on pressure surface?

[Henri Pinheiro]

Dear GrADS users,

I have tried to calculate the PV on a 300 hPa surface using three variables,
u and v wind, and temperature.
The wind data are in a single isobaric level (300 hPa) and temperature data
are in three levels (400, 300 and 250 hPa), that because the calculation of
PV300 is made by the difference between lower and higher levels (400 - 250
hPa).
However, the PV output data is quite strange when it is opened by grads, I
mean I do not think the PV values are correct.
I also noticed that the wind data size are different than the PV output
size, look that:

-rw-r--r-- 1 henri users 27273216 2010-04-26 09:45 PV300_MAM1989.dat
-rw-r--r-- 1 henri users 27131904 2010-04-26 09:45 u300_MAM1989.dat
-rw-r--r-- 1 henri users 27131904 2010-04-26 09:45 v300_MAM1989.dat

Please, could anyone help me with this?
The grads script to compute PV is below...

Thanks in advance
Henri


****** Open the zonal and meridional wind, temperature and it calculates
****** Potential Vorticity (PV) on a pressure surface

'reinit'

'open file_u.ctl'
'open file_v.ctl'
'open file_temp.ctl'

'set gxout fwrite'
'set fwrite PV300_MAM1989.dat'

k=1
while (k<=368)
'set t 'k

'set lev 300'

* ########## Coriolis calculation
'define coriol=2*7.29e-5*sin(lat*3.1415/180)'

* ########## Relative Vorticity calculation
'define vort=hcurl(u300,v300.2)'

* ########## Potential Temperature (theta) calculation
* ########## tmat is Time Mean Absolute Temperature

'define
dt=tmat.3(z-1)*pow(1000/lev(z-1),0.286)-tmat.3(z+1)*pow(1000/lev(z+1),0.286)'
'define dp=100*(lev(z-1)-lev(z+1))'

* ########## differential of theta/pressure

'define dtdp=dt/dp'

* ########### PV calculation on pressure surface

'define pv = -9.8*(coriol+vort)*dtdp'

'd pv'

k= k+1
endwhile

'disable fwrite'

'set gxout shaded'
'set gxout contour'

'quit'



-- 
Henri Rossi Pinheiro
Meteorology Master Student
National Institute for Space Research (INPE)
Weather Forecast and Climate Studies Center (CPTEC)
Cachoeira Paulista, SP/Brazil, 12630-000 Tel.:+55-12-3186-8607
e-mails: henri.pinheiro at cptec.inpe.br
              henri.inpe at gmail.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://gradsusr.org/pipermail/gradsusr/attachments/20100426/5a903ee1/attachment-0003.html