      SUBROUTINE APFORM(N,NW,MU,OMEGA,LOLD,S,T,Y,H,TOLD,KOLD,TOL,W,C,
     +     FCN,SIGMA,YINIT,DYINIT)

C**********************************************************************
C                                                                     *
C     Purpose - This routine integrates one step by RK-methods        *
c               for delay DEs.                                        *
C                                                                     *
C**********************************************************************
C                                                                     *
C     N      Number of equations                                      *
C     NW     First dimension of workspace W & queue KOLD              *
C     MU     Number of solution delay terms                           *
C     OMEGA  Number of derivative delay terms                         *
C     LOLD   The length of queue                                      *
C     S      The number of stages for RK-formula                      *
C     T      Current independent value                                *
C     Y      Current solution values, i.e., at T                      *
C     H      Trial stepsize                                           *
C     TOLD   Queue of T values                                        *
C     KOLD   Queue oF K values and solution values                    *
C     TOL    Tolerance                                                *
C     W      Workspace                                                *
C     C      Communication vector                                     *
C            ...USE C(13) : Minimum stepsize                          *
C                   C(17) : T trial                                   *
C                   C(18) : H trial                                   *
C                   C(24) : # of function evaluations                 *
C                   C(25) : Initial point                             *
C                   C(26) : The first index of queue                  *
C                   C(27) : The last index of queue                   *
C                   C(28) : Flag for computation of K1 in the next    *
C                           step(if = 0.0, need to compute)           *
C                   C(29) : Flag indicating which of the              *
c                           extrapolation or the special interpolant  *
c                           must be used.                             *
C                           If C(29) = 0.0 then use extrapolation     *
C     FCN    Name of derivative routine                               *
C     SIGMA  Name of subroutine for computing delay argument          *
C     YINIT  Name of subroutine for computing initial values at       *
C            initial intervals.                                       *
C     DYINIT Name of subroutine for computing derivative of initial   *
C            values at initial intervals.                             *
C                                                                     *
C**********************************************************************

C**********************************************************************
C Note:                                                               *
C       When using extrapolation, we use global intrp routine.        *
C       And use told(c(27)-1). But c(27) may be one.                  *
C       In this case, we have to use told(lold - 1)                   *
C**********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,MU,OMEGA,LOLD,S
      DOUBLE PRECISION T,H,TOL

C     .. Array Arguments ..
      DOUBLE PRECISION  Y(N),TOLD(LOLD),KOLD(NW,7,LOLD),W(NW,*),C(*)

C     .. Subroutine Arguments ..
      EXTERNAL FCN,SIGMA,YINIT,DYINIT

C     .. Arrays in Common ..
      DOUBLE PRECISION AA(11,11),BB(11),CC(11)

C     .. Local Scalars ..
      INTEGER I,J,K,L,INDEX,IW,MAXNIT,ND,NIT,TMP
      DOUBLE PRECISION ERRKI,SIGN,TEMP,TEMP2

C     .. External Subroutines ..
      EXTERNAL EVALAGS,EVALAGD,INTRPG,DERIVG

C     .. Intrisic Functions ..
      INTRINSIC DMAX1,DFLOAT,IDINT

C     .. Common blocks ..
      COMMON/METH/AA,BB,CC

C     .. Data statements ..
      DATA MAXNIT/5/

C     .. Set up the nubmer of delay terms
      ND = MU + OMEGA

C     ... LOOP for Advanced Delay ...
 1000 CONTINUE

C     Set flag to be 1
      C(28) = 1.D0
C
C***** Beginning of the 1st phase *****

C     Set index to 0
C ... INDEX: The first vanishing delay occures in the INDEX-th
C     stage if not 0.
      INDEX = 0

C     The 2nd stage
C ... Compute Y values
      DO 110 I = 1, N
         W(I,16) = Y(I) + H*AA(2,1)*W(I,1)
 110  CONTINUE

C ... Compute argument for delay values
      TEMP = T + CC(2)*H
      CALL SIGMA(N,ND,TEMP,W(1,16),W(1,15))

C ... Compute solution delay values
      DO 120 I = 1, MU

         IW = 20 + I

         IF (W(I,15) .GT. TEMP) THEN
            PRINT *, '**** WARNING:APFORM *****'
            PRINT *, 'Advanced Delay at Stage  2'
            PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
            PRINT *

C ... Cut the stepseize to half and integrate again
            IF (W(I,15) - TEMP .GE. 1.D0) THEN
               PRINT *, 'ADVANCED TOO MUCH !!!'
               PRINT *
               C(18) = C(18) / 2.D0
               IF (C(18) .LT. C(13)) THEN
                  PRINT *, 'IND = -3,  At T = ', T
                  PRINT *
                  STOP
               END IF
               C(17) = T + C(18)
               H = C(18)
               GOTO 1000
            END IF
         END IF
C
         SIGN = W(I,15) - T
C
C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,YINIT)

C ... In case there is vanishing delay
         ELSE
            INDEX = 2

C ... Compute delay value using extrapolation
            IF (C(29) .EQ. 0.D0) THEN
	       TMP = IDINT(C(27))
               IF (TMP .EQ. 1) TMP = LOLD
               CALL INTRPG(N,NW,TMP,TOLD(TMP-1),W(I,15),W(1,IW),
     +              T-TOLD(TMP-1),KOLD)

C ... Compute delay value using the special interpolant
            ELSE
               DO 130 J = 1, N
                  W(J,IW) = Y(J) +  SIGN*W(J,1)
 130           CONTINUE
            END IF
         END IF

 120  CONTINUE

C ... Compute derivative delay values

      DO 125 I = MU + 1, ND

         IW = 20 + I

         IF (W(I,15) .GT. TEMP) THEN
            PRINT *, '**** WARNING:APFORM *****'
            PRINT *, 'Advanced Delay at Stage  2'
            PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
            PRINT *

C ... Cut the stepseize to half and integrate again
            IF (W(I,15) - TEMP .GE. 1.D0) THEN
               PRINT *, 'ADVANCED TOO MUCH !!!'
               PRINT *
               C(18) = C(18) / 2.D0
               IF (C(18) .LT. C(13)) THEN
                  PRINT *, 'IND = -3,  At T = ', T
                  PRINT *
                  STOP
               END IF
               C(17) = T + C(18)
               H = C(18)
               GOTO 1000
            END IF
         END IF
C
         SIGN = W(I,15) - T
C
C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,DYINIT)

C ... In case there is vanishing delay
         ELSE
            INDEX = 2

C ... Compute delay value using extrapolation
            IF (C(29) .EQ. 0.D0) THEN
	       TMP = IDINT(C(27))
               IF (TMP .EQ. 1) TMP = LOLD
               CALL DERIVG(N,NW,TMP,TOLD(TMP-1),W(I,15),W(1,IW),
     +              T-TOLD(TMP-1),KOLD)

C ... Compute delay value using the special interpolant
            ELSE
               DO 135 J = 1, N
                  W(J,IW) = W(J,1)
 135           CONTINUE
            END IF
         END IF

 125  CONTINUE

C ... Compute K2

      CALL FCN(N,ND,NW,TEMP,W(1,16),W(1,21),W(1,2))

C ... Compute K3 to KS
      DO 140 K = 3, S

         TEMP = T + CC(K)*H

C ... Compute Y values
         DO 150 I = 1, N
            W(I,16) = AA(K,1)*W(I,1)
            DO 160 J = 2, K - 1
               W(I,16) = W(I,16) + AA(K,J)*W(I,J)
 160        CONTINUE
            W(I,16) = Y(I) + H*W(I,16)
 150     CONTINUE

C ... For consistency with asymptotically valid estimate
         if (k. eq. 8) then
            do 165 i = 1, n
               w(i,12) = w(i,16)
 165        continue
         end if

C ... Compute arguments for delay values
         CALL SIGMA(N,ND,TEMP,W(1,16),W(1,15))

C ... Compute solution delay values
         DO 170 I = 1, MU

            IW = 20 + I

            IF (W(I,15) .GT. TEMP) THEN
               PRINT *, '**** WARNING:APFORM *****'
               PRINT *, 'AdvanceD Delay at Stage', K
               PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
               PRINT *

C ... Cut the stepseize to half and integrate again
               IF (W(I,15)-TEMP .GE. 1.D0) THEN
                  PRINT *, 'ADVANCED TOO MUCH !!!'
                  PRINT *
                  C(18) = C(18) / 2.D0
                  IF (C(18) .LT. C(13)) THEN
                     PRINT *, 'IND = -3,  At T = ', T
                     STOP
                  END IF
                  C(17) = T + C(18)
                  C(24) = C(24) + DFLOAT(K - 2)
                  H = C(18)
                  GOTO 1000
               END IF
            END IF

            SIGN = W(I,15) - T
            TEMP2 = SIGN**2 / (CC(K)*H)

C ... In case there is no vanishing delay
            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +              TOLD,KOLD,W,YINIT)

C ... In case there is vanishing delay
            ELSE
               IF (INDEX .EQ. 0) INDEX = K

C ... compute delay values using extrapolation
               IF (C(29) .EQ. 0.D0) THEN
	          TMP = IDINT(C(27))
                  IF (TMP .EQ. 1) TMP = LOLD
                  CALL INTRPG(N,NW,TMP,TOLD(TMP-1),W(I,15),W(1,IW),
     +                 T-TOLD(TMP-1),KOLD)

C ... Compute delay value using the special interpolant
               ELSE
                  DO 180 L = 1, N
                     W(L,IW) = AA(K,1) * W(L,1)
                     DO 190 J = 2, K-1
                        W(L,IW) = W(L,IW) + AA(K,J) * W(L,J)
 190                 CONTINUE
                     W(L,IW) = Y(L) + (SIGN - TEMP2) * W(L,1)
     +                    + TEMP2 * W(L,IW) / CC(K)

 180              CONTINUE
               END IF
            END IF
 170     CONTINUE

C ... Compute derivative delay terms
         DO 175 I = MU + 1, ND

            IW = 20 + I

            IF (W(I,15) .GT. TEMP) THEN
               PRINT *, '**** WARNING:APFORM *****'
               PRINT *, 'AdvanceD Delay at Stage', K
               PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
               PRINT *

C ... Cut the stepseize to half and integrate again
               IF (W(I,15)-TEMP .GE. 1.D0) THEN
                  PRINT *, 'ADVANCED TOO MUCH !!!'
                  PRINT *
                  C(18) = C(18) / 2.D0
                  IF (C(18) .LT. C(13)) THEN
                     PRINT *, 'IND = -3,  At T = ', T
                     STOP
                  END IF
                  C(17) = T + C(18)
                  C(24) = C(24) + DFLOAT(K - 2)
                  H = C(18)
                  GOTO 1000
               END IF
            END IF

            SIGN = W(I,15) - T

C ... In case there is no vanishing delay
            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +              TOLD,KOLD,W,DYINIT)

C ... In case there is vanishing delay
            ELSE
               IF (INDEX .EQ. 0) INDEX = K

C ... compute delay values using extrapolation
               IF (C(29) .EQ. 0.D0) THEN
	          TMP = IDINT(C(27))
                  IF (TMP .EQ. 1) TMP = LOLD
                  CALL DERIVG(N,NW,TMP,TOLD(TMP-1),W(I,15),W(1,IW),
     +                 T-TOLD(TMP-1),KOLD)

C ... Compute delay value using the special interpolant
               ELSE

                  DO 185 L = 1, N
                     W(L,IW) = AA(K,1) * W(L,1)
                     DO 195 J = 2, K-1
                        W(L,IW) = W(L,IW) + AA(K,J) * W(L,J)
 195                 CONTINUE
                     W(L,IW) = W(L,1) - (W(L,1) - W(L,IW) / CC(K)) *
     +                    (2.D0 * SIGN / (H * CC(K)))
 185              CONTINUE
               END IF
            END IF
 175        CONTINUE

C ... Compute K-value

         CALL FCN(N,ND,NW,TEMP,W(1,16),W(1,21),W(1,K))

 140  CONTINUE

C***** End of the 1st phase *****

C***** Beginning of the 2nd phase *****

C ... Set # of iteration to be 0
      NIT = 0

C ... If index is not 0, then there is a vanishg delay
      IF (INDEX .NE. 0) THEN

C ... Loop for iterative scheme ...
 2000    CONTINUE

         NIT = NIT + 1
         ERRKI = 0.D0

C ... Recompute K-values for stage INDEX+1 to s
         DO 210 K = INDEX, S

            DO 220 I = 1, N
               W(I,13) = W(I,K)
 220        CONTINUE

            TEMP = T + CC(K)*H

C ... Recompute y values
            DO 230 I = 1, N
               W(I,16) = AA(K,1)*W(I,1)
               DO 240 J = 2, K-1
                  W(I,16) = W(I,16) + AA(K,J)*W(I,J)
 240           CONTINUE
               W(I,16) = Y(I) + H*W(I,16)
 230        CONTINUE

C ... For consistency with Asymptotically valid estimate
            if (k .eq. 8) then
               do 235 i = 1, n
                  w(i,12) = w(i,16)
 235           continue
            end if

C ... Recompute arguments for delay values
            CALL SIGMA(N,ND,TEMP,W(1,16),W(1,15))

C ... Recompute solution delay values using intrp routine
            DO 250 I = 1, MU

               IW = 20 + I

               IF (W(I,15) .GT. TEMP) THEN
                  PRINT *, '**** WARNING:APFORM *****'
                  PRINT *, '    Advanced Delay at Stage ', K
                  PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
                  PRINT *

C ... Cut the stepseize to half and integrate again
                  IF (W(I,15)-TEMP .GE. 1.D0) THEN
                     PRINT *, 'ADVANCED TOO MUCH !!!'
                     PRINT *
                     C(18) = C(18) / 2.D0
                     IF (C(18) .LT. C(13)) THEN
                        PRINT *, 'IND = -3,  At T = ', T
                        STOP
                     END IF
                     C(17) = T + C(18)
                     C(24) = C(24) + DFLOAT(S-1 + K-INDEX)
                     H = C(18)
                     GOTO 1000
                  END IF
               END IF

               SIGN = W(I,15) - T

C ... In case there is no vanisning delay
               IF (SIGN .LE. 0.D0) THEN
                  CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +                 TOLD,KOLD,W,YINIT)
               ELSE
                  CALL INTRP(N,NW,T,Y,W(I,15),W(1,IW),H,W)
               END IF
 250        CONTINUE

C ... Recompute derivative delay values using intrp routine
            DO 255 I = MU + 1, ND

               IW = 20 + I

               IF (W(I,15) .GT. TEMP) THEN
                  PRINT *, '**** WARNING:APFORM *****'
                  PRINT *, '    Advanced Delay at Stage ', K
                  PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
                  PRINT *

C ... Cut the stepseize to half and integrate again
                  IF (W(I,15)-TEMP .GE. 1.D0) THEN
                     PRINT *, 'ADVANCED TOO MUCH !!!'
                     PRINT *
                     C(18) = C(18) / 2.D0
                     IF (C(18) .LT. C(13)) THEN
                        PRINT *, 'IND = -3,  At T = ', T
                        STOP
                     END IF
                     C(17) = T + C(18)
                     C(24) = C(24) + DFLOAT(S-1 + K-INDEX)
                     H = C(18)
                     GOTO 1000
                  END IF
               END IF

               SIGN = W(I,15) - T

C ... In case there is no vanisning delay
               IF (SIGN .LE. 0.D0) THEN
                  CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +                 TOLD,KOLD,W,DYINIT)
               ELSE

                  CALL DERIV(N,NW,T,Y,W(I,15),W(1,IW),H,W)
               END IF
 255        CONTINUE

C ... Recompute K-value

            CALL FCN(N,ND,NW,TEMP,W(1,16),W(1,21),W(1,K))

            DO 260 I = 1, N
               ERRKI = DMAX1(ERRKI, DABS(W(I,13) - W(I,K)))
 260        CONTINUE

 210     CONTINUE
         C(24) = C(24) + DFLOAT(S - INDEX + 1)

         IF ((ERRKI .GT. TOL) .AND. (NIT .LT. MAXNIT)) GOTO 2000

      END IF

c      PRINT *, 'INDEX =', INDEX
c      PRINT *, 'NIT =', NIT
c      PRINT *, 'ERRKI =', ERRKI
c      PRINT *

