      CHARACTER ME(4)*27,NA*20,CN(4)*37,XS*10,YS*10,YN
      DIMENSION TX(50),CY(50),A(20,21),P(20),CS(50,20),T1(25),C1(25)
      DIMENSION NL(5),X0(20),X1(20),A0(20),A1(20),OC(50)
      DIMENSION XT(50),YT(50),PL(20),PU(20),Q(20)
      COMMON NN,ND,DT,TT(2000),CC(2000,15),R(20)
      ME(1)='GAUSS-NEWTON         METHOD'
      ME(2)='DAMPING GAUSS-NEWTON METHOD'
      ME(3)='MODIFIED MARQUARDT   METHOD'
      ME(4)='SIMPLEX              METHOD'
      CN(1)='NO CONSTRAINT'
      CN(2)='P(I)=Q(I)*Q(I)'
      CN(3)='P(I)=B+(A-B)*SIN^2(Q(I))'
      CN(4)='P(I)=B+(A-B)*(EXP(Q(I))/(1+EXP(Q(I)))'
      RR=RND(0)
      NS=50
      WRITE(*,*)'#######################################'
      WRITE(*,*)'#          MULTI LINES FITTING        #'
      WRITE(*,*)'#      BY RUNGE-KUTTA-GILL METHOD     #'
      WRITE(*,*)'#######################################'
      WRITE(*,*)
      CALL FMA(N,LN,NL,TX,CY,NA)
      DO 10 I=1,4
   10 WRITE(*,'('' ('',I1,'')'',A)')I-1,ME(I)
      WRITE(*,*)
   12 WRITE(*,*)'WHICH ALGORITHM DO YOU SELECT ?'
      READ(*,*)LA
      IF(LA.NE.0.AND.LA.NE.1.AND.LA.NE.2.AND.LA.NE.3)GOTO 12
      WRITE(*,*)'NUMBER OF CONSTANTS(R) ?'
      READ(*,*)MC
      DO 14 I=1,MC
      WRITE(*,'('' R('',I2,'') ?'')') I
      READ(*,*) R(I)
   14 CONTINUE
      WRITE(*,*)'NUMBER OF DIFF.EQUATIONS ?'
      READ(*,*)ND
      WRITE(*,*)'NUMBER OF TIME COURSES(1-5) ?'
      READ(*,*)LN
      WRITE(*,*)'NUMBER OF PARAMETERS ?'
      READ(*,*)M
      NP=M
      T9=-1.E+10
      DO 15 J=1,LN
      DO 27 I=1,NL(J)
      IF(T9.LT.TX(IB+I))T9=TX(IB+I)
   27 CONTINUE
      IB=IB+NL(J)
   15 CONTINUE
   17 WRITE(*,'('' INTEGRATION INTERVAL(>'',E12.5,'') ?'')')T9/2000
      READ(*,*)DT
      IF(DT.LT.T9/2000) GOTO 17
      NN=T9/DT+1
C --EXECUTION OF FITTING JOB-
  840 WRITE(*,*)
      WRITE(*,*)'WEIGHT OF DATA(0,1,2) ?'
      READ(*,*)IW
      WRITE(*,*)' --- CONSTRAINT ON P(I) ---'
      DO 30 I=1,4
      WRITE(*,'('' ('',I1,'') '',A)')I,CN(I)
   30 CONTINUE
   32 WRITE(*,*)'WHICH CONSTRAINT DO YOU SELECT (1,2,3,4)?'
      READ(*,*)NC
      IF(NC.NE.1.AND.NC.NE.2.AND.NC.NE.3.AND.NC.NE.4)GOTO 32
      DO 45 I=1,M
      IF(NC.GT.2)THEN
      WRITE(*,'(''  LOWER, UPPER OF P('',I2,'') ?'')')I
      READ(*,*)PL(I),PU(I)
      END IF
   42 WRITE(*,'(''  INITIAL P('',I2,'') ?'')')I
      READ(*,*)P(I)
      IF(NC.GT.2.AND.(P(I).LE.PL(I).OR.P(I).GE.PU(I)))GOTO 42
   45 CONTINUE
  830 DO 47 I=1,M
      CALL CONST1(I,NC,P,Q,PL,PU)
      A1(I)=Q(I)
   47 CONTINUE
      CALL CALSS(M,SS,LN,NL,TX,CY,IW,P,NC,Q,PL,PU)
      S1=SS
      WRITE(*,'('' INITIAL SS='',E12.5)')SS
