module IVPEPP ! Version 20121207
! This file IVPEPP.F90 is from the EPPEER package for the solution
! of ordinary initial value problems by explicit parallel peer two-step
! methods with automatic stepsize control and continuous output written by
! B.A. Schmitt and R. Weiner available at www.mathematik.uni-marburg.de/~schmitt/peer.
! Parallelization with 2..8 computing cores is provided by the peer methods
! and OpenMP, no parallelization of the right-hand side fcn of the ODE is required.
! In addition to this module also the file IVPRKP.F90 from the package is required.
! This is only a short documentation, a PDF manual man_epp.pdf comes with the package,
! please look there for detailed information. A minimal calling sequence is
!
! call PPSETCOEFF(mnr,stage,methname)
!!-- set start&end times t,tend and initial value y
! call EPPEER(fcn,t,tend,y,cpar)
! irep = PPREPORT(.true.)
! call PPFREECOEFF
!
! where the parameters have the following meaning:
!  MNR: in, integer in the range 2..9, chooses order and degree of parallelization,
!    the number of stages and the order of the method will be MNR
!  STAGE: out, integer returns the chosen number of stages of the method
!  METHNAME: out, string with the name of the method
!  FCN: pure subroutine coding the right-hand side $f(t,u) of the ODE supplied by
!    user, it has to be defined in the form
!     pure subroutine fcn(t,u,udot,par)
!     real(8), intent(in) :: t                   ! time of evaluation
!     real(8), dimension(:), intent(in) :: u     ! actual state of the system
!     real(8), dimension(:), intent(out) :: udot ! result f(t,u)
!     real(8), dimension(:), intent(in) :: par   ! additional parameters of the ODE
!     end subroutine fcn
!  T: in+out, real; start time of integration, on exit: current time of integration
!  TEND: in, real; end time of initital value problem
!  Y: in+out, array; initial value y_0 for IVP, on exit: solution at time TEND,
!     size(Y)=n defines space dimension of ODE
!  CPAR: in, array; additional parameter of the ODE, is passed to calls of FCN through PAR
!  IREP: integer, return code, integration successful for irep==0, error messages
!        are printed to console by call PPREPORT(.true.)
!
! The minimal calling sequence performs integration with automatic stepsize control
! with tolerances 1D-5 and no continuous output. Changes in operation mode are obtained by:
! * call PPSETJOB(stepcon,contout)
!    must be called after PPSETCOEFF, both parameters are input&optional of type logical,
!     STEPCON=.false. switches to fixed stepsize integration --> see PPSETACC
!     CONTOUT=.true.  switches continuous output on --> see PPCONT/SOLOUT
! * call PPSETACC(atol,rtol,hstep)
!    must be called after PPSETCOEFF, changes accuracies or initial stepsize
!     ATOL: in, real; changes absolute error tolerance
!     RTOL: in, real; changes relative error tolerance
!     HSTEP: in, real, optional: sets initial stepsize or constant stepsize for case
!        stepcon=.false.; HSTEP.le.0 triggers stepsize control
! Additional information during/after operation of EPPEER is obtained by
! * call PPGETSTAS(nsteps,nrej,nfcn)
!    to be called after EPPEER gets integration statistics, all parameters are integer&output:
!     NSTEPS: number of all time steps, including rejected steps
!     NREJ:  number of rejected time steps
!     NFCN:  total number of calls to FCN, including startup
! * call PPCONT(ts,tnew,ycon)
!    This module ivpepp saves the stages from the current step and PPCONT evaluates an
!    interpolation polynomial having the approximation order MNR=STAGE. Parameters are:
!     TS: in, real; time, where solution y(ts) is to be computed, should lie in current
!        time interval, TS.le.TNEW
!     TNEW, in, real: end point of current time interval (e.g. TEND after exit from EPPEER)
!     YCON, out, array: value of solution at TS
!-----
! Dense output is computed by setting CONTOUT=.true. and calling EPPEER with an extended
! parameter list
!   call EPPEER(fcn,t,tend,y,cpar,solout)
! Here, solout is a user-supplied subroutine which is called by EPPEER after every
! successful step. The declaration must have the form
!   subroutine SOLOUT(n,ts,tnew)
!     N: in, integer; problem dimension used to allocate local variable(s) YCON(1:N)
!     TS: in/out, real; point of solution evaluation (see below)
!     TNEW: in, real; end of current time interval, passed to PPCONT
! Usage: EPPEER initializes a variable TSOL=T0 and calls SOLOUT after every succesfull step
! in the form CALL SOLOUT(n,TSOL,T+HS), where T+HS is the end of the current time interval.
! SOLOUT may then call PPCONT(ts,tnew,ycon) repeatedly, save or print YCON, and increment TS
! according to the users needs. SOLOUT exits if TS>TNEW (example in driver IVP_PMAIN.F90).
! Note that only for the FSAL methods (mnr=3,5,7,9) dense output is continuous.
!----
! Parallel execution is obtained with OpenMP, recommended command-lines using the GNU gfortran
! compiler (e.g.) for the example ODE in the file MBOD4H.F90 and main program IVP_PMAIN.F90 are:
!    gfortran -c ivprkp.f90 mbod4h.f90
!    gfortran -c -fopenmp ivpepp.f90
!    gfortran -fopenmp ivprkp.o ivpepp.o mbod4h.o ivp_pmain.f90
! The execubale A (A.exe) can now be started. Note, that only the peer integrator IVPEPP and the
! main program IVP_PMAIN.F90 execute parallel.
!
!---------- Begin of IVPEPP module code ------------
  use ivprkp ! Runge-Kutta for startup
