!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
! Added: ODEINT2 (see below) 8/7/2001 DCS.
!
!---------------------------------------------------------------------

module custom_odeint_mod

  USE precision_def
  USE nrtype, except => dp

  type odeint_com_struct
    REAL(rp) :: dxsav 
    REAL(rp), DIMENSION(:), POINTER :: xp
    REAL(rp), DIMENSION(:,:), POINTER :: yp
    INTEGER(I4B) :: nok,nbad,kount
    LOGICAL(LGT) :: save_steps=.false.
    LOGICAL too_many_steps
  end type

  type (odeint_com_struct), save :: ode_com
  
contains

!---------------------------------------------------------------------
!---------------------------------------------------------------------
!---------------------------------------------------------------------

  SUBROUTINE custom_odeint(ystart,x1,x2,eps,h1,hmin,derivs,custom_rkqs)
  USE nrutil, ONLY : nrerror,reallocate
  USE ode_path
  IMPLICIT NONE
  REAL(rp), DIMENSION(:), INTENT(INOUT) :: ystart
  REAL(rp), INTENT(IN) :: x1,x2,eps,h1,hmin
  INTERFACE
    SUBROUTINE derivs(x,y,dydx)
    USE precision_def
    USE nrtype
    IMPLICIT NONE
    REAL(rp), INTENT(IN) :: x
    REAL(rp), DIMENSION(:), INTENT(IN) :: y
    REAL(rp), DIMENSION(:), INTENT(OUT) :: dydx
    END SUBROUTINE derivs
!BL
    SUBROUTINE custom_rkqs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs)
    USE precision_def
    USE nrtype
    IMPLICIT NONE
    REAL(rp), DIMENSION(:), INTENT(INOUT) :: y
    REAL(rp), DIMENSION(:), INTENT(IN) :: dydx,yscal
    REAL(rp), INTENT(INOUT) :: x
    REAL(rp), INTENT(IN) :: htry,eps
    REAL(rp), INTENT(OUT) :: hdid,hnext
    INTERFACE
      SUBROUTINE derivs(x,y,dydx)     
      USE precision_def
      USE nrtype
      IMPLICIT NONE
      REAL(rp), INTENT(IN) :: x
      REAL(rp), DIMENSION(:), INTENT(IN) :: y
      REAL(rp), DIMENSION(:), INTENT(OUT) :: dydx
      END SUBROUTINE derivs
    END INTERFACE
    END SUBROUTINE 
  END INTERFACE
  REAL(rp), PARAMETER :: TINY=1.0e-30_rp
  REAL(rp) :: h,hdid,hnext,x,xsav
  REAL(rp), DIMENSION(size(ystart)) :: dydx,y,yscal
  INTEGER(I4B), PARAMETER :: MAXSTP=10000
  INTEGER(I4B) :: nstp
  x=x1
  h=sign(h1,x2-x1)
  nok=0
  nbad=0
  kount=0
  print *,'ystart ',ystart
  print *,'y ',y
  y(:)=ystart(:)
  if (ode_com%save_steps) then
    xsav=x-2.0_rp*dxsav
    nullify(ode_com%xp, ode_com%yp)
    allocate(ode_com%xp(256))
    allocate(ode_com%yp(size(ystart),size(ode_com%xp)))
  end if
  do nstp=1,MAXSTP
    call derivs(x,y,dydx)
    yscal(:)=abs(y(:))+abs(h*dydx(:))+TINY
    if (ode_com%save_steps .and. (abs(x-xsav) > abs(ode_com%dxsav))) &
      call save_a_step
    if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x
    call custom_rkqs (y,dydx,x,h,eps,yscal,hdid,hnext,derivs)
    if (hdid == h) then
      nok=nok+1
    else
      nbad=nbad+1
    end if
    if ((x-x2)*(x2-x1) >= 0.0) then
      ystart(:)=y(:)
      if (ode_com%save_steps) call save_a_step
      RETURN
    end if
    if (abs(hnext) < hmin)&
      call nrerror('stepsize smaller than minimum in odeint')
    h=hnext
  end do
    call nrerror('too many steps in odeint')
