subroutine direct_vincenty(a,f,glat1,glon1,faz,s,glat2,glon2,baz)

    implicit none

    real(wp),intent(in)  :: a        !! semimajor axis of ellipsoid [m]
    real(wp),intent(in)  :: f        !! flattening of ellipsoid [-]
    real(wp),intent(in)  :: glat1    !! latitude of 1 [rad]
    real(wp),intent(in)  :: glon1    !! longitude of 1 [rad]
    real(wp),intent(in)  :: faz      !! forward azimuth 1->2 [rad]
    real(wp),intent(in)  :: s        !! distance from 1->2 [m]
    real(wp),intent(out) :: glat2    !! latitude of 2 [rad]
    real(wp),intent(out) :: glon2    !! longitude of 2 [rad]
    real(wp),intent(out) :: baz      !! back azimuth 2->1 [rad]

    real(wp) :: r,tu,sf,cf,cu,su,sa,c2a,x,c,d,y,sy,cy,cz,e

    real(wp),parameter :: eps = 0.5e-13_wp

    r  = 1.0_wp - f
    tu = r*sin(glat1)/cos(glat1)
    sf = sin(faz)
    cf = cos(faz)
    if ( cf/=0.0_wp ) then
        baz = atan2(tu,cf)*2.0_wp
    else
        baz = 0.0_wp
    end if
    cu  = 1.0_wp/sqrt(tu*tu+1.0_wp)
    su  = tu*cu
    sa  = cu*sf
    c2a = -sa*sa + 1.0_wp
    x   = sqrt((1.0_wp/r/r-1.0_wp)*c2a+1.0_wp) + 1.0_wp
    x   = (x-2.0_wp)/x
    c   = 1.0_wp - x
    c   = (x*x/4.0_wp+1.0_wp)/c
    d   = (0.375_wp*x*x-1.0_wp)*x
    tu  = s/r/a/c
    y   = tu
    do
        sy = sin(y)
        cy = cos(y)
        cz = cos(baz+y)
        e  = cz*cz*2.0_wp - 1.0_wp
        c  = y
        x  = e*cy
        y  = e + e - 1.0_wp
        y  = (((sy*sy*4.0_wp-3.0_wp)*y*cz*d/6.0_wp+x)*d/4.0_wp-cz)*sy*d + tu
        if ( abs(y-c)<=eps ) exit
    end do
    baz   = cu*cy*cf - su*sy
    c     = r*sqrt(sa*sa+baz*baz)
    d     = su*cy + cu*sy*cf
    glat2 = atan2(d,c)
    c     = cu*cy - su*sy*cf
    x     = atan2(sy*sf,c)
    c     = ((-3.0_wp*c2a+4.0_wp)*f+4.0_wp)*c2a*f/16.0_wp
    d     = ((e*cy*c+cz)*sy*c+y)*sa
    glon2 = glon1 + x - (1.0_wp-c)*d*f
    baz   = atan2(sa,baz) + pi

    end subroutine direct_vincenty