C  echo 'this is not compilable/executable source code see:  nb2.f'
C
C----ROB SALGADO (cps615g)
C prob2.f (rk.f) 10/12/92 (PROBLEM 4.2 -- SEE PROBLEM 5.1 FOR REVISION)
C
C RUNGE-KUTTA (Prob 2)
C
C ---------------------------------------IGNORABLE:  ORIGINAL COMMENTS
C                              
C
C I compiled the program with the mpf alias
C        mpf -nofecommon nb2.f
C
C This is the result of my attempt to run it.
C I have to admit, I don't know what $ Print ctl> is asking for.
C I'll try this again.
C But for now,  here is the output.
C The subroutine follows.
C 
C
C :/\/\/\/\ 
C
C magpie% a.out
C $ Data file>
C f7b16.dat
C $ Print ctl>
C 2
C FIO999 ENDFIL end of file.
C
C MPPE : intercepted signal SIGILL(4) @ FE pc 0x42b170
C on behalf of /u2/cps615g/a.out
C MPPE : You have the DPU
C  ID            Function symbol SRC  SIDE   STAT     FE PC |   ACU PC
C   0               NB2MODE(480):333, BOTH, UFUNC 0x0403f9c | 0x006684
C No MPPEs registered.  Try again? (Y/N)[N]: n
C
C :/\/\/\/ I DON'T UNDERSTAND WHAT THIS MEANS
C --------------------------------------------------------------------


C
C
C------- modeul(h, ns)  --------------------------------

C  Perform ns steps of modified Euler with time-step h.
C  X0, V0 are taken to be current state, and are updated
C  with final state.  X1, V1 are used in the process.
C  Note that X's, V's are calculated redundantly.

      subroutine modeul(h, ns)
      real h
      integer ns

      INCLUDE 'nb2mode.h'

      INTERFACE
        function Grav(Xf)
        INCLUDE 'nbmode.h'
        real, array(1:N, 1:N, 1:3) :: Xf, Grav
CMPF    ONDPU Xf
      END INTERFACE

C----------------------FOR NOW, JUST BETWEEN YOU AND ME,
C                      WE WILL USE THE 4-TH ORDER RUNGE-KUTTA
C                      IF YOU CAN READ THIS MESSAGE, THEN THAT MEANS
C                      I HAVEN'T RENAMED modeul() TO, SAY, rk4()

      real h2
      integer k

      real, array(1:N,1:N,1:3) :: k1,k2,k3,k4

      h2 = 0.5*h        ! h/2, for convenience
      do k = 1, ns      ! ns steps
D     print *, " --- X0  -----------"
D     call print_x(X0, 7, 1)    ! For testing
      k1=h*V0                                    ! X'=V
      k2=h*(V0+k1/2.0)
      k3=h*(V0+k2/2.0)
      k4=h*(V0+k3)
      X1=X0+(k1+k4+2.0*(k2+k3))/6.0              ! this will be the new X

      k1=h*Grav(X0)                              ! V'=Grav(X)
      k2=h*(Grav(X0)+k1/2.0)
      k3=h*(Grav(X0)+k2/2.0)
      k4=h*(Grav(X0)+k3)
      V0=V0+(k1+k4+2.0*(k2+k3))/6.0              ! this is the new V

      X0=X1                                      ! this is the new X


C       (Modified Euler)
C       X1 = X0 + h2*V0         ! Half-way positions
C       V1 = V0 + h2*Grav(X0)   ! Half-way velocities
C       X0 = X0 + h*V1          ! Next positions
C       V0 = V0 + h*Grav(X1)    ! Next velocities

      end do
      end subroutine modeul
 
C----------------
