      IMPLICIT REAL*8(A-H,O,P,R-Z)
	INTEGER COUNT,CONUM,QUOTIENT,ODD
      COMMON /OPT/NT,NOT /TA/T0,T1,TI /OU2/NOUT
	COMMON /CON/CONSTVALUE,QUOTIENT,ODD
	CHARACTER CO*1
C
      CONSTVALUE = 1000000000.0
	CALL OFILE
C
      CALL PARM
C
      CALL INIT

	open(55, file='toleranceFile.dat',form='formatted')
	read(55,*) Q
	close(55)

	ABSERR = Q
      CALL OUT
      CALL PLO(0)
      NUM=1
C
      CALL RG2(T1,RELERR,ABSERR,IER,NUM)
      NOT=NUM
C
  180 CALL CFILE
C
      STOP
      END

C********************************************************
      SUBROUTINE SOLV(T,X)
C     SOLVING CONSTRAINT EQUATIONS
C     IF DEVIATION RATIO IS LESS THAN RELATIVE ERROR(D0),
C     THEN THE PROGRAM OBTAINS SOLUTION.
      DIMENSION DE(40),DF(40),X(130),DE1(40),DF1(40)
      REAL*8 DE,DF,DE1,DF1,D0,D1,D2,FU,T,X
      COMMON /ES/DE,DF,ND
      IF(ND.EQ.0) THEN
        RETURN
      ENDIF
      I0=ND/2
      D0=1.0D-08
C     D0 : RELATIVE ERROR
    1 ICHK=0
C     ICHK : INDICATOR FOR TERMINATION
C          =1  ONE MORE REITERATION
C          =0  TEMINATION
      DO 10 I=1,I0
        D1=DABS(DE(I))*D0
        I1=2*I-1
        DE1(I)=FU(I1,T,X)
        D2=DABS(DE1(I)-DE(I))
        IF(D2.GT.D1) THEN
          ICHK=1
        ENDIF
        DE(I)=DE1(I)
C
        D1=DABS(DF(I))*D0
        I2=2*I
        DF1(I)=FU(I2,T,X)
        D2=DABS(DF1(I)-DF(I))
        IF(D2.GT.D1) THEN
          ICHK=1
        ENDIF
        DF(I)=DF1(I)
C
   10 CONTINUE
C
      IF(ICHK.EQ.0) THEN
        RETURN
      ELSE
        GO TO 1
      ENDIF
C
      END
C
C*********************************
      SUBROUTINE PLO(I)
      CHARACTER LABEL(100)*5, CDIV*1
      DIMENSION OP(100),X(130)
      REAL*8 OP,H,X,T
      COMMON /OU1/OP /OU2/NOUT /OU3/LABEL /RP/H,X,T
      DATA CDIV /','/
C
      IF(I.EQ.0) THEN
        WRITE(10,100) 'TIME',(CDIV, LABEL(J),J=1,NOUT)
        WRITE(*,100) 'TIME',(CDIV, LABEL(J),J=1,NOUT)
  100 FORMAT(A,99(A1,A))
        WRITE(10,200) T,(CDIV, OP(J),J=1,NOUT)
	  WRITE(*,200) T,(CDIV, OP(J),J=1,NOUT)
C  200 FORMAT(D15.8,99(A1,D15.8))
  200 FORMAT(E15.8,99(A1,E15.8))
      ELSE
        WRITE(10,200) T,(CDIV, OP(J),J=1,NOUT)
	  WRITE(*,200) T,(CDIV, OP(J),J=1,NOUT)
      ENDIF

      RETURN
      END


C
      SUBROUTINE RG2(XE,RELERR,ABSERR,IER,NUM)
      IMPLICIT REAL*8(A-H,O-Z)
    
      DIMENSION YH0(130),YHN(130),YL0(130),YLN(130),DX(130)
     & ,ERR(130),ET(130),Y(130)
      COMMON /RP/H,Y,T /BO/NS,ING /TA/T0,T1,TI /OPT/NT,NOT
