      SUBROUTINE KIP (NRW,NSP,NCEQ,RW,TEMP,RK,GAS)
      DOUBLE PRECISION RW(NRW),RK(NCEQ),GAS,TEMP,XL10,RKTEST
      DATA RKTEST/-1.D+30/
      IH1=2*NSP
      IH2=IH1+NCEQ
      IH3=IH2+NCEQ
      RT=GAS*TEMP
      XL10=DLOG(1.D1)
      DO 4 I=1,NCEQ
      I1=IH1+I
      I2=IH2+I
      I3=IH3+I
      IF(RW(I3).GT.RKTEST) GOTO 3
      IF(RW(I2).GT.RKTEST) GOTO 2
      IF(RW(I1).LE.RKTEST) RW(I1)=0.D0
      RK(I)=RW(I1)
      GOTO 4
2     CONTINUE
      RK(I)=DEXP(RW(I1)*XL10-RW(I2)/RT)
      GOTO 4
3     RK(I)=DEXP(RW(I1)*XL10-RW(I2)/RT)*TEMP**RW(I3)
4     CONTINUE
      RETURN
      END
      SUBROUTINE SIMULA (RW,NRW,IW,NIW,KFLAG)
      REAL R,RINT,RTMP
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT,INFDAT
      INTEGER*2 IA,JA,IDA,LR,PLR,IW
      INTEGER*2 ELEM(10,5)
      DOUBLE PRECISION TP(50),ATWS(10),RW
      DOUBLE PRECISION EPS,GAS,HMAX,RK,TEMP
      DIMENSION RW(NRW),IW(NIW)
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
      COMMON / INFORM / IHELP,INFDAT
      COMMON / LARK1 / LR(4000)
      COMMON / LARK2 / PLR(1800)
      COMMON / LARK3 / JA(2800)
      COMMON / LARK4 / IA(401)
      COMMON / LARK5 / IDA(400)
      COMMON / KINPAR / RK(900)
      COMMON / RATE1 / R(900)
      COMMON / RATE2 / RINT(900)
      COMMON / RATE3 / RTMP(900)
      REWIND SIN
      REWIND DOUT
      IF(90.GT.NIW) GOTO 1000
      READ(SIN,9000) IHEAD
      IHEAD=IHEAD*18
      IF(IHEAD.LE.0) GOTO 5
      READ(SIN,9005) (R(L),L=1,IHEAD)
