module rkf45FMM


! Procedura rkf45 
! "Computer Methods for Mathematical  Computations",
! by George E. Forsythe, Michael A. Malcolm, and Cleve B. Moler
! Prentice-Hall, 1977.
!
! Prepis pro Fortran90:
! Last updated: 28 February 2005 by Ralph Carmichael ralph@pdas.com 
! Public Domain Aeronautical Software 
! P.O. Box 1438 Santa Cruz CA 95061 USA 
!
! Drobne upravy pro g95
! Jan Cely, UFKL,Prirodovedecka fakulta MU, leden 2007


use std_type

implicit none
private

public :: rkf45

interface rkf45
	module procedure Rkf45Single,Rkf45Double
end interface



contains

SUBROUTINE Rkf45Single (F,y,t,tout,relerr,abserr,iflag,work,iwork)
! ---------------------------------------------------------------------------
! PURPOSE - Integrate a system of neqn first order differential 
!  equations by the Fehlberg fourth-fifth order Runge-Kutta method.
! NOTES - RKF45 is primarily designed to solve non-stiff and mildly stiff
!   differential equations when derivative evaluations are inexpensive.
!   Rkf45 should generally not be used when the user is demanding 
!	high accuracy.

! ABSTRACT - Subroutine Rkf45 integrates a system of neqn first order
!  ordinary differential equations of the form
!	dy(i)/dt = f(t,y(1),y(2),...,y(neqn))
!	where the y(i) are given at t .
! Typically the subroutine is used to integrate from t to tout but it can
! be used as a one-step integrator to advance the solution a single step
! in the direction of tout. On return, the parameters in the call list are
! set for continuing the integration. The user has only to call Rkf45 again
! (and perhaps define a new value for tout). Actually, Rkf45 is an 
! interfacing routine which calls subroutine Rkfs for the solution. 
! Rkfs in turn calls subroutine Fehl which computes an approximate solution
! over one step. 

! Rkf45 uses the Runge-Kutta-Fehlberg(4,5) method described in the reference
! E. Fehlberg:  Low-order Classical Runge-Kutta formulas with Stepsize 
! Control, NASA TR R-315.

! The performance of Rkf45 is illustrated in the reference
! L.F.Shampine,H.A.Watts,S.Davenport: "Solving Non-stiff Ordinary
! Differential Equations - The State of the Art 
!  Sandia Laboratories Report SAND75-0182, to appear in SIAM Review. 

! First call to Rkf45:

! The user must provide storage in his calling program for the arrays in 
! the call list	y(neqn), work(3+6*neqn), iwork(5), supply 
! subroutine F(t,y,yp) and initialize the following parameters-

! y -- vector of initial conditions
! t -- starting point of integration , must be a variable
! tout -- output point at which solution is desired.
!   t=tout is allowed on the first call only, in which case
!   Rkf45 returns with iflag=2 if continuation is possible.
! relerr,abserr -- relative and absolute local error tolerances
!   which must be non-negative. relerr must be a variable while
!   abserr may be a constant. The code should normally not be
!   used with relative error control smaller than about i.e-8 .
!   to avoid limiting precision difficulties the code requires
!   relerr to be larger than an internally computed relative
!   error parameter which is machine dependent. In particular,
!   pure absolute error is not permitted. If a smaller than
!   allowable value of relerr is attempted, rkf45 increases
!   relerr appropriately and returns control to the user before
!   continuing the integration.
! iflag -- +1,-1 indicator to initialize the code for each new problem. 
!   Normal input is +1. the user should set iflag=-1 only when one-step
!   integrator control is essential. In this case, Rkf45 attempts to 
!   advance the solution a single step in the direction of tout each time
!   it is called. Since this mode of operation results in extra computing 
!   overhead, it should be avoided unless needed. 

! Output from Rkf45:

! y -- solution at t
! t - last point reached in integration.
! iflag = 2 -- integration reached tout. Indicates successful return and is 
!              the normal mode for continuing integration.
!	= -2 -- a single successful step in the direction of tout has been 
!               taken. Normal mode for continuing integration one step 
!               at a time.
!	= 3 -- integration was not completed because relative error tolerance 
!              was too small. relerr has been increased	appropriately for 
!              continuing.
!	= 4 -- integration was not completed because more than 3000 derivative
!              evaluations were needed. This is approximately 500 steps.
!	= 5 -- integration was not completed because solution vanished making
!              a pure relative error test impossible. Must use non-zero abserr
!              to continue. Using the one-step integration mode for one step
!              is a good way to proceed.
!	= 6 -- integration was not completed because requested accuracy could 
!              not be achieved using smallest allowable stepsize. User must
!              increase the error tolerance before continued integration can 
!              be attempted.
!	= 7 -- it is likely that Rkf45 is inefficient for solving this problem.
!              Too much output is restricting the natural stepsize choice. 
!              Use the one-step integrator mode.
!	= 8 -- invalid input parameters. This indicator occurs if any of the 
!              following is satisfied -	SIZE(y) <= 0; 
!              t=tout and iflag /= +1 or -1
!	       relerr or abserr < 0.
!              iflag == 0 or < -2 or > 8
! work,iwork -- information which is usually of no interest to the user but 
!    necessary for subsequent calls.
!    work(1),...,work(neqn) contain the first derivatives of the solution 
!    vector y at t. work(neqn+1) contains the stepsize h to be attempted on 
!    the next step. iwork(1) contains the derivative evaluation counter. 

! Subsequent calls to Rkf45 

! Subroutine Rkf45 returns with all information needed to continue the 
! integration. If the integration reached tout, the user need only define
! a new tout and call Rkf45 again. In the one-step integrator mode (iflag=-2)
! the user must keep in mind that each step taken is in the direction of the
! current tout. Upon reaching tout (indicated by changing iflag to 2),the 
! user must then define a new tout and reset iflag to -2 to continue in the 
! one-step integrator mode. 

! If the integration was not completed but the user still wants to continue 
! (iflag=3,4 cases), he just calls Rkf45 again. With iflag=3, the relerr 
! parameter has been adjusted appropriately for continuing the integration.
! In the case of iflag=4 the function counter will be reset to 0 and 
! another 3000 function evaluations are allowed. 

! However,in the case iflag=5, the user must first alter the error criterion
! to use a positive value of abserr before integration can proceed. If he 
! does not, execution is terminated. 

! Also, in the case iflag=6, it is necessary for the user to reset iflag 
! to 2 (or -2 when the one-step integration mode is being used) as well as 
! increasing either abserr,relerr or both before the integration can be 
! continued. If this is not done, execution will be terminated. The 
! occurrence of iflag=6 indicates a trouble spot (solution is changing 
! rapidly,singularity may be present) and it often is inadvisable to continue.