!--  all variables private
  integer, save, private :: stages ! stage number of peer method
  integer, save, private :: order, ifcn, isteps, irej, maxsteps, irep
  real(8), save, private :: atoli, rtoli, phs, hmin, hmax, sigx
  real(8), save, private :: ptim0  ! start time of current interval, used by ppcont
  logical, save, private :: pfsal, pcout, stpci
  ! matrices
  integer, dimension(20), private :: idx  ! memory pointer in PYY, used idx(1:2*stages)
  real(8), dimension(:), allocatable, save, private :: pc, pe
  real(8), dimension(:,:), allocatable, save, private :: pb, pa0, pcv, ppv, pyy

  contains

!-- EPPEER performs explicit parallel peer integration from t to tend
  subroutine eppeer(fcn,t,tend,y,cpar,solout)
    implicit none
    interface
      pure subroutine fcn(t,u,udot,par)
       real(8), intent(in) :: t
       real(8), dimension(:), intent(in) :: u
       real(8), dimension(:), intent(out) :: udot
       real(8), dimension(:), intent(in) :: par
      end subroutine fcn
      subroutine solout(n,ts,tnew) ! uses private idx, pyy through pcont
       integer, intent(in) :: n
       real(8) :: ts
       real(8), intent(in) :: tnew
      end subroutine solout
    end interface
	optional solout ! if solout is omitted => pcout=.false.
    real(8) :: t,tend
    real(8), dimension(:) :: y
    real(8), dimension(:), intent(in) :: cpar
    integer :: n,stg,ist1,ic,id,new ! local variables
    integer, dimension(2*stages) :: swap ! stage addresses
    real(8) tc,t0,tsol, rer,rc, h0,hn, sig,sighs, ats,rts
    real(8), dimension(stages,stages) :: pa ! coefficient matrix A
    real(8), dimension(size(y)) :: yf
    real(8), dimension(:,:), allocatable :: ff ! stage derivatives
    logical :: sscont ! no stepsize control for stepcon=.false., phs>0
    if (allocated(pyy)) deallocate(pyy)
	if (.not.present(solout)) pcout=.false. ! solout must not be called
    n = size(y) ! problem dimension
	!   -- prepare fsal, indirect addressing
    ist1 = 1; if (pfsal) ist1 = 2
    new = stages+1-ist1 ! new+1 = first new stage
    do stg = 1,2*stages ! old stages first, then new ones with overlap for fsal
      idx(stg) = stg
    end do
    allocate(pyy(n,new+stages)) ! memory for stages
    sscont = (stpci.or.(phs.le.0.D0)) ! switch for stepsize control
    isteps = 0;  irej = 0; irep = 0
    hmax = tend-t ! used by inisteps
	if (phs>hmax) phs = hmax
    if (stages<2.or..not. allocated(pb)) irep=-1 ! ERROR: ppsetcoeff not called ?!
    id = 5  ! startup with DOPRI
    call rksetcoeff(id) ! choose RK method
    rc=1D-1; if (.not.sscont) rc=1D9
    rts = max(rc*rtoli,1D-14) ! higher accuracy for starting proc
    ats = max(rc*atoli,1D-14)
    call rksetparam(rts,ats) ! set RK-parameters
    call fcn(t,y,yf,cpar); ifcn=1
    if (phs<=0.D0) & ! guess initial stepsize
      call inisteps(n,fcn,t,y,yf,phs,cpar)
    irep = 1 ! only one step rkmeth (adjusts stepsize phs for last stage)
    pyy(:,idx(stages)) = y; tc = t; hn = phs
    call rkmeth(n,fcn,tc,t+phs,pyy(:,idx(stages)),hn,irep,cpar); ifcn=ifcn+rkfcn
    call rksetparam(1D9*rtoli,1D9*atoli) ! no stepsize control again
    ! stepsize adjusted, tc = end of start interval:
    h0 = tc-t; t0 = t
    do stg = 1,stages-1 ! startup procedure
      pyy(:,idx(stg)) = y
      if (abs(pc(stg)).ge.1D-12) then
        tc = t0; hn = h0*pc(stg); id = 1
        call rkmeth(n,fcn,tc,t0+h0*pc(stg),pyy(:,idx(stg)),hn,id,cpar)
        ifcn=ifcn+rkfcn
      end if
    end do
    call rkfreecoeff
	allocate(ff(n,new+stages)) ! RK finished, memory for derivatives
    if (ist1>1) ff(:,idx(1)) = yf
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(stg) SCHEDULE(STATIC)
	do stg = ist1,stages
  	  call fcn(t0+h0*pc(stg),pyy(:,idx(stg)),ff(:,idx(stg)),cpar)
	end do
