!> This module packages odeint from numerical recipes.
module mod_odeint
  use mod_comm_lib, only: mpistop
  implicit none
  integer                  :: MAXSTP, NMAX, KMAXX
  PARAMETER (MAXSTP=10000,NMAX=50,KMAXX=200)
  double precision         :: TINY
  PARAMETER (TINY=1.d-30)

contains

  subroutine odeint(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,derivs,rkqs,ierror)

    INTEGER nbad,nok,nvar
    double precision eps,h1,hmin,x1,x2,ystart(nvar)
    INTEGER i,nstp
    double precision h,hdid,hnext,x,xsav,dydx(NMAX),y(NMAX),yscal(NMAX)
    integer kmax, kount
    double precision dxsav, yp(NMAX,KMAXX), xp(KMAXX)
    ! Can be 0: success, 1: hmin too small, 2: MAXSTP exceeded
    integer, intent(out) :: ierror

    EXTERNAL derivs,rkqs

    COMMON /path/ kmax,kount,dxsav,xp,yp

    x     = x1
    h     = sign(h1,x2-x1)
    nok   = 0
    nbad  = 0
    kount = 0
    xsav  = 0.d0

    do i=1,nvar
       y(i)=ystart(i)
    end do

    if (kmax.gt.0) xsav=x-2.d0*dxsav

    do nstp=1,MAXSTP
       call derivs(x,y,dydx)
       do i=1,nvar
          yscal(i)=abs(y(i))+abs(h*dydx(i))+TINY
       end do

       if(kmax.gt.0)then
          if(abs(x-xsav).gt.abs(dxsav)) then
             if(kount.lt.kmax-1)then
                kount=kount+1
                xp(kount)=x
                do i=1,nvar
                   yp(i,kount)=y(i)
                end do
                xsav=x
             end if
          end if
       end if

       if((x+h-x2)*(x+h-x1).gt.0.d0) h=x2-x

       if (any(y(1:nvar) .ne. y(1:nvar)) .or. any(dydx(1:nvar) .ne. &
          dydx(1:nvar)) .or. x .ne. x) then
         write(*,*) "ODEINT BEFORE CALL RKQS: NaN in x, dydx, or y!"
         write(*,*) "x",x
         write(*,*) "y",y
         write(*,*) "dydx",dydx
         write(*,*) "exiting..."
         ierror = 2
         return
       end if

       call rkqs(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs)

       if (any(y(1:nvar) .ne. y(1:nvar)) .or. any(dydx(1:nvar) .ne. &
          dydx(1:nvar)) .or. x .ne. x) then
         write(*,*) "ODEINT AFTER CALL RKQS: NaN in x, dydx, or y!"
         write(*,*) "x",x
         write(*,*) "y",y
         write(*,*) "dydx",dydx
         write(*,*) "exiting..."
         ierror = 2
         return
       end if

       if(hdid.eq.h)then
          nok=nok+1
       else
          nbad=nbad+1
       end if

       if((x-x2)*(x2-x1).ge.0.d0)then
          do i=1,nvar
             ystart(i)=y(i)
          end do
          if(kmax.ne.0)then
             kount=kount+1
             xp(kount)=x
             do i=1,nvar
                yp(i,kount)=y(i)
             end do
          end if

          ierror = 0
          return
       end if

       if(abs(hnext).lt.hmin) then
          ierror = 1
       end if

       h=hnext
    end do

    ierror = 2
  end subroutine odeint

      SUBROUTINE rkqs(y,dydx,n,x,htry,eps,yscal,hdid,hnext,derivs)

      INTEGER n,NMAX
      double precision eps,hdid,hnext,htry,x,dydx(n),y(n),yscal(n)
      EXTERNAL derivs
      PARAMETER (NMAX=50)
!     CU USES derivs,rkck
      INTEGER i
      double precision errmax,h,xnew,yerr(NMAX),ytemp(NMAX),SAFETY,PGROW,&
         PSHRNK,ERRCON
      PARAMETER (SAFETY=0.9d0,PGROW=-.2d0,PSHRNK=-.25d0,ERRCON=1.89d-4)

      h=htry

       if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n)) .or. x &
          .ne. x) then
         write(*,*) "RKQS BEFORE CALL RKCK: NaN in x, dydx, or y!"
         write(*,*) "x",x
         write(*,*) "y",y
         write(*,*) "dydx",dydx
         write(*,*) "exiting..."
         return
       end if

 1    call rkck(y,dydx,n,x,h,ytemp,yerr,derivs)
       if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n)) .or. x &
          .ne. x) then
         write(*,*) "RKQS AFTER CALL RKCK: NaN in x, dydx, or y!"
         write(*,*) "x",x
         write(*,*) "y",y
         write(*,*) "dydx",dydx
         write(*,*) "exiting..."
         return
       end if

      errmax=0.d0
      do 11 i=1,n
         errmax=max(errmax,abs(yerr(i)/yscal(i)))
 11   continue
      errmax=errmax/eps
      if(errmax.gt.1.d0)then
         h=SAFETY*h*(errmax**PSHRNK)
         if(h.lt.0.1d0*h)then
            h=.1d0*h
         end if
         xnew=x+h
         if(xnew.eq.x) call stop('stepsize underflow in rkqs')
         goto 1
      else
         if(errmax.gt.ERRCON)then
            hnext=SAFETY*h*(errmax**PGROW)
         else
            hnext=5.d0*h
         end if
         hdid=h
         x=x+h
         do 12 i=1,n
            y(i)=ytemp(i)
 12      continue
         return
      end if
      END subroutine rkqs