C
      COMMON /FO/DX

      DATA ITEMAX /1000000000/
      DATA EPSMIN / 1.0D-15 /

      TW=(T1-T0)/NOT
      TOUT=T+TW
      N=NS+ING

       DO 10 I=1,N
        YL0(I)=Y(I)
   10  CONTINUE

       CALL SOLV(T,YL0)
       CALL FUNC(T,YL0,N)
       TOLDY=10.0D+74
       DO 20 I=1,N
        IF(DX(I) .NE. 0.0D+0) THEN
         TOL=RELERR*DABS(YL0(I))+ABSERR
         TOLDY=MIN(TOLDY,TOL/DABS(DX(I)))
        END IF
   20  CONTINUE      
       HMIN=EPSMIN*(XE-T)
       IF (TOLDY .EQ. 10.0D+74) THEN
        H=HMIN
       ELSE
        H = MIN((XE - T) / 100.0D+0, TOLDY ** 0.2D+0)
       END IF

       ISTIFF=0
       ICOUNT=0
       IER=0

      DO 30 ITER=1,ITEMAX
        DO 35 I=1,N
          Y(I)=YL0(I)
   35   CONTINUE
        TT=T
        LINDEX = 1
        CALL FEHL(LINDEX,N,YL0,YLN,ERR)

        ERRET=0.0D+0
        DO 40 I=1,N
          ET(I)=RELERR*0.5D+0*(DABS(YL0(I))+DABS(YLN(I)))+ABSERR
          IF (ET(I) .EQ. 0.0D+0) THEN
            GO TO 20000
          ELSE
            ERRET=MAX(ERRET,DABS(ERR(I))/ET(I))
          END IF
   40   CONTINUE
        IF (ERRET .GE. 1.0D+0) THEN
          T=TT
          IF(ERRET .GE. 59049.0D+0) THEN
            H=0.1D+0*H
          ELSE
            H=0.9D+0*H/(ERRET ** 0.2D+0)
          END IF
          IF (H .LE. EPSMIN) GO TO 30000
        ELSEIF ((TT+H) .LT. XE) THEN
          DO 50 I=1,N
            YH0(I)=YL0(I)
            YL0(I)=YLN(I)
   50     CONTINUE
          H1=H
  100     IF(TOUT.LE.(TT+H1)) THEN 
            T=TT
            H=TOUT-T
            LINDEX = 0
            CALL FEHL(LINDEX,N,YH0,YHN,ERR)
            CALL OUT
            CALL PLO(NUM)
            NUM=NUM+1
            TOUT=TOUT+TW
            GO TO 100
          END IF
          H=H1
          T=TT+H

          IF (ERRET.LE.1.889568D-4) THEN
            H = 5.0D+0 * H
          ELSE
            H = 0.9D+0 * H / (ERRET ** 0.2D+0)
          END IF

        ELSE

  150     DO 200 I=1,N
            YH0(I)=YL0(I)
  200     CONTINUE
  300     IF(TOUT.LE.XE) THEN 
            T=TT
            H=TOUT-T
            LINDEX = 0
            CALL FEHL(LINDEX,N,YH0,YHN,ERR)
            CALL OUT
            CALL PLO(NUM)
            NUM=NUM+1
            TOUT=TOUT+TW
            GO TO 300
          ELSE
            T=TT
            H=XE-T
            LINDEX = 0
            CALL FEHL(LINDEX,N,YH0,YHN,ERR)
            CALL OUT
            CALL PLO(NUM)
            RETURN
          END IF
        END IF
   30 CONTINUE

      IER = 10000
      WRITE( * ,10001) X
10001 FORMAT(' ','(SUBR.-RKFDS) TROUBLE(TOO MANY ITERATIONS))',
     &           'AT X = ',1PE15.7)
      RETURN
20000 CONTINUE
      IER = 20000
      WRITE( * ,20001) X
20001 FORMAT(' ','(SUBR.-RKFDS) TROUBLE(INAPPROPRIATE ERROR TOLERANCE)'
     &          ,' AT X = ',1PE15.7)
      RETURN
30000 CONTINUE
      IER = 30000
      WRITE( * ,30001) X
30001 FORMAT(' ','(SUBR.-RKFDS) TROUBLE(TOO SMALL STEP SIZE)',
     &           ' AT X = ',1PE15.7)
      RETURN

40000 CONTINUE

      IER=4000
      WRITE( * , * ) 'THE EQUATIONS ARE STIFF!'
      WRITE( * ,40001) X
40001 FORMAT(' ','STIFFNESS WAS DETECTED AT X = ',1PE15.7,'.')
      WRITE( * , * ) 'TO USE THE SUBROUTINE (KRSNAT OR KRSAT)'
     &,' IS RECOMMENED.'
      RETURN
      END
C