C***** End of the the 2nd phase *****

C***** Beginning of the 3rd phase *****

C ... Improve derivative values on the mesh point for the interpolation
C ... to archieve C1 continuity
      IF ((INDEX .GT. 0) .AND. (INDEX .LE. 8)) THEN
         TEMP = T + CC(8)*H

C ... Compute Y values
         DO 310 I = 1, N
            W(I,12) = AA(8,1)*W(I,1)
            DO 320 J = 2, 7
               W(I,12) = W(I,12) + AA(8,J)*W(I,J)
 320        CONTINUE
            W(I,12) = Y(I) + H*W(I,12)
 310     CONTINUE

C ... Compute arguments for delay values
         CALL SIGMA(N,ND,TEMP,W(1,12),W(1,15))

C ... Compute solution delay values
         DO 330 I = 1, MU

            IW = 20 + I
            SIGN = W(I,15) - T
            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +                      TOLD,KOLD,W,YINIT)
            ELSE
               CALL INTRP(N,NW,T,Y,W(I,15),W(1,IW),H,W)
            END IF

 330     CONTINUE

C ... Compute derivative delay values
         DO 335 I = MU + 1, ND

            IW = 20 + I
            SIGN = W(I,15) - T
            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +                      TOLD,KOLD,W,DYINIT)
            ELSE

               CALL DERIV(N,NW,T,Y,W(I,15),W(1,IW),H,W)
            END IF
 335     CONTINUE

C ... Recompute K-value
         CALL FCN(N,ND,NW,TEMP,W(1,12),W(1,21),W(1,8))

         C(24) = C(24) + 1.D0
         C(28) = 0.D0
      END IF

C***** End of the 3nd phase *****

C ... Finally compute solution at T + H in W(*,12)
      DO 410 I = 1, N
         W(I,12) = BB(1)*W(I,1)
         DO 420 J = 2, S
            W(I,12) = W(I,12) + BB(J)*W(I,J)
 420     CONTINUE
         W(I,12) = Y(I) + H*W(I,12)
 410  CONTINUE

c      print *, 'index =', index

      RETURN
      END

C *** End of APFORM ***

      SUBROUTINE APFORM_ASM(N,NW,MU,OMEGA,LOLD,S,T,Y,H,TOLD,KOLD,TOL,W,C
     +     ,FCN,SIGMA,YINIT,DYINIT)

C**********************************************************************
C                                                                     *
C     Purpose - This routine integrates one step by RK-methods        *
c               for delay DEs.                                        *
C                                                                     *
C**********************************************************************
C                                                                     *
C     N      Number of equations                                      *
C     NW     First dimension of workspace W & queue KOLD              *
C     MU     Number of solution delay terms                           *
C     OMEGA  Number of derivative delay terms                         *
C     LOLD   The length of queue                                      *
C     S      The number of stages for RK-formula                      *
C     T      Current independent value                                *
C     Y      Current solution values, i.e., at T                      *
C     H      Trial stepsize                                           *
C     TOLD   Queue of T values                                        *
C     KOLD   Queue oF K values and solution values                    *
C     TOL    Tolerance                                                *
C     W      Workspace                                                *
C     C      Communication vector                                     *
C            ...USE C(13) : Minimum stepsize                          *
C                   C(17) : T trial                                   *
C                   C(18) : H trial                                   *
C                   C(24) : # of function evaluations                 *
C                   C(25) : Initial point                             *
C                   C(26) : The first index of queue                  *
C                   C(27) : The last index of queue                   *
C                   C(28) : Flag for computation of K1 in the next    *
C                           step(if = 0.0, need to compute)           *
C                   C(29) : Flag indicating which of the              *
c                           extrapolation or the special interpolant  *
c                           must be used.                             *
C                           If C(29) = 0.0 then use extrapolation     *
C     FCN    Name of derivative routine                               *
C     SIGMA  Name of subroutine for computing delay argument          *
C     YINIT  Name of subroutine for computing initial values at       *
C            initial intervals.                                       *
C     DYINIT Name of subroutine for computing derivative of initial   *
C            values at initial intervals.                             *
C                                                                     *
C**********************************************************************

C**********************************************************************
C Note:                                                               *
C       When using extrapolation, we use global intrp routine.        *
C       And use told(c(27)-1). But c(27) may be one.                  *
C       In this case, we have to use told(lold - 1)                   *
C**********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,MU,OMEGA,LOLD,S
      DOUBLE PRECISION T,H,TOL

C     .. Array Arguments ..
      DOUBLE PRECISION  Y(N),TOLD(LOLD),KOLD(NW,7,LOLD),W(NW,*),C(*)

C     .. Subroutine Arguments ..
      EXTERNAL FCN,SIGMA,YINIT,DYINIT

C     .. Arrays in Common ..
      DOUBLE PRECISION AA(11,11),BB(11),CC(11)

C     .. Local Scalars ..
      INTEGER I,J,K,L,INDEX,IW,MAXNIT,ND,NIT,TMP
      DOUBLE PRECISION ERRKI,SIGN,TEMP,TEMP2

C     .. External Subroutines ..
      EXTERNAL EVALAGS,EVALAGD,INTRPG,DERIVG

C     .. Intrisic Functions ..
      INTRINSIC DMAX1,DFLOAT,IDINT

C     .. Common blocks ..
      COMMON/METH/AA,BB,CC

C     .. Data statements ..
      DATA MAXNIT/11/
c      DATA MAXNIT/5/

C     .. Set up the nubmer of delay terms
      ND = MU + OMEGA

C     ... LOOP for Advanced Delay ...
 1000 CONTINUE

C     Set flag to be 1
      C(28) = 1.D0
C
C***** Beginning of the 1st phase *****

C     Set index to 0
C ... INDEX: The first vanishing delay occures in the INDEX-th
C     stage if not 0.
      INDEX = 0

C     The 2nd stage
C ... Compute Y values
      DO 110 I = 1, N
         W(I,16) = Y(I) + H*AA(2,1)*W(I,1)
 110  CONTINUE

C ... Compute argument for delay values
      TEMP = T + CC(2)*H
      CALL SIGMA(N,ND,TEMP,W(1,16),W(1,15))

C ... Compute solution delay values
      DO 120 I = 1, MU

         IW = 20 + I

         IF (W(I,15) .GT. TEMP) THEN
            PRINT *, '**** WARNING:APFORM *****'
            PRINT *, 'Advanced Delay at Stage  2'
            PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
            PRINT *

C ... Cut the stepseize to half and integrate again
            IF (W(I,15) - TEMP .GE. 1.D0) THEN
               PRINT *, 'ADVANCED TOO MUCH !!!'
               PRINT *
               C(18) = C(18) / 2.D0
               IF (C(18) .LT. C(13)) THEN
                  PRINT *, 'IND = -3,  At T = ', T
                  PRINT *
                  STOP
               END IF
               C(17) = T + C(18)
               H = C(18)
               GOTO 1000
            END IF
         END IF
C
         SIGN = W(I,15) - T
C
C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,YINIT)

C ... In case there is vanishing delay
         ELSE
            INDEX = 2

C ... Compute delay value using extrapolation
            IF (C(29) .EQ. 0.D0) THEN
	       TMP = IDINT(C(27))
               IF (TMP .EQ. 1) TMP = LOLD
               CALL INTRPG(N,NW,TMP,TOLD(TMP-1),W(I,15),W(1,IW),
     +              T-TOLD(TMP-1),KOLD)

C ... Compute delay value using the special interpolant
            ELSE
               DO 130 J = 1, N
                  W(J,IW) = Y(J) +  SIGN*W(J,1)
 130           CONTINUE
            END IF
         END IF

 120  CONTINUE

C ... Compute derivative delay values

      DO 125 I = MU + 1, ND

         IW = 20 + I

         IF (W(I,15) .GT. TEMP) THEN
            PRINT *, '**** WARNING:APFORM *****'
            PRINT *, 'Advanced Delay at Stage  2'
            PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
            PRINT *

C ... Cut the stepseize to half and integrate again
            IF (W(I,15) - TEMP .GE. 1.D0) THEN
               PRINT *, 'ADVANCED TOO MUCH !!!'
               PRINT *
               C(18) = C(18) / 2.D0
               IF (C(18) .LT. C(13)) THEN
                  PRINT *, 'IND = -3,  At T = ', T
                  PRINT *
                  STOP
               END IF
               C(17) = T + C(18)
               H = C(18)
               GOTO 1000
            END IF
         END IF
C
         SIGN = W(I,15) - T
C
C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,DYINIT)

C ... In case there is vanishing delay
         ELSE
            INDEX = 2

C ... Compute delay value using extrapolation
            IF (C(29) .EQ. 0.D0) THEN
	       TMP = IDINT(C(27))
               IF (TMP .EQ. 1) TMP = LOLD
               CALL DERIVG(N,NW,TMP,TOLD(TMP-1),W(I,15),W(1,IW),
     +              T-TOLD(TMP-1),KOLD)

C ... Compute delay value using the special interpolant
            ELSE
               DO 135 J = 1, N
                  W(J,IW) = W(J,1)
 135           CONTINUE
            END IF
         END IF

 125  CONTINUE

C ... Compute K2

      CALL FCN(N,ND,NW,TEMP,W(1,16),W(1,21),W(1,2))

C ... Compute K3 to KS
      DO 140 K = 3, S

         TEMP = T + CC(K)*H

C ... Compute Y values
         DO 150 I = 1, N
            W(I,16) = AA(K,1)*W(I,1)
            DO 160 J = 2, K - 1
               W(I,16) = W(I,16) + AA(K,J)*W(I,J)
 160        CONTINUE
            W(I,16) = Y(I) + H*W(I,16)
 150     CONTINUE

c .. for assmyptocally correct defect control
         if (k .eq. 8) then
            do 165 i = 1, n
               w(i,12) = w(i,16)
 165        continue
         end if

C ... Compute arguments for delay values
         CALL SIGMA(N,ND,TEMP,W(1,16),W(1,15))

C ... Compute solution delay values
         DO 170 I = 1, MU

            IW = 20 + I

            IF (W(I,15) .GT. TEMP) THEN
               PRINT *, '**** WARNING:APFORM *****'
               PRINT *, 'AdvanceD Delay at Stage', K
               PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
               PRINT *

C ... Cut the stepseize to half and integrate again
               IF (W(I,15)-TEMP .GE. 1.D0) THEN
                  PRINT *, 'ADVANCED TOO MUCH !!!'
                  PRINT *
                  C(18) = C(18) / 2.D0
                  IF (C(18) .LT. C(13)) THEN
                     PRINT *, 'IND = -3,  At T = ', T
                     STOP
                  END IF
                  C(17) = T + C(18)
                  C(24) = C(24) + DFLOAT(K - 2)
                  H = C(18)
                  GOTO 1000
               END IF
            END IF

            SIGN = W(I,15) - T
            TEMP2 = SIGN**2 / (CC(K)*H)

C ... In case there is no vanishing delay
            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +              TOLD,KOLD,W,YINIT)

C ... In case there is vanishing delay
            ELSE
               IF (INDEX .EQ. 0) INDEX = K

C ... compute delay values using extrapolation
               IF (C(29) .EQ. 0.D0) THEN
	          TMP = IDINT(C(27))
                  IF (TMP .EQ. 1) TMP = LOLD
                  CALL INTRPG(N,NW,TMP,TOLD(TMP-1),W(I,15),W(1,IW),
     +                 T-TOLD(TMP-1),KOLD)

C ... Compute delay value using the special interpolant
               ELSE
                  DO 180 L = 1, N
                     W(L,IW) = AA(K,1) * W(L,1)
                     DO 190 J = 2, K-1
                        W(L,IW) = W(L,IW) + AA(K,J) * W(L,J)
 190                 CONTINUE
                     W(L,IW) = Y(L) + (SIGN - TEMP2) * W(L,1)
     +                    + TEMP2 * W(L,IW) / CC(K)
 180              CONTINUE
               END IF
            END IF
 170     CONTINUE

C ... Compute derivative delay terms
         DO 175 I = MU + 1, ND

            IW = 20 + I

            IF (W(I,15) .GT. TEMP) THEN
               PRINT *, '**** WARNING:APFORM *****'
               PRINT *, 'AdvanceD Delay at Stage', K
               PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
               PRINT *

C ... Cut the stepseize to half and integrate again
               IF (W(I,15)-TEMP .GE. 1.D0) THEN
                  PRINT *, 'ADVANCED TOO MUCH !!!'
                  PRINT *
                  C(18) = C(18) / 2.D0
                  IF (C(18) .LT. C(13)) THEN
                     PRINT *, 'IND = -3,  At T = ', T
                     STOP
                  END IF
                  C(17) = T + C(18)
                  C(24) = C(24) + DFLOAT(K - 2)
                  H = C(18)
                  GOTO 1000
               END IF
            END IF

            SIGN = W(I,15) - T

C ... In case there is no vanishing delay
            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +              TOLD,KOLD,W,DYINIT)

C ... In case there is vanishing delay
            ELSE
               IF (INDEX .EQ. 0) INDEX = K

C ... compute delay values using extrapolation
               IF (C(29) .EQ. 0.D0) THEN
	          TMP = IDINT(C(27))
                  IF (TMP .EQ. 1) TMP = LOLD
                  CALL DERIVG(N,NW,TMP,TOLD(TMP-1),W(I,15),W(1,IW),
     +                 T-TOLD(TMP-1),KOLD)

C ... Compute delay value using the special interpolant
               ELSE

                  DO 185 L = 1, N
                     W(L,IW) = AA(K,1) * W(L,1)
                     DO 195 J = 2, K-1
                        W(L,IW) = W(L,IW) + AA(K,J) * W(L,J)
 195                 CONTINUE
                     W(L,IW) = W(L,1) - (W(L,1) - W(L,IW) / CC(K)) *
     +                    (2.D0 * SIGN / (H * CC(K)))
 185              CONTINUE
               END IF
            END IF
 175        CONTINUE

C ... Compute K-value
         CALL FCN(N,ND,NW,TEMP,W(1,16),W(1,21),W(1,K))

 140  CONTINUE

C***** End of the 1st phase *****

C***** Beginning of the 2nd phase *****

C ... Set # of iteration to be 0
      NIT = 0

C ... If index is not 0, then there is a vanishg delay
      IF (INDEX .NE. 0) THEN

C ... Loop for iterative scheme ...
 2000    CONTINUE

         NIT = NIT + 1
         ERRKI = 0.D0

C ... Recompute K-values for stage INDEX+1 to s
         DO 210 K = INDEX, S

            DO 220 I = 1, N
               W(I,13) = W(I,K)
 220        CONTINUE

            TEMP = T + CC(K)*H

C ... Recompute y values
            DO 230 I = 1, N
               W(I,16) = AA(K,1)*W(I,1)
               DO 240 J = 2, K-1
                  W(I,16) = W(I,16) + AA(K,J)*W(I,J)
 240           CONTINUE
               W(I,16) = Y(I) + H*W(I,16)
 230        CONTINUE

c .. for assmyptocally correct defect control

            if (k .eq. 8) then
               do 235 i = 1, n
                  w(i,12) = w(i,16)
 235           continue
            end if

C ... Recompute arguments for delay values
            CALL SIGMA(N,ND,TEMP,W(1,16),W(1,15))

C ... Recompute solution delay values using intrp routine
            DO 250 I = 1, MU

               IW = 20 + I

               IF (W(I,15) .GT. TEMP) THEN
                  PRINT *, '**** WARNING:APFORM *****'
                  PRINT *, '    Advanced Delay at Stage ', K
                  PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
                  PRINT *

C ... Cut the stepseize to half and integrate again
                  IF (W(I,15)-TEMP .GE. 1.D0) THEN
                     PRINT *, 'ADVANCED TOO MUCH !!!'
                     PRINT *
                     C(18) = C(18) / 2.D0
                     IF (C(18) .LT. C(13)) THEN
                        PRINT *, 'IND = -3,  At T = ', T
                        STOP
                     END IF
                     C(17) = T + C(18)
                     C(24) = C(24) + DFLOAT(S-1 + K-INDEX)
                     H = C(18)
                     GOTO 1000
                  END IF
               END IF

               SIGN = W(I,15) - T

