MODULE DDVERK90_M
!-----------------------------------------------
INTERFACE DDVERK90
	MODULE PROCEDURE DDVERK_1,DDVERK_2,DDVERK_3,DDVERK_4
END INTERFACE
!-----------------------------------------------
TYPE, PUBLIC :: DDVERK_SOL
    INTEGER :: NPTS,FLAG,NE
    INTEGER, DIMENSION(6) :: STATS
    DOUBLE PRECISION, DIMENSION(:), POINTER :: T
    DOUBLE PRECISION, DIMENSION(:,:), POINTER :: Y
    DOUBLE PRECISION, DIMENSION(:), POINTER :: TE
    DOUBLE PRECISION, DIMENSION(:,:), POINTER :: YE
    DOUBLE PRECISION, DIMENSION(:,:,:), POINTER :: QUEUE
    DOUBLE PRECISION, DIMENSION(:), POINTER :: TQUEUE
    INTEGER, DIMENSION(:), POINTER :: IE
    INTEGER, DIMENSION(7):: IPOINT
END TYPE DDVERK_SOL

TYPE, PUBLIC :: DDVERK_OPTS
    LOGICAL, DIMENSION(:), POINTER :: ISTERMINAL
    INTEGER, DIMENSION(:), POINTER :: DIRECTION
    DOUBLE PRECISION :: HINIT,HMAX,RELERR
    DOUBLE PRECISION, DIMENSION(:), POINTER :: ABSERR
    LOGICAL :: NEUTRAL,TRACK_DISCONTINUITIES,INTERPOLATION
    INTEGER :: MAX_EVENTS,MAX_STEPS
    DOUBLE PRECISION, DIMENSION(:), POINTER :: THIT_EXACTLY
	INTEGER :: DEFECT_ESTIMATE_TYPE
	DOUBLE PRECISION :: SCALE
END TYPE DDVERK_OPTS

TYPE, PUBLIC :: DDVERK_INT
    DOUBLE PRECISION, DIMENSION(:), POINTER :: TVALS
    DOUBLE PRECISION, DIMENSION(:,:), POINTER :: YT,DT
    INTEGER, DIMENSION(:), POINTER :: COMPONENTS
END TYPE DDVERK_INT

!-----------------------------------------------
PUBLIC :: DDVERK_SET,DDVERK_PRINT_STATS,DDVERK_VAL, &
			DDVERK_USER,DDVERK_IQUEUE
!-----------------------------------------------
! Set machine characteristic variables using intrinsics.
!							 C(10)			
DOUBLE PRECISION ,PRIVATE :: G_RREB=2.D0*EPSILON(1D0)
DOUBLE PRECISION, PRIVATE :: G_REALMAX=HUGE(1D0),G_UROUND=EPSILON(1D0)
!-----------------------------------------------
!WORK ARRAYS
LOGICAL, PRIVATE :: G_IS_NEUTRAL,G_IS_CONSTANT_HISTORY,G_IS_CONSTANT_DELAYS, &
					G_TRACK_DISCONTINUITIES

LOGICAL, PRIVATE :: G_HAVE_EVENT_FCN,G_HAVE_OUT_FCN,G_HAVE_CHANGE_FCN

DOUBLE PRECISION, DIMENSION(:), POINTER, PRIVATE :: &
					G_PASS_DELAYS,G_PASS_HISTORY


TYPE(DDVERK_SOL), POINTER, PRIVATE :: G_PASS_SOL

INTEGER, PRIVATE :: G_NSPAN,G_NEQN,G_NLAGS,G_NW,G_NW2,G_CL,G_NGUSER,G_NFAIL_STEPS, &
				G_K
DOUBLE PRECISION ,PRIVATE :: G_C6_ONLY_NEXT

DOUBLE PRECISION ,DIMENSION(:),ALLOCATABLE ,PRIVATE:: G_W15

INTEGER,PARAMETER,PRIVATE :: P_INITIAL_G_LOLD_LENGTH=1000
!P_MAX_G_LOLD_LENGTH=P_INITIAL_G_LOLD_LENGTH*2^k
INTEGER,PARAMETER,PRIVATE :: P_MAX_G_LOLD_LENGTH=512000
DOUBLE PRECISION,PARAMETER,PRIVATE :: &
			P_DEFAULT_RELERR=1D-3,P_DEFAULT_ABSERR=1D-6 
	

PRIVATE :: DUMMY_GUSER,DUMMY_CHANGE,SOL_OUT
PRIVATE :: CONSTANT_BETA ,CONSTANT_HISTORY

DOUBLE PRECISION ,DIMENSION(:),ALLOCATABLE, PRIVATE::G_C,G_TOLD
DOUBLE PRECISION ,DIMENSION(:,:,:),ALLOCATABLE, PRIVATE :: G_KOLD
INTEGER ,PRIVATE :: G_LOLD
DOUBLE PRECISION,POINTER,PRIVATE ::G_T
DOUBLE PRECISION,DIMENSION(:),POINTER,PRIVATE ::G_Y
DOUBLE PRECISION,DIMENSION(:,:),POINTER,PRIVATE ::G_W
LOGICAL,PRIVATE::G_LEGAL_W
!-----------------------------------------------
CONTAINS

FUNCTION HZ_ISNAN(X) RESULT(L)
 DOUBLE PRECISION X,Y
 LOGICAL L
 Y=X+(1.D0-2.D0+1.D0)
 !IF(abs(Y-X) < 1.D0) THEN
 IF(Y==X) THEN
   L=.FALSE.
 ELSE
   L=.TRUE.
   !PRINT *,X,Y
 END IF
 !L=ISNAN(X)
END FUNCTION HZ_ISNAN

FUNCTION DDVERK_1(NVAR,NLAGS,DDES,BETA,HISTORY,TSPAN,OPTIONS,EVENT_FCN,     &
                 CHANGE_FCN,OUT_FCN) RESULT(SOL)

	IMPLICIT NONE
    INTEGER, DIMENSION(:) :: NVAR
    INTEGER, DIMENSION(:) :: NLAGS
    TYPE(DDVERK_OPTS), OPTIONAL :: OPTIONS
    TYPE(DDVERK_SOL) :: SOL
    DOUBLE PRECISION, DIMENSION(:) :: TSPAN
    EXTERNAL DDES,BETA,HISTORY ,EVENT_FCN,CHANGE_FCN,OUT_FCN
	OPTIONAL EVENT_FCN,CHANGE_FCN,OUT_FCN ,HISTORY

	PRINT *,"DDVERK_1 CALLED"
    
    NULLIFY(G_PASS_DELAYS,G_PASS_HISTORY)
    G_IS_CONSTANT_DELAYS = .FALSE.
    G_IS_CONSTANT_HISTORY = .FALSE.

	SOL=APPLY_DDVERK1(NVAR,NLAGS,DDES,BETA,HISTORY,TSPAN,OPTIONS,EVENT_FCN,	&
				CHANGE_FCN,OUT_FCN)
END FUNCTION DDVERK_1
!-----------------------------------------------
!-----------------------------------------------
FUNCTION DDVERK_2(NVAR,NLAGS,DDES,DELAYS,HISTORY,TSPAN,OPTIONS,EVENT_FCN,     &
                 CHANGE_FCN,OUT_FCN) RESULT(SOL)

    IMPLICIT NONE
    INTEGER, DIMENSION(:) :: NVAR
    INTEGER, DIMENSION(:) :: NLAGS
    DOUBLE PRECISION, DIMENSION(:), TARGET :: DELAYS
    TYPE(DDVERK_OPTS), OPTIONAL :: OPTIONS
    TYPE(DDVERK_SOL) :: SOL
    DOUBLE PRECISION, DIMENSION(:) :: TSPAN
    EXTERNAL DDES,HISTORY,EVENT_FCN,CHANGE_FCN,OUT_FCN
	OPTIONAL EVENT_FCN,CHANGE_FCN,OUT_FCN

	PRINT *,"DDVERK_2 CALLED"

    NULLIFY(G_PASS_DELAYS,G_PASS_HISTORY)
    G_PASS_DELAYS => DELAYS
    G_IS_CONSTANT_DELAYS = .TRUE.
    G_IS_CONSTANT_HISTORY = .FALSE.
	
	SOL=APPLY_DDVERK1(NVAR,NLAGS,DDES,CONSTANT_BETA,HISTORY,TSPAN,OPTIONS,EVENT_FCN,	&
				CHANGE_FCN,OUT_FCN)
END FUNCTION DDVERK_2
!-----------------------------------------------
!-----------------------------------------------
FUNCTION DDVERK_3(NVAR,NLAGS,DDES,BETA,HISTORY,TSPAN,OPTIONS,EVENT_FCN,     &
                 CHANGE_FCN,OUT_FCN) RESULT(SOL)

    IMPLICIT NONE
    INTEGER, DIMENSION(:) :: NVAR
    INTEGER, DIMENSION(:) :: NLAGS
    DOUBLE PRECISION, DIMENSION(:), TARGET :: HISTORY
    TYPE(DDVERK_OPTS), OPTIONAL :: OPTIONS
    TYPE(DDVERK_SOL) :: SOL
    DOUBLE PRECISION, DIMENSION(:) :: TSPAN
    EXTERNAL DDES,BETA ,EVENT_FCN,CHANGE_FCN,OUT_FCN
    OPTIONAL EVENT_FCN,CHANGE_FCN,OUT_FCN 

	PRINT *,"DDVERK_3 CALLED"

    NULLIFY(G_PASS_DELAYS,G_PASS_HISTORY)
    G_IS_CONSTANT_DELAYS = .FALSE.
    G_PASS_HISTORY => HISTORY
    G_IS_CONSTANT_HISTORY = .TRUE.

	SOL=APPLY_DDVERK1(NVAR,NLAGS,DDES,BETA,CONSTANT_HISTORY,TSPAN,OPTIONS,EVENT_FCN,	&
				CHANGE_FCN,OUT_FCN)
END FUNCTION DDVERK_3
!-----------------------------------------------
!-----------------------------------------------
FUNCTION DDVERK_4(NVAR,NLAGS,DDES,DELAYS,HISTORY,TSPAN,OPTIONS,EVENT_FCN,     &
                 CHANGE_FCN,OUT_FCN) RESULT(SOL)

    IMPLICIT NONE
    INTEGER, DIMENSION(:) :: NVAR
    INTEGER, DIMENSION(:) :: NLAGS
    DOUBLE PRECISION, DIMENSION(:), TARGET :: DELAYS,HISTORY
    TYPE(DDVERK_OPTS), OPTIONAL :: OPTIONS
    TYPE(DDVERK_SOL) :: SOL
    DOUBLE PRECISION, DIMENSION(:) :: TSPAN
    EXTERNAL DDES ,EVENT_FCN,CHANGE_FCN,OUT_FCN
    OPTIONAL EVENT_FCN,CHANGE_FCN,OUT_FCN 

	PRINT *,"DDVERK_4 CALLED"

    NULLIFY(G_PASS_DELAYS,G_PASS_HISTORY)
    G_PASS_DELAYS => DELAYS
    G_IS_CONSTANT_DELAYS = .TRUE.
    G_PASS_HISTORY => HISTORY
    G_IS_CONSTANT_HISTORY = .TRUE.

	SOL=APPLY_DDVERK1(NVAR,NLAGS,DDES,CONSTANT_BETA,CONSTANT_HISTORY,TSPAN,OPTIONS,EVENT_FCN,	&
				CHANGE_FCN,OUT_FCN)
END FUNCTION DDVERK_4

!-----------------------------------------------
!-----------------------------------------------
FUNCTION APPLY_DDVERK1(NVAR,NLAGS,DDES,DELAYS,HISTORY,TSPAN,OPTIONS,EVENT_FCN,     &
                 CHANGE_FCN,OUT_FCN) RESULT(SOL)
	IMPLICIT NONE
    INTEGER :: IFLAG,IER,NOUT
    INTEGER, DIMENSION(:) :: NVAR
    INTEGER, DIMENSION(:) :: NLAGS
    TYPE(DDVERK_OPTS), OPTIONAL :: OPTIONS
    TYPE(DDVERK_OPTS), TARGET :: OPTS
    TYPE(DDVERK_SOL), TARGET :: SOL
    DOUBLE PRECISION, DIMENSION(:) :: TSPAN
    !DOUBLE PRECISION, DIMENSION(2*NVAR(1)) :: YOLD
    EXTERNAL DDES,DELAYS,HISTORY ,EVENT_FCN,CHANGE_FCN,OUT_FCN
	OPTIONAL EVENT_FCN,CHANGE_FCN,OUT_FCN 

	INTEGER ::IFIRST,ILAST,QUEUE_SIZE,QUEUE_INDEX,I
!must be implemented
	!********************************
	!*******
	INTEGER MU,OMEGA
	!***
	PRINT *,"APPLY_DDVERK1 CALLED"

	G_NEQN = NVAR(1)
	G_HAVE_EVENT_FCN = PRESENT(EVENT_FCN)

    IF (G_HAVE_EVENT_FCN) THEN
       IF (SIZE(NVAR) < 2) THEN
          PRINT *, "Must input the number of event functions, NVAR(2)."
          STOP
       ELSE
          G_NGUSER = NVAR(2)
       END IF
       ALLOCATE(SOL%IE(10),SOL%TE(10),SOL%YE(10,G_NEQN),STAT=IER)
       CALL CHECK_STAT(IER,1)
    ELSE
       G_NGUSER = 0
    END IF
    SOL%NE = 0
	

	IF (PRESENT(OPTIONS)) THEN
       OPTS = OPTIONS
    ELSE
       OPTS = DDVERK_SET()
    END IF

    CALL EXPAND_OPTS(G_NEQN,G_NGUSER,OPTS)

    G_HAVE_OUT_FCN = PRESENT(OUT_FCN)

    G_NSPAN = SIZE(TSPAN)
    IF (G_HAVE_OUT_FCN .AND. G_NSPAN > 2) THEN
       ! Change to output at every step.
       TSPAN = (/ TSPAN(1), TSPAN(G_NSPAN) /)
       G_NSPAN = 2
    END IF
    CALL CHECK_TSPAN(TSPAN,NOUT)
    ALLOCATE(SOL%T(NOUT),SOL%Y(NOUT,G_NEQN),STAT=IER)
    CALL CHECK_STAT(IER,2)
    SOL%NPTS = 0
    G_PASS_SOL => SOL

	G_HAVE_CHANGE_FCN=PRESENT(CHANGE_FCN)

	G_IS_NEUTRAL = OPTS%NEUTRAL
	!***********
	MU=NLAGS(1)
	IF(G_IS_NEUTRAL) THEN
		IF (SIZE(NLAGS) < 2) THEN
			PRINT *,' Must input the number of derivative lags, NLAGS(2).'
			STOP
		END IF
		OMEGA=NLAGS(2)
	ELSE
		OMEGA=0
	END IF
	G_NLAGS=MU+OMEGA

	!***********
	IF(G_HAVE_CHANGE_FCN) THEN

       IF (G_HAVE_OUT_FCN) THEN !CHANGE_FCN:yes EVENT_FCN:yes OUT_FCN:yes
			CALL APPLY_DDVERK2(MU,OMEGA,TSPAN,OPTS,DDES,DELAYS,HISTORY,	&
					EVENT_FCN,CHANGE_FCN,OUT_FCN,IFLAG)
       ELSE !CHANGE_FCN:yes EVENT_FCN:yes OUT_FCN:no
			CALL APPLY_DDVERK2(MU,OMEGA,TSPAN,OPTS,DDES,DELAYS,HISTORY,	&
					EVENT_FCN,CHANGE_FCN,SOL_OUT,IFLAG)
       END IF

    ELSEIF (G_HAVE_EVENT_FCN) THEN
  
       IF (G_HAVE_OUT_FCN) THEN !CHANGE_FCN:no EVENT_FCN:yes OUT_FCN:yes
			CALL APPLY_DDVERK2(MU,OMEGA,TSPAN,OPTS,DDES,DELAYS,HISTORY,	&
					EVENT_FCN,DUMMY_CHANGE,OUT_FCN,IFLAG)
       ELSE !CHANGE_FCN:no EVENT_FCN:yes OUT_FCN:no
			CALL APPLY_DDVERK2(MU,OMEGA,TSPAN,OPTS,DDES,DELAYS,HISTORY,	&
					EVENT_FCN,DUMMY_CHANGE,SOL_OUT,IFLAG)
       END IF

    ELSE

       IF (G_HAVE_OUT_FCN) THEN !CHANGE_FCN:no EVENT_FCN:no  OUT_FCN:yes
			CALL APPLY_DDVERK2(MU,OMEGA,TSPAN,OPTS,DDES,DELAYS,HISTORY,	&
					DUMMY_GUSER,DUMMY_CHANGE,OUT_FCN,IFLAG)
	   ELSE !CHANGE_FCN:no EVENT_FCN:no OUT_FCN:no
			PRINT *,"NONE"
			CALL APPLY_DDVERK2(MU,OMEGA,TSPAN,OPTS,DDES,DELAYS,HISTORY,	&
					DUMMY_GUSER,DUMMY_CHANGE,SOL_OUT,IFLAG)
       END IF
	END IF