C ********************************************
      SUBROUTINE FEHL(LINDEX,N,Y0,YN,ERR)
       IMPLICIT REAL*8(A-H,O-Z)
       DIMENSION Y0(N),YN(N),ERR(N)
       DIMENSION AK1(130),AK2(130),AK3(130),AK4(130),AK5(130),AK6(130)
       DIMENSION W(130),DX(130),Y(130)
       COMMON /RP/H,Y,T /FO/DX 
       PARAMETER(ONE = 1, TWO = 2, THR = 3, TWL = 12)
       PARAMETER(AL2 = ONE / 4, AL3 = THR / 8, AL4 = TWL / 13,
     &           AL5 = ONE, AL6 = ONE / 2)
       PARAMETER(B21 = ONE / 4, B31 = THR / 32, B32 = 9.0D+0 / 32,
     &           B41 = 1932.0D+0 / 2197, B42 = -7200.0D+0 / 2197,
     &           B43 = 7296.0D+0 / 2197, B51 = 439.0D+0 / 216,
     &           B52 = -8, B53 = 3680.0D+0 / 513,
     &           B54 = -845.0D+0 / 4104, B61 = -8.0D+0 / 27,
     &           B62 = 2, B63 = -3544.0D+0 / 2565,
     &           B64 = 1859.0D+0 / 4104, B65 = -11.0D+0 / 40)
       PARAMETER(GA1 = 16.0D+0 / 135, GA3 = 6656.0D+0 / 12825,
     &           GA4 = 28561.0D+0 / 56430,
     &           GA5 = -9.0D+0 / 50, GA6 = TWO / 55)
       PARAMETER(DA1 = ONE / 360, DA3 = -128.0D+0 / 4275,
     &           DA4 = -2197.0D+0 / 75240,
     &           DA5 = ONE / 50, DA6 = TWO / 55)
C
      TT=T
      DO 5 I=1,N
        Y(I)=Y0(I)
    5 CONTINUE
      CALL SOLV(T,Y0)
      CALL FUNC(T,Y0,N)

      DO 10 I = 1,N 
        AK1(I)=DX(I)
        W(I) = Y0(I) + H * B21 * AK1(I)
   10 CONTINUE

      T=TT+AL2*H
      CALL SOLV(T,W)
      CALL FUNC(T,W,N)

      DO 20 I = 1,N
        AK2(I)=DX(I)
        W(I) = Y0(I) + H * (B31 * AK1(I) + B32 * AK2(I))
   20 CONTINUE

      T=TT+AL3*H
      CALL SOLV(T,W)
      CALL FUNC(T,W,N)

      DO 30 I = 1,N
         AK3(I)=DX(I)
         W(I) = Y0(I) + H * (B41 * AK1(I) + B42 * AK2(I) + B43 * AK3(I))
   30 CONTINUE

      T=TT+AL4*H
      CALL SOLV(T,W)
      CALL FUNC(T,W,N)

      DO 40 I = 1,N
        AK4(I)=DX(I)
        W(I) = Y0(I) + H * (B51 * AK1(I) + B52 * AK2(I)
     &                  + B53 * AK3(I) + B54 * AK4(I))
   40 CONTINUE

      T=TT+AL5*H
      CALL SOLV(T,W)
      CALL FUNC(T,W,N)

      DO 50 I = 1,N
        AK5(I)=DX(I)
        W(I) = Y0(I) + H * (B61 * AK1(I) + B62 * AK2(I) + B63 * AK3(I)
     &                  + B64 * AK4(I) + B65 * AK5(I))
   50 CONTINUE

      T=TT+AL6*H
      CALL SOLV(T,W)
      CALL FUNC(T,W,N)

      DO 60 I = 1,N
        AK6(I)=DX(I)
        YN(I) = Y(I) + H * (GA1 * AK1(I) + GA3 * AK3(I) + GA4 * AK4(I)
     &                  + GA5 * AK5(I) + GA6 * AK6(I))
       Y(I)=YN(I)
   60 CONTINUE

      T=TT+H
      CALL SOLV(T,Y)
      IF  (LINDEX .EQ. 1) THEN
        DO 70 I = 1,N
          ERR(I) = H * (DA1 * AK1(I) + DA3 * AK3(I) + DA4 * AK4(I)
     &               + DA5 * AK5(I) + DA6 * AK6(I))
   70   CONTINUE
        RETURN
      ELSE
        RETURN
      ENDIF
      END

      SUBROUTINE PARM
      DIMENSION DE(40),DF(40),PT(3),PA(999),X(130)
      REAL*8 T0,T1,TI,DE,DF,PT,PA,V,H,X,T,CONSTVALUE
	INTEGER COUNT,CONUM,QUOTIENT,ODD
      CHARACTER NA*2,NAM*5,LABEL(100)*5,LINE*72
      COMMON /BO/NS,ING /TA/T0,T1,TI /ES/DE,DF,ND /PAT/PT /PAF/PA 
     1       /RP/H,X,T /OPT/NT,NOT /OU2/NOUT /OU3/LABEL
	COMMON /CON/CONSTVALUE,QUOTIENT,ODD	