C ... In case there is no vanisning delay
               IF (SIGN .LE. 0.D0) THEN
                  CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +                 TOLD,KOLD,W,YINIT)
               ELSE
                  CALL INTRP(N,NW,T,Y,W(I,15),W(1,IW),H,W)
               END IF
 250        CONTINUE

C ... Recompute derivative delay values using intrp routine
            DO 255 I = MU + 1, ND

               IW = 20 + I

               IF (W(I,15) .GT. TEMP) THEN
                  PRINT *, '**** WARNING:APFORM *****'
                  PRINT *, '    Advanced Delay at Stage ', K
                  PRINT *, 'Ti =', TEMP, ' Delay Argument =', W(I,15)
                  PRINT *

C ... Cut the stepseize to half and integrate again
                  IF (W(I,15)-TEMP .GE. 1.D0) THEN
                     PRINT *, 'ADVANCED TOO MUCH !!!'
                     PRINT *
                     C(18) = C(18) / 2.D0
                     IF (C(18) .LT. C(13)) THEN
                        PRINT *, 'IND = -3,  At T = ', T
                        STOP
                     END IF
                     C(17) = T + C(18)
                     C(24) = C(24) + DFLOAT(S-1 + K-INDEX)
                     H = C(18)
                     GOTO 1000
                  END IF
               END IF

               SIGN = W(I,15) - T

C ... In case there is no vanisning delay
               IF (SIGN .LE. 0.D0) THEN
                  CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +                 TOLD,KOLD,W,DYINIT)
               ELSE

                  CALL DERIV(N,NW,T,Y,W(I,15),W(1,IW),H,W)
               END IF
 255        CONTINUE

C ... Recompute K-value

            CALL FCN(N,ND,NW,TEMP,W(1,16),W(1,21),W(1,K))

            DO 260 I = 1, N
               ERRKI = DMAX1(ERRKI, DABS(W(I,13) - W(I,K)))
 260        CONTINUE

 210     CONTINUE
         C(24) = C(24) + DFLOAT(S - INDEX + 1)

        IF ((ERRKI .GT. TOL) .AND. (NIT .LT. MAXNIT)) GOTO 2000

      END IF

c      PRINT *, 'INDEX =', INDEX
c      PRINT *, 'NIT =', NIT
c      PRINT *, 'ERRKI =', ERRKI
c      PRINT *

C***** End of the the 2nd phase *****

C***** Beginning of the 3rd phase *****

C ... Improve derivative values on the mesh point for the interpolation
C ... to archieve C1 continuity
      IF ((INDEX .GT. 0) .AND. (INDEX .LE. 8)) THEN
         TEMP = T + CC(8)*H

C ... Compute Y values
         DO 310 I = 1, N
            W(I,12) = AA(8,1)*W(I,1)
            DO 320 J = 2, 7
               W(I,12) = W(I,12) + AA(8,J)*W(I,J)
 320        CONTINUE
            W(I,12) = Y(I) + H*W(I,12)
 310     CONTINUE

C ... Compute arguments for delay values
         CALL SIGMA(N,ND,TEMP,W(1,12),W(1,15))

C ... Compute solution delay values
         DO 330 I = 1, MU

            IW = 20 + I
            SIGN = W(I,15) - T

            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +                      TOLD,KOLD,W,YINIT)
            ELSE
               CALL INTRP(N,NW,T,Y,W(I,15),W(1,IW),H,W)

            END IF
 330     CONTINUE

C ... Compute derivative delay values
         DO 335 I = MU + 1, ND

            IW = 20 + I
            SIGN = W(I,15) - T

            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +                      TOLD,KOLD,W,DYINIT)
            ELSE

               CALL DERIV(N,NW,T,Y,W(I,15),W(1,IW),H,W)
            END IF
 335     CONTINUE

C ... Recompute K-value
         CALL FCN(N,ND,NW,TEMP,W(1,12),W(1,21),W(1,8))

         C(24) = C(24) + 1.D0
         C(28) = 0.D0
      END IF

C***** End of the 3nd phase *****

C***** Beginning of 4th phase (Improve intrpolant for defect) *****
      do 500 k = 9, s
         temp = t + cc(k)*h

C ... Compute Y values
         call intrp(n,nw,t,y,temp,w(1,16),h,w)

C ... Compute arguments for delay values
         call sigma(n,nd,temp,w(1,16),w(1,15))

C ... Compute solution delay values
         do 510 i = 1, mu

            iw = 20 + i
            sign = w(i,15) - t

            if (sign .le. 0.d0) then
               call evalags(n,nw,iw,lold,c(25),c(26),c(27),w(i,15),
     +                      told,kold,w,yinit)
            else
               call intrp(n,nw,t,y,w(i,15),w(1,iw),h,w)
            end if
 510     continue

C ... Compute derivative delay values
         do 520 i = mu + 1, nd

            iw = 20 + i
            sign = w(i,15) - t

            if (sign .le. 0.d0) then
               call evalagd(n,nw,iw,lold,c(25),c(26),c(27),w(i,15),
     +                      told,kold,w,dyinit)
            else
               call deriv(n,nw,t,y,w(i,15),w(1,iw),h,w)
            end if
 520     continue

C ... Recompute K-value
         call fcn(n,nd,nw,temp,w(1,16),w(1,21),w(1,k))

         c(24) = c(24) + 1.d0
 500  continue

C ... Finally compute solution at T + H in W(*,12)
      DO 410 I = 1, N
         W(I,12) = BB(1)*W(I,1)
         DO 420 J = 2, S
            W(I,12) = W(I,12) + BB(J)*W(I,J)
 420     CONTINUE
         W(I,12) = Y(I) + H*W(I,12)
 410  CONTINUE

c      print *, 'index =', index

c      print '(3d25.14)', (w(1,i),i=1,12)
c      print *
c      print '(3d25.14)', (w(2,i),i=1,12)
c      print *
c      print '(3d25.14)', (w(3,i),i=1,12)
c      print *

      RETURN
      END

C *** End of APFORM_ASM ***

      SUBROUTINE DDVERK(N,NW,MU,OMEGA,IND,LOLD,T,TEND,Y,TOL,TOLD,KOLD
     +     ,W,C,FCN,SIGMA,YINIT,DYINIT)

C************************************************************************
C                                                                       *
C            Description of DDVERK (Experimental Code)                  *
C                                                                       *
C     DDVERK is a driver which invokes delay and discontinuity routines.*
C     This driver first initializes parameters and then calls the DVERK *
C     integration routine, RDMETH. After each step of the integration,  *
C     RDMETH returns to this driver and this driver checks if a         *
C     discontinuity exists in the step. If the existence of a           *
C     discontinuity is suspected, this driver calls DSLOC and           *
C     DSLOC detects the discontinuity point. Finally, if the step is    *
C     accepted, DDVERK calls STORE or STORED to store information       *
C     necessary to construct an interpolant in future.                  *
C                                                                       *
C************************************************************************
C                                                                       *
C                          Calling sequence                             *
C                                                                       *
C     CALL DDVERK(N,NW,MU,OMEGA,IND,LOLD,T,TEND,Y,TOL,TOLD,KOLD,W,C,    *
C    +            FCN,SIGMA,YINIT,DYINIT)                               *
C                                                                       *
C     INTEGER N,NW,IND,LOLD,MU,OMEGA                                    *
C     DOUBLE PRECISION T,TEND,Y(N),TOL,TOLD(LOLD),KOLD(NW,7,LOLD)       *
C     DOUBLE PRECISION W(NW,20+MU+OMEGA),C(*)                          *
C     EXTERNAL FCN,SIGMA,YINIT,DYINIT                                   *
C                                                                       *
C************************************************************************
C                                                                       *
C                       Arguments for DDERK                             *
C                                                                       *
C    N      Number of equations                                         *
C    NW     First dimension of workspace W and queue KOLD               *
C           NW must be at least MAX(N,MU+OMEGA)                         *
C    MU     The number of solution delay terms                          *
C    OMEGA  The number of derivative delay terms                        *
C    IND    Indicator --- same as DVERK except for IND = 7:             *
C           IND is set to be 7 after a discontinuity is found.          *
C    LOLD   Length of queue                                             *
C    T      Independent value                                           *
C    TEND   Value of T to which integration is to be carried out        *
C    Y      Solution vector at T                                        *
C    TOL    Tolerance                                                   *
C    TOLD   Queue of independent values T                               *
C    KOLD   Queue of derivative and solution approximations             *
C           The first dimension of KOLD must be NW                      *
C           The second dimension of KOLD must be at least 7             *
C           The third dimension of KOLD must be LOLD                    *
C    W      Workspace array                                             *
C           The first dimension of W must be NW                         *
C           The second dimension of W must be at least 20 + ND          *
C    C      Communication vector - The dimension must be greater than   *
C           or equal to 33. C(1) to C(24) are the same as DVERK.        *
C    FCN    Name of subroutine for computing derivative                 *
C           (must be supplied by a user)                                *
C    SIGMA  Name of subroutine for computing delay argument             *
C           (must be supplied by a user)                                *
C    YINIT  Name of subroutine for computing initial values on the      *
C           initial interval. (must be supplied by a user)              *
C    DYINIT Name of subroutine for computing derivative of initial      *
C           values on the  initial interval. (must be supplied by       *
C           a user)                                                     *
C                                                                       *
C************************************************************************
C                                                                       *
C                          Communication vector                         *
C                                                                       *
C     C(25) Initial time                                                *
C     C(26) The first index of queue                                    *
C     C(27) The last index of queue                                     *
C     C(28) Flag indicating if K1 needs to be computed in the next      *
C           step: if C(28) = 0.0 then we need to compute K1 in the next *
C           step                                                        *
C     C(29) Flag indicating which of the extrapolation or the special   *
C           interpolant must be used: if C(29) = 0.0 then we use        *
C           extrapolation                                               *
C     C(30) TL values used in discontinuity locating routine DSLOC.     *
C     C(31) TH - TL used in DSLOC                                       *
C     C(32) Flag indicating which of the defect estimate must be used:  *
C              1: 1pt, 2: 2pt, 3: asymptotically valid                  *
C              The default value is 1.                                  *
c     C(33) TS value (location of disconinuous point) used in ddverk    *
C                                                                       *
C                                                                       *
C************************************************************************
C                                                                       *
C                               Work array                              *
C                                                                       *
C     W(*, 1:11)           K1 to K11 of CRK formula                     *
C     W(*, 12)             Solution at the right end point of the       *
C                          each step                                    *
C     W(*, 13)             Work space                                   *
C     W(*, 14)             Estimate of Defect                           *
C     W(*, 15)             Delay argument                               *
C     W(*, 16)             Work space                                   *
C     W(*, 17)             Work space                                   *
C     W(*, 18)             Work space                                   *
C     W(*, 21:20+ND)       Delay values                                 *
C                                                                       *
C************************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,MU,OMEGA,IND,LOLD
      DOUBLE PRECISION T,TEND,TOL

C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),TOLD(LOLD),KOLD(NW,7,LOLD)
      DOUBLE PRECISION W(NW,*),C(*)

C     .. Subroutine Arguments ..
      EXTERNAL FCN,SIGMA,YINIT,DYINIT

C     .. Local Scalars ..
      INTEGER I,IND2,P
      DOUBLE PRECISION TS,K, h2, h3

C    .. Local Logicals ..
      LOGICAL GOBACK

C     .. External Subroutines ..
      EXTERNAL RDMETH,DSLOC,INTRP,STORED,STORE

c .. scalr in common
      integer disc_count

C     .. Common blocks ..
      COMMON/GOBACK/GOBACK
      common/disc_c/disc_count

C     .. Data Statements ..
C ... Set the order of methods to be 6
      DATA P/6/
C ... GOBACK: Logical variable
C             --- if true, then return to main routine in every step
      DATA GOBACK/.FALSE./

C ... Initialize communication array

      IF (IND .EQ. 1) THEN
         IND = 2
         DO 110 I = 1, 9
            C(I) = 0.D0
 110     CONTINUE
         C(32) = 1
      END IF

      IF (IND .EQ. 2) THEN
C ... Set up initial value
         CALL YINIT(N,T,Y)

C ... In the first step, the extrapolation cannot be used. Hence set
C     C(29) to be 1.
         C(29) = 1.D0

         GOBACK = C(9) .EQ. 1.D0
C ... Has to be returned to this routine in every step
         C(9) = 1.D0

C ... Substitue initial time to C(25) and TOLD(1)
         C(25) = T
         TOLD(1) = T

C ... Set the first index of the history queue to be 1
         C(26) = 1.D0

C ... Set the last index of the history queue to be 1
         C(27) = 1.D0

C ... Set flag, which determins if K1 (in the current step) is
C     equal to k8 (in the previous step), to 0.d0.
C     If FLAG(=C(28)) is 0, then K1 != K8.
         C(28) = 0.D0

      END IF

C ... Loop for integration till T = Tend
 1000 CONTINUE

C ... Change the information when IND = 7
      IF (IND .EQ. 7) THEN
         T = C(17)
         DO 210 I = 1, N
            Y(I) = W(I,12)
 210     CONTINUE

C ... Add 1 to the number of successful steps
         C(22) = C(22) + 1.D0

C ... Set the number of successive failures to be 0
         C(23) = 0.D0

C ... If T >= Tend, then end the integration
         IF (T .GE. TEND) THEN
            IND = 3
            C(20) = TEND
            C(21) = 1.D0
            RETURN
         END IF
      END IF

 2000 continue
      CALL RDMETH(N,NW,MU,OMEGA,IND,LOLD,T,TEND,Y,TOLD,KOLD,TOL,W,C,
     +            FCN,SIGMA,YINIT,DYINIT)

C ... Exit successfully if IND = 3
      IF (IND .EQ. 3) RETURN

C ... Error exit if IND < 1 or IND > 8
      IF ((IND .LT. 1) .OR. (IND .GT. 9)) GOTO 500

C ... Check for signal of possible discontinuity -- the signal by
C     an attempt to 1.6 times the step when the immediately preceding
C     step was a failure

      IF ((C(23) .GT. 0.D0) .AND. (IND .EQ. 5)
     +     .AND. (TOL .GE. 1.6d0**P*C(19)/0.9D0)) THEN

C ... Use defect to locate and cross discontinuity if one exist
         CALL DSLOC(N,NW,MU,OMEGA,IND2,LOLD,T,TS,Y,W(1,13),TOLD,KOLD,
     +              TOL,W,C,FCN,SIGMA,YINIT,DYINIT)

C ... Now cross the discontinuity if one has been found (IND2 > 0)
C     and prepare for a restart of the integrator.

         print *, 'ind2 =', ind2

         IF (IND2 .GT. 0) THEN
            PRINT *
            PRINT *, 'DETECT DISCONTINUITY AT T = ', TS

            IF (TS .LT. C(17)) THEN
               PRINT *, '***** ERROR AT DDVERK *****'
               PRINT *, 'TS is less than Tn-1 + H'
               PRINT *
               STOP
            END IF

C ... In case the discontinous point is greater than or equal to Tend,
C     we have to inerpolate tend and store informaiton and then return
C     to main routine.

            IF ((TS .GE. TEND) .OR. (DABS(TEND - TS) .LE. C(10))) THEN
               CALL INTRP(N,NW,T,Y,TEND,W(1,13),C(18),W)
               CALL STORED(N,NW,LOLD,T,TEND,Y,W(1,13),TOLD,KOLD,W,C)
               C(17) = TEND

               DO 310 I = 1, N
                  W(I,12) = W(I,13)
 310           CONTINUE

               write(13, 10000) ts
               disc_count = disc_count + 1

               RETURN
            END IF

c ... Try not to use extration, so extend the stepsize and integrate
c ... Note: c(30) = XL and c(33) = XH - XL
            h3 = ts - t
            h2 = dmin1(h3/1.1d0, c(30) - t)

c            print *, 'h3/h1=', h3/c(18), ', h2-h1/h1 =',
c     +           (h2 - c(18))/c(18), ', hl =', c(30)
            if ( (h3/c(18) .gt. 1.2d0) .and. (dabs(h2 - c(18))/c(18)
     +           .gt. 0.1d0)) then

