      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
C
      CALL OUT
      CALL PLO(J)
      J=J+1
      K=0
C
	COUNT = 0
	CONUM = IDNINT((CONSTVALUE*REAL(QUOTIENT)+REAL(ODD))/REAL(NOT))

	DO 155 I=1,ODD
        CALL RUNGE
        IF(CONUM.EQ.COUNT) THEN
          CALL OUT
          CALL PLO(J)
          J=J+1
          COUNT=0
        ENDIF
		COUNT=COUNT+1
  155 CONTINUE
C
	DO 158 K=1,QUOTIENT
		DO 158 I=1,IDNINT(CONSTVALUE)
			CALL RUNGE
			IF(CONUM.EQ.COUNT) THEN
				CALL OUT
				CALL PLO(J)
				J=J+1
				COUNT=0
			ENDIF
		COUNT=COUNT+1
  158 CONTINUE

      CALL CFILE
C
      STOP
      END

C*************************************************
      SUBROUTINE RUNGE
      DIMENSION X(130),X1(130),K0(130),PHI(130),DX(130)
      REAL*8 X,X1,K0,PHI,DX,T,H,CS1,CS2
      COMMON /RP/H,X,T /BO/NS,ING /FO/DX
      N=NS+ING
C
      DO 10 I=1,N
         X1(I)=X(I)
         K0(I)=0.0D-00
         DX(I)=0.0D-00
         PHI(I)=0.0D-00
   10 CONTINUE
      CS1=1.0D-00/DSQRT(2.0D-00)
      CS2=DSQRT(2.0D-00)
C
      CALL SOLV(T,X1)
      CALL FUNC(T,X1,N)
      DO 20 I=1,N
         X1(I)=X(I)+0.5D-00*H*DX(I)
         PHI(I)=PHI(I)+DX(I)
         K0(I)=DX(I)
   20 CONTINUE
C
      T=T+0.5D-00*H
      CALL SOLV(T,X1)
      CALL FUNC(T,X1,N)
      DO 30 I=1,N
         X1(I)=X(I)+(CS1-0.5D-00)*H*K0(I)+(1.0D-00-CS1)*H*DX(I)
         PHI(I)=PHI(I)+(2.0D-00-CS2)*DX(I)
         K0(I)=DX(I)
   30 CONTINUE
C
      CALL SOLV(T,X1)
      CALL FUNC(T,X1,N)
      DO 40 I=1,N
         X1(I)=X(I)-CS1*H*K0(I)+(1.0D-00+CS1)*H*DX(I)
         PHI(I)=PHI(I)+(2.0D-00+CS2)*DX(I)
   40 CONTINUE
C
      T=T+0.5D-00*H
      CALL SOLV(T,X1)
      CALL FUNC(T,X1,N)
      DO 50 I=1,N
         PHI(I)=PHI(I)+DX(I)
         X(I)=X(I)+PHI(I)*H/6.0D-00
   50 CONTINUE
      CALL SOLV(T,X)
C
      RETURN
      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)
  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 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)
      IF(LINE(1:3).NE.'END') THEN
        BACKSPACE(50)
        NA=LINE(1:2)
        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
        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