CONTAINS

  SUBROUTINE save_a_step
  ode_com%kount=ode_com%kount+1
  if (ode_com%kount > size(ode_com%xp)) then
    ode_com%xp=>reallocate(ode_com%xp,2*size(ode_com%xp))
    ode_com%yp=>reallocate(ode_com%yp,size(ode_com%yp,1),size(ode_com%xp))
  end if
  ode_com%xp(ode_com%kount)=x
  ode_com%yp(:,ode_com%kount)=y(:)
  xsav=x
  END SUBROUTINE save_a_step

END SUBROUTINE 

!---------------------------------------------------------------------
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!+
! Subroutine custom_odeint2 (ystart, x1, x2, rel_eps, abs_eps, h1, hmin, &
!                   derivs, custom_rkqs)
!
! This is essentually odeint with a better error limit algorithm.
! Notice that this routine has one more argument (abs_eps) than odeint.
! Also rel_eps (essentually equivalent to eps in odeint) is scalled by the
! step size to to able to relate it to the final accuracy.
!
! Input: [See the odeint writeup for the arguments not shown]
!   rel_eps -- Real(rp): Same as eps for odeint scalled by sqrt(h/(x2-x1))
!               where h is the step size for the current step. rel_eps
!               sets the %error of the result
!   abs_eps -- Real(rp): Sets the absolute error of the result
!
! Essentually (assuming random errors) one of these conditions holds:
!             %error in tracking < rel_eps  
! or
!     absolute error in tracking < abs_eps
!
                           
SUBROUTINE custom_odeint2 (ystart, x1, x2, rel_eps, abs_eps, h1, hmin, derivs, custom_rkqs)

  USE nrutil, ONLY : nrerror,reallocate
  USE ode_path
  IMPLICIT NONE
  REAL(rp), DIMENSION(:), INTENT(INOUT) :: ystart
  REAL(rp), INTENT(IN) :: x1,x2,rel_eps,h1,hmin
  real(rp), intent(in) :: abs_eps

  INTERFACE
    SUBROUTINE derivs(x,y,dydx)
      USE nrtype
      IMPLICIT NONE
      REAL(DP), INTENT(IN) :: x
      REAL(DP), DIMENSION(:), INTENT(IN) :: y
      REAL(DP), DIMENSION(:), INTENT(OUT) :: dydx
    END SUBROUTINE derivs
    SUBROUTINE custom_rkqs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs)
      USE nrtype
      IMPLICIT NONE
      REAL(DP), DIMENSION(:), INTENT(INOUT) :: y
      REAL(DP), DIMENSION(:), INTENT(IN) :: dydx,yscal
      REAL(DP), INTENT(INOUT) :: x
      REAL(DP), INTENT(IN) :: htry,eps
      REAL(DP), INTENT(OUT) :: hdid,hnext
      INTERFACE
        SUBROUTINE derivs(x,y,dydx)
          USE nrtype
          IMPLICIT NONE
          REAL(DP), INTENT(IN) :: x
          REAL(DP), DIMENSION(:), INTENT(IN) :: y
          REAL(DP), DIMENSION(:), INTENT(OUT) :: dydx
        END SUBROUTINE derivs
      END INTERFACE
    END SUBROUTINE 
  END INTERFACE

  REAL(DP), PARAMETER :: TINY=1.0e-30_rp
  REAL(DP) :: h,hdid,hnext,x,xsav
  REAL(DP), DIMENSION(size(ystart)) :: dydx,y,yscal

  real(DP) sqrt_N, rel_eps_eff, abs_eps_eff            ! added
  INTEGER(I4B), PARAMETER :: MAXSTP=10000
  INTEGER(I4B) :: nstp
  ode_com%too_many_steps = .false.