9000  FORMAT(19X,I5)
9005  FORMAT(18A4)
5     IF(KFLAG.GE.1) WRITE(ITOUT,9003)
9003  FORMAT(//,19X,13HACTUAL MODEL:,/,19X,13(1H-),/)
      IF(KFLAG.GE.1.AND.IHEAD.GT.0) WRITE(ITOUT,9006) (R(L),L=1,IHEAD)
      READ(SIN,9000) IEL
      IF(IEL.LE.0) GOTO 11
      DO 10 I=1,IEL
10    READ(SIN,9007) (ELEM(I,L),L=1,5),ATWS(I)
9007  FORMAT(5A1,3X,D21.12)
11    READ(SIN,9000) NSP
      IH=13*NSP
      IF(IH.GE.NIW) GOTO 1000
      IH=IH-NSP
      IF(IH.GE.NRW) GOTO 1001
      IH1=5*NSP
      DO 20 I=1,NSP
      IF(IEL.GT.0) READ(SIN,9009) IH,(IW(5*(IH-1)+L),L=1,5),RW(IH),
     @             RW(NSP+IH),(IW(IH1+IEL*(IH-1)+L),L=1,IEL)
20    IF(IEL.LE.0) READ(SIN,9010) IH,(IW(5*(IH-1)+L),L=1,5),RW(IH)
9009  FORMAT(I4,1X,5A2,1X,D21.12,1X,D21.12,1X,2(I8),/,8(I8))
9010  FORMAT(I5,1X,5A2,2X,D21.12)
      READ(SIN,9015) TEMP
9015  FORMAT(19X,D21.12)
      READ(SIN,9016) GAS
9016  FORMAT(20X,D21.12)
      READ(SIN,9000) NCEQ
      IH1=2*NSP
      IH2=IH1+NCEQ
      IH3=IH2+NCEQ
      DO 30 I=1,NCEQ
      I1=IH1+I
      I2=IH2+I
      I3=IH3+I
30    READ(SIN,9040) RW(I1),RW(I2),RW(I3)
      CALL KIP (NRW,NSP,NCEQ,RW,TEMP,RK,GAS)
      READ(SIN,9000) NLR
      READ(SIN,9500) (PLR(J),J=1,NCEQ)
      NCEQ1=NCEQ+1
      NCEQ2=2*NCEQ
      READ(SIN,9500) (PLR(J),J=NCEQ1,NCEQ2)
      READ(SIN,9500) (LR(J),J=1,NLR)
      READ(SIN,9000) LDIM
      IH=13*NSP+2*LDIM
      IF(IH.GE.NIW) GOTO 1000
      IH=IH-NSP
      IF(IH.GE.NRW) GOTO 1001
      READ(SIN,9500) (JA(J),J=1,LDIM)
      NSP1=NSP+1
      READ(SIN,9500) (IA(J),J=1,NSP1)
      READ(SIN,9500) (IDA(J),J=1,NSP)
      READ(SIN,9015) EPS
      READ(SIN,9000) ITPM
      READ(SIN,9050) (TP(I),I=1,ITPM)
9050  FORMAT(D21.12)
      READ(SIN,9060) IPRINT
      IH=IHEAD/18
      WRITE(DOUT,9001) IH
9001  FORMAT(6H &HEAD,9X,I3)
      IF(IHEAD.GT.0) WRITE(DOUT,9006) (R(L),L=1,IHEAD)
      WRITE(DOUT,9999)
9006  FORMAT(1X,18A4)
      WRITE(DOUT,9011)
9011  FORMAT(5H &DIM)
      WRITE(DOUT,9012) NSP,NCEQ,IEL
9012  FORMAT(5H NSP=,I6,7H, NCEQ=,I6,6H, IEL=,I6,1H,)
      WRITE(DOUT,9999)
      WRITE(DOUT,9052) ITPM
9052  FORMAT(6H &TOUT,/,8H ITPMAX=,I3,5H, TP=)
      WRITE(DOUT,9053) (TP(L),L=1,ITPM)
9053  FORMAT(3(1X,D21.12,1H,))
      WRITE(DOUT,9999)
      IF(IEL.LE.0) GOTO 40
      WRITE(DOUT,9021)
9021  FORMAT(7H &ELEMN)
      WRITE(DOUT,9022) ((ELEM(I,L),L=1,5),I=1,IEL)
9022  FORMAT(7(1X,5A1,4H   ,))
      WRITE(DOUT,9999)
40    CONTINUE
      IF(IEL.LE.0) GOTO 50
      WRITE(DOUT,9025)
9025  FORMAT(6H &EATW,/,4H RK=)
      WRITE(DOUT,9063) (ATWS(L),L=1,IEL)
      WRITE(DOUT,9999)
50    CONTINUE
      WRITE(DOUT,9031)
9031  FORMAT(6H &NAME)
      IH=NSP*5
      WRITE(DOUT,9032) (IW(L),L=1,IH)
9032  FORMAT(5(1X,5A2,3H  ,))
      WRITE(DOUT,9999)
      IF(IEL.LE.0) GOTO 60
      WRITE(DOUT,9035)
9035  FORMAT(6H &SATW,/,4H RK=)
      WRITE(DOUT,9063) (RW(NSP+I),I=1,NSP)
      WRITE(DOUT,9999)
60    CONTINUE
      WRITE(DOUT,9062)
9062  FORMAT(6H &KIPA,/,4H RK=)
      WRITE(DOUT,9063) (RK(L),L=1,NCEQ)
9063  FORMAT(5(1X,D13.7,1H,))
      WRITE(DOUT,9999)
      N=NSP
      NRWH=NRW-N
      HMAX=TP(ITPM)-TP(1)
      CALL META (N,NCEQ,LDIM,NRWH,NIW,RW(1),TP,ITPM,HMAX,EPS,KFLAG,
     @           IPRINT,RW(N+1),IW)
      ENDFILE DOUT
      RETURN
1000  WRITE(ITOUT,9800)
9800  FORMAT(/,41H INTEGER WORK SPACE EXHAUSTED; ENLARGE IW)
      GOTO 1111
1001  WRITE(ITOUT,9801)
9801  FORMAT(/,38H REAL WORK SPACE EXHAUSTED; ENLARGE RW)
1111  KFLAG=-1
      RETURN
9040  FORMAT(3D24.14)
9060  FORMAT(19X,I5)
9500  FORMAT(9I8)
9999  FORMAT(5H &END)
      END
      SUBROUTINE META (N,NCEQ,LDIM,NW,NIW,Y,TP,ITPM,HMAX,EPS,KFLAG,
     @                 IPRINT,W,IW)
      INTEGER*2 IW(NIW)
      DOUBLE PRECISION W(NW)
      DOUBLE PRECISION Y(N),TP(ITPM)
      DOUBLE PRECISION CON1,CON2,EPS,HI,HMAX
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT,INFDAT
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
      COMMON / INFORM / IHELP,INFDAT
      COMMON / COUNT / NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL
      KF=KFLAG
      NSTEP=0
      NFCN=0
      NANFA=0
      NTFAC=0
      NFAC=0
      NSOL=0
      MW=NW-11*N-LDIM
      MWH=MW
      LD=LDIM
      MIW=NIW-13*N
      MIWH=MIW
      M2=N+N
      M4=M2+M2
      M6=M4+M2
      NP1=N+1
      NP2=NP1+N
      NP3=NP2+N
      NP4=NP3+N
      NP5=NP4+N
      NP6=NP5+N
      NP6=NP5+N
      NP7=NP6+N
      NP7=NP6+N
      NP11=11*N+1
      NP13=13*N+1
      NH=11*N+LDIM+1
      HI=0.D0
      CON1=1.D-9
      CON2=1.D0
      IF(EPS.LT.0.D0) CON2=CON1
      IF(EPS.LT.0.D0) EPS=-EPS
      CALL METASC (N,Y,TP,ITPM,EPS,CON1,CON2,HMAX,HI,KFLAG,IPRINT,
     @     W(1),W(NP1),W(NP2),W(NP3),W(NP4),W(NP5),W(NP6),W(NP11),
     @     IW(1),IW(NP1),IW(NP5),IW(NP7),IW(NP13),
     @     M2,M4,M6,MWH,LDIM,MIWH,NCEQ,W(NH))
      IF(KF.LT.2) GOTO 50
      WRITE(ITOUT,8050) NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL
50    RETURN
8050  FORMAT(/,13H ************,/,13H  STATISTICS ,/,13H ************,/,
     @       10H STEPS    ,I10,/,
     @       10H FCN-EVAL ,I10,/,
     @       10H ANAL/FAC ,I10,/,
     @       10H TEST-FAC ,I10,/,
     @       10H FACTOR   ,I10,/,
     @       10H SUBST    ,I10)
      END
      SUBROUTINE METASC (N,Y,TP,ITPM,EPS,CON1,CON2,HMAX,HI,KFLAG,IPRINT,
     @        YM,DY,DZ,DEL,SH,SM,DT,A,IS,IKEEP,IWO,IW1,ILEQ,
     @           M2,M4,M6,NWTMP,LDIM,NLEQ,NCEQ,WTMP)
      DOUBLE PRECISION Y(N),YM(N),DY(N),DZ(N),DEL(N)
      DOUBLE PRECISION SH(N),SM(N),DT(N,5)
      INTEGER*2 IA,IDA,JA,IS(N),IKEEP(M4),IW1(M6),ILEQ(NLEQ),IWO(M2)
      INTEGER IDISP(2)
      DOUBLE PRECISION A(LDIM),WTMP(NWTMP)
      DOUBLE PRECISION TP(ITPM),  D(5,5), AJ(5), AL(5,5)
      REAL R,RINT,RTMP
      INTEGER NJ(5),INCR(5) , NRED(4)
      INTEGER EOUT,MOUT
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT,INFDAT
      DOUBLE PRECISION B1,C,CONTRA,CON1,CON2,COSTF,COSTJ
      DOUBLE PRECISION COSTLR,COSTS,DABS,DEPS,DFLOAT,DM,DMA,DMH
      DOUBLE PRECISION DSQRT,EPH,EPMA,EPMACH,EMIN,EPS,ERR,FC,FCM,FCO
      DOUBLE PRECISION FJ,FJ1,FMIN,FN,G,GD,GH,H,HALF,HH,HMAX,HMAXU
      DOUBLE PRECISION HMIN,HI,HR,H1,OMJ,OMJO,ONE,ONE1,ONE2,QUART
      DOUBLE PRECISION RED,RMIN,ROH,SAFE,SAFEDM,SCALET,SHH,SHI,SMALL
      DOUBLE PRECISION SMALLH,T,TA,TEN,TENTH,TEND,TGROW,TH
      DOUBLE PRECISION THRESH,TN,TOLD,TWO,U,V,W,XKAMIN
      DOUBLE PRECISION XKAACT,YMAX,YMAXH,ZERO,DIFF
      LOGICAL PRERR,PRM,PRM1,PRM2,WINDOW,LPRINT,LRATE
      LOGICAL BD1,BD2,BD3,BANFA
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
      COMMON / INFORM / IHELP,INFDAT
      COMMON / MACHIN / EPMA,SMALLH
      COMMON / COUNT / NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL
      COMMON / LARK3 / JA(2800)
      COMMON / LARK4 / IA(401)
      COMMON / LARK5 / IDA(400)
      COMMON / MA30LE / LP,BD1,BD2,BD3
      COMMON / MA30LF / IRNCP,ICNCP,IRANK,MINIRN,MINICN
      COMMON / MA30LG / DEPS,RMIN
      COMMON / RATE1 / R(900)
      COMMON / RATE2 / RINT(900)
      COMMON / RATE3 / RTMP(900)
      DATA  ZERO/0.D0/ , FMIN/1.D-2/ , QUART/0.25D0/ , HALF/0.5D0/
      DATA  ROH/0.25D0/ , SAFE/0.5D0/ , ONE/1.D0/ , ONE1/1.01D0/
      DATA ONE2 /1.1D0/ , TWO/2.D0/ , TEN/1.D1/ , TENTH/1.D-1/
      DATA SAFEDM/0.8D0/
       EOUT=ITOUT
       MOUT=ITOUT
       LP=EOUT
       EPMACH=EPMA*TEN
       SMALL=SMALLH*TEN**4
      ISMAX=1000
      JRMAX=7
      THRESH=1.D-2
      IF(EPMACH.GT.1.D-10) THRESH=THRESH*5.D0
      IF(EPMACH.GT.1.D-5) THRESH=0.5D0
      TGROW=1.D-6
      IF(EPMACH.GT.1.D-12) TGROW=1.D-5
      IF(EPMACH.GT.1.D-8) TGROW=1.D-3
      IF(EPMACH.GT.1.D-5) TGROW=1.D-1
      NZAMAX=0
      PRERR=.FALSE.
      PRM=.FALSE.
      PRM1=.FALSE.
      PRM2=.FALSE.
      IF(KFLAG.GE.0) PRERR=.TRUE.
      IF(KFLAG.EQ.1) PRM1=.TRUE.
      IF(KFLAG.GE.2) PRM2=.TRUE.
      IF(PRM1.OR.PRM2) PRM=.TRUE.
      IF(.NOT.PRERR) LP=0
      KFLAG=0
      IF(PRM) WRITE(MOUT,9000)
      ISTEP=0
      NP1=N+1
      N2=N+N
      NP2=N2+1
      N3=N2+N
      NP3=N3+1
      N4=N3+N
      NP4=N4+1
      N5=N4+N
      NP5=N5+1
      IELBOW=N
      XKAMIN=1.3D0
      LIRN=LDIM+IELBOW
      LICN=NLEQ-LIRN
      IF(LICN.GT.32000) LICN=32000
      LICNR=NWTMP-N
      NFREE=LICN-LICNR
      IF(NFREE.GT.N2.AND.PRM2) WRITE(MOUT,9300) NFREE
      NFREE=LICNR-LICN
      IF(NFREE.GT.N2.AND.PRM2) WRITE(MOUT,9305) NFREE
      IF(LICNR.LT.LICN) LICN=LICNR
      XKAACT=DFLOAT(LICN)/DFLOAT(LDIM)
      IF(XKAACT.LT.XKAMIN) GOTO 32
      LICNP1=LICN+1
      IDISP(1)=1
      BANFA=.TRUE.
      DEPS=TWO
4322   NJ(1)=2
       NJ(2)=6
       NJ(3)=10
       NJ(4)=14
       NJ(5)=22
      KM=4
      JM=KM+1
      FN=DFLOAT(N)
      FJ1=DFLOAT(NJ(1))
      EPH=ROH*EPS
      AJ(1)=FJ1+ONE
      DO 60 J=2,JM
      J1=J-1
      INCR(J1)=0
      NRED(J1)=0
      FJ=DFLOAT(NJ(J))
      V=AJ(J1)+FJ
      AJ(J)=V
      DO 61 K=1,J1
      W=FJ/DFLOAT(NJ(K))
61    D(J,K)=W*W
      IF(J.EQ.2) GOTO 60
      W=V-FJ1
      DO 62 K1=2,J1
      K=K1-1
      U=(AJ(K1)-V)/(W*DFLOAT(K+K1))
      U=EPH**U
62    AL(J1,K)=U
60    CONTINUE
      COSTF=ONE
      COSTJ=10.*COSTF
      COSTS=ONE
      COSTLR=4.*COSTS
      IF((COSTS+COSTLR+COSTJ).EQ.ZERO) GOTO 64
      AJ(1)=COSTJ+COSTLR+(COSTF+COSTS)*(FJ1+ONE)
      DO 63 J=2,JM
63    AJ(J)=AJ(J-1)+(COSTF+COSTS)*DFLOAT(NJ(J))+COSTS+COSTLR
64    KOH=1
      JOH=2
65    IF(JOH.GE.JM) GOTO 66
      IF(AJ(JOH+1)*ONE1.GT.AJ(JOH)*AL(JOH,KOH)) GOTO 66
      KOH=JOH
      JOH=JOH+1
      GOTO 65
66    K=0
      KM=KOH
      JM=JOH
      INCR(JM)=-1
      OMJO=ZERO
      IF(PRM1) WRITE(MOUT,9100) EPS
      IF(PRM2) WRITE(MOUT,9200) EPS
      IF(PRM2) WRITE(MOUT,9201) KM
      IF(PRM2) WRITE(MOUT,9202) CON1
      IF(PRM2) WRITE(MOUT,9203) CON2
      IF(HMAX.LT.TP(ITPM)-TP(1).AND.PRM) WRITE(MOUT,9001) HMAX
      HMAXU=HMAX
      HMIN=TWO*(TP(ITPM)-TP(1))/DFLOAT(ISMAX)
      FCM=ONE
      IF(HMAXU.GT.HMIN) GOTO 2
      HMAXU=HMIN
      IF(PRERR) WRITE(EOUT,9002) HMAXU
2     EMIN=EPMACH*1.D3
      IF(EPS.GT.EMIN) GOTO 4
      EPS=EMIN
      IF(PRERR) WRITE(EOUT,9004) EPS
4     IF(CON1.GT.SMALL) GOTO 6
      CON1=SMALL
      IF(PRERR) WRITE(EOUT,9006) CON1
6     IF(CON2.GT.EPMACH) GOTO 8
      CON2=EPMACH
      IF(PRERR) WRITE(EOUT,9008) CON2
8     CONTINUE
      IF(PRM1) WRITE(MOUT,9110)
      IF(PRM2) WRITE(MOUT,9210)
      DO 51 I=1,NCEQ
      R(I)=ZERO
51    RINT(I)=ZERO
      YMAX=ZERO
      DO 55 I=1,N
      DO 57 J=1,JM
57    DT(I,J)=ZERO
55    IF(Y(I).GT.YMAX) YMAX=Y(I)
      IF(YMAX.EQ.ZERO) YMAX=ONE
      YMAXH=YMAX*CON1
      SCALET=CON1*YMAX
      IF(SCALET.GT.SMALL) SCALET=SMALL
      DO 59 I=1,N
      U=Y(I)
      IF(U.LT.ZERO) GOTO 33
      IF(U.GE.SCALET) GOTO 58
      SM(I)=YMAX
      SH(I)=YMAXH
      IS(I)=0
      GOTO 59
58    SM(I)=U
      SH(I)=U
      IS(I)=1
59    CONTINUE
      IF(HI.GT.ZERO) GOTO 79
      DM=ZERO
      CALL JACOBI (N,NCEQ,LDIM,TP(1),Y,A)
      IMAX=0
      DO 75 I=1,N
      SHI=ONE/SH(I)
      IMIN=IMAX+1
      IMAX=IA(I+1)-1
      V=ZERO
      DO 70 L=IMIN,IMAX
      IH=JA(L)
      U=SH(IH)*A(L)*SHI
      V=V+DABS(U)
70    CONTINUE
75    IF(V.GT.DM) DM=V
      IF(DM.EQ.ZERO) DM=ONE
      HI=HALF/DM
79    H=HI
      ITP=2
      T=TP(1)
      TEND=TP(2)
      TOLD=T
      H1=TEND-T
      WINDOW=.FALSE.
      LPRINT=.TRUE.
401   CONTINUE
      CALL FCN (N,NCEQ,T,Y,DZ)
      NFCN=NFCN+1
      IF(PRM1) WRITE(MOUT,9120) ISTEP,T
      IF(PRM2) WRITE(MOUT,9220) ISTEP,NFCN,T,K,KOH
      IF(H1.GT.T*EPMACH) GOTO 403
      LPRINT=.TRUE.
      H=HR
      HMAX=HMAXU
      IF(ITP.LT.ITPM) GOTO 405
      CALL SIMDAT (N,NCEQ,Y,TP,ITPM,T,TOLD,IPRINT,LPRINT)
      NWTMP=NZAMAX
      XKAACT=DFLOAT(NZAMAX)/DFLOAT(LDIM)
      IF(ICNCP.GT.10.AND.PRM2) WRITE(MOUT,9310)
      IF(ICNCP.GT.3) GOTO 1000
      NFREE=LICN-NZAMAX
      NFREEH=LICN-(LDIM*14)/10
      IF(NFREE.GT.NFREEH) NFREE=NFREEH
      NFREEH=LICN-(MINICN*12)/10
      IF(NFREE.GT.NFREEH) NFREE=NFREEH
      NFREE=(NFREE*9)/10
      IF(NFREE.LE.10) GOTO 1000
      IF(PRM2) WRITE(MOUT,9320) NFREE
1000  NSTEP=NSTEP+ISTEP
      RETURN
405   ITP=ITP+1
      TEND=TP(ITP)
      H1=TEND-T
      IF(H1.GT.ZERO) GOTO 404
      IF(ITP.LT.ITPM) GOTO 405
      GOTO 35
404   WINDOW=.FALSE.
403   CALL SIMDAT (N,NCEQ,Y,TP,ITPM,T,TOLD,IPRINT,LPRINT)
      LPRINT=.FALSE.
      IF(H1.GE.ONE2*H) GO TO 402
      HR=H
      H=H1
402   JRED=0
      NSTC=0
      DMH=SAFEDM
      DO 409 K=1,KM
409   INCR(K)=INCR(K)+1
      HMAX=H1
      IF(HMAXU.LT.HMAX) HMAX=HMAXU
      CALL JACOBI (N,NCEQ,LDIM,T,Y,A)
      IMAX=0
      DO 100 I=1,N
      SHI=ONE/SH(I)
      IH=IDA(I)
      DZ(I)=DZ(I)*SHI
      IMIN=IMAX+1
      IMAX=IA(I+1)-1
      DO 100 L=IMIN,IMAX
      IH=JA(L)
100   A(L)=-SH(IH)*A(L)*SHI
   10 TN=T+H
      FCM=H/HMAX
      IF(FCM.LT.FMIN) FCM=FMIN
      DM=DMH*TWO
      GH=ZERO
      DO 260 J=1,JM
      LRATE=J.EQ.2
      M=NJ(J)
      G=H/DFLOAT(M)
      GD=GH
      GH=ONE/G
      GD=GH-GD
      DO 2101 I=1,N
      IH=IDA(I)
      DEL(I)=DZ(I)
2101  YM(I)=Y(I)
      IF(.NOT.BANFA) GOTO 2150
      IDISP2=LICN-LDIM
      IDISP(2)=IDISP2+1
      IMAX=0
      DO 2182 I=1,N
      IMIN=IMAX+1
      IMAX=IA(I+1)
      IKEEP(I)=IMAX-IMIN
      IMAX=IMAX-1
      IKEEP(N+I)=I
      IKEEP(N2+I)=I
      DO 2184 II=IMIN,IMAX
      IH=IDISP2+II
      ILEQ(IH)=JA(II)
      WTMP(IH)=A(II)
      IF(JA(II).EQ.I) WTMP(IH)=WTMP(IH)+GH
2184  CONTINUE
2182  CONTINUE
      IDISP2=IDISP2+1
      CALL MA30LA (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),IDISP,IKEEP(NP1),
     @            IKEEP(NP2),ILEQ(LICNP1),LIRN,IW1,IW1(NP1),IW1(NP2),
     @            IW1(NP3),IW1(NP4),IW1(NP5),IWO,IWO(NP1),THRESH,IFLAG)
      IF(IFLAG.LT.0) GOTO 2602
      IELBO2=LICN-IDISP(2)
      CALL MA30LM (N,NP1,IDISP(2),ILEQ,IA,JA,LDIM,ILEQ(LICNP1),IKEEP,
     @            IKEEP(NP1),IKEEP(NP2),IWO)
      NANFA=NANFA+1
      IDISP2=IDISP(2)
      IF(IDISP2.GT.NZAMAX) NZAMAX=IDISP2
      IF(IFLAG.LT.0) GOTO 2602
      BANFA=.FALSE.
      GOTO 2199