C -- JUMP TO SIMPLEX METHOD -
      IF(LA.EQ.3)GOTO 3000
C --GAUSS NEWTON AND MODIFIED MARQUARDT METHOD-
      PC=0.0001
      CF=100.0
      DP=0.001
      DO 65 KK=1,100
      CALL JACOBI(P,NC,Q,PL,PU,DP,CS,M,LN,TX,NL)
      DO 180 I=1,M
      DO 180 J=I,M
      A(I,J)=0.0
      DO 190 L=1,N
      A(I,J)=A(I,J)+CS(L,I)*CS(L,J)/CY(L)**IW
  190 CONTINUE
      A(J,I)=A(I,J)
  180 CONTINUE
      IB=0
      DO 192 J=1,LN
      DO 194 I=1,NL(J)
      T=TX(IB+I)
      CALL TCP(J,T,CP,K)
      OC(IB+I)=CY(IB+I)-CP
  194 CONTINUE
      IB=IB+NL(J)
  192 CONTINUE
      DO 112 I=1,M
      A(I,M+1)=0.0
      DO 112 J=1,N
      A(I,M+1)=A(I,M+1)+CS(J,I)*OC(J)/CY(J)**IW
  112 CONTINUE
      NP=M
      IF (LA.EQ.2)THEN
      DO 120 I=1,M
      A(I,I)=A(I,I)+CF
      X1(I)=A(I,M+1)
  120 CONTINUE
      END IF
      CALL CHORES(A,NP)
      JJ=0
   95 JJ=JJ+1
      IF(JJ.GT.25)THEN
      WRITE(*,*)'OUT OF DAMPING'
      GOTO 730
      END IF
      DO 70 I=1,M
      Q(I)=A1(I)+A(I,M+1)
   70 CONTINUE
      CALL CALSS(M,SS,LN,NL,TX,CY,IW,P,NC,Q,PL,PU)
      DS=ABS(S1-SS)
      IF((LA.NE.2).OR.(SS.EQ.0.0))GOTO 75
C --FLECHER MODIFICATION-
      PW=0.0
      DO 80 I=1,M
      PW=PW+X1(I)*A(I,M+1)+CF*A(I,M+1)*A(I,M+1)
   80 CONTINUE
      IF(DS.GT.0.75*PW)CF=CF/2.0
      IF(DS.LT.0.25*PW)CF=5.0*CF
   75 IF(DS.LE.PC*S1)GOTO 730
C --DAMPING-
      IF((LA.EQ.1).AND.(SS.GT.S1))THEN
      DO 90 I=1,M
      A(I,M+1)=0.5*A(I,M+1)
   90 CONTINUE
      GOTO 95
      END IF
      DO 100 I=1,M
      A1(I)=Q(I)
  100 CONTINUE
      S1=SS
      WRITE(*,*)
      WRITE(*,'('' LOOP='',I4)')KK
      IF(LA.EQ.1)WRITE(*,'('' DAMP='',I2)')JJ
      IF(LA.EQ.2)WRITE(*,'('' FACTOR='',F7.3)')CF
      WRITE(*,'('' P('',I2,'')='',F12.6)')(I,P(I),I=1,M)
      WRITE(*,'('' SS='',E12.5)')SS
   65 CONTINUE