!   SAVE THE NECESSARY INTERPOLATION STRUCTURE IF NECESSARY.
	SOL%IPOINT=0
	SOL%IPOINT(1)=IFLAG
	SOL%IPOINT(3)=G_NEQN
	NULLIFY(SOL%TQUEUE,SOL%QUEUE)
	IF(IFLAG==0) THEN 
		IF(OPTS%INTERPOLATION) THEN
			SOL%IPOINT(2)=1
	
			IFIRST=G_C(26)
			ILAST=G_C(27)
			IF(IFIRST<=ILAST) THEN
				QUEUE_SIZE=ILAST-IFIRST+1
			ELSE
				QUEUE_SIZE=G_LOLD-(IFIRST-ILAST-1)
			END IF
			ALLOCATE(SOL%TQUEUE(QUEUE_SIZE))
			ALLOCATE(SOL%QUEUE(G_NEQN,7,QUEUE_SIZE))
			SOL%IPOINT(4)=QUEUE_SIZE
			QUEUE_INDEX=IFIRST
			DO I=1,QUEUE_SIZE
				SOL%TQUEUE(I)=G_TOLD(QUEUE_INDEX)
				SOL%QUEUE(:,:,I)=G_KOLD(:,:,QUEUE_INDEX)
				QUEUE_INDEX=QUEUE_INDEX+1
				IF(QUEUE_INDEX>G_LOLD) QUEUE_INDEX=1
			END DO
		ELSE
			SOL%IPOINT(2)=0
		END IF
	ELSE
       PRINT *,' *******Error in APPLY_DDVERK1******'
       PRINT *,' One or more errors occurred.'
	   PRINT *,' LAST T= ', G_T
	   PRINT *,' IFLAG= ',IFLAG
	END IF

    CALL CONTRACT_SOL(G_PASS_SOL)

    SOL%FLAG = IFLAG
    SOL%STATS(1) = G_PASS_SOL%NPTS !G_C(22)
    SOL%STATS(2) = G_NFAIL_STEPS
    SOL%STATS(3) = G_C(24)
    SOL%STATS(4) = G_LOLD+(G_NW*7*G_LOLD)+(33+G_NEQN)+G_NGUSER+ &
					G_NEQN*G_NSPAN+G_NGUSER*G_NGUSER+G_NEQN*G_NGUSER+ &
					G_NEQN*G_NLAGS+G_NEQN+(20+G_NLAGS)

	DEALLOCATE(G_TOLD,G_KOLD,G_C,G_W15)

END FUNCTION APPLY_DDVERK1
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE APPLY_DDVERK2(MU,OMEGA,TSPAN,OPTS,DDES,DELAYS,HISTORY,	&
					EVENT_FCN,CHANGE_FCN,OUT_FCN,IFLAG)
	IMPLICIT NONE 
	INTEGER MU,OMEGA,IND
	DOUBLE PRECISION ,TARGET :: T
	DOUBLE PRECISION ::TEND,TOL
	DOUBLE PRECISION ,DIMENSION(:),ALLOCATABLE ,TARGET:: Y
	DOUBLE PRECISION ,DIMENSION(:,:),ALLOCATABLE,TARGET :: W
    DOUBLE PRECISION, DIMENSION(:) :: TSPAN
    TYPE(DDVERK_OPTS), TARGET :: OPTS
    EXTERNAL DDES,DELAYS,HISTORY ,EVENT_FCN,CHANGE_FCN,OUT_FCN
    INTEGER :: IFLAG
	LOGICAL :: RECEIVED_TO_TEND
	INTENT(IN) :: MU,OMEGA,TSPAN
	INTENT(OUT) IFLAG
	
	DOUBLE PRECISION,DIMENSION(:),ALLOCATABLE::GLEFT,GLEFT_NEXT,TEVENT
	DOUBLE PRECISION,DIMENSION(:,:),ALLOCATABLE::YEVENT,DYEVENT
	INTEGER :: NOFTEVENT
	LOGICAL,DIMENSION(:,:),ALLOCATABLE::HAVE_EVENT
	LOGICAL,DIMENSION(G_NGUSER) ::NONE_EVENT
	DOUBLE PRECISION ,DIMENSION(:),ALLOCATABLE :: Y0,DY0
	DOUBLE PRECISION ,DIMENSION(:,:),ALLOCATABLE :: Z0
	DOUBLE PRECISION ,DIMENSION(:),ALLOCATABLE ::YEVENT_CHANGE,DYEVENT_CHANGE

	
	INTEGER I,J

	INTEGER NTHIT
	INTEGER C_THIT_I

	LOGICAL :: PREV_G_EXIST


	DOUBLE PRECISION :: HINIT
	LOGICAL :: QUIT,CHANGED,TERMINATE


	INTEGER :: IT_FIRST,IT_LAST,LAST_IE,TEMP_N,C_TEMP_N,C_LAST_IE

	LOGICAL :: SHOULD_RETURN
	DOUBLE PRECISION ::T_MAX
	DOUBLE PRECISION ,DIMENSION(G_NSPAN) ::TEMP_T
	DOUBLE PRECISION ,DIMENSION(G_NEQN,G_NSPAN) ::TEMP_Y

	DOUBLE PRECISION :: PREV_CRITICAL_EVENT_T


	IF(G_NGUSER>0) THEN
		ALLOCATE(GLEFT(G_NGUSER),GLEFT_NEXT(G_NGUSER),TEVENT(G_NGUSER))
		ALLOCATE(HAVE_EVENT(G_NGUSER,G_NGUSER)) 
		ALLOCATE(YEVENT(G_NEQN,G_NGUSER),DYEVENT(G_NEQN,G_NGUSER))
		ALLOCATE(Y0(G_NEQN),DY0(G_NEQN),Z0(G_NEQN,G_NLAGS))
		ALLOCATE(YEVENT_CHANGE(G_NEQN),DYEVENT_CHANGE(G_NEQN))
	END IF

	PRINT *,"APPLY_DDVERK2 CALLED"
	
	T=TSPAN(1)
	TEND=TSPAN(G_NSPAN)
	G_NW=G_NEQN
	G_NW2=20+MU+OMEGA
	G_CL=33+G_NEQN

	IT_FIRST=1
	IT_LAST=1

	G_LOLD=P_INITIAL_G_LOLD_LENGTH


	ALLOCATE(Y(G_NEQN),G_TOLD(G_LOLD),G_C(G_CL))
	ALLOCATE(W(G_NW,G_NW2))
	ALLOCATE(G_W15(MU+OMEGA))
	ALLOCATE(G_KOLD(G_NW,7,G_LOLD))
	
	G_T => T
	G_Y => Y
	G_W => W
	
    G_LEGAL_W=.FALSE.
	G_K=0
	W=0.D0

	G_C = 0.D0
	G_C(4)=OPTS%HINIT
	G_C(6)=OPTS%HMAX
	
	!G_C(1)=4
	!G_C(1)=1
	!G_C(1)=1
	!G_C(1)=6
	G_C(1)=0


	G_C(5)=OPTS%SCALE
	TOL=OPTS%RELERR
	G_C(33+1:33+G_NEQN)=OPTS%ABSERR



    G_C(32) = OPTS%DEFECT_ESTIMATE_TYPE
    IF (ASSOCIATED(OPTS%THIT_EXACTLY)) THEN
		NTHIT=SIZE(OPTS%THIT_EXACTLY)
	ELSE
		NTHIT=0
	END IF
	G_C(9) =1 !INTERUPT 2
	IND=2
	RECEIVED_TO_TEND = .FALSE.
	IFLAG=0
	PREV_G_EXIST=.FALSE.
	C_THIT_I=1
	
	SHOULD_RETURN=.FALSE.

	NONE_EVENT=.FALSE.

	PREV_CRITICAL_EVENT_T=T

	G_TRACK_DISCONTINUITIES=OPTS%TRACK_DISCONTINUITIES
	G_NFAIL_STEPS=0

	G_C6_ONLY_NEXT=0.D0

	DO WHILE(.NOT. RECEIVED_TO_TEND)
		IF((NTHIT>0) .AND. (C_THIT_I<=NTHIT)) THEN
			G_C(8) =1 !INTERUPT 1
			CALL DDVERK(G_NEQN,G_NW,MU,OMEGA,IND,T,TEND,Y,TOL,&
					  W,G_C,DDES,DELAYS,HISTORY)
			IF((OPTS%THIT_EXACTLY(C_THIT_I)>T) .AND. &
				(OPTS%THIT_EXACTLY(C_THIT_I) < T + G_C(14))) THEN
				G_C(14)=OPTS%THIT_EXACTLY(C_THIT_I)-T
			END IF
			G_C(8) =0
		END IF
		IF((IND/=3) .AND. (IND > 0)) THEN
			CALL DDVERK(G_NEQN,G_NW,MU,OMEGA,IND,T,TEND,Y,TOL, &
					  W,G_C,DDES,DELAYS,HISTORY)
		END IF
		IF(IND < 0 ) THEN
			IFLAG = IND
			EXIT
		ELSEIF((G_C(17)>TEND) .OR. (IND==3) ) THEN
			RECEIVED_TO_TEND = .TRUE.
		END IF
		IF((NTHIT>0) .AND. (C_THIT_I<=NTHIT)) THEN
			IF(ABS(OPTS%THIT_EXACTLY(C_THIT_I)-G_C(17))<=G_UROUND*ABS(G_C(17))) THEN
				C_THIT_I=C_THIT_I+1
			END IF
		END IF
		! Search for Events
		IF((IND==5) .OR. (IND==7)) THEN !accepted step or discountinuity
			LAST_IE=0
			IF(G_NGUSER>0) THEN
				!GLEFT=GLEFT_NEXT
				IF(.NOT. PREV_G_EXIST) THEN
					CALL COMPUTE_YDYZ(T,Y,T,Y0,DY0,Z0,G_C,W,DELAYS,HISTORY, &
							MU,OMEGA)
					CALL EVENT_FCN(T,Y0,DY0,Z0,GLEFT)
				END IF
				IF(IND==7) THEN
					PREV_G_EXIST=.FALSE.
					!TEMP_TRIGHT=0.1D0*T+0.9D0*C(17)
					!TEMP_TRIGHT=C(17)
				ELSE ! IND==5
					PREV_G_EXIST=.TRUE.
					!TEMP_TRIGHT=C(17)
				END IF
				CALL CHECK_FOR_EVENTS(T,Y,G_C(17),GLEFT,EVENT_FCN, &
					OPTS%DIRECTION,HAVE_EVENT,TEVENT,YEVENT,DYEVENT,NOFTEVENT, &
					GLEFT_NEXT,G_C,W,DELAYS,HISTORY,MU,OMEGA,PREV_CRITICAL_EVENT_T)
				GLEFT=GLEFT_NEXT

				IF(NOFTEVENT>0) THEN
					LAST_IE=NOFTEVENT
					DO I=1,NOFTEVENT
						HINIT=0.D0
						QUIT=.FALSE.
						CHANGED=.FALSE.
						YEVENT_CHANGE=YEVENT(:,I)
						DYEVENT_CHANGE=DYEVENT(:,I)
						CALL CHANGE_FCN(HAVE_EVENT(:,I),TEVENT(I),YEVENT_CHANGE, &
							   DYEVENT_CHANGE,HINIT,OPTS%DIRECTION,OPTS%ISTERMINAL, &
								QUIT,CHANGED)
						PREV_CRITICAL_EVENT_T=TEVENT(I)
						IF(CHANGED) THEN
							CALL HANDLE_EVENT_DISC(TEVENT(I),YEVENT_CHANGE, &
									DYEVENT_CHANGE,IND,T,Y,W,G_C)
							IF(HINIT>0.D0) G_C(4)=HINIT
							LAST_IE=I
							!PREV_CRITICAL_EVENT_T=TEVENT(I)
							PREV_G_EXIST=.FALSE.
							EXIT
						ELSE
							TERMINATE=.FALSE.
							DO J=1,G_NGUSER
								IF(HAVE_EVENT(J,I) .AND. OPTS%ISTERMINAL(J)) THEN
									TERMINATE=.TRUE.
									EXIT
								END IF
							END DO
							IF(TERMINATE .OR. QUIT) THEN
								LAST_IE=I
								SHOULD_RETURN=.TRUE.
								EXIT
							ELSE
							END IF
						END IF
					END DO

				END IF
			END IF
			IF(NOFTEVENT>0) THEN
				T_MAX=TEVENT(LAST_IE)
			ELSE
				T_MAX=G_C(17)
			END IF
			IF(G_NSPAN>2) THEN
				IF(TSPAN(IT_FIRST)<T_MAX) THEN
					IT_LAST=IT_FIRST
					DO WHILE((TSPAN(IT_LAST+1) < T_MAX) .AND. (IT_LAST<G_NSPAN))
						IT_LAST=IT_LAST+1
					END DO
					TEMP_N=IT_LAST-IT_FIRST+1
					TEMP_T(1:TEMP_N)=TSPAN(IT_FIRST:IT_LAST)
					DO I=1,TEMP_N
					    CALL INTRP(G_NEQN,G_NW,T,Y,TEMP_T(I),TEMP_Y(:,I),G_C(17)-T,W)
					END DO
					IT_FIRST=IT_LAST+1
				ELSE
					TEMP_N=0
				END IF
			ELSE !Normal Point Output
				TEMP_T(1)=T
				TEMP_Y(:,1)=Y
				TEMP_N=1
			END IF
			C_TEMP_N=1
			C_LAST_IE=1
			DO WHILE((C_TEMP_N<=TEMP_N) .OR. (C_LAST_IE<=LAST_IE))
				IF(LAST_IE >0) THEN
					IF((C_TEMP_N<=TEMP_N) .AND. ((TEMP_T(C_TEMP_N) <= TEVENT(C_LAST_IE)) .OR. &
									(C_LAST_IE>LAST_IE))) THEN
						CALL OUT_FCN(TEMP_T(C_TEMP_N),TEMP_Y(:,C_TEMP_N),.FALSE.,NONE_EVENT)
						C_TEMP_N=C_TEMP_N+1
					ELSEIF(C_LAST_IE <=LAST_IE) THEN
						!PRINT *,'YES TEVENT = ',TEVENT(C_LAST_IE)
						CALL OUT_FCN(TEVENT(C_LAST_IE),YEVENT(:,C_LAST_IE),.TRUE., &
										HAVE_EVENT(:,C_LAST_IE))
						C_LAST_IE=C_LAST_IE+1
					END IF
				ELSE
					CALL OUT_FCN(TEMP_T(C_TEMP_N),TEMP_Y(:,C_TEMP_N),.FALSE.,NONE_EVENT)
					C_TEMP_N=C_TEMP_N+1
				END IF
			END DO
			IF(SHOULD_RETURN) EXIT
			IF(T<10) THEN
				!PRINT *, "T1 = ",T !FMT='( "T1=" E25.12)') T
				!PRINT *, "                    Y(1)= ",Y(1)
				!PRINT *, "IND =  ",IND
				!PRINT *, "T1 =  ",T
				!PRINT *, " C17= ",C(17)
				!PRINT *, "C33= ",C(33)
			END IF
			IF((G_PASS_SOL%NE >0).AND. (LAST_IE>0)) THEN
				!PRINT *,G_PASS_SOL%TE(G_PASS_SOL%NE)
			END IF
		ELSEIF(IND==3) THEN
			CALL OUT_FCN(T,Y,.FALSE.,NONE_EVENT)
		END IF
	END DO
	!PRINT *, "IND = ",IND
	!WRITE(*,FMT='( "C(17)=" E15.8)') C(17)
!	WRITE(*,FMT='( "C(18)=" E15.8)') C(18)
	!PRINT *, "SUCCESSFUL = ",C(22)
	!PRINT *, "C(24) = ",C(24)
	!PRINT *, "T = ",T
	!PRINT *, "Y = ",Y(1)
	!PRINT *, "G_NGUSER = ",G_NGUSER
	RETURN
