! This file ivp_pmain.f90 is the driver program from the EPPEER-package for explicit
! parallel peer methods by B.A.Schmitt and R.Weiner available at
! www.mathematik.uni-marburg.de\~schmitt\peer 
! It solves the ODE from the module ODEPROB for different tolerances and writes
! times and accuracies in a log file; the specific peer method is chosen in line 39
program solve_ivp
  use omp_lib ! OpenMP-Library for parallel execution
  use odeprob ! definition of ODE
  use ivpepp  ! module with eppeer subroutine

  implicit none
  external psolout ! subroutine for dense output, defined below
  integer :: n, ito, stages, irep, nsteps,nrej,nfcn, mnr, nthr
  integer :: i,j
  real(8) :: tm,te, tol, hs, err, ct0,ctm, ot0,otm
  real(8), dimension(:), allocatable :: ym, ue, cpar
  character(16) :: mthn
  character(32) :: logname ! name of log file
  character(9) :: npr = "123456789"
  n = nprob ! ODE problem dimension by global variable
  allocate(ym(n),ue(n),cpar(nparm))
  cpar(1)=1.D0
!--- get solution for error computation:
  call inivals(tm,te,ue,cpar)
  if (exsol) then ! get exact solution ue
    call solution(te,ue,cpar)
  else     ! compute ue with high accuracy
    print "(a)", "-solution approximated"
    mnr = 5 ! choose dopri5
    call rksetcoeff(mnr,mthn)
    call rksetparam(1D-14,1D-14) ! very sharp tolerances
    hs = 0.D0
    call rkmeth(n,fcn,tm,te,ue,hs,irep,cpar) ! reference solution with Runge-Kutta method
    call rkfreecoeff
    if (irep/=0) print "(a)", "*** bad reference solution"
  end if
!--- end test data, initialize method:
  nthr = omp_get_max_threads()
  mnr = 5 ! choose peer method, last digit=#stages
  call ppsetcoeff(mnr,stages,mthn) ! initialize method
  if (stages<1) print *, "peer method ", mnr," not available"
  logname = npr(nthr:nthr)//odename//trim(mthn)
  print "(3a,i2)", "Program IVP, method ",trim(mthn),", Procs=",nthr
  print "(3a,i8,a,i8)", "problem ",odename,", dimension n=",n,", parameters=",nparm
  print "(2a)", ".. writing logfile ", logname
  open(4,file="pout.txt") ! dense output
  open(3,file=trim(logname))
  write (3, "(3a,I2,3a,I8)") "# method ",trim(mthn),", Procs=",nthr, &
    ", problem: ",odename,", dim=", n
  tol = 1.D-1
  do ito =1,11  ! run with different tolerances or stepsizes
    tol = tol*1.D-1
    call cpu_time(ct0)
    ot0 = omp_get_wtime() ! OpenMP start time
	call ppsetacc(tol,tol) ! (re-)set parameters, hs omitted
    call ppsetjob(contout=(ito==3))! dense output for ito==...
    call inivals(tm,te,ym,cpar)   ! get initial values
!    call eppeer(fcn,tm,te,ym,cpar) ! <<---- integration, no dense output
    call eppeer(fcn,tm,te,ym,cpar,psolout) ! <<--------- integration with dense output
    call cpu_time(ctm); ctm = ctm-ct0
    otm = omp_get_wtime()-ot0  ! OpenMP thread time
    irep = ppreport()
	call ppgetstats(nsteps,nrej,nfcn)
    if (irep==0) then ! computation successful
      ym = (ym-ue)**2
      err = sqrt(sum(ym)/n)
      print "(a, G9.2, 3G12.5)", " tol, err, otime, cpu ", tol, err, otm, ctm
      write(3,"(E9.2,2E14.5,I9,E14.5)") tol, err, otm, nfcn, ctm !!otm
    else ! print error message
      print "(a,i4)", "** error in ppmeth, nr:", irep
    end if
    print "(a,2i5,i9)", " steps,rej,nfcn:", nsteps,nrej,nfcn
  end do
  close(3); close(4)
  deallocate(ym,ue,cpar)
  call ppfreecoeff ! release eppeer memory
end program solve_ivp

! subroutine for dense output to file "pout.txt", called by eppeer
  subroutine psolout(n,tsol,tnew) 
    use ivpepp
    implicit none
    integer, intent(in) :: n
    real(8) :: tsol
    real(8), intent(in) :: tnew
    real(8), dimension(n) :: ycon
    do while (tsol<=tnew)
      call ppcont(tsol,tnew,ycon)
      write(4,"(F12.3,2E14.5)") tsol,ycon(1),ycon(2)
      tsol = tsol+0.05D0
    end do
  end subroutine psolout