!

  x=x1
  h=sign(h1,x2-x1)
  nok=0
  nbad=0
  kount=0
  y(:)=ystart(:)
  if (ode_com%save_steps) then
    xsav=x-2.0_rp*dxsav
    nullify(ode_com%xp,ode_com%yp)
    allocate(ode_com%xp(256))
    allocate(ode_com%yp(size(ystart),size(ode_com%xp)))
  end if
  do nstp=1,MAXSTP
    call derivs(x,y,dydx)

! this is the section of executable code that has been changed.

    sqrt_N = sqrt(abs((x2-x1)/h))  ! number of steps we would take with this h
    rel_eps_eff = rel_eps / sqrt_N
    abs_eps_eff = abs_eps / sqrt_N
    yscal(:) = abs(y(:)) + abs(h*dydx(:)) + abs_eps_eff/rel_eps_eff + TINY

!

    if (ode_com%save_steps .and. (abs(x-xsav) > abs(ode_com%dxsav))) &
      call save_a_step
    if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x

    call custom_rkqs(y,dydx,x,h,rel_eps_eff,yscal,hdid,hnext,derivs)  ! changed

    IF(h == 0.)THEN    !dlr 9-19-01
      ode_com%too_many_steps = .true.
      print *,'stepsize underflow in custom_rkqs'
      return
    ENDIF
    if (hdid == h) then
      nok=nok+1
    else
      nbad=nbad+1
    end if
    if ((x-x2)*(x2-x1) >= 0.0) then
      ystart(:)=y(:)
      if (ode_com%save_steps) call save_a_step
      RETURN
    end if
    if (abs(hnext) < hmin)&
      call nrerror('stepsize smaller than minimum in odeint')
    h=hnext
  end do
  ode_com%too_many_steps = .true.  !dlr   9-9-01
  return
  call nrerror('too many steps in odeint')

  CONTAINS

  SUBROUTINE save_a_step
  ode_com%kount=ode_com%kount+1
  if (ode_com%kount > size(ode_com%xp)) then
    ode_com%xp=>reallocate(ode_com%xp,2*size(ode_com%xp))
    ode_com%yp=>reallocate(ode_com%yp,size(ode_com%yp,1),size(ode_com%xp))
  end if
  ode_com%xp(ode_com%kount)=x
  ode_com%yp(:,ode_com%kount)=y(:)
  xsav=x
  END SUBROUTINE save_a_step

