      SUBROUTINE OPTION (MAXSP,MAXEQ,RW,NRW,IW,NIW,KFLAG)
      REAL RW(NRW)
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT,INFDAT
      INTEGER*2 IW(NIW)
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
      COMMON / INFORM / IHELP,INFDAT
      KFIN=KFLAG
      IF(KFLAG.LE.0) RETURN
      IF(KFLAG.NE.1) GOTO 5
1     WRITE(ITOUT,8000)
      WRITE(ITOUT,8001)
      WRITE(ITOUT,8002)
      WRITE(ITOUT,8003)
8000  FORMAT(//,' Choose option:',
     @ /,25H  1  := list this options ,
     @ /,53H  2  := pattern of Jacobian on terminal              ,
     @ /,53H  3  := pattern of Jacobian on printer               ,
     @ /,53H  4  := pattern of Jacobian on plotter               ,
     @ /,53H  5  := pattern of Jacobian on plotter (dialog)      )
8001  FORMAT(
     @   53H  6  := pattern of stoichiometric matrix on terminal ,
     @ /,53H  7  := pattern of stoichiometric matrix on printer  ,
     @ /,53H  8  := pattern of stoichiometric matrix on plotter  ,
     @ /,52H  9  := pattern of stoichiometric matrix on plotter ,
     @  8H(dialog))
8002  FORMAT(
     @   57H 10  := generation of explicit subroutines FCN and JACOBI)
8003  FORMAT(48H 11  := generation of explicit subroutine STRUCT,
     @ /,57H 12  := generation of FCN and JACOBI in full storage mode,
     @ /,57H 13  := list papers related to LARKIN                    ,
     @ /,57H 14  := print papers related to LARKIN                   ,
     @ /,30H  0  := leave subsystem OPTION,/)
      READ(ITIN,*) KFLAG
5     IF(KFLAG.LE.0) GOTO 100
      IF(KFLAG.EQ.1) GOTO 1
      IF(KFLAG.EQ.2) GOTO 10
      IF(KFLAG.EQ.3) GOTO 20
      IF(KFLAG.EQ.4) GOTO 30
      IF(KFLAG.EQ.5) GOTO 30
      IF(KFLAG.GT.5 .AND. KFLAG.LT.10) GOTO 40
      IF(KFLAG.EQ.10) GOTO 50
      IF(KFLAG.EQ.11) GOTO 50
      IF(KFLAG.EQ.12) GOTO 50
      IF(KFLAG.EQ.13) GOTO 60
      IF(KFLAG.EQ.14) GOTO 60
      IF(KFLAG.GT.14) GOTO 100
      GOTO 100
10    LOUT=ITOUT
      ITYPE=1
      GOTO 21
20    LOUT=POUT
      ITYPE=2
21    LIN=SIN
      ITYPE=KFLAG-1
      NH=MAXSP+3
      NH1=NH+1
      NH2=NH1+NH
      NH3=NH2+NH
      NH4=NH3+NH
      NIWH=NIW-NH4+1
      CALL SPAP1 (LIN,LOUT,MAXSP,MAXEQ,IW(NH4),NIWH,IW(1),IW(NH1),
     @            IW(NH2),IW(NH3),NH,ITYPE)
      IF(KFIN.EQ.1) GOTO 1
      GOTO 100
30    LIN=SIN
      LOUT=ITOUT
      ITYPE=KFLAG-1
      NH=MAXSP+3
      NH1=1
      NH2=NH1+NH
      NH3=NH2+NH
      NH4=NH3+NH
      NIWH=NIW-NH4+1
      NHR=NRW/3-2
      NHR2=NHR+1
      NHR3=NHR2+NHR
      CALL SPAP2 (LIN,LOUT,MAXSP,MAXEQ,IW(NH4),NIWH,IW(NH1),
     @            IW(NH2),IW(NH3),NH,ITYPE,RW(1),RW(NHR2),RW(NHR3),NHR)
      IF(KFIN.EQ.1) GOTO 1
      GOTO 100
40    LIN=SIN
      LOUT=ITOUT
      IF(KFLAG.EQ.7) LOUT=POUT
      ITYPE=KFLAG-5
      IF(ITYPE.GT.2) GOTO 45
      NH0=MAXSP+3
      IF(MAXEQ+3.GT.NH0) NH0=MAXEQ+3
      NH1=MAXEQ+NH0
      NH2=NH1+MAXEQ
      NH=NIW-NH2
      CALL STOI1 (LIN,LOUT,MAXSP,MAXEQ,IW(1),IW(NH0),IW(NH1),IW(NH2),
     1            NH,ITYPE)
      IF(KFIN.EQ.1) GOTO 1
      GOTO 100
45    NH0=1
      NH1=MAXEQ+NH0
      NH2=NH1+MAXEQ
      NH=NIW-NH2
      NHR=NRW/3-2
      NHR2=NHR+1
      NHR3=NHR2+NHR
      CALL STOI2 (LIN,LOUT,MAXSP,MAXEQ,IW(NH0),IW(NH1),IW(NH2),
     1            NH,ITYPE,RW(1),RW(NHR2),RW(NHR3),NHR)
      IF(KFIN.EQ.1) GOTO 1
      GOTO 100