c ... Save info in case the next step is rejected

               CALL STORE(N,NW,LOLD,C(17),Y,C(26),C(27),TOLD,KOLD,W)
               c(27) = c(27) - 1.d0
               ind = 8
               c(18) = h2
               c(17) = t + c(18)
               c(14) = c(18)
               c(33) = ts
               go to 2000

            end if

C ... If there is a discontinuity, the extrapolation cannot be used.
C     Hence set C(29) to be 1.
            C(29) = 1.D0

            IND = 7
            K = C(23) - 1.D0
            C(4) = C(18) * 2.D0**K
            C(28) = 0.D0
            c(23) = 0.d0

            CALL STORED(N,NW,LOLD,T,TS,Y,W(1,13),TOLD,KOLD,W,C)
            C(17) = TS

            DO 320 I = 1,N
               W(I,12) = W(I,13)
 320        CONTINUE

            write(13, 10000) ts
            disc_count = disc_count + 1

         ELSE

            CALL STORE(N,NW,LOLD,C(17),Y,C(26),C(27),TOLD,KOLD,W)

         END IF

      ELSE IF (IND .EQ. 5) THEN

         CALL STORE(N,NW,LOLD,C(17),Y,C(26),C(27),TOLD,KOLD,W)

      ELSE IF (IND .EQ. 8) THEN

         c(29) = 1.d0
         ind = 7
         k = c(23) - 1.d0
         c(4) = c(18)
c         c(4) = c(18) * 2.d0 ** k
         c(28) = 0.d0
         c(23) = 0.d0
         call intrp(n,nw,t,y,c(33),w(1,13),c(18),w)
         call stored(n,nw,lold,t,c(33),y,w(1,13),told,kold,w,c)
         c(17) = c(33)
         do 350 i = 1, n
            w(i,12) = w(i,13)
 350     continue

         write(13, 10000) c(33)
         disc_count = disc_count + 1

      else if (ind .eq. 9) then

         c(27) = c(27) + 1.d0
         c(23) = 0.d0
         ind = 5
         do 370 i = 1, n
            w(i,12) = kold(i,3,idint(c(27)))
            w(i,8) = kold(i,4,idint(c(27)))
 370     continue
         c(17) = told(idint(c(27)))
         c(18) = c(17) - told(idint(c(27)) - 1.d0)
         c(14) = c(18)
         c(29) = 0.d0

      end if

 500  CONTINUE

      IF (GOBACK) RETURN

      GO TO 1000

10000 format(f20.12)

      END

C *** End of DDVERK ***

      SUBROUTINE DEFECT(N,NW,MU,OMEGA,LOLD,T,TS,Y,DFCT,TOLD,KOLD,W,C,
     +                  FCN,SIGMA,YINIT,DYINIT)

C********************************************************************
C                                                                   *
C     Purpose - This subroutine computes weighted max norm of       *
c               defect at the point TS.                             *
C                                                                   *
C********************************************************************
C                                                                   *
C    N      Number of equation                                      *
C    NW     First dimension of workspace w & queue kold             *
C    MU     Number of solution delay terms                          *
C    OMEGA  Number of derivative delay terms                        *
C    LOLD   Length of queue                                         *
C    T      Independent value (left end point)                      *
C    TS     Independent value at which defect is to be computed     *
C    Y      Dependent values (at T)                                 *
C    DFCT   Max norm of defect at ts (output)                       *
C    TOLD   Queue of independent values                             *
C    KOLD   Queue of k-values and solutions                         *
C    W      Workspace                                               *
C    C      Communication array                                     *
C    FCN    Name of subroutine for computing derivative             *
C    SIGMA  Name of subroutine for computing delay argument         *
C    YINIT  Name of subroutine for computing initial values at      *
C           the initial intervals.                                  *
C    DYINIT Name of subroutine for computing derivative of initial  *
C           values at the initial intervals.                        *
C                                                                   *
C********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,MU,OMEGA,LOLD
      DOUBLE PRECISION T,TS,DFCT

C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),TOLD(LOLD),KOLD(NW,7,LOLD),W(NW,*),C(*)

C     .. Subroutine Arguments ..
      EXTERNAL FCN,SIGMA,YINIT,DYINIT

C     .. Local Scalars ..
      INTEGER I,K,IW,ND
      DOUBLE PRECISION TEMP,SIGN,H

C     .. External Subroutines ..
      EXTERNAL INTRP,EVALAGS,EVALAGD,DERIV

C     .. Intrinsic Functions ..
      INTRINSIC DMAX1,DABS

c      print *, 't =', t, ' ts =', ts

C ... Error check
      IF (TS .LT. T) THEN
         PRINT *, '***** ERROR AT DEFECT *****'
         PRINT *, 'T=', T, ' TS=', TS
         PRINT *, 'TS is less than T !!!'
         PRINT *
         STOP
      END IF

      ND = MU + OMEGA
      H = C(18)

C ... Compute P(TS)
      CALL INTRP(N,NW,T,Y,TS,W(1,2),H,W)

C ... Compute argument of delay value
      CALL SIGMA(N,ND,TS,W(1,2),W(1,15))

C     Compute P(Delay Arugment)
      DO 110 I = 1, MU
         IW = 20 + I

         IF (W(I,15) .GT. TS) THEN
            PRINT *, '**** Warning:DEFECT *****'
            PRINT *, '    ADVANCED DELAY in DEFECT'
            PRINT *, 'T=',TS,'  Delay Argument=',W(I,15)
            PRINT *
         END IF

         SIGN = W(I,15) - T

C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,YINIT)

C ... In case there is vanishing delay
         ELSE
            CALL INTRP(N,NW,T,Y,W(I,15),W(1,IW),H,W)
         END IF
 110  CONTINUE

C ...  Compute P'(Delay Argument)
      DO 115 I = MU + 1, ND

         IW = 20 + I

         IF (W(I,15) .GT. TS) THEN
            PRINT *, '**** Warning:DEFECT *****'
            PRINT *, '    ADVANCED DELAY in DEFECT'
            PRINT *, 'T=',TS,'  Delay Argument=',W(I,15)
            PRINT *
         END IF

         SIGN = W(I,15) - T

C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,DYINIT)

C ... In case there is vanishing delay
         ELSE
            CALL DERIV(N,NW,T,Y,W(I,15),W(1,IW),H,W)
         END IF
 115  CONTINUE

C ... Compute F(TS, P(TS), P(DEL ARG), P'(DELAY ARG))
      CALL FCN(N,ND,NW,TS,W(1,2),W(1,21),W(1,18))
      C(24) = C(24) + 1.D0

C ... Compute P'(TS) and then defect
      CALL DERIV(N,NW,T,Y,TS,W(1,2),H,W)

c      print *, 'derv =', (w(i,2), i = 1, 3)

      DO 120 I = 1, N
         W(I,14) = W(I,18) - W(I,2)
 120  CONTINUE

c      print *, 'def =', (w(i,14), i=1,3)

C ... Calculate the weighted max norm of W(*,14) as specified by
C     the error control indicator C(1)
      TEMP = 0.D0
      IF (C(1) .NE. 1.D0) GO TO 210

C ... Absolute error control
      DO 200 K = 1, N
         TEMP = DMAX1(TEMP,DABS(W(K,14)))
 200  CONTINUE
      GO TO 500
 210  IF (C(1) .NE. 2.D0) GO TO 240

C ... Relative error control
      DO 230 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)/W(K,18)))
 230  CONTINUE
      GO TO 500
 240  IF (C(1) .NE. 3.D0) GO TO 260

C ... Weights are 1/MAX(C(2),ABS(W(K,18)))
      DO 250 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)) / DMAX1(C(2),
     +        DABS(W(K,18))))
 250  CONTINUE
      GO TO 500
 260  IF (C(1) .NE. 4.D0) GO TO 280

C ... Weights are 1/MAX(C(K+30),ABS(W(K,18)))
      DO 270 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)) / DMAX1(C(K+30),
     +        DABS(W(K,18))))
 270  CONTINUE
      GO TO 500
 280  IF (C(1) .NE. 5.D0) GO TO 300

C ... Weights are 1/C(K+30)
      DO 290 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)/C(K+30)))
 290  CONTINUE
      GO TO 500
 300  CONTINUE

C ... Default case - Weights are 1/MAX(1,ABS(Y'(K)))
      DO 310 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)) / DMAX1(1.D0,
     +        DABS(W(K,18))) )
 310  CONTINUE
 500  CONTINUE

      DFCT = TEMP

      RETURN
      END

C *** End of DEFECT ***

      SUBROUTINE DEFECT2(N,NW,MU,OMEGA,LOLD,T,TS,TS2,Y,DFCT,TOLD,KOLD,W,
     +                  C,FCN,SIGMA,YINIT,DYINIT)

C********************************************************************
C                                                                   *
C     Purpose - This subroutine computes weighted max norm of       *
c               defects at the point TS and TS2.                    *
C                                                                   *
C********************************************************************
C                                                                   *
C    N      Number of equation                                      *
C    NW     First dimension of workspace w & queue kold             *
C    MU     Number of solution delay terms                          *
C    OMEGA  Number of derivative delay terms                        *
C    LOLD   Length of queue                                         *
C    T      Independent value (left end point)                      *
C    TS     Independent value at which defect is to be computed     *
C    TS2    Independent value at which defect is to be computed     *
C    Y      Dependent values (at T)                                 *
C    DFCT   Max norm of defect at ts (output)                       *
C    TOLD   Queue of independent values                             *
C    KOLD   Queue of k-values and solutions                         *
C    W      Workspace                                               *
C    C      Communication array                                     *
C    FCN    Name of subroutine for computing derivative             *
C    SIGMA  Name of subroutine for computing delay argument         *
C    YINIT  Name of subroutine for computing initial values at      *
C           the initial intervals.                                  *
C    DYINIT Name of subroutine for computing derivative of initial  *
C           values at the initial intervals.                        *
C                                                                   *
C********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,MU,OMEGA,LOLD
      DOUBLE PRECISION T,TS,TS2,DFCT

C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),TOLD(LOLD),KOLD(NW,7,LOLD),W(NW,*),C(*)

C     .. Subroutine Arguments ..
      EXTERNAL FCN,SIGMA,YINIT,DYINIT

C     .. Local Scalars ..
      INTEGER I,K,IW,ND
      DOUBLE PRECISION TEMP,SIGN,H

C     .. External Subroutines ..
      EXTERNAL INTRP,EVALAGS,EVALAGD,DERIV

C     .. Intrinsic Functions ..
      INTRINSIC DMAX1,DABS

C ... Error check
      IF ((TS .LT. T) .or. (ts2 .lt. t)) THEN
         PRINT *, '***** ERROR AT DEFECT *****'
         PRINT *, 'T=', T, ' TS=', TS, ' TS2 =', TS2
         PRINT *, 'TS (or TS2) is less than T !!!'
         PRINT *
         STOP
      END IF

      ND = MU + OMEGA
      H = C(18)

C ... Compute P(TS)
      CALL INTRP(N,NW,T,Y,TS,W(1,2),H,W)

C ... Compute argument of delay value
      CALL SIGMA(N,ND,TS,W(1,2),W(1,15))

C     Compute P(Delay Arugment)
      DO 110 I = 1, MU
         IW = 20 + I

         IF (W(I,15) .GT. TS) THEN
            PRINT *, '**** Warning:DEFECT *****'
            PRINT *, '    ADVANCED DELAY in DEFECT'
            PRINT *, 'T=',TS,'  Delay Argument=',W(I,15)
            PRINT *
         END IF

         SIGN = W(I,15) - T

C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,YINIT)

C ... In case there is vanishing delay
         ELSE
            CALL INTRP(N,NW,T,Y,W(I,15),W(1,IW),H,W)
         END IF
 110  CONTINUE

C ...  Compute P'(Delay Argument)
      DO 115 I = MU + 1, ND

         IW = 20 + I

         IF (W(I,15) .GT. TS) THEN
            PRINT *, '**** Warning:DEFECT *****'
            PRINT *, '    ADVANCED DELAY in DEFECT'
            PRINT *, 'T=',TS,'  Delay Argument=',W(I,15)
            PRINT *
         END IF

         SIGN = W(I,15) - T

C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,DYINIT)

C ... In case there is vanishing delay
         ELSE
            CALL DERIV(N,NW,T,Y,W(I,15),W(1,IW),H,W)
         END IF
 115  CONTINUE

C ... Compute F(TS, P(TS), P(DEL ARG), P'(DELAY ARG))
      CALL FCN(N,ND,NW,TS,W(1,2),W(1,21),W(1,18))
      C(24) = C(24) + 1.D0

C ... Compute P'(TS) and then defect
      CALL DERIV(N,NW,T,Y,TS,W(1,2),H,W)

      DO 120 I = 1, N
         W(I,14) = W(I,18) - W(I,2)
 120  CONTINUE

C ... Calculate the weighted max norm of W(*,14) as specified by
C     the error control indicator C(1)
      TEMP = 0.D0
      IF (C(1) .NE. 1.D0) GO TO 210

C ... Absolute error control
      DO 200 K = 1, N
         TEMP = DMAX1(TEMP,DABS(W(K,14)))
 200  CONTINUE
      GO TO 500
 210  IF (C(1) .NE. 2.D0) GO TO 240

C ... Relative error control
      DO 230 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)/W(K,18)))
 230  CONTINUE
      GO TO 500
 240  IF (C(1) .NE. 3.D0) GO TO 260

C ... Weights are 1/MAX(C(2),ABS(W(K,18)))
      DO 250 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)) / DMAX1(C(2),
     +        DABS(W(K,18))))
 250  CONTINUE
      GO TO 500
 260  IF (C(1) .NE. 4.D0) GO TO 280

C ... Weights are 1/MAX(C(K+30),ABS(W(K,18)))
      DO 270 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)) / DMAX1(C(K+30),
     +        DABS(W(K,18))))
 270  CONTINUE
      GO TO 500
 280  IF (C(1) .NE. 5.D0) GO TO 300

C ... Weights are 1/C(K+30)
      DO 290 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)/C(K+30)))
 290  CONTINUE
      GO TO 500
 300  CONTINUE

C ... Default case - Weights are 1/MAX(1,ABS(Y'(K)))
      DO 310 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)) / DMAX1(1.D0,
     +        DABS(W(K,18))) )
 310  CONTINUE
 500  CONTINUE

      DFCT = TEMP

C ... Compute P(TS2)
      CALL INTRP(N,NW,T,Y,TS2,W(1,2),H,W)

C ... Compute argument of delay value
      CALL SIGMA(N,ND,TS2,W(1,2),W(1,15))

C     Compute P(Delay Arugment)
      DO 610 I = 1, MU
         IW = 20 + I

         IF (W(I,15) .GT. TS2) THEN
            PRINT *, '**** Warning:DEFECT *****'
            PRINT *, '    ADVANCED DELAY in DEFECT'
            PRINT *, 'T=',TS2,'  Delay Argument=',W(I,15)
            PRINT *
         END IF

         SIGN = W(I,15) - T

C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,YINIT)

C ... In case there is vanishing delay
         ELSE
            CALL INTRP(N,NW,T,Y,W(I,15),W(1,IW),H,W)
         END IF
 610  CONTINUE

C ...  Compute P'(Delay Argument)
      DO 615 I = MU + 1, ND

         IW = 20 + I

         IF (W(I,15) .GT. TS2) THEN
            PRINT *, '**** Warning:DEFECT *****'
            PRINT *, '    ADVANCED DELAY in DEFECT'
            PRINT *, 'T=',TS2,'  Delay Argument=',W(I,15)
            PRINT *
         END IF

         SIGN = W(I,15) - T

C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,LOLD,C(25),C(26),C(27),W(I,15),
     +           TOLD,KOLD,W,DYINIT)

C ... In case there is vanishing delay
         ELSE
            CALL DERIV(N,NW,T,Y,W(I,15),W(1,IW),H,W)
         END IF
 615  CONTINUE