2150  CONTINUE
      DO 2155 I=1,IDISP2
2155  WTMP(I)=ZERO
      DO 2156 I=1,LDIM
      II=LICN+I
      IH=ILEQ(II)
2156  WTMP(IH)=A(I)
      DO 21560 I=1,N
      IDI=IDA(I)
      IH=ILEQ(LICN+IDI)
21560 WTMP(IH)=WTMP(IH)+GH
      IF(J.EQ.1) DEPS=HALF
      CALL MA30LB (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),IDISP,IKEEP(NP1),
     @            IKEEP(NP2),WTMP(LICNP1),IWO,IFLAG)
      IF(DEPS.NE.TWO) NTFAC=NTFAC+1
      IF(DEPS.EQ.TWO) NFAC=NFAC+1
      DEPS=TWO
      IF(IFLAG.LT.0) GOTO 2160
      IF(J.NE.1.OR.RMIN.GE.TGROW) GOTO 2175
2160  BANFA=.TRUE.
      GOTO 10
2175  CONTINUE
2199  CONTINUE
      CALL MA30LC (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),IDISP,IKEEP(NP1),
     @            IKEEP(NP2),DEL,WTMP(IDISP2+1))
      NSOL=NSOL+1
      M=M-1
      IF(.NOT.LRATE) GOTO 501
      HH=G*HALF
      DO 500 I=1,NCEQ
500   RTMP(I)=HH*R(I)
501   DO 220 K=1,M
      DIFF=ZERO
      IDIFF=1
      DO 2201 I=1,N
      IF(IS(I).LE.0) GOTO 2201
      DIFF=DIFF+DEL(I)**2
      IDIFF=IDIFF+1
2201  YM(I)=YM(I)+DEL(I)*SH(I)
      DIFF=DSQRT(DIFF/DFLOAT(IDIFF))
      IF(DIFF.GT.TEN) GOTO 2603
      TH=T+DFLOAT(K)*G
      CALL FCN (N,NCEQ,TH,YM,DY)
      NFCN=NFCN+1
      DO 2202 I=1,N