END SUBROUTINE APPLY_DDVERK2
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE HANDLE_EVENT_DISC(TEVENT,YEVENT,DYEVENT,IND,T,Y, &
					W,C)
	IMPLICIT NONE
	DOUBLE PRECISION ::TEVENT,T
	DOUBLE PRECISION ,DIMENSION(G_NEQN) ::YEVENT,DYEVENT,Y
	INTEGER ::IND
    DOUBLE PRECISION  C(G_CL)
    DOUBLE PRECISION,DIMENSION(G_NW,G_NW2) :: W


	PRINT *, "HANDLE_EVENT_DISC CALLED"
	PRINT *, 'TEVENT= ',TEVENT
	!PRINT *, 'YEVENT= ',YEVENT(1) ,' ',YEVENT(2)
	!************
	!UNSTORE
	C(27)=C(27)-1
	IF(C(27)==0) C(27)=G_LOLD
	!END UNSTORE
	C(29)=1.D0
	IND=7
	C(28)=0.D0
	G_NFAIL_STEPS=G_NFAIL_STEPS+C(23)
	C(23)=0.D0
	CALL STORED(G_NEQN,G_NW,T,TEVENT,Y,YEVENT,W,C)
	C(17)=TEVENT
	C(18)=TEVENT-T
	W(:,12)=YEVENT
	!************
END SUBROUTINE HANDLE_EVENT_DISC
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE CHECK_FOR_EVENTS(TLEFT,YLEFT,TRIGHT,GLEFT,EVENT_FCN,	&
					REQ_DIRECTION,HAVE_EVENT,TEVENT,YEVENT,DYEVENT,	&
					NOFTEVENT,GLEFT_NEXT,C,W,DELAYS,HISTORY,MU,OMEGA, &
					PREV_CRITICAL_EVENT_T)
	IMPLICIT NONE
	DOUBLE PRECISION :: TLEFT,TRIGHT
	DOUBLE PRECISION ,DIMENSION(G_NEQN) ::YLEFT,DYLEFT
	DOUBLE PRECISION ,DIMENSION(G_NGUSER) :: GLEFT,TEVENT

	DOUBLE PRECISION ,DIMENSION(G_NGUSER),TARGET :: GLEFT_NEXT
	EXTERNAL  EVENT_FCN,DELAYS,HISTORY
	INTEGER ,DIMENSION(G_NGUSER) :: REQ_DIRECTION !, DIRECTION
	INTEGER :: NOFTEVENT
	LOGICAL ,DIMENSION(G_NGUSER,G_NGUSER) :: HAVE_EVENT

	DOUBLE PRECISION,DIMENSION(G_NEQN,G_NGUSER)::YEVENT,DYEVENT

    DOUBLE PRECISION  C(G_CL)
    DOUBLE PRECISION,DIMENSION(G_NW,G_NW2) :: W

	INTEGER MU,OMEGA
	DOUBLE PRECISION  :: PREV_CRITICAL_EVENT_T

	INTENT(IN) :: TLEFT,TRIGHT,GLEFT,REQ_DIRECTION,YLEFT,PREV_CRITICAL_EVENT_T
	
	INTENT(OUT) :: HAVE_EVENT,GLEFT_NEXT,YEVENT ,DYEVENT,NOFTEVENT
	!**************
	DOUBLE PRECISION ,DIMENSION(:),POINTER :: GRIGHT
	DOUBLE PRECISION ,DIMENSION(G_NGUSER) :: ALL_TLEFT,ALL_TRIGHT
	DOUBLE PRECISION ,DIMENSION(G_NGUSER) :: ALL_GLEFT,ALL_GRIGHT
	INTEGER :: I,J
	LOGICAL :: ALL_FINISHED_FIND
	LOGICAL ,DIMENSION(G_NGUSER) :: GHAVE_EVENT ,DONE
	DOUBLE PRECISION :: TMID
	DOUBLE PRECISION,DIMENSION(G_NGUSER) :: GMID
	INTEGER ,DIMENSION(G_NGUSER) :: INDICES
	DOUBLE PRECISION,DIMENSION(G_NEQN,G_NGUSER)::ALL_YTLEFT,ALL_DYTLEFT


	DOUBLE PRECISION ,DIMENSION(G_NEQN) ::YMID,DYMID
	DOUBLE PRECISION ,DIMENSION(G_NEQN) ::YRIGHT,DYRIGHT,YLEFT_TEMP,DYLEFT_TEMP
	DOUBLE PRECISION,DIMENSION(G_NEQN,G_NLAGS)::ZRIGHT,ZMID,ZLEFT_TEMP
	DOUBLE PRECISION :: REL_UROUND

	GRIGHT => GLEFT_NEXT
	!COMPUTE YRIGHT DYRIGHT ZRIGHT
	CALL COMPUTE_YDYZ(TLEFT,YLEFT,TRIGHT,YRIGHT,DYRIGHT,ZRIGHT,C,W,DELAYS, &
			HISTORY,MU,OMEGA)
	!COMPUTE GRIGHT
	CALL EVENT_FCN(TRIGHT,YRIGHT,DYRIGHT,ZRIGHT,GRIGHT)

	ALL_FINISHED_FIND=.TRUE.
	NOFTEVENT=0


	GHAVE_EVENT=.FALSE.
	DO I=1,G_NGUSER
		IF ( (GLEFT(I)*GRIGHT(I) <= 0.D0)) THEN ! .AND. &
				!((GLEFT(I)/=0.D0) .OR.(TLEFT>PREV_CRITICAL_EVENT_T))) THEN
			IF(((REQ_DIRECTION(I)==0).AND.(GLEFT(I)/=GRIGHT(I))) .OR. &
				((REQ_DIRECTION(I)==1) .AND. (GLEFT(I)<GRIGHT(I))) .OR. &
				 ((REQ_DIRECTION(I)==-1) .AND. (GLEFT(I)>GRIGHT(I))) &
					) THEN
				GHAVE_EVENT(I)=.TRUE.

				ALL_TLEFT(I)=TLEFT
				ALL_TRIGHT(I)=TRIGHT		

				ALL_GLEFT(I)=GLEFT(I)
				ALL_GRIGHT(I)=GRIGHT(I)

				ALL_YTLEFT(:,I)=YLEFT
 
 				ALL_FINISHED_FIND=.FALSE.
				NOFTEVENT=1
			END IF
		END IF
	END DO
	IF(NOFTEVENT>0) THEN
		CALL COMPUTE_YDYZ(TLEFT,YLEFT,TLEFT,YLEFT_TEMP,DYLEFT_TEMP, &
				ZLEFT_TEMP,C,W,DELAYS,HISTORY,MU,OMEGA)
		DO I=1,G_NGUSER
			IF(GHAVE_EVENT(I)) THEN
				ALL_DYTLEFT(:,I)=DYLEFT_TEMP
			END IF
		END DO
	END IF

	
	IF(ALL_FINISHED_FIND) THEN
		CALL COMPUTE_YDYZ(TLEFT,YLEFT,TLEFT,ALL_YTLEFT,ALL_DYTLEFT,ZMID,C,W,DELAYS, &
						HISTORY,MU,OMEGA)
	END IF

	REL_UROUND=1.D2*G_UROUND
	DO WHILE (.NOT. ALL_FINISHED_FIND)
		DONE=.FALSE.
		ALL_FINISHED_FIND=.TRUE.
		DO I=1,G_NGUSER
			IF (GHAVE_EVENT(I) .AND. (.NOT. DONE(I))) THEN
				TMID=(ALL_TLEFT(I)+ALL_TRIGHT(I))/2
				!YMID DYMID ZMID
				CALL COMPUTE_YDYZ(TLEFT,YLEFT,TMID,YMID,DYMID,ZMID,C,W,DELAYS, &
								HISTORY,MU,OMEGA)
				CALL EVENT_FCN(TMID,YMID,DYMID,ZMID,GMID)
				DO J=I,G_NGUSER
					IF(GHAVE_EVENT(J) .AND. (.NOT. DONE(J)) .AND. &
						(ALL_TLEFT(J)==ALL_TLEFT(I)) .AND. &
						(ALL_TRIGHT(J)==ALL_TRIGHT(I)) ) THEN
						IF(ALL_GLEFT(J) * GMID(J) <= 0.D0) THEN
							ALL_TRIGHT(J)=TMID
							ALL_GRIGHT(J)=GMID(J)
						ELSE
							ALL_TLEFT(J)=TMID
							ALL_GLEFT(J)=GMID(J)
							ALL_YTLEFT(:,J)=YMID
							ALL_DYTLEFT(:,J)=DYMID
						END IF
						DONE(J)=.TRUE.
					END IF
				END DO
				IF( ALL_TRIGHT(I)-ALL_TLEFT(I) >  &
					REL_UROUND*MAX(ABS(ALL_TRIGHT(I)),ABS(ALL_TLEFT(I))) )THEN
					ALL_FINISHED_FIND=.FALSE.
				END IF
			END IF
		END DO
	END DO
	NOFTEVENT=0
	DO I=1,G_NGUSER
		IF (GHAVE_EVENT(I)) THEN
			IF(ALL_TLEFT(I)-PREV_CRITICAL_EVENT_T <=  &
				REL_UROUND*MAX(ABS(ALL_TLEFT(I)),ABS(PREV_CRITICAL_EVENT_T))) THEN
				GHAVE_EVENT(I)=.FALSE.
			ELSE
				NOFTEVENT=1
			END IF
		END IF
	END DO
	INDICES=1
	DO I=1,G_NGUSER
		IF (GHAVE_EVENT(I)) THEN
			DO J=1,G_NGUSER
				IF (GHAVE_EVENT(J) .AND.		&
					(ALL_TLEFT(J)>ALL_TLEFT(I))) THEN
					INDICES(J)=INDICES(J)+1
					!PRINT *, "YES YES"
					IF(INDICES(J) > NOFTEVENT) THEN
						NOFTEVENT=INDICES(J)
					END IF
				END IF
			END DO
		END IF						
	END DO
	HAVE_EVENT=.FALSE.
	DO I=1,G_NGUSER
		IF (GHAVE_EVENT(I)) THEN
			TEVENT(INDICES(I))=ALL_TLEFT(I)
			HAVE_EVENT(I,INDICES(I))=.TRUE.
			YEVENT(:,INDICES(I))=ALL_YTLEFT(:,I)
			DYEVENT(:,INDICES(I))=ALL_DYTLEFT(:,I)
		END IF
	END DO
	RETURN
END SUBROUTINE CHECK_FOR_EVENTS
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE COMPUTE_YDYZ(T,Y,TOUT,YOUT,DYOUT,ZOUT,C,W,DELAYS,HISTORY,MU,OMEGA)
	IMPLICIT NONE
	DOUBLE PRECISION :: T,TOUT
	DOUBLE PRECISION ,DIMENSION(G_NEQN) :: Y,YOUT,DYOUT
	DOUBLE PRECISION ,DIMENSION(G_NEQN,G_NLAGS) :: ZOUT
    DOUBLE PRECISION  C(G_CL)
    DOUBLE PRECISION,DIMENSION(G_NW,G_NW2) :: W
	EXTERNAL DELAYS,HISTORY
	INTEGER :: MU,OMEGA

	INTENT(IN) :: T,Y
	INTENT(OUT) :: YOUT,DYOUT,ZOUT
	!***************
	DOUBLE PRECISION ,DIMENSION(G_NLAGS) :: BVAL
	INTEGER I
	INTEGER M
	INTEGER MP1
	!COMPUTE YOUT
    CALL INTRP(G_NEQN,G_NW,T,Y,TOUT,YOUT,C(17)-T,W)
	!COMPUTE DYOUT
    CALL DERIV(G_NEQN,G_NW,T,Y,TOUT,DYOUT,C(17)-T,W)
	!COMPUTE ZOUT
	!*****
	CALL DELAYS(TOUT,YOUT,BVAL)
	!
	DO I=1,MU
		IF(BVAL(I)>TOUT) THEN
			PRINT *, "Warning :: Advanced Delay Error in (COMPUTE_YDYZ,1)!"
			STOP !BVAL(I)=TOUT
		END IF
		IF(BVAL(I)>=T) THEN 
		    CALL INTRP(G_NEQN,G_NW,T,Y,BVAL(I),ZOUT(:,I),C(17)-T,W)
		ELSEIF(BVAL(I)>=C(25)) THEN
			CALL TFIND(BVAL(I),C(26),C(27),M)
			IF(M+1 <= G_LOLD) THEN
				MP1=M+1
			ELSE
				MP1=1
			END IF
			CALL INTRPG(G_NEQN,G_NW,MP1,G_TOLD(M),BVAL(I),ZOUT(:,I),G_TOLD(MP1)-G_TOLD(M))
		ELSE
			CALL HISTORY(BVAL(I),ZOUT(:,I),.FALSE.)
		END IF
	END DO
	IF(G_IS_NEUTRAL) THEN
		DO I=OMEGA,G_NLAGS
			IF(BVAL(I)>TOUT) THEN
				PRINT *, "Advanced Delay Error in (COMPUTE_YDYZ,2)!"
				STOP
			ELSEIF(BVAL(I)>=T) THEN 
			    CALL DERIV(G_NEQN,G_NW,T,Y,BVAL(I),ZOUT(:,I),C(17)-T,W)
			ELSEIF(BVAL(I)>=C(25)) THEN
				CALL TFIND(BVAL(I),C(26),C(27),M)
				IF(M+1 <= G_LOLD) THEN
					MP1=M+1
				ELSE
					MP1=1
				END IF
				CALL DERIVG(G_NEQN,G_NW,MP1,G_TOLD(M),BVAL(I),ZOUT(:,I),G_TOLD(MP1)-G_TOLD(M))
			ELSE
				CALL HISTORY(BVAL(I),ZOUT(:,I),.TRUE.)
			END IF
		END DO
	END IF
	!***************
	RETURN
END SUBROUTINE  COMPUTE_YDYZ
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE SOL_OUT(T,Y,IS_EVENT,HAVE_EVENT)

    IMPLICIT NONE
    INTEGER :: NEVENT,NE,NPTS
    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(G_NEQN) :: Y
	LOGICAL :: IS_EVENT
	LOGICAL,DIMENSION(G_NGUSER) :: HAVE_EVENT

	IF (IS_EVENT) THEN
		DO NEVENT=1,G_NGUSER
			IF(HAVE_EVENT(NEVENT)) THEN
				NE = G_PASS_SOL%NE + 1
				IF (NE > SIZE(G_PASS_SOL%TE)) THEN
					CALL EXPAND_SOL_EVENTS(G_PASS_SOL,10)
				END IF
				G_PASS_SOL%NE = NE
				G_PASS_SOL%TE(NE) = T
				G_PASS_SOL%YE(NE,:) = Y(1:G_NEQN)
				G_PASS_SOL%IE(NE) = NEVENT
!				IF((T > 18).AND.(T<22)) THEN
!					T=T
!				END IF
!				PRINT *,'I= ',NEVENT,' TEVENT= ' ,T,' YEVENT= ',Y(1)
			END IF
		END DO
    ELSE
       NPTS = G_PASS_SOL%NPTS + 1
       IF (NPTS > SIZE(G_PASS_SOL%T)) THEN
          CALL EXPAND_SOL(G_PASS_SOL,100)
       END IF
       G_PASS_SOL%NPTS = NPTS
       G_PASS_SOL%T(NPTS) = T
       G_PASS_SOL%Y(NPTS,:) = Y
	   IF( MOD(INT(NPTS),1000) == 0) THEN
	   !IF ((T<1.D-2)) THEN ! .AND. (T <110.000001D0)) THEN
	   !IF ((T>48.4282D0)) THEN ! .AND. (T <110.000001D0)) THEN
       !PRINT *,' T = ' ,T,' Y(1) = ',Y(1),'G_NFAIL_STEPS = ',G_NFAIL_STEPS
       !PRINT *,' T = ' ,T,' Y(1) = ',Y(1),' HTRIAL = ',G_C(18)
       !PRINT *,'T =' ,T,' Y(1)=',Y(1)!,' Y(2)=',Y(2) !,' Y(3)=',Y(3)
	   !IF(T > 92.9525D0) THEN
		PRINT *,'T =' ,T
		!PRINT *,'DEFECT =' ,G_C(19)
		!T=T
	   END IF
    END IF

    RETURN
