mod_odeint.t

Go to the documentation of this file.

!> This module packages odeint from numerical recipes.
module mod_odeint
  use mod_comm_lib, only: mpistop
  implicit none
  double precision         :: tiny
  parameter(tiny=1.d-30)
  integer                  :: maxstp, nmax, kmaxx
  parameter(maxstp=10000,nmax=50,kmaxx=200)
 
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)
    double precision h,hdid,hnext,x,xsav,dydx(NMAX),y(NMAX),yscal(NMAX)
    double precision dxsav, yp(NMAX,KMAXX), xp(KMAXX)
    integer kmax, kount
    INTEGER i,nstp
    ! 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