!$OMP END DO
!$OMP END PARALLEL
	ifcn=ifcn+stages+1-ist1
! end of start interval
    t = t0+h0
	ptim0 = t0; tsol=t0 ! for solout
    if (pcout) call solout(n,tsol,t) ! dense output, may use idx, pyy through ppcont
    y = pyy(:,idx(stages))  ! solution at T
    phs = h0 ! continue with stepsize h0
    if (irep==-2) irep = 0 ! was OK in rkmeth
    do while ((t<tend-atoli).and.(irep==0)) ! time steps, repeat until tend
      isteps = isteps+1
      sig = phs/h0
      do stg = ist1,stages! compute coeff. h0*A row-wise
        sighs = h0
        pa(stg,:) = h0*pa0(stg,:)
        do ic = 1,stages
          sighs = sighs*sig
          pa(stg,:) = pa(stg,:)+pcv(stg,ic)*sighs*ppv(ic,:)
        end do
      end do
      ! compute new stage solutions, summation of small entries first
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(stg,ic) SCHEDULE(STATIC)
      do stg = ist1,stages
        pyy(:,idx(new+stg)) = 0.D0
        do ic = 1,stages
          pyy(:,idx(new+stg)) = pyy(:,idx(new+stg))+pa(stg,ic)*ff(:,idx(ic))
        end do
      end do
!$OMP END DO
!$OMP DO PRIVATE(stg,ic) SCHEDULE(STATIC)
      do stg = ist1,stages
        do ic = 1,stages
          pyy(:,idx(new+stg)) = pyy(:,idx(new+stg))+pb(stg,ic)*pyy(:,idx(ic))
        end do
      end do
!$OMP END DO
     ! function evaluations Fn=fcn(Yn) parallel:
!$OMP DO PRIVATE(stg) SCHEDULE(STATIC)
      do stg = ist1,stages
        call fcn(t+phs*pc(stg),pyy(:,idx(new+stg)),ff(:,idx(new+stg)),cpar)
      end do
!$OMP END DO
!$OMP SINGLE
      ifcn = ifcn+stages+1-ist1
      rer = 0D0
!$OMP END SINGLE
      ! compute error estimate with F(new):