C
      NOUT=0
    1 READ(50,100) LINE
  100 FORMAT(A72)
C
      IF(LINE(1:3).NE.'END') THEN
        BACKSPACE(50)
        NA=LINE(1:2)
C
        IF(NA.EQ.'PA') THEN

          READ(50,110) NO,V
  110 FORMAT(2X,I4,2X,D15.8)
          PA(NO)=V
C
        ELSEIF(NA.EQ.'SU') THEN

          READ(50,120) NO,NAM
C
        ELSEIF(NA.EQ.'LA') THEN
          READ(50,120) NO,NAM
  120 FORMAT(2X,I3,1X,A5)
          NOUT=NOUT+1
          LABEL(NO)=NAM
C
        ELSEIF(NA.EQ.'NS') THEN
          READ(50,130) NO
  130 FORMAT(5X,I11)
          NS=NO
C
        ELSEIF(NA.EQ.'IN') THEN
          READ(50,130) NO
          ING=NO
          N=NS+ING
          DO 10 I=1,N
            X(I)=0.0D-00
   10     CONTINUE
C
        ELSEIF(NA.EQ.'ND') THEN
          READ(50,130) NO
          ND=NO
          IF(ND.GT.0) THEN
            I0=ND/2
            DO 20 I=1,I0
              DE(I)=0.0D-00
              DF(I)=0.0D-00
   20       CONTINUE
          ENDIF
C
        ELSEIF(NA.EQ.'PT') THEN
          READ(50,110) NO,V

          PT(NO)=V
C
        ELSEIF(NA.EQ.'NO') THEN
          READ(50,130) NO

          NOT=NO
C
        ELSEIF(NA.EQ.'OP') THEN
          READ(50,130) NO

          NOUT=NO
C
        ELSEIF(NA.EQ.'ST') THEN
          READ(50,110) NO,V

          X(NO)=V
C
        ELSE
          WRITE(*,140) LINE
  140 FORMAT(1H1,5X,'*** PARAMETER ERROR ***'/
     1           6X,'*** ILLEAGAL INPUT ***'/
     1           6X,'...',A72)
                      RETURN
        ENDIF
        GO TO 1
C
      ELSE
        T0=PT(1)
        T1=PT(2)
        TI=PT(3)
        T=(T1-T0)/TI
	  QUOTIENT = IDNINT(T/CONSTVALUE)
	  ODD = IDNINT(T-CONSTVALUE*REAL(QUOTIENT))

	IF(T.LE.REAL(NOT)) NOT = IDNINT(T);
        H=TI
        T=T0

  185 FORMAT(2D12.5,2I8)
C
        RETURN
      ENDIF
      END

C *********************************************
      INTEGER FUNCTION ICHEK(I,J,K)            
      REAL*8 TC,T0,T1,TI                       
      COMMON /TA/T0,T1,TI /OPT/NT,NOT /ICK/ITC 
      IF(K.EQ.0) THEN                          
        TC=((T1-T0)*DBLE(J)/DBLE(NOT))/TI      
        ITC=IDNINT(TC)                         
        K=1                                    
      ENDIF                                    
CK                                             
      IF(I.LT.NT) THEN                         
        ICHEK=MOD(I,ITC)                       
      ELSE                                     
        ICHEK=0                                
      ENDIF                                    
CK                                             
      RETURN                                   
      END                                      

C **********************************************
      SUBROUTINE INIT
      DIMENSION X(130)
      REAL*8 H,X,T
      COMMON /RP/ H,X,T
      CALL SOLV(T,X)
      RETURN
      END

C **********************************************
      SUBROUTINE OFILE
      CHARACTER NFILE*50
      LOGICAL LCHK

      NFILE='temp'
  260 FORMAT(A50)

        IL=INDEX(NFILE,' ')
        NFILE(IL:IL+3)='.csv'
        OPEN(UNIT=10,FILE=NFILE(1:IL+3), STATUS='UNKNOWN')
        NFILE(IL:IL+3)='.PAR'
        OPEN(UNIT=50,FILE=NFILE(1:IL+3),STATUS='OLD',ERR=300)
        RETURN
  300 WRITE(*,350) NFILE
  350 FORMAT(' ** FILE NAMED ',A50/
     1       '** DOES NOT EXIST. CHECK IT!')
      STOP
C
      END

C **************************
      SUBROUTINE CFILE
C
      CLOSE(50)
      CLOSE(10)
C
      RETURN
      END