END SUBROUTINE SOL_OUT 
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE DUMMY_CHANGE(EVENTS,TEVENT,YEVENT,DYEVENT,HINIT, &
                          DIRECTION,ISTERMINAL,QUIT,CHANGED)
    IMPLICIT NONE
    LOGICAL,DIMENSION(:) :: EVENTS
    LOGICAL, DIMENSION(:) :: ISTERMINAL
    INTEGER, DIMENSION(:) :: DIRECTION
    DOUBLE PRECISION :: TEVENT,HINIT
    DOUBLE PRECISION, DIMENSION(:) :: YEVENT,DYEVENT
    LOGICAL :: QUIT,CHANGED
	CHANGED=.FALSE.
    RETURN
END SUBROUTINE DUMMY_CHANGE
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE DUMMY_GUSER
    IMPLICIT NONE
    RETURN
END SUBROUTINE DUMMY_GUSER
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE CONTRACT_SOL(PASS_SOL)

    IMPLICIT NONE
    TYPE(DDVERK_SOL), POINTER :: PASS_SOL
    INTEGER, DIMENSION(:), ALLOCATABLE :: ITEMP
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TEMP
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: YTEMP
    INTEGER :: NPTS,NE,COLS,IER

    IF (G_HAVE_OUT_FCN) THEN
       DEALLOCATE(PASS_SOL%T,PASS_SOL%Y,STAT=IER)
       CALL CHECK_STAT(IER,36)
       IF (G_HAVE_EVENT_FCN) THEN
          DEALLOCATE(PASS_SOL%TE,PASS_SOL%YE,PASS_SOL%IE,STAT=IER)
          CALL CHECK_STAT(IER,37)
       END IF
       RETURN
    END IF

    NPTS = PASS_SOL%NPTS
    NE = PASS_SOL%NE
    COLS = SIZE(PASS_SOL%Y,2)

    ALLOCATE(TEMP(NPTS),STAT=IER)
    CALL CHECK_STAT(IER,38)
    TEMP = PASS_SOL%T(1:NPTS)
    DEALLOCATE(PASS_SOL%T,STAT=IER)
    CALL CHECK_STAT(IER,39)
    ALLOCATE(PASS_SOL%T(NPTS),STAT=IER)
    CALL CHECK_STAT(IER,40)
    PASS_SOL%T = TEMP
    DEALLOCATE(TEMP,STAT=IER)
    CALL CHECK_STAT(IER,41)

    ALLOCATE(YTEMP(NPTS,COLS),STAT=IER)
    CALL CHECK_STAT(IER,42)
    YTEMP = PASS_SOL%Y(1:NPTS,:)
    DEALLOCATE(PASS_SOL%Y,STAT=IER)
    CALL CHECK_STAT(IER,43)
    ALLOCATE(PASS_SOL%Y(NPTS,COLS),STAT=IER)
    PASS_SOL%Y = YTEMP
    CALL CHECK_STAT(IER,44)
    DEALLOCATE(YTEMP,STAT=IER)
    CALL CHECK_STAT(IER,44)

    IF (NE > 0) THEN

       ALLOCATE(ITEMP(NE),STAT=IER)
       CALL CHECK_STAT(IER,45)
       ITEMP = PASS_SOL%IE(1:NE)
       DEALLOCATE(PASS_SOL%IE,STAT=IER)
       CALL CHECK_STAT(IER,46)
       ALLOCATE(PASS_SOL%IE(NE),STAT=IER)
       CALL CHECK_STAT(IER,47)
       PASS_SOL%IE = ITEMP
       DEALLOCATE(ITEMP,STAT=IER)
       CALL CHECK_STAT(IER,48)

       ALLOCATE(TEMP(NE),STAT=IER)
       CALL CHECK_STAT(IER,49)
       TEMP = PASS_SOL%TE(1:NE)
       DEALLOCATE(PASS_SOL%TE,STAT=IER)
       CALL CHECK_STAT(IER,50)
       ALLOCATE(PASS_SOL%TE(NE),STAT=IER)
       CALL CHECK_STAT(IER,51)
       PASS_SOL%TE = TEMP
       DEALLOCATE(TEMP,STAT=IER)
       CALL CHECK_STAT(IER,52)

       ALLOCATE(YTEMP(NE,COLS),STAT=IER)
       CALL CHECK_STAT(IER,53)
       YTEMP = PASS_SOL%YE(1:NE,:)
       DEALLOCATE(PASS_SOL%YE,STAT=IER)
       CALL CHECK_STAT(IER,54)
       ALLOCATE(PASS_SOL%YE(NE,COLS),STAT=IER)
       CALL CHECK_STAT(IER,55)
       PASS_SOL%YE = YTEMP
       DEALLOCATE(YTEMP,STAT=IER)
       CALL CHECK_STAT(IER,55)

    END IF  

    RETURN
END SUBROUTINE CONTRACT_SOL 

!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE CHECK_TSPAN(TSPAN,NOUT)
    IMPLICIT NONE
    DOUBLE PRECISION, DIMENSION(:) :: TSPAN
    INTEGER :: NOUT
    IF (G_NSPAN < 2) THEN
       PRINT *,' TSPAN must have at least two entries.'
       STOP
    ELSEIF ( ABS(TSPAN(G_NSPAN)-TSPAN(1)) <= 0.0D0 ) THEN
       PRINT *,' The first and last entries of TSPAN must be different.'
       STOP
    END IF
    IF (G_NSPAN == 2) THEN
       NOUT = 100
    ELSE
       NOUT = G_NSPAN
    END IF
    RETURN
END SUBROUTINE CHECK_TSPAN
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE EXPAND_SOL(PASS_SOL,M)

    IMPLICIT NONE
    TYPE(DDVERK_SOL), POINTER :: PASS_SOL
    INTEGER, INTENT(IN) :: M
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TEMP
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: YTEMP
    INTEGER :: ROWS,COLS,IER

    ROWS = SIZE(PASS_SOL%Y,1)
    COLS = SIZE(PASS_SOL%Y,2)

    ALLOCATE(TEMP(ROWS+M),STAT=IER)
    CALL CHECK_STAT(IER,18)
    TEMP(1:ROWS) = PASS_SOL%T
    TEMP(ROWS+1:ROWS+M) = 0D0
    DEALLOCATE(PASS_SOL%T,STAT=IER)
    CALL CHECK_STAT(IER,19)
    ALLOCATE(PASS_SOL%T(ROWS+M),STAT=IER)
    CALL CHECK_STAT(IER,20)
    PASS_SOL%T = TEMP
    DEALLOCATE(TEMP,STAT=IER)
    CALL CHECK_STAT(IER,21)

    ALLOCATE(YTEMP(ROWS+M,COLS),STAT=IER)
    CALL CHECK_STAT(IER,22)
    YTEMP(1:ROWS,1:COLS) = PASS_SOL%Y
    YTEMP(ROWS+1:ROWS+M,1:COLS) = 0D0
    DEALLOCATE(PASS_SOL%Y,STAT=IER)
    CALL CHECK_STAT(IER,23)
    ALLOCATE(PASS_SOL%Y(ROWS+M,COLS),STAT=IER)
    CALL CHECK_STAT(IER,23)
    PASS_SOL%Y = YTEMP
    DEALLOCATE(YTEMP,STAT=IER)
    CALL CHECK_STAT(IER,24)

    RETURN
 END SUBROUTINE EXPAND_SOL 
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE EXPAND_SOL_EVENTS(PASS_SOL,M)

    IMPLICIT NONE
    TYPE(DDVERK_SOL), POINTER :: PASS_SOL
    INTEGER, INTENT(IN) :: M
    INTEGER, DIMENSION(:), ALLOCATABLE :: ITEMP
    DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TEMP
    DOUBLE PRECISION, DIMENSION(:,:), ALLOCATABLE :: YTEMP
    INTEGER :: ROWS,COLS,IER

    ROWS = SIZE(PASS_SOL%YE,1)
    COLS = SIZE(PASS_SOL%YE,2)

    ALLOCATE(ITEMP(ROWS+M),STAT=IER)
    CALL CHECK_STAT(IER,25)
    ITEMP(1:ROWS) = PASS_SOL%IE
    ITEMP(ROWS+1:ROWS+M) = 0
    DEALLOCATE(PASS_SOL%IE,STAT=IER)
    CALL CHECK_STAT(IER,26)
    ALLOCATE(PASS_SOL%IE(ROWS+M),STAT=IER)
    CALL CHECK_STAT(IER,27)
    PASS_SOL%IE = ITEMP
    DEALLOCATE(ITEMP,STAT=IER)
    CALL CHECK_STAT(IER,28)

    ALLOCATE(TEMP(ROWS+M),STAT=IER)
    CALL CHECK_STAT(IER,29)
    TEMP(1:ROWS) = PASS_SOL%TE
    TEMP(ROWS+1:ROWS+M) = 0D0
    DEALLOCATE(PASS_SOL%TE,STAT=IER)
    CALL CHECK_STAT(IER,30)
    ALLOCATE(PASS_SOL%TE(ROWS+M),STAT=IER)
    CALL CHECK_STAT(IER,31)
    PASS_SOL%TE = TEMP
    DEALLOCATE(TEMP,STAT=IER)
    CALL CHECK_STAT(IER,32)

    ALLOCATE(YTEMP(ROWS+M,COLS),STAT=IER)
    CALL CHECK_STAT(IER,33)
    YTEMP(1:ROWS,1:COLS) = PASS_SOL%YE
    YTEMP(ROWS+1:ROWS+M,1:COLS) = 0D0
    DEALLOCATE(PASS_SOL%YE,STAT=IER)
    CALL CHECK_STAT(IER,34)
    ALLOCATE(PASS_SOL%YE(ROWS+M,COLS),STAT=IER)
    CALL CHECK_STAT(IER,35)
    PASS_SOL%YE = YTEMP
    DEALLOCATE(YTEMP,STAT=IER)
    CALL CHECK_STAT(IER,35)

    RETURN
END SUBROUTINE EXPAND_SOL_EVENTS 
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE EXPAND_OPTS(NEQN,NGUSER,OPTS)

    IMPLICIT NONE
    INTEGER :: NEQN,NGUSER,NAE,IER
    TYPE(DDVERK_OPTS) :: OPTS
    DOUBLE PRECISION :: AE

    NAE = SIZE(OPTS%ABSERR)
    IF (NAE == 1) THEN
       AE = OPTS%ABSERR(1)
       DEALLOCATE(OPTS%ABSERR,STAT=IER)
       CALL CHECK_STAT(IER,56)
       ALLOCATE(OPTS%ABSERR(NEQN),STAT=IER)
       CALL CHECK_STAT(IER,57)
       OPTS%ABSERR = AE
    ELSE IF (NAE /= NEQN) THEN
       PRINT *,' AE_VECTOR must have NEQN components.'
       STOP
    END IF

    IF (.NOT. ASSOCIATED(OPTS%ISTERMINAL)) THEN
       ALLOCATE(OPTS%ISTERMINAL(NGUSER),STAT=IER)
       CALL CHECK_STAT(IER,58)
       OPTS%ISTERMINAL = .FALSE.
    END IF
    IF (.NOT. ASSOCIATED(OPTS%DIRECTION)) THEN
       ALLOCATE(OPTS%DIRECTION(NGUSER),STAT=IER)
       CALL CHECK_STAT(IER,59)
       OPTS%DIRECTION = 0
    END IF

    RETURN
END SUBROUTINE EXPAND_OPTS
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE CHECK_STAT(IER,CALLED_FROM)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: IER,CALLED_FROM
    IF (IER /= 0) THEN
       PRINT *,'A storage allocation error occurred.', &
       CALLED_FROM
       STOP
    END IF
    RETURN
END SUBROUTINE CHECK_STAT
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE CONSTANT_HISTORY(T,Y,IS_DERIVATIVE)
	IMPLICIT NONE
	DOUBLE PRECISION :: T
	DOUBLE PRECISION, DIMENSION(G_NEQN) :: Y
	INTENT(IN) :: T
	INTENT(OUT) :: Y
	LOGICAL,INTENT(IN) :: IS_DERIVATIVE
	IF(IS_DERIVATIVE) THEN
		Y=0.D0
	ELSE
		Y=G_PASS_HISTORY
	END IF
	RETURN
END SUBROUTINE CONSTANT_HISTORY
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE CONSTANT_BETA(T,Y,BVAL)
    
    DOUBLE PRECISION :: T
    DOUBLE PRECISION, DIMENSION(G_NEQN) :: Y
    DOUBLE PRECISION, DIMENSION(G_NLAGS) :: BVAL
	INTENT(IN) :: T,Y
	INTENT(OUT) :: BVAL
      BVAL = T-G_PASS_DELAYS

    RETURN
END SUBROUTINE CONSTANT_BETA
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
FUNCTION DDVERK_SET(RE,AE,AE_VECTOR,HINIT,HMAX,ISTERMINAL,DIRECTION,    &
  JUMPS,NEUTRAL,TRACK_DISCONTINUITIES,TRACKING_LEVEL,INTERPOLATION,    &
  THIT_EXACTLY,MAX_EVENTS,MAX_STEPS,DEFECT_ESTIMATE_TYPE,SCALE) RESULT(OPTS)

    IMPLICIT NONE
    LOGICAL, DIMENSION(:), OPTIONAL :: ISTERMINAL
    INTEGER, DIMENSION(:), OPTIONAL :: DIRECTION
    DOUBLE PRECISION, OPTIONAL :: HINIT,HMAX,RE,AE
    DOUBLE PRECISION, DIMENSION(:), OPTIONAL :: AE_VECTOR,JUMPS,       &
    THIT_EXACTLY
    LOGICAL, OPTIONAL :: NEUTRAL,TRACK_DISCONTINUITIES,INTERPOLATION
    INTEGER, OPTIONAL :: TRACKING_LEVEL,MAX_EVENTS,MAX_STEPS
    TYPE(DDVERK_OPTS) :: OPTS
    INTEGER :: IER,NAE
	INTEGER, OPTIONAL :: DEFECT_ESTIMATE_TYPE
	DOUBLE PRECISION, OPTIONAL :: SCALE

	PRINT *,"DDVERK_SET CALLED"

    IF (PRESENT(NEUTRAL)) THEN
      OPTS%NEUTRAL = NEUTRAL
    ELSE
      OPTS%NEUTRAL = .FALSE.
    END IF

    IF (PRESENT(TRACK_DISCONTINUITIES)) THEN
      OPTS%TRACK_DISCONTINUITIES = TRACK_DISCONTINUITIES
    ELSE
      OPTS%TRACK_DISCONTINUITIES = .TRUE.
    END IF

    IF (PRESENT(MAX_EVENTS)) THEN
       OPTS%MAX_EVENTS = MAX_EVENTS
    ELSE
       OPTS%MAX_EVENTS = 0
    END IF

    IF (PRESENT(MAX_STEPS)) THEN
       OPTS%MAX_STEPS = MAX_STEPS
    ELSE
       OPTS%MAX_STEPS = 0
    END IF

    IF (PRESENT(INTERPOLATION)) THEN
      OPTS%INTERPOLATION = INTERPOLATION
    ELSE
      OPTS%INTERPOLATION = .FALSE.
    END IF

    IF (PRESENT(RE)) THEN
       IF (RE <= 0D0) THEN
          PRINT *,' RE must be positive.'
          STOP
       END IF
       OPTS%RELERR = RE
    ELSE
       OPTS%RELERR = P_DEFAULT_RELERR
    END IF

    IF (PRESENT(AE_VECTOR)) THEN
       IF (MINVAL(AE_VECTOR) < 0D0) THEN
          PRINT *,' All components of AE_VECTOR must be non-negative.'
          STOP
       END IF
       NAE = SIZE(AE_VECTOR)
    ELSE
       NAE = 1
    END IF
    ALLOCATE(OPTS%ABSERR(NAE),STAT=IER)
    CALL CHECK_STAT(IER,60)
    IF (PRESENT(AE_VECTOR)) THEN
       OPTS%ABSERR = AE_VECTOR
    ELSE IF (PRESENT(AE)) THEN
       IF (AE <= 0D0) THEN
          PRINT *,' AE must be non-negative.'
          STOP
       END IF
       OPTS%ABSERR = AE
    ELSE
      OPTS%ABSERR = P_DEFAULT_ABSERR
    END IF

    IF (PRESENT(HINIT)) THEN
       OPTS%HINIT = HINIT
    ELSE
       OPTS%HINIT = 0D0
    END IF

    IF (PRESENT(HMAX)) THEN
       OPTS%HMAX = HMAX
    ELSE
       OPTS%HMAX = 0D0
    END IF

    IF (PRESENT(ISTERMINAL)) THEN
       ALLOCATE(OPTS%ISTERMINAL(SIZE(ISTERMINAL)),STAT=IER)
       CALL CHECK_STAT(IER,62)
       OPTS%ISTERMINAL = ISTERMINAL
    ELSE
       NULLIFY(OPTS%ISTERMINAL)
    END IF

    IF (PRESENT(DIRECTION)) THEN
       ALLOCATE(OPTS%DIRECTION(SIZE(DIRECTION)),STAT=IER)
       CALL CHECK_STAT(IER,63)
       OPTS%DIRECTION = DIRECTION
    ELSE
       NULLIFY(OPTS%DIRECTION)
    END IF

    IF (PRESENT(THIT_EXACTLY)) THEN
       ALLOCATE(OPTS%THIT_EXACTLY(SIZE(THIT_EXACTLY)),STAT=IER)
       CALL CHECK_STAT(IER,64)
       OPTS%THIT_EXACTLY = THIT_EXACTLY
    ELSE
       NULLIFY(OPTS%THIT_EXACTLY)
    END IF

	IF (PRESENT(DEFECT_ESTIMATE_TYPE)) THEN
		OPTS%DEFECT_ESTIMATE_TYPE=DEFECT_ESTIMATE_TYPE
	ELSE
		OPTS%DEFECT_ESTIMATE_TYPE=1
	END IF

	IF (PRESENT(SCALE)) THEN
		OPTS%SCALE=SCALE
	ELSE
		OPTS%SCALE=0
	END IF

    RETURN