!$OMP DO PRIVATE(stg,rc) REDUCTION(MAX:rer)
      do id =1,n
        rc = pe(1)*ff(id,idx(new+1))
        do stg = 2,stages
          rc = rc+pe(stg)*ff(id,idx(new+stg))
        end do
        rer = max(rer,abs(rc)/(atoli+rtoli*abs(y(id))))
      end do
!$OMP END DO
!$OMP END PARALLEL
      rer = phs*rer ! quotient of error estimate & tolerance
      ! estimate step ratio:
      hn =exp(-log(1.D-8+rer)/order)
      hn = phs*max(0.2D0,min(0.95D0*hn,sigx)) ! safety interval
      if ((rer > 1.D0).and.sscont) then ! reject step
        irej = irej+1
        if (hn<hmin) irep=1! report failure
      else ! accept step, swap old & new
        swap(1:stages) = idx(new+1:new+stages)
        swap(stages+1:stages+new) = idx(1:new)
        idx = swap
        y = pyy(:,idx(stages))
        ptim0 = t ! save for call of ppcont
        if (pcout) call solout(n,tsol,t+phs) ! dense output, (idx, pyy through ppcont)
        t = t+phs
        h0 = phs
      end if ! of accepted step
      if (isteps > maxsteps) irep=2 ! failure, too many steps
      if (sscont) phs = hn
     if (t+phs>tend) phs = tend-t
    end do ! integration steps
	deallocate(ff)
  end subroutine eppeer

!-- PPSETJOB sets switches (step control etc), optional call after ppsetcoeff
  subroutine ppsetjob(stepcon,contout)
    implicit none
    logical, intent(in), optional :: stepcon, contout
    if (present(stepcon)) stpci = stepcon
	if (present(contout)) pcout = contout
  end subroutine ppsetjob

  !-- PPSETPARAM sets integration parameters, optional call after ppsetcoeff
  subroutine ppsetacc(atol,rtol,hstep)
    implicit none
    real(8), intent(in) :: atol, rtol
    real(8), intent(in), optional :: hstep
    if (atol>0D0) atoli = atol
    if (rtol>0D0) rtoli = rtol
    if (present(hstep)) phs = hstep
  end subroutine ppsetacc

!-- PPREPORT returns status report
  ! if (ppreport==0) : PPMETH successful
  function ppreport(prterr)
    implicit none
	logical, intent(in), optional :: prterr
    integer :: ppreport
    ppreport = irep
	if (present(prterr).and.prterr) then
	  select case(irep)
	  case (-1)
	    print *, "EPPEER: missing inititalization ppsetparam"
	  case (1)
        print *, "EPPEER failure: step size too small (<hmin)"
	  case (2)
	    print *, "EPPEER failure: too many steps (>maxsteps)"
	  end select
	end if
  end function ppreport

!-- PPCONT computes solution y(ts)->ycon, uses IDX and PYY, may be called by SOLOUT
  subroutine ppcont(ts,tnew,ycon)
    implicit none
    real(8), intent(in) :: ts,tnew
    real(8), dimension(:), intent(out) :: ycon
    integer :: j,k,n
    real(8) :: xi, sv
    real(8), dimension(stages):: cd, lag
    xi =(ts-ptim0)/(tnew-ptim0)
    do k=1,stages
      cd(k) = xi-pc(k)
    end do
    ! compute Lagrangian polynomials L_k(xi)
!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(k,j) SCHEDULE(STATIC)
    do k=1,stages
      lag(k) = 1.D0
      do j=1,stages
       if (j/=k) lag(k) = lag(k)*cd(j)/(cd(j)-cd(k))
      end do
    end do
!$OMP END DO
    n = size(ycon)
!$OMP DO PRIVATE(j,k,sv)
    do j = 1,n
	  sv = 0.D0
      do k=1,stages
        sv = sv+lag(k)*pyy(j,idx(k))
      end do
	  ycon(j) = sv
	end do
