anchor.diff:    real(kind=8), intent(OUT) :: lat !< latitude [deg]
     real(kind=8), intent(OUT) :: lon !< longitude [deg]
   
  +  real(kind=8),parameter :: EPSLN=1.d-12
  +!
  +  if ( abs(x) < EPSLN .and. abs(y) < EPSLN ) then
  +    lon = 0.d0
  +    if ( z > 0.d0 ) then
  +      lat = 90.d0
  +    else
  +      lat = -90.d0
  +    end if
  +    return
  +  end if
  +
  +
     lat = atan2(z, sqrt(x*x+y*y))
     !lon = acos(x/cos(lat))
     lon = acos(x/sqrt(x*x+y*y))
  
positive.code:   lon = lon*180.d0/pi
  
    if (y<0) lon = 360.d0-lon
  
  end subroutine xyz2latlon_double
  
  
negative.code:   real(kind=8), intent(OUT) :: y !< y in cartesian
    real(kind=8), intent(OUT) :: z !< z in cartesian
  
  !!$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
  
  