anchor.diff:    real(kind=8), intent(OUT) :: y !< y in cartesian
     real(kind=8), intent(OUT) :: z !< z in cartesian
   
  -  x = cos(lat/180.d0*PI)*cos(lon/180.d0*PI)
  -  y = cos(lat/180.d0*PI)*sin(lon/180.d0*PI)
  -  z = sin(lat/180.d0*PI)
  +!!$d  x = cos(lat/180.d0*PI)*cos(lon/180.d0*PI)
  +!!$d  y = cos(lat/180.d0*PI)*sin(lon/180.d0*PI)
  +!!$d  z = sin(lat/180.d0*PI)
  +  x = cos(lat*D2R)*cos(lon*D2R)
  +  y = cos(lat*D2R)*sin(lon*D2R)
  +  z = sin(lat*D2R)
   
   end subroutine latlon2xyz_double
   
  
positive.code:   real(kind=8), intent(OUT) :: lat !< latitude [deg]
    real(kind=8), intent(OUT) :: lon !< longitude [deg]
  
    lat = atan2(z, sqrt(x*x+y*y))
    !lon = acos(x/cos(lat))
    lon = acos(x/sqrt(x*x+y*y))
  
negative.code:   lon = lon*180.d0/pi
  
    if (y<0) lon = 360.d0-lon
  
    if ( abs(x) < EPSLN ) then
      write(33,*)'xyz2latlon_double:',x,y,z,lon,lat
    end if
  
  end subroutine xyz2latlon_double
  
  