! If (flag=7 is encountered, the user should use the one-step integration mode
! with the stepsize determined by the code or consider switching to the 
! Adams codes DE/STEP,INTRP. If the user insists upon continuing the
! integration with Rkf45, he must reset iflag to 2 before calling Rkf45 again.
! Otherwise, execution will be terminated.

! If iflag=8 is obtained, integration can not be continued unless the invalid
! input parameters are corrected. 

! It should be noted that the arrays work,iwork contain information required
! for subsequent integration. Accordingly, work and iwork should not be 
! altered.

! This interfacing routine merely relieves the user of a long calling list 
! via the splitting apart of two working storage arrays.

INTERFACE
  SUBROUTINE F(arg1,arg2,arg3)
  IMPLICIT NONE
  INTEGER,PARAMETER:: SP = SELECTED_REAL_KIND(6,20)
  REAL(KIND=SP),INTENT(IN):: arg1   ! independent variable
  REAL(KIND=SP),INTENT(IN),DIMENSION(:):: arg2   ! vector of dependent variables
  REAL(KIND=SP),INTENT(OUT),DIMENSION(:):: arg3   ! vector of derivatives
  END Subroutine F
END INTERFACE  


  REAL(KIND=SP), INTENT(IN OUT),DIMENSION(:):: y   ! solution vector at t
  REAL(KIND=SP), INTENT(IN OUT):: t   ! independent variable
  REAL(KIND=SP), INTENT(IN OUT):: tout   ! output point at which solution is desired
  REAL(KIND=SP), INTENT(IN OUT):: relerr   ! relative error tolerance
  REAL(KIND=SP), INTENT(IN):: abserr   !  absolute error tolerance
  INTEGER, INTENT(IN OUT):: iflag   !   indicator for status of work
  REAL(KIND=SP), INTENT(IN OUT),DIMENSION(:):: work   ! work array
  INTEGER, INTENT(IN OUT),DIMENSION(:):: iwork   ! work array

  INTEGER :: k1,k2,k3,k4,k5,k6,k1m
  INTEGER:: neqn
!----------------------------------------------------------------------------
  neqn=SIZE(y)
  k1m=neqn+1   !.....compute indices for the splitting of the work array
  k1=k1m+1
  k2=k1+neqn
  k3=k2+neqn
  k4=k3+neqn
  k5=k4+neqn
  k6=k5+neqn
  
  CALL RkfsSingle(F,neqn,y,t,tout,relerr,abserr,iflag,work(1:neqn),work(neqn+1),  &
    work(k1:k1+neqn-1),work(k2:k2+neqn-1),work(k3:k3+neqn-1),               &
    work(k4:k4+neqn-1), work(k5:k5+neqn-1),work(k6), work(k6+1),            &
    iwork(1),iwork(2),iwork(3),iwork(4),iwork(5))
  RETURN
END Subroutine Rkf45Single   ! ----------------------------------------------------

!+
SUBROUTINE RkfsSingle (F,neqn,y,t,tout,relerr,abserr,iflag,yp,                    &
    h,f1,f2,f3,f4,f5,savre,savae,nfe,kop,init,jflag,kflag)
! ---------------------------------------------------------------------------
! PURPOSE - integrate a system of neqn first order differential equations by 
!   the Fehlberg fourth-fifth order Runge-Kutta method.

! NOTES - Rkfs integrates a system of first order ordinary differential
!  equations as described in the comments for Rkf45. The arrays 
!  yp,f1,f2,f3,f4,and f5 (of dimension at least neqn) and the variables 
!  h,savre,savae,nfe,kop,init,jflag,and kflag are used internally by the 
!  code and appear in the call list to eliminate local retention of variables
!  between calls. Accordingly, they should not be altered. Items of possible 
!  interest are
!   yp - derivative of solution vector at t
!   h - an appropriate stepsize to be used for the next step
!   nfe- counter on the number of derivative function evaluations 

  INTEGER, INTENT(IN):: neqn
  REAL(KIND=SP), INTENT(IN OUT),DIMENSION(:):: y   ! solution vector at t
  REAL(KIND=SP), INTENT(IN OUT):: t   ! independent variable
  REAL(KIND=SP), INTENT(IN):: tout   ! output point at which solution is desired
  REAL(KIND=SP), INTENT(IN OUT):: relerr   ! relative error tolerance
  REAL(KIND=SP), INTENT(IN):: abserr   !  absolute error tolerance
  INTEGER, INTENT(IN OUT):: iflag  !   indicator for status of work
  REAL(KIND=SP), INTENT(IN OUT),DIMENSION(:):: yp   ! derivatives 
  REAL(KIND=SP), INTENT(OUT):: h                    ! step size
  REAL(KIND=SP), INTENT(IN OUT),DIMENSION(:):: f1
  REAL(KIND=SP), INTENT(IN OUT),DIMENSION(:):: f2
  REAL(KIND=SP), INTENT(IN OUT),DIMENSION(:):: f3
  REAL(KIND=SP), INTENT(IN OUT),DIMENSION(:):: f4
  REAL(KIND=SP), INTENT(IN OUT),DIMENSION(:):: f5
  REAL(KIND=SP), INTENT(OUT):: savre
  REAL(KIND=SP), INTENT(OUT):: savae
  INTEGER, INTENT(OUT):: nfe    ! number of function evaluations
  INTEGER, INTENT(OUT):: kop
  INTEGER, INTENT(IN OUT):: init
  INTEGER, INTENT(IN OUT):: jflag
  INTEGER, INTENT(IN OUT):: kflag
  LOGICAL :: hfaild,output

INTERFACE
  SUBROUTINE F(arg1,arg2,arg3)
  IMPLICIT NONE
  INTEGER,PARAMETER:: SP = SELECTED_REAL_KIND(6,20)
  REAL(KIND=SP),INTENT(IN):: arg1
  REAL(KIND=SP),INTENT(IN),DIMENSION(:):: arg2
  REAL(KIND=SP),INTENT(OUT),DIMENSION(:):: arg3
  END Subroutine F
END INTERFACE  

  REAL(KIND=SP):: a,ae,dt,ee,eeoet,esttol,et
  REAL(KIND=SP):: eps
  REAL(KIND=SP):: hmin
  INTEGER:: k,mflag

  REAL(KIND=SP),PARAMETER:: REMIN = 1.0E-12
     ! minimum acceptable value of relerr.  
    ! Attempts to obtain higher accuracy with this routine are usually very
    ! expensive and unsuccessful.

  INTEGER,PARAMETER:: MAXNFE = 3000 ! the expense is controlled by restricting
    ! the number of function evaluations to be approximately MAXNFE. 
    ! As set, this corresponds to about 500 steps.
  REAL(KIND=SP):: rer,s 
  REAL(KIND=SP):: scale                                   
  REAL(KIND=SP):: tol,toln,u26
  REAL(KIND=SP):: ypk
  REAL(KIND=SP),PARAMETER:: ZERO=0.0, ONE=1.0
!-----------------------------------------------------------------------
  eps=EPSILON(ZERO)
  u26=26.0*eps

!.....check input parameters............................................
  IF (neqn < 1 .OR. relerr < ZERO .OR. abserr < ZERO) THEN
    iflag=8
    RETURN
  END IF  
  mflag=ABS(iflag)
  IF (mflag == 0 .OR. mflag > 8) THEN
    iflag=8
    RETURN
  END IF
  
  IF (mflag == 1) GO TO 50

!.....check continuation posibilities...................................
  IF ((t == tout) .AND. (kflag /= 3)) THEN
    iflag=8   ! invalid flag handler
    RETURN
  END IF  
  IF (mflag /= 2) GO TO 25

!...you get here if iflag=+2 or -2
  IF ((kflag == 3) .OR. (init == 0)) GO TO 45
  IF (kflag == 4) GO TO 40
  IF ((kflag == 5) .AND. (abserr == ZERO)) GO TO 30
  IF ((kflag == 6) .AND. (relerr <= savre) .AND. (abserr <= savae))  GO TO 30
  GO TO 50


!...you get here if iflag=3,4,5,6,7, or 8.....
25 continue
   IF (iflag == 3) GO TO 45
   IF (iflag == 4) GO TO 40
   IF ((iflag == 5) .AND. (abserr > ZERO)) GO TO 45


!.....integration cannot be continued since user did not respond to the
!..... instructions pertaining to iflag=5,6,7, or 8
30  continue
	STOP   !!! maybe this is too strong


!.....reset function evaluation counter.................................
40 continue
   nfe=0
   IF (mflag == 2) GO TO 50


!.....reset flag value from previous call...............................
45 continue
   iflag=jflag
   IF (kflag == 3) mflag=ABS(iflag)

50 continue
   jflag=iflag   ! save for subsequent calls
   kflag=0       ! set the continuation flag
  savre=relerr   ! save for subsequent calls
  savae=abserr   ! save for subsequent calls


!.....restrict relative error tolerance to be at least as large as
!       2*eps+remin to avoid limiting precision difficulties arising
!       from impossible accuracy requests.
  rer=eps+eps+remin
  IF (relerr >= rer) GO TO 55


!.....relative error tolerance too small................................
  relerr=rer
  iflag=3
  kflag=3
  RETURN

55 continue
   dt=tout-t
   IF (mflag == 1) GO TO 60
   IF (init == 0) GO TO 65
   GO TO 80

!.....initialization....................................................
60 continue
   init=0
   kop=0
   a=t
   CALL f(a,y,yp)
   nfe=1
   IF (t /= tout) GO TO 65
   iflag=2
   RETURN

65 continue
   init=1
   h=ABS(dt)
   toln=ZERO
   DO  k=1,neqn
     tol=relerr*ABS(y(k))+abserr
     IF (tol <= ZERO) CYCLE
     toln=tol
     ypk=ABS(yp(k))
     IF (ypk*h**5 > tol) h=(tol/ypk)**0.2
   END DO

  IF (toln <= ZERO) h=ZERO
  h=MAX(h, u26*MAX(ABS(t), ABS(dt)))
  jflag=SIGN(2,iflag)

!.....set stepsize for integration in the direction from t to tout
80 continue
   h=SIGN(h,dt)

!.....test to see if rkf45 is being severely impacted by too many outputs
  IF (ABS(h) >= ABS(dt+dt)) kop=kop+1
  IF (kop /= 100) GO TO 85
  kop=0
  iflag=7
  RETURN

85 continue
   IF(ABS(dt) > u26*ABS(t)) GO TO 95

!.....if too close to output point, extrapolate and return
  DO  k=1,neqn
    y(k)=y(k)+dt*yp(k)
  END DO
  a=tout
  CALL F(a,y,yp)
  nfe=nfe+1
  GO TO 300


!.....initialize output point indicator.................................
95 continue
   output=.FALSE.


!.....scale the error tolerances to avoid premature underflow in the
!        error tolerance
  scale=2.0/relerr
  ae=scale*abserr


!.....step by step integration..........................................
N100: DO         !puvodne 100 continue 
   hfaild=.false.
!...set smallest allowable step size...
  hmin=u26*ABS(t)

!.....adjust stepsize if necessary to hit the output point.
!.....look ahead two steps to avoid drastic changes in the stepsize and
!.....thus lessen the impact of output points on the code.
  dt=tout-t
  IF (ABS(dt) >= ABS(h+h)) GO TO 200
  IF (ABS(dt) > ABS(h)) GO TO 150


!.....the next successful step will complete the integration to the
!     output point
  output=.TRUE.
  h=dt
  GO TO 200

150 continue 
   h=0.5*dt

!.....core integrator for taking a single step..........................

200 continue
N200: do  ! 200 continue 
   IF (nfe <= maxnfe) GO TO 220

!.....too much work...
  iflag=4
  kflag=4
  RETURN


!.....advance an approximate solution over one step of length h
220 continue
    CALL FehlSingle(f,y,t,h,yp,f1,f2,f3,f4,f5,f1)
    nfe=nfe+5


!.....compute and test allowable tolerances versus local error estimates
!.....  and remove scaling of tolerances. note that relative error is
!.....  is measured with respect to the average of the magnitudes of the
!.....   solution at the beginning and end of the step.
  eeoet=ZERO
  DO  k=1,neqn
    et=ABS(y(k))+ABS(f1(k))+ae
    IF (et > ZERO) GO TO 240
  
  
!.....Inappropriate error tolerance.....................................
    iflag=5
    RETURN
  
240 continue
    ee=ABS((-2090.0*yp(k)+(21970.0*f3(k)-15048.0*f4(k)))+   &
      (22528.0*f2(k)-27360.0*f5(k)))
    eeoet=MAX(eeoet,ee/et)
  END DO
  esttol=ABS(h)*eeoet*scale/752400.0
  IF (esttol <= ONE) GO TO 260


!.....unsuccessful step. reduce the stepsize and start again............
!     (decrease is limited to a factor of 1/10)
  hfaild=.TRUE.
  output=.FALSE.
  s=0.1
  IF (esttol < 59049.0) s=0.9/esttol**0.2
  h=s*h
  IF (ABS(h) > hmin) cycle N200
 end do N200 


!.....requested error unobtainable at smallest allowable stepsize.......
  iflag=6
  kflag=6
  RETURN


!.....successful step.   store solution at t+h and evaluate derivatives
!     there..........................
260 continue
    t=t+h
    DO  k=1,neqn
      y(k)=f1(k)
    END DO
    a=t
    CALL F(a,y,yp)
    nfe=nfe+1


! Choose next stepsize.  The increase is limited to a factor of 5.
! If step failure has just occured, next stepsize is not allowed
!  to increase.
  s=5.0
  IF (esttol > 1.889568E-4) s=0.9/esttol**0.2
  IF (hfaild) s=MIN(s,ONE)
  h=SIGN(MAX(s*ABS(h),hmin),h)

!.....E N D   O F   C O R E   I N T E G R A T O R   ....................


!.....Should we take another step? .....................................
  IF (output) GO TO 300
  IF (iflag >0) cycle N100
end do N100

!.....integration successfully completed................................
!...one step mode...
  iflag=-2
  RETURN

!...interval mode...
300 continue
    t=tout
    iflag=2
    RETURN
END SUBROUTINE RkfsSingle   ! -----------------------------------------------------

SUBROUTINE FehlSingle (F,y,t,h,yp,f1,f2,f3,f4,f5,s)
! ---------------------------------------------------------------------------
! PURPOSE - Integrate a system of neqn first-order ordinary differential
!  equations of the form dy/dt=F(t,y) [y is a vector], where the initial
!  values of y and the initial values of yp, the derivatives are specified
!  at the starting point t.  Fehl advances the solution over the fixed step h
!  and returns the fifth order (sixth-order locally) solution approximation 
!  at t+h in array s.
INTERFACE
  SUBROUTINE F(arg1,arg2,arg3)
  IMPLICIT NONE
  INTEGER,PARAMETER:: SP = SELECTED_REAL_KIND(6,20)
  REAL(KIND=SP),INTENT(IN):: arg1
  REAL(KIND=SP),INTENT(IN),DIMENSION(:):: arg2
  REAL(KIND=SP),INTENT(OUT),DIMENSION(:):: arg3
  END Subroutine F
END INTERFACE  

!!!  INTEGER, INTENT(IN):: neqn ! number of equations
  REAL(KIND=SP), INTENT(IN),DIMENSION(:):: y ! array of length neqn; function at t
  REAL(KIND=SP), INTENT(IN):: t   ! starting point
  REAL(KIND=SP), INTENT(IN):: h   ! step size
  REAL(KIND=SP), INTENT(IN),DIMENSION(:)::  yp ! array of length neqn; derivatives at t
  REAL(KIND=SP), INTENT(OUT),DIMENSION(:):: f1 ! array of length neqn for internal use
  REAL(KIND=SP), INTENT(OUT),DIMENSION(:):: f2 ! array of length neqn for internal use
  REAL(KIND=SP), INTENT(OUT),DIMENSION(:):: f3 ! array of length neqn for internal use
  REAL(KIND=SP), INTENT(OUT),DIMENSION(:):: f4 ! array of length neqn for internal use
  REAL(KIND=SP), INTENT(OUT),DIMENSION(:):: f5 ! array of length neqn for internal use
  REAL(KIND=SP), INTENT(OUT),DIMENSION(:):: s  ! array of length neqn; the results

  REAL(KIND=SP),PARAMETER:: C1 = 0.25_SP
  REAL(KIND=SP),PARAMETER:: C2 = 3.0_SP/32.0_SP
  REAL(KIND=SP),PARAMETER:: C3 = 3.0_SP
  REAL(KIND=SP),PARAMETER:: C4 = 3.0_SP/8.0_SP
  REAL(KIND=SP),PARAMETER:: C5 = 1.0_SP/2197.0_SP
  REAL(KIND=SP),PARAMETER:: C6 = 1932.0_SP
  REAL(KIND=SP),PARAMETER:: C7 = 7296.0_SP
  REAL(KIND=SP),PARAMETER:: C8 = -7200.0_SP
  REAL(KIND=SP),PARAMETER:: C9 = 12.0_SP/13.0_SP
  REAL(KIND=SP),PARAMETER:: C10 = 1.0_SP/4104.0_SP
  REAL(KIND=SP),PARAMETER:: C11 = 8341.0_SP
  REAL(KIND=SP),PARAMETER:: C12 = -845.0_SP
  REAL(KIND=SP),PARAMETER:: C13 = 29440.0_SP
  REAL(KIND=SP),PARAMETER:: C14 = -32832.0_SP
  REAL(KIND=SP),PARAMETER:: C15 = 1.0_SP/20520.0_SP
  REAL(KIND=SP),PARAMETER:: C16 = -6080.0_SP
  REAL(KIND=SP),PARAMETER:: C17 = 9295.0_SP
  REAL(KIND=SP),PARAMETER:: C18 = -5643.0_SP
  REAL(KIND=SP),PARAMETER:: C19 = 41040.0_SP
  REAL(KIND=SP),PARAMETER:: C20 = -28352.0_SP
  REAL(KIND=SP),PARAMETER:: C21 = 0.5_SP
  REAL(KIND=SP),PARAMETER:: C22 = 1.0_SP/7618050.0_SP
  REAL(KIND=SP),PARAMETER:: C23 = 902880.0_SP
  REAL(KIND=SP),PARAMETER:: C24 = 3855735.0_SP
  REAL(KIND=SP),PARAMETER:: C25 = -1371249.0_SP
  REAL(KIND=SP),PARAMETER:: C26 = 3953664.0_SP
  REAL(KIND=SP),PARAMETER:: C27 = 277020.0_SP

  REAL(KIND=SP):: ch
!-----------------------------------------------------------------------
  ch=C1*h  
  f5=y+ch*yp  
  CALL F(t+ch, f5,f1)

  ch=C2*h  
  f5=y + ch*(yp+C3*f1)
  CALL F(t+C4*h, f5,f2)  

  ch=C5*h
  f5=y + ch*(C6*yp + (C7*f2+C8*f1) )
  CALL F(t+C9*h, f5,f3)  

  ch=C10*h
  f5=y + ch*( (C11*yp + C12*f3) + (C13*f2+C14*f1) )
  CALL F(t+h, f5,f4)

  ch=C15*h
  f1=y + ch*( C16*yp + (C17*f3+C18*f4) + (C19*f1+C20*f2) )  
  CALL F(t+C21*h, f1, f5)

!.....compute approximate solution at t+h
  ch=C22*h
  s=y + ch*( C23*yp + (C24*f3+C25*f4) + (C26*f2+C27*f5) )  
  RETURN
END SUBROUTINE FehlSingle   ! ------------------------------------------------


! ========================= DP ============================




SUBROUTINE Rkf45Double (F,y,t,tout,relerr,abserr,iflag,work,iwork)
! ---------------------------------------------------------------------------
! PURPOSE - Integrate a system of neqn first order differential 
!  equations by the Fehlberg fourth-fifth order Runge-Kutta method.
! NOTES - RKF45 is primarily designed to solve non-stiff and mildly stiff
!   differential equations when derivative evaluations are inexpensive.
!   Rkf45 should generally not be used when the user is demanding 
!	high accuracy.

! ABSTRACT - Subroutine Rkf45 integrates a system of neqn first order
!  ordinary differential equations of the form
!	dy(i)/dt = f(t,y(1),y(2),...,y(neqn))
!	where the y(i) are given at t .
! Typically the subroutine is used to integrate from t to tout but it can
! be used as a one-step integrator to advance the solution a Double step
! in the direction of tout. On return, the parameters in the call list are
! set for continuing the integration. The user has only to call Rkf45 again
! (and perhaps define a new value for tout). Actually, Rkf45 is an 
! interfacing routine which calls subroutine Rkfs for the solution. 
! Rkfs in turn calls subroutine Fehl which computes an approximate solution
! over one step. 

! Rkf45 uses the Runge-Kutta-Fehlberg(4,5) method described in the reference
! E. Fehlberg:  Low-order Classical Runge-Kutta formulas with Stepsize 
! Control, NASA TR R-315.

! The performance of Rkf45 is illustrated in the reference
! L.F.Shampine,H.A.Watts,S.Davenport: "Solving Non-stiff Ordinary
! Differential Equations - The State of the Art 
!  Sandia Laboratories Report SAND75-0182, to appear in SIAM Review. 

! First call to Rkf45:

! The user must provide storage in his calling program for the arrays in 
! the call list	y(neqn), work(3+6*neqn), iwork(5), supply 
! subroutine F(t,y,yp) and initialize the following parameters-

! y -- vector of initial conditions
! t -- starting point of integration , must be a variable
! tout -- output point at which solution is desired.
!   t=tout is allowed on the first call only, in which case
!   Rkf45 returns with iflag=2 if continuation is possible.
! relerr,abserr -- relative and absolute local error tolerances
!   which must be non-negative. relerr must be a variable while
!   abserr may be a constant. The code should normally not be
!   used with relative error control smaller than about i.e-8 .
!   to avoid limiting precision difficulties the code requires
!   relerr to be larger than an internally computed relative
!   error parameter which is machine dependent. In particular,
!   pure absolute error is not permitted. If a smaller than
!   allowable value of relerr is attempted, rkf45 increases
!   relerr appropriately and returns control to the user before
!   continuing the integration.
! iflag -- +1,-1 indicator to initialize the code for each new problem. 
!   Normal input is +1. the user should set iflag=-1 only when one-step
!   integrator control is essential. In this case, Rkf45 attempts to 
!   advance the solution a Double step in the direction of tout each time
!   it is called. Since this mode of operation results in extra computing 
!   overhead, it should be avoided unless needed. 

! Output from Rkf45:

! y -- solution at t
! t - last point reached in integration.
! iflag = 2 -- integration reached tout. Indicates successful return and is 
!              the normal mode for continuing integration.
!	= -2 -- a Double successful step in the direction of tout has been 
!               taken. Normal mode for continuing integration one step 
!               at a time.
!	= 3 -- integration was not completed because relative error tolerance 
!              was too small. relerr has been increased	appropriately for 
!              continuing.
!	= 4 -- integration was not completed because more than 3000 derivative
!              evaluations were needed. This is approximately 500 steps.
!	= 5 -- integration was not completed because solution vanished making
!              a pure relative error test impossible. Must use non-zero abserr
!              to continue. Using the one-step integration mode for one step
!              is a good way to proceed.
!	= 6 -- integration was not completed because requested accuracy could 
!              not be achieved using smallest allowable stepsize. User must
!              increase the error tolerance before continued integration can 
!              be attempted.
!	= 7 -- it is likely that Rkf45 is inefficient for solving this problem.
!              Too much output is restricting the natural stepsize choice. 
!              Use the one-step integrator mode.
!	= 8 -- invalid input parameters. This indicator occurs if any of the 
!              following is satisfied -	SIZE(y) <= 0; 
!              t=tout and iflag /= +1 or -1
!	       relerr or abserr < 0.
!              iflag == 0 or < -2 or > 8
! work,iwork -- information which is usually of no interest to the user but 
!    necessary for subsequent calls.
!    work(1),...,work(neqn) contain the first derivatives of the solution 
!    vector y at t. work(neqn+1) contains the stepsize h to be attempted on 
!    the next step. iwork(1) contains the derivative evaluation counter. 

! Subsequent calls to Rkf45 

! Subroutine Rkf45 returns with all information needed to continue the 
! integration. If the integration reached tout, the user need only define
! a new tout and call Rkf45 again. In the one-step integrator mode (iflag=-2)
! the user must keep in mind that each step taken is in the direction of the
! current tout. Upon reaching tout (indicated by changing iflag to 2),the 
! user must then define a new tout and reset iflag to -2 to continue in the 
! one-step integrator mode. 

! If the integration was not completed but the user still wants to continue 
! (iflag=3,4 cases), he just calls Rkf45 again. With iflag=3, the relerr 
! parameter has been adjusted appropriately for continuing the integration.
! In the case of iflag=4 the function counter will be reset to 0 and 
! another 3000 function evaluations are allowed. 

! However,in the case iflag=5, the user must first alter the error criterion
! to use a positive value of abserr before integration can proceed. If he 
! does not, execution is terminated. 

! Also, in the case iflag=6, it is necessary for the user to reset iflag 
! to 2 (or -2 when the one-step integration mode is being used) as well as 
! increasing either abserr,relerr or both before the integration can be 
! continued. If this is not done, execution will be terminated. The 
! occurrence of iflag=6 indicates a trouble spot (solution is changing 
! rapidly,singularity may be present) and it often is inadvisable to continue.

! If (flag=7 is encountered, the user should use the one-step integration mode
! with the stepsize determined by the code or consider switching to the 
! Adams codes DE/STEP,INTRP. If the user insists upon continuing the
! integration with Rkf45, he must reset iflag to 2 before calling Rkf45 again.
! Otherwise, execution will be terminated.

! If iflag=8 is obtained, integration can not be continued unless the invalid
! input parameters are corrected. 

! It should be noted that the arrays work,iwork contain information required
! for subsequent integration. Accordingly, work and iwork should not be 
! altered.

! This interfacing routine merely relieves the user of a long calling list 
! via the splitting apart of two working storage arrays.

INTERFACE
  SUBROUTINE F(arg1,arg2,arg3)
  use std_type
  IMPLICIT NONE
  REAL(KIND=DP),INTENT(IN):: arg1   ! independent variable
  REAL(KIND=DP),INTENT(IN),DIMENSION(:):: arg2   ! vector of dependent variables
  REAL(KIND=DP),INTENT(OUT),DIMENSION(:):: arg3   ! vector of derivatives
  END Subroutine F
END INTERFACE  


  REAL(KIND=DP), INTENT(IN OUT),DIMENSION(:):: y   ! solution vector at t
  REAL(KIND=DP), INTENT(IN OUT):: t   ! independent variable
  REAL(KIND=DP), INTENT(IN OUT):: tout   ! output point at which solution is desired
  REAL(KIND=DP), INTENT(IN OUT):: relerr   ! relative error tolerance
  REAL(KIND=DP), INTENT(IN):: abserr   !  absolute error tolerance
  INTEGER, INTENT(IN OUT):: iflag   !   indicator for status of work
  REAL(KIND=DP), INTENT(IN OUT),DIMENSION(:):: work   ! work array
  INTEGER, INTENT(IN OUT),DIMENSION(:):: iwork   ! work array

  INTEGER :: k1,k2,k3,k4,k5,k6,k1m
  INTEGER:: neqn
!----------------------------------------------------------------------------
  neqn=SIZE(y)
  k1m=neqn+1   !.....compute indices for the splitting of the work array
  k1=k1m+1
  k2=k1+neqn
  k3=k2+neqn
  k4=k3+neqn
  k5=k4+neqn
  k6=k5+neqn
  
  CALL RkfsDouble(F,neqn,y,t,tout,relerr,abserr,iflag,work(1:neqn),work(neqn+1),  &
    work(k1:k1+neqn-1),work(k2:k2+neqn-1),work(k3:k3+neqn-1),               &
    work(k4:k4+neqn-1), work(k5:k5+neqn-1),work(k6), work(k6+1),            &
    iwork(1),iwork(2),iwork(3),iwork(4),iwork(5))
  RETURN
END Subroutine Rkf45Double
!-----------------------------------------------------------------------------

SUBROUTINE RkfsDouble (F,neqn,y,t,tout,relerr,abserr,iflag,yp,                    &
    h,f1,f2,f3,f4,f5,savre,savae,nfe,kop,init,jflag,kflag)
! ---------------------------------------------------------------------------
! PURPOSE - integrate a system of neqn first order differential equations by 
!   the Fehlberg fourth-fifth order Runge-Kutta method.

! NOTES - Rkfs integrates a system of first order ordinary differential
!  equations as described in the comments for Rkf45. The arrays 
!  yp,f1,f2,f3,f4,and f5 (of dimension at least neqn) and the variables 
!  h,savre,savae,nfe,kop,init,jflag,and kflag are used internally by the 
!  code and appear in the call list to eliminate local retention of variables
!  between calls. Accordingly, they should not be altered. Items of possible 
!  interest are
!   yp - derivative of solution vector at t
!   h - an appropriate stepsize to be used for the next step
!   nfe- counter on the number of derivative function evaluations 

  INTEGER, INTENT(IN):: neqn
  REAL(KIND=DP), INTENT(IN OUT),DIMENSION(:):: y   ! solution vector at t
  REAL(KIND=DP), INTENT(IN OUT):: t   ! independent variable
  REAL(KIND=DP), INTENT(IN):: tout   ! output point at which solution is desired
  REAL(KIND=DP), INTENT(IN OUT):: relerr   ! relative error tolerance
  REAL(KIND=DP), INTENT(IN):: abserr   !  absolute error tolerance
  INTEGER, INTENT(IN OUT):: iflag  !   indicator for status of work
  REAL(KIND=DP), INTENT(IN OUT),DIMENSION(:):: yp   ! derivatives 
  REAL(KIND=DP), INTENT(OUT):: h                    ! step size
  REAL(KIND=DP), INTENT(IN OUT),DIMENSION(:):: f1
  REAL(KIND=DP), INTENT(IN OUT),DIMENSION(:):: f2
  REAL(KIND=DP), INTENT(IN OUT),DIMENSION(:):: f3
  REAL(KIND=DP), INTENT(IN OUT),DIMENSION(:):: f4
  REAL(KIND=DP), INTENT(IN OUT),DIMENSION(:):: f5
  REAL(KIND=DP), INTENT(OUT):: savre
  REAL(KIND=DP), INTENT(OUT):: savae
  INTEGER, INTENT(OUT):: nfe    ! number of function evaluations
  INTEGER, INTENT(OUT):: kop
  INTEGER, INTENT(IN OUT):: init
  INTEGER, INTENT(IN OUT):: jflag
  INTEGER, INTENT(IN OUT):: kflag
  LOGICAL :: hfaild,output

INTERFACE
  SUBROUTINE F(arg1,arg2,arg3)
  use std_type
  IMPLICIT NONE
  REAL(KIND=DP),INTENT(IN):: arg1
  REAL(KIND=DP),INTENT(IN),DIMENSION(:):: arg2
  REAL(KIND=DP),INTENT(OUT),DIMENSION(:):: arg3
  END Subroutine F
END INTERFACE  

  REAL(KIND=DP):: a,ae,dt,ee,eeoet,esttol,et
  REAL(KIND=DP):: eps
  REAL(KIND=DP):: hmin
  INTEGER:: k,mflag

  REAL(KIND=DP),PARAMETER:: REMIN = 1.0E-12
    ! minimum acceptable value of relerr.  
    ! Attempts to obtain higher accuracy with this routine are usually very
    ! expensive and unsuccessful.

  INTEGER,PARAMETER:: MAXNFE = 3000 ! the expense is controlled by restricting
    ! the number of function evaluations to be approximately MAXNFE. 
    ! As set, this corresponds to about 500 steps.
  REAL(KIND=DP):: rer,s 
  REAL(KIND=DP):: scale                                   
  REAL(KIND=DP):: tol,toln,u26
  REAL(KIND=DP):: ypk
  REAL(KIND=DP),PARAMETER:: ZERO=0.0, ONE=1.0
!-----------------------------------------------------------------------
  eps=EPSILON(ZERO)
  u26=26.0*eps

!.....check input parameters............................................
  IF (neqn < 1 .OR. relerr < ZERO .OR. abserr < ZERO) THEN
    iflag=8
    RETURN
  END IF  
  mflag=ABS(iflag)
  IF (mflag == 0 .OR. mflag > 8) THEN
    iflag=8
    RETURN
  END IF
  
  IF (mflag == 1) GO TO 50

!.....check continuation posibilities...................................
  IF ((t == tout) .AND. (kflag /= 3)) THEN
    iflag=8   ! invalid flag handler
    RETURN
  END IF  
  IF (mflag /= 2) GO TO 25

!...you get here if iflag=+2 or -2
  IF ((kflag == 3) .OR. (init == 0)) GO TO 45
  IF (kflag == 4) GO TO 40
  IF ((kflag == 5) .AND. (abserr == ZERO)) GO TO 30
  IF ((kflag == 6) .AND. (relerr <= savre) .AND. (abserr <= savae))  GO TO 30
  GO TO 50


!...you get here if iflag=3,4,5,6,7, or 8.....
25 Continue
   IF (iflag == 3) GO TO 45
   IF (iflag == 4) GO TO 40
   IF ((iflag == 5) .AND. (abserr > ZERO)) GO TO 45


!.....integration cannot be continued since user did not respond to the
!..... instructions pertaining to iflag=5,6,7, or 8
30 Continue 
   STOP   !!! maybe this is too strong


!.....reset function evaluation counter.................................
40 Continue
   nfe=0
   IF (mflag == 2) GO TO 50


!.....reset flag value from previous call...............................
45 Continue
   iflag=jflag
   IF (kflag == 3) mflag=ABS(iflag)

50 Continue
   jflag=iflag   ! save for subsequent calls
   kflag=0       ! set the continuation flag
  savre=relerr   ! save for subsequent calls
  savae=abserr   ! save for subsequent calls


!.....restrict relative error tolerance to be at least as large as
!       2*eps+remin to avoid limiting precision difficulties arising
!       from impossible accuracy requests.
  rer=eps+eps+remin
  IF (relerr >= rer) GO TO 55


!.....relative error tolerance too small................................
  relerr=rer
  iflag=3
  kflag=3
  RETURN

55 Continue
   dt=tout-t
   IF (mflag == 1) GO TO 60
   IF (init == 0) GO TO 65
   GO TO 80

!.....initialization....................................................
60 Continue
   init=0
   kop=0
   a=t
   CALL f(a,y,yp)
   nfe=1
   IF (t /= tout) GO TO 65
   iflag=2
   RETURN

65 Continue
   init=1
   h=ABS(dt)
   toln=ZERO
   DO  k=1,neqn
     tol=relerr*ABS(y(k))+abserr
     IF (tol <= ZERO) CYCLE
     toln=tol
     ypk=ABS(yp(k))
     IF (ypk*h**5 > tol) h=(tol/ypk)**0.2
   END DO

  IF (toln <= ZERO) h=ZERO
  h=MAX(h, u26*MAX(ABS(t), ABS(dt)))
  jflag=SIGN(2,iflag)

!.....set stepsize for integration in the direction from t to tout
80 Continue
   h=SIGN(h,dt)

!.....test to see if rkf45 is being severely impacted by too many outputs
  IF (ABS(h) >= ABS(dt+dt)) kop=kop+1
  IF (kop /= 100) GO TO 85
  kop=0
  iflag=7
  RETURN

85 Continue
   IF(ABS(dt) > u26*ABS(t)) GO TO 95

!.....if too close to output point, extrapolate and return
  DO  k=1,neqn
    y(k)=y(k)+dt*yp(k)
  END DO
  a=tout
  CALL F(a,y,yp)
  nfe=nfe+1
  GO TO 300


!.....initialize output point indicator.................................
95 Continue
   output=.FALSE.


!.....scale the error tolerances to avoid premature underflow in the
!        error tolerance
  scale=2.0/relerr
  ae=scale*abserr


!.....step by step integration..........................................
100 Continue 
N100: DO
   hfaild=.false.
!...set smallest allowable step size...
  hmin=u26*ABS(t)

!.....adjust stepsize if necessary to hit the output point.
!.....look ahead two steps to avoid drastic changes in the stepsize and
!.....thus lessen the impact of output points on the code.
  dt=tout-t
  IF (ABS(dt) >= ABS(h+h)) GO TO 200
  IF (ABS(dt) > ABS(h)) GO TO 150


!.....the next successful step will complete the integration to the
!     output point
  output=.TRUE.
  h=dt
  GO TO 200

150 Continue 
   h=0.5*dt

!.....core integrator for taking a Double step..........................

200 continue
N200:  do
    IF (nfe <= maxnfe) GO TO 220

!.....too much work...
  iflag=4
  kflag=4
  RETURN


!.....advance an approximate solution over one step of length h
220 Continue
    CALL FehlDouble(f,y,t,h,yp,f1,f2,f3,f4,f5,f1)
    nfe=nfe+5


!.....compute and test allowable tolerances versus local error estimates
!.....  and remove scaling of tolerances. note that relative error is
!.....  is measured with respect to the average of the magnitudes of the
!.....   solution at the beginning and end of the step.
  eeoet=ZERO
  DO  k=1,neqn
    et=ABS(y(k))+ABS(f1(k))+ae
    IF (et > ZERO) GO TO 240
  
  
!.....Inappropriate error tolerance.....................................
    iflag=5
    RETURN
  
240 Continue
    ee=ABS((-2090.0*yp(k)+(21970.0*f3(k)-15048.0*f4(k)))+               &
      (22528.0*f2(k)-27360.0*f5(k)))
    eeoet=MAX(eeoet,ee/et)
  END DO
  esttol=ABS(h)*eeoet*scale/752400.0
  IF (esttol <= ONE) GO TO 260


!.....unsuccessful step. reduce the stepsize and start again............
!     (decrease is limited to a factor of 1/10)
  hfaild=.TRUE.
  output=.FALSE.
  s=0.1
  IF (esttol < 59049.0) s=0.9/esttol**0.2
  h=s*h
  IF (ABS(h) <= hmin) exit N200
end do N200

!.....requested error unobtainable at smallest allowable stepsize.......
  iflag=6
  kflag=6
  RETURN


!.....successful step.   store solution at t+h and evaluate derivatives
!     there..........................
260 Continue
    t=t+h
    DO  k=1,neqn
      y(k)=f1(k)
    END DO
    a=t
    CALL F(a,y,yp)
    nfe=nfe+1


! Choose next stepsize.  The increase is limited to a factor of 5.
! If step failure has just occured, next stepsize is not allowed
!  to increase.
  s=5.0
  IF (esttol > 1.889568E-4) s=0.9/esttol**0.2
  IF (hfaild) s=MIN(s,ONE)
  h=SIGN(MAX(s*ABS(h),hmin),h)

!.....E N D   O F   C O R E   I N T E G R A T O R   ....................


!.....Should we take another step? .....................................
  IF (output) GO TO 300
  IF (iflag <= 0) EXIT N100
end do N100


!.....integration successfully completed................................
!...one step mode...
  iflag=-2
  RETURN

!...interval mode...
300 Continue
    t=tout
    iflag=2
    RETURN
END SUBROUTINE RkfsDouble   ! -----------------------------------------------------

!+
SUBROUTINE FehlDouble (F,y,t,h,yp,f1,f2,f3,f4,f5,s)
! ---------------------------------------------------------------------------
! PURPOSE - Integrate a system of neqn first-order ordinary differential
!  equations of the form dy/dt=F(t,y) [y is a vector], where the initial
!  values of y and the initial values of yp, the derivatives are specified
!  at the starting point t.  Fehl advances the solution over the fixed step h
!  and returns the fifth order (sixth-order locally) solution approximation 
!  at t+h in array s.
INTERFACE
  SUBROUTINE F(arg1,arg2,arg3)
  IMPLICIT NONE
  INTEGER,PARAMETER:: DP = SELECTED_REAL_KIND(10,50)
  REAL(KIND=DP),INTENT(IN):: arg1
  REAL(KIND=DP),INTENT(IN),DIMENSION(:):: arg2
  REAL(KIND=DP),INTENT(OUT),DIMENSION(:):: arg3
  END Subroutine F
END INTERFACE  

!!!  INTEGER, INTENT(IN):: neqn ! number of equations
  REAL(KIND=DP), INTENT(IN),DIMENSION(:):: y ! array of length neqn; function at t
  REAL(KIND=DP), INTENT(IN):: t   ! starting point
  REAL(KIND=DP), INTENT(IN):: h   ! step size
  REAL(KIND=DP), INTENT(IN),DIMENSION(:)::  yp ! array of length neqn; derivatives at t
  REAL(KIND=DP), INTENT(OUT),DIMENSION(:):: f1 ! array of length neqn for internal use
  REAL(KIND=DP), INTENT(OUT),DIMENSION(:):: f2 ! array of length neqn for internal use
  REAL(KIND=DP), INTENT(OUT),DIMENSION(:):: f3 ! array of length neqn for internal use
  REAL(KIND=DP), INTENT(OUT),DIMENSION(:):: f4 ! array of length neqn for internal use
  REAL(KIND=DP), INTENT(OUT),DIMENSION(:):: f5 ! array of length neqn for internal use
  REAL(KIND=DP), INTENT(OUT),DIMENSION(:):: s  ! array of length neqn; the results

  REAL(KIND=DP),PARAMETER:: C1 = 0.25_DP
  REAL(KIND=DP),PARAMETER:: C2 = 3.0_DP/32.0_DP
  REAL(KIND=DP),PARAMETER:: C3 = 3.0_DP
  REAL(KIND=DP),PARAMETER:: C4 = 3.0_DP/8.0_DP
  REAL(KIND=DP),PARAMETER:: C5 = 1.0_DP/2197.0_DP
  REAL(KIND=DP),PARAMETER:: C6 = 1932.0_DP
  REAL(KIND=DP),PARAMETER:: C7 = 7296.0_DP
  REAL(KIND=DP),PARAMETER:: C8 = -7200.0_DP
  REAL(KIND=DP),PARAMETER:: C9 = 12.0_DP/13.0_DP
  REAL(KIND=DP),PARAMETER:: C10 = 1.0_DP/4104.0_DP
  REAL(KIND=DP),PARAMETER:: C11 = 8341.0_DP
  REAL(KIND=DP),PARAMETER:: C12 = -845.0_DP
  REAL(KIND=DP),PARAMETER:: C13 = 29440.0_DP
  REAL(KIND=DP),PARAMETER:: C14 = -32832.0_DP
  REAL(KIND=DP),PARAMETER:: C15 = 1.0_DP/20520.0_DP
  REAL(KIND=DP),PARAMETER:: C16 = -6080.0_DP
  REAL(KIND=DP),PARAMETER:: C17 = 9295.0_DP
  REAL(KIND=DP),PARAMETER:: C18 = -5643.0_DP
  REAL(KIND=DP),PARAMETER:: C19 = 41040.0_DP
  REAL(KIND=DP),PARAMETER:: C20 = -28352.0_DP
  REAL(KIND=DP),PARAMETER:: C21 = 0.5_DP
  REAL(KIND=DP),PARAMETER:: C22 = 1.0_DP/7618050.0_DP
  REAL(KIND=DP),PARAMETER:: C23 = 902880.0_DP
  REAL(KIND=DP),PARAMETER:: C24 = 3855735.0_DP
  REAL(KIND=DP),PARAMETER:: C25 = -1371249.0_DP
  REAL(KIND=DP),PARAMETER:: C26 = 3953664.0_DP
  REAL(KIND=DP),PARAMETER:: C27 = 277020.0_DP

  REAL(KIND=DP):: ch
!-----------------------------------------------------------------------
  ch=C1*h  
  f5=y+ch*yp  
  CALL F(t+ch, f5,f1)

  ch=C2*h  
  f5=y + ch*(yp+C3*f1)
  CALL F(t+C4*h, f5,f2)  

  ch=C5*h
  f5=y + ch*(C6*yp + (C7*f2+C8*f1) )
  CALL F(t+C9*h, f5,f3)  

  ch=C10*h
  f5=y + ch*( (C11*yp + C12*f3) + (C13*f2+C14*f1) )
  CALL F(t+h, f5,f4)

  ch=C15*h
  f1=y + ch*( C16*yp + (C17*f3+C18*f4) + (C19*f1+C20*f2) )  
  CALL F(t+C21*h, f1, f5)

!.....compute approximate solution at t+h
  ch=C22*h
  s=y + ch*( C23*yp + (C24*f3+C25*f4) + (C26*f2+C27*f5) )  
  RETURN
END SUBROUTINE FehlDouble
!-----------------------------------------------------------------------------


end module rkf45FMM