50    LIN=SIN
      LOUT=ITOUT
      LOUT1=HDS
      LOUT2=DOUT
      ITYPE=KFLAG-5
      NLRH=4*MAXEQ
      N1=NLRH+1
      NCEQH2=2*MAXEQ
      N2=N1+NCEQH2
      NSPH1=MAXSP+1
      N3=N2+NSPH1
      NSPH=MAXSP
      N4=N3+NSPH
      LDH=NIW-N4+1
      IF(LDH.LE.0) GOTO 510
      CALL EXPLI (LIN,LOUT,LOUT1,LOUT2,MAXSP,MAXEQ,IW(1),IW(N1),IW(N2),
     @            IW(N3),IW(N4),NLRH,NCEQH2,NSPH1,NSPH,LDH,ITYPE)
      IF(KFIN.EQ.1) GOTO 1
      GOTO 100
60    CONTINUE
      ITINS=ITIN
      IF(KFLAG.EQ.14) ITIN=POUT
      CALL HELP (23)
      ITIN=ITINS
      IF(KFIN.EQ.1) GOTO 1
      GOTO 100
100   RETURN
510   WRITE(ITOUT,9510)
9510  FORMAT(/,41H INTEGER work space exhausted; enlarge IW,/)
      IF(KFIN.EQ.1) GOTO 1
      RETURN
      END
      SUBROUTINE SPAP1 (LIN,LOUT,MAXSP,MAXEQ,JA,NIW,ROW,IA,IR,IC,
     @                  NH,ITYPE)
      INTEGER LIN,LOUT
      INTEGER*2 JA(NIW),IA(NH),ROW(NH),IR(NH),IC(NH)
      INTEGER*2 STAR,BLANK,SIGN,RW,RS
      DATA  STAR/1H*/ , BLANK/1H / , RW/1H-/ , RS/1H^/
      DATA  CHART/4H*TEM/
      REWIND LIN
      MAXP=120
      MAXT=76
      READ(LIN,9000) IHEAD
      IHEAD=IHEAD*36
      IF(IHEAD.EQ.0) GOTO 5
      READ(LIN,9005) (JA(L),L=1,IHEAD)
9000  FORMAT(19X,I5)
9005  FORMAT(36A2)
      WRITE(LOUT,9003)
