[ncl-talk] Problem with divergence calculation

[Mateus da Silva Teixeira]

Hi,

I'm trying to reproduce with *center_finite_difference* function the
computation of wind divergence, to compare with *uv2dv_cfd* function.
However, I'm getting some strange values for divergence.

For *Dv/Dy*, I'm using:

dv_dy = center_finite_diff_n( v850, lat*g2r, False, 0, 0 )/rTerra

For Du/Dx, I'm using:






*dlon = (lon(2)-lon(1))*g2rdu_dx = new( (/dimsizes(u850)/), typeof( u850 ),
u850 at _FillValue )term3 = new( (/dimsizes(v850)/), typeof(v850),
v850 at _FillValue )do i = 0, dimsizes(lat)-1      dx = rTerra * cos(
g2r*lat(i) ) * dlon      du_dx(i,:) = center_finite_diff( u850(i,:), dx,
True, 0 )*

*      term3(i,:) = ( v850(i,:)/rTerra ) * tan( lat(i)*g2r )  end do*

with *rTerra* being the mean radius of the Earth, and finally, for
divergence (using the same equation calculated by *uv2dv_cfd*):

*div = du_dx + dv_dy - term3*

With *uv2dv_cfd*, I'm using:

*divV = uv2dv_cfd( u850, v850, lat, lon, 3 )*

Below, a comparison between the results for latitude 90, for both
calculations:

(0) div manual      div function
(0) -0.435886     8.54325e-07
(1) 0.282385     8.54325e-07
(2) 0.641552     8.54325e-07
(3) 0.640924     8.54325e-07
(4) -3.11505     8.54325e-07
...
(138) 0.0877819     8.54325e-07
(139) 0.806053     8.54325e-07
(140) -2.59076     8.54325e-07
(141) -1.87243     8.54325e-07
(142) 2.96093     8.54325e-07
(143) -0.795052     8.54325e-07

And below, the comparison for latitude 87.5, for both calculations:

(0) div manual      div function
(0) 2.54982e-06     2.54933e-06
(1) 3.44865e-06     3.44799e-06
(2) 4.34874e-06     4.34791e-06
(3) 1.30406e-06     1.30381e-06
(4) 2.92225e-06     2.9217e-06
...
(138) 2.74009e-06     2.73957e-06
(139) 1.3016e-06     1.30136e-06
(140) 3.98753e-06     3.98678e-06
(141) 7.39285e-06     7.39144e-06
(142) 2.72861e-06     2.72809e-06
(143) 2.54953e-06     2.54905e-06

I've noted this problem only in the poles and seems to be related with *cos*
function in the computation of gradient of zonal wind over longitude, since
it doesn't give zero (or a very small value) for *cos(90)* or *cos(-90).*

I'm using NCEP Reanalysis 1 data and NCL 6.4.0.

Am I doing something wrong in these calculations?

Thanks,

Mateus
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190710/3774e6af/attachment.html>