C --OUTPUT-
  730 WRITE(*,*)
      WRITE(*,'('' **** '',A,'' ****'')')NA
      WRITE(*,*)'BY ',ME(LA+1)
      WRITE(*,'('' WEIGHT =1/CP^('',I1,'')'')')IW
      WRITE(*,'('' CONSTRAINT ON P(I): '',A)')CN(NC)
      WRITE(*,'('' LOOP = '',I4)')KK
      IF ((LA.NE.3).AND.(N.GT.M))THEN
      CALL VARIAN(X1,A,NP,NL,TX,CS,CY,LN,IW,M,N,P,NC,Q,PL,PU,DP)
      END IF
      IF(SS.EQ.0.0)THEN
      WRITE(*,*)'AIC=-INFINITE'
      ELSE
      WRITE(*,'('' AIC = '',F12.6)')N*ALOG(SS)+2.*M
      END IF
      WRITE(*,*)
      DO 110 I=1,M
      IF((LA.NE.3).AND.(X1(I).GT.0.0).AND.(N.GT.M))THEN
      IF (NC.EQ.1)SD=1.0
      IF (NC.EQ.2)SD=2.0*Q(I)
      IF (NC.EQ.3)SD=2.0*(PL(I)-PU(I))*SIN(Q(I))*COS(Q(I))
      IF (NC.EQ.4)SD=(PL(I)-PU(I))*EXP(Q(I))/(1.0+EXP(Q(I)))**2
      WRITE(*,508)I,P(I),ABS(SD)*SQRT(X1(I)*SS/FLOAT(N-M))
  508 FORMAT(' FILAL P(',I2,')=',F12.6,'   S.D.=',F12.6)
      ELSE
      WRITE(*,'('' FINAL P('',I2,'')='',F12.6)')I,P(I)
      END IF
  110 CONTINUE
      WRITE(*,'('' SS='',E12.5)')SS
      WRITE(*,*)
      CALL TMP
      IB=0
      DO 115 J=1,LN
      Y9=-1.E+10
      WRITE(*,'('' ## TIME COURSE'',I2,'' ##'')')J
      NE=NL(J)
      DO 122 I=1,NE
      T=TX(IB+I)
      CALL TCP(J,T,CP,K)
      WRITE(*,512)J,T,J,CP,CY(IB+I)
  512 FORMAT(' T',I2,'=',F12.6,' ,CP',I2,' = ',F12.6,'(',F12.6,')')
      IF(Y9.LT.CY(IB+I))Y9=CY(IB+I)
      T1(I)=T
      C1(I)=CY(IB+I)
  122 CONTINUE
      D1=T9/NS
      DO 800 I=1,NS
      XT(I)=D1*(I-1)
      T=XT(I)
      CALL TCP(J,T,CP,K)
      YT(I)=CP
      IF(Y9.LT.YT(I))Y9=YT(I)
  800 CONTINUE
      WRITE(*,*)'NAME OF T-AXIS ?'
      READ(*,'(A)')XS
      WRITE(*,*)'NAME OF CP-AXIS ?'
      READ(*,'(A)')YS
      CALL PLOT(NE,T1,C1,NS,XT,YT,0.,T9,0.,Y9,XS,YS)
      WRITE(*,*) 'SAVE(Y/N) ?'
      READ(*,'(A)') YN
      IF(YN.EQ.'Y') CALL SAVE(NE,T1,C1,NS,XT,YT,M,P,MC,R,N,SS)
      IB=IB+NL(J)
  115 CONTINUE
C--SELECTION OF ALGORITHM-
      DO 140 I=1,4
      WRITE(*,'('' ( '',I2,'') '',A)')I-1,ME(I)
  140 CONTINUE
      WRITE(*,*)'( -1) END'
      WRITE(*,*)
  142 WRITE(*,*)' WHICH ALGORITHM DO YOU SELECT (0-3, OR -1) ? '
      READ(*,*)LA
      IF(LA.LT.0)STOP
      IF(LA.NE.0.AND.LA.NE.1.AND.LA.NE.2.AND.LA.NE.3)GOTO 142
      WRITE(*,*) 'RE-INPUT INITIAL VALUES (Y/N) ?'
      READ(*,'(A)')YN
      IF(YN.EQ.'Y')GOTO 840
      GOTO 830