C ... Compute F(TS2, P(TS2), P(DEL ARG), P'(DELAY ARG))
      CALL FCN(N,ND,NW,TS2,W(1,2),W(1,21),W(1,18))
      C(24) = C(24) + 1.D0

C ... Compute P'(TS2) and then defect
      CALL DERIV(N,NW,T,Y,TS2,W(1,2),H,W)

      DO 620 I = 1, N
         W(I,14) = W(I,18) - W(I,2)
 620  CONTINUE

C ... Calculate the weighted max norm of W(*,14) as specified by
C     the error control indicator C(1)
      TEMP = 0.D0
      IF (C(1) .NE. 1.D0) GO TO 710

C ... Absolute error control
      DO 700 K = 1, N
         TEMP = DMAX1(TEMP,DABS(W(K,14)))
 700  CONTINUE
      GO TO 900
 710  IF (C(1) .NE. 2.D0) GO TO 740

C ... Relative error control
      DO 730 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)/W(K,18)))
 730  CONTINUE
      GO TO 900
 740  IF (C(1) .NE. 3.D0) GO TO 760

C ... Weights are 1/MAX(C(2),ABS(W(K,18)))
      DO 750 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)) / DMAX1(C(2),
     +        DABS(W(K,18))))
 750  CONTINUE
      GO TO 900
 760  IF (C(1) .NE. 4.D0) GO TO 780

C ... Weights are 1/MAX(C(K+30),ABS(W(K,18)))
      DO 770 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)) / DMAX1(C(K+30),
     +        DABS(W(K,18))))
 770  CONTINUE
      GO TO 900
 780  IF (C(1) .NE. 5.D0) GO TO 800

C ... Weights are 1/C(K+30)
      DO 790 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)/C(K+30)))
 790  CONTINUE
      GO TO 900
 800  CONTINUE

C ... Default case - Weights are 1/MAX(1,ABS(Y'(K)))
      DO 810 K = 1, N
         TEMP = DMAX1(TEMP, DABS(W(K,14)) / DMAX1(1.D0,
     +        DABS(W(K,18))) )
 810  CONTINUE
 900  CONTINUE

c      print *, 'ts =', ts, ' ts2 =', ts2
c      print *, 'dfct1 =', dfct, ', defct2 =', temp

      dfct = dmax1(dfct,temp)

      RETURN
      END

C *** End of DEFECT2 ***

      SUBROUTINE DEFECTG(N,NW,MU,OMEGA,LOLD,M,T1,TS,H,TOLD,KOLD,W,C,
     +                   FCN,SIGMA,YINIT,DYINIT)

C**********************************************************************
C                                                                     *
C     Purpose - This subroutine computes defect at XS and put them in *
C               W(*,14). The DEFECT computes the defect only on the   *
C               current step using INTRP and DERIV while DEFECTG      *
C               computes the defect in any step using INTRPG and      *
C               DERIVG.                                               *
C                                                                     *
C**********************************************************************
C                                                                     *
C    N      Number of equation                                        *
C    NW     First dimension of workspace W & queue KOLD               *
C    MU     Number of solution delay terms                            *
C    OMEGA  Number of derivative delay terms                          *
C    LOLD   Length of queue                                           *
C    M      Index for the queue, i.e, TOLD(M-1) <= TS <= TOLD(M)      *
C    T1     Independent value (left end point)                        *
C    TS     Independent value at which defect is to be computed       *
C    H      Stepsize, i.e, TOLD(M) - TOLD(M-1)                        *
C    W      Work space                                                *
C    C      Communication array                                       *
C    TOLD   Queue of independent values                               *
C    KOLD   Queue of k-values and solutions                           *
C    FCN    Name of subroutine for computing derivative               *
C    SIGMA  Name of subroutine for computing delay argument           *
C    YINIT  Name of subroutine for computing initial values at        *
C           initial intervals                                         *
C    DYINIT Name of subroutine for computing derivative of initial    *
C           values at initial intervals                               *
C                                                                     *
C**********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,LOLD,M,MU,OMEGA
      DOUBLE PRECISION T1,TS,H

C     .. Array Arguments ..
      DOUBLE PRECISION TOLD(LOLD),KOLD(NW,7,LOLD),W(NW,*),C(*)

C     .. Subroutine Arguments ..
      EXTERNAL FCN,SIGMA,YINIT,DYINIT

C     .. Local Scalars ..
      INTEGER I,ND

C     .. External Subroutines ..
      EXTERNAL INTRPG,EVALAGS,EVALAGD,DERIVG

C ... Error check
      IF (TS .LT. T1) THEN
         PRINT *, '***** ERROR AT DEFECT *****'
         PRINT *, 'X=', T1, ' TS=', TS
         PRINT *, 'TS is less than X !!!'
         PRINT *
         STOP
      END IF

      ND = MU + OMEGA

C ... Compute P(TS)
      CALL INTRPG(N,NW,M,T1,TS,W(1,2),H,KOLD)

C ... Compute argument of delay value
      CALL SIGMA(N,ND,TS,W(1,2),W(1,15))

C ... Compute P(Delay Arguemnt)
      DO 20 I = 1, MU

         IF (W(I,15) .GT. TS) THEN
            PRINT *, '**** WARNING AT DEFECTG *****'
            PRINT *, '    ADVANCED DELAY in DEFECTG'
            PRINT *, 'X=',TS,' Delay Argument =', W(I,15)
            PRINT *
         END IF

         CALL EVALAGS(N,NW,20+I,LOLD,C(25),C(26),C(27),W(I,15),
     +        TOLD,KOLD,W,YINIT)

 20   CONTINUE

C ... Compute P'(Delay Arugment)
      DO 25 I = MU + 1, ND

         IF (W(I,15) .GT. TS) THEN
            PRINT *, '**** WARNING AT DEFECTG *****'
            PRINT *, '    ADVANCED DELAY in DEFECTG'
            PRINT *, 'X=',TS,' Delay Argument =', W(I,15)
            PRINT *
         END IF

         CALL EVALAGD(N,NW,20+I,LOLD,C(25),C(26),C(27),W(I,15),
     +        TOLD,KOLD,W,DYINIT)

 25   CONTINUE

C ... Compute F(TS, P(TS), P(DEL ARG), P'(DELAY ARG))
      CALL FCN(N,ND,NW,TS,W(1,2),W(1,21),W(1,18))
      C(24) = C(24) + 1.D0

C ... Compute P'(TS) and then defect
      CALL DERIVG(N,NW,M,T1,TS,W(1,2),H,KOLD)

      DO 30 I = 1, N
         W(I,14) = W(I,18) - W(I,2)
 30   CONTINUE

      RETURN
      END

C *** End of DEFECTG ***

      SUBROUTINE DERIV(N,NW,T,Y,TOUT,YPOUT,HTRIAL,W)

C**********************************************************************
C                                                                     *
C   Purpose - This subroutine computes derivatives at TOUT            *
C                                                                     *
C**********************************************************************
C                                                                     *
C     N      Number of equations                                      *
C     NW     First dimension of workspace W & queue KOLD              *
C     T      Current independent value                                *
C     Y      Current solution values, i.e., at T                      *
C     TOUT   Independent value at which derivative is computed        *
C     YPOUT  Derivative value at TOUT (OUTPUT)                        *
C     HTRIAL Stepsize                                                 *
C     W      Workspace                                                *
C                                                                     *
C**********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW
      DOUBLE PRECISION T,TOUT,HTRIAL

C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),YPOUT(N),W(NW,12)

C     .. Local Scalars ..
      INTEGER I
      DOUBLE PRECISION TAU,TAU2,TAU3,TAU4,TAU5
      double precision DD0,DD1,DD2,DD3,DD9,DD10,DD11

      TAU = (TOUT - T)/HTRIAL
      TAU2 = TAU * TAU
      TAU3 = TAU2 * TAU
      TAU4 = TAU3 * TAU
      TAU5 = TAU4 * TAU

      DD0 = 4.285714285716394d0 * TAU - 57.85714285714626d0 * TAU2 +
     +     248.5714285713638d0 * TAU3 - 404.9999999997968d0 * TAU4  +
     +     209.9999999998538d0 * TAU5
      DD1 = 1.d0 - 14.53571428571327d0 * TAU + 71.98214285713236d0 *
     +     TAU2 - 154.5714285713923d0 * TAU3  + 146.8749999999520d0 *
     +     TAU4 - 50.74999999997890d0 * TAU5
      DD2 = - 4.285714285716394d0 * TAU + 57.85714285714626d0 * TAU2 -
     +     248.5714285713638d0 * TAU3  + 404.9999999997968d0 * TAU4 -
     +     209.9999999998538d0 * TAU5
      DD3 = 0.5309523809522572D0 * TAU - 7.234523809521525D0
     +  * TAU2  + 31.62857142856110D0 * TAU3 - 53.20833333331760D0
     +  * TAU4  + 29.28333333332581D0 * TAU5
      DD9 =  6.857142857142692D0 * TAU - 88.30476190475958D0
     +  * TAU2  + 348.6476190476093D0 * TAU3 - 498.6666666666516D0
     +  * TAU4  + 231.4666666666593D0 * TAU5
      DD10 =  28.58333333332976D0 * TAU - 225.8083333333001D0
     +  * TAU2  + 617.3999999998936D0 * TAU3 - 700.2916666665310D0
     +  * TAU4  + 280.1166666666077D0 * TAU5
      DD11 = - 17.14999999999986D0 * TAU + 191.5083333333260D0
     +  * TAU2  - 594.5333333332940D0 * TAU3 + 700.2916666666015D0
     +  * TAU4  - 280.1166666666346D0 * TAU5

      DO 110 I = 1, N
         YPOUT(I) = (DD0 * y(i) + DD2 * w(i,12)) / HTRIAL +
     +        DD1 * w(i,1) + DD3 * w(i,8) + DD9 *
     +        w(i,9) + DD10 * w(i,10) + DD11 * w(i,11)
 110  CONTINUE

      RETURN
      END

C *** End of DERIV ***

      SUBROUTINE DERIVG(N,NW,M,T,TOUT,YPOUT,HTRIAL,KOLD)

C**********************************************************************
C                                                                     *
C   Purpose - This subroutine computes derivatives at TOUT            *
C             The DERIV computes the derivative only on the current   *
C             step while DERIVG computes the derivative in any step   *
C                                                                     *
C**********************************************************************
C                                                                     *
C     N      Number of equations                                      *
C     NW     First dimension of workspace W & queue KOLD              *
C     M      Index for the queue, i.e. TOLD(M-1) <= TOUT <= TOLD(M)   *
C     T      Current independent value                                *
C     TOUT   Independent value at which derivative is computed        *
C     YPOUT  Derivative value at TOUT (OUTPUT)                        *
C     HTRIAL Stepsize                                                 *
C     KOLD   Queue of K values and solution values                    *
C                                                                     *
C**********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,M
      DOUBLE PRECISION T,TOUT,HTRIAL

C     .. Array Arguments ..
      DOUBLE PRECISION YPOUT(N),KOLD(NW,7,*)

C     .. Local Scalars ..
      INTEGER I
      DOUBLE PRECISION TAU,TAU2,TAU3,TAU4,TAU5
      DOUBLE PRECISION DD0,DD1,DD2,DD3,DD9,DD10,DD11

      TAU = (TOUT - T)/HTRIAL
      TAU2 = TAU * TAU
      TAU3 = TAU2 * TAU
      TAU4 = TAU3 * TAU
      TAU5 = TAU4 * TAU

      DD0 = 4.285714285716394d0 * TAU - 57.85714285714626d0 * TAU2 +
     +     248.5714285713638d0 * TAU3 - 404.9999999997968d0 * TAU4  +
     +     209.9999999998538d0 * TAU5
      DD1 = 1.d0 - 14.53571428571327d0 * TAU + 71.98214285713236d0 *
     +     TAU2 - 154.5714285713923d0 * TAU3  + 146.8749999999520d0 *
     +     TAU4 - 50.74999999997890d0 * TAU5
      DD2 = - 4.285714285716394d0 * TAU + 57.85714285714626d0 * TAU2 -
     +     248.5714285713638d0 * TAU3  + 404.9999999997968d0 * TAU4 -
     +     209.9999999998538d0 * TAU5
      DD3 = 0.5309523809522572D0 * TAU - 7.234523809521525D0
     +  * TAU2  + 31.62857142856110D0 * TAU3 - 53.20833333331760D0
     +  * TAU4  + 29.28333333332581D0 * TAU5
      DD9 =  6.857142857142692D0 * TAU - 88.30476190475958D0
     +  * TAU2  + 348.6476190476093D0 * TAU3 - 498.6666666666516D0
     +  * TAU4  + 231.4666666666593D0 * TAU5
      DD10 =  28.58333333332976D0 * TAU - 225.8083333333001D0
     +  * TAU2  + 617.3999999998936D0 * TAU3 - 700.2916666665310D0
     +  * TAU4  + 280.1166666666077D0 * TAU5
      DD11 = - 17.14999999999986D0 * TAU + 191.5083333333260D0
     +  * TAU2  - 594.5333333332940D0 * TAU3 + 700.2916666666015D0
     +  * TAU4  - 280.1166666666346D0 * TAU5

      DO 110 I = 1, N
         YPOUT(I) = (DD0 * KOLD(I,1,M) + DD2 * KOLD(I,3,M)) / HTRIAL +
     +        DD1 * KOLD(I,2,M) + DD3 * KOLD(I,4,M) + DD9 *
     +        KOLD(I,5,M) + DD10 * KOLD(I,6,M) + DD11 * KOLD(I,7,M)
 110  CONTINUE

      RETURN
      END

C *** End of DERIVG ***

      SUBROUTINE DSLOC(N,NW,MU,OMEGA,IND2,LOLD,T,TS,Y,YS,TOLD,KOLD,
     +                 TOL,W,C,FCN,SIGMA,YINIT,DYINIT)

C********************************************************************
C                                                                   *
C     Purpose - This routine searches discontinuity point           *
c               using the techiniques proposed by                   *
c               enright et al                                       *
C                                                                   *
C********************************************************************
C                                                                   *
C    N      Number of equation                                      *
C    NW     First dimension of workspace W & queue KOLD             *
C    MU     Number of solution delay terms                          *
C    OMEGA  Number of derivative delay terms                        *
C    IND2   Flag for the discontinuity                              *
C           ... if found IND2 > 0                                   *
C    LOLD   Length of queue                                         *
C    T      Independent value (Tn-1)                                *
C    TS     Discontinous point to be found (output)                 *
C    Y      Solution at T                                           *
C    YS     Solution at TS (Output)                                 *
C    TOLD   Queue of independent values                             *
C    KOLD   Queue of K-values and solutions                         *
C    TOL    Tolerance                                               *
C    W      Workspace                                               *
C    C      Communication array                                     *
C    FCN    Name of subroutine for computing derivative             *
C    SIGMA  Name of subroutine for computing delay argument         *
C    YINIT  Name of subroutine for computing initial values at      *
C           the initial interval                                    *
C    DYINIT Name of subroutine for computing initial values at      *
C           the initial interval                                    *
C                                                                   *
C********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,IND2,LOLD,MU,OMEGA
      DOUBLE PRECISION T,TS,TOL

C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),YS(N),TOLD(LOLD),KOLD(NW,7,LOLD)
      DOUBLE PRECISION W(NW,*),C(*)

C     .. Subroutine Arguments ..
      EXTERNAL FCN,SIGMA,YINIT,DYINIT

C     .. Local Scalars ..
      DOUBLE PRECISION K2,K21,K22,DFCT,EPS,DELTA,TL,TH

C     .. External Subroutines ..
      EXTERNAL DEFECT,INTRP

C     .. Intrisic Functions ..
      INTRINSIC DMAX1,DMIN1

C     .. Data statements ..
C        SET THE VALUES FOR DECIDING IF WE ARE TO THE RIGHT OF Td.
      DATA K21,K22/400.D0, 50.D0/

C ... Set Tl to be Tn-1 + H and Th to be Tn-1 + 2H
      TL = C(17)
      TH = TL + C(18)

C ... To check if there is discontinuity, first we have to check
C     The scale of defect at TH. Thus, seT K2 = 400.
      K2 = K21
      DFCT = C(19)

C ... Set DELTA = 400 * DEFECT at Tstar, WHERE Tn-1 < Tstar < Tn-1 + H
C     DELTA SHOULD BE 10*TOL < DELTA < 1
      DELTA = DMAX1(10.D0*TOL, K2*DFCT)
      DELTA = DMIN1(1.D0, DELTA)

      TS = TH
      CALL DEFECT(N,NW,MU,OMEGA,LOLD,T,TS,Y,DFCT,TOLD,KOLD,W,C,
     +            FCN,SIGMA,YINIT,DYINIT)

C ... Set EPS to be TOL / DEFECT at Tn-1 + 2H for stopping criterion
c      EPS = DMIN1(C(18), TOL / DFCT)
      EPS = DMIN1(C(18)/20.d0, TOL / DFCT)

C ... If DEFECT at Tn-1 + 2H is less than 400 * DEFECT at Tstar,
C     where Tn-1 < Tstar < Tn-1 + H, there is no discontinuity.
      IND2 = -1
      IF (DFCT .LT. DELTA) RETURN

C ... Check the DEFECT at Tn-1 + 3/2 H and then decide which K2 value
C     should be used.
      TS = TL + 0.5D0 * C(18)
      CALL DEFECT(N,NW,MU,OMEGA,LOLD,T,TS,Y,DFCT,TOLD,KOLD,W,C,
     +            FCN,SIGMA,YINIT,DYINIT)

C ... If DEFECT at Tn-1 + 3/2 H > 400 * DEFECT at Tstar,
C     then Tn-1 + H < Td < Tn-1 + 3/2 H
      IF (DFCT .GT. DELTA) THEN
         TH = TS
         EPS = DMIN1(C(18)/20.d0, TOL / DFCT)

C     If DEFECT at Tn-1 + 3/2 H <= 400 * DEFECT at Tstar,
C     then Tn-1 + 3/2 H < Td < Tn-1 + 2H
      ELSE
         TL = TS
         K2 = K22
         DELTA = DMAX1(10.D0*TOL, K2*DFCT)
         DELTA = DMIN1(1.D0, DELTA)
      END IF

      IND2 = 2

C ... Loop --- Bisection search
 1000 CONTINUE

         IF ((TH - TL) .LT. EPS) THEN
            TS = TH
            CALL INTRP(N,NW,T,Y,TS,YS,C(18),W)
            C(30) = TL
            C(31) = TH - TL
            RETURN
         END IF

         TS = (TH + TL) / 2.D0
         CALL DEFECT(N,NW,MU,OMEGA,LOLD,T,TS,Y,DFCT,TOLD,KOLD,W,C,
     +               FCN,SIGMA,YINIT,DYINIT)

         IND2 = IND2 + 1

         IF (DFCT .LT. DELTA) THEN
            TL = TS
         ELSE
            TH = TS
            EPS = DMIN1(C(18)/20.d0, TOL / DFCT)
         END IF

         GOTO 1000

C     End loop --- Bisection search

      END

C *** End of DSLOC ***

      SUBROUTINE EVALAGS(N,NW,IW,LOLD,X0,IFIRST,ILAST,DELARG,TOLD,KOLD,
     +                  W,YINIT)

C********************************************************************
C                                                                   *
C     Purpose - This is the routine that computes delay values      *
c               by using interpolation routine.                     *
C                                                                   *
C********************************************************************
C                                                                   *
C     N       Number of equations                                   *
C     NW      First dimension of workspace W & queue KOLD           *
C     IW      Index for strage space                                *
C     LOLD    Length of the queue                                   *
C     X0      Initial independet value                              *
C     IFIRST  The first index of the queue                          *
C     ILAST   The last index of the queue                           *
C     DELARG  Delay argement                                        *
C     TOLD    Queue of the past independent variables               *
C     KOLD    Queue of the past K-values and soutions at mesh       *
C     W       Workspace matrix (W(*,21:20+I) is used to store delay *
C             values and W(*,21+I:20+ND+I) is used to store         *
C             derivative delay values)                              *
C     YINIT   Name of subroutine for computing initial values       *
C             at initial interval.                                  *
C                                                                   *
C********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,IW,LOLD
      DOUBLE PRECISION X0,IFIRST,ILAST,DELARG

C     .. Array Arguments ..
      DOUBLE PRECISION TOLD(LOLD),KOLD(NW,7,LOLD),W(NW,*)

C     .. Subroutine Arguments ..
      EXTERNAL YINIT

C     .. Local Scalars ..
      INTEGER J,M
      DOUBLE PRECISION H

C     .. External Subroutines ..
      EXTERNAL TFIND,INTRPG

C ... Check whthere the delay is before initial point
      IF (DELARG .LE. X0) THEN
         CALL YINIT(N,DELARG,W(1,IW))

C ... Otherwise find the index in the queue for the delay function
      ELSE

         CALL TFIND(LOLD,TOLD,DELARG,IFIRST,ILAST,M)

C ... In case the data is on the mesh
         IF (TOLD(M) .EQ. DELARG) THEN
            DO 120 J = 1, N
               W(J,IW) = KOLD(J,3,M)
 120        CONTINUE

C ... Otherwise evaluate delay value using intrp routine
C ... In case the index is equal to the last index
         ELSE IF (M .EQ. LOLD) THEN
c         IF (M .EQ. LOLD) THEN
            H = TOLD(1) - TOLD(M)
            CALL INTRPG(N,NW,1,TOLD(M),DELARG,W(1,IW),H,KOLD)

C ... Otherwise
         ELSE
            H = TOLD(M+1) - TOLD(M)
            CALL INTRPG(N,NW,M+1,TOLD(M),DELARG,W(1,IW),H,KOLD)
         END IF
      END IF

      RETURN
      END

C *** End of EVALAGS ***

      SUBROUTINE EVALAGD(N,NW,IW,LOLD,X0,IFIRST,ILAST,DELARG,TOLD,KOLD,
     +                  W,DYINIT)

C********************************************************************
C                                                                   *
C     Purpose - This is the routine that computes delay values      *
c               by using interpolation routine.                     *
C                                                                   *
C********************************************************************
C                                                                   *
C     N       Number of equations                                   *
C     NW      First dimension of workspace W & queue KOLD           *
C     IW      Index of storage space                                *
C     LOLD    Length of the queue                                   *
C     X0      Initial independet value                              *
C     IFIRST  The first index of the queue                          *
C     ILAST   The last index of the queue                           *
C     DELARG  Delay argement                                        *
C     TOLD    Queue of the past independent variables               *
C     KOLD    Queue of the past K-values and soutions at mesh       *
C     W       Workspace matrix (W(*,21:20+I) is used to store delay *
C             values and W(*,21+I:20+ND+I) is used to store         *
C             derivative delay values)                              *
C     DYINIT  Name of subroutine for computing derivative initial   *
C             values at initial interval.                           *
C                                                                   *
C********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,IW,LOLD
      DOUBLE PRECISION X0,IFIRST,ILAST,DELARG

C     .. Array Arguments ..
      DOUBLE PRECISION TOLD(LOLD),KOLD(NW,7,LOLD),W(NW,*)

C     .. Subroutine Arguments ..
      EXTERNAL DYINIT

C     .. Local Scalars ..
      INTEGER J,M
      DOUBLE PRECISION H

c     .. arrays in common
c      double precision c(40)
c      common /comm/c

C     .. External Subroutines ..
      EXTERNAL TFIND,DERIVG

C ... Check whthere the delay is before initial point
      IF (DELARG .LE. X0) THEN
         CALL DYINIT(N,DELARG,W(1,IW))

c      if ((c(40) .eq. 1.d0) .and. (delarg .le. x0)) then
c         call dyinit(n,delarg,w(1,iw))
c      else if ((c(40) .eq. 2.d0) .and. (delarg .lt. x0)) then
c         call dyinit(n,delarg,w(1,iw))

C ... Otherwise find the index in the queue for the delay function
      ELSE

         CALL TFIND(LOLD,TOLD,DELARG,IFIRST,ILAST,M)

C ... In case the data is on the mesh
         IF (TOLD(M) .EQ. DELARG) THEN
            DO 120 J = 1, N
               W(J,IW) = KOLD(J,4,M)
 120        CONTINUE

C ... Otherwise evaluate delay value using intrp routine
c ... In case the index is equal to the last index
         ELSE IF (M .EQ. LOLD) THEN
c         IF (M .EQ. LOLD) THEN
            H = TOLD(1) - TOLD(M)
            CALL DERIVG(N,NW,1,TOLD(M),DELARG,W(1,IW),H,KOLD)

C ... Otherwise
         ELSE
            H = TOLD(M+1) - TOLD(M)
            CALL DERIVG(N,NW,M+1,TOLD(M),DELARG,W(1,IW),H,KOLD)
         END IF
      END IF

      RETURN
      END

C *** End of EVALAGD ***

      SUBROUTINE INTRP(N,NW,T,Y,TOUT,YPOUT,HTRIAL,W)

C**********************************************************************
C                                                                     *
C   Purpose - This subroutine computes solutions at TOUT              *
C                                                                     *
C**********************************************************************
C                                                                     *
C     N      Number of equations                                      *
C     NW     First dimension of workspace W & queue KOLD              *
C     T      Current independent value                                *
C     Y      Current solution values, i.e., at T                      *
C     TOUT   Independent value at which solution is computed          *
C     YPOUT  Solution  at TOUT (Output)                               *
C     HTRIAL Stepsize                                                 *
C     W      Workspace                                                *
C                                                                     *
C**********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW
      DOUBLE PRECISION T,TOUT,HTRIAL

C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),YPOUT(N),W(NW,12)

