From: David Brown <dbrown_at_nyahnyahspammersnyahnyah>

Date: Thu Oct 01 2009 - 18:29:21 MDT

Date: Thu Oct 01 2009 - 18:29:21 MDT

Hi Mateus,

I think I have to disagree with Mary on this. I do not think it is a

bug, but a natural consequence of NCL's rules for the propagation of

missing values. The only way i think it can be considered a bug is in

that the initial choices of default missing values for NCL values were

a bit myopic. The default missing values should have been outside the

range of "normal" calculations. Likewise we should be encouraging

users to set their own missing values outside the range they expect

their computations to encompass, as Mary does here when she recommends

you use 1e20 rather than -32767 as the missing value. If you look at

this section of the NCL reference manual:

http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclExpressions.shtml#MissingValues

you will see that in the evaluation of any expression if any variable

has a _FillValue attribute it is propagated to the result of the

expression with elements on the left side of each operation prevailing

if there are different fill values. NCL has no way to know that in

this case an intermediate result is just accidently equal to the

_FillValue of the left hand variable in the calculation. And of course

once you have a missing value it is "viral" -- all subsequent

calculations are short-circuited and the missing value is propagated

as the result.

The bottom line is that this behavior results from a fundamental

property of NCL's data model that is built into all the arithmetic

operations. I don't think it could be changed without a radical

reworking of much significant code.

-dave

On Oct 1, 2009, at 3:17 PM, Mary Haley wrote:

*> Mateus,
*

*>
*

*> I've tracked down the source of the problem, and I think it's a bug
*

*> on our end, which I've reported.
*

*>
*

*> The issue is that when dh=2114, and the two v's equal -10.7 and -4.8,
*

*> then this part of the calcuation is equal to exactly -32767, which is
*

*> the missing value for "dh":
*

*>
*

*> dh * ( v(c,t,k-1,j,i) + v(c,t,k,j,i) )
*

*>
*

*> For some reason this makes the rest of the calculation missing.
*

*>
*

*> When you changed it so that it was multiplying by 0.5 first, then 0.5
*

*> got multiplied to dh before the rest of the calculation was done, so
*

*> you didn't end up with -32767.
*

*>
*

*> You've already found one work-around. Another one is to change the
*

*> missing value of "dh" to something other than -32767, like 1e20.
*

*>
*

*> I'll let you know if I hear about a fix, or whether this is a not a
*

*> bug.
*

*>
*

*> --Mary
*

*>
*

*> On Oct 1, 2009, at 10:18 AM, Mateus Teixeira wrote:
*

*>
*

*>> Hi Mary,
*

*>>
*

*>> Here are the results of your codes:
*

*>>
*

*>> ncl 0> su = 0.
*

*>> ncl 1> dh = 1271
*

*>> ncl 2> v1 = 0.699997
*

*>> ncl 3> v2 = -4.8
*

*>> ncl 4> su1 = su + dh * (v1 + v2) / 2
*

*>> ncl 5> su2 = su + 0.5 * dh * (v1 + v2)
*

*>> ncl 6> print("su1 = " + su1)
*

*>> (0) su1 = -2605.55
*

*>> ncl 7> print("su2 = " + su2)
*

*>> (0) su2 = -2605.55
*

*>>
*

*>> ncl 10> su = -2605.55
*

*>> ncl 11> dh = 2114
*

*>> ncl 12> v1 = -4.8
*

*>> ncl 13> v2 = -10.7
*

*>> ncl 14> su1 = su + dh * (v1 + v2) / 2
*

*>> ncl 15> su2 = su + 0.5 * dh * (v1 + v2)
*

*>> ncl 16> print("su1 = " + su1)
*

*>> (0) su1 = -18989.1
*

*>> ncl 17> print("su2 = " + su2)
*

*>> (0) su2 = -18989.1
*

*>>
*

*>> As I mentioned in the previous email, with both formulae the result
*

*>> is
*

*>> the same, when computed by hand, but in the loop it goes to a missing
*

*>> value.
*

*>>
*

*>> And here are the outputs from printVarSummary for all variables used
*

*>> in the formula:
*

*>>
*

*>>
*

*>> Variable: dh
*

*>> Type: float
*

*>> Total Size: 4 bytes
*

*>> 1 values
*

*>> Number of Dimensions: 1
*

*>> Dimensions and sizes: [1]
*

*>> Coordinates:
*

