!_______________________________________________________________________
! This module contains IVP integration routines
! from the FORTRAN 77 library RKSUITE
! Aleksandar Donev, November 2000
!_______________________________________________________________________
MODULE ODE
      ! use RKSUITE
      IMPLICIT NONE
      PUBLIC
!
      INTEGER, PARAMETER :: wp = KIND (0.0D0)
      INTEGER, SAVE :: flag = 1
      EXTERNAL :: SETUP, CT, INTRP
!
CONTAINS
      SUBROUTINE ODESolve (f, n_eq, x_range, n_points, y0, x, y, relerr, abserr)
         INTEGER, INTENT (IN) :: n_eq
         REAL (KIND=wp), DIMENSION (2), INTENT (IN) :: x_range
         INTEGER, INTENT (IN) :: n_points
         REAL (KIND=wp), DIMENSION (n_eq), INTENT (IN) :: y0
         REAL (KIND=wp), DIMENSION (n_eq), INTENT (OUT), OPTIONAL :: x
         REAL (KIND=wp), DIMENSION (n_eq, n_points) :: y
         REAL (KIND=wp), INTENT (IN) :: relerr
         REAL (KIND=wp), DIMENSION (n_eq), INTENT (IN), OPTIONAL :: abserr
         INTERFACE
            SUBROUTINE f (x, y, dy)
               INTEGER, PARAMETER :: wp = KIND (0.0D0)
               REAL (KIND=wp), INTENT (IN) :: x
               REAL (KIND=wp), DIMENSION (*), INTENT (IN) :: y
               REAL (KIND=wp), DIMENSION (*), INTENT (OUT) :: dy
            END SUBROUTINE f
         END INTERFACE
!
         REAL (KIND=wp), DIMENSION (n_eq) :: tolerr, y_now, dy_now, dy, y_out ! y_got, , dy_got, dy_max
         REAL (KIND=wp), DIMENSION (32*n_eq) :: workspace
         REAL (KIND=wp), DIMENSION (6*n_eq) :: workint
         REAL (KIND=wp) :: dx, x_want, x_now ! , x_got
         INTEGER :: i
!
         IF (present(abserr)) THEN
            tolerr = abserr
         ELSE
            tolerr = 0.01_wp * relerr
         END IF
!
         ! call SETUP(NEQ=n_eq,TSTART=x_range(1),YSTART=y0,TEND=x_range(2),TOL=relerr, &
         !            THRES=tolerr,method=2,task='C',ERRASS=.FALSE.,HSTART=0.0D0, &
         !            WORK=workspace,LENWRK=32*n_eq,MESAGE=.TRUE.)
         CALL SETUP (n_eq, x_range(1), y0, x_range(2)+(x_range(2)-x_range(1))/10, relerr, tolerr, 2, 'C', .FALSE., 0.0D0, &
                     & workspace, 32*n_eq, .TRUE.)
!
         dx = (x_range(2)-x_range(1)) / (n_points-1)
         x_want = x_range (1)
         x_now = x_range (1)
         i = 1
         MainLoop: DO
            TakeSteps: DO
               ! call CT(F=f,TNOW=x_now,YNOW=y_now,YPNOW=dy_now, &
               !        WORK=workspace,CFLAG=flag)
               IF (x_now > x_want) EXIT TakeSteps
               CALL CT (f, x_now, y_now, dy_now, workspace, flag)
            END DO TakeSteps
            InterpResults: DO
               IF (present(x)) x (i) = x_want
               ! call INTRP(TWANT=x_want,REQEST='S',NWANT=n_eq,YWANT=y(:,i),YPWANT=dy,F=f, &
               !            WORK=workspace,WRKINT=workint,LENINT=6*n_eq)
               CALL INTRP (x_want, 'S', n_eq, y_out, dy, f, workspace, workint, 6*n_eq)
               x_want = x_want + dx
               y (:, i) = y_out
               i = i + 1
               IF (i > n_points) EXIT MainLoop
               IF (x_want > x_now) EXIT InterpResults
            END DO InterpResults
         END DO MainLoop
      END SUBROUTINE ODESolve
END MODULE ODE
!