Google
 

Trailing-Edge - PDP-10 Archives - decus_20tap3_198111 - decus/20-0079/sinc.for
There are no other files named sinc.for in the archive.
C     PROGRAM SINC (INPUT,OUTPUT,TAPE1=INPUT,TAPE6=OUTPUT)
C
C     .... A SIMULATOR FOR NON-LINEAR ELECTRONIC CIRCUITS
C
      COMMON  YETC(5930)
      COMMON /SPARSE/ IORD(300),NUR(300),NVV(16),IPOS(1600),IUC(600),
     1 IV(200)
      COMMON /NODES/ NREN(100),LNOD,NX,MDIM,NNOD,NUT,NLT,NVT
      COMMON /TRANS/ MREF(100),NQ(100,17),VT,TOLER,NXTR,ITR,IND
      COMMON /MODEL/ NTYP(30),BETA(30,8),CSAT(30,1),RBS(30,3),RCL(30,3),
     1 ROU(30,2),FT(30,4),CJ1(30,4),CJ2(30,6),TEMP(30)
      COMMON /MODPR/ BEF(30),BER(30),CSA(30),TCF(30),TCR(30),GBS(30),
     1 GCL(30),BWF(30),BC1(30),BC2(30),BC3(30),CJE(30),CJC(30),CT1(30),
     2 CT2(30)
      COMMON /ELMTS/ KIND(100),ND(100,4),VALU(100),NELT
      COMMON /VSRCE/ VSRC(10),NAMV(10),TVLT(10,10),VVAL(10,10),
     1 NVVAL(10),KVOL(10),NV(10,2),NVSC
      COMMON /CSRCE/ CSRC(10),NAMC(10),TCUR(10,10),CVAL(10,10),
     1 NCVAL(10),KCUR(10),NCS(10,2),NCSC
      COMMON /TIMEV/ TTIME,DELT,TIMP,DELP,TMAX,TPRT,DPRT,ITER,NCUT,ICONV
      COMMON /NAMES/ NAME(100),NAMQ(100),MODN(100),NAMB(30),NMOD
      COMMON /OUTPT/ KOUT(10),NAMO(10),NODO(10,2),NOUT
      COMMON /PLOTS/ PMIN(6),PMAX(6),KPLO(2,6),NPVAL(6),NAMP(2,6),
     1 NP(2,2,6),NPLT,NPTS,FSTR,FSTP,NHAR
      COMMON /STORE/ LQ(100,3),LD(100,2),LC(10,2),XALU(100),LCL(10),
     1 LVL(10),TC1(100),TC2(100),TEMPE(10),NTMP,KTMP,NOIS,NPOS,NNEG
      COMMON /CARD/ KARD(80),IPNT(80),IPT,NERR,KFLG,NAM
      DOUBLE PRECISION TITLE
      COMMON /JOBCNT/ TITLE(9),PRGVER,NALTR,IREQ,ISTOP
      DATA PRGVER /3HD4 /                                               10/13/72
      DIMENSION ADATE(2),ATIME(2)
C
      NALTR=0
      KFLG=0
C
C.... GET PROGRAM STARTING TIME
 100  IF (NALTR.NE.0) GO TO 120
      IDED=1
      CALL CLKST(IDED)
C
C.... READ AND PRINT TITLE
      IF (KFLG.EQ.0) GO TO 105
      KFLG=0
      K=-7
      DO 102 I=1,9
      K=K+8
      L=K+7
      ENCODE (8,901,TITLE(I)) (KARD(J),J=K,L)
102   CONTINUE
      GO TO 110
105   READ (1,900,END=800) TITLE
110   CALL DATE(ADATE)
      CALL TIME(ATIME(1),ATIME(2))
      WRITE (6,905) PRGVER,ADATE,ATIME,TITLE
C
C     ... READ CIRCUIT DATA
120   ISTOP=0
      CALL READIN
      IF (ISTOP.EQ.0) GO TO 200
      WRITE (6,1130)
      GO TO 999
C
C     ... INITIALIZE CIRCUIT
 200  IF (NALTR.GT.3) GO TO 999
      IDED=2
      CALL CLKST(IDED)
      CALL INITL
      IF (ISTOP.NE.0) GO TO 999
C
C     ... PERFORM CIRCUIT ANALYSIS
250   DO 260 KTMP=1,NTMP
      IDED=2
      CALL CLKST(IDED)
      IF (TEMPE(KTMP).LT.100.) GO TO 260
C
C     ... CALCULATE DC SOLUTION
      CALL DCSOL
      IF (ISTOP.NE.0) GO TO 260
      IF (TMAX.EQ.0.) GO TO 260
      IDED=2
      CALL CLKST(IDED)
C
C     ... PERFORMS TRANSIENT ANALYSIS
      CALL CONTL
C
C     ... PRINT PLOTS IF REQUESTED
      IF (NPLT.EQ.0) GO TO 260
      IDED=2
      CALL CLKST(IDED)
      CALL PLOT
260   CONTINUE
C
C.... CALCULATE AND PRINT RUN TIME
999   IDED=2
      CALL CLKST(IDED)
      GO TO 100
 800  CALL EXIT
