[ncl-talk] Problem with divergence calculation
Dennis Shea
shea at ucar.edu
Mon Jul 15 09:57:43 MDT 2019
Oi Mateus,
Your data is global. Hence, to derive divergence you can use either CFD or,
better, SPH (Spherical Harmonics) to derive divergence. The latter are more
accurate. However, visually, it is hard notice the difference. [See
attached].
See Example 5 @* http://www.ncl.ucar.edu/Applications/wind.shtml*
<http://www.ncl.ucar.edu/Applications/wind.shtml>
Example 6 is a more 'advanced' application.
If you want to plot the divergence over a specified region, use the
*{latS:latN},{lon{W:lonE}
*syntax.
---
That said ..... the [edited] output from the attached script indicates some
very subtle issue for the CFD gradients in the 'y' [latitude] direction.
The same base function is used ti calculate the gradients
Variable: divs
Dimensions and sizes: [lat | 73] x [lon | 144]
Coordinates:
lat: [-90..90]
lon: [ 0..357.5]
===============
(0) divergence (1/s) : min=-2.68841e-05 max=3.66626e-05
(0) ---
(0) du/dx , du/dlon (1/s) : min=-4.13793e-05 max=5.20984e-05 <==
spherical harmonics
(0) du/dy , du/dlat (1/s) : min=-7.50804e-05 max=6.49166e-05
(0) ---
(0) dv/dx , dv/dlon (1/s) : min=-6.04759e-05 max=6.54951e-05 ,==
spherical harmonics
(0) dv/dy , dv/dlat (1/s) : min=-4.19307e-05 max=4.15874e-05
(0) ---
(0) divf: min=-2.48427e-05 max=3.53151e-05 <== CFD
(0) dudx: cfd (1/s) : min=-6.51086e-05 max=5.64755e-05 <=== good
(*0) dudy: cfd (1/s) : min=-32.9175 max=32.9176 <==================
???*
(0) ---
(0) dvdx: cfd (1/s) : min=-3.72306e-05 max=2.75183e-05 <=== good
*(0) dvdy: cfd (1/s) : min=-32.9176 max=32.9176 <==================
???*
*HTH*
*D*
On Wed, Jul 10, 2019 at 1:03 PM Mateus da Silva Teixeira <
mateusstex at gmail.com> wrote:
> Hi Dennis,
>
> Thank you for your answer, but this didn't make difference in the
> calculations. As I could read in documentation, there is no warning about
> latitude order.
>
> Trying more comparisons, I've calculated the wind gradients with
> *grad_latlon_cfd* function to get *du_dx* and *dv_dy* and the results are
> the same, as shown below:
>
> *For latitude -90*
> (0) du_dx manual du_dx grad funct
> (0) -32.9238 -32.9175
> (1) -32.9238 -32.9175
> (2) -28.8081 -28.8026
> (3) -28.8081 -28.8026
>
> (0) dv_dy manual dv_dy grad funct
> (0) 5.39676e-06 5.39574e-06
> (1) 5.39676e-06 5.39574e-06
> (2) 5.39676e-06 5.39574e-06
> (3) 5.39676e-06 5.39574e-06
>
> *For latitude -87.5*
> (0) du_dx manual du_dx grad funct
> (0) 2.88688e-05 2.88633e-05
> (1) 2.88688e-05 2.88633e-05
> (2) 2.47451e-05 2.47403e-05
> (3) 2.88688e-05 2.88633e-05
>
> (0) dv_dy manual dv_dy grad funct
> (0) 2.69839e-06 2.69787e-06
> (1) 2.51851e-06 2.51803e-06
> (2) 2.33858e-06 2.33813e-06
> (3) 2.33861e-06 2.33816e-06
>
> And below the divergence, but along with the calculation using gradients
> obtained with *grad_latlon_cfd*:
>
> *For latitude -90*
> (0) div manual div function div grad funct
> (0) -0.24189 5.72546e-06 -0.235588
> (1) -0.960165 5.72546e-06 -0.953863
> (2) 2.43722 5.72546e-06 2.44274
> (3) 1.71895 5.72546e-06 1.72446
>
> *For latitude -87.5*
> (0) div manual div function div grad funct
> (0) 4.24092e-06 4.24011e-06 4.23487e-06
> (1) 4.78015e-06 4.77923e-06 4.77414e-06
> (2) 1.19561e-06 1.19538e-06 1.19042e-06
> (3) 6.0385e-06 6.03735e-06 6.03253e-06
>
> The data used for these calculations has its latitude ordered
> south-to-north. I think that I'm missing something, but I can't figure it
> out.
>
> Best regards,
>
> Mateus
>
>
> Em qua, 10 de jul de 2019 às 13:23, Dennis Shea <shea at ucar.edu> escreveu:
>
>> NCEP reanalyses are ordered North-to-South. The function requires
>> South-to-North. Assuming (time,lev,Lat,Lon)
>>
>> V = V(:,:,::-1,:)
>> printVarSummary(V)
>>
>> Look at the lat coordinate.
>>
>> Sent from my iPhone
>>
>> On Jul 10, 2019, at 8:33 AM, Mateus da Silva Teixeira via ncl-talk <
>> ncl-talk at ucar.edu> wrote:
>>
>> 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
>>
>>
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk at ucar.edu
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190715/475475f4/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: div.ncl_talk.ncl
Type: application/octet-stream
Size: 5222 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190715/475475f4/attachment-0001.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mateus.png
Type: image/png
Size: 370678 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20190715/475475f4/attachment-0001.png>
More information about the ncl-talk
mailing list