*>> Number Of Attributes: 1
*

*>> _FillValue : -32767
*

*>>
*

*>> Variable: sdh
*

*>> Type: float
*

*>> Total Size: 4 bytes
*

*>> 1 values
*

*>> Number of Dimensions: 1
*

*>> Dimensions and sizes: [1]
*

*>> Coordinates:
*

*>> Number Of Attributes: 1
*

*>> _FillValue : -32767
*

*>>
*

*>>
*

*>> Variable: sv
*

*>> Type: float
*

*>> Total Size: 4 bytes
*

*>> 1 values
*

*>> Number of Dimensions: 1
*

*>> Dimensions and sizes: [1]
*

*>> Coordinates:
*

*>> Number Of Attributes: 1
*

*>> _FillValue : -32767
*

*>>
*

*>>
*

*>> Variable: v (subsection)
*

*>> Type: float
*

*>> Total Size: 4 bytes
*

*>> 1 values
*

*>> Number of Dimensions: 1
*

*>> Dimensions and sizes: [1]
*

*>> Coordinates:
*

*>> Number Of Attributes: 8
*

*>> lon : 302.5
*

*>> lat : -35
*

*>> nivel : 250
*

*>> tempo : 1
*

*>> padrao : 1
*

*>> _FillValue : -32767
*

*>> long_name : vento meridional
*

*>> units : m/s
*

*>>
*

*>>
*

*>> Variable: hgt (subsection)
*

*>> Type: float
*

*>> Total Size: 4 bytes
*

*>> 1 values
*

*>> Number of Dimensions: 1
*

*>> Dimensions and sizes: [1]
*

*>> Coordinates:
*

*>> Number Of Attributes: 8
*

*>> lon : 302.5
*

*>> lat : -35
*

*>> nivel : 300
*

*>> tempo : 1
*

*>> padrao : 1
*

*>> _FillValue : -32767
*

*>> long_name : altura geopotencial
*

*>> units : m
*

*>>
*

*>>
*

*>> And here is the loop
*

*>>
*

*>> sdh = 0.
*

*>> su = 0.
*

*>> sv = 0.
*

*>> dh = 0.
*

*>> do k = k10000, k3000
*

*>> dh = hgt(c,t,k-1,j,i) - hgt(c,t,k,j,i)
*

*>> sdh = sdh + dh
*

*>> su = su + 0.5 * dh * ( u(c,t,k-1,j,i) + u(c,t,k,j,i) )
*

*>> sv = sv + 0.5 * dh * ( v(c,t,k-1,j,i) + v(c,t,k,j,i) )
*

*>> end do
*

*>>
*

*>> And all these lines are inside other loops for c, t, i, and j. I
*

*>> searched for missing values in the data (hgt, u, and v), founding
*

*>> nothing. It is very strange. I thought it was due a error
*

*>> accumulations, but the variables involved are equaled to zero
*

*>> before k
*

*>> loop.
*

*>>
*

*>> Best regards,
*

*>>
*

*>> Mateus
*

*>>
*

*>>
*

*>>
*

*>> 2009/10/1 Mary Haley <haley@ucar.edu>:
*

*>>> Hi Mateus,
*

*>>>
*

*>>> I was unable to reproduce your problem. Can you try this code and
*

*>>> let me
*

*>>> know what you get:
*

*>>>
*

*>>> su = 0.
*

*>>> dh = 1271
*

*>>> v1 = 0.699997
*

*>>> v2 = -4.8
*

*>>> su1 = su + dh * (v1 + v2) / 2
*

*>>> su2 = su + 0.5 * dh * (v1 + v2)
*

*>>> print("su1 = " + su1)
*

*>>> print("su2 = " + su2)
*

*>>>
*

*>>> su = -2605.55
*

*>>> dh = 2114
*

*>>> v1 = -4.8
*

*>>> v2 = -10.7
*

*>>> su1 = su + dh * (v1 + v2) / 2
*

*>>> su2 = su + 0.5 * dh * (v1 + v2)
*

*>>> print("su1 = " + su1)
*

*>>> print("su2 = " + su2)
*

*>>>
*

*>>> My first thought is that this has something to do with the types
*

*>>> of your
*

*>>> variables, and the order of which the arithmetic is being done
*

*>>> on them.
*

*>>>
*

*>>> Can you do a "printVarSummary" on all your variables and send me
*

*>>> the output
*

*>>> so I have the correct types?
*

*>>>
*