2202  DY(I)=DY(I)/SH(I)-DEL(I)*GH
      CALL MA30LC (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),IDISP,IKEEP(NP1),
     @            IKEEP(NP2),DY,WTMP(IDISP2+1))
      NSOL=NSOL+1
      DO 2203 I=1,N
2203  DEL(I)=DEL(I)+TWO*DY(I)
      IF(.NOT.LRATE) GOTO 220
      DO 502 I=1,NCEQ
502   RTMP(I)=RTMP(I)+G*R(I)
220   CONTINUE
      DO 2204 I=1,N
2204  YM(I)=YM(I)+DEL(I)*SH(I)
      CALL FCN (N,NCEQ,TN,YM,DY)
      NFCN=NFCN+1
      DO 2205 I=1,N
2205  DEL(I)=DY(I)/SH(I)-DEL(I)*GH
      CALL MA30LC (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),IDISP,IKEEP(NP1),
     @            IKEEP(NP2),DEL,WTMP(IDISP2+1))
      NSOL=NSOL+1
      DMA=DM
      DM=ZERO
      DO 2206 I=1,N
      U=DEL(I)*SH(I)
      DEL(I)=U
      YM(I)=YM(I)+U
      IF(J.GT.3) GOTO 2206
      U=DABS(U)
      V=DABS(YM(I))
      IF(V.LT.SH(I)) V=SH(I)
      U=U/V
      IF(U.GT.DM) DM=U
2206  CONTINUE
      IF(J.GT.3) GOTO 2209
      IF(DM.LT.DMA*HALF) GOTO 2207
      IF(DM.LE.TENTH.AND.J.GT.1) GOTO 2209
      NSTC=NSTC+1
      GOTO 2601
2207  IF(J.EQ.1) DMH=DM*HALF
      IF(DM.LT.HALF) GOTO 2209
      DO 2208 L=JOH,JM
      IF(INCR(L).GT.0) INCR(L)=0
      INCR(L)=INCR(L)-2
2208  CONTINUE
2209  IF(.NOT.LRATE) GOTO 504
      DO 503 I=1,NCEQ
503   RTMP(I)=RTMP(I)+HH*R(I)
504   ERR=ZERO
      DO 234 I=1,N
      C=YM(I)
      V=DT(I,1)
      DT(I,1)=C
      IF(J.EQ.1) GOTO 234
      TA=C
      DO 231 K=2,J
      JK=J-K+1
      B1=D(J,JK)
      W=C-V
      U=W/(B1-ONE)
      C=B1*U
      V=DT(I,K)
      DT(I,K)=U
231   TA=U+TA
      YM(I)=TA
      TA=DABS(TA)
      CONTRA=CON1
      IF(IS(I).GE.1) CONTRA=CON2
      SHH=SM(I)*CONTRA
      IF(TA.LT.SHH) TA=SHH
      IF(TA.LT.SMALL) TA=SMALL
      U=U/TA
      ERR=ERR+U*U
234   CONTINUE
      IF(J.EQ.1) GOTO 260
      ERR=DSQRT(ERR/FN)
      KONV=0
      IF(ERR.LT.EPS) KONV=1
      ERR=ERR/EPH
      K=J-1
      FC=ERR**(ONE/DFLOAT(K+J))
      IF(FC.LT.FCM) FC=FCM
      OMJ=FC*AJ(J)
      IF(J.GT.2.AND.OMJ*ONE1.GT.OMJO.OR.K.GT.JOH) GOTO 235
      KO=K
      JO=J
      OMJO=OMJ
      FCO=FC
235   IF(J.LT.KOH.AND.WINDOW) GOTO 260
      IF(KONV.EQ.0) GOTO 236
      IF(KO.LT.K.OR.INCR(J).LT.0) GOTO 20
      IF(NRED(KO).GT.0) NRED(KO)=NRED(KO)-1
      FC=FCO/AL(J,K)
      IF(FC.LT.FCM) FC=FCM
      J1=J+1
      IF(AJ(J+1)*FC*ONE1.GT.OMJO) GOTO 20
      FCO=FC
      KO=JO
      JO=JO+1
      GOTO 20
236   RED=ONE/FCO
      JK=KM
      IF(JOH.LT.KM) JK=JOH
      IF(K.GE.JK) GOTO 239
      IF(KO.LT.KOH) RED=AL(KOH,KO)*RED
      IF(AL(JK,KO).LT.FCO) GOTO 239
260   CONTINUE
239   RED=RED*SAFE
      H=H*RED
2392  IF(.NOT.WINDOW) GOTO 2390
      NRED(KOH)=NRED(KOH)+1
      DO 2391 L=KOH,KM
2391  INCR(L)=-2-NRED(KOH)
2390  JRED=JRED+1
      IF(PRM2) WRITE(MOUT,9015) JRED,RED,KOH
      IF(JRED.GT.JRMAX) GOTO 34
      GOTO 10
2601  HMAX=G*FJ1*HALF
      RED=HMAX/H
      H=HMAX
      IF(PRM2) WRITE(MOUT,9010)
      IF(JRED.GT.0.OR.NSTC.GT.0) GOTO 2390
      GOTO 2392
2602  HMAX=G*FJ1*QUART
      RED=HMAX/H
      H=HMAX
      IF(PRM2) WRITE(MOUT,9011)
      IF(JRED.LE.1) GOTO 2392
      IF(IFLAG.GT.-3) GOTO 2392
      IF(IFLAG.NE.-3) GOTO 37
      IELBOW=MINIRN+2*N
      LIRN=LDIM+IELBOW
      LICN=NLEQ-LIRN
      LICNR=NWTMP-N
      IF(LICNR.LT.LICN) LICN=LICNR
      XKAACT=DFLOAT(LICN)/DFLOAT(LDIM)
      IF(XKAACT.LT.XKAMIN) GOTO 37
      LICNP1=LICN+1
      IDISP(1)=1
      BANFA=.TRUE.
      DEPS=TWO
      GOTO 2392
2603  HMAX=G*FJ1*QUART
      RED=HMAX/H
      H=HMAX
      IF(PRM2) WRITE(MOUT,9012)
      IF(JRED.GT.0) GOTO 2390
      GOTO 2392
20    TOLD=T
      T=TN
      H1=TEND-T
      ISTEP=ISTEP+1
      WINDOW=.TRUE.
      DO 510 I=1,NCEQ
510   RINT(I)=RINT(I)+RTMP(I)
      DO 2606 I=1,N
      U=YM(I)
      Y(I)=U
      IF(U.LT.ZERO) Y(I)=ZERO
      U=DABS(U)
      IF(U.GT.SH(I)) GOTO 2607
      IF(IS(I).LE.0) GOTO 2606
      IF(IS(I).EQ.2) GOTO 2608
      IF(IS(I).EQ.3) GOTO 2605
      IF(U.GE.SM(I)*CON2) GOTO 2609
      IS(I)=2
2608  IF(U.LT.SM(I)*CON2*EPS) IS(I)=3
      GOTO 2606
2605  IF(U.LT.SM(I)*CON2) GOTO 2606
      IS(I)=2
      GOTO 2606
2607  IF(IS(I).EQ.3 .AND. PRM1) WRITE(MOUT,9020) I
      IF(U.GT.SM(I).OR.IS(I).LE.0) SM(I)=U
2609  SH(I)=U
      IS(I)=1
2606  CONTINUE
      IF(ISTEP.GT.ISMAX) GOTO 31
      H=H/FCO
      KOH=KO
      JOH=KOH+1
      IF(H.GT.T*EPMACH) GO TO 401
      IF(PRERR) WRITE(EOUT,9030) H,T
      KFLAG=-1
      GOTO 39
31    IF(PRERR) WRITE(EOUT,9031) ISMAX
      KFLAG=-2
      GOTO 39
32    IF(PRERR) WRITE(EOUT,9032)
      KFLAG=-3
      T=TP(1)
      GOTO 39
33    IF(PRERR) WRITE(EOUT,9033)
      KFLAG=-4
      GOTO 39
34    IF(PRERR) WRITE(EOUT,9034) JRMAX
      KFLAG=-5
      GOTO 39
35    IF(PRERR) WRITE(EOUT,9035) (TP(I),I=1,ITPM)
      KFLAG=-6
      GOTO 39
37    IF(PRERR) WRITE(EOUT,9032)
      KFLAG=-3
      GOTO 39
39    IF(PRERR) WRITE(EOUT,9039)
      NSTEP=NSTEP+ISTEP
      TP(ITPM)=T
      HMAX=HMAXU
      RETURN
