!  Adams-Bashforth/Adams-Moulton predictor/corrector for use with
!  nbmstep.f90.

!  Advance ns steps of h.  Current state is in X0, V0, which will be
!  updated to new state.  History is maintained in VH, AH, HI

      subroutine Multistep(h, ns)
      USE nbmsmod
      USE nbrmod
      real h
      integer ns

      INTERFACE
	subroutine Grav(X)
        USE nbrmod
        real, dimension(1:N, 1:3) :: X
	end subroutine Grav
       END INTERFACE

      real, dimension(1:N, 1:3) :: A0
      real h24
      integer u, ti

      h24 = h/24.0
      do, u = 1, ns
	call Grav(X0);  A0 = Ac
!  Predictor
        X1 = X0 + h24*(55.0*V0 - 59.0*VH(HI(1),:,:)  &
     			       + 37.0*VH(HI(2),:,:)  &
     			       -  9.0*VH(HI(3),:,:))
	V1 = V0 + h24*(55.0*A0 - 59.0*AH(HI(1),:,:)  &
     			       + 37.0*AH(HI(2),:,:)  &
     			       -  9.0*AH(HI(3),:,:))
	VH(HI(3),:,:) = V0				! For future steps
	AH(HI(3),:,:) = A0
!  Corrector
	X0 = X0 + h24*(9.0*V1 + 19.0*V0 - 5.0*VH(HI(1),:,:)  &
                                        +     VH(HI(2),:,:))
        call Grav(X1)
	V0 = V0 + h24*(9.0*Ac + 19.0*A0 - 5.0*AH(HI(1),:,:)  &
     					+     AH(HI(2),:,:))
!  Adjust indices
        ti    = HI(3)
	HI(3) = HI(2)
	HI(2) = HI(1)
	HI(1) = ti
      end do
      end subroutine Multistep