!$OMP END DO
!$OMP END PARALLEL
  end subroutine ppcont

  !-- VANDM solves Vandermonde systems PM_new*V=PM
  subroutine vandm(n,pc,pm) ! [Björck&Pereyra/1970]
    implicit none
    integer :: n, j,k
    real(8), dimension(n), intent(in) :: pc
    real(8), dimension(:,:) :: pm
    do k=1,n-1
      do j=n,k+1,-1
        pm(:,j) = pm(:,j)-pm(:,j-1)*pc(k)
      end do
    end do
    do k=n-1,1,-1
      do j=k+1,n
        pm(:,j) = pm(:,j)/(pc(j)-pc(j-k))
      end do
      do j=k,n-1
        pm(:,j) = pm(:,j)-pm(:,j+1)
      end do
    end do
  end subroutine vandm

  !-- PPSETCOEFF creates coefficient matrices, sets defaults
  subroutine ppsetcoeff(mnr,pstage,methname)
    implicit none
    integer :: mnr, i,j
    integer, intent(out) :: pstage
    character(len=16), intent(out) :: methname
    real(8) :: s
    stages = mnr-(mnr/10)*10 ! last digit of MNR
	pstage = stages
    allocate(pb(stages,stages),pa0(stages,stages),pcv(stages,stages))
    allocate(ppv(stages,stages),pc(stages),pe(stages))
    pfsal = .false.
    select case(mnr)
    case(3) ! EPP3f4
      methname = "EPP3f4"
      order = 4 ! stages = 3
      sigx = 1.8D0
	  pfsal = .true.
      pc = (/ 0.D0, 0.75D0, 1.D0 /)
      pb(1,:)=(/0.D0, 0.D0, 1.D0 /)
      pb(2,:)=(/-12163/17521920D0,12163/72D2,-60328969/876096D2 /)
	  pb(3,:)=(/-1/144D1,169/1D2,-4963/72D2/)
    case (4) ! EPP4y3
      methname = "EPP4y3"
      order = 5 ! stages = 4
      sigx = 1.6D0 ! 20060411
      pc = (/ 1.33880820864483004D0,1.70380840062134099D0,1.86823097147835395D0, 1.D0 /)
      pb(1,:)=(/ 0.00534248151938605D0,-0.12787556805087963D0,-0.01346305501819259D0,&
       1.13599614154968617D0/)
      pb(2,:)=(/ 0.00291846709251725D0,-0.00273308699096256D0,-0.07182299833148626D0,&
       1.07163761822993158D0 /)
      pb(3,:)=(/ 0.00028874763346411D0,-0.00495214157424223D0,-0.00260939452842348D0,&
       1.00727278846920161D0 /)
      pb(4,:)=(/ 0.D0,0.D0,0.D0,1.D0 /)
    case(5) ! EPP5f3
      methname = "EPP5f3"
      order = 6 ! stages = 5
      sigx = 1.5D0
	  pfsal = .true.
      pc = (/ 0.D0, 0.39884436163035D0, 1.69682706703097D0, 2.22894777271564D0, 1D0 /)
      pb(1,:)=(/ 0.D0, 0.D0, 0.D0, 0.D0, 1.D0 /)
      pb(2,:)=(/0.241672476438073D0, 0.666924367890126D0,-0.0885084327153958D0,&
        0.021731059295136D0, 0.158180529092061D0/)
      pb(3,:)=(/0.0350785127485454D0,-0.0487307782453271D0, 0.0296222650530156D0,&
       -0.0810700078774731D0, 1.06510000832124D0/)
      pb(4,:)=(/-0.0850957114286093D0, 0.0590239723048774D0, 0.0590155919885721D0,&
        0.033557017674838D0, 0.933499129460322D0/)
      pb(5,:)=(/0.200970258904312D0, 0.584979404360599D0,-0.0151944433990113D0,&
       -0.0406515692479206D0, 0.269896349382021D0/)
    case (6) ! EPP6j1
      methname = "EPP6j1"
      order = 7 ! stages = 6
      sigx = 1.5D0
      pc = (/6.1182488158460324D-1,1.0734784354567433D0,1.7733348046756701D0,&
        1.9723174701317718D0, 1.4155260278449762D0, 1.D0 /)
      pb(1,:)=(/-2.0180146181687607D-4, 1.6304614802106130D-2,-1.2851544816318155D-2,&
        2.6210658256845634D-3, 3.7912816867901877D-3, 9.9033638396355417D-1/)
      pb(2,:)=(/5.4402447198783130D-5, 2.4618095745274586D-4, 4.1283950478615055D-3,&
        1.0819600004066808D-3, -6.4960842108610951D-3, 1.0009851457579413D0/)
      pb(3,:)=(/-1.3283089416477777D-4, 2.6585250025434079D-4,-3.9985754350926959D-4,&
       -1.4512932179574525D-2, 1.0221144035648493D-2, 1.0045586240813458D0/)
      pb(4,:)=(/1.9579814800354238D-4, -1.4979121212198304D-4, 1.4147303648950299D-4,&
       -1.8142965943514442D-4, -1.7806845742653309D-2, 1.0178007954297175D0/)
      pb(5,:)=(/-7.6319822223659301D-6, 1.8177963233112666D-4,-1.7553812394822423D-4,&
       -4.0614157234655195D-5, 5.3690770730854430D-4, 9.9950509692376555D-1/)
      pb(6,:)=(/0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1.D0 /)
    case (7) ! EPP7f4
	  methname = "EPP7f4"
	  order = 8
	  sigx = 1.5D0
      pfsal = .true.
	  pc=(/0.D0,-1.0887681117114D0, 0.679619857341234D0, -0.801197239808466D0,&
	   1.72075750720492D0, 1.45437989720083D0, 1.D0/)
	  pb(1,:)=(/ 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1.D0 /)
      pb(2,:)=(/-0.0658766764985973D0, 0.100370063156494D0, 0.296507780215623D0,&
 	   -0.0812832310498541D0,-0.041347305457385D0, 0.0129705580263485D0, 0.77865881160737D0/)
      pb(3,:)=(/-0.10465106060643D0, 0.268362172668672D0, 0.194763281593793D0,&
	   -0.0959611977243247D0, 0.0113615560477841D0,-0.0927079454884518D0, 0.818833193508957D0/)
      pb(4,:)=(/-0.0755612754400796D0,-0.0892441539408943D0, 0.153325223233285D0,&
	   0.00681557654247543D0, 0.00854034318702452D0, 0.0255227436206913D0, 0.970601542797498D0/)
      pb(5,:)=(/0.0736022934008542D0, 0.155460799331655D0,-0.192113386065735D0,&
	   -0.0136767232186511D0,-0.0798789827494036D0, 0.0368174391741583D0, 1.01978856012712D0/)
      pb(6,:)=(/-0.0491366058290124D0,-0.098740769585523D0, 0.416720303081099D0,&
	   -0.0290143944255234D0,-0.0455618741851349D0, 0.0126240801815893D0, 0.793109260762506D0/)
      pb(7,:)=(/-0.166222868645607D0, 0.19236125941836D0, 0.421175782008803D0,&
	   -0.0853717083407181D0, 0.0624408466481007D0,-0.18968929236399D0, 0.765305981275052D0/)
    case (8) ! EPP8-d
      methname = "EPP8d"
      order = 9
      sigx=1.4D0
      pc = (/0.26041740957753135D0, 0.52923626823623069D0, 1.54653689839871537D0,&
          1.77514379674033294D0, 0.090107569010498D0, 0.19995116126552217D0,&
          1.32615523512850911D0, 1.D0 /)
      pb(1,:)=(/ 0.00218540643609441D0, 0.00964151773849017D0, 0.02740905027639759D0,&
         -0.00757481375423608D0, -0.00627931764830473D0, -0.00651613228833201D0,&
         -0.00103882981570947D0,  0.98217311905560012D0 /)
      pb(2,:)=(/ 0.0026911855247508D0, -0.00999004293842317D0, 0.01299877362586784D0,&
         -0.00362117297559769D0, -0.01393348294082885D0, 0.00206411701387576D0,&
         -0.00794032918151776D0,  1.01773095187187307D0 /)
      pb(3,:)=(/ 0.00199608410884834D0, -0.01002762418952864D0, 0.02034106943279384D0,&
         -0.00932460631409766D0, -0.00937443352246159D0, 0.00649579387949097D0,&
         -0.02592014846435581D0,  1.02581386506931055D0  /)
      pb(4,:)=(/ 0.00103036221311367D0, -0.01190896158535306D0, 0.02063699940945861D0,&
         -0.00672959115904196D0, -0.00583850891833535D0, 0.03351758481419382D0,&
         -0.03846471330760932D0,  1.0077568285335736D0  /)
      pb(5,:)=(/ 0.00148756342489522D0, -0.01047547599327585D0, 0.0208838878315454D0,&
         -0.00663219960933854D0, -0.00563964345840616D0, 0.02716167598460901D0,&
         -0.00664272798306554D0, 0.97985691980303647D0  /)
      pb(6,:)=(/ 0.00205269876381091D0, -0.00960860565429818D0, 0.0205427263285296D0,&
         -0.00625381166734364D0, -0.00577560681234781D0, 0.01293837851580094D0,&
          0.00660793064341287D0,  0.97949628988243531D0 /)
      pb(7,:)=(/ 0.00184988219096286D0, -0.01079724930569055D0, 0.02051636130170474D0,&
         -0.00632847064898024D0, -0.00550013664675566D0, 0.01390490619566548D0,&
         -0.01323688831455224D0,  0.99959159522764561D0 /)
      pb(8,:)=(/ 0.00186731024859287D0, -0.01061165690522434D0, 0.02042910717986102D0,&
         -0.00630049368165016D0, -0.00568571957918147D0, 0.01348102696866063D0,&
         -0.01331088571679290D0,  1.00013131148573434D0 /)
    case (9) ! EPP9f2
      methname = "EPP9f2"
      order = 10
      sigx = 1.4D0
	  pfsal = .true.
      pc = (/0.D0,-1.45065446034357D0,-1.14222878247125D0, 0.380810721807532D0,&
       1.40211112656114D0,-0.681879393558872D0, 1.8663786044472D0,&
       1.75672620008449D0,1.D0/)
      pb(1,:)=(/ 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 0.D0, 1.D0 /)
      pb(2,:)=(/-0.533252224247118D0,-0.0975501226629606D0, 0.116672069301833D0,&
       0.44697499897106D0, 0.0340845693880289D0, 0.135289891492841D0,&
       0.00622389138487414D0, 0.0325928274134378D0, 0.858964098958003D0/)
      pb(3,:)=(/-0.301652204208198D0, 0.0520009020704396D0,-0.0464702777791137D0,&
       0.280682940598104D0, 0.0671718625949898D0, 0.135006432610688D0,&
       -0.046367639817997D0, 0.0345253672919836D0, 0.825102616639103D0/)
      pb(4,:)=(/0.297011977623876D0,-0.162918534218843D0, 0.0861949381787168D0,&
       -0.22497850740489D0,-0.160217066458187D0,-0.116187126867031D0,&
       -0.0279445699477793D0, 0.0483634067976424D0, 1.26067548229649D0/)
      pb(5,:)=(/-0.0700215072654026D0,-0.0871778303108089D0, 0.018588359758533D0,&
       0.0493308059368063D0, 0.0386560715601355D0, 0.00206062231467469D0,&
       0.0203705003455984D0,-0.0164980878257943D0, 1.04469106548626D0/)
      pb(6,:)=(/-0.377028842757212D0, 0.220063630217656D0,-0.129814104899979D0,&
       0.27827383620039D0, 0.186943412912638D0, 0.0936142543473008D0,&
       0.04468074060207D0,-0.0421130346891518D0, 0.725380108066289D0/)
      pb(7,:)=(/-0.210835638176618D0,-0.0796361355345654D0, 0.0735641672208091D0,&
       0.241807873963969D0,-0.103620229495115D0,-0.0163192383349393D0,&
       -0.00975174063171978D0, 0.082602741614847D0, 1.02218819937333D0/)
      pb(8,:)=(/0.0761288961852642D0,-0.0447817877333645D0,-0.0309993817579791D0,&
       -0.00589722673837671D0,-0.0760908310760073D0,-0.0992689288448062D0,&
       -0.00659612005493727D0, 0.0536747338124468D0, 1.13383064620776D0/)
      pb(9,:)=(/0.0451553425438244D0,-0.162203290825595D0, 0.0723244975922378D0,&
       -0.0184649776742533D0,-0.0998398851224182D0,-0.0683944910444858D0,&
       -0.0184731703513451D0, 0.0570903861232346D0, 1.1928055887588D0/)
    case default !
      methname = "error"
      pstage = -1
    end select