END FUNCTION DDVERK_SET
!-----------------------------------------------
!-----------------------------------------------
SUBROUTINE DDVERK_PRINT_STATS(SOL)
    IMPLICIT NONE
    TYPE(DDVERK_SOL) :: SOL
	PRINT *,"DDVERK_PRINT_STATS CALLED"
!must be implemented   
    PRINT *,' '
    PRINT *,' Integration statistics: '
    PRINT *,' '
    PRINT *,' Number of successful steps                  = ',  &
      SOL%STATS(1)
    PRINT *,' Number of failed steps                      = ',  &
      SOL%STATS(2)
    PRINT *,' Number of derivative evaluations            = ',  &
      SOL%STATS(3)
    PRINT *,' Total array storage used                    = ',  &
      SOL%STATS(4)
    PRINT *,' '
END SUBROUTINE DDVERK_PRINT_STATS
!-----------------------------------------------
!-----------------------------------------------
FUNCTION DDVERK_VAL(T,SOL,COMPONENTS,DERIVATIVES) RESULT(YINT)
	IMPLICIT NONE
	TYPE(DDVERK_SOL), TARGET :: SOL
	TYPE(DDVERK_INT), TARGET :: YINT
	DOUBLE PRECISION, DIMENSION(:) :: T
	INTEGER, DIMENSION(:), OPTIONAL :: COMPONENTS
	LOGICAL, OPTIONAL :: DERIVATIVES
	!
	INTEGER ::IFLAG,NEQN,QSIZE
	LOGICAL ::INTERPOLATION,POLYD
	INTEGER :: NT,NC,I,J
	INTEGER :: IER
	INTEGER :: M,GUESS
	DOUBLE PRECISION D0,D1,D2,D3,D9,D10,D11
	DOUBLE PRECISION DD0,DD1,DD2,DD3,DD9,DD10,DD11
	DOUBLE PRECISION ::HTRIAL
	PRINT *,"DDVERK_VAL CALLED"
!   CHECK DIMENSIONS.
	IFLAG=SOL%IPOINT(1)
	INTERPOLATION=(SOL%IPOINT(2)==1)
	NEQN=SOL%IPOINT(3)
	QSIZE=SOL%IPOINT(4)

	NULLIFY(YINT%TVALS)
	NULLIFY(YINT%YT,YINT%DT)
	NULLIFY(YINT%COMPONENTS)
	IF(IFLAG /= 0) THEN
		!
		STOP
	END IF
	IF(.NOT. INTERPOLATION) THEN
		!
		PRINT *,' To use DDVERK_VAL the problem should be solved with INTERPOLATION'
		STOP
	END IF

    NT = SIZE(T)
    IF (NT < 1) THEN
		PRINT *,' Error in DDVERK_VAL.'
        PRINT *,' The length of T must be positive.'
        STOP
    END IF

!    PROCESS THE COMPONENTS ARRAY, YINT%COMPONENTS.
	IF (PRESENT(COMPONENTS)) THEN
		NC = SIZE(COMPONENTS)
        IF (NC < 1 .OR. NC > NEQN) THEN
			PRINT *,' Error in DDVERK_VAL.'
            PRINT *,'The size of COMPONENTS must be between 1 AND NEQN.'
            STOP
        ELSE
			DO I = 1, NC
				IF (COMPONENTS(I) < 1 .OR. COMPONENTS(I) > NEQN) THEN
					PRINT *,' Error in DDVERK_VAL.'
					PRINT *,' COMPONENTS(I) must be between 1 and NEQN'
					STOP
				END IF
			END DO
        END IF
        ALLOCATE(YINT%COMPONENTS(NC),STAT=IER)
        CALL CHECK_STAT(IER,65)
        YINT%COMPONENTS = COMPONENTS
	ELSE
		NC = NEQN
        ALLOCATE(YINT%COMPONENTS(NC),STAT=IER)
        CALL CHECK_STAT(IER,66)
        DO I = 1, NC
			YINT%COMPONENTS(I) = I
        END DO
	END IF

!   SET THE LOCAL CALCULATE DERIVATIVES FLAG.
    IF (PRESENT(DERIVATIVES)) THEN
		POLYD = DERIVATIVES
    ELSE
		POLYD = .FALSE.
    END IF

!   PROCESS T ARRAY, YINT%TVALS.
    ALLOCATE(YINT%TVALS(NT),STAT=IER)
    CALL CHECK_STAT(IER,67)
    YINT%TVALS(1:NT) = T(1:NT)

!   PROCESS SOLUTION ARRAY, YINT%YT.
    ALLOCATE(YINT%YT(NT,NC),STAT=IER)
    CALL CHECK_STAT(IER,68)
    !POLYS = .TRUE.

!   PROCESS DERIVATIVE ARRAY, YINT%DT.
    IF (POLYD) THEN
		ALLOCATE(YINT%DT(NT,NC),STAT=IER)
        CALL CHECK_STAT(IER,69)
    END IF

	GUESS=1
	DO I=1,NT
		M=FIND_LEFT_T_INDEX(QSIZE,SOL%TQUEUE,GUESS,T(I))
		IF((M>=1).AND.(M<=QSIZE-1)) THEN
			HTRIAL=SOL%TQUEUE(M+1)-SOL%TQUEUE(M)
			CALL COMPUTE_INTRP_COEF(SOL%TQUEUE(M),T(I),	&
					HTRIAL,D0,D1,D2,D3,D9,D10,D11)
			DO J=1,NC
				CALL APLLY_INTRP_FORMULA(HTRIAL,SOL%QUEUE(YINT%COMPONENTS(J),:,M+1), &
					D0,D1,D2,D3,D9,D10,D11,YINT%YT(I,J))
			END DO
			IF(POLYD) THEN
				CALL COMPUTE_DERIV_COEF(SOL%TQUEUE(M),T(I),	&
						HTRIAL,DD0,DD1,DD2,DD3,DD9,DD10,DD11)
				DO J=1,NC
					CALL APLLY_DERIV_FORMULA(HTRIAL,SOL%QUEUE(YINT%COMPONENTS(J),:,M+1), &
						DD0,DD1,DD2,DD3,DD9,DD10,DD11,YINT%DT(I,J))
				END DO
			END IF
		ELSE
            PRINT *,' Error in DDVERK_VAL '
			PRINT *,' T must falls before the initial time and the end time.'
		END IF
		GUESS=M
	END DO
END FUNCTION DDVERK_VAL

SUBROUTINE DDVERK_USER(TVAL_P,YOFT,DYOFT)
	IMPLICIT NONE
    DOUBLE PRECISION ::TVAL_P,TVAL
    DOUBLE PRECISION, DIMENSION(G_NEQN) :: YOFT,DYOFT
	DOUBLE PRECISION :: T
    INTENT(IN) :: TVAL_P
    INTENT(OUT) :: YOFT,DYOFT
	!must be implemented
	INTEGER M
	INTEGER MP1
	DOUBLE PRECISION :: H,SIGN,TEMP2,ALAKI
	INTEGER :: J,L

    DOUBLE PRECISION AA(11,11),BB(11),CC(11)
	!C     .. Common blocks ..
    COMMON/METH/AA,BB,CC

	TVAL=TVAL_P

	IF (.NOT.(ALLOCATED(G_KOLD))) THEN
		PRINT *, ' DDVERK_USER may only be called by the user routines'
        PRINT *, ' DDES, BETA, HISTORY, EF, and CHNG. Use DDVERK_VAL'
        PRINT *, ' to perform the requested interpolation.'
        RETURN
    END IF

	T=G_T

	IF(TVAL>G_C(17)) THEN
		
		PRINT *, 'Warning in DDVERK_USER'
        PRINT *, ' TVAL is larger the current value of TRIGHT'
		PRINT *, ' TVAL= ',TVAL,' TRIGHT= ',G_C(17)

		!TVAL=G_C(17)
		
	END IF

	IF(TVAL>=G_T) THEN
		!IF(G_LEGAL_W) THEN
		!	CALL INTRP(G_NEQN,G_NW,G_T,G_Y,TVAL,YOFT,G_C(17)-G_T,G_W)
		!	CALL DERIV(G_NEQN,G_NW,G_T,G_Y,TVAL,DYOFT,G_C(17)-G_T,G_W)
		!ELSE
		!	IF(G_T>G_C(25)) THEN
		!		M=G_C(27)-1
		!		IF(M==0) M=G_LOLD
		!		IF(M+1 <= G_LOLD) THEN
		!			MP1=M+1
		!		ELSE
		!			MP1=1
		!		END IF
		!		CALL INTRPG(G_NEQN,G_NW,MP1,G_TOLD(M),TVAL,YOFT,G_TOLD(MP1)-G_TOLD(M))
		!		CALL DERIVG(G_NEQN,G_NW,MP1,G_TOLD(M),TVAL,DYOFT,G_TOLD(MP1)-G_TOLD(M))
		!	ELSE
		!		YOFT=G_Y
		!		DYOFT=0
		!	END IF
		!END IF
!		IF(.FALSE.) THEN 
		ALAKI=G_C(25)
		IF(G_T>G_C(25)) THEN
			M=G_C(27)-1
			IF(M==0) M=G_LOLD
			IF(M+1 <= G_LOLD) THEN
				MP1=M+1
			ELSE
				MP1=1
			END IF
			CALL INTRPG(G_NEQN,G_NW,MP1,G_TOLD(M),TVAL,YOFT,G_TOLD(MP1)-G_TOLD(M))
			CALL DERIVG(G_NEQN,G_NW,MP1,G_TOLD(M),TVAL,DYOFT,G_TOLD(MP1)-G_TOLD(M))
		ELSE
			!----------------------------
			SIGN=TVAL-G_T
			IF(G_K<=2) THEN
				DO J = 1, G_NEQN
					YOFT(J) = G_Y(J) +  SIGN*G_W(J,1)
					DYOFT(J) = G_W(J,1)
				END DO
				!-------------
			ELSE
				H=G_C(18)
				TEMP2=SIGN**2/(CC(G_K)*H)
				DO  L = 1, G_NEQN
					YOFT(L) = AA(G_K,1) * G_W(L,1)
					DYOFT(L) = AA(G_K,1) * G_W(L,1)
					DO J = 2, G_K-1
						YOFT(L) = YOFT(L) + AA(G_K,J) * G_W(L,J)
						DYOFT(L) =DYOFT(L) + AA(G_K,J) * G_W(L,J)
					END DO
					YOFT(L) = G_Y(L) + (SIGN - TEMP2) * G_W(L,1) &
						+ TEMP2 * YOFT(L) / CC(G_K)
					DYOFT(L) = G_W(L,1) - (G_W(L,1) - DYOFT(L) / CC(G_K)) *  &
						(2.D0 * SIGN / (H * CC(G_K)))
				END DO
			END IF
			!----------------------------
		END IF
	ELSEIF(TVAL>=G_C(25)) THEN
		CALL TFIND(TVAL,G_C(26),G_C(27),M)
		IF(M+1 <= G_LOLD) THEN
			MP1=M+1
		ELSE
			MP1=1
		END IF
		CALL INTRPG(G_NEQN,G_NW,MP1,G_TOLD(M),TVAL,YOFT,G_TOLD(MP1)-G_TOLD(M))
		CALL DERIVG(G_NEQN,G_NW,MP1,G_TOLD(M),TVAL,DYOFT,G_TOLD(MP1)-G_TOLD(M))
	ELSE
        PRINT *,' Error. T falls before the initial time.'
        PRINT *,' There is insufficient information in the solution'
        PRINT *,' queue to perform the interpolation requested in'
        PRINT *,' DDVERK_USER. Use your history function to obtain the'
        PRINT *,' interpolated solution requested in DDVERK_USER.'
	END IF

	RETURN
END SUBROUTINE DDVERK_USER

SUBROUTINE DDVERK_IQUEUE(IQ)
	IMPLICIT NONE
	INTEGER IQ
	INTENT(OUT) :: IQ
	IF (.NOT.(ALLOCATED(G_C))) THEN
		PRINT *, ' DDVERK_USER may only be called by the user routines'
        PRINT *, ' DDES, BETA, HISTORY, EF, and CHNG. Use DDVERK_VAL'
        PRINT *, ' to perform the requested interpolation.'
        RETURN
    END IF
	IQ=G_C(27)-1
	RETURN
END SUBROUTINE DDVERK_IQUEUE


FUNCTION FIND_LEFT_T_INDEX(QSIZE,TQUEUE,GUESS,TOUT) RESULT(M)
	IMPLICIT NONE
	INTEGER :: QSIZE,GUESS
	DOUBLE PRECISION :: TOUT
	DOUBLE PRECISION,DIMENSION(QSIZE)::TQUEUE
	INTEGER :: M
	INTEGER :: LEFT_INDEX,RIGHT_INDEX,MID_INDEX

	IF((GUESS<1).OR.(GUESS>QSIZE-1)) GUESS=1
	IF((TOUT >=TQUEUE(GUESS)).AND.(TOUT < TQUEUE(GUESS+1))) THEN
		M=GUESS
	ELSE
		IF(TOUT==TQUEUE(QSIZE)) THEN
			M=QSIZE-1
		ELSE
			IF((TOUT>=TQUEUE(1)).AND.(TOUT<TQUEUE(QSIZE))) THEN
				LEFT_INDEX=1
				RIGHT_INDEX=QSIZE
				DO WHILE(RIGHT_INDEX-LEFT_INDEX > 1)
					MID_INDEX=(LEFT_INDEX+RIGHT_INDEX)/2
					IF(TOUT < TQUEUE(MID_INDEX)) THEN
						RIGHT_INDEX=MID_INDEX
					ELSE
						LEFT_INDEX=MID_INDEX
					END IF
				END DO
				M=LEFT_INDEX
			ELSE
				M=0
			END IF
		END IF
	END IF
END FUNCTION FIND_LEFT_T_INDEX
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
!-----------------------------------------------
!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------
      SUBROUTINE APFORM(N,NW,MU,OMEGA,S,T,Y,H,TOL,W,C, &
          FCN,SIGMA,YINIT)

!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 G_KOLD              *
!C     MU     Number of solution delay terms                           *
!C     OMEGA  Number of derivative delay terms                         *
!C     G_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     G_TOLD   Queue of T values                                        *
!C     G_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 G_TOLD(c(27)-1). But c(27) may be one.                  *
!C       In this case, we have to use G_TOLD(G_LOLD - 1)                   *
!C**********************************************************************

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

!C     .. Array Arguments ..
      DOUBLE PRECISION  Y(N),C(G_CL)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

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

!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 MAX,DFLOAT,IDINT

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

!C     .. Data statements ..
      DATA MAXNIT/5/
!C     .. Set up the nubmer of delay terms
	  G_LEGAL_W=.FALSE.
	  G_K=0
	  !W(:,2:11)=0.D0
      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(TEMP,W(1,16),G_W15(1))

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

         IW = 20 + I

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

!C ... Cut the stepseize to half and integrate again
            IF (G_W15(I) - TEMP .GE. 1.D0) THEN
               PRINT *, 'ADVANCED TOO MUCH (1)!!!'
               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 = G_W15(I) - T
!C
!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                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 = G_LOLD
               CALL INTRPG(N,NW,TMP,G_TOLD(TMP-1),G_W15(I),W(1,IW), &
                   T-G_TOLD(TMP-1))