C     .. Local Scalars ..
      INTEGER I
      DOUBLE PRECISION TAU,TAU2,TAU3,TAU4,TAU5,TAU6
      double precision d0,d1,d2,d3,d9,d10,d11
c     double precision bb(11)

      TAU = (TOUT - T)/HTRIAL
      TAU2 = TAU * TAU
      TAU3 = TAU2 * TAU
      TAU4 = TAU3 * TAU
      TAU5 = TAU4 * TAU
      TAU6 = TAU5 * TAU

      D0 = 1.D0 + 2.142857142858197D0 * TAU2 - 19.28571428571542D0
     + * TAU3 + 62.14285714284096D0 * TAU4 - 80.99999999995936D0
     + * TAU5 + 34.99999999997563D0 * TAU6
      D1 = TAU - 7.267857142856635D0 * TAU2 + 23.99404761904412D0
     + * TAU3 - 38.64285714284807D0 * TAU4 + 29.37499999999040D0
     + * TAU5 - 8.458333333329817D0 * TAU6
      D2 =  - 2.142857142858197D0 * TAU2 + 19.28571428571542D0 * TAU3
     + - 62.14285714284096D0 * TAU4 + 80.99999999995936D0 * TAU5 -
     + 34.99999999997563D0 * TAU6
      D3 = 0.2654761904761286D0 * TAU2 - 2.411507936507175D0 * TAU3
     + + 7.907142857140276D0 * TAU4 - 10.64166666666352D0 * TAU5
     + + 4.880555555554301D0 * TAU6
      D9 = 3.428571428571346D0 * TAU2  - 29.43492063491986D0
     + * TAU3  + 87.16190476190233D0 * TAU4 - 99.73333333333032D0
     + * TAU5  + 38.57777777777655D0 * TAU6
      D10 = 14.29166666666488D0 * TAU2  - 75.26944444443335D0
     + * TAU3  + 154.3499999999734D0 * TAU4 - 140.0583333333062D0
     + * TAU5  + 46.68611111110129D0 * TAU6
      D11 = - 8.574999999999928D0 * TAU2  + 63.83611111110867D0
     + * TAU3  - 148.6333333333235D0 * TAU4 + 140.0583333333203D0
     + * TAU5  - 46.68611111110577D0 * TAU6

      DO 110 I = 1, N
         YPOUT(I) = D0 * y(i) + D2 * w(i,12) + HTRIAL * (D1
     +        * w(i,1) + D3 * w(i,8) + D9 * w(i,9) + D10
     +        * w(i,10) + D11 * w(i,11))
 110  CONTINUE

      RETURN
      END

C *** End of INTRP ***

      SUBROUTINE INTRPG(N,NW,M,T,TOUT,YPOUT,HTRIAL,KOLD)

C**********************************************************************
C                                                                     *
C   Purpose - This subroutine computes solutions at TOUT              *
C             The INTRP computes the solution only on the current     *
C             step while INTRPG computes the solution in any step     *
C                                                                     *
C**********************************************************************
C                                                                     *
C     N      Number of equations                                      *
C     NW     First dimension of workspace W & queue KOLD              *
C     M      Index for the queue, i.e. TOLD(M-1) <= TOUT <= TOLD(M)   *
C     T      Current independent value                                *
C     TOUT   Independent value at which solution is computed          *
C     YPOUT  Solution at TOUT (OUTPUT)                                *
C     HTRIAL Stepsize                                                 *
C     KOLD   Queue of K values and solution values                    *
C                                                                     *
C**********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,M
      DOUBLE PRECISION T,TOUT,HTRIAL

C     .. Array Arguments ..
      DOUBLE PRECISION YPOUT(N),KOLD(NW,7,*)

C     .. Local Scalars ..
      INTEGER I
      DOUBLE PRECISION TAU,TAU2,TAU3,TAU4,TAU5,TAU6
      DOUBLE PRECISION D0,D1,D2,D3,D9,D10,D11
C
      TAU = (TOUT - T)/HTRIAL
      TAU2 = TAU * TAU
      TAU3 = TAU2 * TAU
      TAU4 = TAU3 * TAU
      TAU5 = TAU4 * TAU
      TAU6 = TAU5 * TAU

      D0 = 1.D0 + 2.142857142858197D0 * TAU2 - 19.28571428571542D0
     + * TAU3 + 62.14285714284096D0 * TAU4 - 80.99999999995936D0
     + * TAU5 + 34.99999999997563D0 * TAU6
      D1 = TAU - 7.267857142856635D0 * TAU2 + 23.99404761904412D0
     + * TAU3 - 38.64285714284807D0 * TAU4 + 29.37499999999040D0
     + * TAU5 - 8.458333333329817D0 * TAU6
      D2 =  - 2.142857142858197D0 * TAU2 + 19.28571428571542D0 * TAU3
     + - 62.14285714284096D0 * TAU4 + 80.99999999995936D0 * TAU5 -
     + 34.99999999997563D0 * TAU6
      D3 = 0.2654761904761286D0 * TAU2 - 2.411507936507175D0 * TAU3
     + + 7.907142857140276D0 * TAU4 - 10.64166666666352D0 * TAU5
     + + 4.880555555554301D0 * TAU6
      D9 = 3.428571428571346D0 * TAU2  - 29.43492063491986D0
     + * TAU3  + 87.16190476190233D0 * TAU4 - 99.73333333333032D0
     + * TAU5  + 38.57777777777655D0 * TAU6
      D10 = 14.29166666666488D0 * TAU2  - 75.26944444443335D0
     + * TAU3  + 154.3499999999734D0 * TAU4 - 140.0583333333062D0
     + * TAU5  + 46.68611111110129D0 * TAU6
      D11 = - 8.574999999999928D0 * TAU2  + 63.83611111110867D0
     + * TAU3  - 148.6333333333235D0 * TAU4 + 140.0583333333203D0
     + * TAU5  - 46.68611111110577D0 * TAU6

      DO 110 I = 1, N
         YPOUT(I) = D0 * KOLD(I,1,M) + D2 * KOLD(I,3,M) + HTRIAL * (D1
     +        * KOLD(I,2,M) + D3 * KOLD(I,4,M) + D9 * KOLD(I,5,M) + D10
     +        * KOLD(I,6,M) + D11 * KOLD(I,7,M))
 110  CONTINUE

      RETURN
      END

      SUBROUTINE RDMETH(N,NW,MU,OMEGA,IND,LOLD,T,TEND,Y,TOLD,KOLD,TOL,
     +                  W,C,FCN,SIGMA,YINIT,DYINIT)