!    compute other coefficients
    do i = 1,stages ! row sums of B must be 1 exactly
      s = 1.D0
      do j = 1,stages
        s = s-pb(i,j)
      end do
!	  if (abs(s)>1D-12) print *, "B!!", i, s
      pb(i,stages) = pb(i,stages)+s
    end do
    ! compute matrix PCV =diag(C)*VdM
    do i = 1,stages
      pcv(i,1) = pc(i)
      do j = 2,stages
        pcv(i,j) = pc(i)*pcv(i,j-1)
      end do
    end do
   ! compute matrix PA0, ppv used as temporary memory
    ppv = -pb
    ppv(:,stages) = ppv(:,stages)+1.D0 ! now \ones*e_s^T-B
    pa0 = matmul(ppv,pcv)
    do j=2,stages  ! scale columns
      pa0(:,j) = pa0(:,j)/j
    end do
    call vandm(stages,pc,pa0) ! now PA0=(\ones*e_s^T-B)*C*V*(V*D)^(-1)
    ! compute matrix PPV:
    ppv(1,:) = 1.D0 ! start Pascal matrix
    do i = 2,stages
      ppv(i,1) = 0.D0
      do j = 1,stages-1
        ppv(i,j+1) = ppv(i,j)+ppv(i-1,j)
      end do
    end do ! end Pascal matrix
    do i = 2,stages ! scale rows
      ppv(i,:) = ppv(i,:)/i
    end do
    call vandm(stages,pc,ppv) ! now PPV=D^(-1)*P*V^(-1)
    pe = ppv(stages,:) ! error estimate with last row of PPV
    ! set default parameters:
    stpci = .true. ! step control on
	pcout = .false.! dense output off
    maxsteps = 1000000 ! maximal number of time steps
	phs = -1D0 ! initial stepsize (none => own estimate)
    hmin = 1.D-10 ! minimal allowed stepsize
    atoli = 1.D-5; rtoli=atoli ! absolut and relative tolerances
  end subroutine ppsetcoeff

  !-- PPGETSTATS returns integration statistics
  subroutine ppgetstats(nsteps,nrej,nfcn)
    implicit none
    integer, intent(out) :: nsteps, nrej, nfcn
    integer :: i
    nsteps = isteps
    nrej = irej
    nfcn = ifcn
  end subroutine ppgetstats

