center_finite_diff
Performs a centered finite difference operation on the rightmost dimension.
Prototype
function center_finite_diff ( q : numeric, r : numeric, rCyclic : logical, Option ) return_val [dimsizes(q)] : numeric
Arguments
qA multi-dimensional array.
rA scalar, one-dimensional, or multi-dimensional array containing the coordinates along which q is to be differenced. Does need not be equally spaced from a computational point of view.
- scalar: r assumed to be the (constant) distance between adjacent points.
- one-dimensional (and the same size as the rightmost dimension of q): applied to all dimensions of q.
- multi-dimensional: then it must be the same size as q.
True: q treated as cyclic in r and the end values and all the returned values will be calculated via centered differences. False: q NOT treated as cyclic in r and the end values will use a one-sided difference scheme for the end points. q should not include a cyclic point.
result(n) = (q(n+1)-q(n))/(r(n+1)-r(n)) for the initial value result(m) = (q(m)-q(m-1))/(r(m)-r(m-1)) for the last valueOption
Reserved for future use. Currently not used. Set to an integer.
Description
Performs a centered finite difference operation on the rightmost dimension. If missing values are present, the calculation will occur at all points possible, but coordinates which could not be used will set to missing.
result(n) = (q(n+1)-q(n-1))/(r(n+1)-r(n-1))
Examples
Example 1
q = (/30,33,39,36,41,37/) r = 2. ; constant dqdr = center_finite_diff (q,r,False,0) ; result dqdr(0) = (33-30)/2 = 1.5 dqdr(1) = (39-30)/4 = 2.25 dqdr(2) = (36-33)/4 = 0.75 dqdr(3) = (41-39)/4 = 0.5 dqdr(4) = (37-36)/4 = 0.25 dqdr(5) = (37-41)/2 = -2.0Example 2:
theta = (/ 298,299,300,302,...,345,355,383/) ; potential temp. p = (/1000,950,900,850,...,200,150,100/) ; pressure (hPa) dtdp = center_finite_diff (theta,p,False,0)Example 3
p is not equally spaced. Demonstrates defining the appropriate pressure levels when the r coordinate is not equally spaced. The end points are computed via a one-sided difference.
theta = (/ 297, 298,300,302,...,345,350,355,383/) ; potential temp. p = (/1013,1000,900,850,...,200,175,150,100/) ; pressure (hPa) dtdp = center_finite_diff (theta,p,False,0) np = dimsizes(p) pPlot = new ( np, "float" ) ; arithmetic mean pPlot(0) = (p(0)+p(1))*0.5 ; set bottom pPlot(np-1) = (p(np-1)+p(np-2))*0.5 ; set top pPlot(1:np-2) = (p(0:np-3) + p(2:np-1))*0.5 ; mid points ; or ; log (mass) wgted pPlot(0) = exp((log(p(0))+log(p(1)))*0.5) ; set bottom pPlot(np-1) = exp((log(p(np-1))+log(p(np-2)))*0.5) ; set top pPlot(1:np-2) = exp((log(p(0:np-3)) + log(p(2:np-1)))*0.5) ; mid pointsExample 4
Let T be four-dimensional with dimensions time,level,lat,lon. P is one-dimensional. Perform the finite differencing only in the vertical (level) dimension. This requires that T be reordered to put level in the rightmost dimension.
dTdP = center_finite_diff(T(time|:,lat|:,lon|:,lev|:),P,False,0) ; returns dTdP (time,lat,lon,lev). this variable can be reordered ; again to place it back in the original orderExample 5
Now P is also four-dimensional and requires reordering:
dTdP = center_finite_diff (T(time|:,lat|:,lon|:,lev|:) \
,P(time|:,lat|:,lon|:,lev|:), False,0)
; returns dTdP(time,lat,lon,lev)
Example 6
Assume that the longitude
coordinate variable
associated with T in the examples above is cyclic and is equally spaced
in degrees but not in physical space.
dlon = (lon(2)-lon(1))*0.0174533 ; convert to radians
; pre-allocate space
dTdX = new ( (/ntim,klev,nlat,mlon/), typeof(T), T@_FillValue)
do nl=0,nlat-1 ; loop over each latitude
dX = 6378388.*cos(0.0174533*lat(nl))*dlon ; constant at this latitude
dTdX(:,:,nl,:) = center_finite_diff (T(:,:,nl;,:), dX , True,0)
end do
; result: dTdX(time,lev,lat,lon)