!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 (G_W15(I) .GT. TEMP) THEN
            PRINT *, '**** WARNING:APFORM *****'
            PRINT *, 'Advanced Delay at Stage  2'
            PRINT *, 'Ti =', TEMP, ' Delay Argument =', G_W15(I)
            PRINT *

!C ... Cut the stepseize to half and integrate again
            IF (G_W15(I) - TEMP .GE. 1.D0) THEN
               PRINT *, 'ADVANCED TOO MUCH (2)!!!'
               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 = G_W15(I) - T
!C
!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                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 = G_LOLD
               CALL DERIVG(N,NW,TMP,G_TOLD(TMP-1),G_W15(I),W(1,IW), &
                   T-G_TOLD(TMP-1))

!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(TEMP,W(1,16),W(1,21),W(1,2))

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

         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(TEMP,W(1,16),G_W15(1))

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

            IW = 20 + I

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

!C ... Cut the stepseize to half and integrate again
               IF (G_W15(I)-TEMP .GE. 1.D0) THEN
                  PRINT *, 'ADVANCED TOO MUCH (3)!!!'
                  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 = G_W15(I) - 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,C(25),C(26),C(27),G_W15(I), &
                   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 = G_LOLD
                  CALL INTRPG(N,NW,TMP,G_TOLD(TMP-1),G_W15(I),W(1,IW), &
                      T-G_TOLD(TMP-1))

!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 (G_W15(I) .GT. TEMP) THEN
               PRINT *, '**** WARNING:APFORM *****'
               PRINT *, 'AdvanceD Delay at Stage', K
               PRINT *, 'Ti =', TEMP, ' Delay Argument =', G_W15(I)
               PRINT *

!C ... Cut the stepseize to half and integrate again
               IF (G_W15(I)-TEMP .GE. 1.D0) THEN
                  PRINT *, 'ADVANCED TOO MUCH (4)!!!'
                  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 = G_W15(I) - T

!C ... In case there is no vanishing delay
            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                   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 = G_LOLD
                  CALL DERIVG(N,NW,TMP,G_TOLD(TMP-1),G_W15(I),W(1,IW), &
                      T-G_TOLD(TMP-1))

!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(TEMP,W(1,16),W(1,21),W(1,K))
         !CALL FCN(TEMP,W(:,16),W(1:G_NEQN,21:20+G_NLAGS),W(:,K))

 140  CONTINUE
	  G_K=S

!C***** End of the 1st phase *****
	
	!G_LEGAL_W=.TRUE.

!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
			G_K = K

            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(TEMP,W(1,16),G_W15(1))

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

               IW = 20 + I

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

!C ... Cut the stepseize to half and integrate again
                  IF (G_W15(I)-TEMP .GE. 1.D0) THEN
                     PRINT *, 'ADVANCED TOO MUCH (5)!!!'
                     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 = G_W15(I) - T

!C ... In case there is no vanisning delay
               IF (SIGN .LE. 0.D0) THEN
                  CALL EVALAGS(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                      W,YINIT)
               ELSE 
                  CALL INTRP(N,NW,T,Y,G_W15(I),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 (G_W15(I) .GT. TEMP) THEN
                  PRINT *, '**** WARNING:APFORM *****'
                  PRINT *, '    Advanced Delay at Stage ', K
                  PRINT *, 'Ti =', TEMP, ' Delay Argument =', G_W15(I)
                  PRINT *

!C ... Cut the stepseize to half and integrate again
                  IF (G_W15(I)-TEMP .GE. 1.D0) THEN
                     PRINT *, 'ADVANCED TOO MUCH (6)!!!'
                     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 = G_W15(I) - T

!C ... In case there is no vanisning delay
               IF (SIGN .LE. 0.D0) THEN
                  CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                      W,YINIT)
               ELSE 

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

!C ... Recompute K-value

            CALL FCN(TEMP,W(1,16),W(1,21),W(1,K))

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

 210     CONTINUE
		 G_K = S
         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(TEMP,W(1,12),G_W15(1))

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

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

 330     CONTINUE

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

            IW = 20 + I
            SIGN = G_W15(I) - T
            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                           W,YINIT)
            ELSE 

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

!C ... Recompute K-value
         CALL FCN(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 SUBROUTINE APFORM


!C *** End of APFORM ***






      SUBROUTINE APFORM_ASM(N,NW,MU,OMEGA,S,T,Y,H,TOL,W,C &
          ,FCN,SIGMA,YINIT)

!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 G_KOLD              *
!C     MU     Number of solution delay terms                           *
!C     OMEGA  Number of derivative delay terms                         *
!C     G_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     G_TOLD   Queue of T values                                        *
!C     G_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 G_TOLD(c(27)-1). But c(27) may be one.                  *
!C       In this case, we have to use G_TOLD(G_LOLD - 1)                   *
!C**********************************************************************

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

!C     .. Array Arguments ..
      DOUBLE PRECISION  Y(N),C(G_CL)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

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

!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 MAX,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
	  G_LEGAL_W=.FALSE.
	  !W(:,2:11)=0.D0
      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(TEMP,W(1,16),G_W15(1))

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

         IW = 20 + I

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

!C ... Cut the stepseize to half and integrate again
            IF (G_W15(I) - TEMP .GE. 1.D0) THEN
               PRINT *, 'ADVANCED TOO MUCH (7)!!!'
               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 = G_W15(I) - T
!C
!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                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 = G_LOLD
               CALL INTRPG(N,NW,TMP,G_TOLD(TMP-1),G_W15(I),W(1,IW), &
                   T-G_TOLD(TMP-1))

!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 (G_W15(I) .GT. TEMP) THEN
            PRINT *, '**** WARNING:APFORM *****'
            PRINT *, 'Advanced Delay at Stage  2'
            PRINT *, 'Ti =', TEMP, ' Delay Argument =', G_W15(I)
            PRINT *

!C ... Cut the stepseize to half and integrate again
            IF (G_W15(I) - TEMP .GE. 1.D0) THEN
               PRINT *, 'ADVANCED TOO MUCH (8)!!!'
               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 = G_W15(I) - T
!C
!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                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 = G_LOLD
               CALL DERIVG(N,NW,TMP,G_TOLD(TMP-1),G_W15(I),W(1,IW), &
                   T-G_TOLD(TMP-1))

!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(TEMP,W(1,16),W(1,21),W(1,2))

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

         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(TEMP,W(1,16),G_W15(1))

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

            IW = 20 + I

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

!C ... Cut the stepseize to half and integrate again
               IF (G_W15(I)-TEMP .GE. 1.D0) THEN
                  PRINT *, 'ADVANCED TOO MUCH (9)!!!'
                  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 = G_W15(I) - 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,C(25),C(26),C(27),G_W15(I), &
                   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 = G_LOLD
                  CALL INTRPG(N,NW,TMP,G_TOLD(TMP-1),G_W15(I),W(1,IW), &
                      T-G_TOLD(TMP-1))

!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 (G_W15(I) .GT. TEMP) THEN
               PRINT *, '**** WARNING:APFORM *****'
               PRINT *, 'AdvanceD Delay at Stage', K
               PRINT *, 'Ti =', TEMP, ' Delay Argument =', G_W15(I)
               PRINT *

!C ... Cut the stepseize to half and integrate again
               IF (G_W15(I)-TEMP .GE. 1.D0) THEN
                  PRINT *, 'ADVANCED TOO MUCH (10)!!!'
                  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 = G_W15(I) - T

!C ... In case there is no vanishing delay
            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                   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 = G_LOLD
                  CALL DERIVG(N,NW,TMP,G_TOLD(TMP-1),G_W15(I),W(1,IW), &
                      T-G_TOLD(TMP-1))

!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(TEMP,W(1,16),W(1,21),W(1,K))

 140  CONTINUE
	  G_K = S

!C***** End of the 1st phase *****
	  !G_LEGAL_W=.TRUE.
!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
			G_K = K

            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(TEMP,W(1,16),G_W15(1))

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

               IW = 20 + I

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

!C ... Cut the stepseize to half and integrate again
                  IF (G_W15(I)-TEMP .GE. 1.D0) THEN
                     PRINT *, 'ADVANCED TOO MUCH (11)!!!'
                     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 = G_W15(I) - T

!C ... In case there is no vanisning delay
               IF (SIGN .LE. 0.D0) THEN
                  CALL EVALAGS(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                      W,YINIT)
               ELSE 
                  CALL INTRP(N,NW,T,Y,G_W15(I),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 (G_W15(I) .GT. TEMP) THEN
                  PRINT *, '**** WARNING:APFORM *****'
                  PRINT *, '    Advanced Delay at Stage ', K
                  PRINT *, 'Ti =', TEMP, ' Delay Argument =', G_W15(I)
                  PRINT *

!C ... Cut the stepseize to half and integrate again
                  IF (G_W15(I)-TEMP .GE. 1.D0) THEN
                     PRINT *, 'ADVANCED TOO MUCH (12)!!!'
                     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 = G_W15(I) - T

!C ... In case there is no vanisning delay
               IF (SIGN .LE. 0.D0) THEN
                  CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                      W,YINIT)
               ELSE 

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

!C ... Recompute K-value

            CALL FCN(TEMP,W(1,16),W(1,21),W(1,K))

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

 210     CONTINUE
		 G_K = S
         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(TEMP,W(1,12),G_W15(1))

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

            IW = 20 + I
            SIGN = G_W15(I) - T

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

            END IF
 330     CONTINUE

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

            IW = 20 + I
            SIGN = G_W15(I) - T

            IF (SIGN .LE. 0.D0) THEN
               CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                           W,YINIT)
            ELSE 

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

!C ... Recompute K-value
         CALL FCN(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
		 G_K = K
         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(temp,w(1,16),G_W15(1))

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

            iw = 20 + i
            sign = G_W15(i) - t

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

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

            iw = 20 + i
            sign = G_W15(i) - t

            if (sign .le. 0.d0) then
               call evalagd(n,nw,iw,c(25),c(26),c(27),G_W15(i), &
                           w,yinit)
            else 
               call deriv(n,nw,t,y,G_W15(i),w(1,iw),h,w)
            end if
 520     continue

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

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

!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 SUBROUTINE APFORM_ASM


!C *** End of APFORM_ASM ***








      SUBROUTINE DDVERK(N,NW,MU,OMEGA,IND,T,TEND,Y,TOL, &
          W,C,FCN,SIGMA,YINIT)

!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,G_LOLD,T,TEND,Y,TOL,G_TOLD,G_KOLD,W,C,    *
!C    +            FCN,SIGMA,YINIT,DYINIT)                               *
!C                                                                       *
!C     INTEGER N,NW,IND,G_LOLD,MU,OMEGA                                    *
!C     DOUBLE PRECISION T,TEND,Y(N),TOL,G_TOLD(G_LOLD),G_KOLD(NW,7,G_LOLD)       *
!C     DOUBLE PRECISION W(NW,20+MU+OMEGA),C(G_CL)                          *
!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 G_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    G_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    G_TOLD   Queue of independent values T                               *
!C    G_KOLD   Queue of derivative and solution approximations             *
!C           The first dimension of G_KOLD must be NW                      *
!C           The second dimension of G_KOLD must be at least 7             *
!C           The third dimension of G_KOLD must be G_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 DSLO!C                                       *
!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
      DOUBLE PRECISION T,TEND,TOL

!C     .. Array Arguments ..
      DOUBLE PRECISION Y(N)
      DOUBLE PRECISION C(G_CL)
  	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W


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

!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(T,Y,.FALSE.)

!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 G_TOLD(1)
         C(25) = T
         G_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
		 G_NFAIL_STEPS=G_NFAIL_STEPS+C(23)
         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,T,TEND,Y,TOL,W,C, &
                 FCN,SIGMA,YINIT) 

!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) &
		  .AND. G_TRACK_DISCONTINUITIES) THEN

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

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

         !print *, 'ind2 =', ind2 !HZ

         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
			!G_C6_ONLY_NEXT=TS-C(17) !HZ
			!CALL STORE(N,NW,C(17),Y,C(26),C(27),W) !HZ
			!GO TO 500 ! HZ

!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. DABS(TEND)*C(10))) THEN
               CALL INTRP(N,NW,T,Y,TEND,W(1,13),C(18),W)
               CALL STORED(N,NW,T,TEND,Y,W(1,13),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 = MIN(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,C(17),Y,C(26),C(27),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         
			G_NFAIL_STEPS=G_NFAIL_STEPS+C(23)
            c(23) = 0.d0
   
            CALL STORED(N,NW,T,TS,Y,W(1,13),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,C(17),Y,C(26),C(27),W) 

         END IF

      ELSE IF (IND .EQ. 5) THEN

         CALL STORE(N,NW,C(17),Y,C(26),C(27),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
		 G_NFAIL_STEPS=G_NFAIL_STEPS+C(23)
         c(23) = 0.d0
         call intrp(n,nw,t,y,c(33),w(1,13),c(18),w)
         call stored(n,nw,t,c(33),y,w(1,13),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
		 G_NFAIL_STEPS=G_NFAIL_STEPS+C(23)
         c(23) = 0.d0
         ind = 5
         do 370 i = 1, n
            w(i,12) = G_KOLD(i,3,idint(c(27)))
            w(i,8) = G_KOLD(i,4,idint(c(27)))
 370     continue
         c(17) = G_TOLD(idint(c(27)))
         c(18) = c(17) - G_TOLD(idint(c(27)) - 1)
         c(14) = c(18)
         c(29) = 0.d0

      end if

 500  CONTINUE
      IF (GOBACK) RETURN

      GO TO 1000

10000 format(f20.12)

      END SUBROUTINE DDVERK


!C *** End of DDVERK ***





      SUBROUTINE DEFECT(N,NW,MU,OMEGA,T,TS,Y,DFCT,W,C,	&
                       FCN,SIGMA,YINIT)

!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 G_KOLD             *
!C    MU     Number of solution delay terms                          *
!C    OMEGA  Number of derivative delay terms                        *
!C    G_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    G_TOLD   Queue of independent values                             *
!C    G_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********************************************************************

	  IMPLICIT NONE
!C     .. Scalar Arguments ..
      INTEGER N,NW,MU,OMEGA
      DOUBLE PRECISION T,TS,DFCT

!C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),C(G_CL)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

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

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

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

!C     .. Intrinsic Functions ..
      INTRINSIC MAX,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)
	  !PRINT *,'DiffrenceY(1)=',Y(1)-W(1,2)
!C ... Compute argument of delay value
      CALL SIGMA(TS,W(1,2),G_W15(1))

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

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

         SIGN = G_W15(I) - T

!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                W,YINIT)

!C ... In case there is vanishing delay
         ELSE
            CALL INTRP(N,NW,T,Y,G_W15(I),W(1,IW),H,W)
         END IF
 110  CONTINUE
       
!C ...  Compute P'(Delay Argument)
      DO 115 I = MU + 1, ND

         IW = 20 + I

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

         SIGN = G_W15(I) - T

!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                W,YINIT)

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

!C ... Compute F(TS, P(TS), P(DEL ARG), P'(DELAY ARG))
      CALL FCN(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              1
      DO 200 K = 1, N
		IF(HZ_ISNAN(W(K,14)) .OR. (HZ_ISNAN(W(K,14)-W(K,14)))) THEN
			 TEMP = HUGE(1.D0)
		ELSE
		     TEMP = MAX(TEMP,DABS(W(K,14)))
		END IF
 200  CONTINUE
      GO TO 500
 210  IF (C(1) .NE. 2.D0) GO TO 240

!C ... Relative error control              2
      DO 230 K = 1, N
		IF(HZ_ISNAN(W(K,14)) .OR. (HZ_ISNAN(W(K,14)-W(K,14)))) THEN
			 TEMP = HUGE(1.D0)
		ELSE
         TEMP = MAX(TEMP, DABS(W(K,14)/W(K,18)))
		END IF
 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)))  3
      DO 250 K = 1, N
         TEMP = MAX(TEMP, DABS(W(K,14)) / MAX(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+33),ABS(W(K,18)))   4
      DO 270 K = 1, N
		IF(HZ_ISNAN(W(K,14)) .OR. (HZ_ISNAN(W(K,14)-W(K,14)))) THEN
			 TEMP = HUGE(1.D0)
		ELSE
		     TEMP = MAX(TEMP, DABS(W(K,14)) / MAX(C(K+33), &
			     DABS(W(K,18))))
		END IF
 270  CONTINUE
      GO TO 500
 280  IF (C(1) .NE. 5.D0) GO TO 300

