dstode Subroutine

subroutine dstode(Neq,Y,Yh,Nyh,Yh1,Ewt,Savf,Acor,Wm,Iwm,f,jac,pjac,slvs)
use M_odepack
implicit none
integer,parameter :: dp=kind(0.0d0)

integer                     :: Neq(*)
real(kind=dp),intent(inout) :: Y(*)
integer                     :: Nyh
real(kind=dp),intent(inout) :: Yh(Nyh,*)
real(kind=dp),intent(inout) :: Yh1(*)
real(kind=dp)               :: Ewt(*)
real(kind=dp),intent(inout) :: Savf(*)
real(kind=dp),intent(inout) :: Acor(*)
real(kind=dp)               :: Wm(*)
integer                     :: Iwm(*)
external f
external jac
external pjac
external slvs

real(kind=dp) :: dcon, ddn, del, delp, dsm, dup, exdn, exsm, exup, r, rh, rhdn, rhsm, rhup, told
integer       :: i, i1, iredo, iret, j, jb, m, ncf, newq

dls1%kflag = 0
told = dls1%tn
ncf = 0
dls1%ierpj = 0
dls1%iersl = 0
dls1%jcur = 0
dls1%icf = 0
delp = 0.0D0
if ( dls1%jstart>0 ) goto 400
if ( dls1%jstart==-1 ) then
!-----------------------------------------------------------------------
!  The following block handles preliminaries needed when dls1%JSTART = -1.
!  IPUP is set to dls1%MITER to force a matrix update.
!  If an order increase is about to be considered (IALTH = 1),
!  IALTH is reset to 2 to postpone consideration one more step.
!  If the caller has changed dls1%METH, DCFODE is called to reset
!  the coefficients of the method.
!  If the caller has changed dls1%MAXORD to a value less than the current
!  order NQ, NQ is reduced to dls1%MAXORD, and a new dls1%H chosen accordingly.
!  If dls1%H is to be changed, YH must be rescaled.
!  If dls1%H or dls1%METH is being changed, IALTH is reset to L = NQ + 1
!  to prevent further changes in dls1%H for that many steps.
!-----------------------------------------------------------------------
   dls1%ipup = dls1%miter
   dls1%lmax = dls1%maxord + 1
   if ( dls1%ialth==1 ) dls1%ialth = 2
   if ( dls1%meth/=dls1%meo ) then
      call dcfode(dls1%meth,dls1%elco,dls1%tesco)
      dls1%meo = dls1%meth
      if ( dls1%nq<=dls1%maxord ) then
         dls1%ialth = dls1%l
         iret = 1
         goto 100
      endif
   elseif ( dls1%nq<=dls1%maxord ) then
      goto 200
   endif
   dls1%nq = dls1%maxord
   dls1%l = dls1%lmax
   do i = 1, dls1%l
      dls1%el(i) = dls1%elco(i,dls1%nq)
   enddo
   dls1%nqnyh = dls1%nq*Nyh
   dls1%rc = dls1%rc*dls1%el(1)/dls1%el0
   dls1%el0 = dls1%el(1)
   dls1%conit = 0.5D0/(dls1%nq+2)
   ddn = dvnorm(dls1%n,Savf,Ewt)/dls1%tesco(1,dls1%l)
   exdn = 1.0D0/dls1%l
   rhdn = 1.0D0/(1.3D0*ddn**exdn+0.0000013D0)
   rh = min(rhdn,1.0D0)
   iredo = 3
   if ( dls1%h==dls1%hold ) then
      rh = max(rh,dls1%hmin/abs(dls1%h))
   else
      rh = min(rh,abs(dls1%h/dls1%hold))
      dls1%h = dls1%hold
   endif
   goto 300
else
   if ( dls1%jstart==-2 ) goto 200