*>>> --Mary
*

*>>>
*

*>>>
*

*>>> On Oct 1, 2009, at 8:13 AM, Mateus Teixeira wrote:
*

*>>>
*

*>>>> Dear NCL users,
*

*>>>>
*

*>>>> I'm doing some algebric computations that are having strange
*

*>>>> behavior
*

*>>>> depending on way the formula is written. When written this way
*

*>>>>
*

*>>>> su = su + dh * ( u(c,t,k-1,j,i) + u(c,t,k,j,i) ) / 2
*

*>>>>
*

*>>>> I'm getting a missing value (the formula is inside a do loop), as
*

*>>>> shown
*

*>>>> below:
*

*>>>>
*

*>>>> (0) sv(before)=0; dh=1271; v(k-1)=0.699997; v(k)=-4.8
*

*>>>> (0) sv(after)=-2605.55
*

*>>>> (0) sv(before)=-2605.55; dh=2114; v(k-1)=-4.8; v(k)=-10.7
*

*>>>> (0) sv(after)=-32767 <= MISSING VALUE
*

*>>>>
*

*>>>> But when the formula is written as
*

*>>>>
*

*>>>> su = su + 0.5 * dh * ( u(c,t,k-1,j,i) + u(c,t,k,j,i) )
*

*>>>>
*

*>>>> I haven't a missingl value, as shown below for the same stop
*

*>>>> point showed
*

*>>>> above
*

*>>>>
*

*>>>> (0) sv(before)=0; dh=1271; v(k-1)=0.699997; v(k)=-4.8
*

*>>>> (0) sv(after)=-2605.55
*

*>>>> (0) sv(before)=-2605.55; dh=2114; v(k-1)=-4.8; v(k)=-10.7
*

*>>>> (0) sv(after)=-18989.1 <=CORRECT VALUE
*

*>>>>
*

*>>>> I made these same computations by hand in NCL with both formulae
*

*>>>> and
*

*>>>> get the correct value given by the second way of the formula.
*

*>>>>
*

*>>>> Can anyone imagine what is causing this? I thought that the float
*

*>>>> representation of the 1/2 could be the source of the error, but the
*

*>>>> difference between missing value and the correct value is not
*

*>>>> little.
*

*>>>>
*

*>>>> Best regards,
*

*>>>>
*

*>>>> --
*

*>>>> Mateus da Silva Teixeira
*

*>>>> Meteorologista
*

*>>>> Instituto de Pesquisas Meteorológicas - IPMet
*

*>>>> Universidade Estadual Paulista - UNESP
*

*>>>> Av. Luis Edmundo Carrijo Coube, 14-01 - Cx. Postal 281 - CEP
*

*>>>> 17033-360
*

*>>>> Bauru - SP - Brasil
*

*>>>> Fone: +55 14 3103-6030 / fax: 3203-3649
*

*>>>>
*

*>>>> Registered Linux User #466740 (http://counter.li.org/)
*

*>>>> _______________________________________________
*

*>>>> ncl-talk mailing list
*

*>>>> List instructions, subscriber options, unsubscribe:
*

*>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
*

*>>>
*

*>>>
*

*>>
*

*>>
*

*>>
*

*>> --
*

*>> Mateus da Silva Teixeira
*

*>> Meteorologista
*

*>> Instituto de Pesquisas Meteorológicas - IPMet
*

*>> Universidade Estadual Paulista - UNESP
*

*>> Av. Luis Edmundo Carrijo Coube, 14-01 - Cx. Postal 281 - CEP
*

*>> 17033-360
*

*>> Bauru - SP - Brasil
*

*>> Fone: +55 14 3103-6030 / fax: 3203-3649
*

*>>
*

*>> Registered Linux User #466740 (http://counter.li.org/)
*

*>> _______________________________________________
*

*>> ncl-talk mailing list
*

*>> List instructions, subscriber options, unsubscribe:
*

*>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
*

*>
*

*> _______________________________________________
*

*> ncl-talk mailing list
*

*> List instructions, subscriber options, unsubscribe:
*

*> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
*

_______________________________________________

ncl-talk mailing list

List instructions, subscriber options, unsubscribe:

http://mailman.ucar.edu/mailman/listinfo/ncl-talk

Received on Thu Oct 1 18:29:26 2009

*
This archive was generated by hypermail 2.1.8
: Mon Oct 26 2009 - 15:01:33 MDT
*