!C ... Weights are 1/C(K+33)            5
      DO 290 K = 1, N
		IF(HZ_ISNAN(W(K,14)) .OR. (HZ_ISNAN(W(K,14)-W(K,14)))) THEN
			 TEMP = HUGE(1.D0)
		ELSE
         TEMP = MAX(TEMP, DABS(W(K,14)/C(K+33)))
		END IF
 290  CONTINUE
      GO TO 500

300  CONTINUE
!C ... Default case - Weights are 1/MAX(1,ABS(Y'(K)))    other > 6 or 0
      DO 310 K = 1, N
		IF(HZ_ISNAN(W(K,14)) .OR. (HZ_ISNAN(W(K,14)-W(K,14)))) THEN
			 TEMP = HUGE(1.D0)
		ELSE
         TEMP = MAX(TEMP, DABS(W(K,14)) / MAX(1.D0, &
             DABS(W(K,18))) )
		END IF
 310  CONTINUE
 500  CONTINUE

      DFCT = TEMP


      RETURN
      END SUBROUTINE DEFECT


!C *** End of DEFECT ***





      SUBROUTINE DEFECT2(N,NW,MU,OMEGA,T,TS,TS2,Y,DFCT,W, &
                       C,FCN,SIGMA,YINIT)

!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 G_KOLD             *
!C    MU     Number of solution delay terms                          *
!C    OMEGA  Number of derivative delay terms                        *
!C    G_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    G_TOLD   Queue of independent values                             *
!C    G_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
      DOUBLE PRECISION T,TS,TS2,DFCT

!C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),C(G_CL)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

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

!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 MAX,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(TS,W(1,2),G_W15(1))

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

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

         SIGN = G_W15(I) - T

!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                W,YINIT)

!C ... In case there is vanishing delay
         ELSE
            CALL INTRP(N,NW,T,Y,G_W15(I),W(1,IW),H,W)
         END IF
 110  CONTINUE
       
!C ...  Compute P'(Delay Argument)
      DO 115 I = MU + 1, ND

         IW = 20 + I

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

         SIGN = G_W15(I) - T

!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                W,YINIT)

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