C --SIMPLEX METHOD-
 3000 AA=1.0
      BB=0.5
      CD=2.0
      SG=1.0E+10
      PC=0.00001
      DO 200 I=1,M
      A(I,1)=Q(I)
  200 CONTINUE
      DO 210 J=2,M+1
      DO 210 I=1,M
      A(I,J)=2.0*RND(1)*A(I,1)+0.01*(RND(1)-0.5)
  210 CONTINUE
      DO 220 K=1,M+1
      DO 230 I=1,M
      Q(I)=A(I,K)
  230 CONTINUE
      CALL CALSS(M,SS,LN,NL,TX,CY,IW,P,NC,Q,PL,PU)
      A0(K)=SS
  220 CONTINUE
      KK=0
 1000 KK=KK+1
      IF(KK.GE.500)GOTO 730
      WRITE(*,*)
      WRITE(*,'('' LOOP='',I4)')KK
      DO 225 I=1,M
      WRITE(*,'('' P('',I2,'')='',F10.5)')I,P(I)
  225 CONTINUE
      WRITE(*,'('' SS='',E12.5)')SS
      SR=0.0
      DO 235 I=1,M+1
      SR=SR+A0(I)
  235 CONTINUE
      IF (ABS(SR-SG).LT.PC*SG)THEN
      DO 236 I=1,M
      Q(I)=A(I,JL)
  236 CONTINUE
      SS=A0(JL)
      GOTO 730
      END IF
      SG=SR
      SR=0.0
      SL=1.0E+10
      DO 255 J=1,M+1
      IF(SR.LT.A0(J))THEN
      JH=J
      SR=A0(J)
      END IF
      IF(SL.GT.A0(J))THEN
      JL=J
      SL=A0(J)
      END IF
  255 CONTINUE
      SR=0.0
      DO 260 J=1,M+1
      IF((J.NE.JH).AND.(SR.LT.A0(J)))THEN
      JS=J
      SR=A0(J)
      END IF
  260 CONTINUE
      DO 265 I=1,M
      X0(I)=0.0
      DO 270 J=1,M+1
      IF (J.NE.JH)X0(I)=X0(I)+A(I,J)
  270 CONTINUE
      X0(I)=X0(I)/M
  265 CONTINUE
      DO 275 I=1,M
      A1(I)=(1.0+AA)*X0(I)-AA*A(I,JH)
      Q(I)=A1(I)
  275 CONTINUE
      CALL CALSS(M,SS,LN,NL,TX,CY,IW,P,NC,Q,PL,PU)
      SR=SS
      IF(SR.LE.A0(JS))GOTO 1500
      IF(SR.LT.A0(JH))THEN
      DO 285 I=1,M
      A(I,JH)=A1(I)
  285 CONTINUE
      A0(JH)=SR
      END IF
      DO 290 I=1,M
      A1(I)=BB*A(I,JH)+(1.0-BB)*X0(I)
      Q(I)=A1(I)
  290 CONTINUE
      CALL CALSS(M,SS,LN,NL,TX,CY,IW,P,NC,Q,PL,PU)
      SR=SS
      IF(SR.LE.A0(JH))THEN
      DO 295 I=1,M
      A(I,JH)=A1(I)
  295 CONTINUE
      A0(JH)=SR
      ELSE
      DO 300 K=1,M+1
      DO 305 I=1,M
      A(I,K)=(A(I,K)+A(I,JL))/2.0
      Q(I)=A(I,K)
  305 CONTINUE
      CALL CALSS(M,SS,LN,NL,TX,CY,IW,P,NC,Q,PL,PU)
      A0(K)=SS
  300 CONTINUE
      END IF
      GOTO 1000
 1500 IF(SR.LT.A0(JL))THEN
      DO 320 I=1,M
      X1(I)=CD*A1(I)+(1.0-CD)*X0(I)
      Q(I)=X1(I)
  320 CONTINUE
      CALL CALSS(M,SS,LN,NL,TX,CY,IW,P,NC,Q,PL,PU)
      SL=SS
      IF(SL.LT.SR)THEN
      DO 325 I=1,M
      A(I,JH)=X1(I)
  325 CONTINUE
      A0(JH)=SL
      GOTO 1000
      END IF
      END IF
      DO 315 I=1,M
      A(I,JH)=A1(I)
  315 CONTINUE
      A0(JH)=SR
      GOTO 1000
      END
