c **********************************************************************
      subroutine solvexf(v0,theta,ga,dt,xf,tf)
c **********************************************************************
      implicit real*8 (a-h,o-z)
c Constant
      pi    = 4.0d0*atan(1.0d0)
c
      t  = 0.0d0
      x  = 0.0d0
      y  = 0.0d0
      vx = v0*cos(theta*pi/180)
      vy = v0*sin(theta*pi/180)
c
 10   continue
        call onestep(x,y,vx,vy,t,dt,ga)
      if( y .ge. 0.0d0 ) goto 10
c
      call solvey0(x,y,vx,vy,t,ga)
      xf = x
      tf = t
c
      end
c **********************************************************************
      subroutine solvey0(x,y,vx,vy,t,ga)
c **********************************************************************
c   This subroutine gives the position where y = 0.
      implicit real*8 (a-h,o-z)
c
      eps = 1.0d-12
      do n = 1, 1000
        if(abs(y).lt.eps) goto 20
        dt = - y / vy
        call onestep(x,y,vx,vy,t,dt,ga)
      enddo
        write(*,*) " I cannot find the position where y = 0."
20    continue
c
      end