Re: Divergence on a pressure/latitude grid

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
         rad = 4.*atan(1.)/180. ; radians
         rcon = re*rad
                                 ; 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

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
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,
Adam

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

-- 
--------------------------------------------------------------
Adam Phillips			             asphilli_at_ucar.edu
National Center for Atmospheric Research   tel: (303) 497-1726
ESSL/CGD/CAS                               fax: (303) 497-1333
P.O. Box 3000				
Boulder, CO 80307-3000	  http://www.cgd.ucar.edu/cas/asphilli
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk