      SUBROUTINE Euler_Richardson(x,y,vx,vy,t,C,g,dt,dt_2)
      IMPLICIT NONE
      REAL*8 ax, axmid, ay, aymid, v, v2, vmid, vmid2, vxmid, vymid,
     +       xmid, ymid, x, y, vx, vy, t, C, g, dt, dt_2
      v2 = vx*vx + vy*vy
      v = sqrt(v2)
      ax = -C*v*vx
      ay = -g - C*v*vy
      vxmid = vx + ax*dt_2                  ! velocity at midpoint
      vymid = vy + ay*dt_2
      xmid = x + vx*dt_2                    ! position at midpoint
      ymid = y + vy*dt_2
      vmid2 = vxmid*vxmid + vymid*vymid
      vmid = sqrt(vmid2)
      axmid = -C*vmid*vxmid                 ! acceleration at midpoint
      aymid = -g - C*vmid*vymid
      vx = vx + axmid*dt
      vy = vy + aymid*dt
      x = x + vxmid*dt
      y = y + vymid*dt
      t = t + dt
      END