9000  FORMAT(///,12H M E T A S C,/,12H -----------)
9001  FORMAT(40H USER PRESCRIBED MAXIMUM STEPSIZE: HMAX=,D10.3,/)
9002  FORMAT(/,43H PRESCRIBED STEPSIZE HMAX TOO SMALL; SET TO,D10.3)
9004  FORMAT(/,42H  PRESCR. REL. PREC. EPS TOO SMALL; SET TO,D10.3)
9006  FORMAT(/,27H CON1 TOO SMALL; REDUCED TO,D10.3)
9008  FORMAT(/,27H CON2 TOO SMALL; REDUCED TO,D10.3)
9010  FORMAT(/,22H STABILITY TEST FAILED)
9011  FORMAT(/,24H LU-DECOMPOSITION FAILED)
9012  FORMAT(/,17H UNBOUNDED GROWTH)
9015  FORMAT(19H STEPSIZE REDUCTION,I7,13H.RED. FACTOR:,D10.3,I5)
9020  FORMAT(/,15H COMPONENT NO.:,I4,28H WHICH HAS POSSIBLY NO RIGHT,
     @         23H DIGIT BEGINNS TO ARISE,/,
     @         50H IF THE O.D.E. IS NOT STABLE, RESULTS MAY BE WRONG,/,
     @         50H PROBLEM SEEMS TO BE OSCILLARITY,   REDUCE CON2   )
9030  FORMAT(/,23H STEPSIZE TOO SMALL; H=,D18.10,6H AT T=,D18.10)
9031  FORMAT(/,17H MORE THAN ISMAX=,I5,12H BASIC STEPS)
9032  FORMAT(/,38H NOT ENOUGH SPACE FOR LU-DECOMPOSITION,/,
     @         36H ENLARGE INTEGER AND REAL WORK SPACE)
9033  FORMAT(/,38H NEGATIVE INITIAL CONCENTRATION GIVEN:,D21.12)
9034  FORMAT(/,17H MORE THAN JRMAX=,I2,20H STEPSIZE REDUCTIONS,
     @         15H PER BASIC STEP)
9035  FORMAT(/,34H PRESCRIBED OUTPUT POINTS INVALID:,17(/,3D21.12))
9039  FORMAT(///,19H INTEGRATION FAILED,/)
9100  FORMAT(/,36H PRESCRIBED RELATIVE PRECISION: EPS=,D10.3,/)
9110  FORMAT(//,8X,4HSTEP,13X,1HT)
9120  FORMAT(/,I10,7X,5H   T=,D18.10)
9200  FORMAT(/,36H PRESCRIBED RELATIVE PRECISION: EPS=,D10.3)
9201  FORMAT(35H MAXIMUM COLUMN NUMBER:         KM=,I2)
9202  FORMAT(37H SCALING PARAMETERS:            CON1=,D10.3)
9203  FORMAT(32X,5HCON2=,D10.3,/)
9210  FORMAT(//,8X,4HSTEP,6X,3HFCN,10X,1HT,20X,7HK   KOH)
9220  FORMAT(/,2I10,7H     T=,D18.10,3X,2I5)
9300  FORMAT(/,41H REAL AND INTEGER WORK SPACE NOT BALANCED,/,
     @37H INTEGER WORK SPACE MAY BE REDUCED BY,I6,18H STORAGE LOCATIONS)
9305  FORMAT(/,41H REAL AND INTEGER WORK SPACE NOT BALANCED,/,
     @34H REAL WORK SPACE MAY BE REDUCED BY,I6,18H STORAGE LOCATIONS)
9310  FORMAT(/,47H REAL AND INTEGER WORK SPACE MAY BE ENLARGED TO,
     @24H SPEED UP LINEAR ALGEBRA)
9320  FORMAT(/,46H REAL AND INTEGER WORK SPACE MAY BE REDUCED BY,I6,
     @18H STORAGE LOCATIONS)
      END
      SUBROUTINE SIMDAT (NSP,NCEQ,CONC,TOUT,IOUT,T,TOLD,IPRINT,LPRINT)
      DOUBLE PRECISION CONC(NSP),TOUT(IOUT),P,RK,T,TOLD,ZERO
      REAL R,RINT,RTMP
      INTEGER ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT,INFDAT
      LOGICAL LPRINT
      COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT
      COMMON / INFORM / IHELP,INFDAT
      COMMON / KINPAR / RK(900)
      COMMON / RATE1 / R(900)
      COMMON / RATE2 / RINT(900)
      COMMON / RATE3 / RTMP(900)
      DATA ZERO/0.D0/
      WRITE(DOUT,9001)
      WRITE(DOUT,9005) T
      WRITE(DOUT,9010) (CONC(I),I=1,NSP)
      WRITE(DOUT,9000)
      P=ZERO
      DO 11 I=1,NSP
11    P=P+CONC(I)
      WRITE(DOUT,9100)
      WRITE(DOUT,9105) T,P
      WRITE(DOUT,9000)
      IF(IPRINT.EQ.0) RETURN
      IF(IPRINT.LE.3.AND..NOT.LPRINT) RETURN
      IPR=IPRINT
      IF(IPRINT.GT.3) IPR=IPR-3
      IF(IPR.NE.1.AND.IPR.NE.3) GOTO 20
      WRITE(DOUT,9400)
      WRITE(DOUT,9405) T
      WRITE(DOUT,9010) (RINT(I),I=1,NCEQ)
      WRITE(DOUT,9000)
20    IF(IPR.EQ.1) RETURN
      WRITE(DOUT,9300)
      WRITE(DOUT,9305) T
      WRITE(DOUT,9010) (R(I),I=1,NCEQ)
      WRITE(DOUT,9000)
      RETURN
9000  FORMAT(5H &END)
9001  FORMAT(6H &CONC)
9005  FORMAT(3H T= ,D21.12,4H, C=)
9010  FORMAT(5(1X,D13.7,1H,))
9100  FORMAT(6H &SUMC)
9105  FORMAT(3H T= ,D21.12,6H, SUM=,D13.7)
9300  FORMAT(6H &RATE)
9305  FORMAT(3H T= ,D21.12,5H, RK=)
9400  FORMAT(7H &IRATE)
9405  FORMAT(3H T= ,D21.12,5H, RK=)
      END
      SUBROUTINE FCN (NSP,NCEQ,T,Y,DY)
      DOUBLE PRECISION Y(NSP),DY(NSP),T,RT,RK
      INTEGER*2 LR,PLR
      REAL R
      COMMON / LARK1 / LR(4000)
      COMMON / LARK2 / PLR(1800)
      COMMON / KINPAR / RK(900)
      COMMON / RATE1 / R(900)
      DO 1 I=1,NSP
      DY(I)=0.D0
1     CONTINUE
      IR2=0
      NK=NCEQ
      DO 10 K=1,NCEQ
      RT=RK(K)
      IL1=IR2+1
      IL2=PLR(K)
      NK=NK+1
      IR2=PLR(NK)
      IF(IL2.LT.IL1) GOTO 3
      DO 2 I=IL1,IL2
      LRI=LR(I)
2     RT=RT*Y(LRI)
3     CONTINUE
      DO 6 I=IL1,IR2
      LRI=LR(I)
      IF(I.GT.IL2) GOTO 5
      DY(LRI)=DY(LRI) - RT
      GOTO 6
5     DY(LRI)=DY(LRI) + RT
6     CONTINUE
10    R(K)=RT
      RETURN
      END
      SUBROUTINE JACOBI (NSP,NCEQ,LDIM,T,Y,A)
      DOUBLE PRECISION A(LDIM),Y(NSP),T,TERM,RKK,RK
      INTEGER*2 LR,PLR,JA,IDA
      INTEGER IN(10),I1(10)
      COMMON / LARK1 / LR(4000)
      COMMON / LARK2 / PLR(1800)
      COMMON / LARK3 / JA(2800)
      COMMON / LARK5 / IDA(400)
      COMMON / KINPAR / RK(900)
      DO 1 L=1,LDIM
      A(L)=0.D0
1     CONTINUE
      IR2=0
      NK=NCEQ
      DO 100 K=1,NCEQ
      IL1=IR2+1
      IL2=PLR(K)
      IK=IL2-IR2
      NK=NK+1
      IR2=PLR(NK)
      IF(IK.LE.0) GOTO 100
      IM=0
      DO 2 M=IL1,IR2
      IM=IM+1
      LRI=LR(M)
      I1(IM)=IDA(LRI)
2     IN(IM)=LRI
      RKK=RK(K)
      DO 40 J=1,IK
      TERM=RKK
      JJ=IN(J)
      DO 42 I=1,IK
      L=IN(I)
42    IF(J.NE.I) TERM=TERM*Y(L)
      DO 50 I=1,IM
      L=I1(I)
      IF(JA(L)-JJ) 51,54,52
51    L=L+1
      IF(JA(L).NE.JJ) GOTO 51
      GOTO 54
52    L=L-1
      IF(JA(L).NE.JJ) GOTO 52
54    IF(I.GT.IK) GOTO 56
      A(L)=A(L) - TERM
      GOTO 50
56    A(L)=A(L) + TERM
50    CONTINUE
40    CONTINUE
100   CONTINUE
      RETURN
      END
      SUBROUTINE MA30LA(N,ICN,A,LICN,LENR,LENRL,IDISP,IP,IQ,IRN,
     @LIRN,LENC,IFIRST,LASTR,NEXTR,LASTC,NEXTC,IPTR,IPC,UIN,IFLAG)
      DOUBLE PRECISION A(LICN),U,UIN,AU,UMAX,AMAX,ZERO
      DOUBLE PRECISION DABS,DMAX1,DMIN1
      INTEGER PIVOT,PIVEND,DISPC,OLDPIV,OLDEND,PIVROW
      INTEGER ROWI
      INTEGER IDISP(2)
      LOGICAL ABORT1,ABORT2,ABORT3
      INTEGER*2 ICN(LICN),LENR(N),LENRL(N),IP(N),IQ(N),LENC(N)
      INTEGER*2 IRN(LIRN),IPC(N),IPTR(N)
      INTEGER*2 IFIRST(N),LASTR(N),NEXTR(N),LASTC(N),NEXTC(N)
      COMMON /MA30LE/ LP,ABORT1,ABORT2,ABORT3
      COMMON /MA30LF/ IRNCP,ICNCP,IRANK,MINIRN,MINICN
      DATA UMAX/.999999999D0/
      DATA ZERO/0.0D0/
      MINIRN=0
      MINICN=IDISP(1)-1
      MOREI=0
      IRANK=N
      IRNCP=0
      ICNCP=0
      IFLAG=0
      U=DMIN1(UIN,UMAX)
      U=DMAX1(U,ZERO)
      IBEG=IDISP(1)
      IACTIV=IDISP(2)
      NZROW=LICN-IACTIV+1
      MINICN=NZROW+MINICN
      IPTR(1)=IACTIV
      IF (N.EQ.1) GO TO 20
      NM1=N-1
      DO 10 I=1,NM1
   10 IPTR(I+1)=IPTR(I)+LENR(I)
  20  ITOP=LICN
      DO 110 I=1,N
      LENRL(I)=0
  110 LENC(I)=0
      IF (ITOP-IACTIV.LT.LIRN) GO TO 120
      MINIRN=ITOP-IACTIV+1
      PIVOT=0
      GO TO 1050
  120 DO 130 II=IACTIV,ITOP
      I=ICN(II)
  130 LENC(I)=LENC(I)+1
      IPC(N)=LIRN+1
      IF(N.EQ.1) GOTO 142
      DO 140 JJ=2,N
      J=N-JJ+1
  140 IPC(J)=IPC(J+1)-LENC(J+1)
  142 DO 160 INDROW=1,N
      J1=IPTR(INDROW)
      J2=J1+LENR(INDROW)-1
      IF (J1.GT.J2) GO TO 160
      DO 150 JJ=J1,J2
      J=ICN(JJ)
      IPOS=IPC(J)-1
      IRN(IPOS)=INDROW
      IPC(J)=IPOS
  150 CONTINUE
  160 CONTINUE
      DISPC=IPC(1)
      NZCOL=LIRN-DISPC+1
      MINIRN=MAX0(NZCOL,MINIRN)
      NZMIN=1
      DO 170 I=1,N
  170 IFIRST(I)=0
      DO 190 JJ=1,N
      J=N-JJ+1
      NZ=LENC(J)
      IF (NZ.NE.0) GO TO 180
      IPC(J)=0
      LASTC(J)=0
      GO TO 190
  180 ISW=IFIRST(NZ)
      IFIRST(NZ)=-J
      LASTC(J)=0
      NEXTC(J)=-ISW
      ISW1=IABS(ISW)
      IF (ISW.NE.0) LASTC(ISW1)=J
  190 CONTINUE
      DO 210 II=1,N
      I=N-II+1
      NZ=LENR(I)
      IF (NZ.NE.0) GO TO 200
      IPTR(I)=0
      LASTR(I)=0
      GO TO 210
  200 ISW=IFIRST(NZ)
      IFIRST(NZ)=I
      IF (ISW.GT.0) GO TO 205
      NEXTR(I)=0
      LASTR(I)=ISW
      GO TO 210
 205  NEXTR(I)=ISW
      LASTR(I)=LASTR(ISW)
      LASTR(ISW)=I
  210 CONTINUE
      DO 930 PIVOT=1,N
      NZ2=NZMIN
      JCOST=N*N
      DO 290 L=1,2
      LL=L
      DO 280 NZ=NZ2,N
      IF (JCOST.LE.(NZ-1)**2) GO TO 380
      IJFIR=IFIRST(NZ)
      IF (IJFIR) 212,211,215
 211  IF (LL.EQ.1) NZMIN=NZ+1
      GO TO 280
 212  LL=2
      IJFIR=-IJFIR
      GO TO 245
 215  LL=2
      DO 235 IDUMMY=1,N
      IF (IJFIR.EQ.0) GO TO 240
      I=IJFIR
      IJFIR=NEXTR(I)
      AMAX=ZERO
      J1=IPTR(I)+LENRL(I)
      J2=IPTR(I)+LENR(I)-1
      DO 220 JJ=J1,J2
 220  AMAX=DMAX1(AMAX,DABS(A(JJ)))
      AU=AMAX*U
      DO 230 JJ=J1,J2
      IF (DABS(A(JJ)).LE.AU.AND.L.EQ.1) GO TO 230
      J=ICN(JJ)
      KCOST=(NZ-1)*(LENC(J)-1)
      IF (KCOST.GE.JCOST) GO TO 230
      JCOST=KCOST
      IJPOS=JJ
      IPIV=I
      JPIV=J
      IF (JCOST.LE.(NZ-1)**2) GO TO 380
  230 CONTINUE
 235  CONTINUE
 240  IJFIR=IFIRST(NZ)
      IJFIR=-LASTR(IJFIR)
 245  IF (JCOST.LE.NZ*(NZ-1)) GO TO 380
      DO 270 IDUMMY=1,N
      IF (IJFIR.EQ.0) GO TO 280
      J=IJFIR
      IJFIR=NEXTC(IJFIR)
      I1=IPC(J)
      I2=I1+NZ-1
      DO 260 II=I1,I2
      I=IRN(II)
      KCOST=(NZ-1)*(LENR(I)-LENRL(I)-1)
      IF (KCOST.GE.JCOST) GO TO 260
      J1=IPTR(I)+LENRL(I)
      J2=IPTR(I)+LENR(I)-1
      AMAX=ZERO
      DO 250 JJ=J1,J2
      AMAX=DMAX1(AMAX,DABS(A(JJ)))
  250 IF (ICN(JJ).EQ.J) JPOS=JJ
      IF (DABS(A(JPOS)).LE.AMAX*U.AND.L.EQ.1) GO TO 260
      JCOST=KCOST
      IPIV=I
      JPIV=J
      IJPOS=JPOS
      IF (JCOST.LE.NZ*(NZ-1)) GO TO 380
  260 CONTINUE
  270 CONTINUE
  280 CONTINUE
      IRANK=IRANK-1
  290 CONTINUE
      IF (IFLAG.NE.2.AND.IFLAG.NE.-5) IFLAG=1
      IRANK=IRANK-N+PIVOT+1
      IF (.NOT.ABORT1) GO TO 300
      IDISP(2)=IACTIV
      IFLAG=-1
      IF (LP.NE.0) WRITE(LP,50)
 50   FORMAT(65H ERROR RETURN FROM MA30LA BECAUSE MATRIX IS STRUCTURALLY
     @ SINGULAR)
      GO TO 1110
  300 K=PIVOT-1
      DO 350 I=1,N
      IF (LASTR(I).NE.0) GO TO 350
      K=K+1
      LASTR(I)=K
      IF (LENRL(I).EQ.0) GO TO 340
      MINICN=MAX0(MINICN,NZROW+IBEG-1+MOREI+LENRL(I))
      IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 320
      CALL MA30LD(A,ICN,IPTR(1),N,IACTIV,ITOP,.TRUE.)
      IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 320
      MOREI=MOREI+IBEG-IDISP(1)
      IBEG=IDISP(1)
      IF (LP.NE.0) WRITE(LP,310)
      IFLAG=-5
  310 FORMAT(48H LU DECOMPOSITION DESTROYED TO CREATE MORE SPACE)
      IF (ABORT3) GO TO 1030
  320 J1=IPTR(I)
      J2=J1+LENRL(I)-1
      IPTR(I)=0
      DO 330 JJ=J1,J2
      A(IBEG)=A(JJ)
      ICN(IBEG)=ICN(JJ)
      ICN(JJ)=0
  330 IBEG=IBEG+1
      NZROW=NZROW-LENRL(I)
  340 IF (K.EQ.N) GO TO 360
  350 CONTINUE
  360 K=PIVOT-1
      DO 370 I=1,N
      IF (LASTC(I).NE.0) GO TO 370
      K=K+1
      LASTC(I)=-K
      IF (K.EQ.N) GO TO 940
  370 CONTINUE
  380 ISING=PIVOT
      IF (A(IJPOS).NE.ZERO) GO TO 390
      ISING=-ISING
      IF (IFLAG.NE.-5) IFLAG=2
      IF (.NOT.ABORT2) GO TO 390
      IDISP(2)=IACTIV
      IFLAG=-2
      IF (LP.NE.0) WRITE(LP,70)
 70   FORMAT(65H ERROR RETURN FROM MA30LA BECAUSE MATRIX IS NUMERICALLY
     @ SINGULAR)
      GO TO 1110
  390 OLDPIV=IPTR(IPIV)+LENRL(IPIV)
      OLDEND=IPTR(IPIV)+LENR(IPIV)-1
      DO 460 JJ=OLDPIV,OLDEND
      J=ICN(JJ)
      LC=LASTC(J)
      NC=NEXTC(J)
      IF (NC.NE.0) LASTC(NC)=LC
      IF (LC.EQ.0) GO TO 440
      NEXTC(LC)=NC
      GO TO 460
 440  NZ=LENC(J)
      ISW=IFIRST(NZ)
      IF (ISW.GT.0) LASTR(ISW)=-NC
      IF (ISW.LT.0) IFIRST(NZ)=-NC
  460 CONTINUE
      I1=IPC(JPIV)
      I2=I1+LENC(JPIV)-1
      DO 530 II=I1,I2
      I=IRN(II)
      LR=LASTR(I)
      NR=NEXTR(I)
      IF (NR.NE.0) LASTR(NR)=LR
      IF (LR.LE.0) GO TO 500
      NEXTR(LR)=NR
      GO TO 530
 500  NZ=LENR(I)-LENRL(I)
      IF (NR.NE.0) IFIRST(NZ)=NR
      IF (NR.EQ.0) IFIRST(NZ)=LR
  530 CONTINUE
      LASTC(JPIV)=ISING
      LASTR(IPIV)=PIVOT
      IF (OLDPIV.EQ.IJPOS) GO TO 540
      AU=A(OLDPIV)
      A(OLDPIV)=A(IJPOS)
      A(IJPOS)=AU
      ICN(IJPOS)=ICN(OLDPIV)
      ICN(OLDPIV)=JPIV
  540 MINICN=MAX0(MINICN,NZROW+IBEG-1+MOREI+LENR(IPIV))
      IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 550
      CALL MA30LD(A,ICN,IPTR(1),N,IACTIV,ITOP,.TRUE.)
      OLDPIV=IPTR(IPIV)+LENRL(IPIV)
      OLDEND=IPTR(IPIV)+LENR(IPIV)-1
      IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 550
      MOREI=MOREI+IBEG-IDISP(1)
      IBEG=IDISP(1)
      IF (LP.NE.0) WRITE(LP,310)
      IFLAG=-5
      IF (ABORT3) GO TO 1030
      IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 550
      IFLAG=-4
      GO TO 1030
  550 IJPOS=0
      J1=IPTR(IPIV)
      DO 570 JJ=J1,OLDEND
      A(IBEG)=A(JJ)
      ICN(IBEG)=ICN(JJ)
      IF (IJPOS.NE.0) GO TO 560
      IF (ICN(JJ).EQ.JPIV) IJPOS=IBEG
      ICN(JJ)=0
      GO TO 570
  560 K=IBEG-IJPOS
      J=ICN(JJ)
      ICN(JJ)=IQ(J)
      IQ(J)=-K
  570 IBEG=IBEG+1
      IJP1=IJPOS+1
      PIVEND=IBEG-1
      LENPIV=PIVEND-IJPOS
      NZROW=NZROW-LENRL(IPIV)-1
      IPTR(IPIV)=OLDPIV+1
      IF (LENPIV.EQ.0) IPTR(IPIV)=0
      DO 600 JJ=IJPOS,PIVEND
      J=ICN(JJ)
      I1=IPC(J)
      LENC(J)=LENC(J)-1
      I2=IPC(J)+LENC(J)-1
      IF (I2.LT.I1) GO TO 590
      DO 580 II=I1,I2
      IF (IRN(II).NE.IPIV) GO TO 580
      IRN(II)=IRN(I2+1)
      GO TO 590
  580 CONTINUE
  590 IRN(I2+1)=0
  600 CONTINUE
      NZCOL=NZCOL-LENPIV-1
      NZPC=LENC(JPIV)
      IF (NZPC.EQ.0) GO TO 870
      DO 820 III=1,NZPC
      II=IPC(JPIV)+III-1
      I=IRN(II)
      J1=IPTR(I)+LENRL(I)
      IEND=IPTR(I)+LENR(I)-1
      DO 610 JJ=J1,IEND
      IF (ICN(JJ).NE.JPIV) GO TO 610
      AU=ZERO
      IF (A(IJPOS).NE.ZERO) AU=-A(JJ)/A(IJPOS)
      A(JJ)=A(J1)
      A(J1)=AU
      ICN(JJ)=ICN(J1)
      ICN(J1)=JPIV
      LENRL(I)=LENRL(I)+1
      GO TO 620
  610 CONTINUE
  620 IF (LENPIV.EQ.0) GO TO 820
      ROWI=J1+1
      IOP=0
      IF (ROWI.GT.IEND) GO TO 640
      DO 630 JJ=ROWI,IEND
      J=ICN(JJ)
      IF (IQ(J).GT.0) GO TO 630
      IOP=IOP+1
      PIVROW=IJPOS-IQ(J)
      A(JJ)=A(JJ)+AU*A(PIVROW)
      ICN(PIVROW)=-ICN(PIVROW)
  630 CONTINUE
  640 IFILL=LENPIV-IOP
      IF (IFILL.EQ.0) GO TO 740
      MINICN=MAX0(MINICN,MOREI+IBEG-1+NZROW+IFILL+LENR(I))
      DO 650 JDIFF=1,IFILL
      JNPOS=IEND+JDIFF
      IF (JNPOS.GT.LICN) GO TO 660
      IF (ICN(JNPOS).NE.0) GO TO 660
  650 CONTINUE
      IEND=IEND+1
      GO TO 740
  660 JMORE=IFILL-JDIFF+1
      I1=IPTR(I)
      DO 670 JDIFF=1,JMORE
      JNPOS=I1-JDIFF
      IF (JNPOS.LT.IACTIV) GO TO 680
      IF (ICN(JNPOS).NE.0) GO TO 690
  670 CONTINUE
  680 JNPOS=I1-JMORE
      GO TO 700
  690 JNPOS=IACTIV-LENR(I)-IFILL
  700 IF (JNPOS.GE.IBEG) GO TO 720
      CALL MA30LD(A,ICN,IPTR(1),N,IACTIV,ITOP,.TRUE.)
      I1=IPTR(I)
      IEND=I1+LENR(I)-1
      JNPOS=IACTIV-LENR(I)-IFILL
      IF (JNPOS.GE.IBEG) GO TO 720
      MOREI=MOREI+IBEG-IDISP(1)-LENPIV-1
      IF (LP.NE.0) WRITE(LP,310)
      IFLAG=-5
      IF (ABORT3) GO TO 1030
      IBEG=IDISP(1)
      ICN(IBEG)=JPIV
      A(IBEG)=A(IJPOS)
      IJPOS=IBEG
      DO 710 JJ=IJP1,PIVEND
      IBEG=IBEG+1
      A(IBEG)=A(JJ)
  710 ICN(IBEG)=ICN(JJ)
      IJP1=IJPOS+1
      PIVEND=IBEG
      IBEG=IBEG+1
      IF (JNPOS.GE.IBEG) GO TO 720
      IFLAG=-4
      GO TO 1030
  720 IACTIV=MIN0(IACTIV,JNPOS)
      IPTR(I)=JNPOS
      DO 730 JJ=I1,IEND
      A(JNPOS)=A(JJ)
      ICN(JNPOS)=ICN(JJ)
      JNPOS=JNPOS+1
  730 ICN(JJ)=0
      IEND=JNPOS
  740 NZROW=NZROW+IFILL
      DO 810 JJ=IJP1,PIVEND
      J=ICN(JJ)
      IF (J.LT.0) GO TO 800
      A(IEND)=AU*A(JJ)
      ICN(IEND)=J
      IEND=IEND+1
      MINIRN=MAX0(MINIRN,NZCOL+LENC(J)+1)
      JEND=IPC(J)+LENC(J)
      JROOM=NZPC-III+1+LENC(J)
      IF (JEND.GT.LIRN) GO TO 750
      IF (IRN(JEND).EQ.0) GO TO 790
  750 IF (JROOM.LT.DISPC) GO TO 760
      CALL MA30LD(A,IRN,IPC(1),N,DISPC,LIRN,.FALSE.)
      IF (JROOM.LT.DISPC) GO TO 760
      JROOM=DISPC-1
      IF (JROOM.GE.LENC(J)+1) GO TO 760
      GO TO 1050
  760 JBEG=IPC(J)
      JEND=IPC(J)+LENC(J)-1
      JZERO=DISPC-1
      DISPC=DISPC-JROOM
      IDISPC=DISPC
      DO 770 II=JBEG,JEND
      IRN(IDISPC)=IRN(II)
      IRN(II)=0
  770 IDISPC=IDISPC+1
      IPC(J)=DISPC
      JEND=IDISPC
      DO 780 II=JEND,JZERO
  780 IRN(II)=0
  790 IRN(JEND)=I
      NZCOL=NZCOL+1
      LENC(J)=LENC(J)+1
      GO TO 810
  800 ICN(JJ)=-J
  810 CONTINUE
      LENR(I)=LENR(I)+IFILL
  820 CONTINUE
      I1=IPC(JPIV)
      I2=IPC(JPIV)+LENC(JPIV)-1
      NZCOL=NZCOL-LENC(JPIV)
      DO 860 II=I1,I2
      I=IRN(II)
      IRN(II)=0
      NZ=LENR(I)-LENRL(I)
      IF (NZ.NE.0) GO TO 830
      LASTR(I)=0
      GO TO 860
  830 IFIR=IFIRST(NZ)
      IFIRST(NZ)=I
      IF (IFIR) 840,855,850
 840  LASTR(I)=IFIR
      NEXTR(I)=0
      GO TO 860
 850  LASTR(I)=LASTR(IFIR)
      NEXTR(I)=IFIR
      LASTR(IFIR)=I
      GO TO 860
 855  LASTR(I)=0
      NEXTR(I)=0
      NZMIN=MIN0(NZMIN,NZ)
 860  CONTINUE
  870 IPC(JPIV)=0
      IF (LENPIV.EQ.0) GO TO 930
      NZROW=NZROW-LENPIV
      JVAL=IJP1
      JZER=IPTR(IPIV)
      IPTR(IPIV)=0
      DO 880 JCOUNT=1,LENPIV
      J=ICN(JVAL)
      IQ(J)=ICN(JZER)
      ICN(JZER)=0
      JVAL=JVAL+1
  880 JZER=JZER+1
      DO 920 JJ=IJP1,PIVEND
      J=ICN(JJ)
      NZ=LENC(J)
      IF (NZ.NE.0) GO TO 890
      LASTC(J)=0
      GO TO 920
  890 IFIR=IFIRST(NZ)
      LASTC(J)=0
      IF (IFIR) 900,910,915
 900  IFIRST(NZ)=-J
      IFIR=-IFIR
      LASTC(IFIR)=J
      NEXTC(J)=IFIR
      GO TO 920
 910  IFIRST(NZ)=-J
      NEXTC(J)=0
      NZMIN=MIN0(NZMIN,NZ)
      GO TO 920
 915  LC=-LASTR(IFIR)
      LASTR(IFIR)=-J
      NEXTC(J)=LC
      IF (LC.NE.0) LASTC(LC)=J
  920 CONTINUE
  930 CONTINUE
  940 IF (N.NE.N) IACTIV=IPTR(N+1)
  950 CONTINUE
      IEND=IBEG-1
      DO 980 JJ=1,IEND
      JOLD=ICN(JJ)
  980 ICN(JJ)=LASTC(JOLD)
      DO 990 II=1,N
      I=LASTR(II)
      NEXTR(I)=LENR(II)
  990 NEXTC(I)=LENRL(II)
      DO 1000 I=1,N
      LENRL(I)=NEXTC(I)
 1000 LENR(I)=NEXTR(I)
      DO 1010 II=1,N
      I=LASTR(II)
      J=LASTC(II)
      NEXTR(I)=IP(II)
 1010 NEXTC(J)=IQ(II)
      DO 1020 I=1,N
      IP(I)=NEXTR(I)
 1020 IQ(I)=NEXTC(I)
      IDISP(2)=IEND
      GO TO 1110
 1030 IDISP(2)=IACTIV
      IF (LP.EQ.0) GO TO 1110
      WRITE(LP,1040)
 1040 FORMAT(53H ERROR RETURN FROM MA30LA BECAUSE LICN NOT BIG ENOUGH)
      GO TO 1080
 1050 IF (IFLAG.EQ.-5) IFLAG=-6
      IF (IFLAG.NE.-6) IFLAG=-3
      IDISP(2)=IACTIV
      IF (LP.EQ.0) GO TO 1110
      IF (IFLAG.EQ.-3) WRITE(LP,1060)
      IF (IFLAG.EQ.-6) WRITE(LP,1070)
 1060 FORMAT(53H ERROR RETURN FROM MA30LA BECAUSE LIRN NOT BIG ENOUGH)
 1070 FORMAT(49H ERROR RETURN FROM MA30LA LIRN AND LICN TOO SMALL)
 1080 WRITE(LP,1090) PIVOT
 1090 FORMAT(10H AT STAGE ,I5)
      IF (PIVOT.EQ.0) WRITE(LP,1100) MINIRN
 1100 FORMAT(34H TO CONTINUE SET LIRN TO AT LEAST ,I8)
 1110 RETURN
      END
      SUBROUTINE MA30LD(A,ICN,IPTR,N,IACTIV,ITOP,REALS)
      DOUBLE PRECISION A(ITOP)
      LOGICAL REALS
      INTEGER*2 ICN(ITOP),IPTR(N)
      COMMON /MA30LF/ IRNCP,ICNCP,IRANK,MINIRN,MINICN
      IF (REALS) ICNCP=ICNCP+1
      IF (.NOT.REALS) IRNCP=IRNCP+1
      DO 10 J=1,N
      K=IPTR(J)
      IF (K.LT.IACTIV) GO TO 10
      IPTR(J)=ICN(K)
      ICN(K)=-J
   10 CONTINUE
      KN=ITOP+1
      KL=ITOP-IACTIV+1
      DO 30 K=1,KL
      JPOS=ITOP-K+1
      IF (ICN(JPOS).EQ.0) GO TO 30
      KN=KN-1
      IF (REALS) A(KN)=A(JPOS)
      IF (ICN(JPOS).GE.0) GO TO 20
      J=-ICN(JPOS)
      ICN(JPOS)=IPTR(J)
      IPTR(J)=KN
   20 ICN(KN)=ICN(JPOS)
   30 CONTINUE
      IACTIV=KN
      RETURN
      END
      SUBROUTINE MA30LB(N,ICN,A,LICN,LENR,LENRL,IDISP,IP,IQ,W,IW,
     @IFLAG)
      DOUBLE PRECISION A(LICN),W(N),AU,EPS,ROWMAX,ZERO,ONE,RMIN
      DOUBLE PRECISION DABS,DMAX1
      INTEGER IDISP(2),PIVPOS
      LOGICAL ABORT1,ABORT2,ABORT3,STAB
      INTEGER*2 ICN(LICN),LENR(N),LENRL(N),IP(N),IQ(N),IW(N)
      COMMON /MA30LE/ LP,ABORT1,ABORT2,ABORT3
      COMMON /MA30LG/ EPS,RMIN
      DATA ZERO/0.0D0/,ONE/1.0D0/
      STAB=EPS.LE.ONE
      RMIN=EPS
      ISING=0
      IFLAG=0
      IF (N.EQ.1) GO TO 170
      DO 10 I=1,N
   10 W(I)=ZERO
      IW(1)=IDISP(1)
      DO 20 I=2,N
   20 IW(I)=IW(I-1)+LENR(I-1)
      DO 140 I=1,N
      ISTART=IW(I)
      IFIN=ISTART+LENR(I)-1
      ILEND=ISTART+LENRL(I)-1
      IF (ISTART.GT.ILEND) GO TO 70
      DO 30 JJ=ISTART,IFIN
      J=ICN(JJ)
   30 W(J)=A(JJ)
      DO 50 JJ=ISTART,ILEND
      J=ICN(JJ)
      IPIVJ=IW(J)+LENRL(J)
      AU=-W(J)/A(IPIVJ)
      W(J)=AU
      IPIVJ=IPIVJ+1
      JFIN=IW(J)+LENR(J)-1
      IF (IPIVJ.GT.JFIN) GO TO 50
      DO 40 JAYJAY=IPIVJ,JFIN
      JAY=ICN(JAYJAY)
   40 W(JAY)=W(JAY)+AU*A(JAYJAY)
   50 CONTINUE
      DO 60 JJ=ISTART,IFIN
      J=ICN(JJ)
      A(JJ)=W(J)
   60 W(J)=ZERO
   70 PIVPOS=ILEND+1
      IF (PIVPOS.GT.IFIN) GO TO 150
      IF (A(PIVPOS).EQ.ZERO) GO TO 150
      IF (.NOT.STAB) GO TO 140
      ROWMAX=ZERO
      DO 130 JJ=PIVPOS,IFIN
  130 ROWMAX=DMAX1(ROWMAX,DABS(A(JJ)))
      IF (DABS(A(PIVPOS))/ROWMAX.GE.RMIN) GO TO 140
      IFLAG=I
      RMIN=DABS(A(PIVPOS))/ROWMAX
  140 CONTINUE
      GO TO 170
  150 IF (LP.NE.0) WRITE(LP,160) I
  160 FORMAT(53H ERROR RETURN FROM MA30LB SINGULARITY DETECTED IN ROW,
     1I8)
      IFLAG=-I
  170 RETURN
      END
      SUBROUTINE MA30LC(N,ICN,A,LICN,LENR,LENRL,IDISP,IP,IQ,X,W)
      DOUBLE PRECISION A(LICN),X(N),W(N),WII,WI,RESID,ZERO
      INTEGER IDISP(2)
      INTEGER*2 ICN(LICN),LENR(N),LENRL(N),IP(N),IQ(N)
      COMMON /MA30LH/ RESID
      DATA ZERO/0.0D0/
      RESID=ZERO
      DO 10 II=1,N
   10 W(II)=X(IP(II))
      IBLOCK=1
      DO 120 I=1,N
      IF (LENRL(I).EQ.0) GO TO 50
      WI=W(I)
      IEND=IBLOCK+LENRL(I)-1
      DO 40 JJ=IBLOCK,IEND
      J=ICN(JJ)
   40 WI=WI+A(JJ)*W(J)
      W(I)=WI
   50 IBLOCK=IBLOCK+LENR(I)
 120  CONTINUE
      J1=IBLOCK
      DO 100 III=1,N
      II=N-III+1
      J2=J1-1
      J1=J1-LENR(II)
      JPIV=J1+LENRL(II)
      JPIVP1=JPIV+1
      IF (J2.LT.JPIVP1) GO TO 90
      WII=W(II)
      DO 80 JJ=JPIVP1,J2
      J=ICN(JJ)
   80 WII=WII-A(JJ)*W(J)
      W(II)=WII
   90 W(II)=W(II)/A(JPIV)
  100 CONTINUE
      DO 130 II=1,N
      I=IQ(II)
  130 X(I)=W(II)
      RETURN
      END
      SUBROUTINE MA30LM (N,NP1,NZLU,ICN,IA,JA,NZ,IMAP,LENR,IP,IQ,IW)
      INTEGER*2 IW(N),IQ(N)
      INTEGER*2 ICN(NZLU),IA(NP1),JA(NZ),IMAP(NZ),LENR(N),IP(N)
      J1=1
      DO 100 INEW=1,N
      J2=J1+LENR(INEW)-1
      DO 30 JJ=J1,J2
      JNEW=ICN(JJ)
      JOLD=IQ(JNEW)
30    IW(JOLD)=JJ
      IOLD=IP(INEW)
      JAY1=IA(IOLD)
      JAY2=IA(IOLD+1)-1
      DO 40 JJ=JAY1,JAY2
      J=JA(JJ)
      IMAP(JJ)=IW(J)
40    CONTINUE
      J1=J2+1
100   CONTINUE
      RETURN
      END
 