!C (C) Copr. 1986-92 Numerical Recipes Software Vs1&v%1jw#<?4210(9p#.

      SUBROUTINE rkck(y,dydx,n,x,h,yout,yerr,derivs)

      INTEGER n,NMAX
      double precision  h,x,dydx(n),y(n),yerr(n),yout(n)
      EXTERNAL derivs
      PARAMETER (NMAX=50)
!CU    USES derivs
      INTEGER i
      double precision ak2(NMAX),ak3(NMAX),ak4(NMAX),ak5(NMAX),ak6(NMAX),&
         ytemp(NMAX),A2,A3,A4,A5,A6,B21,B31,B32,B41,B42,B43,B51, B52,B53,B54,&
         B61,B62,B63,B64,B65,C1,C3,C4,C6,DC1,DC3,DC4,DC5,DC6
      PARAMETER (A2=.2d0,A3=.3d0,A4=.6d0,A5=1.d0, A6=.875d0,B21=.2d0,&
         B31=3.d0/40.d0,B32=9.d0/40.d0,B41=.3d0,B42=-.9d0,B43=1.2d0,&
         B51=-11.d0/54.d0,B52=2.5d0,B53=-70.d0/27.d0,B54=35.d0/27.d0,&
         B61=1631.d0/55296.d0, B62=175.d0/512.d0, B63=575.d0/13824.d0,&
         B64=44275.d0/110592.d0, B65=253.d0/4096.d0,C1=37.d0/378.d0,&
         C3=250.d0/621.d0,C4=125.d0/594.d0,C6=512.d0/1771.d0,&
         DC1=C1-2825.d0/27648.d0,DC3=C3-18575.d0/48384.d0,&
         DC4=C4-13525.d0/55296.d0,DC5=-277.d0/14336.d0,DC6=C6-.25d0)
      if (any(y(1:n) .ne. y(1:n)) .or. any(dydx(1:n) .ne. dydx(1:n))) then
        write(*,*) "NaNs IN RKCK, STEP 0!"
        write(*,*) "y0",y(1:n)
        write(*,*) "derivs",dydx(1:n)
        write(*,*) "ABORTING..."
        call mpistop("")
      end if
      do 11 i=1,n
         ytemp(i)=y(i)+B21*h*dydx(i)
 11   continue
      call derivs(x+A2*h,ytemp,ak2)
      if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak2(1:n) .ne. ak2(1:n))) &
         then
        write(*,*) "NaNs IN RKCK, STEP 1!"
        write(*,*) "y1",ytemp(1:n)
        write(*,*) "derivs",ak2(1:n)
        write(*,*) "ABORTING..."
        call mpistop("")
      end if
      do 12 i=1,n
         ytemp(i)=y(i)+h*(B31*dydx(i)+B32*ak2(i))
 12   continue
      call derivs(x+A3*h,ytemp,ak3)
      if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak3(1:n) .ne. ak3(1:n))) &
         then
        write(*,*) "NaNs IN RKCK, STEP 2!"
        write(*,*) "y2",ytemp(1:n)
        write(*,*) "derivs",ak3(1:n)
        write(*,*) "ABORTING..."
        call mpistop("")
      end if
      do 13 i=1,n
         ytemp(i)=y(i)+h*(B41*dydx(i)+B42*ak2(i)+B43*ak3(i))
 13   continue
      call derivs(x+A4*h,ytemp,ak4)
      if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak4(1:n) .ne. ak4(1:n))) &
         then
        write(*,*) "NaNs IN RKCK, STEP 3!"
        write(*,*) "y3",ytemp(1:n)
        write(*,*) "derivs",ak4(1:n)
        write(*,*) "ABORTING..."
        call mpistop("")
      end if
      do 14 i=1,n
         ytemp(i)=y(i)+h*(B51*dydx(i)+B52*ak2(i)+B53*ak3(i)+B54*ak4(i))
 14   continue
      call derivs(x+A5*h,ytemp,ak5)
      if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak5(1:n) .ne. ak5(1:n))) &
         then
        write(*,*) "NaNs IN RKCK, STEP 4!"
        write(*,*) "y4",ytemp(1:n)
        write(*,*) "derivs",ak5(1:n)
        write(*,*) "ABORTING..."
        call mpistop("")
      end if
      do 15 i=1,n
         ytemp(i)=y(i)+h*(B61*dydx(i)+B62*ak2(i)+B63*ak3(i)+B64*ak4(i)+&
            B65*ak5(i))
 15   continue
      call derivs(x+A6*h,ytemp,ak6)
      if (any(ytemp(1:n) .ne. ytemp(1:n)) .or. any(ak6(1:n) .ne. ak6(1:n))) &
         then
        write(*,*) "NaNs IN RKCK, STEP 1!"
        write(*,*) "y15",ytemp(1:n)
        write(*,*) "derivs",ak6(1:n)
        write(*,*) "ABORTING..."
        call mpistop("")
      end if
      do 16 i=1,n
         yout(i)=y(i)+h*(C1*dydx(i)+C3*ak3(i)+C4*ak4(i)+C6*ak6(i))
 16   continue
      do 17 i=1,n
         yerr(i)=h*(DC1*dydx(i)+DC3*ak3(i)+DC4*ak4(i)+DC5*ak5(i)+DC6*ak6(i))
 17   continue
      return
      END subroutine rkck
!C  (C) Copr. 1986-92 Numerical Recipes Software Vs1&v%1jw#<?4210(9p#.
  subroutine stop(text)
    character(len=*), intent(in)   :: text
    print*, text
    STOP
  end subroutine stop

end module mod_odeint