!-----------------------------------------------------------------------
!  On the first call, the order is set to 1, and other variables are
!  initialized.  RMAX is the maximum ratio by which dls1%H can be increased
!  in a single step.  It is initially 1.E4 to compensate for the small
!  initial dls1%H, but then is normally equal to 10.  If a failure
!  occurs (in corrector convergence or error test), RMAX is set to 2
!  for the next increase.
!-----------------------------------------------------------------------
   dls1%lmax = dls1%maxord + 1
   dls1%nq = 1
   dls1%l = 2
   dls1%ialth = 2
   dls1%rmax = 10000.0D0
   dls1%rc = 0.0D0
   dls1%el0 = 1.0D0
   dls1%crate = 0.7D0
   dls1%hold = dls1%h
   dls1%meo = dls1%meth
   dls1%nslp = 0
   dls1%ipup = dls1%miter
   iret = 3
!-----------------------------------------------------------------------
!  DCFODE is called to get all the integration coefficients for the
!  current dls1%METH.  Then the EL vector and related constants are reset
!  whenever the order NQ is changed, or at the start of the problem.
!-----------------------------------------------------------------------
   call dcfode(dls1%meth,dls1%elco,dls1%tesco)
endif
 100  continue
do i = 1, dls1%l
   dls1%el(i) = dls1%elco(i,dls1%nq)
enddo
dls1%nqnyh = dls1%nq*Nyh
dls1%rc = dls1%rc*dls1%el(1)/dls1%el0
dls1%el0 = dls1%el(1)
dls1%conit = 0.5D0/(dls1%nq+2)
select case (iret)
case (2)
   rh = max(rh,dls1%hmin/abs(dls1%h))
   goto 300
case (3)
   goto 400
case default
endselect
!-----------------------------------------------------------------------
!  If dls1%H is being changed, the dls1%H ratio RH is checked against
!  RMAX, dls1%HMIN, and dls1%HMXI, and the YH array rescaled.  IALTH is set to
!  L = NQ + 1 to prevent a change of dls1%H for that many steps, unless
!  forced by a convergence or error test failure.
!-----------------------------------------------------------------------
 200  continue
if ( dls1%h==dls1%hold ) goto 400
rh = dls1%h/dls1%hold
dls1%h = dls1%hold
iredo = 3
 300  continue
rh = min(rh,dls1%rmax)
rh = rh/max(1.0D0,abs(dls1%h)*dls1%hmxi*rh)
r = 1.0D0
do j = 2, dls1%l
   r = r*rh
   do i = 1, dls1%n
      Yh(i,j) = Yh(i,j)*r
   enddo
enddo
dls1%h = dls1%h*rh
dls1%rc = dls1%rc*rh
dls1%ialth = dls1%l
if ( iredo==0 ) then
   dls1%rmax = 10.0D0
   goto 1200
endif
!-----------------------------------------------------------------------
!  This section computes the predicted values by effectively
!  multiplying the YH array by the Pascal Triangle matrix.
!  RC is the ratio of new to old values of the coefficient  dls1%H*EL(1).
!  When RC differs from 1 by more than dls1%CCMAX, IPUP is set to dls1%MITER
!  to force PJAC to be called, if a Jacobian is involved.
!  In any case, PJAC is called at least every dls1%MSBP steps.
!-----------------------------------------------------------------------
 400  continue
if ( abs(dls1%rc-1.0D0)>dls1%ccmax ) dls1%ipup = dls1%miter
if ( dls1%nst>=dls1%nslp+dls1%msbp ) dls1%ipup = dls1%miter
dls1%tn = dls1%tn + dls1%h
i1 = dls1%nqnyh + 1
do jb = 1, dls1%nq
   i1 = i1 - Nyh
! dir$ ivdep
   do i = i1, dls1%nqnyh
      Yh1(i) = Yh1(i) + Yh1(i+Nyh)
   enddo
enddo
!-----------------------------------------------------------------------
!  Up to dls1%MAXCOR corrector iterations are taken.  A convergence test is
!  made on the R.M.S. norm of each correction, weighted by the error
!  weight vector EWT.  The sum of the corrections is accumulated in the
!  vector ACOR(i).  The YH array is not altered in the corrector loop.
!-----------------------------------------------------------------------
 500  continue
m = 0
do i = 1, dls1%n
   Y(i) = Yh(i,1)