C=====  SUBROUTINES  =====
C ----  MODIFIED CHORESKY METHOD ---
      SUBROUTINE CHORES(A,NP)
      DIMENSION A(20,21),SL(20,21)
      IF (NP.EQ.1)THEN
      A(1,2)=A(1,2)/A(1,1)
      RETURN
      END IF
      DO 10 JS=1,NP
      SL(1,JS)=A(1,JS)
   10 CONTINUE
      DO 20 IS=2,NP
      DO 30 JS=IS,NP
      DS=0
      DO 40 KS=1,IS-1
      DS=DS+SL(KS,IS)*SL(KS,JS)/SL(KS,KS)
   40 CONTINUE
      SL(IS,JS)=A(IS,JS)-DS
   30 CONTINUE
   20 CONTINUE
      SL(1,NP+1)=A(1,NP+1)
      DO 50 IS=2,NP
      DS=0
      DO 60 KS=1,IS-1
      DS=DS+SL(KS,IS)*SL(KS,NP+1)/SL(KS,KS)
   60 CONTINUE
      SL(IS,NP+1)=A(IS,NP+1)-DS
   50 CONTINUE
      A(NP,NP+1)=SL(NP,NP+1)/SL(NP,NP)
      DO 70 IS=NP-1,1,-1
      DS=0
      DO 80 KS=IS+1,NP
      DS=DS+SL(IS,KS)*A(KS,NP+1)
   80 CONTINUE
      A(IS,NP+1)=(SL(IS,NP+1)-DS)/SL(IS,IS)
   70 CONTINUE
      END
C --CALCULATION OF SS-
      SUBROUTINE CALSS(M,SS,LN,NL,TX,CY,IW,P,NC,Q,PL,PU)
      DIMENSION NL(5),TX(50),CY(50),P(20),Q(20),PL(20),PU(20)
      CALL RUNGE(M,P,NC,Q,PL,PU)
      SS=0.
      IB=0
      DO 10 J=1,LN
      DO 20 I=1,NL(J)
      T=TX(IB+I)
      CALL TCP(J,T,CP,K)
      SS=SS+(CY(IB+I)-CP)*(CY(IB+I)-CP)/CY(IB+I)**IW
   20 CONTINUE
      IB=IB+NL(J)
   10 CONTINUE
      END
C --JACOBIAN-
      SUBROUTINE JACOBI(P,NC,Q,PL,PU,DP,CS,M,LN,TX,NL)
      DIMENSION P(20),Q(20),CS(50,20),PL(20),PU(20),CT(2000,15)
      DIMENSION NL(5),TX(50)
      COMMON NN,ND,DT,TT(2000),CC(2000,15),R(20)
      DO 10 JS=1,M
      QT=P(JS)
      PT=Q(JS)
      Q(JS)=PT+DP
      CALL RUNGE(M,P,NC,Q,PL,PU)
      DO 20 J=1,NN
      DO 20 I=1,LN
      CT(J,I)=CC(J,I)
   20 CONTINUE
      Q(JS)=PT-DP
      CALL RUNGE(M,P,NC,Q,PL,PU)
      IB=0
      DO 30 J=1,LN
      DO 40 I=1,NL(J)
      T=TX(IB+I)
      CALL TCP(J,T,CP,K)
      CS(IB+I,JS)=(CT(K-1,J)-CP)/2./DP
   40 CONTINUE
      IB=IB+NL(J)
   30 CONTINUE
      P(JS)=QT
      Q(JS)=PT
   10 CONTINUE
      END
