Re: possible bug with gc_latlon

This is a numerical rounding error that is occuring in
the underlying fortran code

       gcdist = acos( sin(rlat1r)*sin(rlat2r) +
      * cos(rlat1r)*cos(rlat2r)*cos(dlonr) )

so if
x = sin(rlat1r)*sin(rlat2r) +cos(rlat1r)*cos(rlat2r)*cos(dlonr)

ncl 0> x = 0.0
ncl 1> acx = acos(x)
ncl 2> print(acx) ; ==> 1.570796

ncl 4> x = 1.0
ncl 5> acx = acos(x)
ncl 6> print(acx) ; ==> 0.0

ncl 7> x = 1.0000001
ncl 8> acx = acos(x)
ncl 9> print(acx) ; ==> nan

On 02/11/2011 09:27 AM, Mary Haley wrote:
> Hi Matt,
>
> Thanks so much for providing the files and script up front.
>
> I ran this on my 64-bit Mac and a 64-bit LINUX box, and got "Nan" on both of them.
>
> On a 32-bit Linux box, I got 0.
>
> I've filed a ticket for this (NCL-995) and will look at it shortly.
>
> --Mary
>
> On Feb 9, 2011, at 6:15 PM, Matthew Higgins wrote:
>
>> Hi all,
>>
>> I think I may have found a platform-dependent bug with gc_latlon. I've
>> found that when my input latitude and longitudes are the identical,
>> rather than having a 0 returned (as in a distance of 0 km between the
>> two points), occasionally I will get an NaN returned. I only see this
>> issue if I'm running on a Cray XE6 - when I run the same code on my
>> local linux box, I get a 0. In both cases I am running NCL 5.2.1.
>>
>> Here is some example code. NCAR NCL folks - I have uploaded the two
>> .nc files referenced in the code below to the CGD ftp site if you'd
>> like to try to replicate what I am seeing.
>>
>> Matt
>>
>> Here's the script:
>>
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>> f1 = addfile("geo_em.d01.50km.nc", "r")
>> lat2d_50 = f1->XLAT_M(0,:,:)
>> lon2d_50 = f1->XLONG_M(0,:,:)
>> lat1d_50 = ndtooned(lat2d_50)
>> lon1d_50 = ndtooned(lon2d_50)
>> f2 = addfile("geo_em.d01.100km.nc", "r")
>> lat2d_100 = f2->XLAT_M(0,:,:)
>> lon2d_100 = f2->XLONG_M(0,:,:)
>> lat1d_100 = ndtooned(lat2d_100)
>> lon1d_100 = ndtooned(lon2d_100)
>>
>> i = 5832
>> q = 23532
>>
>> print("example lat/lon in 100km grid "+lat1d_100(i)+" "+lon1d_100(i))
>> print("example lat/lon in 50km grid "+lat1d_50(q)+" "+lon1d_50(q))
>> dist = gc_latlon(lat1d_100(i),lon1d_100(i),lat1d_50(q),lon1d_50(q),2,4)
>> print("distance in km between these points is "+dist)
>>
>>
>>
>> --
>> Matt Higgins
>> Cooperative Institute for Research in Environmental Sciences
>> University of Colorado - Boulder
>> (303) 735-1636
>> _______________________________________________
>> 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

-- 
======================================================
Dennis J. Shea                  tel: 303-497-1361    |
P.O. Box 3000                   fax: 303-497-1333    |
Climate Analysis Section                             |
Climate & Global Dynamics Div.                       |
National Center for Atmospheric Research             |
Boulder, CO  80307                                   |
USA                        email: shea 'at' ucar.edu |
======================================================