Calculating vertical gradient of theta...???

[Thomas Spengler]

Dear all,

I am interested in calculating the brunt vaisala frequency and hence
need the vertical gradient of potential temperature.
My dataset is in pressure coordinates and I can simply define theta
(potential temperature) as

'define th=t*pow((1000/lev),(287/1004))'

In order to calculate dth/dp you would simply use

'define dthdp=(th(z-1)-th)/(100*(lev(z-1)-lev))'

But that's not working because GrADS has problems with shifting the
levels of a self-defined variable!!!
(problem with th(z-1))

To be able to calculate dthdp you need to define a th1 in the following way

'define th1=t(z-1)*pow((1000/lev(z-1)),(287/1004))'

to be able to finally calculate dthdp

'define dthdp=(th1-th)/(100*(lev(z-1)-lev))'

My question now is:
Why is GrADS not able to shift coordinates in self-defined variables?
It would be much more convenient to be able to simply use th(z-1)
instead of a new defined variable th1.

A simple example illustrating the problem of GrADS is defining t1 as

'define t1=t(z-1)'

and then defining t2 in a way so that it should represent t again

'define t2=t1(z+1)'

You can easily convince yourself with this simple test that there is a
big problem within GrADS since t2 will not be equal to t.
(I am working with the 1.8 S11 version of GrADS)

Would be great if someone has a solution to the problem.

Best regards,
Thomas