enddo
call f(Neq,dls1%tn,Y,Savf)
dls1%nfe = dls1%nfe + 1
if ( dls1%ipup>0 ) then
!-----------------------------------------------------------------------
!  If indicated, the matrix P = I - dls1%h*dls1%el(1)*J is reevaluated and
!  preprocessed before starting the corrector iteration.  IPUP is set
!  to 0 as an indicator that this has been done.
!-----------------------------------------------------------------------
   call pjac(Neq,Y,Yh,Nyh,Ewt,Acor,Savf,Wm,Iwm,f,jac)
   dls1%ipup = 0
   dls1%rc = 1.0D0
   dls1%nslp = dls1%nst
   dls1%crate = 0.7D0
   if ( dls1%ierpj/=0 ) goto 800
endif
do i = 1, dls1%n
   Acor(i) = 0.0D0
enddo
 600  continue
if ( dls1%miter/=0 ) then
!-----------------------------------------------------------------------
!  In the case of the chord method, compute the corrector error,
!  and solve the linear system with that as right-hand side and
!  P as coefficient matrix.
!-----------------------------------------------------------------------
   do i = 1, dls1%n
      Y(i) = dls1%h*Savf(i) - (Yh(i,2)+Acor(i))
   enddo
   call slvs(Wm,Iwm,Y,Savf)
   if ( dls1%iersl<0 ) goto 800
   if ( dls1%iersl>0 ) goto 700
   del = dvnorm(dls1%n,Y,Ewt)
   do i = 1, dls1%n
      Acor(i) = Acor(i) + Y(i)
      Y(i) = Yh(i,1) + dls1%el(1)*Acor(i)
   enddo
else
!-----------------------------------------------------------------------
!  In the case of functional iteration, update Y directly from
!  the result of the last function evaluation.
!-----------------------------------------------------------------------
   do i = 1, dls1%n
      Savf(i) = dls1%h*Savf(i) - Yh(i,2)
      Y(i) = Savf(i) - Acor(i)
   enddo
   del = dvnorm(dls1%n,Y,Ewt)
   do i = 1, dls1%n
      Y(i) = Yh(i,1) + dls1%el(1)*Savf(i)
      Acor(i) = Savf(i)
   enddo
endif
!-----------------------------------------------------------------------
!  Test for convergence.  If M.gt.0, an estimate of the convergence
!  rate constant is stored in CRATE, and this is used in the test.
!-----------------------------------------------------------------------
if ( m/=0 ) dls1%crate = max(0.2D0*dls1%crate,del/delp)
dcon = del*min(1.0D0,1.5D0*dls1%crate)/(dls1%tesco(2,dls1%nq)*dls1%conit)
if ( dcon<=1.0D0 ) then
!-----------------------------------------------------------------------
!  The corrector has converged.  JCUR is set to 0
!  to signal that the Jacobian involved may need updating later.
!  The local error test is made and control passes to statement 500
!  if it fails.
!-----------------------------------------------------------------------
   dls1%jcur = 0
   if ( m==0 ) dsm = del/dls1%tesco(2,dls1%nq)
   if ( m>0 ) dsm = dvnorm(dls1%n,Acor,Ewt)/dls1%tesco(2,dls1%nq)
   if ( dsm>1.0D0 ) then
!-----------------------------------------------------------------------
!  The error test failed.  dls1%KFLAG keeps track of multiple failures.
!  Restore dls1%TN and the YH array to their previous values, and prepare
!  to try the step again.  Compute the optimum step size for this or
!  one lower order.  After 2 or more failures, dls1%H is forced to decrease
!  by a factor of 0.2 or less.
!-----------------------------------------------------------------------
      dls1%kflag = dls1%kflag - 1
      dls1%tn = told
      i1 = dls1%nqnyh + 1
      do jb = 1, dls1%nq
         i1 = i1 - Nyh
! dir$ ivdep
         do i = i1, dls1%nqnyh
            Yh1(i) = Yh1(i) - Yh1(i+Nyh)
         enddo
      enddo
      dls1%rmax = 2.0D0
      if ( abs(dls1%h)<=dls1%hmin*1.00001D0 ) then