C**********************************************************************
C                                                                     *
C     Purpose - This routine sets up necessary data such as           *
C               stepsize, minimum stepsize, etc. Then call APFORM to  *
C               integrates one step by RK-methods for delay DEs.      *
C               This subroutine is based on DVERK.                    *
C                                                                     *
C**********************************************************************
C                                                                     *
C     N      Number of equations                                      *
C     NW     First dimension of workspace W & queue KOLD              *
C     MU     Number of solution delay terms                           *
C     OMEGA  Number of derivative delay terms                         *
C     IND    Indicator --- same as DVERK except IND = 7:              *
C            IND is set to be 7 after a discontinuity is found.       *
C     LOLD   The length of queue                                      *
C     T      Current independent value                                *
C     Y      Current solution, i.e., at T                             *
C     TOLD   Queue of T values                                        *
C     KOLD   Queue oF K values and solution values                    *
C     TOL    Tolerance                                                *
C     W      Workspace                                                *
C     C      Communication vector                                     *
C     FCN    Name of derivative routine                               *
C     SIGMA  Name of subroutine for computing delay argument          *
C     YINIT  Name of subroutine for computing initial values at       *
C            initial intervals.                                       *
C     DYINIT Name of subroutine for computing derivative of initial   *
C            values at initial intervals.                             *
C                                                                     *
C**********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,IND,LOLD,MU,OMEGA
      DOUBLE PRECISION T,TEND,TOL

C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),TOLD(LOLD),KOLD(NW,7,LOLD),W(NW,*),C(*)

C     .. Subroutine Arguments ..
      EXTERNAL FCN,SIGMA,YINIT,DYINIT

C     .. Arrays in Common ..
      DOUBLE PRECISION AA(11,11),BB(11),CC(11)

C     .. Local Scalars ..
      INTEGER I,J,K,P,SSTAR,ND
      DOUBLE PRECISION H,HSCALE,TAUS,TAUS2,TEMP

C     .. External Subroutines ..
      EXTERNAL APFORM,DEFECT,EVALAGS,EVALAGD

C     .. Intrisic Functions ..
      INTRINSIC DABS,DFLOAT,DMIN1,DMAX1,DSIGN

C     .. Common blocks ..
      COMMON/METH/AA,BB,CC

C     ******************************************************************
C     * BEGIN INITIALIZATION, PARAMETER CHECKING, INTERRUPT RE-ENTRIES *
C     ******************************************************************

C ... Abort if IND out of range 1 to 7
      IF (IND.LT.1 .OR. IND.GT.8) GO TO 3000

      SSTAR = 11
      P = 6
      TAUS = .64D0
      TAUS2 = 0.85d0
      if (C(32) .eq. 3) TAUS = .859d0

      ND = MU + OMEGA

C ... Cases - Initial entry, normal re-entry, interrupt re-entries
      GO TO (100, 100, 200, 1111, 2222, 2222, 200, 2500), IND

C ... Case 1 - Initial entry (IND .EQ. 1 OR 2)
C ... Abort if N  > NW or TOL <= 0
 100  IF (N.GT.NW .OR. TOL.LE.0.D0) GO TO 3000

C ... Set up coefficients of RK-formula
      CC(1) = 0.D0
      CC(2) = 0.1D0
      CC(3) = 0.16D0
      CC(4) = 0.2040816326530613D0
      CC(5) = 0.6111D0
      CC(6) = 0.7778D0
      CC(7) = 0.1D1
      CC(8) = 0.1D1
      CC(9) = 0.5D0
      CC(10) = 0.1428571428571428D0
      CC(11) = 0.2857142857142857D0

      AA(2,1) = 0.1D0
      AA(3,1) = 0.32D-1
      AA(3,2) = 0.128D0
      AA(4,1) = 0.388656087174562D-1
      AA(4,2) = 0.9349845727545482D-1
      AA(4,3) = 0.7171756666015026D-1
      AA(5,1) = 0.9182106096978206D0
      AA(5,2) = -0.2806349313292671D0
      AA(5,3) = -0.4995012311203205D1
      AA(5,4) = 0.4968536632834652D1
      AA(6,1) = -0.1172309621097296D1
      AA(6,2) = 0.2577443286829264D0
      AA(6,3) = 0.7019866335902897D1
      AA(6,4) = -0.5919067674431333D1
      AA(6,5) = 0.5915666309428068D0
      AA(7,1) = 0.2052791398807182D1
      AA(7,2) = -0.2063792647412788D0
      AA(7,3) = -0.1319762903121892D2
      AA(7,4) = 0.1236821894967594D2
      AA(7,5) = -0.7279931113442874D0
      AA(7,6) = 0.7109910588213573D0
      AA(8,1) = 0.9752588552523605D-1
      AA(8,2) = 0.D0
      AA(8,3) = -0.4526249411455864D0
      AA(8,4) = 0.8022821731202383D0
      AA(8,5) = 0.2301388037778522D0
      AA(8,6) = 0.24584593549251D0
      AA(8,7) = 0.7683214322974987D-1
      AA(9,1) = 0.1089434604424868D0
      AA(9,2) = 0.D0
      AA(9,3) = -0.5148608705527948D0
      AA(9,4) = 0.8462848943497683D0
      AA(9,5) = 0.6712861232697209D-1
      AA(9,6) = -0.5132033903415767D-2
      AA(9,7) = -0.3361406266302458D-1
      AA(9,8) = 0.312500000000111D-1
      AA(10,1) = 0.8592255366477703D-1
      AA(10,2) = 0.D0
      AA(10,3) = -0.1524656894033379D0
      AA(10,4) = 0.2543389082003915D0
      AA(10,5) = 0.3082416724957052D-1
      AA(10,6) = 0.1449348604485854D-1
      AA(10,7) = -0.3149720157394562D-2
      AA(10,8) = -0.1427976438387382D-2
      AA(10,9) = -0.8567858630333475D-1
      AA(11,1) = 0.4859699015987465D-1
      AA(11,2) = 0.D0
      AA(11,3) = -0.573655693811402D-1
      AA(11,4) = 0.9353712191970651D-1
      AA(11,5) = 0.5261500935769923D-2
      AA(11,6) = -0.3816683388062108D-2
      AA(11,7) = -0.5124089070117709D-2
      AA(11,8) = 0.5665928935000797D-2
      AA(11,9) = 0.4744196280052546D-2
      AA(11,10) = 0.1942148893232018D0

      DO 110 J = 1, 7
         BB(J) = AA(8,J)
 110  CONTINUE
      DO 120 J = 8, 11
         BB(J) = 0.D0
 120  CONTINUE

      IF (IND .EQ. 2) GO TO 130

C ... Initial entry without options (IND .EQ. 1)
C ... Set C(1) to C(9) equal to 0
      DO 140 K = 1, 9
         C(K) = 0.D0
 140  CONTINUE
      GO TO 150

 130  CONTINUE
C ... Initial entry with options (IND .EQ. 2)
C ... Make C(1) to C(9) non-negative
      DO 160 K = 1, 9
         C(K) = DABS(C(K))
 160   CONTINUE

C ... Make floor values non-negative if they are to be used
      IF (C(1).NE.4.D0 .AND. C(1).NE.5.D0) GO TO 180
      DO 170 K = 1, N
         C(K+30) = DABS(C(K+30))
 170  CONTINUE
 180  CONTINUE
 150  CONTINUE

C ... Initialize rreb, dwarf, prev tend, flag, counts
      C(10) = 16.D0**(-13)
      C(11) = 1.D-28

C ... Set previous TEND initially to initial value of T
      C(20) = T

      DO 190 K = 21, 24
         C(K) = 0.D0
 190  CONTINUE
      GO TO 300

C ... Case 2 - Normal re-entry (IND .EQ. 3) or restart after
C              discontinuity (IND .EQ. 7)
C ....Abort if TEND reached, and either T changed or TEND not
 200  IF (C(21) .NE. 0.D0 .AND.
     +     (T. NE. C(20) .OR. TEND .EQ. C(20))) GO TO 3000

C ... Re-initialize flag
      C(21) = 0.D0
      GO TO 300

C ... Case 3 - Re-entry following an interrupt (IND .EQ. 4 TO 6)
C              transfer control to the appropriate re-entry point
C              This has already been handled by the computed go to
C              end case
 300   CONTINUE

C ... End initialization, etc.

C     ******************************************************************
C     * LOOP THROUGH THE FOLLOWING 4 STAGES, ONCE FOR EACH TRIAL  STEP *
C     * UNTIL THE OCCURRENCE OF ONE OF THE FOLLOWING		       *
C     *    (A) THE NORMAL RETURN (WITH IND .EQ. 3) ON REACHING TEND IN *
C     *        STAGE 4						       *
C     *    (B) AN ERROR RETURN (WITH IND .LT. 0) IN STAGE 1 OR STAGE 4 *
C     *    (C) AN INTERRUPT RETURN (WITH IND  .EQ.  4,	5  OR  6),  IF *
C     *        REQUESTED, IN STAGE 1 OR STAGE 4 		       *
C     ******************************************************************
C
99999 CONTINUE

C
C	 ***************************************************************
C	 * STAGE 1 - PREPARE - DO CALCULATIONS OF  HMIN,  HMAX,  ETC., *
C	 * AND SOME PARAMETER  CHECKING,  AND  END  UP	WITH  SUITABLE *
C	 * VALUES OF HMAG, TTRIAL AND HTRIAL IN PREPARATION FOR TAKING *
C	 * AN INTEGRATION STEP. 				       *
C	 ***************************************************************
C
C***********ERROR RETURN (WITH IND=-1) IF NO OF FCN EVALS TOO GREAT
      IF (C(7).EQ.0.D0 .OR. C(24).LT.C(7)) GO TO 400
      IND = -1
      RETURN
 400  CONTINUE

C ... Calculate slope (adding 1 to no of fcn evals) if IND .NE. 6
      IF (IND .EQ. 6) GO TO 410

C ... If C(28) is equal to 1.D0, then K1(in this step) is equal to K8
C     (in the previous step). Thus, no need to compute K1 here.

      IF (C(28) .EQ. 1.D0) THEN
         DO 420 I = 1, N
            W(I,1) = W(I,8)
 420     CONTINUE
      ELSE

         CALL SIGMA(N,ND,T,Y,W(1,15))

         DO 430 I = 1, MU
            IF (W(I,15) .GT. T) THEN
               PRINT *, '***** Warning: RDMETH *****'
               PRINT *, ' ADVANCED DELAY AT STAGE  1'
               PRINT *, 'T= ', T, ' Delay Argument =', W(I,15)

            ELSE
               CALL EVALAGS(N,NW,20+I,LOLD,C(25),C(26),C(27),W(I,15),
     +                     TOLD,KOLD,W,YINIT)
            END IF

 430     CONTINUE

         DO 435 I = MU + 1, ND
            IF (W(I,15) .GT. T) THEN
               PRINT *, '***** Warning: RDMETH *****'
               PRINT *, ' ADVANCED DELAY AT STAGE  1'
               PRINT *, 'T= ', T, ' Delay Argument =', W(I,15)

            ELSE

               CALL EVALAGD(N,NW,20+I,LOLD,C(25),C(26),C(27),W(I,15),
     +                     TOLD,KOLD,W,DYINIT)

            END IF

 435     CONTINUE

         CALL FCN(N,ND,NW,T,Y,W(1,21),W(1,1))
         C(24) = C(24) + 1.D0
      END IF

 410  CONTINUE

C ... Calculate HMIN - Use default unless value prescribed
      C(13) = C(3)
      IF (C(3) .NE. 0.D0) GO TO 440
      C(13) = 1.D-8
c      C(13) = 1.D-6
 440  CONTINUE

C ... Calculate SCALE - Use default unless value prescribed
      C(15) = C(5)
      IF (C(5) .EQ. 0.D0) C(15) = 1.D0

C ... Calculate HMAX - Consider 4 cases
C ... Case 1 both HMAX and SCALE prescribed
      IF (C(6).NE.0.D0 .AND. C(5).NE.0.D0)
     +     C(16) = DMIN1(C(6), 2.D0/C(5))

C ... Case 2 - HMAX prescribed, but SCALE not
      IF (C(6).NE.0.D0 .AND. C(5).EQ.0.D0) C(16) = C(6)

C ... Case 3 - HMAX not prescribed, but SCALE is
      IF (C(6).EQ.0.D0 .AND. C(5).NE.0.D0) C(16) = 2.D0/C(5)

C ... Case 4 - Neither HMAX nor SCALE is provided
      IF (C(6).EQ.0.D0 .AND. C(5).EQ.0.D0) C(16) = 1.5D0

C ... Error return (with IND=-2) if HMIN > HMAX
      IF (C(13) .LE. C(16)) GO TO 450
      IND = -2
      RETURN
 450  CONTINUE

C ... Calculate preliminary HMAG - Consider 3 cases
      IF (IND .GT. 2 .AND. IND .LT. 7) GO TO 460

C ... Case 1 - Initial entry - use prescribed value of HSTART, if
C              any, else default
      C(14) = C(4)

C ... If the delay vanishes at the initial time, then set the stepsize
C     small.
      if (omega .ge. 1) then
         CALL SIGMA(N,ND,T,Y,W(1,15))
         DO 470 I = 1, ND
            IF (((W(I,15) - T) .GT. -C(10))) THEN
               IF ( C(4) .EQ.  0.D0)  C(14) = 1.D2*C(13)
            END IF
 470     CONTINUE
      end if

      IF (C(14) .EQ. 0.D0) C(14) = C(16)*TOL**(1.D0/DFLOAT(P))

      if ((C(14) .gt. 0.d0) .and. (c(14) .le. 1.d-3)) then
         hscale = 50.d0
      else
         HSCALE = 2.D0
      end if

      GO TO 480

 460  IF (C(23) .GT. 1.D0) GO TO 490
      IF (C(23) .EQ. 1.D0) HSCALE = 2.D0

C ... Case 2 - After a successful step, or at most  one  failure,
C              Use MIN(2, .9*(TOL/EST)**(1/P))*HMAG, but avoid possible
C              overflow. Then avoid reduction by more than half.

      TEMP = HSCALE * C(14)
      IF (TOL .LT. (HSCALE/.9D0)**P*C(19))
     +     TEMP = .9D0*(TOL/C(19))**(1.D0/DFLOAT(P))*C(14)

      C(14) = DMAX1(TEMP, .5D0*C(14))
      GO TO 480
 490  CONTINUE

C ... Case 3 - After two or more successive failures
      C(14) = .5D0*C(14)
 480  CONTINUE

C ... Check against HMAX
      C(14) = DMIN1(C(14), C(16))

C ... Check against HMIN
      C(14) = DMAX1(C(14), C(13))

C***********INTERRUPT NO 1 (WITH IND=4) IF REQUESTED
      IF (C(8) .EQ. 0.D0) GO TO 1111
      IND = 4
      RETURN

C ... Resume here on re-entry with IND = 4   ........RE-ENTRY..
 1111 CONTINUE

C ... Calculate HMAG, TTRIAL - Depending on preliminary HMAG, TEND
      IF (C(14) .GE. DABS(TEND - T)) GO TO 500

C ... Do not step more than half way to TEND
      C(14) = DMIN1(C(14), .5D0*DABS(TEND - T))
      C(17) = T + DSIGN(C(14), TEND - T)
      GO TO 510
 500  CONTINUE

C ... Hit TEND exactly
      C(14) = DABS(TEND - T)
      C(17) = TEND
 510  CONTINUE

C ... Calculate HTRIAL
      C(18) = C(17) - T

C ... END STAGE 1

C*********************************************************************
C   Stage 2 - Calculate YTRIAL using the subroutine APFORM
 2500 continue

      H = C(18)

      if ((C(32) .eq. 1) .or. (C(32) .eq. 2)) then

         CALL APFORM(N,NW,MU,OMEGA,LOLD,SSTAR,T,Y,H,TOLD,KOLD,TOL,W,C,
     +            FCN,SIGMA,YINIT,DYINIT)

      else if (C(32) .eq. 3) then

         CALL APFORM_ASM(N,NW,MU,OMEGA,LOLD,SSTAR,T,Y,H,TOLD,KOLD,TOL,W
     +        ,C,FCN,SIGMA,YINIT,DYINIT)

      else

         print *, 'Parameter Error in C(32)'
         stop

      end if

C ... Add SSTAR to the no of fcn evals
          C(24) = C(24) + DFLOAT(SSTAR - 1)

C ... End stage 2
C*********************************************************************