END SUBROUTINE 

	SUBROUTINE custom_rkqs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs)
	USE nrtype; USE nrutil, ONLY : assert_eq,nrerror
	IMPLICIT NONE
	REAL(DP), DIMENSION(:), INTENT(INOUT) :: y
	REAL(DP), DIMENSION(:), INTENT(IN) :: dydx,yscal
	REAL(DP), INTENT(INOUT) :: x
	REAL(DP), INTENT(IN) :: htry,eps
	REAL(DP), INTENT(OUT) :: hdid,hnext
	INTERFACE
		SUBROUTINE derivs(x,y,dydx)
		USE nrtype
		IMPLICIT NONE
		REAL(DP), INTENT(IN) :: x
		REAL(DP), DIMENSION(:), INTENT(IN) :: y
		REAL(DP), DIMENSION(:), INTENT(OUT) :: dydx
		END SUBROUTINE derivs
	END INTERFACE
	INTEGER(I4B) :: ndum
	REAL(DP) :: errmax,h,htemp,xnew
	REAL(DP), DIMENSION(size(y)) :: yerr,ytemp
	REAL(DP), PARAMETER :: SAFETY=0.9_rp,PGROW=-0.2_rp,PSHRNK=-0.25_rp,&
		ERRCON=1.89e-4
	ndum=assert_eq(size(y),size(dydx),size(yscal),'custom_rkqs')
	h=htry
	do
		call custom_rkck(y,dydx,x,h,ytemp,yerr,derivs)
		errmax=maxval(abs(yerr(:)/yscal(:)))/eps
		if (errmax <= 1.0) exit
		htemp=SAFETY*h*(errmax**PSHRNK)
		h=sign(max(abs(htemp),0.1_rp*abs(h)),h)
		xnew=x+h
		if (xnew == x) return       !dlr 9-19-01
		if (xnew == x) call nrerror('stepsize underflow in custom_rkqs')
	end do
	if (errmax > ERRCON) then
		hnext=SAFETY*h*(errmax**PGROW)
	else
		hnext=5.0_rp*h
	end if
	hdid=h
	x=x+h
	y(:)=ytemp(:)
	END SUBROUTINE 


	SUBROUTINE custom_rkck(y,dydx,x,h,yout,yerr,derivs)
	USE nrtype; USE nrutil, ONLY : assert_eq
	IMPLICIT NONE
	REAL(DP), DIMENSION(:), INTENT(IN) :: y,dydx
	REAL(DP), INTENT(IN) :: x,h
	REAL(DP), DIMENSION(:), INTENT(OUT) :: yout,yerr
	INTERFACE
		SUBROUTINE derivs(x,y,dydx)
		USE nrtype
		IMPLICIT NONE
		REAL(DP), INTENT(IN) :: x
		REAL(DP), DIMENSION(:), INTENT(IN) :: y
		REAL(DP), DIMENSION(:), INTENT(OUT) :: dydx
		END SUBROUTINE derivs
	END INTERFACE
	INTEGER(I4B) :: ndum
	REAL(DP), DIMENSION(size(y)) :: ak2,ak3,ak4,ak5,ak6,ytemp
	REAL(DP), PARAMETER :: A2=0.2_rp,A3=0.3_rp,A4=0.6_rp,A5=1.0_rp,&
		A6=0.875_rp,B21=0.2_rp,B31=3.0_rp/40.0_rp,B32=9.0_rp/40.0_rp,&
		B41=0.3_rp,B42=-0.9_rp,B43=1.2_rp,B51=-11.0_rp/54.0_rp,&
		B52=2.5_rp,B53=-70.0_rp/27.0_rp,B54=35.0_rp/27.0_rp,&
		B61=1631.0_rp/55296.0_rp,B62=175.0_rp/512.0_rp,&
		B63=575.0_rp/13824.0_rp,B64=44275.0_rp/110592.0_rp,&
		B65=253.0_rp/4096.0_rp,C1=37.0_rp/378.0_rp,&
		C3=250.0_rp/621.0_rp,C4=125.0_rp/594.0_rp,&
		C6=512.0_rp/1771.0_rp,DC1=C1-2825.0_rp/27648.0_rp,&
		DC3=C3-18575.0_rp/48384.0_rp,DC4=C4-13525.0_rp/55296.0_rp,&
		DC5=-277.0_rp/14336.0_rp,DC6=C6-0.25_rp
	ndum=assert_eq(size(y),size(dydx),size(yout),size(yerr),'custom_rkck')
	ytemp=y+B21*h*dydx
	call derivs(x+A2*h,ytemp,ak2)
	ytemp=y+h*(B31*dydx+B32*ak2)
	call derivs(x+A3*h,ytemp,ak3)
	ytemp=y+h*(B41*dydx+B42*ak2+B43*ak3)
	call derivs(x+A4*h,ytemp,ak4)
	ytemp=y+h*(B51*dydx+B52*ak2+B53*ak3+B54*ak4)
	call derivs(x+A5*h,ytemp,ak5)
	ytemp=y+h*(B61*dydx+B62*ak2+B63*ak3+B64*ak4+B65*ak5)
	call derivs(x+A6*h,ytemp,ak6)
	yout=y+h*(C1*dydx+C3*ak3+C4*ak4+C6*ak6)
	yerr=h*(DC1*dydx+DC3*ak3+DC4*ak4+DC5*ak5+DC6*ak6)
	END SUBROUTINE custom_rkck

end module