9003  FORMAT(//,19X,14HActual System:,/)
      WRITE(LOUT,9005) (JA(L),L=1,IHEAD)
5     READ(LIN,9000) IEL
      IF(IEL.EQ.0) GOTO 10
      DO 9 I=1,IEL
9     READ(LIN,9001) CHAR
9001  FORMAT(A4)
10    READ(LIN,9000) NSP
      IF(NSP.EQ.0) GOTO 12
      IF(NSP.GE.MAXSP) GOTO 510
      DO 11 I=1,10000
      READ(LIN,9001) CHAR
11    IF(CHAR.EQ.CHART) GOTO 33
12    READ(LIN,9001) CHAR
33    READ(LIN,9001) CHAR
      READ(LIN,9000) NCEQ
      IF(NCEQ.EQ.0) GOTO 14
      IF(NCEQ.GT.MAXEQ) GOTO 520
      DO 13 I=1,NCEQ
13    READ(LIN,9001) CHAR
14    READ(LIN,9000) NLR
      IF(NCEQ.GT.NIW.OR.NLR.GT.NIW) GOTO 530
      READ(LIN,9500) (JA(J),J=1,NCEQ)
      NCEQ1=NCEQ+1
      NCEQ2=2*NCEQ
      READ(LIN,9500) (JA(J),J=NCEQ1,NCEQ2)
      READ(LIN,9500) (JA(J),J=1,NLR)
      READ(LIN,9000) LDIM
      READ(LIN,9500) (JA(J),J=1,LDIM)
      IF(LDIM.GT.NIW) GOTO 530
      NSP1=NSP+1
      READ(LIN,9500) (IA(J),J=1,NSP1)
      READ(LIN,9500) (IR(J),J=1,NSP)
9500  FORMAT(9I8)
      IF(ITYPE.EQ.1.AND.NSP.GT.MAXT) GOTO 550
      IF(ITYPE.EQ.2.AND.NSP.GT.MAXP) GOTO 550
      N=NSP
      N2=N+2
      DO 80 I=1,N2
      IR(I)=I
      IC(I)=I
80    CONTINUE
      WRITE(LOUT,1000)
      DO 100 L=1,N2
100   ROW(L)=RW
      WRITE(LOUT,1001) (ROW(L),L=1,N2)
      DO 200 I=1,N
      DO 222 K=1,N
222   ROW(K+1)=BLANK
      ROW(1)=RS
      IMIN=IA(IR(I))
      IMAX=IA(IR(I)+1)-1
      DO 202 J=IMIN,IMAX
      SIGN=STAR
      ROW(IC(JA(J))+1)=SIGN
202   CONTINUE
      ROW(N+2)=RS
      WRITE(LOUT,1001) (ROW(L),L=1,N2)
200   CONTINUE
      DO 300 L=1,N2
300   ROW(L)=RW
      WRITE(LOUT,1001) (ROW(L),L=1,N2)
      REWIND LIN
      RETURN
510   WRITE(LOUT,8110) NSP,MAXSP
8110  FORMAT(/,52H Actual dimension NSP greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
520   WRITE(LOUT,8120) NCEQ,MAXEQ
8120  FORMAT(/,53H Actual dimension NCEQ greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
530   WRITE(LOUT,8130)
8130  FORMAT(/,41H INTEGER work space exhausted; enlarge IW)
      GOTO 599
550   WRITE(LOUT,8150)
8150  FORMAT(/,44H Problem size too large for this output type)
599   REWIND LIN
      ITYPE=-1
      RETURN
1000  FORMAT(1H ,//,10X,22HStructure of Jacobian:,/)
1001  FORMAT(1X,132A1)
      END
      SUBROUTINE SPAP2 (LIN,LOUT,MAXSP,MAXEQ,JA,NIW,IA,IR,IC,
     @                  NH,ITYPE,XP,YP,INDEX,NRDIM)
      INTEGER LIN,LOUT
      REAL XP(*),YP(*),TEXT(8)
      INTEGER*2 JA(NIW),IA(NH),IR(NH),IC(NH),INDEX(*)
      INTEGER IFLAG(10)
      DATA  CHART/4H*TEM/ , CHARB/4H    /
      REWIND LIN
      DO 1 I=1,8
1     TEXT(I)=CHARB
      READ(LIN,9000) IHEAD
      IHEAD=IHEAD*36
      IF(IHEAD.EQ.0) GOTO 5
      READ(LIN,9011) (TEXT(L),L=1,8)
      BACKSPACE LIN
      READ(LIN,9005) (JA(L),L=1,IHEAD)
9000  FORMAT(19X,I5)
9011  FORMAT(8A4)
9005  FORMAT(36A2)
5     READ(LIN,9000) IEL
      IF(IEL.EQ.0) GOTO 10
      DO 9 I=1,IEL
9     READ(LIN,9001) CHAR
9001  FORMAT(A4)
10    READ(LIN,9000) NSP
      IF(NSP.EQ.0) GOTO 12
      IF(NSP.GE.MAXSP) GOTO 510
      DO 11 I=1,10000
      READ(LIN,9001) CHAR
11    IF(CHAR.EQ.CHART) GOTO 33
12    READ(LIN,9001) CHAR
33    READ(LIN,9001) CHAR
      READ(LIN,9000) NCEQ
      IF(NCEQ.EQ.0) GOTO 14
      IF(NCEQ.GT.MAXEQ) GOTO 520
      DO 13 I=1,NCEQ
13    READ(LIN,9001) CHAR
14    READ(LIN,9000) NLR
      IF(NCEQ.GT.NIW.OR.NLR.GT.NIW) GOTO 530
      READ(LIN,9500) (JA(J),J=1,NCEQ)
      NCEQ1=NCEQ+1
      NCEQ2=2*NCEQ
      READ(LIN,9500) (JA(J),J=NCEQ1,NCEQ2)
      READ(LIN,9500) (JA(J),J=1,NLR)
      READ(LIN,9000) LDIM
      READ(LIN,9500) (JA(J),J=1,LDIM)
      IF(LDIM.GT.NIW) GOTO 530
      NSP1=NSP+1
      READ(LIN,9500) (IA(J),J=1,NSP1)
      READ(LIN,9500) (IR(J),J=1,NSP)
9500  FORMAT(9I8)
      IF(ITYPE.GT.2.AND.NRDIM.LT.LDIM) GOTO 550
      N=NSP
      N2=N+2
      DO 80 I=1,N2
      IR(I)=I
      IC(I)=I
80    CONTINUE
      DO 200 I=1,N
      IMIN=IA(IR(I))
      IMAX=IA(IR(I)+1)-1
      DO 202 J=IMIN,IMAX
      INDEX(J)=0
      ICOL=IC(JA(J))
      IROW=NSP+1-I
      XP(J)=FLOAT(ICOL)
      YP(J)=FLOAT(IROW)
202   CONTINUE
200   CONTINUE
      ICASE=0
      H=.3
      XLEN=FLOAT(2*NSP+2)*H
      YLEN=XLEN
      RED=2*H
      DO 300 I=1,LDIM
      XP(I)=XP(I)*RED
      YP(I)=YP(I)*RED
300   CONTINUE
      IFLAG(1)=1
      IFLAG(2)=1
      IFLAG(3)=1
      IF(ITYPE.EQ.4)
     1        CALL DLG (XLEN,YLEN,H,XP,YP,LDIM,TEXT,IFLAG)
      CALL PLOMAGR (NSP,NSP,XLEN,YLEN,H,XP,YP,INDEX,LDIM,ICASE,LOUT,
     1            TEXT,IFLAG)
      REWIND LIN
      RETURN
510   WRITE(LOUT,8110) NSP,MAXSP
8110  FORMAT(/,52H Actual dimension NSP greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
520   WRITE(LOUT,8120) NCEQ,MAXEQ
8120  FORMAT(/,53H Actual dimension NCEQ greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
530   WRITE(LOUT,8130)
8130  FORMAT(/,41H INTEGER work space exhausted; enlarge IW)
      GOTO 599
550   WRITE(LOUT,8150)
8150  FORMAT(/,38H REAL work space exhausted; enlarge RW)
599   REWIND LIN
      ITYPE=-1
      RETURN
      END
      SUBROUTINE STOI1 (LIN,LOUT,MAXSP,MAXEQ,ROW,PLS,PRS,LR,NIW,ITYPE)
      INTEGER LIN,LOUT
      INTEGER*2 PLS(MAXEQ),PRS(MAXEQ),LR(NIW)
      INTEGER*2 ROW(*)
      INTEGER*2 PLUS,MINUS,BLANK,SIGN,RW,RS
      DATA  PLUS/1H+/ , MINUS/1H-/ , BLANK/1H / , RW/1H-/ , RS/1H^/
      DATA  CHART/4H*TEM/
      REWIND LIN
      MAXP=120
      MAXT=76
      READ(LIN,9000) IHEAD
      IHEAD=IHEAD*36
      IF(IHEAD.EQ.0) GOTO 5
      READ(LIN,9005) (LR(L),L=1,IHEAD)
9000  FORMAT(19X,I5)
9005  FORMAT(36A2)
      IF(ITYPE.LE.2) WRITE(LOUT,9003)
9003  FORMAT(//,19X,14HActual System:,/)
      IF(ITYPE.LE.2) WRITE(LOUT,9005) (LR(L),L=1,IHEAD)
5     READ(LIN,9000) IEL
      IF(IEL.EQ.0) GOTO 10
      DO 9 I=1,IEL
9     READ(LIN,9001) CHAR
9001  FORMAT(A4)
10    READ(LIN,9000) NSP
      IF(NSP.EQ.0) GOTO 12
      IF(NSP.GE.MAXSP) GOTO 510
      DO 11 I=1,10000
      READ(LIN,9001) CHAR
11    IF(CHAR.EQ.CHART) GOTO 33
12    READ(LIN,9001) CHAR
33    READ(LIN,9001) CHAR
      READ(LIN,9000) NCEQ
      IF(NCEQ.EQ.0) GOTO 14
      IF(NCEQ.GT.MAXEQ) GOTO 520
      DO 13 I=1,NCEQ
13    READ(LIN,9001) CHAR
14    READ(LIN,9000) NLR
      IF(NCEQ.GT.NIW.OR.NLR.GT.NIW) GOTO 530
      READ(LIN,9500) (PLS(J),J=1,NCEQ)
      READ(LIN,9500) (PRS(J),J=1,NCEQ)
      READ(LIN,9500) (LR(J),J=1,NLR)
9500  FORMAT(9I8)
      IF(ITYPE.EQ.1.AND.NSP.GT.MAXT) GOTO 550
      IF(ITYPE.EQ.2.AND.NSP.GT.MAXP) GOTO 550
      N=NSP
      N2=N+2
      IF(ITYPE.LE.2) WRITE(LOUT,1000)
      DO 100 L=1,N2
100   ROW(L)=RW
      IF(ITYPE.LE.2) WRITE(LOUT,1001) (ROW(L),L=1,N2)
      IR2=0
      DO 250 K=1,NCEQ
      DO 220 L=1,N
220   ROW(L+1)=BLANK
      ROW(N+2)=RS
      ROW(1)=RS
      IL1=IR2+1
      IL2=PLS(K)
      IR1=IL2+1
      IR2=PRS(K)
      IF(IL2.LT.IL1) GOTO 232
      DO 231 J=IL1,IL2
      SIGN=MINUS
      ROW(LR(J)+1)=SIGN
231   CONTINUE
232   CONTINUE
      IF(IR2.LT.IR1) GOTO 234
      DO 233 J=IR1,IR2
      SIGN=PLUS
      ROW(LR(J)+1)=SIGN
233   CONTINUE
234   CONTINUE
      IF(ITYPE.LE.2) WRITE(LOUT,1001) (ROW(L),L=1,N2)
250   CONTINUE
      DO 300 L=1,N2
300   ROW(L)=RW
      IF(ITYPE.LE.2) WRITE(LOUT,1001) (ROW(L),L=1,N2)
      REWIND LIN
      RETURN
510   WRITE(LOUT,8110) NSP,MAXSP
8110  FORMAT(/,52H Actual dimension NSP greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
520   WRITE(LOUT,8120) NCEQ,MAXEQ
8120  FORMAT(/,53H Actual dimension NCEQ greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
530   WRITE(LOUT,8130)
8130  FORMAT(/,41H INTEGER work space exhausted; enlarge IW)
      GOTO 599
550   WRITE(LOUT,8150)
8150  FORMAT(/,44H Problem size too large for this output type)
599   REWIND LIN
      ITYPE=-1
      RETURN
1000  FORMAT(1H ,//,2X,35HStructure of Stoichiometric Matrix:,/)
1001  FORMAT(1X,132A1)
      END
      SUBROUTINE STOI2 (LIN,LOUT,MAXSP,MAXEQ,PLS,PRS,LR,NIW,ITYPE,
     1                  XP,YP,INDEX,NRDIM)
      REAL TEXT(8)
      INTEGER LIN,LOUT
      REAL XP(*),YP(*)
      INTEGER*2 PLS(MAXEQ),PRS(MAXEQ),LR(NIW)
      INTEGER*2 INDEX(*)
      INTEGER IFLAG(10)
      DATA  CHART/4H*TEM/ , CHARB/4H    /
      REWIND LIN
      DO 1 I=1,8
1     TEXT(I)=CHARB
      READ(LIN,9000) IHEAD
      IHEAD=IHEAD*36
      IF(IHEAD.EQ.0) GOTO 5
      READ(LIN,9011) (TEXT(L),L=1,8)
      BACKSPACE LIN
      READ(LIN,9005) (LR(L),L=1,IHEAD)
9000  FORMAT(19X,I5)
9011  FORMAT(8A4)
9005  FORMAT(36A2)
5     READ(LIN,9000) IEL
      IF(IEL.EQ.0) GOTO 10
      DO 9 I=1,IEL
9     READ(LIN,9001) CHAR
9001  FORMAT(A4)
10    READ(LIN,9000) NSP
      IF(NSP.EQ.0) GOTO 12
      IF(NSP.GE.MAXSP) GOTO 510
      DO 11 I=1,10000
      READ(LIN,9001) CHAR
11    IF(CHAR.EQ.CHART) GOTO 33
12    READ(LIN,9001) CHAR
33    READ(LIN,9001) CHAR
      READ(LIN,9000) NCEQ
      IF(NCEQ.EQ.0) GOTO 14
      IF(NCEQ.GT.MAXEQ) GOTO 520
      DO 13 I=1,NCEQ
13    READ(LIN,9001) CHAR
14    READ(LIN,9000) NLR
      IF(NCEQ.GT.NIW.OR.NLR.GT.NIW) GOTO 530
      READ(LIN,9500) (PLS(J),J=1,NCEQ)
      READ(LIN,9500) (PRS(J),J=1,NCEQ)
      READ(LIN,9500) (LR(J),J=1,NLR)
9500  FORMAT(9I8)
      IF(ITYPE.GT.2.AND.NLR.GT.NRDIM) GOTO 550
      N=NSP
      IR2=0
      DO 250 K=1,NCEQ
      IL1=IR2+1
      IL2=PLS(K)
      IR1=IL2+1
      IR2=PRS(K)
      IF(IL2.LT.IL1) GOTO 232
      DO 231 J=IL1,IL2
      ISP=LR(J)
      INDEX(J)=2
      XP(J)=FLOAT(K)
      YP(J)=FLOAT(NSP+1-ISP)
231   CONTINUE
232   CONTINUE
      IF(IR2.LT.IR1) GOTO 234
      DO 233 J=IR1,IR2
      ISP=LR(J)
      INDEX(J)=1
      XP(J)=FLOAT(K)
      YP(J)=FLOAT(NSP+1-ISP)
233   CONTINUE
234   CONTINUE
250   CONTINUE
      ICASE=1
      H=.3
      XLEN=FLOAT(2*NCEQ+2)*H
      YLEN=FLOAT(2*NSP+2)*H
      RED=2*H
      DO 300 I=1,NLR
      XP(I)=XP(I)*RED
      YP(I)=YP(I)*RED
300   CONTINUE
      IFLAG(1)=1
      IFLAG(2)=1
      IFLAG(3)=1
      IF(ITYPE.EQ.4)
     1        CALL DLG (XLEN,YLEN,H,XP,YP,NLR,TEXT,IFLAG)
      CALL PLOMAGR (NSP,NCEQ,XLEN,YLEN,H,XP,YP,INDEX,NLR,ICASE,LOUT,
     1            TEXT,IFLAG)
      REWIND LIN
      RETURN
510   WRITE(LOUT,8110) NSP,MAXSP
8110  FORMAT(/,52H Actual dimension NSP greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
520   WRITE(LOUT,8120) NCEQ,MAXEQ
8120  FORMAT(/,53H Actual dimension NCEQ greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
530   WRITE(LOUT,8130)
8130  FORMAT(/,41H INTEGER work space exhausted; enlarge IW)
      GOTO 599
550   WRITE(LOUT,8150)
8150  FORMAT(/,38H REAL work space exhausted; enlarge RW)
599   REWIND LIN
      ITYPE=-1
      RETURN
      END
      SUBROUTINE EXPLI (LIN,LOUT,LOUT1,LOUT2,MAXSP,MAXEQ,LR,PLR,IA,
     @                  IDA,JA,NLRH,NCEQH2,NSPH1,NSPH,LDH,KFLAG)
      IMPLICIT INTEGER (A-Z)
      INTEGER*2 LR(NLRH),PLR(NCEQH2),IA(NSPH1),IDA(NSPH),JA(LDH)
      INTEGER*2 KEY(80),SIG(80),LQ(80),VZ(542),PDY(542)
      INTEGER*2 PLUS,PARL,PARR,STAR,MINUS,CC,YY
      LOGICAL   REP
      LOGICAL*1 IJ(542)
      DATA PLUS/1H+/ , PARL/1H(/ , PARR/1H)/ ,
     @     STAR/1H*/ , MINUS/1H-/ , CC/1HC/ , YY/1HY/
      DATA POWER/1HD/
      DATA  CHART/4H*TEM/
      REWIND LIN
      REWIND LOUT1
      REWIND LOUT2
      READ(LIN,9000) IHEAD
      IH=IHEAD*36
      IF(IH.EQ.0) GOTO 5
      READ(LIN,9005) (LR(L),L=1,IH)
9000  FORMAT(19X,I5)
9005  FORMAT(36A2)
      IF(KFLAG.EQ.5.OR.KFLAG.EQ.7) WRITE(LOUT1,9200)
      IF(KFLAG.EQ.6) WRITE(LOUT1,9400)
      IF(KFLAG.EQ.5) WRITE(LOUT2,9300)
      IF(KFLAG.EQ.7) WRITE(LOUT2,9700)
      WRITE(LOUT1,9100)
      IF(KFLAG.EQ.5.OR.KFLAG.EQ.7) WRITE(LOUT2,9100)
      DO 4 I=1,IHEAD
      L1=(I-1)*36+1
      L2=L1+35
      WRITE(LOUT1,9102) (LR(L),L=L1,L2)
      IF(KFLAG.EQ.5.OR.KFLAG.EQ.7) WRITE(LOUT2,9102) (LR(L),L=L1,L2)
4     CONTINUE
      WRITE(LOUT1,9100)
      IF(KFLAG.EQ.5.OR.KFLAG.EQ.7) WRITE(LOUT2,9100)
5     READ(LIN,9000) IEL
      IF(IEL.EQ.0) GOTO 10
      DO 9 I=1,IEL
9     READ(LIN,9001) CHAR
9001  FORMAT(A4)
10    READ(LIN,9000) NSP
      IF(NSP.EQ.0) GOTO 12
      IF(NSP.GT.MAXSP) GOTO 510
      DO 11 I=1,10000
      READ(LIN,9001) CHAR
11    IF(CHAR.EQ.CHART) GOTO 33
12    READ(LIN,9001) CHAR
33    READ(LIN,9001) CHAR
      READ(LIN,9000) NCEQ
      IF(NCEQ.EQ.0) GOTO 14
      IF(NCEQ.GT.MAXEQ) GOTO 520
      DO 13 I=1,NCEQ
13    READ(LIN,9001) CHAR
14    READ(LIN,9000) NLR
      IF(NLR.GT.NLRH) GOTO 540
      READ(LIN,9500) (PLR(J),J=1,NCEQ)
      NCEQ1=NCEQ+1
      NCEQ2=2*NCEQ
      READ(LIN,9500) (PLR(J),J=NCEQ1,NCEQ2)
      READ(LIN,9500) (LR(J),J=1,NLR)
      READ(LIN,9000) LDIM
      IF(LDIM.GT.LDH) GOTO 530
      READ(LIN,9500) (JA(J),J=1,LDIM)
      NSP1=NSP+1
      READ(LIN,9500) (IA(J),J=1,NSP1)
      READ(LIN,9500) (IDA(J),J=1,NSP)
9500  FORMAT(9I8)
      IF(KFLAG.NE.6) GOTO 16
      WRITE(LOUT1,9402)
      WRITE(LOUT1,9404) LDIM
      WRITE(LOUT1,9405) NSP1
      WRITE(LOUT1,9406) NSP
      WRITE(LOUT1,9411) (J,IDA(J),J=1,NSP)
      WRITE(LOUT1,9414) (J,IA(J),J=1,NSP1)
      WRITE(LOUT1,9412) (J,JA(J),J=1,LDIM)
      WRITE(LOUT1,9100)
      WRITE(LOUT1,9450)
      RETURN
16    WRITE(LOUT1,9203) NCEQ
      WRITE(LOUT1,9100)
      WRITE(LOUT1,9207) NCEQ
      WRITE(LOUT1,9208) NCEQ
      WRITE(LOUT1,9100)
      IF(KFLAG.EQ.5) WRITE(LOUT2,9303)
      IF(KFLAG.EQ.7) WRITE(LOUT2,9703)
      WRITE(LOUT2,9100)
      IF(KFLAG.EQ.5) WRITE(LOUT2,9307) NCEQ
      IF(KFLAG.EQ.7) WRITE(LOUT2,9707) NCEQ
      IF(KFLAG.EQ.7) WRITE(LOUT2,9708) POWER
      WRITE(LOUT2,9100)
      IR2=0
      IA2=0
      DO 1100 I=1,NCEQ
      IL1=IR2+1
      IL2=PLR(I)
      IR2=PLR(NCEQ+I)
      IF(IL2.LT.IL1) GOTO 1101
      WRITE(LOUT1,9210) I,I,(STAR,LR(L),L=IL1,IL2)
      GOTO 1100
1101  WRITE(LOUT1,9211) I,I
9211  FORMAT(T7,2HR(,I3,4H)=C(,I3,1H))
1100  CONTINUE
      DO 1200 I=1,NSP
      PS=0
      IR2=0
      R2=NCEQ
      DO 1210 J=1,NCEQ
      IORD=0
      IL1=IR2+1
      IL2=PLR(J)
      IR1=IL2+1
      IR2=PLR(NCEQ+J)
      IF(IL2.LT.IL1) GOTO 1221
      DO 1220 L=IL1,IL2
      IF(LR(L).NE.I) GOTO 1220
      IORD=IORD+1
      PS=PS+1
      PDY(PS)=J
      VZ(PS)=MINUS
1220  CONTINUE
1221  CONTINUE
      IF(IR2.LT.IR1) GOTO 1232
      DO 1230 R=IR1,IR2
      IF(LR(R).NE.I) GOTO 1230
      IORD=IORD-1
      IF(IORD.LT.0) GOTO 1231
      PS=PS-1
      GOTO 1230
1231  PS=PS+1
      PDY(PS)=J
      VZ(PS)=PLUS
1230  CONTINUE
1232  CONTINUE
      IF(IORD.EQ.0) GOTO 1210
      DO 1241 L=IL1,IL2
1241  IJ(LR(L))=.FALSE.
1210  CONTINUE
      IF(PS.GT.0) GOTO 1250
      I1=I+1
      WRITE(LOUT1,9218) I,POWER
      IF(KFLAG.EQ.5) WRITE(LOUT2,9310) IDA(I),POWER
      IF(KFLAG.EQ.7) WRITE(LOUT2,9710) I,I,POWER
      GOTO 1200
1250  L2=PS
      IF(L2.GT.53) L2=53
      WRITE(LOUT1,9220) I,(VZ(L),PDY(L),L=1,L2)
1251  L1=L2+1
      L2=L2+52
      IF(L2.GT.PS) L2=PS
      IF(L2.LE.L1) GOTO 1260
      WRITE(LOUT1,9219) I,I,(VZ(L),PDY(L),L=L1,L2)
      GOTO 1251
1260  CONTINUE
      IA1=IA2+1
      IA2=IA(I+1)-1
      DO 1300 JJ=IA1,IA2
      J=JA(JJ)
      P1=0
      IR2J=0
      REP=.FALSE.
      DO 1350 K=1,PS
      IF(PDY(K).EQ.1) IL1J=1
      IF(PDY(K).GE.2) IL1J=PLR(PDY(K)-1+NCEQ)+1
      IL2J=PLR(PDY(K))
      IORD=0
      DO 1351 L=IL1J,IL2J
      IF(LR(L).EQ.J) GOTO 1352
1351  CONTINUE
      GOTO 1350
1352  P1=P1+1
      R1=P1
      LQ(P1)=PDY(K)
      KEY(P1)=CC
      SIG(P1)=VZ(K)
      DO 1355 L=IL1J,IL2J
      LSL=LR(L)
      P1=P1+1
      LQ(P1)=LSL
      KEY(P1)=YY
      SIG(P1)=STAR
      IF(LSL.NE.J) GOTO 1355
      IF(IORD.EQ.0) P1=P1-1
      IORD=IORD+1
1355  CONTINUE
1356  IF(IORD.LE.1) GOTO 1358
      IORD=IORD-1
      R2=P1
      DO 1357 R=R1,R2
      P1=P1+1
      LQ(P1)=LQ(R)
      KEY(P1)=KEY(R)
1357  SIG(P1)=SIG(R)
      GOTO 1356
1358  IF(P1.LE.52) GOTO 1350
      IF(KFLAG.EQ.7) GOTO 13000
      IF(.NOT.REP) WRITE(LOUT2,9330) JJ,(SIG(L),KEY(L),PARL,
     @                                    LQ(L),PARR,L=1,P1A)
      IF(REP) WRITE(LOUT2,9331) JJ,JJ,(SIG(L),KEY(L),PARL,
     @                                   LQ(L),PARR,L=1,P1A)
      GOTO 13001
13000 CONTINUE
      ICOL=JA(JJ)
      IF(.NOT.REP) WRITE(LOUT2,9730) I,ICOL,(SIG(L),KEY(L),PARL,
     @                                    LQ(L),PARR,L=1,P1A)
      IF(REP) WRITE(LOUT2,9731) I,ICOL,I,ICOL,(SIG(L),KEY(L),PARL,
     @                                   LQ(L),PARR,L=1,P1A)
13001 REP=.TRUE.
      R=0
      P1A=P1A+1
      DO 1359 L=P1A,P1
      R=R+1
      LQ(R)=LQ(L)
      KEY(R)=KEY(L)
1359  SIG(R)=SIG(L)
      P1=R
1350  P1A=P1
      IF(P1A.EQ.0) GOTO 1390
      IF(KFLAG.EQ.7) GOTO 13005
      IF(.NOT.REP) WRITE(LOUT2,9330) JJ,(SIG(L),KEY(L),PARL,
     @                                    LQ(L),PARR,L=1,P1A)
      IF(REP) WRITE(LOUT2,9331) JJ,JJ,(SIG(L),KEY(L),PARL,
     @                                   LQ(L),PARR,L=1,P1A)
      GOTO 1300
13005 CONTINUE
      ICOL=JA(JJ)
      IF(.NOT.REP) WRITE(LOUT2,9730) I,ICOL,(SIG(L),KEY(L),PARL,
     @                                    LQ(L),PARR,L=1,P1A)
      IF(REP) WRITE(LOUT2,9731) I,ICOL,I,ICOL,(SIG(L),KEY(L),PARL,
     @                                   LQ(L),PARR,L=1,P1A)
      GOTO 1300
1390  IF(KFLAG.EQ.5) WRITE(LOUT2,9310) JJ,POWER
      ICOL=JA(JJ)
      IF(KFLAG.EQ.7) WRITE(LOUT2,9710) I,ICOL,POWER
1300  CONTINUE
1200  CONTINUE
      WRITE(LOUT1,9100)
      WRITE(LOUT1,9240)
      WRITE(LOUT1,9100)
      WRITE(LOUT1,9100)
      WRITE(LOUT1,9250)
      WRITE(LOUT2,9100)
      WRITE(LOUT2,9350)
      REWIND LIN
      ENDFILE LOUT1
      ENDFILE LOUT2
      REWIND LOUT1
      REWIND LOUT2
      RETURN
510   WRITE(LOUT,8110) NSP,MAXSP
8110  FORMAT(/,52H Actual dimension NSP greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
520   WRITE(LOUT,8120) NCEQ,MAXEQ
8120  FORMAT(/,53H Actual dimension NCEQ greater than defined max value,
     @       /,I4,6H .GT. ,I4)
      GOTO 599
530   WRITE(LOUT,8130)
8130  FORMAT(/,41H INTEGER work space exhausted; enlarge IW)
      GOTO 599
540   WRITE(LOUT,8140)
8140  FORMAT(/,28H Enlarge definition of MAXEQ)
599   REWIND LIN
      KFLAG=-1
      RETURN
9100  FORMAT(1HC)
9102  FORMAT(2HC ,35A2)
9200  FORMAT(6X,32HSUBROUTINE FCN (NSP,NCEQ,T,Y,DY))
9203  FORMAT(6X,38HDOUBLE PRECISION Y(NSP),DY(NSP),T,C,R(,I3,1H))
9207  FORMAT(6X,20HCOMMON / KINPAR / C(,I3,1H))
9208  FORMAT(6X,20HCOMMON / RATE1 / RH(,I3,1H))
9210  FORMAT(6X,2HR(,I3,4H)=C(,I3,1H),6(A1,2HY(,I3,1H)))
9220  FORMAT(6X,3HDY(,I3,4H) = ,8(A1,2HR(,I3,1H)),
     @       5(/,5X,4H@   ,9(A1,2HR(,I3,1H))))
9219  FORMAT(6X,3HDY(,I3,7H) = DY(,I3,1H),7(A1,2HR(,I3,1H)),
     @       5(/,5X,4H@   ,9(A1,2HR(,I3,1H))))
9218  FORMAT(6X,3HDY(,I3,6H) = 0.,A1,1H0)
9240  FORMAT(6X,13HDO 1 I=1,NCEQ,/,1H1,5X,10HRH(I)=R(I))
9250  FORMAT(23HC    END SUBROUTINE FCN,/,1HC,/,6X,6HRETURN,
     @       /,6X,3HEND)
9300  FORMAT(6X,39HSUBROUTINE JACOBI (NSP,NCEQ,LDIM,T,Y,A))
9303  FORMAT(6X,35HDOUBLE PRECISION A(LDIM),Y(NSP),T,C)
9307  FORMAT(6X,20HCOMMON / KINPAR / C(,I3,1H))
9350  FORMAT(26HC    END SUBROUTINE JACOBI,/,1HC,/,6X,6HRETURN,
     @       /,6X,3HEND)
9310  FORMAT(6X,2HA(,I5,6H) = 0.,A1,1H0)
9330  FORMAT(6X,2HA(,I5,2H)=,8(3A1,I3,A1),
     @              6(/,5X,3H@  ,9(3A1,I3,A1)))
9331  FORMAT(6X,2HA(,I5,4H)=A(,I5,1H),7(3A1,I3,A1),
     @                      6(/,5X,4H@   ,9(3A1,I3,A1)))
9400  FORMAT(6X,17HSUBROUTINE STRUCT)
9402  FORMAT(6X,18HINTEGER*2 JA,IA,DP)
9404  FORMAT(6X,20HCOMMON / LARK3 / JA(,I4,1H))
9405  FORMAT(6X,20HCOMMON / LARK4 / IA(,I4,1H))
9406  FORMAT(6X,20HCOMMON / LARK5 / DP(,I3,1H))
9411  FORMAT(6X,3HDP(,I3,4H) = ,I5)
9412  FORMAT(6X,3HJA(,I5,4H) = ,I3)
9414  FORMAT(6X,3HIA(,I4,4H) = ,I5)
9450  FORMAT(26HC    END SUBROUTINE STRUCT,/,1HC,/,6X,6HRETURN,
     @       /,6X,3HEND)
9700  FORMAT(6X,39HSUBROUTINE JACOBI (NDCL,NSP,NCEQ,T,Y,A))
9703  FORMAT(6X,39HDOUBLE PRECISION A(NDCL,NSP),Y(NSP),T,C)
9707  FORMAT(6X,20HCOMMON / KINPAR / C(,I3,1H))
9708  FORMAT(1HC,/,6X,12HDO 1 I=1,NSP,/,6X,12HDO 1 J=1,NSP,
     @       /,1H1,5X,9HA(I,J)=0.,A1,1H0)
9710  FORMAT(6X,2HA(,I3,1H,,I3,6H) = 0.,A1,1H0)
9730  FORMAT(6X,2HA(,I3,1H,,I3,2H)=,8(3A1,I3,A1),
     @              6(/,5X,3H@  ,9(3A1,I3,A1)))
9731  FORMAT(6X,2HA(,I3,1H,,I3,4H)=A(,I3,1H,,I3,1H),7(3A1,I3,A1),
     @                      6(/,5X,4H@   ,9(3A1,I3,A1)))
      END
      SUBROUTINE DLG (XLEN,YLEN,H,XP,YP,NDIM,TEXT,IFLAG)
      REAL XP(*),YP(*),TEXT(8)
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT,INFDAT
      INTEGER IFLAG(10)
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
      COMMON / INFORM / IHELP,INFDAT
C
5     LIN=ITIN
      LOUT=ITOUT
C
10    CONTINUE
      WRITE(LOUT,9000) XLEN,YLEN,(TEXT(L),L=1,8)
9000  FORMAT(/,22H  Enter Option Number:,/
     1        ,13H  0:     Quit,/
     2        ,17H  1:     XLENGTH=,F7.2,11H   YLENGTH=,F7.2,/
     3        ,15H  2:     Text:",8A4,1H")
      IF(IFLAG(2).EQ.1) WRITE(LOUT,9021)
      IF(IFLAG(2).EQ.0) WRITE(LOUT,9020)
9020  FORMAT(26H  3:     Title not plotted)
9021  FORMAT(22H  3:     Title plotted)
      IF(IFLAG(3).EQ.1) WRITE(LOUT,9031)
      IF(IFLAG(3).EQ.0) WRITE(LOUT,9030)
9030  FORMAT(27H  4:     Legend not plotted)
9031  FORMAT(23H  4:     Legend plotted)
      READ(LIN,*) IOPT
      IF(IOPT.EQ.0) RETURN
      IF(IOPT.LT.0.OR.IOPT.GT.4) GOTO 10
      IF(IOPT.EQ.4.AND.IFLAG(3).EQ.0) IFH=1
      IF(IOPT.EQ.4.AND.IFLAG(3).EQ.1) IFH=0
      IF(IOPT.EQ.4) IFLAG(3)=IFH
      IF(IOPT.EQ.3.AND.IFLAG(2).EQ.0) IFH=1
      IF(IOPT.EQ.3.AND.IFLAG(2).EQ.1) IFH=0
      IF(IOPT.EQ.3) IFLAG(2)=IFH
      IF(IOPT.EQ.2) WRITE(LOUT,9200)
9200  FORMAT(35H  Enter Text (up to 32 characters):)
      IF(IOPT.EQ.2) READ(LIN,9210) (TEXT(L),L=1,8)
9210  FORMAT(8A4)
      IF(IOPT.NE.1) GOTO 10
      WRITE(LOUT,9100)
9100  FORMAT(25H  Enter rescaling factor:)
      READ(LIN,*) RED
      H=RED*H
      XLEN=RED*XLEN
      YLEN=RED*YLEN
      DO 50 I=1,NDIM
      XP(I)=RED*XP(I)
      YP(I)=RED*YP(I)
50    CONTINUE
      GOTO 10
      END
 