C*********************************************************************
C   Stage 3 - Calculate the error estimate EST. First calculate
C   the DEFECT at TAUS and use the norm of this vector to control
C   the stepsize. Note that W(*,2) and W(*,14) are used for temp
C   storage.

          IF ((C(32) .eq. 1) .or. (C(32) .eq. 3)) then

             CALL DEFECT(N,NW,MU,OMEGA,LOLD,T,T+TAUS*H,Y,C(19),TOLD,KOLD
     +            ,W,C,FCN,SIGMA,YINIT,DYINIT)

          else if (C(32) .eq. 2) then

             CALL DEFECT2(N,NW,MU,OMEGA,LOLD,T,T+TAUS*H,T+TAUS2*H,Y,
     +            C(19),TOLD,KOLD,W,C,FCN,SIGMA,YINIT,DYINIT)

          else

             print *, 'Parameter Error in C(32)'
             stop

          end if

C ... END Stage 3
C*********************************************************************

C*********************************************************************
C   Stage 4 - Make decisions.

C ... Set IND=5 if step acceptable, else set IND=6

         if (c(19) .le. tol) then
c         if ((dabs(w(1,14)) .le. tol) .and.
c     +        (dmax1(dabs(w(2,14)), dabs(w(3,14))) .le. 1.d-4)) then

            if (ind .ne. 8) ind = 5
         else
            if (ind .ne. 8) then
               ind = 6
            else
               ind = 9
            end if
         end if

C ... If accepted then, we can use extrap in the next step,so C(29) = 0
      if ((ind .eq. 8) .or. (ind .eq. 9)) return
      IF (IND .EQ. 5)  C(29) = 0.D0

C***********INTERRUPT NO 2 IF REQUESTED
      IF (C(9) .EQ. 0.D0) GO TO 2222
      RETURN

C ... Resume here on re-entry with IND = 5 OR 6   ...RE-ENTRY..
 2222 CONTINUE

      IF (IND .EQ. 6) GO TO 600

C ... Step accepted (IND .EQ. 5), so update T, Y from TTRIAL,
C     YTRIAL, add 1 to the no of successful steps, and set the no of
C     successive failures to zero
      T = C(17)
      DO 610 K = 1, N
         Y(K) = W(K,12)
 610  CONTINUE
      C(22) = C(22) + 1.D0
      C(23) = 0.D0

C**************RETURN(WITH IND=3, TEND SAVED, FLAG SET) IF T .EQ. TEND
      IF (T .NE. TEND) GO TO 620
      IND = 3
      C(20) = TEND
      C(21) = 1.D0
      RETURN
 620  CONTINUE
      GO TO 630
 600  CONTINUE

C ... Step not accepted (IND .EQ. 6), so add 1 to the no of
C     successive failures
      C(23) = C(23) + 1.D0

C**************ERROR RETURN (WITH IND=-3) IF HMAG .LE. HMIN
      IF (C(14) .GT. C(13)) GO TO 640
      IND = -3
      RETURN
 640  CONTINUE
 630  CONTINUE

C END STAGE 4
C***********************************************************************

      GO TO 99999
C ... End loop

C ... Begin abort action
 3000  CONTINUE

       WRITE(6,3010) IND, TOL, T, N, C(13), TEND, NW, C(16), C(20),
     +      C(22), C(23), C(24), (Y(K), K = 1, N)
 3010  FORMAT( /// 1H ,58HCOMPUTATION STOPPED IN DNORK WITH THE FOLLOWIN
     +G VALUES -
     +   / 1H , 5HIND =, I4, 5X, 6HTOL  =, 1PD13.6, 5X, 11HX         =,
     +          1PD22.15
     +   / 1H , 5HN   =, I4, 5X, 6HHMIN =, 1PD13.6, 5X, 11HTEND      =,
     +          1PD22.15
     +   / 1H , 5HNW  =, I4, 5X, 6HHMAX =, 1PD13.6, 5X, 11HPREV TEND =,
     +          1PD22.15
     +   / 1H , 14X, 27HNO OF SUCCESSFUL STEPS    =, 0PF8.0
     +   / 1H , 14X, 27HNO OF SUCCESSIVE FAILURES =, 0PF8.0
     +   / 1H , 14X, 27HNO OF FUNCTION EVALS      =,0PF8.0
     +   / 1H , 23HTHE COMPONENTS OF Y ARE
     +   // (1H , 1P5D24.15))

      STOP

C ... End abort action

      END

C *** End RDMETH ***

      SUBROUTINE STORE(N,NW,LOLD,T,Y,IFIRST,ILAST,TOLD,KOLD,W)

C*********************************************************************
C                                                                    *
C     Purpose - This routine is for storing T value, K-values, and   *
C               solutions to TOLD and KOLD queue.                    *
C               When K-values for Tn-1 to Tn and Yn have             *
C               computed, store                                      *
C               TOLD(ILAST+1) <--- Tn                                *
C               KOLD(*,1,ILAST+1) <--- Yn-1                          *
C               KOLD(*,2,ILAST+1) <--- K1                            *
C               KOLD(*,3,ILAST+1) <--- Yn                            *
C               KOLD(*,4,ILAST+1) <--- Ks+1(K8)                      *
C               KOLD(*,5,ILAST+1) <--- Ks+2(K9)                      *
C               KOLD(*,6,ILAST+1) <--- Ks+3(K10)                     *
C               KOLD(*,7,ILAST+1) <--- Ks+4(K11)                     *
C                                                                    *
C*********************************************************************
C                                                                    *
C    N      Number of equation                                       *
C    NW     First dimension of workspace W & queue KOLD              *
C    LOLD   Length of queue                                          *
C    T      Independent value (Tn)                                   *
C    Y      Solutions at Tn-1                                        *
C    IFIRST First index of queue                                     *
C    ILAST  Last index of queue                                      *
C    TOLD   Queue of independent values                              *
C    KOLD   Queue of k-values and solutions                          *
C    W      Workspace                                                *
C                                                                    *
C*********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,LOLD
      DOUBLE PRECISION T,IFIRST,ILAST

C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),TOLD(LOLD),KOLD(NW,7,LOLD),W(NW,*)

C     .. Local Scalars ..
      INTEGER I,M,IFST,ILST

C     .. Intrisic Functions ..
      INTRINSIC IDINT,DFLOAT

      IFST = IDINT(IFIRST)
      ILST = IDINT(ILAST)

C ... CASE 1:  IFST <= ILST & ILST < LOLD
      IF ((ILST .GE. IFST) .AND. (ILST .LT. LOLD)) THEN
         M = ILST + 1
         ILAST = ILAST + 1.D0

C ... CASE 2: IFST < ILST = LOLD
      ELSE IF (ILST .EQ. LOLD) THEN
         print *, 'Warning:the length of the queue is short!'
         M = 1
         IFIRST = 2.D0
         ILAST = 1.D0

C ... CASE 3: ILST < IFST < LOLD
      ELSE IF (IFST .LT. LOLD) THEN
         M = ILST + 1
         IFIRST = IFIRST + 1.D0
         ILAST = ILAST + 1.D0

C ... CASE 4: ILST < IFST = LOLD
      ELSE
         M = LOLD
         IFIRST = 1.D0
         ILAST = DFLOAT(LOLD)
      END IF

C ... Store T value
      TOLD(M) = T

C ... Store K values and soltuion
      DO 100 I = 1, N
         KOLD(I,1,M) = Y(I)
         KOLD(I,2,M) = W(I,1)
         KOLD(I,3,M) = W(I,12)
         KOLD(I,4,M) = W(I,8)
         KOLD(I,5,M) = W(I,9)
         KOLD(I,6,M) = W(I,10)
         KOLD(I,7,M) = W(I,11)
 100  CONTINUE

c      print *, 'm =', m
c      print *, 'told =', told(m)
c      print *, 'kold1 =', (kold(1,i,m), i= 1, 7)
c      print *, 'kold2 =', (kold(2,i,m), i= 1, 7)
c      print *, 'kold3 =', (kold(3,i,m), i= 1, 7)

      RETURN
      END

C *** End of STORE ***

      SUBROUTINE STORED(N,NW,LOLD,T,TS,Y,YS,TOLD,KOLD,W,C)

C*********************************************************************
C                                                                    *
C     Purpose - This routine is for storing T value, derivatives,    *
c               and solutions to TOLD and KOLD queue just after      *
c               finding discontinuity.                               *
C                                                                    *
C               TOLD(ILAST+1) <--- Ts                                *
C               KOLD(*,1,ILAST+1) <--- Yn-1                          *
C               KOLD(*,2,ILAST+1) <--- K1                            *
C               KOLD(*,3,ILAST+1) <--- Ys                            *
C               KOLD(*,4,ILAST+1) <--- Z'(Ts)                        *
C               KOLD(*,5,ILAST+1) <--- Z'(T+C9*Hs)                   *
C               KOLD(*,6,ILAST+1) <--- Z'(T+C10*Hs)                  *
C               KOLD(*,7,ILAST+1) <--- Z'(T+C11*Hs)                  *
C                                                                    *
C*********************************************************************
C                                                                    *
C    N      Number of equation                                       *
C    NW     First dimension of workspace W & queue KOLD              *
C    LOLD   Length of queue                                          *
C    T      Independent value (left mesh point)                      *
C    TS     Independent value (right mesh point)                     *
C    Y      Solution at T                                            *
C    YS     Solution at TS                                           *
C    TOLD   QUEUE OF INDEPENDENT VALUES                              *
C    KOLD   QUEUE OF K-VALUES AND SOLUTIONS                          *
C    W      WORKSPACE                                                *
C    C      COMMUNICATION ARRAY                                      *
C                                                                    *
C*********************************************************************

C     .. Scalar Arguments ..
      INTEGER N,NW,LOLD
      DOUBLE PRECISION T,TS

C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),YS(N),TOLD(LOLD),KOLD(NW,7,LOLD)
      DOUBLE PRECISION W(NW,*),C(*)

C     .. Local Scalars ..
      INTEGER I,IFST,ILST,M
      DOUBLE PRECISION  C9,C10,C11,H,HS

C     .. External Subroutines ..
      EXTERNAL DERIV

C     .. Intrisic Functions ..
      INTRINSIC IDINT,DFLOAT

C     .. Data statements ..
      DATA C9/0.5D0/
      DATA C10/0.1428571428571428D0/
      DATA C11/0.2857142857142857D0/

      H = C(18)
      HS = TS - T

      IFST = IDINT(C(26))
      ILST = IDINT(C(27))

C ... CASE 1:  IFST <= ILST & ILST < LOLD
      IF ((ILST .GE. IFST) .AND. (ILST .LT. LOLD)) THEN
         M = ILST + 1
         C(27) = C(27) + 1.D0

C ... CASE 2: IFST < ILST = LOLD
      ELSE IF (ILST .EQ. LOLD) THEN
         print *, 'Warning:the length of the queue is short!'
         M = 1
         C(26) = 2.D0
         C(27) = 1.D0

C ... CASE 3: ILST < IFST < LOLD
      ELSE IF (IFST .LT. LOLD) THEN
         M = ILST + 1
         C(26) = C(26) + 1.D0
         C(27) = C(27) + 1.D0

C ... CASE 4: ILST < IFST = LOLD
      ELSE
         M = LOLD
         C(26) = 1.D0
         C(27) = DFLOAT(LOLD)
      END IF

C ... Store T value
      TOLD(M) = TS

C ... Store Yn-1, K1, Ys, F(Ts, Ys, DEL) to KOLD
      DO 20 I = 1, N
         KOLD(I,1,M) = Y(I)
         KOLD(I,2,M) = W(I,1)
         KOLD(I,3,M) = YS(I)
 20   CONTINUE

C ... Store Z'(T+C9*Hs), Z'(T+C10*Hs), Z'(T+C11*Hs) to KOLD
      CALL DERIV(N,NW,T,Y,TS,KOLD(1,4,M),H,W)
      CALL DERIV(N,NW,T,Y,T+C9*HS,KOLD(1,5,M),H,W)
      CALL DERIV(N,NW,T,Y,T+C10*HS,KOLD(1,6,M),H,W)
      CALL DERIV(N,NW,T,Y,T+C11*HS,KOLD(1,7,M),H,W)

      RETURN
      END

C *** End of STORED ***

      SUBROUTINE TFIND(N,TOLD,T,IFIRST,ILAST,M)

C****************************************************************
C                                                               *
C     Purpose - This subroutine finds the index for the         *
C               independent valriable from the queue of the     *
C               past values by using binary search.             *
C                                                               *
C     If TOLD(M) <= T < TOLD(M+1), M is returned.               *
C                                                               *
C****************************************************************
C                                                               *
C     N       The length of queue                               *
C     TOLD    Queue of past independent variable                *
C     T       Independent variavle - wish to find the indet     *
C     IFIRST  The first index of the queue                      *
C     ILAST   The last index of the queue                       *
C     M       Index found (Output)                              *
C                                                               *
C****************************************************************

C     .. Scalar Arguments ..
      INTEGER N,M
      DOUBLE PRECISION T,IFIRST,ILAST

C     .. Array Arguments ..
      DOUBLE PRECISION TOLD(N)

C     .. Local Scalars ..
      INTEGER I,J,K,IFST,ILST,REALK,RK,RKP,RKM

C     .. Intrisic Functions ..
      INTRINSIC IDINT,DFLOAT

C ... Initialize variables
      IFST = IDINT(IFIRST)
      ILST = IDINT(ILAST)
      I = IFST

C ... Calculate the last index J for binary search
      IF (ILST .EQ. IFST) THEN
         PRINT *, '***** ERROR:TFIND ***** '
         PRINT *, 'No data is stored in the queue.:'
         Print *, 'T =', t
         PRINT *
         STOP
      ELSE IF (IFST .LT. ILST) THEN
         J = ILST
      ELSE
         J = N + ILST
      END IF

C ... Loop - Binary search
 100  CONTINUE
C ... Calculate the middel of I and J
         K = (I + J) / 2

C ... Calculate real index for K, K+1, and K-1
         RK  = REALK(K, N)
         RKP = REALK(K+1, N)
         RKM = REALK(K-1, N)

C ... CASE 1:  TOLD(RK) = T
         IF (TOLD(RK) .EQ. T) THEN
            M = RK
            RETURN

C ... CASE 2:  TOLD(RK) < T
         ELSE IF (TOLD(RK) .LT. T) THEN

C ... CASE 2.1: IF RK IS ILST THEN MUST STOP
            IF (RK .EQ. ILST) THEN
               PRINT *, '***** ERROR:TFIND ***** '
               PRINT *, 'The delay is advanced than the history'
               print *, 'T =', t
               PRINT *
               STOP

C ... CASE 2.2: TOLD(RK) < T < TOLD(RK+1) THEN M = RK
            ELSE IF (T .LT. TOLD(RKP)) THEN
               M = RK
               RETURN

C ... CASE 2.3: TOLD(RK+1) < T , CHANGE I TO K + 1
            ELSE
               I = K + 1
            END IF

C ... CASE 3: NOW T < TOLD(RK)
C ... CASE 3.1: IF K IS IFST THEN MUST STOP
         ELSE IF (K .EQ. IFST) THEN
            PRINT *, ' ***** ERROR:TFIND ***** '
            PRINT *, 'The length of the queue is short'
            PRINT *, 'T =', t, ', k =', k, ', Told(k) =', told(k)
            PRINT *
            STOP

C ... CASE 3.2: TOLD(RK-1) < T < TOLD(RK) THEM M = RK - 1
         ELSE IF (TOLD(RKM) .LT. T) THEN
            M = RKM
            RETURN

C ... CASE 3.3: T < TOLD(RK-1)
         ELSE
            J = K - 1
         END IF
      GO TO 100
C ... END LOOP - BINARY SEARCH

      END

C *** End of TFIND ***

      FUNCTION REALK(K,N)

C*****************************************************************
C                                                                *
C     Purpose - This function computes real index of the queue.  *
C                                                                *
C*****************************************************************
C                                                                *
C   K      The index                                             *
C   N      The length of queue                                   *
C   REALK  Real index (Output)                                   *
C                                                                *
C*****************************************************************

C     .. Scalar Arguments ..
      INTEGER K, N, REALK

C ... If the K is between 1 to N then K is the real index
      IF (K .LE. N) THEN
         REALK = K

C ... If the K is bigger than N then return K MOD N
      ELSE
         REALK = MOD(K, N)
      END IF

      RETURN
      END

C *** End of REALK