!C ... Compute F(TS, P(TS), P(DEL ARG), P'(DELAY ARG))
      CALL FCN(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 = MAX(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 = MAX(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 = MAX(TEMP, DABS(W(K,14)) / MAX(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+33),ABS(W(K,18)))
      DO 270 K = 1, N
		IF(HZ_ISNAN(W(K,14)) .OR. (HZ_ISNAN(W(K,14)-W(K,14)))) THEN
			 TEMP = HUGE(1.D0)
		ELSE
		     TEMP = MAX(TEMP, DABS(W(K,14)) / MAX(C(K+33), &
			     DABS(W(K,18))))
		END IF
 270  CONTINUE
      GO TO 500
 280  IF (C(1) .NE. 5.D0) GO TO 300

!C ... Weights are 1/C(K+33)
      DO 290 K = 1, N
		IF(HZ_ISNAN(W(K,14)) .OR. (HZ_ISNAN(W(K,14)-W(K,14)))) THEN
			 TEMP = HUGE(1.D0)
		ELSE
         TEMP = MAX(TEMP, DABS(W(K,14)/C(K+33)))
		 END IF
 290  CONTINUE
      GO TO 500
 300  CONTINUE

!C ... Default case - Weights are 1/MAX(1,ABS(Y'(K)))
      DO 310 K = 1, N
		IF(HZ_ISNAN(W(K,14)) .OR. (HZ_ISNAN(W(K,14)-W(K,14)))) THEN
			 TEMP = HUGE(1.D0)
		ELSE
         TEMP = MAX(TEMP, DABS(W(K,14)) / MAX(1.D0, &
             DABS(W(K,18))) )
		END IF
 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(TS2,W(1,2),G_W15(1))

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

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

         SIGN = G_W15(I) - T

!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGS(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                W,YINIT)

!C ... In case there is vanishing delay
         ELSE
            CALL INTRP(N,NW,T,Y,G_W15(I),W(1,IW),H,W)
         END IF
 610  CONTINUE
       
!C ...  Compute P'(Delay Argument)
      DO 615 I = MU + 1, ND

         IW = 20 + I

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

         SIGN = G_W15(I) - T

!C ... In case there is no vanishing delay
         IF (SIGN .LE. 0.D0) THEN
            CALL EVALAGD(N,NW,IW,C(25),C(26),C(27),G_W15(I), &
                W,YINIT)

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

!C ... Compute F(TS2, P(TS2), P(DEL ARG), P'(DELAY ARG))
      CALL FCN(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 = MAX(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 = MAX(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 = MAX(TEMP, DABS(W(K,14)) / MAX(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+33),ABS(W(K,18)))
      DO 770 K = 1, N
         TEMP = MAX(TEMP, DABS(W(K,14)) / MAX(C(K+33), &
             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+33)
      DO 790 K = 1, N
         TEMP = MAX(TEMP, DABS(W(K,14)/C(K+33)))
 790  CONTINUE
      GO TO 900
 800  CONTINUE

!C ... Default case - Weights are 1/MAX(1,ABS(Y'(K))) >6 or 0
      DO 810 K = 1, N
         TEMP = MAX(TEMP, DABS(W(K,14)) / MAX(1.D0, &
             DABS(W(K,18))) )
 810  CONTINUE
 900  CONTINUE

!C      print *, 'ts =', ts, ' ts2 =', ts2
!C      print *, 'dfct1 =', dfct, ', defct2 =', temp

      dfct = MAX(dfct,temp)

      RETURN
      END SUBROUTINE DEFECT2


!C *** End of DEFECT2 ***




      SUBROUTINE DEFECTG(N,NW,MU,OMEGA,M,T1,TS,H,W,C, &
                        FCN,SIGMA,YINIT)


!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 G_KOLD               *
!C    MU     Number of solution delay terms                            *
!C    OMEGA  Number of derivative delay terms                          *
!C    G_LOLD   Length of queue                                           *
!C    M      Index for the queue, i.e, G_TOLD(M-1) <= TS <= G_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, G_TOLD(M) - G_TOLD(M-1)                        *
!C    W      Work space                                                *
!C    !C      Communication array                                       *
!C    G_TOLD   Queue of independent values                               *
!C    G_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,M,MU,OMEGA
      DOUBLE PRECISION T1,TS,H

!C     .. Array Arguments ..
      DOUBLE PRECISION C(G_CL)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

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

!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)

!C ... Compute argument of delay value
      CALL SIGMA(TS,W(1,2),G_W15(1))

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

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

         CALL EVALAGS(N,NW,20+I,C(25),C(26),C(27),G_W15(I), &
             W,YINIT)

 20   CONTINUE

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

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

         CALL EVALAGD(N,NW,20+I,C(25),C(26),C(27),G_W15(I), &
             W,YINIT)

 25   CONTINUE

!C ... Compute F(TS, P(TS), P(DEL ARG), P'(DELAY ARG))
      CALL FCN(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)

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

      RETURN
      END SUBROUTINE DEFECTG


!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 G_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)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

!C     .. Local Scalars ..
      INTEGER I

      double precision DD0,DD1,DD2,DD3,DD9,DD10,DD11


	  CALL COMPUTE_DERIV_COEF(T,TOUT,HTRIAL,DD0,DD1,DD2,DD3,DD9,DD10,DD11)

      DO I = 1, N
		  CALL APLLY_DERIV_FORMULA(HTRIAL,&
			(/y(i),w(i,1),w(i,12),w(i,8),w(i,9),w(i,10),w(i,11)/), &
			DD0,DD1,DD2,DD3,DD9,DD10,DD11 ,YPOUT(I))
	  END DO

      RETURN
      END SUBROUTINE DERIV


!C *** End of DERIV ***






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

!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 G_KOLD              *
!C     M      Index for the queue, i.e. G_TOLD(M-1) <= TOUT <= G_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     G_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)

!C     .. Local Scalars ..
      INTEGER I
      DOUBLE PRECISION DD0,DD1,DD2,DD3,DD9,DD10,DD11

	  CALL COMPUTE_DERIV_COEF(T,TOUT,HTRIAL,DD0,DD1,DD2,DD3,DD9,DD10,DD11)
      DO I = 1, N
		CALL APLLY_DERIV_FORMULA(HTRIAL,G_KOLD(I,:,M),   &
				DD0,DD1,DD2,DD3,DD9,DD10,DD11,YPOUT(I))
	  END DO

      RETURN
      END SUBROUTINE DERIVG


!C *** End of DERIVG ***

SUBROUTINE COMPUTE_DERIV_COEF(T,TOUT,HTRIAL,DD0,DD1,DD2,DD3,DD9,DD10,DD11)
      DOUBLE PRECISION ,INTENT(IN)::T,TOUT,HTRIAL
	  DOUBLE PRECISION ,INTENT(OUT) ::DD0,DD1,DD2,DD3,DD9,DD10,DD11
      DOUBLE PRECISION :: TAU,TAU2,TAU3,TAU4,TAU5
	  DOUBLE PRECISION :: T_RIGHT,T_LEFT
	  DOUBLE PRECISION  C9,C10,C11
      DATA C9/0.5D0/
      DATA C10/0.1428571428571428D0/
      DATA C11/0.2857142857142857D0/


      TAU = (TOUT - T)/HTRIAL
!IF(.FALSE.) THEN
!IF(HTRIAL > 1.D-7) THEN
!IF((.NOT. G_IS_NEUTRAL) .AND. (HTRIAL > 1.D-1)) THEN
!IF(.TRUE.) THEN
IF(.NOT. G_IS_NEUTRAL) THEN
      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
ELSE
!	PRINT *,'HTRIAL=',HTRIAL
		DD0=TAU
		DD1=T
END IF
		 	  
	  RETURN
END SUBROUTINE COMPUTE_DERIV_COEF

SUBROUTINE APLLY_DERIV_FORMULA(HTRIAL,K,DD0,DD1,DD2,DD3,DD9,DD10,DD11,YPOUT_I)
	  IMPLICIT NONE
      DOUBLE PRECISION ,INTENT(IN)::HTRIAL
	  DOUBLE PRECISION ,DIMENSION(7),INTENT(IN):: K
	  DOUBLE PRECISION ,INTENT(IN)::DD0,DD1,DD2,DD3,DD9,DD10,DD11
      DOUBLE PRECISION ,INTENT(OUT)::YPOUT_I
	  
	  DOUBLE PRECISION ,DIMENSION(5)::MY_K
	  INTEGER :: MY_DD0

	  DOUBLE PRECISION ,DIMENSION(5)::Z
	  DOUBLE PRECISION ,DIMENSION(5)::x
	  DOUBLE PRECISION ::x0,T
	  DOUBLE PRECISION  C9,C10,C11
	  INTEGER J,L

      DATA C9/0.5D0/
      DATA C10/0.1428571428571428D0/
      DATA C11/0.2857142857142857D0/

    	T=DD1

!IF(.FALSE.) THEN
!IF(HTRIAL > 1.D-7) THEN
!IF((.NOT. G_IS_NEUTRAL) .AND. (HTRIAL > 1.D-1)) THEN
!IF(.TRUE.) THEN
IF(.NOT. G_IS_NEUTRAL) THEN
      YPOUT_I = (DD0 * K(1) + DD2 * K(3)) / HTRIAL +	&
             DD1 * K(2) + DD3 * K(4) + DD9 *				&
             K(5) + DD10 * K(6) + DD11 * K(7)
ELSE
	  x=(/0.D0,C10,C11,C9,1.D0/)*HTRIAL+T
	  Z=(/K(2),K(6),K(7),K(5),K(4)/)
	  DO J=1,4
		DO L=J+1,5
			Z(L)=(Z(L)-Z(J))/(x(L)-x(J))
		END DO
	  END DO
	  
	  x0=DD0*HTRIAL+T
	  YPOUT_I=Z(5)
	  DO J=1,4
		YPOUT_I=YPOUT_I*(x0-x(5-J))+Z(5-J)
	  END DO
END IF

	  
END SUBROUTINE APLLY_DERIV_FORMULA


      SUBROUTINE DSLOC(N,NW,MU,OMEGA,IND2,T,TS,Y,YS,	&
                      TOL,W,C,FCN,SIGMA,YINIT)

!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 G_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    G_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    G_TOLD   Queue of independent values                             *
!C    G_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,MU,OMEGA
      DOUBLE PRECISION T,TS,TOL

!C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),YS(N)
      DOUBLE PRECISION C(G_CL)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

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

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

!C     .. External Subroutines ..
!      EXTERNAL DEFECT,INTRP

!C     .. Intrisic Functions ..
      INTRINSIC MAX,MIN

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

	  !PRINT *, "DSLOC CALLED"

!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 = MAX(10.D0*TOL, K2*DFCT)
      DELTA = MIN(1.D0, DELTA)

      !TS = TL !HZ

	  !CALL DEFECT(N,NW,MU,OMEGA,T,TS,Y,DFCT_L,W,C,	& !HZ
       !          FCN,SIGMA,YINIT) !ZH
	  
	  TS = TH
      CALL DEFECT(N,NW,MU,OMEGA,T,TS,Y,DFCT,W,C,	&
                 FCN,SIGMA,YINIT)

!C ... Set EPS to be TOL / DEFECT at Tn-1 + 2H for stopping criterion
!C      EPS = MIN(C(18), TOL / DFCT)
      EPS = MIN(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,T,TS,Y,DFCT,W,C,	&
                 FCN,SIGMA,YINIT)

!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 = MIN(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 = MAX(10.D0*TOL, K2*DFCT)
         DELTA = MIN(1.D0, DELTA)
      END IF

      IND2 = 2
		
!C ... Loop --- Bisection search

 1000 CONTINUE

         IF ((TH - TL) .LT. EPS) THEN
         !IF ((TH - TL) .LT. G_UROUND*1.D2*MAX(ABS(TH),ABS(TL))) THEN !HZ
            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,T,TS,Y,DFCT,W,C,	&
                    FCN,SIGMA,YINIT)

         IND2 = IND2 + 1

         IF (DFCT .LT. DELTA) THEN
         !IF (DFCT .LT. 1D2*(MAX(DFCT_L,TINY(1.D0)))) THEN !HZ
            TL = TS
         ELSE
            TH = TS
            EPS = MIN(C(18)/20.d0, TOL / DFCT)
         END IF

         GOTO 1000

!C     End loop --- Bisection search

      END SUBROUTINE DSLOC


!C *** End of DSLOC ***





      SUBROUTINE EVALAGS(N,NW,IW,X0,IFIRST,ILAST,DELARG,	&
                       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 G_KOLD           *
!C     IW      Index for stroge space                                *
!C     G_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     G_TOLD    Queue of the past independent variables               *
!C     G_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
      DOUBLE PRECISION X0,IFIRST,ILAST,DELARG

!C     .. Array Arguments ..
      DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

!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(DELARG,W(1,IW), .FALSE.)
!		PRINT *, "SEE CALLED YINIT"		
!C ... Otherwise find the index in the queue for the delay function
      ELSE 

         CALL TFIND(DELARG,IFIRST,ILAST,M)

!C ... In case the data is on the mesh
         IF (G_TOLD(M) .EQ. DELARG) THEN
            DO 120 J = 1, N
               W(J,IW) = G_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. G_LOLD) THEN
!C         IF (M .EQ. G_LOLD) THEN
            H = G_TOLD(1) - G_TOLD(M)
            CALL INTRPG(N,NW,1,G_TOLD(M),DELARG,W(1,IW),H)

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

      RETURN
      END SUBROUTINE EVALAGS


!C *** End of EVALAGS ***


      SUBROUTINE EVALAGD(N,NW,IW,X0,IFIRST,ILAST,DELARG, &
                       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 G_KOLD           *
!C     IW      Index of storage space                                *
!C     G_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     G_TOLD    Queue of the past independent variables               *
!C     G_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
      DOUBLE PRECISION X0,IFIRST,ILAST,DELARG

!C     .. Array Arguments ..
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W


!C     .. Subroutine Arguments ..
      EXTERNAL YINIT

!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 YINIT(N,DELARG,W(1,IW),.TRUE.)
!		PRINT *, "SEE CALLED DYINIT"		

!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(DELARG,IFIRST,ILAST,M)


!C ... In case the data is on the mesh
         IF (G_TOLD(M) .EQ. DELARG) THEN
            DO 120 J = 1, N
               W(J,IW) = G_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. G_LOLD) THEN
!C         IF (M .EQ. G_LOLD) THEN
            H = G_TOLD(1) - G_TOLD(M)
            CALL DERIVG(N,NW,1,G_TOLD(M),DELARG,W(1,IW),H)

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

      RETURN
      END SUBROUTINE EVALAGD


!C *** End of EVALAGD ***




      SUBROUTINE INTRP(N,NW,T,Y,TOUT,YOUT,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 G_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     YOUT  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),YOUT(N)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W
	  DOUBLE PRECISION D0,D1,D2,D3,D9,D10,D11

!C     .. Local Scalars ..
      INTEGER I
	  CALL COMPUTE_INTRP_COEF(T,TOUT,HTRIAL,D0,D1,D2,D3,D9,D10,D11)
      DO I = 1, N
		  CALL APLLY_INTRP_FORMULA(HTRIAL,&
			(/y(i),w(i,1),w(i,12),w(i,8),w(i,9),w(i,10),w(i,11)/), &
			D0,D1,D2,D3,D9,D10,D11 ,YOUT(I))
	  END DO

      RETURN
      END SUBROUTINE INTRP


!C *** End of INTRP ***






      SUBROUTINE INTRPG(N,NW,M,T,TOUT,YOUT,HTRIAL)
		IMPLICIT NONE

!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 G_KOLD              *
!C     M      Index for the queue, i.e. G_TOLD(M-1) <= TOUT <= G_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     G_KOLD   Queue of K values and solution values                    *
!C                                                                     *
!C**********************************************************************

!C     .. Scalar Arguments ..
      INTEGER N,NW,M
      DOUBLE PRECISION T,TOUT,HTRIAL
	  DOUBLE PRECISION D0,D1,D2,D3,D9,D10,D11

!C     .. Array Arguments ..
      DOUBLE PRECISION YOUT(N)

	
!C     .. Local Scalars ..
      INTEGER I
	  CALL COMPUTE_INTRP_COEF(T,TOUT,HTRIAL,D0,D1,D2,D3,D9,D10,D11)

      DO I = 1, N
		  CALL APLLY_INTRP_FORMULA(HTRIAL,G_KOLD(I,:,M), &
			D0,D1,D2,D3,D9,D10,D11,YOUT(I))
	  END DO
	
      RETURN
      END  SUBROUTINE INTRPG

SUBROUTINE COMPUTE_INTRP_COEF(T,TOUT,HTRIAL,D0,D1,D2,D3,D9,D10,D11)
      DOUBLE PRECISION ,INTENT(IN)::T,TOUT,HTRIAL
	  DOUBLE PRECISION ,INTENT(OUT) ::D0,D1,D2,D3,D9,D10,D11
      DOUBLE PRECISION TAU,TAU2

      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

END SUBROUTINE COMPUTE_INTRP_COEF

SUBROUTINE APLLY_INTRP_FORMULA(HTRIAL,K,D0,D1,D2,D3,D9,D10,D11,YOUT_I)
	  IMPLICIT NONE
      DOUBLE PRECISION ,INTENT(IN)::HTRIAL
	  DOUBLE PRECISION ,DIMENSION(7),INTENT(IN):: K
	  DOUBLE PRECISION ,INTENT(IN)::D0,D1,D2,D3,D9,D10,D11
      DOUBLE PRECISION ,INTENT(OUT)::YOUT_I

      YOUT_I = D0 * K(1)+ D2 *K(3) + HTRIAL * (D1	&
             * K(2) + D3 * K(4) + D9 * K(5) + D10       &
             * K(6) + D11 * K(7))
END SUBROUTINE APLLY_INTRP_FORMULA

      SUBROUTINE RDMETH(N,NW,MU,OMEGA,IND,T,TEND,Y,TOL,	&
                       W,C,FCN,SIGMA,YINIT)
		IMPLICIT NONE
		INTEGER :: BAD_I

!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 G_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     G_LOLD   The length of queue                                      *
!C     T      Current independent value                                *
!C     Y      Current solution, i.e., at T                             *
!C     G_TOLD   Queue of T values                                        *
!C     G_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,MU,OMEGA
      DOUBLE PRECISION T,TEND,TOL

!C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),C(G_CL)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

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

!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
	  LOGICAL :: ACCEPTED_DEFECT

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

!C     .. Intrisic Functions ..
      INTRINSIC DABS,DFLOAT,MIN,MAX,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+33) = DABS(C(K+33))
 170  CONTINUE
 180  CONTINUE
 150  CONTINUE

!C ... Initialize rreb, dwarf, prev tend, flag, counts
      !C(10) = 16.D0**(-13)
	  C(10) = G_RREB
      !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(T,Y,G_W15(1))

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

            ELSE 
               CALL EVALAGS(N,NW,20+I,C(25),C(26),C(27),G_W15(I),	&
                          W,YINIT)
            END IF

 430     CONTINUE

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

            ELSE 
               
               CALL EVALAGD(N,NW,20+I,C(25),C(26),C(27),G_W15(I),	&
                          W,YINIT)
               
            END IF

 435     CONTINUE

         CALL FCN(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) = MIN(1.D-8,TOL/(C(22)+100))
!      C(13) = MAX(1.D-10,1.D4*T*G_UROUND)
      !C(13) = 1.D-24
      !C(13) = 1.D-8
	  C(13)=MAX(T,1.D1*TINY(1.D0))*G_UROUND*1.D2
      !C(13) = MIN(1.D-8,1.D4*T*G_UROUND)
!C      C(13) = 1.D-6
 440  CONTINUE
		!PRINT *,'ALAKI =',T
!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) = MIN(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
	  !IF (C(6).EQ.0.D0 .AND. C(5).EQ.0.D0) C(16) = TEND-T
	  !IF(G_C6_ONLY_NEXT /=0.D0) C(16)=MIN(C(16),G_C6_ONLY_NEXT) !HZ
	  !G_C6_ONLY_NEXT=0.D0 !HZ

!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(T,Y,G_W15(1))
         DO 470 I = 1, ND
            IF (((T-G_W15(I) ) .LT. DABS(T)*C(10))) THEN
				PRINT *, "T-G_W15(I) = ",T-G_W15(I)
               IF ( C(4) .EQ.  0.D0)  C(14) = 1.D2*C(13) 
               !IF ( C(4) .EQ.  0.D0)  C(14) = 1.D8*C(13) 
            END IF
 470     CONTINUE
      end if

      !IF (C(14) .EQ. 0.D0) C(14) = MIN(1.D0,C(16))*TOL**(1.D0/DFLOAT(P))
      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.
	  !IF(IND == 5) THEN
	!	C(14) = 1.5*
	!	GO TO 480
	 ! END IF
      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) = MAX(TEMP, 0.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
	  !IF(C(16) > 0.D0 ) THEN 
	  C(14) = MIN(C(14), C(16))
	  !END IF
!C ... Check against HMIN
      C(14) = MAX(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) = MIN(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,SSTAR,T,Y,H,TOL,W,C, &
                 FCN,SIGMA,YINIT) 

      else if (C(32) .eq. 3) then
         
         CALL APFORM_ASM(N,NW,MU,OMEGA,SSTAR,T,Y,H,TOL,W	&
             ,C,FCN,SIGMA,YINIT) 

      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,T,T+TAUS*H,Y,C(19),&
                 W,C,FCN,SIGMA,YINIT)
			 !TEMPTEMP=C(19) !HZ
             !CALL DEFECT(N,NW,MU,OMEGA,T,T+H,Y,C(19),& !HZ
                 !W,C,FCN,SIGMA,YINIT) !HZ
				 !C(19)=MAX(C(19),TEMPTEMP) !HZ

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

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

          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
    	ACCEPTED_DEFECT=.TRUE.
	    DO I=1,N
			!if((ABS(W(I,14)) > ABS(W(I,2))*tol) .OR. &
			if((ABS(W(I,14)) > ABS(W(I,2))*tol+C(I+33)) .OR. &
			!if((H*ABS(W(I,14)) > ABS(W(I,2))*tol+C(I+33)*H) .OR. &
			!if((H*ABS(W(I,14)) > ABS(W(I,2))*tol+C(I+33)) .OR. &
				HZ_ISNAN(W(I,14)) .OR. HZ_ISNAN(W(I,18))) then
				!.FALSE. ) then
				ACCEPTED_DEFECT=.FALSE.
				BAD_I=I
				EXIT
			end if
		END DO

         !if (c(19) .le. tol/H) then 
         !if (c(19) .le. tol/1.D0) then 
         if (ACCEPTED_DEFECT) then 
		 !.OR. (DABS(C(14)-C(13)) .LE. DABS(C(14))*1D-8)) then
!C         if ((dabs(w(1,14)) .le. tol) .and. 
!C     +        (MAX(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
				!w(:,2:11)=0.D0
            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
	  G_NFAIL_STEPS=G_NFAIL_STEPS+C(23)
      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
	  IF(IND==-3) THEN
		PRINT *,'BAD_I=', BAD_I
		PRINT *,'C(19)=',C(19),'H= ',C(18)
		PRINT *, 'dY=',W(BAD_I,2)
		PRINT *, 'Defect=',W(BAD_I,14)
	      IND = -3
		  !IND = 5
	  END IF

      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 FOLLOWING	&
		& 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 SUBROUTINE RDMETH


!C *** End RDMETH ***





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

!C*********************************************************************
!C                                                                    *
!C     Purpose - This routine is for storing T value, K-values, and   *
!C               solutions to G_TOLD and G_KOLD queue.                    *
!C               When K-values for Tn-1 to Tn and Yn have             *
!C               computed, store                                      *
!C               G_TOLD(ILAST+1) <--- Tn                                *
!C               G_KOLD(*,1,ILAST+1) <--- Yn-1                          *
!C               G_KOLD(*,2,ILAST+1) <--- K1                            *
!C               G_KOLD(*,3,ILAST+1) <--- Yn                            *
!C               G_KOLD(*,4,ILAST+1) <--- Ks+1(K8)                      *
!C               G_KOLD(*,5,ILAST+1) <--- Ks+2(K9)                      *
!C               G_KOLD(*,6,ILAST+1) <--- Ks+3(K10)                     *
!C               G_KOLD(*,7,ILAST+1) <--- Ks+4(K11)                     *
!C                                                                    *
!C*********************************************************************
!C                                                                    *
!C    N      Number of equation                                       *
!C    NW     First dimension of workspace W & queue G_KOLD              *
!C    G_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    G_TOLD   Queue of independent values                              *
!C    G_KOLD   Queue of k-values and solutions                          *
!C    W      Workspace                                                *
!C                                                                    *
!C*********************************************************************


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

!C     .. Array Arguments ..
      DOUBLE PRECISION Y(N)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

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

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


      IFST = IDINT(IFIRST)
      ILST = IDINT(ILAST)
	  !WRITE(*,*) 'STORE = ',T

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

!C ... CASE 2: IFST < ILST = G_LOLD
      ELSE IF (ILST .EQ. G_LOLD) THEN
		 IF(G_LOLD<P_MAX_G_LOLD_LENGTH) THEN
			CALL EXTEND_QUEUE
	         M = ILST + 1
		     ILAST = ILAST + 1.D0
		 ELSE
	         print *, '***Warning: STORE***' 
		     print *, '    the length of the queue is short!'
			 print *, ' T = ',T
	         M = 1
		     IFIRST = 2.D0
			 ILAST = 1.D0
		 END IF
!C ... CASE 3: ILST < IFST < G_LOLD
      ELSE IF (IFST .LT. G_LOLD) THEN
         M = ILST + 1
         IFIRST = IFIRST + 1.D0
         ILAST = ILAST + 1.D0

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

!C ... Store T value
      G_TOLD(M) = T

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

!C      print *, 'm =', m
!C      print *, 'G_TOLD =', G_TOLD(m)
!C      print *, 'G_KOLD1 =', (G_KOLD(1,i,m), i= 1, 7)
!C      print *, 'G_KOLD2 =', (G_KOLD(2,i,m), i= 1, 7)
!C      print *, 'G_KOLD3 =', (G_KOLD(3,i,m), i= 1, 7) 

      RETURN
      END SUBROUTINE STORE



!C *** End of STORE ***





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

!C*********************************************************************
!C                                                                    *
!C     Purpose - This routine is for storing T value, derivatives,    *
!C               and solutions to G_TOLD and G_KOLD queue just after      *
!C               finding discontinuity.                               *
!C                                                                    *
!C               G_TOLD(ILAST+1) <--- Ts                                *
!C               G_KOLD(*,1,ILAST+1) <--- Yn-1                          *
!C               G_KOLD(*,2,ILAST+1) <--- K1                            *
!C               G_KOLD(*,3,ILAST+1) <--- Ys                            *
!C               G_KOLD(*,4,ILAST+1) <--- Z'(Ts)                        *
!C               G_KOLD(*,5,ILAST+1) <--- Z'(T+C9*Hs)                   *
!C               G_KOLD(*,6,ILAST+1) <--- Z'(T+C10*Hs)                  *
!C               G_KOLD(*,7,ILAST+1) <--- Z'(T+C11*Hs)                  *
!C                                                                    *
!C*********************************************************************
!C                                                                    *
!C    N      Number of equation                                       *
!C    NW     First dimension of workspace W & queue G_KOLD              *
!C    G_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    G_TOLD   QUEUE OF INDEPENDENT VALUES                              *
!C    G_KOLD   QUEUE OF K-VALUES AND SOLUTIONS                          *
!C    W      WORKSPACE                                                *
!C    !C      COMMUNICATION ARRAY                                      *
!C                                                                    *
!C*********************************************************************

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

!C     .. Array Arguments ..
      DOUBLE PRECISION Y(N),YS(N)
      DOUBLE PRECISION C(G_CL)
	  DOUBLE PRECISION,DIMENSION(NW,G_NW2) :: W

!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 < G_LOLD
      IF ((ILST .GE. IFST) .AND. (ILST .LT. G_LOLD)) THEN
         M = ILST + 1
         C(27) = C(27) + 1.D0

!C ... CASE 2: IFST < ILST = G_LOLD
      ELSE IF (ILST .EQ. G_LOLD) THEN
		 IF(G_LOLD<P_MAX_G_LOLD_LENGTH) THEN
			CALL EXTEND_QUEUE
	         M = ILST + 1
		     ILAST = ILAST + 1.D0
		 ELSE
	         print *, '***Warning: STORED***' 
		     print *, '    the length of the queue is short!' 
			 M = 1
			C(26) = 2.D0
			C(27) = 1.D0
		END IF

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

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

!C ... Store T value
      G_TOLD(M) = TS

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

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

      RETURN
      END SUBROUTINE STORED


!C *** End of STORED ***

SUBROUTINE EXTEND_QUEUE
	IMPLICIT NONE
	DOUBLE PRECISION , DIMENSION(:),POINTER :: G_TOLD_TEMP
	DOUBLE PRECISION , DIMENSION(:,:,:),POINTER ::G_KOLD_TEMP
	ALLOCATE(G_KOLD_TEMP(G_NW,7,G_LOLD),G_TOLD_TEMP(G_LOLD))

	G_TOLD_TEMP=G_TOLD(1:G_LOLD)
	G_KOLD_TEMP=G_KOLD(1:G_NW,1:7,1:G_LOLD)
	DEALLOCATE(G_TOLD)
	DEALLOCATE(G_KOLD)
	ALLOCATE(G_TOLD(2*G_LOLD))
	ALLOCATE(G_KOLD(G_NW,7,2*G_LOLD))
	G_TOLD(1:G_LOLD)=G_TOLD_TEMP
	G_KOLD(1:G_NW,1:7,1:G_LOLD)=G_KOLD_TEMP
	!
	G_LOLD=2*G_LOLD

	DEALLOCATE(G_KOLD_TEMP,G_TOLD_TEMP)
END SUBROUTINE EXTEND_QUEUE



      SUBROUTINE TFIND(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 G_TOLD(M) <= T < G_TOLD(M+1), M is returned.               *
!C                                                               *
!C****************************************************************
!C                                                               *
!C     N       The length of queue                               *
!C     G_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     .. Local Scalars ..
      INTEGER I,J,K,IFST,ILST,RK,RKP,RKM

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


!C ... Initialize variables
	  N=G_LOLD
      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:  G_TOLD(RK) = T
         IF (G_TOLD(RK) .EQ. T) THEN
            M = RK
            RETURN

!C ... CASE 2:  G_TOLD(RK) < T
         ELSE IF (G_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: G_TOLD(RK) < T < G_TOLD(RK+1) THEN M = RK
            ELSE IF (T .LT. G_TOLD(RKP)) THEN
               M = RK
               RETURN

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

!C ... CASE 3: NOW T < G_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, ', G_TOLD(k) =', G_TOLD(k)
            PRINT *
            STOP

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

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

      END SUBROUTINE TFIND



!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 FUNCTION REALK
!C *** End of REALK
!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------
END MODULE DDVERK90_M