!  -- PPFREECOEFF releases all memory
  subroutine ppfreecoeff
    deallocate(pb,pa0,pcv,ppv,pc,pe)
 	if (allocated(pyy)) deallocate(pyy)
 end subroutine ppfreecoeff

!  -- inisteps computes initial stepsize [dopri/Hairer-Wanner]
  subroutine inisteps(n,fcn,t,y,fv1,hs,par)
    implicit none
    interface
      subroutine fcn(t,u,udot,par)
      real(8), intent(in) :: t
      real(8), dimension(:), intent(in) :: u
      real(8), dimension(:), intent(out) :: udot
      real(8), dimension(:), intent(in) :: par
      end subroutine fcn
    end interface
    integer, intent(in) :: n
    real(8), dimension(:), intent(in) :: y, fv1, par
    real(8) :: t,hs, h,h1,rc,dnf,dny,der2,der12
    real(8), dimension(n) :: y2, f2
    integer :: id
    dnf = 0.D0; dny = 0.D0
    do id = 1,n
      rc = atoli+rtoli*abs(y(id))
      dnf = dnf +(fv1(id)/rc)**2
      dny = dny +(y(id)/rc)**2
    end do
    if (min(dnf,dny)<1D-10) then
      h = 1.D-6
    else
      h =1.D-2*sqrt(dny/dnf)
    end if
    h = min(h,hmax)
    y2 = y+h*fv1 ! explicit Euler step
    call fcn(t+h,y2,f2,par)
    der2 = 0.D0
    do id = 1,n
      rc = atoli+rtoli*abs(y(id))
      der2 = der2 +((f2(id)-fv1(id))/rc)**2
    end do
    der2 = sqrt(der2)/h ! estimate for second derivative
    der12 = max(der2,sqrt(dnf))
    if (der12<=1.D-15) then
      h1 = max(1.D-6,h*1.D-3)
    else
      h1 = (1.D-2/der12)**(1.D0/order)
    end if
    hs = min(1D2*h,h1,hmax)
  end subroutine inisteps

end module IVPEPP