C --VARIANCE-
      SUBROUTINE VARIAN(X1,A,NP,NL,TX,CS,CY,LN,IW,M,N,P,NC,Q,
     +PL,PU,DP)
      DIMENSION NL(5),X1(20),A(20,21),TX(50),CS(50,20)
     +,CY(50),P(20),Q(20),PL(20),PU(20)
      CALL JACOBI(P,NC,Q,PL,PU,DP,CS,M,LN,TX,NL)
      DO 180 I=1,M
      DO 180 J=1,M
      A(I,J)=0.0
      DO 190 L=1,N
      A(I,J)=A(I,J)+CS(L,I)*CS(L,J)/CY(L)**IW
  190 CONTINUE
      A(J,I)=A(I,J)
  180 CONTINUE
      DO 10 K=1,M
      DO 20 I=1,M
      A(I,M+1)=0
   20 CONTINUE
      A(K,M+1)=1
      CALL CHORES(A,NP)
      X1(K)=A(K,M+1)
   10 CONTINUE
      END
C --RANDOM NUMBER GENERATION-
      FUNCTION RND(IN)
      IF(IN.EQ.1)GOTO 10
      X=50249347.
      A=23.
      FM=(1.E+8)+1
   10 W=X*A
      V=AMOD(W,FM)
      RND=V/FM
      X=V
      END
C =====GRAPH SUBROUTINES=====
      SUBROUTINE PLOT(NE,TX,CY,NS,XT,YT,X0,X9,Y0,Y9,XS,YS)
      DIMENSION TX(NE),CY(NE),XT(NS),YT(NS),XG(6)
      CHARACTER XS*10,YS*10,IP(51,17)
      XM=50.
      YM=16.
      DO 100 I=1,51
      DO 100 J=1,17
      IP(I,J)=' '
  100 CONTINUE
      DO 110 I=1,NS
      IX=INT(XM*(XT(I)-X0)/(X9-X0)+0.5)+1
      IY=INT(YM*(Y9-YT(I))/(Y9-Y0)+0.5)+1
      IP(IX,IY)='*'
  110 CONTINUE
      DO 120 I=1,NE
      IX=INT(XM*(TX(I)-X0)/(X9-X0)+0.5)+1
      IY=INT(YM*(Y9-CY(I))/(Y9-Y0)+0.5)+1
      IP(IX,IY)='O'
  120 CONTINUE
      WRITE(*,'(5X,A10)')YS
      DO 130 K=1,17
      JJ=K-4*INT(K/4)
      IF(JJ.EQ.1)THEN
      YY=(Y9-Y0)*(16-K+1)/16.0+Y0
      IF((YY.GE.10000.0).OR.((YY.LT.0.1).AND.(YY.NE.0.0)))THEN
      WRITE(*,'(1X,E12.4,''+'',51A1)')YY,(IP(I,K),I=1,51)
      ELSE
      WRITE(*,'(5X,F8.3,''+'',51A1)')YY,(IP(I,K),I=1,51)
      ENDIF
      ELSE
      WRITE(*,'(13X,''I'',51A1)')(IP(I,K),I=1,51)
      ENDIF
  130 CONTINUE
      WRITE(*,'(13X,''++'',5A10)')('---------+',I=1,5)
      DO 140 I=1,6
      XG(I)=(X9-X0)*(I-1)/5.0+X0
  140 CONTINUE
      WRITE(*,'(7X,6F10.2)')(XG(I),I=1,6)
      WRITE(*,'(35X,A10)')XS
      END
      SUBROUTINE TMP
      CHARACTER TP
      WRITE(*,*)'DO YOU CONTINUE(Y/N)?'
      READ(*,'(A1)')TP
      IF(TP.EQ.'N')STOP
      END
