! ************************************************************
! Runge-Kutta-Dormand-Prince integrator capable of integrating
! forwards only (t0<tn, pas>0)
! ************************************************************
SUBROUTINE RKDP54(EDO_DYDT, t0, y0, ndim_y, pas, tn, y, &
                  tol_abs, tol_rel, EDO_TRACE, n_traceunit)

IMPLICIT NONE

! Parametres l'adjustement du pas
DOUBLE PRECISION, PARAMETER :: max_fac = 2.0D0
DOUBLE PRECISION, PARAMETER :: min_fac = 0.5D0
DOUBLE PRECISION, PARAMETER :: sec_fac = 0.8D0


! Variables (muettes) de la liste d'arguments
DOUBLE PRECISION, INTENT(IN)        :: t0, tn
DOUBLE PRECISION, DIMENSION(ndim_y), INTENT(IN) :: y0
DOUBLE PRECISION, DIMENSION(ndim_y), INTENT(OUT) :: y
DOUBLE PRECISION, INTENT(INOUT)     :: pas
EXTERNAL                            :: EDO_DYDT
INTEGER, INTENT(IN)                 :: ndim_y
DOUBLE PRECISION, INTENT(IN)        :: tol_abs, tol_rel
EXTERNAL                            :: EDO_TRACE
INTEGER, INTENT(IN)                 :: n_traceunit


! Variables locales:
INTEGER :: j, nombre_pas_acc, nombre_pas_tot

DOUBLE PRECISION, DIMENSION(ndim_y) :: k1, k2, k3, k4, k5, k6, k7
DOUBLE PRECISION, DIMENSION(ndim_y) :: y1, y2, y3, y4, y5, y6, y7
DOUBLE PRECISION, DIMENSION(ndim_y) :: y_tilde, y_test
DOUBLE PRECISION                    :: delta, erreur_estimee
DOUBLE PRECISION                    :: t
DOUBLE PRECISION                    :: t1, t2, t3, t4, t5, t6, t7

DOUBLE PRECISION                    :: stiffness_diag
INTEGER                             :: n_stiffness_diag = 0

! Coefficients de la méthode de Dormand-Prince 5(4)
! Tableau a_ij
DOUBLE PRECISION, PARAMETER :: a21 =      1D0/     5D0
DOUBLE PRECISION, PARAMETER :: a31 =      3D0/    40D0, &
                               a32 =      9D0/    40D0
DOUBLE PRECISION, PARAMETER :: a41 =     44D0/    45D0, &
                               a42 =    -56D0/    15D0, &
                               a43 =     32D0/     9D0
DOUBLE PRECISION, PARAMETER :: a51 =  19372D0/  6561D0, &
                               a52 = -25360D0/  2187D0, &
                               a53 =  64448D0/  6561D0, &
                               a54 =   -212D0/   729D0
DOUBLE PRECISION, PARAMETER :: a61 =   9017D0/  3168D0, &
                               a62 =   -355D0/    33D0, &
                               a63 =  46732D0/  5247D0, &
                               a64 =     49D0/   176D0, &
                               a65 =  -5103D0/ 18656D0
DOUBLE PRECISION, PARAMETER :: a71 =     35D0/   384D0, &
                               a72 =      0D0         , &
                               a73 =    500D0/  1113D0, &
                               a74 =    125D0/   192D0, &
                               a75 =  -2187D0/  6784D0, &
                               a76 =     11D0/    84D0

! Tableaux b_j et b t_j
!  - Tableau b_j = a_7j
!  - Tableau b t_j
DOUBLE PRECISION, PARAMETER :: bt1 =   5179D0/ 57600D0, &
                               bt2 =      0D0         , &
                               bt3 =   7571D0/ 16695D0, &
                               bt4 =    393D0/   640D0, &
                               bt5 = -92097D0/339200D0, &
                               bt6 =    187D0/  2100D0, &
                               bt7 =      1D0/    40D0

! Tableaux c_i
DOUBLE PRECISION, PARAMETER ::  c1 =      0D0         , &
                                c2 =      1D0/     5D0, &
                                c3 =      3D0/    10D0, &
                                c4 =      4D0/     5D0, &
                                c5 =      8D0/     9D0, &
                                c6 =      1D0         , &
                                c7 =      1D0

! Initialiser les compteurs a 0
nombre_pas_acc = 0
nombre_pas_tot = 0