!-----------------------------------------------------------------------
!  All returns are made through this section.  dls1%H is saved in HOLD
!  to allow the caller to change dls1%H on the next step.
!-----------------------------------------------------------------------
         dls1%kflag = -1
         goto 1300
      elseif ( dls1%kflag<=-3 ) then
!-----------------------------------------------------------------------
!  Control reaches this section if 3 or more failures have occured.
!  If 10 failures have occurred, exit with dls1%KFLAG = -1.
!  It is assumed that the derivatives that have accumulated in the
!  YH array have errors of the wrong order.  Hence the first
!  derivative is recomputed, and the order is set to 1.  Then
!  dls1%H is reduced by a factor of 10, and the step is retried,
!  until it succeeds or dls1%H reaches dls1%HMIN.
!-----------------------------------------------------------------------
         if ( dls1%kflag==-10 ) then
            dls1%kflag = -1
            goto 1300
         else
            rh = 0.1D0
            rh = max(dls1%hmin/abs(dls1%h),rh)
            dls1%h = dls1%h*rh
            do i = 1, dls1%n
               Y(i) = Yh(i,1)
            enddo
            call f(Neq,dls1%tn,Y,Savf)
            dls1%nfe = dls1%nfe + 1
            do i = 1, dls1%n
               Yh(i,2) = dls1%h*Savf(i)
            enddo
            dls1%ipup = dls1%miter
            dls1%ialth = 5
            if ( dls1%nq==1 ) goto 400
            dls1%nq = 1
            dls1%l = 2
            iret = 3
            goto 100
         endif
      else
         iredo = 2
         rhup = 0.0D0
         goto 900
      endif
   else
!-----------------------------------------------------------------------
!  After a successful step, update the YH array.
!  Consider changing dls1%H if IALTH = 1.  Otherwise decrease IALTH by 1.
!  If IALTH is then 1 and NQ .lt. dls1%MAXORD, then ACOR is saved for
!  use in a possible order increase on the next step.
!  If a change in dls1%H is considered, an increase or decrease in order
!  by one is considered also.  A change in dls1%H is made only if it is by a
!  factor of at least 1.1.  If not, IALTH is set to 3 to prevent
!  testing for that many steps.
!-----------------------------------------------------------------------
      dls1%kflag = 0
      iredo = 0
      dls1%nst = dls1%nst + 1
      dls1%hu = dls1%h
      dls1%nqu = dls1%nq
      do j = 1, dls1%l
         do i = 1, dls1%n
            Yh(i,j) = Yh(i,j) + dls1%el(j)*Acor(i)
         enddo
      enddo
      dls1%ialth = dls1%ialth - 1
      if ( dls1%ialth==0 ) then
!-----------------------------------------------------------------------
!  Regardless of the success or failure of the step, factors
!  RHDN, RHSM, and RHUP are computed, by which dls1%H could be multiplied
!  at order NQ - 1, order NQ, or order NQ + 1, respectively.
!  In the case of failure, RHUP = 0.0 to avoid an order increase.
!  The largest of these is determined and the new order chosen
!  accordingly.  If the order is to be increased, we compute one
!  additional scaled derivative.
!-----------------------------------------------------------------------
         rhup = 0.0D0
         if ( dls1%l/=dls1%lmax ) then
            do i = 1, dls1%n
               Savf(i) = Acor(i) - Yh(i,dls1%lmax)
            enddo
            dup = dvnorm(dls1%n,Savf,Ewt)/dls1%tesco(3,dls1%nq)
            exup = 1.0D0/(dls1%l+1)
            rhup = 1.0D0/(1.4D0*dup**exup+0.0000014D0)
         endif
         goto 900
      else
         if ( dls1%ialth<=1 ) then
            if ( dls1%l/=dls1%lmax ) then
               do i = 1, dls1%n
                  Yh(i,dls1%lmax) = Acor(i)
               enddo
            endif
         endif
         goto 1200
      endif
   endif