C
900   FORMAT (9A8)
901   FORMAT (8A1)
 905  FORMAT (1H1,11H- S I N C -,2X,5HVER. ,A3,2XA5,A4,2XA5,A4,2X,
     1 45(1HX),///1X9A8///1X,90(1HX)//)
 1130 FORMAT (1X,28H****DATA ERROR - JOB ABORTED)
      END
C
      SUBROUTINE CONTL
C
C     ... CONTROLS TRANSIENT ANALYSES
C
      COMMON Y(600),YL(600),YV(200),YD(250),V(250),V1(250),BCUR(100),
     1 CBC(100),BCIP(100),CBE(100),BEIP(100),VBC(100),VBE(100),GBE(100),
     2 GBC(100),CJEP(100),CJCP(100),PLTP(2,200,6),OUPUT(10),VSLP(10),
     3 CSLP(10),VPI(250)
      COMMON /TRANS/ MREF(100),NQ(100,17),VT,TOLER,NXTR,ITR,IND
      COMMON /NODES/ NREN(100),LNOD,NX,MDIM,NNOD,NUT,NLT,NVT
      COMMON /ELMTS/ KIND(100),ND(100,4),VALU(100),NELT
      COMMON /VSRCE/ VSRC(10),NAMV(10),TVLT(10,10),VVAL(10,10),
     1 NVVAL(10),KVOL(10),NV(10,2),NVSC
      COMMON /CSRCE/ CSRC(10),NAMC(10),TCUR(10,10),CVAL(10,10),
     1 NCVAL(10),KCUR(10),NCS(10,2),NCSC
      COMMON /PLOTS/ PMIN(6),PMAX(6),KPLO(2,6),NPVAL(6),NAMP(2,6),
     1 NP(2,2,6),NPLT,NPTS,FSTR,FSTP,NHAR
      COMMON /OUTPT/ KOUT(10),NAMO(10),NODO(10,2),NOUT
      COMMON /TIMEV/ TTIME,DELT,TIMP,DELP,TMAX,TPRT,DPRT,ITER,NCUT,ICONV
      DOUBLE PRECISION TITLE
      COMMON /JOBCNT/ TITLE(9),PRGVER,KBREF,IREQ,ISTOP
C
C     ... START TRANSIENT ANALYSIS
      IF (NOUT.EQ.0) GO TO 130
      WRITE (6,900) TITLE,PRGVER
      WRITE (6,910) (NAMO(I),I=1,NOUT)
      WRITE (6,920) ((NODO(I,J),J=1,2),I=1,NOUT)
      WRITE (6,930)
C
C.... INITIALIZE TIME VALUES
130   TTIME=0.
      TPRT=0.
      DELIM=DPRT*1.E-12
      DELT=DPRT
      NCUT=1
      NPTS=0
      ITER=0
      JTER=0
      NSOL=0
      DO 132 I=1,NVSC
      VSLP(I)=0.0
132   NVVAL(I)=0
      DO 135 I=1,NCSC
      CSLP(I)=0.0
135   NCVAL(I)=0
C
C.... PRINT OUTPUTS IF TIME
140   NOTXAT=0
      NSOL=NSOL+1
      IF (TTIME.LT.TPRT-1.E-15) GO TO 175
      IF (NOUT.EQ.0) GO TO 164
      DO 160 I=1,NOUT
      IF (KOUT(I)) 155,145,150
C.... VOLTAGE OUTPUT
145   N1=NODO(I,1)+1
      N2=NODO(I,2)+1
      N1=NREN(N1)
      N2=NREN(N2)
      OUPUT(I)=V(N1)-V(N2)
      GO TO 160
C.... CURRENT OUTPUT ACROSS ELEMENT
150   N=KOUT(I)
      OUPUT(I)=BCUR(N)*1.E3
      GO TO 160
C.... CURRENT OUTPUT ACROSS C-SOURCE
155   N=-KOUT(I)
      OUPUT(I)=CSRC(N)*1.E3
160   CONTINUE
      WRITE (6,930) TTIME,(OUPUT(I),I=1,NOUT)
C
C.... STORE PLOT DATA
164   IF (NPLT.EQ.0) GO TO 170
      NPTS=NPTS+1
      IF (NPTS.GT.200) GO TO 170
      DO 169 K=1,NPLT
      N=NPVAL(K)
      DO 168 I=1,N
      IF (KPLO(I,K)) 167,165,166
C.... VOLTAGE OUTPUT PLOT
165   N1=NP(1,I,K)+1
      N2=NP(2,I,K)+1
      N1=NREN(N1)
      N2=NREN(N2)
      PLTP(I,NPTS,K)=V(N1)-V(N2)
      GO TO 168
C.... CURRENT OUTPUT PLOT
166   KP=KPLO(I,K)
      PLTP(I,NPTS,K)=BCUR(KP)*1.E3
      GO TO 168
C.... CURRENT OUT ACROSS C-SOURCE
167   KP=-KPLO(I,K)
      PLTP(I,NPTS,K)=CSRC(KP)*1.E3
168   CONTINUE
169   CONTINUE
C
C.... ADVANCE TIME
170   TPRT=TPRT+DPRT
      IF (TPRT.GT.TMAX) GO TO 550
C     ... INCREMENT TIME VALUE
 175  IF (NCUT.NE.0) GO TO 185
      IF (DELT.LT.DELA) GO TO 180
      IF (DELP.LE.DELA) GO TO 180
      DELT=DELP
      GO TO 185
180   IF (ITR.GT.3) GO TO 185
      DELT=DELT+DELT
185   TIMP=TTIME
      TTIME=TIMP+DELT
      DELP=DELT
C     ... CUT TIME TO MATCH PRINT INTERVAL IF NECESSARY
      IF (TTIME.LE.TPRT) GO TO 190
      TTIME=TPRT
      DELT=TTIME-TIMP
 190  DELA=DELT
      ITRLIM=8
      NCUT=0
      ITR=0
      DO 195 I=1,NNOD
195   VPI(I)=V(I)
C
C     ... ADJUST TIME FOR BREAK POINTS
 200  IF (NVSC.EQ.0) GO TO 280
      DO 250 I=1,NVSC
      IF (KVOL(I).NE.1) GO TO 250
      N=NVVAL(I)
      IF (TTIME-TVLT(N+1,I)) 250,250,225
 225  N=N+1
      DTEM=TVLT(N+1,I)-TVLT(N,I)
      SLP=VSLP(I)
      VSLP(I)=(VVAL(N+1,I)-VVAL(N,I))/DTEM
      IF (ABS(VSLP(I))-ABS(SLP)) 245,245,230
 230  DTEM=DTEM*0.1
      IF (TVLT(N,I)-TIMP-DTEM) 240,240,235
 235  TTIME=TVLT(N,I)
      DELT=TTIME-TIMP
      VSLP(I)=SLP
      GO TO 250
 240  DELT=AMIN1(DELT,DTEM)
      TTIME=TIMP+DELT
 245  NVVAL(I)=N
 250  CONTINUE
C
 280  IF (NCSC.EQ.0) GO TO 400
      DO 330 I=1,NCSC
      IF (KCUR(I).NE.1) GO TO 330
      N=NCVAL(I)
      IF (TTIME-TCUR(N+1,I)) 330,330,305
 305  N=N+1
      DTEM=TCUR(N+1,I)-TCUR(N,I)
      SLP=CSLP(I)
      CSLP(I)=(CVAL(N+1,I)-CVAL(N,I))/DTEM
      IF (ABS(CSLP(I))-ABS(SLP)) 325,325,310
 310  DTEM=DTEM*0.1
      IF (TCUR(N,I)-TIMP-DTEM) 320,320,315
 315  TTIME=TCUR(N,I)
      DELT=TTIME-TIMP
      CSLP(I)=SLP
      GO TO 330
 320  DELT=AMIN1(DELT,DTEM)
      TTIME=TIMP+DELT
 325  NCVAL(I)=N
 330  CONTINUE
C
C     ... CALCULATE NEW VOLTAGE VALUES
400   CALL TIMSOL (ITRLIM)
      IF (ICONV.EQ.0) GO TO 140
C
C.... NO CONVERGENCE - ADJUST TIME AND RETRY
      JTER=JTER+ITR
      DELT=DELT*.125
      IF (DELT.GT.DELIM) GO TO 480
      WRITE (6,940)
      GO TO 550
480   TTIME=TIMP+DELT
      NCUT=NCUT+1
      IF (NCUT.GE.7) GO TO 520
      ITR=0
      DO 490 I=1,NNOD
490   V(I)=VPI(I)
      GO TO 200
C.... TOO MANY CUTS.  PRINT WARNING AND CONTINUE
 520  IF (NOUT.NE.0) WRITE (6,950)
      NOTXAT=NOTXAT+1
      IF (NOTXAT.LT.2) GO TO 140
      WRITE (6,970)
C
550   WRITE (6,980) ITER,JTER,NSOL
      RETURN
C
 900  FORMAT (1H1,9A8,8X,5HSINC-,A3/)
910   FORMAT (/5X,4HTIME,1X,9(9X,A4))
920   FORMAT (10X,9(7X,2I3))
930   FORMAT (1X,1PE11.3,1X,0P9F13.6)
940   FORMAT (//1X,24HSOLUTION TOO SLOW - STOP/)
950   FORMAT (1X,20H**SOLUTION NOT EXACT)
 970  FORMAT (//1X,29H****UNSTABLE SOLUTION - ABORT/)
980   FORMAT (//1X,20HTOTAL NEWTON ITER. =,I5,2H -,I4,3X,11HNO. SOLTN =,
     1 I4/)
      END
C
      SUBROUTINE TIMSOL (ITRLIM)
C
C     .... SUBROUTINE TO CALCULATE VOLTAGES AT A GIVEN TIME
C
      COMMON Y(600),YL(600),YV(200),YD(250),V(250),V1(250),BCUR(100),
     1 CBC(100),BCIP(100),CBE(100),BEIP(100),VBC(100),VBE(100),GBE(100),
     2 GBC(100),CJEP(100),CJCP(100),PLTP(2,200,6),OUPUT(10),VSLP(10),
     3 CSLP(10),VPI(250)
      COMMON /TRANS/ MREF(100),NQ(100,17),VT,TOLER,NXTR,ITR,IND
      COMMON /MODEL/ NTYP(30),BETA(30,8),CSAT(30,1),RBS(30,3),RCL(30,3),
     1 ROU(30,2),FT(30,4),CJ1(30,4),CJ2(30,6),TEMP(30)
      COMMON /ELMTS/ KIND(100),ND(100,4),VALU(100),NELT
      COMMON /VSRCE/ VSRC(10),NAMV(10),TVLT(10,10),VVAL(10,10),
     1 NVVAL(10),KVOL(10),NV(10,2),NVSC
      COMMON /CSRCE/ CSRC(10),NAMC(10),TCUR(10,10),CVAL(10,10),
     1 NCVAL(10),KCUR(10),NCS(10,2),NCSC
      COMMON /NODES/ NREN(100),LNOD,NX,MDIM,NNOD,NUT,NLT,NVT
      COMMON /TIMEV/ TTIME,DELT,TIMP,DELP,TMAX,TPRT,DPRT,ITER,NCUT,ICONV
C
C     ... ASSIGN SOURCE VALUES FOR NEW TIME
      IF (NVSC.EQ.0) GO TO 370
      DO 360 I=1,NVSC
      IF (KVOL(I)-1) 360,350,345
 345  VSRC(I)=VVAL(1,I)*SIN(VVAL(2,I)*TTIME+VVAL(3,I))+VVAL(4,I)
      GO TO 360
 350  N=NVVAL(I)
      IF (VSLP(I).EQ.0.0) GO TO 355
      VSRC(I)=VVAL(N,I)+VSLP(I)*(TTIME-TVLT(N,I))
      GO TO 360
 355  VSRC(I)=VVAL(N,I)
 360  CONTINUE
C
 370  IF (NCSC.EQ.0) GO TO 400
      DO 390 I=1,NCSC
      IF (KCUR(I)-1) 390,380,375
 375  CSRC(I)=CVAL(1,I)*SIN(CVAL(2,I)*TTIME+CVAL(3,I))
      GO TO 390
 380  N=NCVAL(I)
      IF (CSLP(I).EQ.0.0) GO TO 385
      CSRC(I)=CVAL(N,I)+CSLP(I)*(TTIME-TCUR(N,I))
      GO TO 390
 385  CSRC(I)=CVAL(N,I)
 390  CONTINUE
C
C     ... CLEAR Y-MATRIX AND C-VECTOR
400   DO 105 I=1,NUT
105   Y(I)=0.0
      DO 110 I=1,NLT
110   YL(I)=0.0
      DO 115 I=1,NVT
115   YV(I)=0.0
      DO 120 I=1,NNOD
      V1(I)=V(I)
      V(I)=0.0
120   YD(I)=0.0
C
C     .... ADD R,C,L ELEMENT CONTRIBUTION
      IF (NELT.EQ.0) GO TO 155
      DO 150 N=1,NELT
      N1=ND(N,1)
      N2=ND(N,2)
      N3=ND(N,3)
      N4=ND(N,4)
      DELVP=VPI(N1)-VPI(N2)
      IF (KIND(N)-2) 125,130,135
C     .... RESISTOR
125   G=VALU(N)
      GO TO 145
C     ... CAPACITOR
130   G=VALU(N)/DELT
      CT=G*DELVP+BCUR(N)
      GO TO 140
C     ... INDUCTOR
135   G=VALU(N)*DELT
      CT=-G*DELVP-BCUR(N)
140   V(N1)=V(N1)+CT
      V(N2)=V(N2)-CT
145   YD(N1)=YD(N1)+G
      YD(N2)=YD(N2)+G
      Y(N3)=Y(N3)-G
      Y(N4)=Y(N4)-G
150   CONTINUE
C
C     ... ADD TRANSISTORS AND CURRENT SOURCES
155   CALL CKTMOD (1)
C
C     .... SOLVE EQUATIONS FOR NEW VOLTAGES
      ITR=ITR+1
      ITER=ITER+1
      CALL SOLUTN
C
C     .... CHECK FOR CONVERGENCE OF NON-LINEAR MODELS
      ICONV=0
      IF (NXTR.EQ.0) GO TO 170
      TOL=VT*TOLER
      DO 165 N=1,NXTR
      NC=NQ(N,10)
      NB=NQ(N,11)
      NE=NQ(N,3)
      PE=V(NB)-V(NE)
      PC=V(NB)-V(NC)
      K=MREF(N)
      IF (NTYP(K).NE.0) GO TO 160
      PE=-PE
      PC=-PC
160   IF (ABS(PE-VBE(N)).GT.TOL) GO TO 470
      IF (ABS(PC-VBC(N)).GT.TOL) GO TO 470
165   CONTINUE
      ITR=0
C
C     ... CALCULATE BRANCH CURRENTS
170   IF (NELT.EQ.0) GO TO 465
      DO 460 N=1,NELT
      N1=ND(N,1)
      N2=ND(N,2)
      DELV=V(N1)-V(N2)
      DELVP=VPI(N1)-VPI(N2)
      IF (KIND(N)-2) 440,445,450
C.... RESISTOR
440   BCUR(N)=DELV*VALU(N)
      GO TO 460
C     ... CAPACITOR
 445  G=VALU(N)/DELT
      BCUR(N)=G*(DELV-DELVP)-BCUR(N)
      GO TO 460
C     ... INDUCTOR
 450  G=VALU(N)*DELT
      BCUR(N)=G*(DELV+DELVP)+BCUR(N)
460   CONTINUE
465   RETURN
C
470   ICONV=1
      IF (ITR-ITRLIM) 400,500,500
500   RETURN
      END
C
      SUBROUTINE DCSOL
C
C     .... SUBROUTINE TO CALCULATE THE DC OPERATING POINT
C
      COMMON Y(600),YL(600),YV(200),YD(250),V(250),V1(250),BCUR(100),
     1 CBC(100),BCIP(100),CBE(100),BEIP(100),VBC(100),VBE(100),GBE(100),
     2 GBC(100),CJEP(100),CJCP(100),PLTP(2,200,6),OUPUT(10),VSLP(10),
     3 CSLP(10),VPI(250)
      COMMON /NODES/ NREN(100),LNOD,NX,MDIM,NNOD,NUT,NLT,NVT
      COMMON /TRANS/ MREF(100),NQ(100,17),VT,TOLER,NXTR,ITR,IND
      COMMON /MODEL/ NTYP(30),BETA(30,8),CSAT(30,1),RBS(30,3),RCL(30,3),
     1 ROU(30,2),FT(30,4),CJ1(30,4),CJ2(30,6),TEMP(30)
      COMMON /MODPR/ BEF(30),BER(30),CSA(30),TCF(30),TCR(30),GBS(30),
     1 GCL(30),BWF(30),BC1(30),BC2(30),BC3(30),CJE(30),CJC(30),CT1(30),
     2 CT2(30)
      COMMON /ELMTS/ KIND(100),ND(100,4),VALU(100),NELT
      COMMON /VSRCE/ VSRC(10),NAMV(10),TVLT(10,10),VVAL(10,10),
     1 NVVAL(10),KVOL(10),NV(10,2),NVSC
      COMMON /CSRCE/ CSRC(10),NAMC(10),TCUR(10,10),CVAL(10,10),
     1 NCVAL(10),KCUR(10),NCS(10,2),NCSC
      COMMON /TIMEV/ TTIME,DELT,TIMP,DELP,TMAX,TPRT,DPRT,ITER,NCUT,ICONV
      COMMON /NAMES/ NAME(100),NAMQ(100),MODN(100),NAMB(30),NMOD
      COMMON /STORE/ LQ(100,3),LD(100,2),LC(10,2),XALU(100),LCL(10),
     1 LVL(10),TC1(100),TC2(100),TEMPE(10),NTMP,KTMP,NOIS,NPOS,NNEG
      DOUBLE PRECISION TITLE
      COMMON /JOBCNT/ TITLE(9),PRGVER,KBREF,IREQ,ISTOP
      DIMENSION NDM(4),VDM(4)
      DATA LPAR /1H(/
C
      IF (NELT.EQ.0) GO TO 125
      DO 120 I=1,NELT
      DT=TEMPE(KTMP)-TEMPE(1)
      TEMF=1.+DT*(TC1(I)+DT*TC2(I))
      IF (KIND(I)-2) 105,110,115
105   VALU(I)=1./(XALU(I)*TEMF)
      GO TO 120
110   VALU(I)=(XALU(I)+XALU(I))*TEMF
      GO TO 120
115   VALU(I)=.5/(XALU(I)*TEMF)
120   CONTINUE
125   DO 130 I=1,NNOD
130   V(I)=0.0
      IF (NVSC.EQ.0) GO TO 145
      DO 140 I=1,NVSC
      IF (KVOL(I).EQ.2) GO TO 135
      VSRC(I)=VVAL(1,I)
      GO TO 140
135   VSRC(I)=VVAL(1,I)*SIN(VVAL(3,I))+VVAL(4,I)
140   CONTINUE
145   IF (NCSC.EQ.0) GO TO 160
      DO 155 I=1,NCSC
      IF (KCUR(I).EQ.2) GO TO 150
      CSRC(I)=CVAL(1,I)
      GO TO 155
150   CSRC(I)=CVAL(1,I)*SIN(CVAL(3,I))
155   CONTINUE
160   IF (NXTR.EQ.0) GO TO 170
      VT=TEMPE(KTMP)*8.6164E-5
      IF (KTMP.GT.1) GO TO 170
      DO 165 I=1,NXTR
      VBE(I)=.6
165   VBC(I)=0.
C
C     .... ITERATE TO SOLUTION
170   DO 225 ITER=1,100
      ITR=ITER-1
      IND=ITR
C
C     ... CLEAR MATRIX AND UPDATE VOLTAGES
      DO 175 I=1,NUT
175   Y(I)=0.0
      DO 180 I=1,NLT
180   YL(I)=0.0
      DO 185 I=1,NVT
185   YV(I)=0.0
      DO 190 I=1,NNOD
      V1(I)=V(I)
      V(I)=0.0
190   YD(I)=0.0
C
C     .... ADD RESISTORS AND INDUCTORS TO MATRIX
      IF (NELT.EQ.0) GO TO 210
      DO 205 I=1,NELT
      IF (KIND(I).EQ.2) GO TO 205
      N1=ND(I,1)
      N2=ND(I,2)
      N3=ND(I,3)
      N4=ND(I,4)
      IF (KIND(I).EQ.1) GO TO 195
C     ... INDUCTOR - 1 OHM RESISTANCE
      G=1.0
      GO TO 200
C     ... RESISTOR
195   G=VALU(I)
200   YD(N1)=YD(N1)+G
      YD(N2)=YD(N2)+G
      Y(N3)=Y(N3)-G
      Y(N4)=Y(N4)-G
205   CONTINUE
C
C     ... ADD TRANSISTORS AND CURRENT SOURCES
210   CALL CKTMOD (0)
C
C     .... SOLVE FOR NEW VOLTAGES
      CALL SOLUTN
C
C     .... CHECK FOR CONVERGENCE
      IF (ITR.EQ.0) GO TO 220
      DO 215 I=1,MDIM
      IF (ABS(V(I)-V1(I)).GT.1.E-4) GO TO 220
215   CONTINUE
      IF (NXTR.EQ.0) GO TO 425
      DO 410 N=1,NXTR
      K=MREF(N)
      NC=NQ(N,10)
      NB=NQ(N,11)
      NE=NQ(N,3)
      PE=V(NB)-V(NE)
      PC=V(NB)-V(NC)
      IF (NTYP(K).NE.0) GO TO 405
      PE=-PE
      PC=-PC
405   IF (ABS(PE-VBE(N)).GT.1.E-4) GO TO 220
      IF (ABS(PC-VBC(N)).GT.1.E-4) GO TO 220
410   CONTINUE
425   ICONV=ICONV+1
      IF (ICONV-1) 225,225,235
220   ICONV=0
225   CONTINUE
C
C     .... ITERATION LIMIT EXCEEDED - PRINT FINAL VALUES
      WRITE (6,260)
      DO 232 N=2,LNOD
      N1=N-1
      I=NREN(N)
      IF (I.EQ.0) GO TO 232
      PCTERR=0.
      IF (V(I).EQ.0.) GO TO 230
      PCTERR=ABS((V1(I)-V(I))/V(I))*100.
230   WRITE (6,265) N1,V(I),PCTERR
232   CONTINUE
      ISTOP=1
      RETURN
C
C     .... CONVERGENCE OBTAINED - CALCULATE BRANCH CURRENTS
235   ITR=0
      IF (NELT.EQ.0) GO TO 255
      DO 250 I=1,NELT
      IF (KIND(I).EQ.2) GO TO 245
      N1=ND(I,1)
      N2=ND(I,2)
      IF (KIND(I).EQ.1) GO TO 240
C     ... OBTAIN INDUCTOR CURRENT
      BCUR(I)=V(N1)-V(N2)
      GO TO 250
C     .... CALCULATE RESISTOR CURRENT
240   BCUR(I)=(V(N1)-V(N2))*VALU(I)
      GO TO 250
C     .... ZERO CAPACITOR CURRENT
245   BCUR(I)=0.
250   CONTINUE
C
C     ... PRINT NODE VOLTAGES
255   N=0
      WRITE (6,270) TITLE,PRGVER
      WRITE (6,275) TEMPE(KTMP),ITER
300   DO 310 K=1,4
305   N=N+1
      IF (N.GT.LNOD) GO TO 315
      I=NREN(N)
      IF (I.EQ.0) GO TO 305
      NDM(K)=N-1
      VDM(K)=V(I)
310   CONTINUE
      K=5
315   K1=K-1
      IF (K1.EQ.0) GO TO 350
      WRITE (6,280) (LPAR,NDM(K),VDM(K),K=1,K1)
      IF (N.LT.LNOD) GO TO 300
C
C     ... PRINT TRANSISTOR AND RESISTOR CURRENTS
350   IF (NXTR.EQ.0) GO TO 370
      WRITE (6,940)
      DO 365 N=1,NXTR
      K=MREF(N)
      CF=AMAX1(CBE(N),1.E-9)
      BF=BEF(K)/(BC1(K)+BC2(K)/SQRT(CF)+BC2(K)*CF)
      CEB=CBE(N)/BF
      CCB=CBC(N)/BER(K)
      VCE=VBE(N)-VBC(N)
      FAC=1.-BWF(K)*VBC(N)
      CC=((CBE(N)-CBC(N))*FAC-CCB)*1.E3
      CB=(CEB+CCB)*1.E3
      VBET=VBE(N)
      NTP=1H 
      IF (NTYP(K).NE.0) GO TO 360
      CC=-CC
      CB=-CB
      VBET=-VBET
      VCE=-VCE
      NTP=1HP
360   WRITE (6,945) NAMQ(N),NTP,CC,CB,VBET,VCE,BF
365   CONTINUE
370   IF (NELT.EQ.0) GO TO 390
      WRITE (6,950)
      DO 380 N=1,NELT
      IF (KIND(N).EQ.2) GO TO 380
      N1=ND(N,1)
      N2=ND(N,2)
      CURR=BCUR(N)*1.E3
      RPW=CURR*(V(N1)-V(N2))
      WRITE (6,955) NAME(N),CURR,RPW
380   CONTINUE
390   WRITE (6,955)
C
450   RETURN
C
260   FORMAT (/1X,35HNO CONVERGENCE AFTER 100 ITERATIONS/1X,35(1H*)//
     1 1X,4HNODE,4X,7HVOLTAGE,4X,'PER CENT ERROR'/)
265   FORMAT (1X,I3,F12.5,F11.3)
270   FORMAT (1H1,9A8,8X,5HSINC-,A3)
275   FORMAT (/1X,6HTEMP =,F9.2//1X,14HNO. OF ITER. =,I4//1X,
     1 16HNODE VOLTAGES --)
280   FORMAT (2X,4(A1,I2,1H),F11.5,4X))
940   FORMAT (//1X,21HTRANSISTOR OP. PT. --/1X,7HNAME TP,5X,7HIC (MA),
     1 4X,7HIB (MA),6X,3HVBE,8X,3HVCE,6X,2HBF)
945   FORMAT (1X,A5,A2,4F11.5,F9.2)
950   FORMAT (//1X,20HRESISTOR CURRENTS --/1X,4HNAME,4X,12HCURRENT (MA),
     1 3X,10HPOWER (MW))
955   FORMAT (1X,A5,2X,F12.5,F13.4)
      END
C
      SUBROUTINE CKTMOD (MODE)
C
C     .... SET UP EQUIVALENT CIRCUITS
C
      COMMON Y(600),YL(600),YV(200),YD(250),V(250),V1(250),BCUR(100),
     1 CBC(100),BCIP(100),CBE(100),BEIP(100),VBC(100),VBE(100),GBE(100),
     2 GBC(100),CJEP(100),CJCP(100),PLTP(2,200,6),OUPUT(10),VSLP(10),
     3 CSLP(10),VPI(250)
      COMMON /TRANS/ MREF(100),NQ(100,17),VT,TOLER,NXTR,ITR,IND
      COMMON /MODEL/ NTYP(30),BETA(30,8),CSAT(30,1),RBS(30,3),RCL(30,3),
     1 ROU(30,2),FT(30,4),CJ1(30,4),CJ2(30,6),TEMP(30)
      COMMON /CSRCE/ CSRC(10),NAMC(10),TCUR(10,10),CVAL(10,10),
     1 NCVAL(10),KCUR(10),NCS(10,2),NCSC
      COMMON /TIMEV/ TTIME,DELT,TIMP,DELP,TMAX,TPRT,DPRT,ITER,NCUT,ICONV
      COMMON /NAMES/ NAME(100),NAMQ(100),MODN(100),NAMB(30),NMOD
      COMMON /STORE/ LQ(100,3),LD(100,2),LC(10,2),XALU(100),LCL(10),
     1 LVL(10),TC1(100),TC2(100),TEMPE(10),NTMP,KTMP,NOIS,NPOS,NNEG
      COMMON /MODPR/ BEF(30),BER(30),CSA(30),TCF(30),TCR(30),GBS(30),
     1 GCL(30),BWF(30),BC1(30),BC2(30),BC3(30),CJE(30),CJC(30),CT1(30),
     2 CT2(30)
C
      IF (MODE.EQ.0) GO TO 100
      DELTI=1.0/DELT
100   IF (NXTR.EQ.0) GO TO 450
      IF (IND.NE.0) GO TO 120
C
C     ... INITIALISE MODEL PARAMETERS
      DO 440 K=1,NMOD
      IF (KTEMP.GT.1) GO TO 410
C
C     ... TEMP. DEPEND. PARAMETERS
410   CSA(K)=CSAT(K,1)*(TEMPE(KTMP)/TEMP(K))**3
      CSA(K)=CSA(K)*EXP(13920.*(1./TEMP(K)-1./TEMPE(KTMP)))
      DT=TEMPE(KTMP)-TEMP(K)
      TEMF=1.+DT*(BETA(K,6)+DT*BETA(K,7))
      BEF(K)=TEMF
      BER(K)=TEMF*BETA(K,8)
      IF (RBS(K,1).EQ.0.) GO TO 420
      GBS(K)=1./(RBS(K,1)*(1.+DT*(RBS(K,2)+DT*RBS(K,3))))
420   IF (RCL(K,1).EQ.0.) GO TO 440
      GCL(K)=1./(RCL(K,1)*(1.+DT*(RCL(K,2)+DT*RCL(K,3))))
440   CONTINUE
      VT10=VT*10.
      VT2=VT*2.3
C
C     ... PROCESS TRANSISTORS
120   DO 340 N=1,NXTR
      K=MREF(N)
      NT=NTYP(K)
      NC=NQ(N,10)
      NB=NQ(N,11)
      N1=NQ(N,1)
      N2=NQ(N,2)
      NE=NQ(N,3)
      N4=NQ(N,4)
      N5=NQ(N,5)
      N6=NQ(N,6)
      N7=NQ(N,7)
      N8=NQ(N,8)
      N9=NQ(N,9)
C
C     ... DC TRANSISTOR MODEL
      IF (IND.NE.0) GO TO 125
      PE=VBE(N)
      PC=VBC(N)
      GO TO 130
125   PE=V1(NB)-V1(NE)
      PC=V1(NB)-V1(NC)
      IF (NT.NE.0) GO TO 130
      PE=-PE
      PC=-PC
C
C     ... DETERMINE JUNCTION VOLTAGES
130   CS=CSA(K)
      IF (ITR.EQ.0) GO TO 150
      IF (VBE(N).LT.VT10) GO TO 137
      BVT=PE-VBE(N)
      IF (ABS(BVT).LE.VT2) GO TO 140
      IF (BVT) 132,132,135
132   PE=VBE(N)-VT2
      GO TO 140
135   PE=VBE(N)+VT2
      GO TO 140
137   IF (PE.LE.VT10) GO TO 140
      PE=VT10
140   IF (VBC(N).LT.VT10) GO TO 147
      BVT=PC-VBC(N)
      IF (ABS(BVT).LE.VT2) GO TO 150
      IF (BVT) 142,142,145
142   PC=VBC(N)-VT2
      GO TO 150
145   PC=VBC(N)+VT2
      GO TO 150
147   IF (PC.LE.VT10) GO TO 150
      PC=VT10
150   VBE(N)=PE
      VBC(N)=PC
C
C     ... CALCULATE CURRENTS AND CONDUCTANCES
      IF (PE.LE.0.0) GO TO 160
      CF=CS*(EXP(PE/VT)-1.0)
155   GMF=(CF+CS)/VT
      CE=CF-GMF*PE
      GO TO 170
160   GMF=CS/VT
      CF=GMF*PE
      CE=0.0
170   IF (PC.LE.0.0) GO TO 180
      CR=CS*(EXP(PC/VT)-1.0)
175   GMR=(CR+CS)/VT
      CC=CR-GMR*PC
      GO TO 190
180   GMR=CS/VT
      CR=GMR*PC
      CC=0.0
C
C     ... DC BASE CIRCUIT
190   DF=AMAX1(CF,1.E-9)
      BF=BEF(K)/(BC1(K)+BC2(K)/SQRT(DF)+BC3(K)*DF)
      BR=BER(K)
      CEB=CE/BF
      CCB=CC/BR
      GPIF=GMF/BF+1.E-9
      GPIR=GMR/BR+1.E-9
      CBE(N)=CF
      CBC(N)=CR
      IF (MODE.EQ.0) GO TO 250
C
C     ... TRANSIENT TRANSISTOR MODEL
      IF (ITR.NE.0) GO TO 230
C     ... SET PAST CURRENTS
      BEIP(N)=CBE(N)
      BCIP(N)=CBC(N)
C     ... JUNCTION CAPACITANCES
      CJEP(N)=0.
      IF (CJE(K).EQ.0) GO TO 215
      VQT=PE/CJ1(K,3)
      IF (VQT.GE..5) GO TO 210
      CJEP(N)=CJE(K)/EXP(CJ1(K,4)*ALOG(1.-VQT))
      GO TO 215
210   CJEP(N)=CJE(K)*(1.+VQT)/CT1(K)
215   CJCP(N)=0.
      IF (CJC(K).EQ.0) GO TO 230
      VQT=PC/CJ2(K,3)
      IF (VQT.GE..5) GO TO 220
      CJCP(N)=CJC(K)/EXP(CJ2(K,4)*ALOG(1.-VQT))
      GO TO 230
220   CJCP(N)=CJC(K)*(1.+VQT)/CT2(K)
C
C     ... TRAN. BASE AND COLLECTOR CIRCUIT
230   TCE=TCF(K)*DELTI
      TCC=TCR(K)*DELTI
      BEC=CJEP(N)*DELTI
      BCC=CJCP(N)*DELTI
C     ... SPLIT CAPACITANCE
      IF (CJ2(K,5).EQ.0.) GO TO 235
      IF (N2.EQ.NB) GO TO 235
      GCST=BCC*CJ2(K,5)
      BCC=BCC-GCST
      CTB=GCST*(VPI(N2)-VPI(NC))
      N16=NQ(N,16)
      N17=NQ(N,17)
      YD(N2)=YD(N2)+GCST
      YD(NC)=YD(NC)+GCST
      Y(N16)=Y(N16)-GCST
      Y(N17)=Y(N17)-GCST
      V(N2)=V(N2)+CTB
      V(NC)=V(NC)-CTB
C     ... SUBSTRATE CAPACITANCE
235   IF (CJ2(K,6).EQ.0.) GO TO 237
      GCST=CJ2(K,6)*DELTI
      CTC=GCST*VPI(NC)
      YD(NC)=YD(NC)+GCST
      V(NC)=V(NC)+CTC
C
C     ... DIFFUSION AND JUNCTION CAPACITANCES
237   GPIF=GPIF+GMF*TCE+BEC
      GPIR=GPIR+GMR*TCC+BCC
      BEC=BEC*(VPI(NB)-VPI(NE))
      BCC=BCC*(VPI(NB)-VPI(NC))
      IF (NT.NE.0) GO TO 240
      BEC=-BEC
      BCC=-BCC
240   CEB=CEB+(CE-BEIP(N))*TCE-BEC
      CCB=CCB+(CC-BCIP(N))*TCC-BCC
C
C     ... SET CM-EMT 2-PORT PARAMETERS
250   CTB=CEB+CCB
C
C     ... BASE WIDTH MODULATION
      VAI=BWF(K)
      IF (VAI.EQ.0.) GO TO 255
      FAC=1.-VAI*PC
      GMF=GMF*FAC
      GMR=GMR*FAC+(CF-CR)*VAI
      CTC=(CF-CR)*FAC-GMF*PE+GMR*PC-CCB
      GO TO 257
255   CTC=CE-CC-CCB
257   IF (NT.NE.0) GO TO 260
      CTB=-CTB
      CTC=-CTC
C
C     ... INCORPORATE SERIES RESISTANCES
260   IF (N1.EQ.NC) GO TO 270
      GC=GCL(K)
      N12=NQ(N,12)
      N13=NQ(N,13)
      YD(N1)=YD(N1)+GC
      YD(NC)=YD(NC)+GC
      Y(N12)=Y(N12)-GC
      Y(N13)=Y(N13)-GC
270   IF (N2.EQ.NB) GO TO 300
      GB=GBS(K)
      N14=NQ(N,14)
      N15=NQ(N,15)
      YD(N2)=YD(N2)+GB
      YD(NB)=YD(NB)+GB
      Y(N14)=Y(N14)-GB
      Y(N15)=Y(N15)-GB
C
C     ... ADD TO MATRIX
300   YD(NC)=YD(NC)+(GPIR+GMR)
      YD(NB)=YD(NB)+(GPIF+GPIR)
      YD(NE)=YD(NE)+(GPIF+GMF)
      Y(N4)=Y(N4)-(GPIR+GMR-GMF)
      Y(N5)=Y(N5)-GMF
      Y(N6)=Y(N6)-GPIR
      Y(N7)=Y(N7)-GPIF
      Y(N8)=Y(N8)-GMR
      Y(N9)=Y(N9)-(GPIF+GMF-GMR)
      V(NC)=V(NC)-CTC
      V(NB)=V(NB)-CTB
      V(NE)=V(NE)+(CTC+CTB)
340   CONTINUE
C
C. .. ADD CURRENT SOURCE CONTRIBUTION
450   IF (NCSC.EQ.0) GO TO 500
      DO 460 I=1,NCSC
      N1=NCS(I,1)
      N2=NCS(I,2)
      V(N1)=V(N1)-CSRC(I)
460   V(N2)=V(N2)+CSRC(I)
C
500   RETURN
      END
C
      SUBROUTINE SOLUTN
C
C     .... SOLVES EQUATIONS WITH SPECIAL MATRIX INDICATORS
C
      COMMON Y(600),YL(600),YV(200),YD(250),V(250)
      COMMON /SPARSE/ IORD(300),NUR(300),NVV(16),IPOS(1600),IUC(600),
     1 IV(200)
      DIMENSION NLC(300),ILR(600)
      EQUIVALENCE (NUR,NLC),(IUC,ILR)
      COMMON /NODES/ NREN(100),LNOD,NX,MDIM,NNOD,NUT,NLT,NVT
      COMMON /VSRCE/ VSRC(10),NAMV(10),TVLT(10,10),VVAL(10,10),
     1 NVVAL(10),KVOL(10),NV(10,2),NVSC
C
C     ... PROCESS VOLTAGE SOURCES
      IF (NVSC.EQ.0) GO TO 120
      DO 115 I=1,NVSC
      TEM=VSRC(I)
      V(MDIM+I)=TEM
      IF (TEM.EQ.0.0) GO TO 115
      NL=NVV(I)
      NLE=NVV(I+1)
105   IF (NL-NLE) 110,115,115
110   IR=IV(NL)
      V(IR)=V(IR)-YV(NL)*TEM
      NL=NL+1
      GO TO 105
115   CONTINUE
C
120   V(NNOD)=0.
      YD(250)=0.
      MD1=MDIM-1
      IF (MD1) 190,125,130
125   V(1)=V(1)/YD(1)
      GO TO 190
C
C     ... L*U DECOMPOSITION
130   KNT=0
      DO 155 I=1,MD1
      L=IORD(I)
      NUS=NUR(I)
      NUE=NUR(I+1)
      NL=NLC(I)
      NLE=NLC(I+1)
135   IF (NL-NLE) 140,155,155
140   TEM=YL(NL)/YD(L)
      YL(NL)=TEM
      NL=NL+1
      NU=NUS
145   IF (NU-NUE) 150,135,135
150   KNT=KNT+1
      K=IPOS(KNT)
      Y(K)=Y(K)-Y(NU)*TEM
      NU=NU+1
      GO TO 145
155   CONTINUE
C
C     ... FORWARD SUBSTITUTION
      DO 170 I=1,MD1
      L=IORD(I)
      NL=NLC(I)
      NLE=NLC(I+1)
160   IF (NL-NLE) 165,170,170
165   IR=ILR(NL)
      V(IR)=V(IR)-YL(NL)*V(L)
      NL=NL+1
      GO TO 160
170   CONTINUE
C
C     ... BACK SUBSTITUTION
      L=IORD(MDIM)
      V(L)=V(L)/YD(L)
      DO 185 I=1,MD1
      J=MDIM-I
      L=IORD(J)
      NU=NUR(J)
      NUE=NUR(J+1)
175   IF (NU-NUE) 180,185,185
180   IC=IUC(NU)
      V(L)=V(L)-Y(NU)*V(IC)
      NU=NU+1
      GO TO 175
185   V(L)=V(L)/YD(L)
C
190   RETURN
      END                                                               -
C
      SUBROUTINE INITL
C
C     .... RENUMBER NODES AND SET UP SPARSE MATRIX INDICATORS
C
      COMMON NROW(300),NCON(300),NSR(300),ISC(1600),NLC(300),ILR(600)
      COMMON /SPARSE/ IORD(300),NUR(300),NVV(16),IPOS(1600),IUC(600),
     1 IV(200)
      COMMON /NODES/ NREN(100),LNOD,NX,MDIM,NNOD,NUT,NLT,NVT
      COMMON /TRANS/ MREF(100),NQ(100,17),VT,TOLER,NXTR,ITR,IND
      COMMON /MODEL/ NTYP(30),BETA(30,8),CSAT(30,1),RBS(30,3),RCL(30,3),
     1 ROU(30,2),FT(30,4),CJ1(30,4),CJ2(30,6),TEMP(30)
      COMMON /ELMTS/ KIND(100),ND(100,4),VALU(100),NELT
      COMMON /VSRCE/ VSRC(10),NAMV(10),TVLT(10,10),VVAL(10,10),
     1 NVVAL(10),KVOL(10),NV(10,2),NVSC
      COMMON /CSRCE/ CSRC(10),NAMC(10),TCUR(10,10),CVAL(10,10),
     1 NCVAL(10),KCUR(10),NCS(10,2),NCSC
      DOUBLE PRECISION TITLE
      COMMON /JOBCNT/ TITLE(9),PRGVER,KBREF,IREQ,ISTOP
      DIMENSION IFM(2)
C
C     .... SET UP NODE CHECKING AND RENUMBERING ARRAYS
      LNOD=0
      DO 105 I=1,100
      NREN(I)=0
105   NCON(I)=0
C     .... SET UP NODE INDICATOR FOR GROUND AND V-SOURCES
      NREN(1)=1
      IF (NVSC.EQ.0) GO TO 115
      DO 110 I=1,NVSC
      N1=NV(I,1)+1
      NCON(N1)=NCON(N1)+1
      NCON(1)=NCON(1)+1
      IF (NREN(N1).EQ.0) GO TO 110
      WRITE (6,420) NV(I,1)
      ISTOP=1
110   NREN(N1)=1
C     .... SET UP INDICATOR FOR CURRENT SOURCES
115   IF (NCSC.EQ.0) GO TO 125
      DO 120 I=1,NCSC
      N1=NCS(I,1)+1
      N2=NCS(I,2)+1
      NCON(N1)=NCON(N1)+1
      NCON(N2)=NCON(N2)+1
120   LNOD=MAX0(LNOD,N1,N2)
C     .... SET UP INDICATORS FOR ELEMENTS
125   IF (NELT.EQ.0) GO TO 135
      DO 130 I=1,NELT
      N1=ND(I,1)+1
      N2=ND(I,2)+1
      NCON(N1)=NCON(N1)+1
      NCON(N2)=NCON(N2)+1
130   LNOD=MAX0(LNOD,N1,N2)
C     .... SET UP INDICATORS FOR TRANSISTORS
135   IF (NXTR.EQ.0) GO TO 145
      DO 140 I=1,NXTR
      N2=NQ(I,2)+1
      N3=NQ(I,3)+1
      NCON(N2)=NCON(N2)+1
      NCON(N3)=NCON(N3)+1
      N1=NQ(I,1)+1
      NCON(N1)=NCON(N1)+1
140   LNOD=MAX0(LNOD,N1,N2,N3)
C
C     ... ESTABLISH COMPACT NODE NUMBERS
145   MDIM=0
      IF (NCON(1).LE.1) GO TO 390
      DO 160 I=2,LNOD
      IF (NREN(I).NE.0) GO TO 160
      IF (NCON(I)-1) 160,150,155
150   J=I-1
      WRITE (6,410) J
      ISTOP=1
155   MDIM=MDIM+1
      NREN(I)=MDIM
160   CONTINUE
C     ... CREATE TRANSISTOR INTERNAL NODES
      NX=MDIM
      IF (NXTR.EQ.0) GO TO 185
      DO 180 I=1,NXTR
      K=MREF(I)
      IF (RCL(K,1).EQ.0.) GO TO 165
      MDIM=MDIM+1
      NQ(I,10)=MDIM
      GO TO 170
165   NQ(I,10)=0
170   IF (RBS(K,1).EQ.0.) GO TO 175
      MDIM=MDIM+1
      NQ(I,11)=MDIM
      GO TO 180
175   NQ(I,11)=0
180   CONTINUE
C
C     ... VOLTAGE SOURCE NODES
185   NNOD=MDIM
      IF (NVSC.EQ.0) GO TO 195
      DO 190 I=1,NVSC
      N1=NV(I,1)+1
      NNOD=NNOD+1
190   NREN(N1)=NNOD
C     ... GROUND NODE
195   NNOD=NNOD+1
      NREN(1)=NNOD
      IF (NNOD.GT.250) GO TO 395
C
C     ... RENUMBER NODES
      IF (NELT.EQ.0) GO TO 205
      DO 200 I=1,NELT
      N1=ND(I,1)+1
      N2=ND(I,2)+1
      ND(I,1)=NREN(N1)
      ND(I,2)=NREN(N2)
200   CONTINUE
205   IF (NXTR.EQ.0) GO TO 220
      DO 215 I=1,NXTR
      N2=NQ(I,2)+1
      N3=NQ(I,3)+1
      NQ(I,2)=NREN(N2)
      NQ(I,3)=NREN(N3)
      N1=NQ(I,1)+1
      NQ(I,1)=NREN(N1)
      IF (NQ(I,10).NE.0) GO TO 210
      NQ(I,10)=NQ(I,1)
210   IF (NQ(I,11).NE.0) GO TO 215
      NQ(I,11)=NQ(I,2)
215   CONTINUE
220   IF (NCSC.EQ.0) GO TO 230
      DO 225 I=1,NCSC
      N1=NCS(I,1)+1
      N2=NCS(I,2)+1
      NCS(I,1)=NREN(N1)
225   NCS(I,2)=NREN(N2)
C
C     ... SET UP MATRIX STRUCTURE INDICATORS
230   ND1=NNOD-1
      NR=0
      DO 365 I=1,ND1
      NSR(I)=NR+1
      IF (NELT.EQ.0) GO TO 270
      DO 265 J=1,NELT
      N1=ND(J,1)
      N2=ND(J,2)
      IF (N1.EQ.I) GO TO 235
      IF (N2.NE.I) GO TO 265
      IFN=N1
      GO TO 240
235   IFN=N2
240   IF (IFN-MDIM) 245,245,265
245   IF (IFN.EQ.I) GO TO 265
      N=NSR(I)
250   IF (N-NR) 255,255,260
255   IF (IFN.EQ.ISC(N)) GO TO 265
      N=N+1
      GO TO 250
260   NR=NR+1
      ISC(NR)=IFN
265   CONTINUE
C
270   IF (NXTR.EQ.0) GO TO 365
      DO 360 J=1,NXTR
      N1=NQ(J,1)
      N2=NQ(J,2)
      N3=NQ(J,3)
      NC=NQ(J,10)
      NB=NQ(J,11)
      NX=-1
      IF (N1.EQ.I) GO TO 275
271   NX=0
      IF (N2.EQ.I) GO TO 285
      IF (NC.EQ.I) GO TO 295
      IF (NB.EQ.I) GO TO 300
273   NX=1
      IF (N3.NE.I) GO TO 360
      IFN=NC
      IFM(1)=NB
      GO TO 310
275   IF (NC.EQ.N1) GO TO 280
      IFN=NC
      GO TO 315
280   IFN=NB
      IFM(1)=N3
      GO TO 310
285   IF (NB.EQ.N2) GO TO 290
      IFN=NB
      GO TO 315
290   IFN=NC
      IFM(1)=N3
      GO TO 310
295   IFN=N1
      IFM(1)=NB
      GO TO 305
300   IFN=N2
      IFM(1)=NC
305   IFM(2)=N3
      M=-3
      GO TO 320
310   M=-2
      GO TO 320
315   M=0
      GO TO 335
320   K=0
325   M=M+1
      IF (IFN-MDIM) 330,330,355
330   IF (IFN.EQ.I) GO TO 355
335   N=NSR(I)
340   IF (N-NR) 345,345,350
345   IF (IFN.EQ.ISC(N)) GO TO 355
      N=N+1
      GO TO 340
350   NR=NR+1
      ISC(NR)=IFN
355   IF (M.EQ.0) GO TO 357
      K=K+1
      IFN=IFM(K)
      GO TO 325
357   IF (NX) 271,273,360
360   CONTINUE
365   CONTINUE
      NSR(NNOD)=NR+1
C
C     ... DETERMINE INDICATORS FOR MATRIX REDUCTION
      CALL INDMT (ISTOP)
      IF (ISTOP.NE.0) GO TO 405
C
C     ... ASSIGN MATRIX LOCATION NUMBERS
      IF (NELT.EQ.0) GO TO 375
      DO 370 I=1,NELT
      N1=ND(I,1)
      N2=ND(I,2)
      ND(I,3)=INDX(N1,N2)
      ND(I,4)=INDX(N2,N1)
370   CONTINUE
375   IF (NXTR.EQ.0) GO TO 385
      DO 380 I=1,NXTR
C     ... INTRINSIC XTR
      NC=NQ(I,10)
      NB=NQ(I,11)
      N3=NQ(I,3)
      NQ(I,4)=INDX(NC,NB)
      NQ(I,5)=INDX(NC,N3)
      NQ(I,6)=INDX(NB,NC)
      NQ(I,7)=INDX(NB,N3)
      NQ(I,8)=INDX(N3,NC)
      NQ(I,9)=INDX(N3,NB)
C     ... RESISTANCES
      N1=NQ(I,1)
      N2=NQ(I,2)
      NQ(I,12)=INDX(N1,NC)
      NQ(I,13)=INDX(NC,N1)
      NQ(I,14)=INDX(N2,NB)
      NQ(I,15)=INDX(NB,N2)
C     ... SPLIT CAPACITANCE
      NQ(I,16)=INDX(N2,NC)
      NQ(I,17)=INDX(NC,N2)
380   CONTINUE
385   RETURN
C
390   WRITE (6,415)
      GO TO 400
395   WRITE (6,425)
400   ISTOP=1
405   RETURN
C
410   FORMAT (/1X,8H****NODE,I3,12H IS SINGULAR/)
415   FORMAT (/1X,28H****CIRCUIT NOT GROUNDED****/)
420   FORMAT (/1X,30H****MULTIPLE V-SOURCES AT NODE,I3/)
425   FORMAT (/1X,31H****NUMBER OF NODES EXCEEDS 250/)
      END
C
      SUBROUTINE INDMT (ISTOP)
C
C     .... ESTABLISH MATRIX REDUCTION INDICATORS
C
      COMMON NROW(300),NCON(300),NSR(300),ISC(1600),NLC(300),ILR(600)
      COMMON /SPARSE/ IORD(300),NUR(300),NVV(16),IPOS(1600),IUC(600),
     1 IV(200)
      COMMON /NODES/ NREN(100),LNOD,NX,MDIM,NNOD,NUT,NLT,NVT
      COMMON /TRANS/ MREF(100),NQ(100,17),VT,TOLER,NXTR,ITR,IND
      COMMON /ELMTS/ KIND(100),ND(100,4),VALU(100),NELT
      COMMON /VSRCE/ VSRC(10),NAMV(10),TVLT(10,10),VVAL(10,10),
     1 NVVAL(10),KVOL(10),NV(10,2),NVSC
C
C     ... SET UP INDICATORS FOR V-SOURCES
      NVT=0
      IF (NVSC.EQ.0) GO TO 120
      NL=NSR(MDIM+1)-1
      NU=NVSC+1
      DO 105 I=1,NU
105   NVV(I)=NSR(MDIM+I)-NL
      N=NL+1
      NLE=NSR(NNOD)-1
110   IF (N-NLE) 115,115,120
115   NVT=NVT+1
      IV(NVT)=ISC(N)
      N=N+1
      GO TO 110
C
C     ... COUNT OFF DIAGONAL ELEMENTS
120   DO 125 I=1,MDIM
      IORD(I)=I
      NCON(I)=I
      NROW(I)=NSR(I+1)-NSR(I)
125   CONTINUE
C
C     ... COLUMN AND ROW RENUMBERING AND INDICATOR SET-UP
      NU=1
      NL=1
      NPS=0
      MP1=MDIM+1
      MD1=MDIM-1
      IF (MD1) 215,215,130
130   DO 210 I=1,MD1
      NUR(I)=NU
      NLC(I)=NL
      JS=I+1
      L=IORD(I)
C     ... SEARCH FOR MINIMUM  NROW
      DO 135 J=JS,MDIM
      IR=IORD(J)
      IF (NROW(IR).GE.NROW(L)) GO TO 135
      IORD(J)=L
      NCON(L)=J
      IORD(I)=IR
      NCON(IR)=I
      L=IR
135   CONTINUE
      NPAS=0
      DO 205 J=JS,MDIM
      IR=IORD(J)
      IF (LOC(IR,L).EQ.0) GO TO 205
      ILR(NL)=IR
      NROW(IR)=NROW(IR)-1
      NL=NL+1
      IF (NPAS.NE.0) GO TO 165
      DO 160 K=JS,MDIM
      IC=IORD(K)
      IF (LOC(L,IC).EQ.0) GO TO 160
      IUC(NU)=IC
      IF (IR.EQ.IC) GO TO 155
      IF (LOC(IR,IC).NE.0) GO TO 155
      NROW(IR)=NROW(IR)+1
      IRS=IR+1
      DO 140 IT=IRS,MP1
140   NSR(IT)=NSR(IT)+1
      NS=NSR(IRS)
      N=NSR(MP1)
143   N=N-1
      IF (N-NS) 150,145,145
145   ISC(N)=ISC(N-1)
      GO TO 143
150   ISC(N)=IC
155   NU=NU+1
      NPS=NPS+1
160   CONTINUE
      NPAS=1
      GO TO 205
165   NCT=NUR(I)
170   IF (NCT-NU) 175,205,205
175   IC=IUC(NCT)
      IF (IR.EQ.IC) GO TO 200
      IF (LOC(IR,IC).NE.0) GO TO 200
      NROW(IR)=NROW(IR)+1
      IRS=IR+1
      DO 180 IT=IRS,MP1
180   NSR(IT)=NSR(IT)+1
      NS=NSR(IRS)
      N=NSR(MP1)
185   N=N-1
      IF (N-NS) 195,190,190
190   ISC(N)=ISC(N-1)
      GO TO 185
195   ISC(N)=IC
200   NCT=NCT+1
      NPS=NPS+1
      GO TO 170
205   CONTINUE
210   CONTINUE
C
215   NLC(MDIM)=NL
      NUR(MDIM)=NU
      NUT=NU
      NLT=NL
      WRITE (6,260) NUT,NLT,NVT,NPS
      IF (NUT.GT.600) GO TO 250
      IF (NPS.GT.1600) GO TO 250
C
C     ... ASSIGN OPERATION NUMBERS
      IF (MD1.EQ.0) GO TO 245
      NPS=0
      DO 240 I=1,MD1
      NUS=NUR(I)
      NUE=NUR(I+1)
      NL=NLC(I)
      NLE=NLC(I+1)
220   IF (NL-NLE) 225,240,240
225   IR=ILR(NL)
      NL=NL+1
      NU=NUS
230   IF (NU-NUE) 235,220,220
235   IC=IUC(NU)
      NPS=NPS+1
      IPOS(NPS)=INDX(IR,IC)
      NU=NU+1
      GO TO 230
240   CONTINUE
245   RETURN
C
250   WRITE (6,255)
      ISTOP=1
      RETURN
C
C
255   FORMAT (/1X,24H****MATRIX TOO DENSE****/)
260   FORMAT (/1X,4HNU =,I4,6H  NL =,I4,7H   NV =,I4,7H  NPS =,I5/)
      END                                                               -
C
      FUNCTION INDX (IR,IC)
C
C     .... OBTAINS MATRIX LOCATIONS
C
      COMMON NROW(300),NCON(300),NSR(300),ISC(1600),NLC(300),ILR(600)
      COMMON /SPARSE/ IORD(300),NUR(300),NVV(16),IPOS(1600),IUC(600),
     1 IV(200)
      COMMON /NODES/ NREN(100),LNOD,NX,MDIM,NNOD,NUT,NLT,NVT
C
      IF (IR.EQ.IC) GO TO 175
      IF (IR-MDIM) 50,50,180
50    IF (IC-NNOD) 60,180,180
60    IF (IC-MDIM) 70,70,150
70    IS=NCON(IR)
      JS=NCON(IC)
      IF (JS.GE.IS) GO TO 120
C     ... LOWER TRIANGLE
      N=NLC(JS)
      NE=NLC(JS+1)
100   IF (N-NE) 105,180,180
105   IF (IR.EQ.ILR(N)) GO TO 110
      N=N+1
      GO TO 100
110   INDX=600+N
      RETURN
C     ... UPPER TRIANGLE
120   N=NUR(IS)
      NE=NUR(IS+1)
130   IF (N-NE) 135,180,180
135   IF (IC.EQ.IUC(N)) GO TO 140
      N=N+1
      GO TO 130
140   INDX=N
      RETURN
C     ... VOLTAGE SOURCE NODES
150   K=IC-MDIM
      N=NVV(K)
      NE=NVV(K+1)
160   IF (N-NE) 165,180,180
165   IF (IR.EQ.IV(N)) GO TO 170
      N=N+1
      GO TO 160
170   INDX=1200+N
      RETURN
C     ... DIAGONAL ELEMENT
175   INDX=1400+IR
      RETURN
C     ... UNUSED LOCATION
180   INDX=1400+NNOD
      RETURN
      END
      FUNCTION LOC (IR,IC)
C
C     ... TEST FOR NON-ZERO MATRIX ELEMENT
C
      COMMON NROW(300),NCON(300),NSR(300),ISC(1600)
      N=NSR(IR)
      NE=NSR(IR+1)
200   IF (N-NE) 210,220,220
210   IF (IC.EQ.ISC(N)) GO TO 290
      N=N+1
      GO TO 200
220   LOC=0
      RETURN
290   LOC=1
      RETURN
      END
C
      SUBROUTINE PLOT
C
C.... SUBROUTINE TO PRINT TIME RESPONSE PLOTS
C
      COMMON Y(600),YL(600),YV(200),YD(250),V(250),V1(250),BCUR(100),
     1 FA(50),FB(50),FC(50),FFA(50),FCR(50),VFNR(200),LINE(101),DM(449),
     1 PLTP(2,200,6)
      COMMON /PLOTS/ PMIN(6),PMAX(6),KPLO(2,6),NPVAL(6),NAMP(2,6),
     1 NP(2,2,6),NPLT,NPTS,FSTR,FSTP,NHAR
      COMMON /TIMEV/ TTIME,DELT,TIMP,DELP,TMAX,TPRT,DPRT,ITER,NCUT,ICONV
      DOUBLE PRECISION TITLE
      COMMON /JOBCNT/ TITLE(9),PRGVER,KBREF,IREQ,ISTOP
      DIMENSION LIN1(50),LIN2(51)
      EQUIVALENCE (LIN1(1),LINE(1)),(LIN2(1),LINE(51))
      DIMENSION VAL(11),NUM(5),ISTOR(5)
      DATA NUM/1H+,1H*,1H3,1H4,1H5/, KBL,IPL/1H ,1H./
      PLTIM(N)=DPRT*FLOAT(N-1)
C
      IF (NPTS.GT.200) NPTS=200
      IF (IREQ.EQ.0) GO TO 102
      IF (NHAR.GT.49) NHAR=49
      IRM=1
      INRI=0
      DPRD=DPRT/37.3
      FSTP=FSTP+DPRD
      FSTR=FSTR+DPRD
C     LIMITS OF PLOT PROGRAM FOR FOURIER ANALYSIS
      DO 407 KPT=1,NPTS
      IF(PLTIM(KPT).GT.FSTR) GO TO 420
  407 CONTINUE
      GO TO 560
  420 KMIN=KPT-1
      DO 430 KPT=KMIN,NPTS
      IF (PLTIM(KPT).GT.FSTP) GO TO 440
430   CONTINUE
      GO TO 560
  440 KMAY=KPT-1
      KMAN=KMAY-1
      GO TO 102
  560 INRI=1
      KPTF=KPT
      FKTP=PLTIM(KPT)
C
C.... SEQUENCE THROUGH PLOTS
  102 DO 210 KPLT=1,NPLT
C.... PRINT OUTPUTS
      WRITE (6,250) TITLE,PRGVER
      NVL=NPVAL(KPLT)
      NL=MIN0(2,NVL)
      DO 100 I=1,NVL
100   WRITE (6,260) NUM(I),NAMP(I,KPLT),(NP(J,I,KPLT),J=1,2)
      WRITE (6,260)
C.... CALCULATE AND PRINT COORDINATE VALUES
      PMN=PMIN(KPLT)
      PMX=PMAX(KPLT)
      WID=(PMX-PMN)/100.
      IF (WID.GT.0.0) GO TO 105
      PMN=-10.0
      WID=0.2
 105  DO 110 I=1,11
110   VAL(I)=FLOAT(I-1)*(WID*10.)+PMN
      WRITE (6,270) VAL
      DO 120 I=2,101
120   LINE(I)=KBL
C
C.... PRINT PLOT
      PRT=0.
      DO 200 KPT=1,NPTS
      MIN=101
      MAX=1
      LINE(1)=IPL
      LINE(51)=IPL
      DO 140 I=1,NVL
      IPOS=INT((PLTP(I,KPT,KPLT)-PMN)/WID+1.5)
      IF (IPOS.GT.101) GO TO 130
      IF (IPOS.LT.1) GO TO 130
      LINE(IPOS)=NUM(I)
      ISTOR(I)=IPOS
      MIN=MIN0(MIN,IPOS)
      MAX=MAX0(MAX,IPOS)
      GO TO 140
 130  ISTOR(I)=1
140   CONTINUE
C
C.... PRINT LINE
      IF (MIN.GT.50) GO TO 150
      IF (MAX.LE.50) GO TO 160
C.... PRINT ALL ACROSS
      WRITE (6,280) PRT,LINE
      GO TO 180
C.... PRINT LAST 50 POSITIONS
150   WRITE (6,290) PRT,LIN2
      GO TO 180
C.... PRINT FIRST 50 POSITIONS
160   WRITE (6,300) PRT,LIN1
 180  WRITE (6,305) (PLTP(I,KPT,KPLT),I=1,NL)
      DO 190 I=1,NVL
      N=ISTOR(I)
190   LINE(N)=KBL
      PRT=PRT+DPRT
200   CONTINUE
      WRITE (6,310)
C
C     ... PROGRAM TO PERFORM FOURIER ANALYSIS
      IF (IREQ.EQ.0) GO TO 209
      IF(INRI.EQ.0) GO TO 445
      PRINT 508,KPTF,FKTP
      GO TO 209
  445 INU=0
      DO 205 II=1,NL
C        NUFU = 2N + 1    IS ODD
      NUFU=KMAY-KMIN
      IF((NUFU-NUFU/2-NUFU/2).EQ.0)GO TO 470
      FDIV=(FSTP-FSTR)/NUFU
      DO 460 KPT=KMIN,KMAN
      DPPI=FSTR+(KPT-KMIN)*FDIV-PLTIM(KPT)
      INU=INU+1
  460 VFNR(INU)=PLTP(II,KPT,KPLT)+(PLTP(II,KPT+1,KPLT)-PLTP(II,KPT,KPLT)
     1)*DPPI/DPRT
      GO TO 480
C
  470 NUFU=NUFU-1
      FDIV=(FSTP-FSTR)/NUFU
      DPPI=FSTR-FDIV-PLTIM(KMIN-1)
      FDIX=DPRT+DPRT-FDIV
      KMO=0
      KMQ=0
      DO 462 KPT=KMIN,KMAN
      IF(KMO.EQ.1)GO TO 466
      DPPI=DPPI+FDIV-DPRT*(1+KMQ)
      KMQ=0
      INU=INU+1
      VFNR(INU)=PLTP(II,KPT,KPLT)+(PLTP(II,KPT+1,KPLT)-PLTP(II,KPT,KPLT)
     1)*DPPI/DPRT
      DIFX=FDIX-DPPI
      IF (DIFX.LE.0.0) KMO=1
      GO TO 467
  466 KMO=0
      KMQ=1
  467 CONTINUE
  462 CONTINUE
  480 IF(NUFU.EQ.INU)GO TO 469
      PRINT 502,NUFU,INU
      GO TO 209
  469 CONTINUE
      PRINT 501,NHAR,FSTR,FSTP,NUFU
C
C     ... FOURIER ANALYSIS AND OUTPUT CALCULATIONS
      NUFV=(NUFU-1)/2
      CALL FORIT(VFNR,NUFV,NHAR,FA,FB)
      PRINT 504
      NHAS=NHAR+1
      DO 510 I=1,NHAS
      FC(I)=SQRT(FA(I)*FA(I)+FB(I)*FB(I))
  510 FFA(I)=ATAN2(-FB(I),FA(I))*180./3.1415926
      DO 520 I=1,NHAS
      IPOR=I-1
      FCR(I)=100.*FC(I)/FC(IRM+1)
  520 PRINT 506,IPOR,FB(I),FA(I),FC(I),FCR(I),FFA(I)
 205  CONTINUE
  209 CONTINUE
210   CONTINUE
      RETURN
C
 250  FORMAT (1H7/1H1,9A8,8X,5HSINC-,A3//1HX)
 260  FORMAT (1X,A1,1H%,1X,A4,2I3)
 270  FORMAT (3X,11F10.4/10X,1HI,10(9X,1HI)/10X,101(1H.))
 280  FORMAT (1X,1PE8.2,1X,101A1)
 290  FORMAT (1X,1PE8.2,1X,1H.,49X,51A1)
 300  FORMAT (1X,1PE8.2,1X,50A1,1H.)
 305  FORMAT (1H+,110X,2F10.5)
 310  FORMAT (1HY)
  501 FORMAT(1X,' FOURIERANALYSIS IN ',I4,' HARMONICS FROM',E11.4,' TILL
     1',E11.4,' SECONDS CALCULATED WITH',I5,' POINTS '/)
  502 FORMAT(5X,' ERROR = GIVEN NUMBER OF INTERPOLATIONPOINTS',I6,' DIFF
     1ERS FROM COMPUTED NUMBER ',I6/)
  504 FORMAT(1X,'ORDER HARM',3X,'SINE',8X,'COSINE',15X,'MAGNITUDE  RELAT
     1 MAG   PHASE'/)
  506 FORMAT(7X,I2,3X,2E12.4,9X,1E11.4,2F11.4)
  508 FORMAT(5X,'ERROR = TIME LIMITS OF FOUR OUTSIDE COMPUTED OUTPUTTIME
     1          POINT IS ',I3,' TIME IS ',1PE8.2)
      END

      SUBROUTINE FORIT (FNT,N,M,A,B)
C
      DIMENSION FNT(1),A(M),B(M)
C
C     ... SET CONSTANTS
      COEF=1.0/(FLOAT(N)+0.5)
      CONST=3.141593*COEF
      S1=SIN(CONST)
      C1=COS(CONST)
      C=1.0
      S=0.0
      FNTZ=FNT(1)
      J=1
      M1=M+1
      N2=N+N+1
C
C     ... COMPUTE FOURIER COEFFICIENTS RECURSIVELY
   70 U2=0.0
      U1=0.0
      I=N2
   75 U0=FNT(I)+2.0*C*U1-U2
      U2=U1
      U1=U0
      I=I-1
      IF (I-1) 80,80,75
   80 A(J)=COEF*(FNTZ+C*U1-U2)
      B(J)=COEF*S*U1
      IF (J-M1) 90,100,100
   90 Q=C1*C-S1*S
      S=C1*S+S1*C
      C=Q
      J=J+1
      GO TO 70
  100 A(1)=A(1)*0.5
      RETURN
      END
C
      SUBROUTINE MODELL(IDEC)
C
C     .... READS TRANSISTOR MODEL CARD
C
      COMMON /MODEL/ NTYP(30),BETA(30,8),CSAT(30,1),RBS(30,3),RCL(30,3),
     1 ROU(30,2),FT(30,4),CJ1(30,4),CJ2(30,6),TEMP(30)
      COMMON /MODPR/ BEF(30),BER(30),CSA(30),TCF(30),TCR(30),GBS(30),
     1 GCL(30),BWF(30),BC1(30),BC2(30),BC3(30),CJE(30),CJC(30),CT1(30),
     2 CT2(30)
      COMMON /NAMES/ NAME(100),NAMQ(100),MODN(100),NAMB(30),NMOD
      DOUBLE PRECISION TITLE
      COMMON /JOBCNT/ TITLE(9),PRGVER,NALTR,IREQ,ISTOP
      COMMON /CARD/ KARD(80),IPNT(80),IPT,NERR,KFLG,NAM
      DOUBLE PRECISION NAMP
      DIMENSION MPL(20),DFVAL(32),LOCAT(13),NDUM(5),NAMP(35)
      DATA MPL /3HNPN,3HPNP,2HBF,2HBR,3HISS,2HRB,2HRC,2HRO,2HFT,3HTSA,
     1 3HCJE,3HCJC,3HCSU,3HTEM,2HTF,2HTR,2HVA,2HIS,2*0/
      DATA DFVAL /100.,6*0.,1.,1.E-14,6*0.,0.,1.E-3,4*0.,
     1 0.,0.,.7,.33333,0.,0.,.5,.33333,0.,0.,300./
      DATA LOCAT /1,8,9,10,13,16,18,21,22,26,31,32,33/
      DATA KBL,KPL /1H ,1H+/
      DATA NAMP /8HBF (MAX),8H  /ICMAX,6H  /BF ,6H  /IC ,6H  /VCE,
     1 6H  /TC1,6H  /TC2,6HBR    ,6HISS   ,6HRB    ,6H TC1  ,6H TC2  ,
     16HRC    ,6H TC1  ,
     2 6H TC2  ,6HRO    ,6H   IC ,6HFT    ,6H   IC ,6H   VCE,6HTSAT  ,
     2 6HCJE   ,6H   VBE,
     3 7H   PHIE,6H   NE ,6HCJC   ,6H   VBC,7H   PHIC,6H   NC ,
     3 8H   RATIO,
     4 6HCSUB  ,6HTEMP  ,6HTF    ,6HTR    ,6H1/VA   /
C
      GO TO (1,2),IDEC
1     MFLG=0
      IF (NALTR.EQ.0) GO TO 230
      DO 225 I=1,NMOD
      IF (NAM.EQ.NAMB(I)) GO TO 250
225   CONTINUE
230   NMOD=NMOD+1
      IF (NMOD.GT.30) GO TO 600
      I=NMOD
      NAMB(I)=NAM
      DO 245 J=1,32
245   BETA(I,J)=DFVAL(J)
      NTYP(I)=1
      TCF(I)=-1.
      TCR(I)=-1.
      BWF(I)=-1.
C
C     ... INTERPRET CARD
250   N=0
265   N=N+1
267   NAM=KBL
270   VAL=VALUE(2)
      IF (NAM.EQ.0) GO TO 265
      IF (NAM.EQ.1) GO TO 335
      IF (NAM.EQ.KBL) GO TO 332
C     ... SET CODE NUMBER
300   DO 310 J=1,18
      IF (NAM.EQ.MPL(J)) GO TO 320
310   CONTINUE
      NAM=KBL
      VAL=VALUE(2)
      IF (NAM.EQ.1) GO TO 335
      IF (NAM.EQ.KBL) GO TO 270
      GO TO 300
320   KODE=J-2
      IF (KODE) 321,321,323
321   NTYP(I)=IABS(KODE)
      GO TO 265
323   VAL=VALUE(1)
C
C     ... SET MODEL VALUES
      IF (KODE-13) 330,325,326
325   TCF(I)=VAL
      GO TO 265
326   IF (KODE-15) 327,328,329
327   TCR(I)=VAL
      GO TO 265
328   BWF(I)=VAL
      GO TO 265
329   CSAT(I,1)=VAL
      GO TO 265
C
C     ... SET PARAMETER VALUE
330   N=LOCAT(KODE)
      GO TO 333
332   IF (N.GE.33) GO TO 265
333   BETA(I,N)=VAL
      GO TO 265
335   IF (NERR.EQ.0) GO TO 345
      MFLG=1
      WRITE (6,920) IPNT
      WRITE (6,925) NERR
      DO 340 J=1,80
340   IPNT(J)=KBL
      NERR=0
345   READ (1,905,END=347) KARD
      GO TO 410
347   KFLG=-1
      GO TO 450
410   IF (KARD(1).NE.KPL) GO TO 420
      IPT=2
      WRITE (6,920) KARD
      GO TO 265
420   KFLG=1
450   IF (MFLG.EQ.0) GO TO 460
      NMOD=NMOD-1
      GO TO 600
C
C     ... CHECK PARAME3ER VALUES
460   IF (BETA(I,1).LE.0.) BETA(I,1)=DFVAL(1)
      IF (BETA(I,8).LE.0.) BETA(I,1)=DFVAL(8)
      IF (CSAT(I,1).LE.0.) CSAT(I,1)=DFVAL(9)
      RBS(I,1)=ABS(RBS(I,1))
      RCL(I,1)=ABS(RCL(I,1))
      ROU(I,1)=ABS(ROU(I,1))
      CJ1(I,1)=ABS(CJ1(I,1))
      CJ2(I,1)=ABS(CJ2(I,1))
      CJ2(I,6)=ABS(CJ2(I,6))
      IF (CJ1(I,3).LE.0.) CJ1(I,3)=DFVAL(24)
      IF (CJ1(I,4).LE.0.) CJ1(I,4)=DFVAL(25)
      IF (CJ2(I,3).LE.0.) CJ2(I,3)=DFVAL(28)
      IF (CJ2(I,4).LE.0.) CJ2(I,4)=DFVAL(29)
      IF (TEMP(I).LT.100.) TEMP(I)=100.
      K=I
C     ... JUNSTION CAPACITANCES
      CT1(K)=1.5*(.5**CJ1(K,4))
      VQT=CJ1(K,2)/CJ1(K,3)
      IF (VQT.GE..5) GO TO 350
      CJE(K)=CJ1(K,1)*(1.-VQT)**CJ1(K,4)
      GO TO 355
350   CJE(K)=CJ1(K,1)*CT1(K)/(1.+VQT)
355   CT2(K)=1.5*(.5**CJ2(K,4))
      VQT=CJ2(K,2)/CJ2(K,3)
      IF (VQT.GE..5) GO TO 360
      CJC(K)=CJ2(K,1)*(1.-VQT)**CJ2(K,4)
      GO TO 365
360   CJC(K)=CJ2(K,1)*CT2(K)/(1.+VQT)
C     ... TRANSIT TIMES
365   IF (TCF(K).GE.0.) GO TO 375
      TCF(K)=0.
      IF (FT(K,1).EQ.0.) GO TO 375
      TF=.159/FT(K,1)
      IF (FT(K,2).EQ.0.) GO TO 370
      GM=FT(K,2)/(8.6164E-5*TEMP(K))
      TF=TF-(CJE(K)+CJC(K)*(1.+GM*RCL(K,1)))/GM
370   IF (TF.LT.0.) TF=0.
      TCF(K)=TF
375   IF (TCR(K).GE.0.) GO TO 377
      BF=BETA(K,1)
      BR=BETA(K,8)
      TR=((1.+BF+BR)*FT(K,4)/BF-(1.+BR)*TCF(K))/BR
      IF (TR.LT.0.) TR=0.
      TCR(K)=TR
C     ... BETA COEFFICIENTS
377   IF (BETA(K,4).EQ.0.) GO TO 380
      IF (BETA(K,3).EQ.0.) GO TO 380
      BC2(K)=SQRT(BETA(K,4))/BETA(K,3)
      BC3(K)=BC2(K)/(2.*BETA(K,2)**1.5)
      BC1(K)=1./BETA(K,1)-1.89*(BC2(K)**.667)*(BC3(K)**.333)
      GO TO 400
380   BC1(K)=1./BETA(K,1)
      BC2(K)=0.
      BC3(K)=0.
C     ... SET EARLY VOLTAGE
400   IF (BWF(K).GT.0.) GO TO 405
      BWF(K)=0.
      IF (ROU(K,1).EQ.0.) GO TO 600
      BWF(K)=1./(ROU(K,1)*ROU(K,2))
      GO TO 600
405   BWF(K)=1./BWF(K)
600   RETURN
C
C     ... PRINT MODEL PARAMETERS
2     IF (NMOD.EQ.0) GO TO 850
      WRITE (6,980)
      NE=0
605   NS=NE+1
      NE=MIN0(NMOD,NS+4)
      IF (NE.GT.30) GO TO 635
      IF (NS-NE) 610,610,850
610   WRITE (6,985) (NAMB(J),J=NS,NE)
      N=0
      DO 620 K=NS,NE
      N=N+1
      IF (NTYP(K).NE.0) GO TO 615
      NDUM(N)=MPL(2)
      GO TO 620
615   NDUM(N)=MPL(1)
620   CONTINUE
      WRITE (6,990) (NDUM(J),J=1,N)
      DO 625 KODE=1,12
      KS=LOCAT(KODE)
      KE=LOCAT(KODE+1)
      N=KS
621   DO 622 J=NS,NE
      IF (BETA(J,N).NE.0.) GO TO 623
622   CONTINUE
      IF (N-KS) 625,625,624
623   WRITE (6,995) NAMP(N),(BETA(J,N),J=NS,NE)
624   N=N+1
      IF (N-KE) 621,625,625
625   CONTINUE
      WRITE (6,995)
      WRITE (6,995) NAMP(33),(TCF(J),J=NS,NE)
      WRITE (6,995) NAMP(34),(TCR(J),J=NS,NE)
      WRITE (6,995) NAMP(35),(BWF(J),J=NS,NE)
      GO TO 605
C
635   ISTOP=1
      WRITE (6,1080)
850   RETURN
C
905   FORMAT (80A1)
920   FORMAT (3X,80A1)
925   FORMAT (1X,2H**,I3,23H ERROR(S) ON ABOVE CARD )
980   FORMAT (/1X,11HBJT MODELS-/1X,9HPARAMETER,4X,5HVALUE)
985   FORMAT (/1X,4HNAME,4X,5(4X,A5,3X))
990   FORMAT (1X,4HTYPE,4X,5(4X,A3,5X))
995   FORMAT (1X,A8,1P5E12.3)
1080  FORMAT (/1X,28H****PROGRAM STORAGE EXCEEDED/)
      END
C
      SUBROUTINE READIN
C
C.... SUBROUTINE TO READ AND PROCESS CIRCUIT DATA
C
      COMMON /TRANS/ MREF(100),NQ(100,17),VT,TOLER,NXTR,ITR,IND
      COMMON /ELMTS/ KIND(100),ND(100,4),VALU(100),NELT
      COMMON /VSRCE/ VSRC(10),NAMV(10),TVLT(10,10),VVAL(10,10),
     1 NVVAL(10),KVOL(10),NV(10,2),NVSC
      COMMON /CSRCE/ CSRC(10),NAMC(10),TCUR(10,10),CVAL(10,10),
     1 NCVAL(10),KCUR(10),NC(10,2),NCSC
      COMMON /PLOTS/ PMIN(6),PMAX(6),KPLO(2,6),NPVAL(6),NAMP(2,6),
     1 NP(2,2,6),NPLT,NPTS,FSTR,FSTP,NHAR
      COMMON /OUTPT/ KOUT(10),NAMO(10),NODO(10,2),NOUT
      COMMON /TIMEV/ TTIME,DELT,TIMP,DELP,TMAX,TPRT,DPRT,ITER,NCUT,ICONV
      COMMON /NAMES/ NAME(100),NAMQ(100),MODN(100),NAMB(30),NMOD
      COMMON /STORE/ LQ(100,3),LD(100,2),LC(10,2),XALU(100),LCL(10),
     1 LVL(10),TC1(100),TC2(100),TEMPE(10),NTMP,KTMP,NOIS,NPOS,NNEG
      DOUBLE PRECISION TITLE
      COMMON /JOBCNT/ TITLE(9),PRGVER,NALTR,IREQ,ISTOP
      COMMON /CARD/ KARD(80),IPNT(80),IPT,NERR,KFLG,NAM
      DIMENSION LETR(20)
      DATA LETR/1HR,1HC,1HL,1HQ,1HB,1HV,1HI,1H*,3HVOU,3HIOU,3HTIM,
     1 3HPLO,3HTOL,3HFOR,3HTEM,3HPRI,3HMOD,3HTRA,2*0/
      DATA KBL,IHI,IHV,ISN,KDT /1H ,1HI,1HV,1H$,1H./
      DATA KVOUT,IOUT,IHN,IHB /4HVOUT,4HIOUT,1HN,1HB/
C
      IF (NALTR.EQ.0) GO TO 100
      WRITE (6,90) TITLE,PRGVER
 90   FORMAT (1H1,9A8,8X,5HSINC-,A3///,2X,20HCONTINUATION OF DATA //)
      DO 20 I=1,NXTR
      NQ(I,1)=LQ(I,1)
      NQ(I,2)=LQ(I,2)
20    NQ(I,3)=LQ(I,3)
      DO 30 I=1,NELT
      ND(I,1)=LD(I,1)
30    ND(I,2)=LD(I,2)
      IF (NCSC.EQ.0) GO TO 45
      DO 40 I=1,NCSC
      NCVAL(I)=LCL(I)
      IF (KCUR(I).NE.2) GO TO 35
      CVAL(2,I)=CVAL(5,I)
      CVAL(3,I)=CVAL(6,I)
35    NC(I,1)=LC(I,1)
40    NC(I,2)=LC(I,2)
45    IF (NVSC.EQ.0) GO TO 110
      DO 50 I=1,NVSC
      NVVAL(I)=LVL(I)
      IF (KVOL(I).NE.2) GO TO 50
      VVAL(2,I)=VVAL(5,I)
      VVAL(3,I)=VVAL(6,I)
50    CONTINUE
      GO TO 110
C
C     ... INITIALIZE
100   NMOD=0
      NXTR=0
      NELT=0
      NVSC=0
      NCSC=0
      NSRC=0
      NSIN=0
      NOUT=0
      NPLT=0
      IREQ=0
      NOIS=0
      NTMP=0
      TMAX=0.
      TOLER=.1
      DPRT=1.
      DO 105 I=1,80
105   IPNT(I)=KBL
      KFLG=0
C
C.... READ CARD AND TEST TYPE
110   NERR=0
      IF (KFLG) 630,111,113
111   READ (1,910,END=630) KARD
113   KFLG=0
      IPT=1
      IF (KARD(1).EQ.KDT) IPT=2
      NAM=NVAL(3)
      IF (NAM.EQ.3HEEC) GO TO 629
      WRITE (6,1020) KARD
      IF (NAM.EQ.3HEND) GO TO 635
      IF (NAM.EQ.3HALT) GO TO 640
      IF (NAM.EQ.LETR(17)) GO TO 325
      IF (NAM.EQ.LETR(18)) GO TO 400
      DO 120 KODE=9,16
      IF (NAM.EQ.LETR(KODE)) GO TO 160
120   CONTINUE
      DO 125 KODE=1,8
      IF (KARD(1).EQ.LETR(KODE)) GO TO 130
125   CONTINUE
      WRITE (6,1010)
      GO TO 110
C
C     ... GET ELEMENT NAME
 130  IPT=1
      NAM=NVAL(5)
C
C.... BRANCH TO TYPE
160   CONTINUE
      GO TO (290,290,290,300,330,170,230,110,420,420,400,460,510,520,
     1 530,550), KODE
C
C.... VOLTAGE SOURCE LINE
170   IF (NALTR.EQ.0) GO TO 174
      DO 172 L=1,NVSC
      IF (NAM.EQ.NAMV(L)) GO TO 176
172   CONTINUE
174   NVSC=NVSC+1
      IF (NVSC.GT.10) GO TO 110
      L=NVSC
      NAMV(L)=NAM
176   NV(L,1)=VALUE(0)
      NV(L,2)=VALUE(0)
C
C.... CHECK FOR TYPE OF V-SOURCE
      IF (KARD(3).EQ.IHI.AND.KARD(4).EQ.IHN) GO TO 200
C.... CONVENTIONAL V-SOURCE
      KVOL(L)=1
      NSRC=NSRC+1
      NAM=0
      DO 180 I=1,9
      VVAL(I,L)=VALUE(1)
      TVLT(I,L)=VALUE(1)
      IF (NAM.NE.0) GO TO 185
180   CONTINUE
      I=10
185   NVV=I-1
      IF (NVV.GE.1) GO TO 190
      KVOL(L)=0
      NVV=1
190   NVVAL(L)=NVV
      TVLT(NVV+1,L)=1.E10
      VVAL(NVV+1,L)=VVAL(NVV,L)
      IF (NERR) 110,110,600
C.... SINUSOIDAL V-SOURCE
200   NSIN=NSIN+1
      KVOL(L)=2
      DO 210 I=1,4
 210  VVAL(I,L)=VALUE(1)
      IF (NERR) 110,110,600
C
C.... CURRENT SOURCE LINE
230   IF (NALTR.EQ.0) GO TO 234
      DO 232 L