module piab
    implicit none
    real(8),parameter::pi=3.14159265358979323846D0
    real(8),parameter::a=6378137.000000D0 !Equatorial radius
    real(8),parameter::b=6356752.314140D0 !Polar radius
end module piab


program f90_DIST1
    use piab
    implicit none
    character(len=100)::strcom
    integer::ndata
    character(len=50),allocatable::locname(:)
    real(8),allocatable::xlon(:),ylat(:),dist(:,:)
    real(8)::d1,d2
    real(8)::lat1,lon1
    real(8)::lat2,lon2
    character(len=50)::fnameR,fnameW
    integer::i,j
    character(len=100)::dummy
    character(len=1000)::linebuf

    call getarg(1,fnameR)
    call getarg(2,fnameW)

    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,*) ndata
    allocate(locname(1:ndata))
    allocate(xlon(1:ndata),ylat(1:ndata))
    allocate(dist(1:ndata,1:ndata))
    do i=1,ndata
        read(11,*) locname(i),xlon(i),ylat(i)
    end do
    close(11)

    do j=1,ndata
        do i=1,ndata
            dist(i,j)=0.0D0
        end do
    end do
    do j=2,ndata
        do i=1,j-1
            dist(i,j)=DLA(xlon(i),ylat(i),xlon(j),ylat(j))
        end do
    end do

    open(12,file=fnameW,status='replace')
    write(12,'(a)') trim(adjustl(strcom))
    write(12,'(i5)') ndata
    do j=1,ndata
        write(dummy,*) locname(j)
        linebuf=trim(dummy)
        do i=1,j
            write(dummy,'(f10.3)') dist(i,j)/1000.0D0
            linebuf=trim(linebuf)//','//trim(dummy)
        end do
        call del_spaces(linebuf)
        write(12,'(a)') trim(adjustl(linebuf))
    end do
    close(12)

stop

contains

    real(8) function DLA(lon1,lat1,lon2,lat2)
        !Lambert-Andyer's formula
        real(8),intent(in)::lat1,lon1,lat2,lon2
        real(8)::IA,LA
        real(8)::IB,LB
        real(8)::phiA,phiB,XX,delta
        real(8)::rho
        real(8)::FF

        IA=lat1/180.0*pi
        LA=lon1/180.0*pi
        IB=lat2/180.0*pi
        LB=lon2/180.0*pi

        FF=(a-b)/a
        phiA=atan((1.0-FF)*tan(IA))
        phiB=atan((1.0-FF)*tan(IB))
        XX=acos(sin(phiA)*sin(phiB)+cos(phiA)*cos(phiB)*cos(LA-LB))
        delta=FF/8.0*&
        (&
         (sin(XX)-XX)*(sin(phiA)+sin(phiB))**2.0D0/(cos(0.5D0*XX))**2.0D0&
        -(sin(XX)+XX)*(sin(phiA)-sin(phiB))**2.0D0/(sin(0.5D0*XX))**2.0D0&
        )
        DLA=a*(XX+delta)
    end function DLA
    !**************************************************
    ! Deleate space
    !**************************************************
    subroutine del_spaces(s)
        character (*), intent (inout) :: s
        character (len=len(s)) tmp
        integer i, j
        j = 1
        do i = 1, len(s)
            if (s(i:i)==' ') cycle
            tmp(j:j) = s(i:i)
            j = j + 1
        end do
        s = tmp(1:j-1)
    end subroutine del_spaces

end program f90_DIST1