C -- RUNGE-KUTTA-GILL METHOD -
      SUBROUTINE RUNGE(M,P,NC,Q,PL,PU)
      DIMENSION P(20),Q(20),PL(20),PU(20)
      DIMENSION DC(15),C0(15),C(15),QQ(15),RR(15),AK(15)
      COMMON NN,ND,DT,TT(2000),CC(2000,15),R(20)
      SQ=SQRT(.5)
      DO 10 I=1,M
      CALL CONST2(I,NC,P,Q,PL,PU)
   10 CONTINUE
      CALL INI(C0,P)
      DO 110 I=1,ND
      QQ(I)=0.
      C(I)=C0(I)
  110 CONTINUE
      DO 150 J=1,NN
      T=DT*(J-1)
      TT(J)=T
      DO 155 I=1,ND
      CC(J,I)=C(I)
  155 CONTINUE
      CALL DIFF(DC,C,T,P)
      DO 120 I=1,ND
      AK(I)=DT*DC(I)
      RR(I)=(AK(I)-2.*QQ(I))/2.
      C(I)=C(I)+RR(I)
      QQ(I)=QQ(I)+3.*RR(I)-AK(I)/2.
  120 CONTINUE
      CALL DIFF(DC,C,T,P)
      DO 130 I=1,ND
      AK(I)=DT*DC(I)
      RR(I)=(1.-SQ)*(AK(I)-QQ(I))
      C(I)=C(I)+RR(I)
      QQ(I)=QQ(I)+3.*RR(I)-(1.-SQ)*AK(I)
  130 CONTINUE
      CALL DIFF(DC,C,T,P)
      DO 140 I=1,ND
      AK(I)=DT*DC(I)
      RR(I)=(1.+SQ)*(AK(I)-QQ(I))
      C(I)=C(I)+RR(I)
      QQ(I)=QQ(I)+3.*RR(I)-(1.+SQ)*AK(I)
  140 CONTINUE
      CALL DIFF(DC,C,T,P)
      DO 150 I=1,ND
      AK(I)=DT*DC(I)
      RR(I)=(AK(I)-2.*QQ(I))/6.
      C(I)=C(I)+RR(I)
      QQ(I)=QQ(I)+3.*RR(I)-AK(I)/2.
  150 CONTINUE
      END
      SUBROUTINE TCP(J,T,CP,K)
      COMMON NN,ND,DT,TT(2000),CC(2000,15),R(20)
      SP=1.E+10
      DO 10 K=1,NN
      S2=(TT(K)-T)*(TT(K)-T)
      IF(SP.LE.S2)GOTO 20
      SP=S2
   10 CONTINUE
   20 CP=CC(K-1,J)
      END
C=====CONSTRAINTS=====
      SUBROUTINE CONST1(I,NC,P,Q,PL,PU)
      DIMENSION P(20),Q(20),PL(20),PU(20)
      GOTO (10,20,30,40)NC
   10 Q(I)=P(I)
      RETURN
   20 Q(I)=SQRT(P(I))
      RETURN
   30 TM=SQRT((P(I)-PU(I))/(PL(I)-PU(I)))
      Q(I)=ATAN(TM/SQRT(1-TM*TM))
      RETURN
   40 Q(I)=ALOG((P(I)-PU(I))/(PL(I)-P(I)))
      RETURN
      END
      SUBROUTINE CONST2(I,NC,P,Q,PL,PU)
      DIMENSION P(20),Q(20),PL(20),PU(20)
      GOTO(10,20,30,40)NC
   10 P(I)=Q(I)
      RETURN
   20 P(I)=Q(I)*Q(I)
      RETURN
   30 P(I)=PU(I)+(PL(I)-PU(I))*SIN(Q(I))*SIN(Q(I))
      RETURN
   40 P(I)=PU(I)+(PL(I)-PU(I))*EXP(Q(I))/(1+EXP(Q(I)))
      RETURN
      END
