module m_spherdist
contains
real function spherdist(lon1,lat1,lon2,lat2)
! --- -----------------------------------------
! --- Computes the distance between geo. pos.
! --- lon1,lat1 and lon2,lat2. 
! --- INPUT is in degrees.
! --- -----------------------------------------

   implicit none
   REAL, intent(in) :: lon1,lat1,lon2,lat2 ! Pos. in degrees

   ! KAL -- The routines below lead to lack of precision
   !     -- so we hardcode reals to real*8 in temp calcs
   real*8, parameter :: invradian=0.017453292
   real*8, parameter :: rearth=6371001.0     ! Radius of earth

   real*8  rlon1,rlat1,rlon2,rlat2           ! Pos. in radians
   real*8  x1,y1,z1,x2,y2,z2                 ! Cartesian position
   real*8  dx,dy,dz,dr,dott                  ! Cartesian distances


   rlon1=lon1*invradian             !lon1 in rad
   rlat1=(90.-lat1)*invradian       !90-lat1 in rad 

   rlon2=lon2*invradian             !lon2 in rad
   rlat2=(90.-lat2)*invradian       !90-lat2 in rad 

   x1= SIN(rlat1)*COS(rlon1)        !x,y,z of pos 1.
   y1= SIN(rlat1)*SIN(rlon1)
   z1= COS(rlat1) 

   x2= SIN(rlat2)*COS(rlon2)        !x,y,z of pos 2.
   y2= SIN(rlat2)*SIN(rlon2)
   z2= COS(rlat2) 

   dx=x2-x1                         !distances in x, y, z 
   dy=y2-y1
   dz=z2-z1

   !print '(9f14.8,e14.5)',z1-z2,x1-x2,y1-y2,x1,x2,y1,y2,z1,z2,&
   !x1*x2+y1*y2+z1*z2
   dott=min(x1*x2+y1*y2+z1*z2,1.)
   !dr=SQRT(dx*dx+dy*dy+dz*dz)       !distance pytagaros
   !dr=acos(x1*x2+y1*y2+z1*z2)       ! Arc length
   dr=acos(dott)       ! Arc length

   spherdist=dr*rearth

end function spherdist
end module m_spherdist