56 real :: rlat1,rlon1,rlat2,rlon2,
distance
57 real(kind=dp) :: clat1,clat2,slat1,slat2,cdlon,crd
58 real(kind=dp),
parameter :: rerth=6.3712e6_dp
59 real(kind=dp),
parameter :: pi=3.14159265358979_dp, dpr=180.0_dp/pi
61 if ((abs(rlat1-rlat2).lt.0.03).and. &
62 (abs(rlon1-rlon2).lt.0.03))
then
65 clat1=cos(
real(rlat1,kind=dp)/dpr)
66 slat1=sin(
real(rlat1,kind=dp)/dpr)
67 clat2=cos(
real(rlat2,kind=dp)/dpr)
68 slat2=sin(
real(rlat2,kind=dp)/dpr)
69 cdlon=cos(
real((rlon1-rlon2),kind=dp)/dpr)
70 crd=slat1*slat2+clat1*clat2*cdlon
71 distance=
real(rerth*acos(crd)/1000.0_dp)
real function distance(rlat1, rlon1, rlat2, rlon2)