C   ====== FILE I/O ======= 
      SUBROUTINE FMA(N,LN,NL,X,Y,FN)
      DIMENSION NL(5),X(50),Y(50)
      CHARACTER FN*10,IC*1
  100 WRITE(*,*) 'E)XIT, I)NPUT, R)EAD, W)RITE, S)TOP ?'
      READ(*,'(A)') IC
      IF(IC.EQ.'S') STOP
      IF(IC.EQ.'E'.AND.N.GT.0) RETURN
      IF(IC.EQ.'I') GOTO 1100
      IF(IC.EQ.'R') GOTO 1200
      IF(IC.EQ.'W') GOTO 1300
      GOTO 100
 1100 WRITE(*,*) 'NUMBER OF LINES (1-5) ?'
      READ(*,*) LN
      N=0
      IB=0
      DO 1115 J=1,LN
      WRITE(*,'('' NUMBER OF POINTS ('',I2,'') ?'')') J
      READ(*,*) NL(J)
      N=N+NL(J)
      DO 1127 I=1,NL(J)
      WRITE(*,1125) J,I,J,I
 1125 FORMAT(' T',I2,'(',I2,')  , CP',I2,'(',I2,') ?')
      READ(*,*) X(IB+I),Y(IB+I)
 1127 CONTINUE
      IB=IB+NL(J)
 1115 CONTINUE 
      GOTO 100
 1200 WRITE(*,*) 'FILE NAME TO READ ?'
      READ(*,'(A)') FN
      OPEN(UNIT=5,FILE=FN)
      READ(5,*) LN
      N=0
      IB=0
      DO 1220 J=1,LN
      READ(5,*) NL(J)
      N=N+NL(J)
      WRITE(*,*)
      DO 1210 I=1,NL(J)
      READ(5,*) X(IB+I),Y(IB+I)
      WRITE(*,1225) J,I,X(IB+I),J,I,Y(IB+I)
 1225 FORMAT(' T',I2,'(',I2,')=',F10.5,'  , CP',I2,'(',I2,')=',F10.5)
 1210 CONTINUE
      IB=IB+NL(J)
 1220 CONTINUE
      CLOSE(UNIT=5)
      GOTO 100
 1300 WRITE(*,*) 'FILE NAME TO WRITE ?'
      READ(*,'(A)') FN
      OPEN(UNIT=6,FILE=FN,STATUS='NEW')
      WRITE(6,*) LN
      N=0
      IB=0
      DO 1315 J=1,LN
      WRITE(6,*) NL(J)
      N=N+NL(J)
      DO 1327 I=1,NL(J)
      WRITE(6,'(F12.5,'','',F12.5)') X(IB+I),Y(IB+I)
 1327 CONTINUE
      IB=IB+NL(J)
 1315 CONTINUE 
      CLOSE(UNIT=6)
      GOTO 100
      END
C  ======  GRAPHIC DATA SAVE =======
      SUBROUTINE SAVE(NE,TX,CY,NS,XT,YT,M,P,MC,R,N,SS)
      DIMENSION TX(50),CY(50),XT(50),YT(50),P(20),R(20)
      CHARACTER NA*12
      WRITE(*,*)'FILE NAME ?'
      READ(*,'(A)') NA
      OPEN(UNIT=6,FILE=NA,STATUS='NEW')
      WRITE(6,*) NE
      DO 10 I=1,NE
      WRITE(6,'(F12.5,'','',F12.5)') TX(I),CY(I)
   10 CONTINUE
      WRITE(6,*) NS
      DO 20 I=1,NS
      WRITE(6,'(F12.5,'','',F12.5)') XT(I),YT(I)
   20 CONTINUE
      DO 30 I=1,M
      WRITE(6,'(''P('',I2,'')=,'',F12.5)')I,P(I)
   30 CONTINUE
      DO 40 I=1,MC
      WRITE(6,'(''R('',I2,'')=,'',F12.5)')I,R(I)
   40 CONTINUE
      IF(SS.GT.0)THEN
      WRITE(6,'(''AIC='',F12.5)')N*ALOG(SS)+2.*M
      ENDIF
      CLOSE(UNIT=6)
      END
C  =====  DIFF. EQUATIONS ======
      SUBROUTINE DIFF(DC,C,T,P)
      DIMENSION C0(15),DC(15),C(15),P(20)
      COMMON NN,ND,DT,TT(2000),CC(2000,15),R(20)
      DC(1)=-P(2)*C(1)
      DC(2)=-P(2)*C(2)+P(1)*P(4)*P(3)*EXP(-P(4)*T)
      END
      SUBROUTINE INI(C0,P)
      DIMENSION C0(15),P(20)
      COMMON NN,ND,DT,TT(2000),CC(2000,15),R(20)
      C0(1)=P(1)
      C0(2)=0.
      END
