re: Potential Vorticity and spiciness

[Jonathan Callahan <callahan@xxxxxxxxxxxxx>]

Paulo,

I contacted Jennifer Miletta Adams about this and she sent me back an
old journal script she had used.  Here's her response.  I hope this is a
good starting place.


-- Jonathan Callahan


>  What I calculated was Isentropic
> Potential Vorticity -- not sure whether this is exactly what
> Paolo wants, but that's what I've got. I have no idea what
> spiciness means...
> 
> Anyway, I've copied  the .jnl file below. Looks like the
> only argument is a value for the theta surface. I'm sure
> that temperature data file is no longer around, but
> otherwise it should work. Feel free to edit and post to
> ferret-users as you see fit. 
> 
> --Jennifer

----------------------------------------
Jennifer Miletta Adams
Center for Ocean-Land-Atmosphere Studies
4041 Powder Mill Road, Suite 302
Calverton, MD 20705-3106
Phone: (301)902-1278 or (301)595-7000
Fax:   (301)595-9793
Email: jma@cola.iges.org


\can mode verify
use tmpu_winter8586.cdf
let/qui g = 9.81
let/qui p0 = 1000
let/qui omega = 7.292e-5
let/qui pi = 3.14159
let/qui deg2rad = pi/180.

let/qui f = 2 * omega *
sin(deg2rad*y[g=tmpu,d=tmpu_winter8586])

let/qui theta =
tmpu[d=tmpu_winter8586]*((p0/z[g=tmpu,d=tmpu_winter8586])^0.286)
set var/title="Potential Temperature"/unit=K theta
let/qui theta_$1 = theta[z=@weq:$1]

use uwnd_winter8586.cdf
let/qui uwnd_$1 = uwnd[d=uwnd_winter8586] * theta_$1
set var/tit="U-Wind on the $1K Isentrope" uwnd_$1

use vwnd_winter8586.cdf
let/qui vwnd_$1 = vwnd[d=vwnd_winter8586] * theta_$1
set var/tit="V-Wind on the $1K Isentrope" vwnd_$1

let/qui dthetadp = theta[z=@ddc]
let/qui dthetadp_$1 = dthetadp * theta_$1
let/qui utheta = uwnd_$1[z=@sum]
let/qui vtheta = vwnd_$1[z=@sum]
let/qui zeta_$1 = vtheta[x=@ddc] - utheta[y=@ddc]

! add the 1e4 factor to correct the units to PV Units
let/qui ipv_$1 = 1e4 * (-1)*g * (f + zeta_$1) *
dthetadp_$1[z=@sum]

fill/x=-180:180/tit="Potential Vorticity on the $1 K Theta
Surface" IPV_$1
cont/over/lev=(2,10,1,-3)/nolab IPV_$1
cont/over/lev=(2,2,1,-3)/pen=7/nolab IPV_$1