# Re: Divergence on a pressure/latitude grid

Date: Thu, 13 Sep 2007 09:44:30 -0600

Hi Erwan,

Presumably you want something like:

ep_div = dv/dlat + dw/dp

This seems to be a rather straight forward programming exercise.

=====

function dv_lat_p(v[*][*]:numeric, w[*][*]:numeric, lat[*], p[*])
; ============
; UNTESTED: user must verify that units and finite diff are correct
; ===========
; nomenclature:
; v,w - meridional and vertical components
; These MUST have the "_FillValue" assigned and
; coordinate variables "lat" and "p"
;
; v(klvl,nlat) , w(klvl,nlat)
; lat - south to north order
; p - bottom to top [Pa]
;
; Note: ; the "top" and "bottom" lwvels will be returned as missing

local re, rad, rcon, dlat, dims, nlat, klvl, dv, kl, nl, dp, dlat
begin
dims = dimsizes(v)
klvl = dims(0)
nlat = dims(1)

dv = v ; create a return array the same size as v
; also this transfers all coordinate information
dv = v@_FillValue ; initialize all to missing value

re = 6.37122e06 ; radius of earth
; delta latitude [array syntax]
; 1/(2*dlat)
dlat = new ( nlat, typeof(lat), -999)
dlat(1:nlat-2) = rcon*(lat(2:nlat-1)-lat(0:nlat-3)) ; meters
dlat_at_units = "m"
;print ( dlat )

; delta pressure
dp = new (klvl, typeof(p), -999)
dp(1:klvl-2) = (p(0:klvl-3) - p(2:klvl-1))
;dp_at_units = p_at_units
;print ( dp )
; compute divergence using centered differences
; classic do notation [clear]
do kl=1,klvl-2
do nl=1,nlat-2
dv(kl,nl) = (v(kl,nl+1)-v(kl,nl-1))/dlat(nl) \
+ (w(kl+1,nl)-w(kl-1,nl))/dp(kl)
end do
end do
; array notation [slightly faster]
;;do kl=1,klvl-2
;; dv(kl,1:nlat-2) = (v(kl,2:nlat-1) -
v(kl,0:nlat-3))/dlat(1:nlat-2) \
;; + (w(kl+1,1:nlat-2) - w(kl-1,1:nlat-2))/dp(kl)
;;end do

dv_at_long_name = "Monier Divergence"
dv_at_units = "?"

return (dv)
end

begin
nlat = 64
lat = latGau(nlat, "lat", "latitude", "degrees_north")
p = (/1000,925,850,700,600,500,300,200,100,70,50,20,10/) ; hPa
p = p*100 ; Pa
klvl = dimsizes(p)

v = random_normal ( 10, 15, (/klvl,nlat/))
w = random_normal ( 2, 10, (/klvl,nlat/))
v@_FillValue = 1.e20
w@_FillValue = 1.e20

EPDIV= dv_lat_p(v, w, lat, p)
printVarSummary( EPDIV )
end

Good luck,

Erwan Monier wrote:
> Hi,
>
> I would like to calculate the divergence of the Eliassen-Palm flux
> vector, which has vertical and meridional coordinates. I have tried to
> use the function uv2dvF and trick the script to think that the vertical
> coordinate corresponds to the longitudinal coordinate. However, since
> the vertical coordinate is on pressure levels
> (1000,925,850,700...,100,70,50,20,10,...), the grid is neither fixed nor
> gaussian and therefore the divergence is not properly calculated. Is
> there anything I can do to properly calculate the divergence in a
> height/pressure vs latitude cross section vector ?
>
> Please let me know if you have any advice on this issue. Thank you,
>
> Erwan
> _______________________________________________
> ncl-talk mailing list
> ncl-talk_at_ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

```--
--------------------------------------------------------------