! Transcrire les conditions initiales dans les variables
! qui sont utilisees pour la recurrence
t = t0
y = y0

! EDO_TRACE n'est appele que si n_traceunit>0
IF(n_traceunit>0) THEN
  CALL EDO_TRACE(t, y, ndim_y, n_traceunit)
ENDIF


iterations: DO WHILE(t < tn)

  pas = MIN(pas, (tn-t)) ! limiter le pas pour ne pas dépasser tn

  t1 = t + pas*c1
  y1 = y
  CALL EDO_DYDT(t1, y1, ndim_y, k1)

  t2 = t + pas*c2
  y2 = y + pas*a21*k1
  CALL EDO_DYDT(t2, y2, ndim_y, k2)

  t3 = t + pas*c3
  y3 = y + pas*(a31*k1 + a32*k2)
  CALL EDO_DYDT(t3, y3, ndim_y, k3)

  t4 = t + pas*c4
  y4 = y + pas*(a41*k1 + a42*k2 + a43*k3)
  CALL EDO_DYDT(t4, y4, ndim_y, k4)

  t5 = t + pas*c5
  y5 = y + pas*(a51*k1 + a52*k2 + a53*k3 + a54*k4)
  CALL EDO_DYDT(t5, y5, ndim_y, k5)

  t6 = t + pas*c6
  y6 = y + pas*(a61*k1 + a62*k2 + a63*k3 + a64*k4 + a65*k5)
  CALL EDO_DYDT(t6, y6, ndim_y, k6)

  t7 = t + pas*c7
  y7 = y + pas*(a71*k1 + a72*k2 + a73*k3 + a74*k4 + a75*k5 + a76*k6)
  CALL EDO_DYDT(t7, y7, ndim_y, k7)

  y_test  = y7
  y_tilde = y + pas*(bt1*k1 + bt2*k2 + bt3*k3 + bt4*k4 + bt5*k5 + bt6*k6 + bt7*k7)

  stiffness_diag =   SUM((k7 - k6)*(k7 - k6)) &
                   / SUM((y7 - y6)*(y7 - y6))

  IF(stiffness_diag > 3.25D0) THEN
    n_stiffness_diag = n_stiffness_diag+1
    IF(n_stiffness_diag == 1) THEN
      WRITE(*,'("RKDP54: Stiffness first encountered at t = ", E13.6)') t
      WRITE(*,'("        stiffness diagnostic value: ", E13.6)') stiffness_diag
    ENDIF
!    IF(MOD(n_stiffness_diag,1000) == 0) THEN
!      WRITE(*,'("RKDP54: Until now, ",I0, " steps had stiffness diagnosed")') n_stiffness_diag
!    ENDIF
  ENDIF

  delta  = MAXVAL(ABS(y_test+y))/2D0*tol_rel + tol_abs
  erreur_estimee = MAXVAL(ABS(y_tilde - y_test))

  IF(erreur_estimee > delta*pas) THEN
    nombre_pas_tot=nombre_pas_tot+1

    IF(erreur_estimee >= ((delta*pas)*(sec_fac/min_fac)**4)) THEN
      pas = pas * min_fac
    ELSE
      pas = pas * sec_fac*((delta*pas)/erreur_estimee)**(1D0/4D0)
    ENDIF
  ELSE
    t = t + pas
    y = y_test
    nombre_pas_acc=nombre_pas_acc+1
    nombre_pas_tot=nombre_pas_tot+1

    IF(n_traceunit>0) THEN
      CALL EDO_TRACE(t, y, ndim_y, n_traceunit)
    ENDIF
    
    IF(erreur_estimee <= ((delta*pas)*(sec_fac/max_fac)**4)) THEN
      pas = pas * max_fac
    ELSE
      pas = pas * sec_fac*((delta*pas)/erreur_estimee)**(1D0/4D0)
    ENDIF
  ENDIF


ENDDO iterations


WRITE(*,'("RKDP54: number of steps accepted = ", I0)')   nombre_pas_acc
WRITE(*,'("RKDP54: number of steps rejected = ", I0)')   nombre_pas_tot-nombre_pas_acc
WRITE(*,'("RKDP54: number of steps calculated = ", I0)') nombre_pas_tot
WRITE(*,'("RKDP54: number of steps (total) with diagnosed stiffness ", I0)') n_stiffness_diag

RETURN
END SUBROUTINE RKDP54