else
   m = m + 1
   if ( m/=dls1%maxcor ) then
      if ( m<2 .or. del<=2.0D0*delp ) then
         delp = del
         call f(Neq,dls1%tn,Y,Savf)
         dls1%nfe = dls1%nfe + 1
         goto 600
      endif
   endif
endif
!-----------------------------------------------------------------------
!  The corrector iteration failed to converge.
!  If dls1%MITER .ne. 0 and the Jacobian is out of date, PJAC is called for
!  the next try.  Otherwise the YH array is retracted to its values
!  before prediction, and dls1%H is reduced, if possible.  If dls1%H cannot be
!  reduced or dls1%MXNCF failures have occurred, exit with dls1%KFLAG = -2.
!-----------------------------------------------------------------------
 700  continue
if ( dls1%miter/=0 .and. dls1%jcur/=1 ) then
   dls1%icf = 1
   dls1%ipup = dls1%miter
   goto 500
endif
 800  continue
dls1%icf = 2
ncf = ncf + 1
dls1%rmax = 2.0D0
dls1%tn = told
i1 = dls1%nqnyh + 1
do jb = 1, dls1%nq
   i1 = i1 - Nyh
! dir$ ivdep
   do i = i1, dls1%nqnyh
      Yh1(i) = Yh1(i) - Yh1(i+Nyh)
   enddo
enddo
if ( dls1%ierpj<0 .or. dls1%iersl<0 ) then
   dls1%kflag = -3
   goto 1300
elseif ( abs(dls1%h)<=dls1%hmin*1.00001D0 ) then
   dls1%kflag = -2
   goto 1300
elseif ( ncf==dls1%mxncf ) then
   dls1%kflag = -2
   goto 1300
else
   rh = 0.25D0
   dls1%ipup = dls1%miter
   iredo = 1
   rh = max(rh,dls1%hmin/abs(dls1%h))
   goto 300
endif
 900  continue
exsm = 1.0D0/dls1%l
rhsm = 1.0D0/(1.2D0*dsm**exsm+0.0000012D0)
rhdn = 0.0D0
if ( dls1%nq/=1 ) then
   ddn = dvnorm(dls1%n,Yh(1,dls1%l),Ewt)/dls1%tesco(1,dls1%nq)
   exdn = 1.0D0/dls1%nq
   rhdn = 1.0D0/(1.3D0*ddn**exdn+0.0000013D0)
endif
if ( rhsm>=rhup ) then
   if ( rhsm>=rhdn ) then
      newq = dls1%nq
      rh = rhsm
      goto 1000
   endif
elseif ( rhup>rhdn ) then
   newq = dls1%l
   rh = rhup
   if ( rh<1.1D0 ) then
      dls1%ialth = 3
      goto 1200
   else
      r = dls1%el(dls1%l)/dls1%l
      do i = 1, dls1%n
         Yh(i,newq+1) = Acor(i)*r
      enddo
      goto 1100
   endif
endif
newq = dls1%nq - 1
rh = rhdn
if ( dls1%kflag<0 .and. rh>1.0D0 ) rh = 1.0D0
 1000 continue
if ( (dls1%kflag==0) .and. (rh<1.1D0) ) then
   dls1%ialth = 3
   goto 1200
else
   if ( dls1%kflag<=-2 ) rh = min(rh,0.2D0)
!-----------------------------------------------------------------------
!  If there is a change of order, reset NQ, dls1%l, and the coefficients.
!  In any case dls1%H is reset according to RH and the YH array is rescaled.
!  Then exit from 690 if the step was OK, or redo the step otherwise.
!-----------------------------------------------------------------------
   if ( newq==dls1%nq ) then
      rh = max(rh,dls1%hmin/abs(dls1%h))
      goto 300
   endif
endif
 1100 continue
dls1%nq = newq
dls1%l = dls1%nq + 1
iret = 2
goto 100
 1200 continue
r = 1.0D0/dls1%tesco(2,dls1%nqu)
do i = 1, dls1%n
   Acor(i) = Acor(i)*r
enddo
 1300 continue
dls1%hold = dls1%h
dls1%jstart = 1
end subroutine dstode