inverse_vincenty Subroutine

    subroutine inverse_vincenty(a,rf,b1,l1,b2,l2,faz,baz,s,it,sig,lam,kind)

    implicit none

    real(wp),intent(in)  :: a     !! Equatorial semimajor axis
    real(wp),intent(in)  :: rf    !! reciprocal flattening (1/f)
    real(wp),intent(in)  :: b1    !! latitude of point 1 (rad, positive north)
    real(wp),intent(in)  :: l1    !! longitude of point 1 (rad, positive east)
    real(wp),intent(in)  :: b2    !! latitude of point 2 (rad, positive north)
    real(wp),intent(in)  :: l2    !! longitude of point 2 (rad, positive east)
    real(wp),intent(out) :: faz   !! Forward azimuth (rad, clockwise from north)
    real(wp),intent(out) :: baz   !! Back azimuth (rad, clockwise from north)
    real(wp),intent(out) :: s     !! Ellipsoidal distance
    integer,intent(out)  :: it    !! iteration count
    real(wp),intent(out) :: sig   !! spherical distance on auxiliary sphere
    real(wp),intent(out) :: lam   !! longitude difference on auxiliary sphere
    integer,intent(out)  :: kind  !! solution flag: kind=1, long-line; kind=2, antipodal

    real(wp) :: beta1,beta2,biga,bigb,bige,bigf,boa,c,cosal2,coslam,&
                cossig,costm,costm2,cosu1,cosu2,d,dsig,ep2,l,prev,&
                sinal,sinlam,sinsig,sinu1,sinu2,tem1,tem2,temp,test,z

    real(wp),parameter :: tol = 1.0e-14_wp   !! convergence tolerance
    real(wp),parameter :: eps = 1.0e-15_wp   !! tolerance for zero

    boa = 1.0_wp - 1.0_wp/rf   ! b/a

    beta1 = atan(boa*tan(b1))  ! better reduced latitude
    sinu1 = sin(beta1)
    cosu1 = cos(beta1)
    beta2 = atan(boa*tan(b2))  ! better reduced latitude
    sinu2 = sin(beta2)
    cosu2 = cos(beta2)

    l = l2 - l1  ! longitude difference [-pi,pi]
    if ( l>pi ) l = l - twopi
    if ( l<-pi ) l = l + twopi
    prev = l
    test = l
    it   = 0
    kind = 1
    lam  = l

    longline : do  ! long-line loop (kind=1)

        sinlam = sin(lam)
        coslam = cos(lam)
        temp   = cosu1*sinu2 - sinu1*cosu2*coslam
        sinsig = sqrt((cosu2*sinlam)**2+temp**2)
        cossig = sinu1*sinu2 + cosu1*cosu2*coslam
        sig    = atan2(sinsig,cossig)
        if ( abs(sinsig)<eps ) then
            sinal = cosu1*cosu2*sinlam/sign(eps,sinsig)
        else
            sinal = cosu1*cosu2*sinlam/sinsig
        endif
        cosal2 = -sinal**2 + 1.0_wp
        if ( abs(cosal2)<eps ) then
            costm = -2.0_wp*(sinu1*sinu2/sign(eps,cosal2)) + cossig
        else
            costm = -2.0_wp*(sinu1*sinu2/cosal2) + cossig
        endif
        costm2 = costm*costm
        c = ((-3.0_wp*cosal2+4.0_wp)/rf+4.0_wp)*cosal2/rf/16.0_wp

        antipodal : do  ! antipodal loop (kind=2)

            it = it + 1
            d = (((2.0_wp*costm2-1.0_wp)*cossig*c+costm)*sinsig*c+sig)*(1.0_wp-c)/rf
            if ( kind==1 ) then
                lam = l + d*sinal
                if ( abs(lam-test)>=tol ) then
                    if ( abs(lam)>pi ) then
                        kind   = 2
                        lam    = pi
                        if ( l<0.0_wp ) lam = -lam
                        sinal  = 0.0_wp
                        cosal2 = 1.0_wp
                        test   = 2.0_wp
                        prev   = test
                        sig    = pi - abs(atan(sinu1/cosu1)+atan(sinu2/cosu2))
                        sinsig = sin(sig)
                        cossig = cos(sig)
                        c      = ((-3.0_wp*cosal2+4.0_wp)/rf+4.0_wp)*cosal2/rf/16.0_wp
                        if ( abs(sinal-prev)<tol ) exit longline
                        if ( abs(cosal2)<eps ) then
                            costm = -2.0_wp*(sinu1*sinu2/sign(eps,cosal2)) + cossig
                        else
                            costm = -2.0_wp*(sinu1*sinu2/cosal2) + cossig
                        endif
                        costm2 = costm*costm
                        cycle antipodal
                    endif
                    if ( ((lam-test)*(test-prev))<0.0_wp .and. it>5 ) &
                         lam = (2.0_wp*lam+3.0_wp*test+prev)/6.0_wp  ! refined converge.
                    prev = test
                    test = lam
                    cycle longline
                endif
            else
                sinal  = (lam-l)/d
                if ( ((sinal-test)*(test-prev))<0.0_wp .and. it>5 ) &
                       sinal = (2.0_wp*sinal+3.0_wp*test+prev)/6.0_wp  ! refined converge.
                prev   = test
                test   = sinal
                cosal2 = -sinal**2 + 1.0_wp
                sinlam = sinal*sinsig/(cosu1*cosu2)
                coslam = -sqrt(abs(-sinlam**2+1.0_wp))
                lam    = atan2(sinlam,coslam)
                temp   = cosu1*sinu2 - sinu1*cosu2*coslam
                sinsig = sqrt((cosu2*sinlam)**2+temp**2)
                cossig = sinu1*sinu2 + cosu1*cosu2*coslam
                sig    = atan2(sinsig,cossig)
                c      = ((-3.0_wp*cosal2+4.0_wp)/rf+4.0_wp)*cosal2/rf/16.0_wp
                if ( abs(sinal-prev)>=tol ) then
                    if ( abs(cosal2)<eps ) then
                        costm = -2.0_wp*(sinu1*sinu2/sign(eps,cosal2)) + cossig
                    else
                        costm = -2.0_wp*(sinu1*sinu2/cosal2) + cossig
                    endif
                    costm2 = costm*costm
                    cycle antipodal
                endif
            endif

            exit longline  !finished

        end do antipodal

    end do longline

    ! convergence

    if ( kind==2 ) then  ! antipodal
        faz  = sinal/cosu1
        baz  = sqrt(-faz**2+1.0_wp)
        if ( temp<0.0_wp ) baz = -baz
        faz  = atan2(faz,baz)
        tem1 = -sinal
        tem2 = sinu1*sinsig - cosu1*cossig*baz
        baz  = atan2(tem1,tem2)
    else  ! long-line
        tem1 = cosu2*sinlam
        tem2 = cosu1*sinu2 - sinu1*cosu2*coslam
        faz  = atan2(tem1,tem2)
        tem1 = -cosu1*sinlam
        tem2 = sinu1*cosu2 - cosu1*sinu2*coslam
        baz  = atan2(tem1,tem2)
    endif
    if ( faz<0.0_wp ) faz = faz + twopi
    if ( baz<0.0_wp ) baz = baz + twopi

    ! helmert 1880 from vincenty "geodetic inverse solution between antipodal points"

    ep2  = 1.0_wp/(boa*boa) - 1.0_wp
    bige = sqrt(1.0_wp+ep2*cosal2)
    bigf = (bige-1.0_wp)/(bige+1.0_wp)
    biga = (1.0_wp+bigf*bigf/4.0_wp)/(1.0_wp-bigf)
    bigb = bigf*(1.0_wp-0.375_wp*bigf*bigf)
    z    = bigb/6.0_wp*costm*(-3.0_wp+4.0_wp*sinsig**2)*(-3.0_wp+4.0_wp*costm2)
    dsig = bigb*sinsig*(costm+bigb/4.0_wp*(cossig*(-1.0_wp+2.0_wp*costm2)-z))
    s    = (boa*a)*biga*(sig-dsig)

    end subroutine inverse_vincenty