FTN66,Y,J,E 
*************************************************************** 
* (C) COPYRIGHT HEWLETT-PACKARD COMPANY 1980.  ALL RIGHTS     * 
* RESERVED.  NO PART OF THIS PROGRAM MAY BE PHOTOCOPIED, RE-  * 
* REPRODUCED, OR TRANSLATED TO ANOTHER PROGRAM LANGUAGE WITH- * 
* OUT THE PRIOR WRITTEN CONSENT OF HEWLETT-PACKARD COMPANY.   * 
*************************************************************** 
*  NAME:   DCTR1                      * 
*  SOURCE: 92091-18005                * 
*  RELOC:  92091-16005                * 
*  PGMR:   C. GEBER                   * 
*************************************** 
$PAGE 
$EMA /BLANK/
$CDS ON 
      SUBROUTINE BJT
     + ,92091-16005 REV.5000 <890418.1140>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE PROCESSES BJTS FOR DC AND TRANSIENT ANALYSES.
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      DIMENSION VBEO(1),VBCO(1),CCO(1),CBO(1),GPIO(1),GMUO(1),GMO(1), 
     1   GOO(1),QBE(1),CQBE(1),QBC(1),CQBC(1),QCS(1),CQCS(1),QBX(1),
     2   CQBX(1),GXO(1),CEXBC(1),GEQCBO(1)
      EQUIVALENCE (VBEO(1),VALUE(1)),(VBCO(1),VALUE(2)),
     1   (CCO(1),VALUE(3)),(CBO(1),VALUE(4)),(GPIO(1),VALUE(5)),
     2   (GMUO(1),VALUE(6)),(GMO(1),VALUE(7)),(GOO(1),VALUE(8)),
     3   (QBE(1),VALUE(9)),(CQBE(1),VALUE(10)),(QBC(1),VALUE(11)),
     4   (CQBC(1),VALUE(12)),(QCS(1),VALUE(13)),(CQCS(1),VALUE(14)),
     5   (QBX(1),VALUE(15)),(CQBX(1),VALUE(16)),(GXO(1),VALUE(17)), 
     6   (CEXBC(1),VALUE(18)),(GEQCBO(1),VALUE(19)) 
$PAGE 
      LOC=LOCATE(12)
   10 IF (LOC.EQ.0) RETURN
      LOCV=NODPLC(LOC+1)
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NODE3=NODPLC(LOC+4) 
      NODE4=NODPLC(LOC+5) 
      NODE5=NODPLC(LOC+6) 
      NODE6=NODPLC(LOC+7) 
      NODE7=NODPLC(LOC+30)
      LOCM=NODPLC(LOC+8)
      IOFF=NODPLC(LOC+9)
      TYPE=NODPLC(LOCM+2) 
      LOCM=NODPLC(LOCM+1) 
      LOCT=NODPLC(LOC+22) 
      GCCS=0.0D0
      CEQCS=0.0D0 
      GEQBX=0.0D0 
      CEQBX=0.0D0 
      GEQCB=0.0D0 
C 
C  DC MODEL PARAMTERS 
C 
      AREA=VALUE(LOCV+1)
      BFM=VALUE(LOCM+2) 
      BRM=VALUE(LOCM+8) 
      CSAT=VALUE(LOCM+1)*AREA 
      RBPR=VALUE(LOCM+18)/AREA
      RBPI=VALUE(LOCM+16)/AREA-RBPR 
      GCPR=VALUE(LOCM+20)*AREA
      GEPR=VALUE(LOCM+19)*AREA
      OVA=VALUE(LOCM+4) 
      OVB=VALUE(LOCM+10)
      OIK=VALUE(LOCM+5)/AREA
      C2=VALUE(LOCM+6)*AREA 
      VTE=VALUE(LOCM+7)*VT
      OIKR=VALUE(LOCM+11)/AREA
      C4=VALUE(LOCM+12)*AREA
      VTC=VALUE(LOCM+13)*VT 
      VCRIT=VALUE(LOCM+54)
      TD=VALUE(LOCM+28) 
      XJRB=VALUE(LOCM+17)*AREA
$PAGE 
C 
C  INITIALIZATION 
C 
      ICHECK=1
      GO TO (100,20,30,50,60,70),INITF
   20 IF (MODE.NE.1.OR.MODEDC.NE.2.OR.NOSOLV.EQ.0) GO TO 25 
      VBE=TYPE*VALUE(LOCV+2)
      VCE=TYPE*VALUE(LOCV+3)
      VBC=VBE-VCE 
      VBX=VBC 
      VCS=0.0D0 
      GO TO 300 
   25 IF (IOFF.NE.0) GO TO 40 
      VBE=VCRIT 
      VBC=0.0D0 
      GO TO 300 
   30 IF (IOFF.EQ.0) GO TO 100
   40 VBE=0.0D0 
      VBC=0.0D0 
      GO TO 300 
   50 VBE=VBEO(LX0+LOCT)
      VBC=VBCO(LX0+LOCT)
      VBX=TYPE*(VALUE(LVNIM1+NODE2)-VALUE(LVNIM1+NODE4))
      VCS=TYPE*(VALUE(LVNIM1+NODE7)-VALUE(LVNIM1+NODE4))
      GO TO 300 
   60 VBE=VBEO(LX1+LOCT)
      VBC=VBCO(LX1+LOCT)
      VBX=TYPE*(VALUE(LVNIM1+NODE2)-VALUE(LVNIM1+NODE4))
      VCS=TYPE*(VALUE(LVNIM1+NODE7)-VALUE(LVNIM1+NODE4))
      IF (MODE.NE.2.OR.NOSOLV.EQ.0) GO TO 300 
      VBX=TYPE*(VALUE(LOCV+2)-VALUE(LOCV+3))
      VCS=0.0D0 
      GO TO 300 
   70 XFACT=DELTA/DELOLD(2) 
      VBEO(LX0+LOCT)=VBEO(LX1+LOCT) 
      VBE=(1.0D0+XFACT)*VBEO(LX1+LOCT)-XFACT*VBEO(LX2+LOCT) 
      VBCO(LX0+LOCT)=VBCO(LX1+LOCT) 
      VBC=(1.0D0+XFACT)*VBCO(LX1+LOCT)-XFACT*VBCO(LX2+LOCT) 
      CCO(LX0+LOCT)=CCO(LX1+LOCT) 
      CBO(LX0+LOCT)=CBO(LX1+LOCT) 
      GPIO(LX0+LOCT)=GPIO(LX1+LOCT) 
      GMUO(LX0+LOCT)=GMUO(LX1+LOCT) 
      GMO(LX0+LOCT)=GMO(LX1+LOCT) 
      GOO(LX0+LOCT)=GOO(LX1+LOCT) 
      GXO(LX0+LOCT)=GXO(LX1+LOCT) 
      GO TO 110 
$PAGE 
C 
C  COMPUTE NEW NONLINEAR BRANCH VOLTAGES
C 
  100 VBE=TYPE*(VALUE(LVNIM1+NODE5)-VALUE(LVNIM1+NODE6))
      VBC=TYPE*(VALUE(LVNIM1+NODE5)-VALUE(LVNIM1+NODE4))
  110 DELVBE=VBE-VBEO(LX0+LOCT) 
      DELVBC=VBC-VBCO(LX0+LOCT) 
      VBX=TYPE*(VALUE(LVNIM1+NODE2)-VALUE(LVNIM1+NODE4))
      VCS=TYPE*(VALUE(LVNIM1+NODE7)-VALUE(LVNIM1+NODE4))
      CCHAT=CCO(LX0+LOCT)+(GMO(LX0+LOCT)+GOO(LX0+LOCT))*DELVBE
     1   -(GOO(LX0+LOCT)+GMUO(LX0+LOCT))*DELVBC 
      CBHAT=CBO(LX0+LOCT)+GPIO(LX0+LOCT)*DELVBE+GMUO(LX0+LOCT)*DELVBC 
C 
C   BYPASS IF SOLUTION HAS NOT CHANGED
C 
      IF (INITF.EQ.6) GO TO 200 
      TOL=RELTOL*DMAX1(DABS(VBE),DABS(VBEO(LX0+LOCT)))+VNTOL
      IF (DABS(DELVBE).GE.TOL) GO TO 200
      TOL=RELTOL*DMAX1(DABS(VBC),DABS(VBCO(LX0+LOCT)))+VNTOL
      IF (DABS(DELVBC).GE.TOL) GO TO 200
      TOL=RELTOL*DMAX1(DABS(CCHAT),DABS(CCO(LX0+LOCT)))+ABSTOL
      IF (DABS(CCHAT-CCO(LX0+LOCT)).GE.TOL) GO TO 200 
      TOL=RELTOL*DMAX1(DABS(CBHAT),DABS(CBO(LX0+LOCT)))+ABSTOL
      IF (DABS(CBHAT-CBO(LX0+LOCT)).GE.TOL) GO TO 200 
      VBE=VBEO(LX0+LOCT)
      VBC=VBCO(LX0+LOCT)
      CC=CCO(LX0+LOCT)
      CB=CBO(LX0+LOCT)
      GPI=GPIO(LX0+LOCT)
      GMU=GMUO(LX0+LOCT)
      GM=GMO(LX0+LOCT)
      GO=GOO(LX0+LOCT)
      GX=GXO(LX0+LOCT)
      GEQCB=GEQCBO(LX0+LOCT)
      IF (MODE.NE.1) GO TO 800
      GO TO 900 
C 
C  LIMIT NONLINEAR BRANCH VOLTAGES
C 
  200 ICHK1=1 
      CALL PNJLIM(VBE,VBEO(LX0+LOCT),VT,VCRIT,ICHECK) 
      CALL PNJLIM(VBC,VBCO(LX0+LOCT),VT,VCRIT,ICHK1)
      IF (ICHK1.EQ.1) ICHECK=1
$PAGE 
C 
C  DETERMINE DC CURRENT AND DERIVITIVES 
C 
  300 VTN=VT*VALUE(LOCM+3)
      IF (VBE.LE.-5.0D0*VTN) GO TO 320
      EVBE=DEXP(VBE/VTN)
      CBE=CSAT*(EVBE-1.0D0)+GMIN*VBE
      GBE=CSAT*EVBE/VTN+GMIN
      IF (C2.NE.0.0D0) GO TO 310
      CBEN=0.0D0
      GBEN=0.0D0
      GO TO 350 
  310 EVBEN=DEXP(VBE/VTE) 
      CBEN=C2*(EVBEN-1.0D0) 
      GBEN=C2*EVBEN/VTE 
      GO TO 350 
  320 GBE=-CSAT/VBE+GMIN
      CBE=GBE*VBE 
      GBEN=-C2/VBE
      CBEN=GBEN*VBE 
  350 VTN=VT*VALUE(LOCM+9)
      IF (VBC.LE.-5.0D0*VTN) GO TO 370  
      EVBC=DEXP(VBC/VTN)
      CBC=CSAT*(EVBC-1.0D0)+GMIN*VBC
      GBC=CSAT*EVBC/VTN+GMIN
      IF (C4.NE.0.0D0) GO TO 360
      CBCN=0.0D0
      GBCN=0.0D0
      GO TO 400 
  360 EVBCN=DEXP(VBC/VTC) 
      CBCN=C4*(EVBCN-1.0D0) 
      GBCN=C4*EVBCN/VTC 
      GO TO 400 
  370 GBC=-CSAT/VBC+GMIN
      CBC=GBC*VBC 
      GBCN=-C4/VBC
      CBCN=GBCN*VBC 
C 
C  DETERMINE BASE CHARGE TERMS
C 
  400 Q1=1.0D0/(1.0D0-OVA*VBC-OVB*VBE)
      IF (OIK.NE.0.0D0) GO TO 405 
      IF (OIKR.NE.0.0D0) GO TO 405
      QB=Q1 
      DQBDVE=Q1*QB*OVB
      DQBDVC=Q1*QB*OVA
      GO TO 410 
  405 Q2=OIK*CBE+OIKR*CBC 
      ARG=DMAX1(0.0D0,1.0D0+4.0D0*Q2) 
      SQARG=1.0D0 
      IF (ARG.NE.0.0D0) SQARG=DSQRT(ARG)
      QB=Q1*(1.0D0+SQARG)/2.0D0 
      DQBDVE=Q1*(QB*OVB+OIK*GBE/SQARG)
      DQBDVC=Q1*(QB*OVA+OIKR*GBC/SQARG) 
$PAGE 
C 
C  WEIL'S APPROX. FOR EXCESS PHASE APPLIED WITH BACKWARD- 
C  EULER INTEGRATION
C 
  410 CC=0.0D0
      CEX=CBE 
      GEX=GBE 
      IF (MODE.EQ.1) GO TO 420
      IF (TD.EQ.0.0D0) GO TO 420
      ARG1=DELTA/TD 
      ARG2=3.0D0*ARG1 
      ARG1=ARG2*ARG1
      DENOM=1.0D0+ARG1+ARG2 
      ARG3=ARG1/DENOM 
      IF (INITF.NE.5) GO TO 411 
      CEXBC(LX1+LOCT)=CBE/QB
      CEXBC(LX2+LOCT)=CEXBC(LX1+LOCT) 
  411 CC=(CEXBC(LX1+LOCT)*(1.0D0+DELTA/DELOLD(2)+ARG2)
     1  -CEXBC(LX2+LOCT)*DELTA/DELOLD(2))/DENOM 
      CEX=CBE*ARG3
      GEX=GBE*ARG3
      CEXBC(LX0+LOCT)=CC+CEX/QB 
C 
C  DETERMINE DC INCREMENTAL CONDUCTANCES
C 
  420 CC=CC+(CEX-CBC)/QB-CBC/BRM-CBCN 
      CB=CBE/BFM+CBEN+CBC/BRM+CBCN
      GX=RBPR+RBPI/QB 
      IF (XJRB.EQ.0.0D0) GO TO 430
      ARG1=DMAX1(CB/XJRB,1.0D-9)
      ARG2=(-1.0D0+DSQRT(1.0D0+14.59025D0*ARG1))/2.4317D0/DSQRT(ARG1) 
      ARG1=DTAN(ARG2) 
      GX=RBPR+3.0D0*RBPI*(ARG1-ARG2)/ARG2/ARG1/ARG1 
  430 IF (GX.NE.0.0D0) GX=1.0D0/GX
      GPI=GBE/BFM+GBEN
      GMU=GBC/BRM+GBCN
      GO=(GBC+(CEX-CBC)*DQBDVC/QB)/QB 
      GM=(GEX-(CEX-CBC)*DQBDVE/QB)/QB-GO
      IF (MODE.NE.1) GO TO 500
      IF ((MODEDC.EQ.2).AND.(NOSOLV.NE.0)) GO TO 500
      IF (INITF.EQ.4) GO TO 500 
      GO TO 700 
$PAGE 
C 
C  CHARGE STORAGE ELEMENTS
C 
  500 TF=VALUE(LOCM+24) 
      TR=VALUE(LOCM+33) 
      CZBE=VALUE(LOCM+21)*AREA
      PE=VALUE(LOCM+22) 
      XME=VALUE(LOCM+23)
      CDIS=VALUE(LOCM+32) 
      CTOT=VALUE(LOCM+29)*AREA
      CZBC=CTOT*CDIS
      CZBX=CTOT-CZBC
      PC=VALUE(LOCM+30) 
      XMC=VALUE(LOCM+31)
      FCPE=VALUE(LOCM+46) 
      CZCS=VALUE(LOCM+38)*AREA
      PS=VALUE(LOCM+39) 
      XMS=VALUE(LOCM+40)
      XTF=VALUE(LOCM+25)
      OVTF=VALUE(LOCM+26) 
      XJTF=VALUE(LOCM+27)*AREA
      IF (TF.EQ.0.0D0) GO TO 505  
      IF (VBE.LE.0.0D0) GO TO 505 
      ARGTF=0.0D0 
      ARG2=0.0D0
      ARG3=0.0D0
      IF (XTF.EQ.0.0D0) GO TO 504 
      ARGTF=XTF 
      IF (OVTF.NE.0.0D0) ARGTF=ARGTF*DEXP(VBC*OVTF) 
      ARG2=ARGTF
      IF (XJTF.EQ.0.0D0) GO TO 503  
      TEMP=CBE/(CBE+XJTF) 
      ARGTF=ARGTF*TEMP*TEMP 
      ARG2=ARGTF*(3.0D0-TEMP-TEMP)
  503 ARG3=CBE*ARGTF*OVTF 
  504 CBE=CBE*(1.0D0+ARGTF)/QB
      GBE=(GBE*(1.0D0+ARG2)-CBE*DQBDVE)/QB
      GEQCB=TF*(ARG3-CBE*DQBDVC)/QB 
  505 IF (VBE.GE.FCPE) GO TO 510
      ARG=1.0D0-VBE/PE
      SARG=DEXP(-XME*DLOG(ARG)) 
      QBE(LX0+LOCT)=TF*CBE+PE*CZBE*(1.0D0-ARG*SARG)/(1.0D0-XME) 
      CAPBE=TF*GBE+CZBE*SARG
      GO TO 520 
  510 F1=VALUE(LOCM+47) 
      F2=VALUE(LOCM+48) 
      F3=VALUE(LOCM+49) 
      CZBEF2=CZBE/F2
      QBE(LX0+LOCT)=TF*CBE+CZBE*F1+CZBEF2*(F3*(VBE-FCPE)
     1   +(XME/(PE+PE))*(VBE*VBE-FCPE*FCPE))
      CAPBE=TF*GBE+CZBEF2*(F3+XME*VBE/PE) 
  520 FCPC=VALUE(LOCM+50) 
      F1=VALUE(LOCM+51) 
      F2=VALUE(LOCM+52) 
      F3=VALUE(LOCM+53) 
      IF (VBC.GE.FCPC) GO TO 530
      ARG=1.0D0-VBC/PC
      SARG=DEXP(-XMC*DLOG(ARG)) 
      QBC(LX0+LOCT)=TR*CBC+PC*CZBC*(1.0D0-ARG*SARG)/(1.0D0-XMC) 
      CAPBC=TR*GBC+CZBC*SARG
      GO TO 540 
  530 CZBCF2=CZBC/F2
      QBC(LX0+LOCT)=TR*CBC+CZBC*F1+CZBCF2*(F3*(VBC-FCPC)
     1   +(XMC/(PC+PC))*(VBC*VBC-FCPC*FCPC))
      CAPBC=TR*GBC+CZBCF2*(F3+XMC*VBC/PC) 
  540 IF (VBX.GE.FCPC) GO TO 550  
      ARG=1.0D0-VBX/PC
      SARG=DEXP(-XMC*DLOG(ARG)) 
      QBX(LX0+LOCT)=PC*CZBX*(1.0D0-ARG*SARG)/(1.0D0-XMC)
      CAPBX=CZBX*SARG 
      GO TO 560 
  550 CZBXF2=CZBX/F2
      QBX(LX0+LOCT)=CZBX*F1+CZBXF2*(F3*(VBX-FCPC)+(XMC/(PC+PC))*
     1   (VBX*VBX-FCPC*FCPC)) 
      CAPBX=CZBXF2*(F3+XMC*VBX/PC)
  560 IF (VCS.GE.0.0D0) GO TO 570 
      ARG=1.0D0-VCS/PS
      SARG=DEXP(-XMS*DLOG(ARG)) 
      QCS(LX0+LOCT)=PS*CZCS*(1.0D0-ARG*SARG)/(1.0D0-XMS)
      CAPCS=CZCS*SARG 
      GO TO 580 
  570 QCS(LX0+LOCT)=VCS*CZCS*(1.0D0+XMS*VCS/(2.0D0*PS)) 
      CAPCS=CZCS*(1.0D0+XMS*VCS/PS) 
C 
C  STORE SMALL-SIGNAL PARAMETERS
C 
  580 IF ((MODE.EQ.1).AND.(MODEDC.EQ.2).AND.(NOSOLV.NE.0)) GO TO 700
      IF (INITF.NE.4) GO TO 600 
      VALUE(LX0+LOCT+9)=CAPBE 
      VALUE(LX0+LOCT+11)=CAPBC
      VALUE(LX0+LOCT+13)=CAPCS
      VALUE(LX0+LOCT+15)=CAPBX
      VALUE(LX0+LOCT+17)=GEQCB
      GO TO 1000
$PAGE 
C 
C  TRANSIENT ANALYSIS 
C 
  600 IF (INITF.NE.5) GO TO 610 
      QBE(LX1+LOCT)=QBE(LX0+LOCT) 
      QBC(LX1+LOCT)=QBC(LX0+LOCT) 
      QBX(LX1+LOCT)=QBX(LX0+LOCT) 
      QCS(LX1+LOCT)=QCS(LX0+LOCT) 
  610 CALL INTGR8(GEQ,CEQ,CAPBE,LOCT+8) 
      GEQCB=GEQCB*AG(1) 
      GPI=GPI+GEQ 
      CB=CB+CQBE(LX0+LOCT)
      CALL INTGR8(GEQ,CEQ,CAPBC,LOCT+10)
      GMU=GMU+GEQ 
      CB=CB+CQBC(LX0+LOCT)
      CC=CC-CQBC(LX0+LOCT)
      IF (INITF.NE.5) GO TO 700 
      CQBE(LX1+LOCT)=CQBE(LX0+LOCT) 
      CQBC(LX1+LOCT)=CQBC(LX0+LOCT) 
C 
C  CHECK CONVERGENCE
C 
  700 IF (INITF.NE.3) GO TO 710 
      IF (IOFF.EQ.0) GO TO 710
      GO TO 750 
  710 IF (ICHECK.EQ.1) GO TO 720
      TOL=RELTOL*DMAX1(DABS(CCHAT),DABS(CC))+ABSTOL 
      IF (DABS(CCHAT-CC).GT.TOL) GO TO 720
      TOL=RELTOL*DMAX1(DABS(CBHAT),DABS(CB))+ABSTOL 
      IF (DABS(CBHAT-CB).LE.TOL) GO TO 750
  720 NONCON=NONCON+1 
  750 VBEO(LX0+LOCT)=VBE
      VBCO(LX0+LOCT)=VBC
      CCO(LX0+LOCT)=CC
      CBO(LX0+LOCT)=CB
      GPIO(LX0+LOCT)=GPI
      GMUO(LX0+LOCT)=GMU
      GMO(LX0+LOCT)=GM
      GOO(LX0+LOCT)=GO
      GXO(LX0+LOCT)=GX
      GEQCBO(LX0+LOCT)=GEQCB
      IF (MODE.EQ.1) GO TO 900
C 
C     CHARGE STORAGE FOR C-S AND B-X JUNCTIONS
C 
  800 CALL INTGR8(GCCS,CEQ,CAPCS,LOCT+12) 
      CEQCS=TYPE*(CQCS(LX0+LOCT)-VCS*GCCS)
      CALL INTGR8(GEQBX,CEQ,CAPBX,LOCT+14)
      CEQBX=TYPE*(CQBX(LX0+LOCT)-VBX*GEQBX) 
      IF (INITF.NE.5) GO TO 900 
      CQBX(LX1+LOCT)=CQBX(LX0+LOCT) 
      CQCS(LX1+LOCT)=CQCS(LX0+LOCT) 
$PAGE 
C 
C  LOAD CURRENT EXCITATION VECTOR 
C 
  900 CEQBE=TYPE*(CC+CB-VBE*(GM+GO+GPI)+VBC*(GO-GEQCB)) 
      CEQBC=TYPE*(-CC+VBE*(GM+GO)-VBC*(GMU+GO)) 
      VALUE(LVN+NODE2)=VALUE(LVN+NODE2)-CEQBX 
      VALUE(LVN+NODE4)=VALUE(LVN+NODE4)+CEQCS+CEQBX+CEQBC 
      VALUE(LVN+NODE5)=VALUE(LVN+NODE5)-CEQBE-CEQBC 
      VALUE(LVN+NODE6)=VALUE(LVN+NODE6)+CEQBE 
      VALUE(LVN+NODE7)=VALUE(LVN+NODE7)-CEQCS 
C 
C  LOAD Y MATRIX
C 
      LOCY=LVN+NODPLC(LOC+24) 
      VALUE(LOCY)=VALUE(LOCY)+GCPR
      LOCY=LVN+NODPLC(LOC+25) 
      VALUE(LOCY)=VALUE(LOCY)+GX+GEQBX
      LOCY=LVN+NODPLC(LOC+26) 
      VALUE(LOCY)=VALUE(LOCY)+GEPR
      LOCY=LVN+NODPLC(LOC+27) 
      VALUE(LOCY)=VALUE(LOCY)+GMU+GO+GCPR+GCCS+GEQBX
      LOCY=LVN+NODPLC(LOC+28) 
      VALUE(LOCY)=VALUE(LOCY)+GX  +GPI+GMU+GEQCB
      LOCY=LVN+NODPLC(LOC+29) 
      VALUE(LOCY)=VALUE(LOCY)+GPI+GEPR+GM+GO
      LOCY=LVN+NODPLC(LOC+10) 
      VALUE(LOCY)=VALUE(LOCY)-GCPR
      LOCY=LVN+NODPLC(LOC+11) 
      VALUE(LOCY)=VALUE(LOCY)-GX
      LOCY=LVN+NODPLC(LOC+12) 
      VALUE(LOCY)=VALUE(LOCY)-GEPR
      LOCY=LVN+NODPLC(LOC+13) 
      VALUE(LOCY)=VALUE(LOCY)-GCPR
      LOCY=LVN+NODPLC(LOC+14) 
      VALUE(LOCY)=VALUE(LOCY)-GMU+GM
      LOCY=LVN+NODPLC(LOC+15) 
      VALUE(LOCY)=VALUE(LOCY)-GM-GO 
      LOCY=LVN+NODPLC(LOC+16) 
      VALUE(LOCY)=VALUE(LOCY)-GX
      LOCY=LVN+NODPLC(LOC+17) 
      VALUE(LOCY)=VALUE(LOCY)-GMU-GEQCB 
      LOCY=LVN+NODPLC(LOC+18) 
      VALUE(LOCY)=VALUE(LOCY)-GPI 
      LOCY=LVN+NODPLC(LOC+19) 
      VALUE(LOCY)=VALUE(LOCY)-GEPR
      LOCY=LVN+NODPLC(LOC+20) 
      VALUE(LOCY)=VALUE(LOCY)-GO+GEQCB
      LOCY=LVN+NODPLC(LOC+21) 
      VALUE(LOCY)=VALUE(LOCY)-GPI-GM-GEQCB
      LOCY=LVN+NODPLC(LOC+31) 
      VALUE(LOCY)=VALUE(LOCY)+GCCS
      LOCY=LVN+NODPLC(LOC+32) 
      VALUE(LOCY)=VALUE(LOCY)-GCCS
      LOCY=LVN+NODPLC(LOC+33) 
      VALUE(LOCY)=VALUE(LOCY)-GCCS
      LOCY=LVN+NODPLC(LOC+34) 
      VALUE(LOCY)=VALUE(LOCY)-GEQBX 
      LOCY=LVN+NODPLC(LOC+35) 
      VALUE(LOCY)=VALUE(LOCY)-GEQBX 
 1000 LOC=NODPLC(LOC) 
      GO TO 10
      END 
$EMA /BLANK/                            
      SUBROUTINE CMEYER ( VGS0 , VGD0 , VGB0 , VON0 , VDSAT0 , VGS1 , 
     1    VGD1 , VGB1 , COVLGS , COVLGD , COVLGB , CGS0 , CGD0 , CGB0 , 
     2    CGS1 , CGD1 , CGB1 )
     + ,92091-16005 REV.5000 <890418.1140>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE COMPUTES THE MOSFET OVERLAP CAPACITANCES AS FUNCTIONS
C OF THE DEVICE TERMINAL VOLTAGES.
C 
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /MOSARG/ VTO,BETA,GAMMA,PHI,PHIB,COX,XNSUB,XNFS,XD,XJ,XLD, 
     1   XLAMDA,UO,UEXP,VBP,UTRA,VMAX,XNEFF,XL,XW,VBI,VON,VDSAT,QSPOF,
     2   BETA0,BETA1,CDRAIN,XQCO,XQC,FNARRW,FSHORT,LEV
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
C 
      INDAX=1 
      VGS=VGS1
      VGD=VGD1
      VGB=VGB1
      VONS=VON
      VBS=VGS-VGB 
      VDBSAT=VDSAT-VBS
      VDB=VGB-VGD 
   10 VDS=VGS-VGD 
      VGBT=VGS-VONS 
      IF (VGBT.GT.-PHI) GO TO 100 
      CGB=COX+COVLGB
      CGD=COVLGD
      CGS=COVLGS
      GO TO 430 
C 
  100 IF (VGBT.GT.-PHI/2.0D0) GO TO 200 
      CGB=-VGBT*COX/PHI+COVLGB
      CGD=COVLGD
      CGS=COVLGS
      GO TO 430 
C 
  200 IF (VGBT.GT.0.0D0) GO TO 300
      CGB=-VGBT*COX/PHI+COVLGB
      CGD=COVLGD
      CGS=COX/(7.5D-1*PHI)*VGBT+COX/1.5D0+COVLGS
      GO TO 430 
C 
  300 IF (VDBSAT.GT.VDB) GO TO 400
      CGB=COVLGB
      CGD=COVLGD
      CGS=COX/1.5D0+COVLGS
      GO TO 430 
$PAGE 
  400 VDDIF=2.0D0*VDBSAT-VDB
      VDDIF1=VDBSAT-VDB-1.0D-12 
      VDDIF2=VDDIF*VDDIF
      CGD=COX*(1.0D0-VDBSAT*VDBSAT/VDDIF2)/1.5D0+COVLGD 
      CGS=COX*(1.0D0-VDDIF1*VDDIF1/VDDIF2)/1.5D0+COVLGS 
      CGB=COVLGB
C 
  430 GO TO (440,560), INDAX
  440 INDAX=2 
      CGS1=CGS
      CGD1=CGD
      CGB1=CGB
      VGS=VGS0
      VGD=VGD0
      VGB=VGB0
      VONS=VON0 
      VBS=VGS-VGB 
      VDBSAT=VDSAT0-VBS 
      VDB=VGB-VGD 
      GO TO 10
C 
  560 CGS0=CGS
      CGD0=CGD
      CGB0=CGB
C 
C  FINISHED 
C 
 1000 RETURN
      END 
$EMA /BLANK/
      SUBROUTINE COMCOF 
     + ,92091-16005 REV.5000 <890418.1140>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE CALCULATES THE TIMESTEP-DEPENDENT TERMS USED IN THE
C NUMERICAL INTEGRATION.
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /FLAGS/ IPRNTA,IPRNTL,IPRNTM,IPRNTN,IPRNTO,LIMTIM,LIMPTS,
     1   LVLCOD,LVLTIM,ITL1,ITL2,ITL3,ITL4,ITL5,IGOOF,NOGO,KEOF 
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
      DIMENSION GMAT(7,7) 
C 
C  COMPUTE COEFFICIENTS FOR PARTICULAR INTEGRATION METHOD 
C 
      IF (METHOD.NE.1) GO TO 5
      IF (IORD.EQ.1) GO TO 5
C 
C...  TRAPEZOIDAL METHOD
C 
      AG(1)=1.0D0/DELTA/(1.0D0-XMU) 
      AG(2)=XMU/(1.0D0-XMU) 
      GO TO 200 
$PAGE 
C 
C  CONSTRUCT GEAR COEFFICIENT MATRIX
C 
    5 ISTOP=IORD+1
      CALL ZERO8(AG,ISTOP)
      AG(2)=-1.0D0
      DO 10 I=1,ISTOP 
      GMAT(1,I)=1.0D0 
   10 CONTINUE
      DO 20 I=2,ISTOP 
      GMAT(I,1)=0.0D0 
   20 CONTINUE
      ARG=0.0D0 
      DO 40 I=2,ISTOP 
      ARG=ARG+DELOLD(I-1) 
      ARG1=1.0D0
      DO 30 J=2,ISTOP 
      ARG1=ARG1*ARG 
      GMAT(J,I)=ARG1
   30 CONTINUE
   40 CONTINUE
C 
C  SOLVE FOR GEAR COEFFICIENTS AG(*)
C 
C 
C  LU DECOMPOSITION 
C 
      DO 70 I=2,ISTOP 
      JSTART=I+1
      IF (JSTART.GT.ISTOP) GO TO 70 
      DO 60 J=JSTART,ISTOP
      GMAT(J,I)=GMAT(J,I)/GMAT(I,I) 
      DO 50 K=JSTART,ISTOP
      GMAT(J,K)=GMAT(J,K)-GMAT(J,I)*GMAT(I,K) 
   50 CONTINUE
   60 CONTINUE
   70 CONTINUE
C 
C  FORWARD SUBSTITUTION 
C 
      DO 90 I=2,ISTOP 
      JSTART=I+1
      IF (JSTART.GT.ISTOP) GO TO 90 
      DO 80 J=JSTART,ISTOP
      AG(J)=AG(J)-GMAT(J,I)*AG(I) 
   80 CONTINUE
   90 CONTINUE
$PAGE 
C 
C  BACKWARD SUBSTITUTION
C 
      AG(ISTOP)=AG(ISTOP)/GMAT(ISTOP,ISTOP) 
      IR=ISTOP
      DO 110 I=2,ISTOP
      JSTART=IR 
      IR=IR-1 
      DO 100 J=JSTART,ISTOP 
      AG(IR)=AG(IR)-GMAT(IR,J)*AG(J)
  100 CONTINUE
      AG(IR)=AG(IR)/GMAT(IR,IR) 
  110 CONTINUE
C 
C  FINISHED 
C 
  200 RETURN
      END 
$EMA /BLANK/                            
      SUBROUTINE DCDCMP 
     + ,92091-16005 REV.5000 <890418.1140>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE SWAPS ROWS AND COLUMNS IN THE COEFFICIENT MATRIX ACCOR-
C DING TO THE NUMERICAL PIVOTING AND MINIMUM FILLIN TERMS.IT THEN PERFORMS
C AN IN-PLACE LU FACTORIZATION OF THE COEFFICIENT MATRIX. 
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /MISCEL/ ATIME,APROG(3),ADATE,ATITLE(10),DEFL,DEFW,DEFAD,
     1  DEFAS,RSTATS(50),IWIDTH,LWIDTH,NOPAGE 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /FLAGS/ IPRNTA,IPRNTL,IPRNTM,IPRNTN,IPRNTO,LIMTIM,LIMPTS,
     1   LVLCOD,LVLTIM,ITL1,ITL2,ITL3,ITL4,ITL5,IGOOF,NOGO,KEOF 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      COMMON/DEBUG/ IDEBUG(20)
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
$PAGE 
      CALL SECOND(T1) 
      IF (IPIV.LE.0) GO TO 12 
      IF (IDEBUG(11).LE.0) GO TO  3 
      CALL DMPMAT(6HDCDCMP) 
      IDEBUG(11)=IDEBUG(11)-1 
    3 DO 10 I=2,NSTOP 
      NO=0
      LOC=NODPLC(JCPT+I)
    5 IF (LOC.EQ.0) GO TO 7 
      NO=NO+1 
      LOC=NODPLC(JCPT+LOC)
      GO TO 5 
    7 NODPLC(NUMOFF+I)=NO 
      NO=0
      LOC=NODPLC(IRPT+I)
    8 IF (LOC.EQ.0) GO TO 9 
      NO=NO+1 
      LOC=NODPLC(IRPT+LOC)
      GO TO 8 
    9 NODPLC(NMOFFC+I)=NO 
   10 CONTINUE
   12 N=1 
C 
C     FIND NEXT PIVOT 
C 
   15 N=N+1 
      IF (IPIV.GT.0) GO TO 20 
      IF (IDEBUG(13).LE.0) GO TO 17 
      CALL DMPMAT(6HDCDCM2) 
      IDEBUG(13)=IDEBUG(13)-1 
   17 IF (IDEBUG(14).LE.0) GO TO 18 
      IF (MODE.NE.2) GO TO 18 
      CALL DMPMAT(6HDCDCM3) 
      IDEBUG(14)=IDEBUG(14)-1 
   18 IF (IDEBUG(15).LE.0.OR.ICALC.LE.10) GO TO 19
      CALL DMPMAT(6HDCDCM4) 
      IDEBUG(15)=IDEBUG(15)-1 
   19 NXTI=N
      NXTJ=N
      GO TO 120 
$PAGE 
C 
C     SEARCH THE CORESPONDING COLUMN FOR MAX ENTRY
C 
   20 VMAX=0.0D0
      LOCI=N
   25 LOCI=NODPLC(IRPT+LOCI)
      IF (LOCI.EQ.0) GO TO 50 
      I=NODPLC(IROWNO+LOCI) 
      IF (I.LT.N) GO TO 25
   30 IF (DABS(VALUE(LVN+LOCI)).LE.VMAX) GO TO 25 
      VMAX=DABS(VALUE(LVN+LOCI))
      GO TO 25
   50 IF (VMAX.GT.PIVTOL) GO TO 60
      WRITE(6,51) N,VMAX
   51 FORMAT('0*ERROR*:  MAXIMUM ENTRY IN THIS COLUMN AT STEP ',I4,' (',
     1   1PD12.6,') IS LESS THAN PIVTOL') 
      IGOOF=1 
      RETURN
   60 EPSREL=DMAX1(PIVREL*VMAX,PIVTOL)
      IF (N.GE.NSTOP) GO TO 200 
      IF (IPIV.LE.0) GO TO 120
C 
C     PIVOTING ON THE DIAGONAL
C 
      MINOP=100000
      NXTI=0
      DO 70 I=N,NSTOP 
      I1=NODPLC(IRSWPF+I) 
      J1=NODPLC(ICSWPF+I) 
      ISPOT=INDXX(I1,J1)
      IF (ISPOT.EQ.1) GO TO 70
      IF (DABS(VALUE(LVN+ISPOT)).LT.EPSREL) GO TO 70
      NOP=(NODPLC(NUMOFF+I)-1)*(NODPLC(NMOFFC+I)-1) 
      IF (NOP.GE.MINOP) GO TO 70
      MINOP=NOP 
      NXTI=I
      NXTJ=I
      IF (MINOP.LE.0) GO TO 95
   70 CONTINUE
      IF (NXTI.LE.0) GO TO 75 
      IF (NXTI-N) 120,120,100 
$PAGE 
C 
C     PIVOTING ON THE ENTIRE MATRIX 
C 
   75 DO 90 I=N,NSTOP 
      LOC=I 
   80 LOC=NODPLC(JCPT+LOC)
      IF (LOC.EQ.0) GO TO 90
      J=NODPLC(JCOLNO+LOC)
      IF (J.LT.N) GO TO 80
      IF (DABS(VALUE(LVN+LOC)).LT.EPSREL) GO TO 80
      NOP=(NODPLC(NUMOFF+I)-1)*(NODPLC(NMOFFC+J)-1) 
      IF (NOP.GE.MINOP) GO TO 80
      MINOP=NOP 
      NXTI=I
      NXTJ=J
      IF (MINOP.LE.0) GO TO 95
   90 CONTINUE
      IF (NXTI.GT.0) GO TO 95 
      WRITE (6,92)
   92 FORMAT('0*ABORT*:  PIVOT NOT IN DCDCMP')
      IGOOF=1 
      GO TO 200 
   95 IF (NXTI.EQ.N.AND.NXTJ.EQ.N) GO TO 120
      IF (NXTI.EQ.N) GO TO 105
C 
C     A PIVOT HAS BEEN FOUND
C 
  100 LOAD=NODPLC(IRSWPF+NXTI)
      LR=NODPLC(IRSWPF+N) 
      NODPLC(IRSWPF+NXTI)=LR
      NODPLC(IRSWPR+LR)=NXTI
      NODPLC(IRSWPF+N)=LOAD 
      NODPLC(IRSWPR+LOAD)=N 
      NOFF=NODPLC(NUMOFF+NXTI)
      NODPLC(NUMOFF+NXTI)=NODPLC(NUMOFF+N)
      NODPLC(NUMOFF+N)=NOFF 
      IF (NXTJ.EQ.N) GO TO 110
  105 LOAD=NODPLC(ICSWPF+NXTJ)
      LC=NODPLC(ICSWPF+N) 
      NODPLC(ICSWPF+NXTJ)=LC
      NODPLC(ICSWPR+LC)=NXTJ
      NODPLC(ICSWPF+N)=LOAD 
      NODPLC(ICSWPR+LOAD)=N 
      NOFF=NODPLC(NMOFFC+NXTJ)
      NODPLC(NMOFFC+NXTJ)=NODPLC(NMOFFC+N)
      NODPLC(NMOFFC+N)=NOFF 
  110 CALL SWAPIJ(NXTI,N,NXTJ,N)
      NXTI=N
      NXTJ=N
$PAGE 
C 
C     CALCULATE CONTRIBUTION FROM NXTI, NXTJ AND FIND FILL-INS
C 
  120 IF (N.GE.NSTOP) GO TO 200 
      N1=NODPLC(IRSWPF+NXTI)
      N2=NODPLC(ICSWPF+NXTJ)
      LOCNN=INDXX(N1,N2)
      IF (IPIV.LE.0 .AND. DABS(VALUE(LVN+LOCNN)).LT.PIVTOL) GO TO 220 
C 
C     DOWN COL J
C 
      LOCR=NODPLC(IRPT+LOCNN) 
  125 IF (LOCR.EQ.0) GO TO 180
      I=NODPLC(IROWNO+LOCR) 
      VALUE(LVN+LOCR)=VALUE(LVN+LOCR)/VALUE(LVN+LOCNN)
      LOCC=NODPLC(JCPT+LOCNN) 
C 
C     FOR EACH COLUMN ELEMENT LOOK UP ROW NXTI
C 
  130 IF (LOCC.EQ.0) GO TO 170
      J=NODPLC(JCOLNO+LOCC) 
C 
C     CHECK FOR FILL-IN (I,J) 
C 
      IF (IPIV.LE.0) GO TO 135
      CALL SIZMEM(JCPT,ISIZE1)
      CALL RESERV(I,J)
      CALL SIZMEM(JCPT,ISIZE2)
      IF (ISIZE1.EQ.ISIZE2) GO TO 135 
      CALL EXTMEM(LVN,1)
      NTTBR=NTTBR+1 
      VALUE(LVN+NSTOP+NTTBR)=0.0D0
C 
C     LOCATE ELEMENT (I,J)
C 
  135 IF (J.LT.I) GO TO 145 
      LOCIJ=LOCC
  140 LOCIJ=NODPLC(IRPT+LOCIJ)
      IF (NODPLC(IROWNO+LOCIJ).EQ.I) GO TO 155
      GO TO 140 
  145 LOCIJ=LOCR
  150 LOCIJ=NODPLC(JCPT+LOCIJ)
      IF (NODPLC(JCOLNO+LOCIJ).EQ.J) GO TO 155
      GO TO 150 
  155 VALUE(LVN+LOCIJ)=VALUE(LVN+LOCIJ)-
     1                  VALUE(LVN+LOCC)*VALUE(LVN+LOCR) 
  160 LOCC=NODPLC(JCPT+LOCC)
      GO TO 130 
  170 LOCR=NODPLC(IRPT+LOCR)
      IF (IPIV.LE.0) GO TO 125
      NODPLC(NUMOFF+I)=NODPLC(NUMOFF+I)-1 
      GO TO 125 
$PAGE 
C 
C     REDUCE NMOFFC FOR EACH ELEMENT IN COL NXTI
C 
  180 IF (IPIV.LE.0) GO TO 15 
      LOCC=NODPLC(JCPT+LOCNN) 
  185 IF (LOCC.EQ.0) GO TO 15 
      J=NODPLC(JCOLNO+LOCC) 
      NODPLC(NMOFFC+J)=NODPLC(NMOFFC+J)-1 
      LOCC=NODPLC(JCPT+LOCC)
      GO TO 185 
C 
C     DONE
C 
  200 IF (IPIV.EQ.0) GO TO 210
      IF (IDEBUG(17).LE.0) GO TO 202
      CALL DMPMAT(6HDCDCM5) 
      IDEBUG(17)=IDEBUG(17)-1 
  202 CALL MATLOC 
      RSTATS(49)=RSTATS(49)+1.0D0 
      IPIV=0
      IF (LVLCOD.EQ.2) LVLCOD=3 
      IFILL=NTTBR-NTTAR 
      PERSPA=100.0D0*(1.0D0-DFLOAT(NTTBR)/DFLOAT(NSTOP*NSTOP))
$PAGE 
C 
C  CALCULATION OF OPERATION COUNT (OPERATION := `*' or `/'):
C 
C     NOFFR := OFF-DIAGONAL ELEMENTS IN ROW, NOT INCLUDING DIAGONAL,
C                COUNTING ONLY THOSE ELEMENTS IN THE REMAINDER MATRIX 
C     NOFFC := OFF-DIAGONAL ELEMENTS IN COLUMN, NOT INCLUDING DIAGONAL, 
C                COUNTING ONLY THOSE ELEMENTS IN THE REMAINDER MATRIX 
C 
C     THEN WE HAVE
C 
C        LU DECOMPOSITION     REQUIRES SIGMA(2,NSTOP-1) {NOFFC + NOFFC*NOFFR} 
C        FORWARD SUBSTITUTION          SIGMA(2,NSTOP-1) {NOFFC + 1}   +   1 
C        BACKWARD SUBSTITUTION         SIGMA(2,NSTOP-1) {NOFFR} 
C 
C     WHICH SUMS TO 
C 
C               SIGMA(2,NSTOP-1) {NOFFC + NOFFC*NOFFR + (NOFFC+1) + NOFFR} + 1
C         OR
C               SIGMA(2,NSTOP-1) {NOFFC*(NOFFR+2) + NOFFR + 1}   +   1
C 
      IOPS=1
      NSTOP1=NSTOP-1
      DO 205 I=2,NSTOP1 
      NOFFR=NODPLC(NUMOFF+I)-1
      NOFFC=NODPLC(NMOFFC+I)-1
      IOPS=IOPS+NOFFR+NOFFC*(NOFFR+2)+1 
  205 CONTINUE
      RSTATS(20)=NSTOP
      RSTATS(21)=NTTAR
      RSTATS(22)=NTTBR
      RSTATS(23)=IFILL
      RSTATS(24)=0.0D0
      RSTATS(25)=NTTBR
      RSTATS(26)=IOPS 
      RSTATS(27)=PERSPA 
      GO TO 215 
  210 IF (IDEBUG(18).LE.0) GO TO 212
      CALL DMPMAT(6HDCDCM6) 
      IDEBUG(18)=IDEBUG(18)-1 
  212 IF (IDEBUG(19).LE.0.OR.ICALC.LE.10) GO TO 215 
      CALL DMPMAT(6HDCDCM7) 
      IDEBUG(19)=IDEBUG(19)-1 
  215 CALL SECOND(T2) 
      RSTATS(50)=RSTATS(50)+T2-T1 
       RETURN 
  220 IPIV=1
      WRITE(6,221) N,NXTI,NXTJ,ITERNO,TIME
  221 FORMAT(' PIVOT CHANGE ON FLY: N= ',I5,' NXTI= ',I5,' NXTJ= ', 
     1       I5,' ITERNO= ',I5,' TIME= ',1PD12.5) 
      RSTATS(49)=RSTATS(49)+1.0D0 
      GO TO 20
      END 
$EMA /BLANK/                            
      SUBROUTINE DCSOL
     + ,92091-16005 REV.5000 <890418.1140>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE SOLVES THE SYSTEM OF CIRCUIT EQUATIONS BY PERFORMING 
C A FORWARD AND BACKWARD SUBSTITUTION STEP USING THE PREVIOUSLY-COMPUTED
C LU FACTORS. 
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /MISCEL/ ATIME,APROG(3),ADATE,ATITLE(10),DEFL,DEFW,DEFAD,
     1  DEFAS,RSTATS(50),IWIDTH,LWIDTH,NOPAGE 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
C  FORWARD SUBSTITUTION 
C 
      CALL SECOND(T1) 
      DO 20 I=2,NSTOP 
      LOC=I 
      IORD=NODPLC(IRSWPF+I) 
   10 LOC=NODPLC(JCPT+LOC)
      IF (NODPLC(JCOLNO+LOC).GE.I) GO TO 20 
      J=NODPLC(JCOLNO+LOC)
      JORD=NODPLC(IRSWPF+J) 
      VALUE(LVN+IORD)=VALUE(LVN+IORD)-
     1           VALUE(LVN+LOC)*VALUE(LVN+JORD) 
      GO TO 10
   20 CONTINUE
$PAGE 
C 
C     BACK SUBSTITUTION 
C 
      I=NSTOP 
      IORD=NODPLC(IRSWPF+I) 
      JORD=NODPLC(ICSWPF+I) 
      LOCNN=INDXX(IORD,JORD)
   30 VALUE(LVN+IORD)=VALUE(LVN+IORD)/VALUE(LVN+LOCNN)
      I=I-1 
      IF (I.LE.1) GO TO 100 
      IORD=NODPLC(IRSWPF+I) 
      LOC=I 
   35 LOC=NODPLC(JCPT+LOC)
   40 IF (NODPLC(JCOLNO+LOC).NE.I) GO TO 35 
      LOCNN=LOC 
   50 LOC=NODPLC(JCPT+LOC)
      IF (LOC.EQ.0) GO TO 30
      J=NODPLC(JCOLNO+LOC)
      JORD=NODPLC(IRSWPF+J) 
      VALUE(LVN+IORD)=VALUE(LVN+IORD)-
     1           VALUE(LVN+LOC)*VALUE(LVN+JORD) 
      GO TO 50
  100 CALL SECOND(T2) 
      RSTATS(46)=RSTATS(46)+T2-T1 
      RETURN
      END 
$EMA /BLANK/                            
$ALIAS MOSFET='MOSFET_LINK' 
      SUBROUTINE DCTRAN 
     + ,92091-16005 REV.5000 <890418.1140>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C 
C     THIS ROUTINE CONTROLS THE DC TRANSFER CURVE, DC OPERATING POINT,
C AND * ANALYSES.  THE VARIABLES MODE AND MODEDC (DEFINED BELOW)
C DETERMINE EXACTLY WHICH ANALYSIS IS PERFORMED.
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /MISCEL/ ATIME,APROG(3),ADATE,ATITLE(10),DEFL,DEFW,DEFAD,
     1  DEFAS,RSTATS(50),IWIDTH,LWIDTH,NOPAGE 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /FLAGS/ IPRNTA,IPRNTL,IPRNTM,IPRNTN,IPRNTO,LIMTIM,LIMPTS,
     1   LVLCOD,LVLTIM,ITL1,ITL2,ITL3,ITL4,ITL5,IGOOF,NOGO,KEOF 
      COMMON /DC/ TCSTAR(2),TCSTOP(2),TCINCR(2),ICVFLG,ITCELM(2),KSSOP, 
     1   KINEL,KIDIN,KOVAR,KIDOUT 
      COMMON /TRAN/ TSTEP,TSTOP,TSTART,DELMAX,TDMAX,FORFRE,JTRFLG 
      COMMON /CJE/ MAXTIM,ITIME,ICOST 
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
      LOGICAL MEMPTR
C 
      DIMENSION SUBTIT(4,2) 
      DIMENSION AVHDR(3),AVFRM(4) 
      DATA AVHDR / 8H( (2X,A4, 8H,3X,A7,3, 5HX)//) /
      DATA AVFRM / 8H( (1H ,A, 8H1,I3,1H), 8H,F10.4,3, 4HX)/) / 
      DATA ANODE, AVLTG / 4HNODE, 7HVOLTAGE / 
      DATA SUBTIT / 8HSMALL SI, 8HGNAL BIA, 8HS SOLUTI, 8HON      , 
     1              8HINITIAL , 8HTRANSIEN, 8HT SOLUTI, 8HON      / 
      DATA LPRN /1H(/ 
$PAGE 
C 
C      THE VARIABLES *MODE*, *MODEDC*, AND *INITF* ARE USED BY SPICE TO 
C KEEP TRACK OF THE STATE OF THE ANALYSIS.  THE VALUES OF THESE FLAGS 
C (AND THE CORRESPONDING MEANINGS) ARE AS FOLLOWS:
C 
C        FLAG    VALUE    MEANING 
C        ----    -----    ------- 
C 
C        MODE      1      DC ANALYSIS (SUBTYPE DEFINED BY *MODEDC*) 
C                  2      TRANSIENT ANALYSIS
C                  3      AC ANALYSIS (SMALL SIGNAL)
C 
C        MODEDC    1      DC OPERATING POINT
C                  2      INITIAL OPERATING POINT FOR TRANSIENT ANALYSIS
C                  3      DC TRANSFER CURVE COMPUTATION 
C 
C        INITF     1      CONVERGE WITH 'OFF' DEVICES ALLOWED TO FLOAT
C                  2      INITIALIZE JUNCTION VOLTAGES
C                  3      CONVERGE WITH 'OFF' DEVICES HELD 'OFF'
C                  4      STORE SMALL-SIGNAL PARAMETERS AWAY
C                  5      FIRST TIMEPOINT IN TRANSIENT ANALYSIS 
C                  6      PREDICTION STEP 
C 
C NOTE:  *MODEDC* IS ONLY SIGNIFICANT IF *MODE* = 1.
C 
C 
C  INITIALIZE 
C 
      CALL SECOND(T1) 
C 
C.. DON'T TAKE ANY CHANCES WITH LX3, SET TO LARGE NUMBER
C 
      LX3=20000000
      LX2=20000000
C 
C.. SEE IF LX3 AND LX2 TABLES ARE NEEDED
C 
      NOLX2=0 
      NOLX3=0 
   20 LOCTIM=5
C 
C.. POST-PROCESSING INITIALIZATION
C 
      IF (IPOSTP.EQ.0) GO TO 25 
      NUMCUR=JELCNT(9)
      NUMPOS=NUNODS+NUMCUR
      CALL GETM8(IBUFF,NUMPOS)
      NUMPOS=NUMPOS*4 
      IF (NUMCUR.EQ.0) GO TO 25 
      LOC=LOCATE(9) 
      LOCCUR=NODPLC(LOC+6)-1
$PAGE 
C 
C...  SET UP FORMAT 
C 
   25 NVPRLN=4+(LWIDTH-72)/19 
      NVPRLN=MIN0(NVPRLN,NCNODS-1)
      IPOS=2
      CALL ALFNUM(NVPRLN,AVFRM,IPOS)
      IPOS=2
      CALL ALFNUM(NVPRLN,AVHDR,IPOS)
C 
C...  ALLOCATE STORAGE
C 
      IF (MODE.EQ.2) GO TO 35 
      NEED=4*NSTOP+NTTBR+NXTRM
      CALL AVLM8(NAVL)
      IF (NEED.LE.NAVL) GO TO 30
C 
C...  NOT ENOUGH MEMORY FOR DC OPERATING POINT ANALYSIS 
C 
      WRITE(6,26) NEED,NAVL 
   26 FORMAT('0INSUFFICIENT MEMORY AVAILABLE FOR DC ANALYSIS.',/
     1' MEMORY REQUIRED ',I6,', MEMORY AVAILABLE ',I6,'.')
      NOGO=1
      GO TO 1100
   30 CALL GETM8(LVNIM1,NSTOP)
      CALL GETM8(LVN,NSTOP+NTTBR) 
      CALL SLPMEM(LVN,NSTOP)
      CALL GETM8(LX0,NXTRM) 
      CALL GETM8(LVNTMP,NSTOP)
      IF (MODEDC.NE.3) GO TO 45 
   35 CALL GETM8(LX1,NXTRM) 
      IF (NOLX2.EQ.0) CALL GETM8(LX2,NXTRM) 
      IF (MODE.NE.2) GO TO 40 
      IF (NOLX3.EQ.0) CALL GETM8(LX3,NXTRM) 
      CALL GETM8(LTD,0) 
   40 CALL GETM8(LOUTPT,0)
   45 CALL CRUNCH 
   50 IF (MODE.EQ.2) GO TO 500
      TIME=0.0D0
      AG(1)=0.0D0 
      CALL SORUPD 
      IF (MODEDC.EQ.3) GO TO 300
$PAGE 
C 
C  ....  SINGLE POINT DC ANALYSIS 
C 
C 
C  COMPUTE DC OPERATING POINT 
C 
  100 INITF=2 
      CALL ITER8(ITL1)
      RSTATS(6)=RSTATS(6)+ITERNO
      IF (IGOOF.NE.0) GO TO 150 
      IF (MODEDC.NE.1) GO TO 120
      INITF=4 
      CALL DIODE
      CALL BJT
      CALL JFET 
      CALL MOSFET 
C 
C  PRINT OPERATING POINT
C 
  120 IF ((MODE.EQ.1).AND.(MODEDC.EQ.2).AND.(NOSOLV.NE.0)) GO TO 1000 
      CALL TITLE(-1,LWIDTH,1,SUBTIT(1,MODEDC))
      WRITE (6,AVHDR) (ANODE,AVLTG,I=1,NVPRLN)
      WRITE (6,AVFRM) (LPRN,NODPLC(JUNODE+I),VALUE(LVNIM1+I),I=2,NCNODS)
      GO TO 1000
C 
C  NO CONVERGENCE 
C 
  150 NOGO=1
      WRITE (6,151) 
  151 FORMAT('1*ERROR*:  NO CONVERGENCE IN DC ANALYSIS'/'0LAST NODE VOL'
     1   ,'TAGES:'/)
      WRITE (6,AVHDR) (ANODE,AVLTG,I=1,NVPRLN)
      WRITE (6,AVFRM) (LPRN,NODPLC(JUNODE+I),VALUE(LVNIM1+I),I=2,NCNODS)
      GO TO 1000
$PAGE 
C 
C  ....  DC TRANSFER CURVES 
C 
  300 NUMOUT=JELCNT(41)+1 
      ITEMP=ITCELM(1) 
      LOCS=NODPLC(ITEMP+1)
      TEMVAL=VALUE(LOCS+1)
      ICVFL2=1
      IF (ITCELM(2).EQ.0) GO TO 310 
      ITEMP=ITCELM(2) 
      LOCS2=NODPLC(ITEMP+1) 
      TEMV2=VALUE(LOCS2+1)
      VALUE(LOCS2+1)=TCSTAR(2)
      TEMP=DABS((TCSTOP(2)-TCSTAR(2))/TCINCR(2))+0.5D0
      ICVFL2=IDINT(TEMP)+1
      ICVFL2=MAX0(ICVFL2,1) 
  310 DELTA=TCINCR(1) 
      DO 320 I=1,7
      DELOLD(I)=DELTA 
  320 CONTINUE
      ICVFL1=ICVFLG/ICVFL2
      VALUE(LOCS+1)=TCSTAR(1) 
      ICALC=0 
      ICAL2=0 
      LOCTIM=3
  340 INITF=2 
      CALL ITER8(ITL1)
      RSTATS(4)=RSTATS(4)+ITERNO
      CALL COPY8(VALUE(LX0+1),VALUE(LX1+1),NXTRM) 
      IF (NOLX2.EQ.0) CALL COPY8(VALUE(LX0+1),VALUE(LX2+1),NXTRM) 
      IF (IGOOF.NE.0) GO TO 450 
      GO TO 360 
  350 CALL GETCJE 
      IF ((MAXTIM-ITIME).LE.LIMTIM) GO TO 460 
      INITF=6 
      CALL ITER8(ITL2)
      RSTATS(4)=RSTATS(4)+ITERNO
      IF (IGOOF.NE.0) GO TO 340 
$PAGE 
C 
C  STORE OUTPUTS
C 
  360 CALL EXTMEM(LOUTPT,NUMOUT)
      LOCO=LOUTPT+ICALC*NUMOUT
      ICALC=ICALC+1 
      ICAL2=ICAL2+1 
      VALUE(LOCO+1)=VALUE(LOCS+1) 
      LOC=LOCATE(41)
  370 IF (LOC.EQ.0) GO TO 400 
      IF (NODPLC(LOC+5).NE.0) GO TO 380 
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      ISEQ=NODPLC(LOC+4)
      VALUE(LOCO+ISEQ)=VALUE(LVNIM1+NODE1)-VALUE(LVNIM1+NODE2)
      LOC=NODPLC(LOC) 
      GO TO 370 
  380 IPTR=NODPLC(LOC+2)
      IPTR=NODPLC(IPTR+6) 
      ISEQ=NODPLC(LOC+4)
      VALUE(LOCO+ISEQ)=VALUE(LVNIM1+IPTR) 
      LOC=NODPLC(LOC) 
      GO TO 370 
C 
C  INCREMENT SOURCE VALUE 
C 
  400 IF (IPOSTP.EQ.0) GO TO 410
      VALUE(IBUFF+1)=VALUE(LOCS+1)
      CALL COPY8(VALUE(LVNIM1+2),VALUE(IBUFF+2),NUNODS-1) 
      IF (NUMCUR.NE.0) CALL COPY8(VALUE(LVNIM1+LOCCUR+1), 
     1  VALUE(IBUFF+NUNODS+1),NUMCUR) 
  410 IF (ICALC.GE.ICVFLG) GO TO 490
      IF (ICAL2.GE.ICVFL1) GO TO 480
      IF (NOLX2.NE.0) GO TO 420 
      CALL PTRMEM(LX2,ITEMP)
      CALL PTRMEM(LX1,LX2)
      GO TO 430 
  420 CALL PTRMEM(LX1,ITEMP)
  430 CALL PTRMEM(LX0,LX1)
      CALL PTRMEM(ITEMP,LX0)
      VALUE(LOCS+1)=TCSTAR(1)+DFLOAT(ICAL2)*DELTA 
      GO TO 350 
$PAGE 
C 
C  NO CONVERGENCE 
C 
  450 ITEMP=ITCELM(1) 
      LOCE=NODPLC(ITEMP+1)
      WRITE (6,451) VALUE(LOCE),VALUE(LOCS+1) 
  451 FORMAT('1*ERROR*:  NO CONVERGENCE IN DC TRANSFER CURVES AT ',A8,
     1   ' = ',1PD10.3/'0LAST NODE VOLTAGES:'/) 
      WRITE (6,AVHDR) (ANODE,AVLTG,I=1,NVPRLN)
      WRITE (6,AVFRM) (LPRN,NODPLC(JUNODE+I),VALUE(LVNIM1+I),I=2,NCNODS)
      GO TO 470 
  460 WRITE (6,461) 
  461 FORMAT('0*ERROR*:  CPU TIME LIMIT EXCEEDED ... ANALYSIS STOPPED'/)
  470 NOGO=1
      GO TO 490 
C 
C... RESET FIRST SWEEP VARIABLE ... STEP SECOND 
C 
  480 ICAL2=0 
      VALUE(LOCS+1)=TCSTAR(1) 
      VALUE(LOCS2+1)=VALUE(LOCS2+1)+TCINCR(2) 
      GO TO 340 
C 
C  FINISHED WITH DC TRANSFER CURVES 
C 
  490 VALUE(LOCS+1)=TEMVAL
      IF (ITCELM(2).NE.0) VALUE(LOCS2+1)=TEMV2
      IF (IPOSTP.EQ.0) GO TO 1000 
      GO TO 1000
$PAGE 
C 
C  ....  TRANSIENT ANALYSIS 
C 
  500 NUMOUT=JELCNT(42)+1 
C 
C...  LIMIT DELMAX IF NO ENERGY-STORAGE ELEMENTS
C 
      NUMESE=JELCNT(2)+JELCNT(3)+JELCNT(11)+JELCNT(12)+JELCNT(13) 
     1   +JELCNT(14)
      IF (NUMESE.EQ.0) DELMAX=DMIN1(DELMAX,TSTEP) 
      INITF=5 
      IORD=1
      LOCTIM=9
      ICALC=0 
      NUMTP=0 
      NUMRTP=0
      NUMNIT=0
      TIME=0.0D0
      IBKFLG=1
      DELBKP=DELMAX 
      NBKPT=1 
      DELTA=DELMAX
      DO 510 I=1,7
      DELOLD(I)=DELTA 
  510 CONTINUE
      DELNEW=DELTA
      DELMIN=1.0D-9*DELMAX
      GO TO 650 
C 
C  INCREMENT TIME, UPDATE SOURCES, AND SOLVE NEXT TIMEPOINT 
C 
  600 TIME=TIME+DELTA 
      CALL SORUPD 
      IF (NOGO.NE.0) GO TO 950
      CALL GETCJE 
      IF ((MAXTIM-ITIME).LE.LIMTIM) GO TO 920 
      IF ((ITL5.NE.0).AND.(NUMNIT.GE.ITL5)) GO TO 905 
      CALL COMCOF 
      IF (INITF.NE.5) INITF=6 
      ITRLIM=ITL4 
      IF ((NUMTP.EQ.0).AND.(NOSOLV.NE.0)) ITRLIM=ITL1 
      CALL ITER8(ITRLIM)
      NUMNIT=NUMNIT+ITERNO
      NUMTP=NUMTP+1 
      IF (NUMTP.NE.1) GO TO 605 
      IF (NOLX2.EQ.0) CALL COPY8(VALUE(LX1+1),VALUE(LX2+1),NXTRM) 
      IF (NOLX3.EQ.0) CALL COPY8(VALUE(LX1+1),VALUE(LX3+1),NXTRM) 
$PAGE 
C 
C.. NOTE THAT TIME-POINT IS CUT WHEN ITRLIM EXCEEDED REGARDLESS 
C.. OF WHICH TIME-STEP CONTOL IS SPECIFIED THRU 'LVLTIM'. 
C 
  605 IF (IGOOF.EQ.0) GO TO 610 
      JORD=IORD 
      IORD=1
      IF (JORD.GE.5) CALL CLRMEM(LX7) 
      IF (JORD.GE.4) CALL CLRMEM(LX6) 
      IF (JORD.GE.3) CALL CLRMEM(LX5) 
      IF ((JORD.GE.2).AND.(METHOD.NE.1)) CALL CLRMEM(LX4) 
      IGOOF=0 
      TIME=TIME-DELTA 
      DELTA=DELTA/8.0D0 
      GO TO 620 
  610 DELNEW=DELTA
      IF (NUMTP.EQ.1) GO TO 630 
      CALL TRUNC(DELNEW)
      IF (DELNEW.GE.(0.9D0*DELTA)) GO TO 630
      TIME=TIME-DELTA 
      DELTA=DELNEW
  620 NUMRTP=NUMRTP+1 
      IBKFLG=0
      DELOLD(1)=DELTA 
      IF (DELTA.GE.DELMIN) GO TO 600
      TIME=TIME+DELTA 
      GO TO 900 
$PAGE 
C 
C  DETERMINE ORDER OF INTEGRATION METHOD
C 
C...  SKIP IF TRAPEZOIDAL ALGORITHM USED
C 
  630 IF ((METHOD.EQ.1).AND.(IORD.EQ.2)) GO TO 650
      IF (NUMTP.EQ.1) GO TO 650 
      ORDRAT=1.05D0 
      IF (IORD.GT.1) GO TO 635
      IORD=2
      CALL TRUNC(DELNEW)
      IORD=1
      IF ((DELNEW/DELTA).LE.ORDRAT) GO TO 650 
      IF (MAXORD.LE.1) GO TO 650
      IORD=2
      IF (METHOD.EQ.1) GO TO 650
      CALL GETM8(LX4,NXTRM) 
      GO TO 650 
  635 IF (IORD.LT.MAXORD) GO TO 640 
      IORD=IORD-1 
      CALL TRUNC(DELNEW)
      IORD=IORD+1 
      IF ((DELNEW/DELTA).LE.ORDRAT) GO TO 650 
      GO TO 642 
  640 IORD=IORD-1 
      CALL TRUNC(DELNEW)
      IORD=IORD+1 
      IF ((DELNEW/DELTA).LE.ORDRAT) GO TO 645 
  642 IORD=IORD-1 
      IF (IORD.EQ.1) CALL CLRMEM(LX4) 
      IF (IORD.EQ.2) CALL CLRMEM(LX5) 
      IF (IORD.EQ.3) CALL CLRMEM(LX6) 
      IF (IORD.EQ.4) CALL CLRMEM(LX7) 
      GO TO 650 
  645 IORD=IORD+1 
      CALL TRUNC(DELNEW)
      IORD=IORD-1 
      IF ((DELNEW/DELTA).LE.ORDRAT) GO TO 650 
      IORD=IORD+1 
      IF (IORD.EQ.2) CALL GETM8(LX4,NXTRM)
      IF (IORD.EQ.3) CALL GETM8(LX5,NXTRM)
      IF (IORD.EQ.4) CALL GETM8(LX6,NXTRM)
      IF (IORD.EQ.5) CALL GETM8(LX7,NXTRM)
$PAGE 
C 
C  STORE OUTPUTS
C 
  650 IF ((TIME+DELTA).LE.TSTART) GO TO 685 
      IF ((NUMTP.EQ.0).AND.(NOSOLV.NE.0)) GO TO 685 
      CALL EXTMEM(LOUTPT,NUMOUT)
      LOCO=LOUTPT+ICALC*NUMOUT
      ICALC=ICALC+1 
      VALUE(LOCO+1)=TIME
      LOC=LOCATE(42)
  670 IF (LOC.EQ.0) GO TO 682 
      IF (NODPLC(LOC+5).NE.0) GO TO 680 
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      ISEQ=NODPLC(LOC+4)
      VALUE(LOCO+ISEQ)=VALUE(LVNIM1+NODE1)-VALUE(LVNIM1+NODE2)
      LOC=NODPLC(LOC) 
      GO TO 670 
  680 IPTR=NODPLC(LOC+2)
      IPTR=NODPLC(IPTR+6) 
      ISEQ=NODPLC(LOC+4)
      VALUE(LOCO+ISEQ)=VALUE(LVNIM1+IPTR) 
      LOC=NODPLC(LOC) 
      GO TO 670 
  682 IF (IPOSTP.EQ.0) GO TO 684  
      VALUE(IBUFF+1)=TIME 
      CALL COPY8(VALUE(LVNIM1+2),VALUE(IBUFF+2),NUNODS-1) 
      IF (NUMCUR.NE.0) CALL COPY8(VALUE(LVNIM1+LOCCUR+1), 
     1  VALUE(IBUFF+NUNODS+1),NUMCUR) 
  684 CONTINUE
C 
C  UPDATE TRANSMISSION LINE DELAY TABLE 
C 
  685 IF (JELCNT(17).EQ.0) GO TO 694
      CALL SIZMEM(LTD,LTDSIZ) 
      NUMTD=LTDSIZ/NTLIN
      IF (NUMTD.LE.3) GO TO 689 
      BAKTIM=TIME-TDMAX 
      IF (BAKTIM.LT.0.0D0) GO TO 689
      LCNTR=0 
      LTEMP=LTD 
      DO 686 I=1,NUMTD
      IF (VALUE(LTEMP+1).GE.BAKTIM) GO TO 687 
      LTEMP=LTEMP+NTLIN 
      LCNTR=LCNTR+1 
  686 CONTINUE
      GO TO 689 
$PAGE 
  687 IF (LCNTR.LE.2) GO TO 689 
      LCNTR=LCNTR-2 
      NWORDS=LCNTR*NTLIN
      LTEMP=LTEMP-NTLIN-NTLIN 
      CALL COPY8(VALUE(LTEMP+1),VALUE(LTD+1),LTDSIZ-NWORDS) 
      CALL RELMEM(LTD,NWORDS) 
      CALL SIZMEM(LTD,LTDSIZ) 
  689 CALL EXTMEM(LTD,NTLIN)
      LTDPTR=LTD+LTDSIZ 
      VALUE(LTDPTR+1)=TIME
      LOC=LOCATE(17)
  690 IF (LOC.EQ.0) GO TO 693 
      LOCV=NODPLC(LOC+1)
      Z0=VALUE(LOCV+1)
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NODE3=NODPLC(LOC+4) 
      NODE4=NODPLC(LOC+5) 
      IBR1=NODPLC(LOC+8)
      IBR2=NODPLC(LOC+9)
      LSPOT=NODPLC(LOC+30)+LTDPTR 
      IF ((INITF.EQ.5).AND.(NOSOLV.NE.0)) GO TO 691 
      VALUE(LSPOT)=VALUE(LVNIM1+NODE3)-VALUE(LVNIM1+NODE4)
     1   +VALUE(LVNIM1+IBR2)*Z0 
      VALUE(LSPOT+1)=VALUE(LVNIM1+NODE1)-VALUE(LVNIM1+NODE2)
     1   +VALUE(LVNIM1+IBR1)*Z0 
      GO TO 692 
  691 VALUE(LSPOT)=VALUE(LOCV+7)+VALUE(LOCV+8)*Z0 
      VALUE(LSPOT+1)=VALUE(LOCV+5)+VALUE(LOCV+6)*Z0 
  692 LOC=NODPLC(LOC) 
      GO TO 690 
C 
C  ADD TWO *FAKE* BACKPOINTS TO LTD FOR INTERPOLATION NEAR TIME=0.0D0 
C 
  693 IF (NUMTD.NE.0) GO TO 694 
      CALL EXTMEM(LTD,NTLIN+NTLIN)
      CALL COPY8(VALUE(LTD+1),VALUE(LTD+NTLIN+1),NTLIN) 
      CALL COPY8(VALUE(LTD+1),VALUE(LTD+2*NTLIN+1),NTLIN) 
      VALUE(LTD+2*NTLIN+1)=TIME 
      VALUE(LTD+NTLIN+1)=TIME-DELTA 
      VALUE(LTD+1)=TIME-DELTA-DELTA 
$PAGE 
C 
C  ROTATE STATE VECTOR STORAGE
C 
C.. TIME-POINT ACCEPTED 
C 
  694 CALL COPY8(DELOLD(1),DELOLD(2),6) 
      DELTA=DELNEW
      DELOLD(1)=DELTA 
      GO TO (710,706,702,698,696,696), IORD 
  696 CALL PTRMEM(LX7,ITEMP)
      CALL PTRMEM(LX6,LX7)
      GO TO 700 
  698 CALL PTRMEM(LX6,ITEMP)
  700 CALL PTRMEM(LX5,LX6)
      GO TO 704 
  702 CALL PTRMEM(LX5,ITEMP)
  704 CALL PTRMEM(LX4,LX5)
      GO TO 708 
  706 IF (METHOD.EQ.1) GO TO 710
      CALL PTRMEM(LX4,ITEMP)
  708 CALL PTRMEM(LX3,LX4)
      GO TO 713 
  710 IF (NOLX3.EQ.0) GO TO 712 
      IF (NOLX2.EQ.0) GO TO 711 
      CALL PTRMEM(LX1,ITEMP)
      GO TO 714 
  711 CALL PTRMEM(LX2,ITEMP)
      CALL PTRMEM(LX1,LX2)
      GO TO 714 
  712 CALL PTRMEM(LX3,ITEMP)
  713 CALL PTRMEM(LX2,LX3)
      CALL PTRMEM(LX1,LX2)
  714 CALL PTRMEM(LX0,LX1)
      CALL PTRMEM(ITEMP,LX0)
$PAGE 
C 
C  CHECK BREAKPOINTS
C 
  750 IF (IBKFLG.EQ.0) GO TO 760
C 
C.. JUST ACCEPTED ANALYSIS AT BREAKPOINT
C 
      JORD=IORD 
      IORD=1
      IF (JORD.GE.5) CALL CLRMEM(LX7) 
      IF (JORD.GE.4) CALL CLRMEM(LX6) 
      IF (JORD.GE.3) CALL CLRMEM(LX5) 
      IF ((JORD.GE.2).AND.(METHOD.NE.1)) CALL CLRMEM(LX4) 
      IBKFLG=0
      NBKPT=NBKPT+1 
      IF (NBKPT.GT.NUMBKP) GO TO 950
      TEMP=DMIN1(DELBKP,VALUE(LSBKPT+NBKPT)-TIME) 
      DELTA=DMIN1(DELTA,0.1D0*TEMP,DELMAX)
      IF (NUMTP.EQ.0) DELTA=DELTA/10.0D0
      DELOLD(1)=DELTA 
      GO TO 600 
  760 DEL1=VALUE(LSBKPT+NBKPT)-TIME 
      IF ((1.01D0*DELTA).LE.DEL1) GO TO 600 
      IBKFLG=1
      DELBKP=DELTA
      DELTA=DEL1
      DELOLD(1)=DELTA 
      GO TO 600 
C 
C  TRANSIENT ANALYSIS FAILED
C 
  900 WRITE (6,901) 
  901 FORMAT('1*ERROR*:  INTERNAL TIMESTEP TOO SMALL IN TRANSIENT ANALYS
     1IS'/) 
      GO TO 910 
  905 WRITE (6,906) ITL5
  906 FORMAT('1*ERROR*:  TRANSIENT ANALYSIS ITERATIONS EXCEED LIMIT OF '
     1,I5,/'0THIS LIMIT MAY BE OVERRIDDEN USING THE ITL5 PARAMETER ON TH
     2E .OPTION CARD')
  910 WRITE (6,911) TIME,DELTA,NUMNIT 
  911 FORMAT(1H0,10X,'TIME = ',1PD12.5,';  DELTA = ',D12.5,';  NUMNIT = 
     1',I6/)
      WRITE (6,916) 
  916 FORMAT(1H0/'0LAST NODE VOLTAGES:'/) 
      WRITE (6,AVHDR) (ANODE,AVLTG,I=1,NVPRLN)
      WRITE (6,AVFRM) (LPRN,NODPLC(JUNODE+I),VALUE(LVNIM1+I),I=2,NCNODS)
      GO TO 930 
  920 WRITE (6,921) TIME
  921 FORMAT('0*ERROR*:  CPU TIME LIMIT EXCEEDED IN TRANSIENT ANALYSIS '
     1   ,'AT TIME = ',1PD13.6/)
  930 NOGO=1
$PAGE 
C 
C  FINISHED WITH TRANSIENT ANALYSIS 
C 
  950 RSTATS(10)=RSTATS(10)+NUMNIT
      RSTATS(30)=RSTATS(30)+NUMTP 
      RSTATS(31)=RSTATS(31)+NUMRTP
      RSTATS(32)=RSTATS(32)+NUMNIT
      IF (IPOSTP.EQ.0) GO TO 1000 
C 
C  RETURN UNNEEDED MEMORY 
C 
 1000 IF (MODE.EQ.2) GO TO 1010 
      IF (MODEDC.NE.3) GO TO 1100 
 1010 CALL CLRMEM(LVNIM1) 
      CALL CLRMEM(LX0)
      CALL CLRMEM(LVN)
      CALL CLRMEM(LX1)
      IF (MEMPTR(MACINS)) CALL CLRMEM(MACINS) 
      IF (NOLX2.EQ.0) CALL CLRMEM(LX2)  
      CALL CLRMEM(LVNTMP) 
      IF ((MODE.EQ.1).AND.(MODEDC.EQ.3)) GO TO 1020 
      IF (NOLX3.EQ.0) CALL CLRMEM(LX3)  
      IF (MODE.EQ.1) GO TO 1020 
      CALL CLRMEM(LTD)
      IF (IORD.EQ.1) GO TO 1020 
      IF (METHOD.EQ.1) GO TO 1020 
      CALL CLRMEM(LX4)
      IF (IORD.EQ.2) GO TO 1020 
      CALL CLRMEM(LX5)
      IF (IORD.EQ.3) GO TO 1020 
      CALL CLRMEM(LX6)
      IF (IORD.EQ.4) GO TO 1020 
      CALL CLRMEM(LX7)
 1020 CALL EXTMEM(LOUTPT,2*NUMOUT)
 1100 IF (IPOSTP.NE.0) CALL CLRMEM(IBUFF) 
      CALL SECOND(T2) 
      RSTATS(LOCTIM)=RSTATS(LOCTIM)+T2-T1 
      RETURN
      END 
$EMA /BLANK/                            
      SUBROUTINE DIODE
     + ,92091-16005 REV.5000 <890418.1140>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE PROCESSES DIODES FOR DC AND TRANSIENT ANALYSES.
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      DIMENSION VDO(1),CDO(1),GDO(1),QD(1),CQD(1) 
      EQUIVALENCE (VDO(1),VALUE(1)),(CDO(1),VALUE(2)),
     1   (GDO(1),VALUE(3)),(QD(1),VALUE(4)),(CQD(1),VALUE(5)) 
C 
      LOC=LOCATE(11)
   10 IF (LOC.EQ.0) RETURN
      LOCV=NODPLC(LOC+1)
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NODE3=NODPLC(LOC+4) 
      LOCM=NODPLC(LOC+5)
      IOFF=NODPLC(LOC+6)
      LOCM=NODPLC(LOCM+1) 
      LOCT=NODPLC(LOC+11) 
C 
C  DC MODEL PARAMETERS
C 
      AREA=VALUE(LOCV+1)
      CSAT=VALUE(LOCM+1)*AREA 
      GSPR=VALUE(LOCM+2)*AREA 
      VTE=VALUE(LOCM+3)*VT
      BV=VALUE(LOCM+13) 
      VCRIT=VALUE(LOCM+18)
$PAGE 
C 
C  INITIALIZATION 
C 
      ICHECK=1
      GO TO (100,20,30,50,60,70),INITF
   20 IF (MODE.NE.1.OR.MODEDC.NE.2.OR.NOSOLV.EQ.0) GO TO 25 
      VD=VALUE(LOCV+2)
      GO TO 300 
   25 IF (IOFF.NE.0) GO TO 40 
      VD=VCRIT
      GO TO 300 
   30 IF (IOFF.EQ.0) GO TO 100
   40 VD=0.0D0
      GO TO 300 
   50 VD=VDO(LX0+LOCT)
      GO TO 300 
   60 VD=VDO(LX1+LOCT)
      GO TO 300 
   70 XFACT=DELTA/DELOLD(2) 
      VDO(LX0+LOCT)=VDO(LX1+LOCT) 
      VD=(1.0D0+XFACT)*VDO(LX1+LOCT)-XFACT*VDO(LX2+LOCT)
      CDO(LX0+LOCT)=CDO(LX1+LOCT) 
      GDO(LX0+LOCT)=GDO(LX1+LOCT) 
      GO TO 110 
C 
C  COMPUTE NEW NONLINEAR BRANCH VOLTAGE 
C 
  100 VD=VALUE(LVNIM1+NODE3)-VALUE(LVNIM1+NODE2)
  110 DELVD=VD-VDO(LX0+LOCT)
      CDHAT=CDO(LX0+LOCT)+GDO(LX0+LOCT)*DELVD 
C 
C  BYPASS IF SOLUTION HAS NOT CHANGED 
C 
      IF (INITF.EQ.6) GO TO 200 
      TOL=RELTOL*DMAX1(DABS(VD),DABS(VDO(LX0+LOCT)))+VNTOL
      IF (DABS(DELVD).GE.TOL) GO TO 200 
      TOL=RELTOL*DMAX1(DABS(CDHAT),DABS(CDO(LX0+LOCT)))+ABSTOL
      IF (DABS(CDHAT-CDO(LX0+LOCT)).GE.TOL) GO TO 200 
      VD=VDO(LX0+LOCT)
      CD=CDO(LX0+LOCT)
      GD=GDO(LX0+LOCT)
      GO TO 800 
C 
C  LIMIT NEW JUNCTION VOLTAGE 
C 
  200 VLIM=VTE+VTE
      IF (BV.EQ.0.0D0) GO TO 205  
      IF (VD.LT.DMIN1(0.0D0,-BV+10.0D0*VTE)) GO TO 210
  205 CALL PNJLIM(VD,VDO(LX0+LOCT),VTE,VCRIT,ICHECK)
      GO TO 300 
  210 VDTEMP=-(VD+BV) 
      CALL PNJLIM(VDTEMP,-(VDO(LX0+LOCT)+BV),VTE,VCRIT,ICHECK)
      VD=-(VDTEMP+BV) 
$PAGE 
C 
C  COMPUTE DC CURRENT AND DERIVITIVES 
C 
  300 IF (VD.LT.-5.0D0*VTE) GO TO 310 
      EVD=DEXP(VD/VTE)
      CD=CSAT*(EVD-1.0D0)+GMIN*VD 
      GD=CSAT*EVD/VTE+GMIN
      GO TO 330 
  310 IF (BV.EQ.0.0D0) GO TO 315
      IF (VD.LT.-BV) GO TO 320
  315 GD=-CSAT/VD+GMIN
      CD=GD*VD
      GO TO 330 
  320 EVREV=DEXP(-(BV+VD)/VT) 
      CD=-CSAT*(EVREV-1.0D0+BV/VT)
      GD=CSAT*EVREV/VT
  330 IF (MODE.NE.1) GO TO 500
      IF ((MODEDC.EQ.2).AND.(NOSOLV.NE.0)) GO TO 500
      IF (INITF.EQ.4) GO TO 500 
      GO TO 700 
C 
C  CHARGE STORAGE ELEMENTS
C 
  500 TAU=VALUE(LOCM+4) 
      CZERO=VALUE(LOCM+5)*AREA
      PB=VALUE(LOCM+6)
      XM=VALUE(LOCM+7)
      FCPB=VALUE(LOCM+12) 
      IF (VD.GE.FCPB) GO TO 510 
      ARG=1.0D0-VD/PB 
      SARG=DEXP(-XM*DLOG(ARG))
      QD(LX0+LOCT)=TAU*CD+PB*CZERO*(1.0D0-ARG*SARG)/(1.0D0-XM)
      CAPD=TAU*GD+CZERO*SARG
      GO TO 520 
  510 F1=VALUE(LOCM+15) 
      F2=VALUE(LOCM+16) 
      F3=VALUE(LOCM+17) 
      CZOF2=CZERO/F2
      QD(LX0+LOCT)=TAU*CD+CZERO*F1+CZOF2*(F3*(VD-FCPB)
     1   +(XM/(PB+PB))*(VD*VD-FCPB*FCPB)) 
      CAPD=TAU*GD+CZOF2*(F3+XM*VD/PB) 
C 
C  STORE SMALL-SIGNAL PARAMETERS
C 
  520 IF ((MODE.EQ.1).AND.(MODEDC.EQ.2).AND.(NOSOLV.NE.0)) GO TO 700
      IF (INITF.NE.4) GO TO 600 
      VALUE(LX0+LOCT+4)=CAPD
      GO TO 1000
$PAGE 
C 
C  TRANSIENT ANALYSIS 
C 
  600 IF (INITF.NE.5) GO TO 610 
      QD(LX1+LOCT)=QD(LX0+LOCT) 
  610 CALL INTGR8(GEQ,CEQ,CAPD,LOCT+3)
      GD=GD+GEQ 
      CD=CD+CQD(LX0+LOCT) 
      IF (INITF.NE.5) GO TO 700 
      CQD(LX1+LOCT)=CQD(LX0+LOCT) 
C 
C  CHECK CONVERGENCE
C 
  700 IF (INITF.NE.3) GO TO 710 
      IF (IOFF.EQ.0) GO TO 710
      GO TO 750 
  710 IF (ICHECK.EQ.1) GO TO 720
      TOL=RELTOL*DMAX1(DABS(CDHAT),DABS(CD))+ABSTOL 
      IF (DABS(CDHAT-CD).LE.TOL) GO TO 750
  720 NONCON=NONCON+1 
  750 VDO(LX0+LOCT)=VD
      CDO(LX0+LOCT)=CD
      GDO(LX0+LOCT)=GD
C 
C  LOAD CURRENT VECTOR
C 
  800 CDEQ=CD-GD*VD 
      VALUE(LVN+NODE2)=VALUE(LVN+NODE2)+CDEQ
      VALUE(LVN+NODE3)=VALUE(LVN+NODE3)-CDEQ
C 
C  LOAD MATRIX
C 
      LOCY=LVN+NODPLC(LOC+13) 
      VALUE(LOCY)=VALUE(LOCY)+GSPR
      LOCY=LVN+NODPLC(LOC+14) 
      VALUE(LOCY)=VALUE(LOCY)+GD
      LOCY=LVN+NODPLC(LOC+15) 
      VALUE(LOCY)=VALUE(LOCY)+GD+GSPR 
      LOCY=LVN+NODPLC(LOC+7)
      VALUE(LOCY)=VALUE(LOCY)-GSPR
      LOCY=LVN+NODPLC(LOC+8)
      VALUE(LOCY)=VALUE(LOCY)-GD
      LOCY=LVN+NODPLC(LOC+9)
      VALUE(LOCY)=VALUE(LOCY)-GSPR
      LOCY=LVN+NODPLC(LOC+10) 
      VALUE(LOCY)=VALUE(LOCY)-GD
 1000 LOC=NODPLC(LOC) 
      GO TO 10
      END 
$EMA /BLANK/
      SUBROUTINE EVPOLY( RESULT , ITYPE , LCOEF , NCOEF , LARG ,
     1  NARG,LEXP)
     + ,92091-16005 REV.5000 <890418.1140>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE EVALUATES A POLYNOMIAL.  LCOEF POINTS TO THE COEF- 
C FICIENTS, AND LARG POINTS TO THE VALUES OF THE POLYNOMIAL ARGUMENT(S).
C 
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      IF (ITYPE) 100,200,300
C 
C  INTEGRATION (POLYNOMIAL *MUST* BE ONE-DIMENSIONAL) 
C 
  100 RESULT=0.0D0
      ARG=1.0D0 
      ARG1=VALUE(LARG+1)
      DO 110 I=1,NCOEF
      ARG=ARG*ARG1
      RESULT=RESULT+VALUE(LCOEF+I)*ARG/DFLOAT(I)
  110 CONTINUE
      GO TO 1000
C 
C  EVALUATION OF THE POLYNOMIAL 
C 
  200 RESULT=VALUE(LCOEF+1) 
      IF (NCOEF.EQ.1) GO TO 1000
      CALL ZERO4(NODPLC(LEXP+1),NARG) 
      DO 220 I=2,NCOEF
      CALL NXTPWR(NODPLC(LEXP+1),NARG)
      IF (VALUE(LCOEF+I).EQ.0.0D0) GO TO 220
      ARG=1.0D0 
      DO 210 J=1,NARG 
      CALL EVTERM(VAL,VALUE(LARG+J),NODPLC(LEXP+J)) 
      ARG=ARG*VAL 
  210 CONTINUE
      RESULT=RESULT+VALUE(LCOEF+I)*ARG
  220 CONTINUE
      GO TO 1000
$PAGE 
C 
C  PARTIAL DERIVATIVE WITH RESPECT TO THE ITYPE*TH VARIABLE 
C 
  300 RESULT=0.0D0
      IF (NCOEF.EQ.1) GO TO 1000
      CALL ZERO4(NODPLC(LEXP+1),NARG) 
      DO 330 I=2,NCOEF
      CALL NXTPWR(NODPLC(LEXP+1),NARG)
      IF (NODPLC(LEXP+ITYPE).EQ.0) GO TO 330
      IF (VALUE(LCOEF+I).EQ.0.0D0) GO TO 330
      ARG=1.0D0 
      DO 320 J=1,NARG 
      IF (J.EQ.ITYPE) GO TO 310 
      CALL EVTERM(VAL,VALUE(LARG+J),NODPLC(LEXP+J)) 
      ARG=ARG*VAL 
      GO TO 320 
  310 CALL EVTERM(VAL,VALUE(LARG+J),NODPLC(LEXP+J)-1) 
      ARG=ARG*DFLOAT(NODPLC(LEXP+J))*VAL
  320 CONTINUE
      RESULT=RESULT+VALUE(LCOEF+I)*ARG
  330 CONTINUE
C 
C  FINISHED 
C 
 1000 RETURN
      END 
      SUBROUTINE EVTERM( VAL , ARG , IEXP ) 
     + ,92091-16005 REV.5000 <890418.1140>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE EVALUATES ONE TERM OF A POLYNOMIAL.
C 
      JEXP=IEXP+1 
      IF (JEXP.GE.6) GO TO 60 
      GO TO (10,20,30,40,50), JEXP
   10 VAL=1.0D0 
      GO TO 100 
   20 VAL=ARG 
      GO TO 100 
   30 VAL=ARG*ARG 
      GO TO 100 
   40 VAL=ARG*ARG*ARG 
      GO TO 100 
   50 VAL=ARG*ARG 
      VAL=VAL*VAL 
      GO TO 100 
   60 IF (ARG.EQ.0.0D0) GO TO 70
      ARGEXP=DFLOAT(IEXP)*DLOG(DABS(ARG)) 
      IF (ARGEXP.LT.-200.0D0) GO TO 70
      VAL=DEXP(ARGEXP)
      IF ((IEXP/2)*2.EQ.IEXP) GO TO 100 
      VAL=DSIGN(VAL,ARG)
      GO TO 100 
   70 VAL=0.0D0 
C 
C  FINISHED 
C 
  100 RETURN
      END 
      SUBROUTINE FETLIM( VNEW , VOLD , VTO )
     + ,92091-16005 REV.5000 <890418.1140>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE LIMITS THE PER-ITERATION CHANGE OF FET VOLTAGES. 
C 
C 
C      THREE REGIONS OF OPERATION ARE IDENTIFIED: 
C 
C                  V < VTO        DEFINITELY OFF
C            VTO <= V <= VTO+3.5D0    OFF OR ON DEPENDING ON VBS
C        VTO+3.5D0 < V              DEFINITELY ON 
C 
      VTSTHI=DABS(2.0D0*(VOLD-VTO))+2.0D0 
      VTSTLO=VTSTHI/2.0D0+2.0D0 
      VTOX=VTO+3.5D0
      DELV=VNEW-VOLD
C 
      IF (VOLD.LT.VTO) GO TO 300
      IF (VOLD.LT.VTOX) GO TO 200 
C 
C  ON ... 
C 
      IF (DELV.GT.0.0D0) GO TO 120
C 
C...  GOING OFF 
C 
      IF (VNEW.LT.VTOX) GO TO 110 
      IF (-DELV.LE.VTSTLO) GO TO 500
      VNEW=VOLD-VTSTLO
      GO TO 500 
  110 VNEW=DMAX1(VNEW,VTO+2.0D0)
      GO TO 500 
C 
C...  STAYING ON
C 
  120 IF (DELV.LT.VTSTHI) GO TO 500 
      VNEW=VOLD+VTSTHI
      GO TO 500 
$PAGE 
C 
C  MIDDLE REGION ...
C 
  200 IF (DELV.GT.0.0D0) GO TO 210
C 
C...  DECREASING
C 
      VNEW=DMAX1(VNEW,VTO-0.5D0)
      GO TO 500 
C 
C...  INCREASING
C 
  210 VNEW=DMIN1(VNEW,VTO+4.0D0)
      GO TO 500 
C 
C  OFF ...
C 
  300 IF (DELV.GT.0.0D0) GO TO 310
      IF (-DELV.LE.VTSTHI) GO TO 500
      VNEW=VOLD-VTSTHI
      GO TO 500 
  310 VTEMP=VTO+0.5D0 
      IF (VNEW.GT.VTEMP) GO TO 320
      IF (DELV.LE.VTSTLO) GO TO 500 
      VNEW=VOLD+VTSTLO
      GO TO 500 
  320 VNEW=VTEMP
C 
C  FINISHED 
C 
  500 RETURN
      END 
$EMA /BLANK/                            
      SUBROUTINE INTGR8( GEQ , CEQ , CAPVAL , LOCT )
     + ,92091-16005 REV.5000 <890418.1140>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE PERFORMS THE ACTUAL NUMERICAL INTEGRATION FOR EACH 
C CIRCUIT ELEMENT.
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      DIMENSION QCAP(1),CCAP(1) 
      EQUIVALENCE (QCAP(1),VALUE(1)),(CCAP(1),VALUE(2)) 
C 
      IF (METHOD.EQ.2) GO TO 100
C 
C  TRAPEZOIDAL ALGORITHM
C 
      IF (IORD.EQ.1) GO TO 100
      CCAP(LX0+LOCT)=-CCAP(LX1+LOCT)*AG(2)
     1   +AG(1)*(QCAP(LX0+LOCT)-QCAP(LX1+LOCT)) 
      GO TO 190 
$PAGE 
C 
C  GEARS ALGORITHM
C 
  100 GO TO (110,120,130,140,150,160), IORD 
  110 CCAP(LX0+LOCT)=AG(1)*QCAP(LX0+LOCT)+AG(2)*QCAP(LX1+LOCT)
      GO TO 190 
  120 CCAP(LX0+LOCT)=AG(1)*QCAP(LX0+LOCT)+AG(2)*QCAP(LX1+LOCT)
     1              +AG(3)*QCAP(LX2+LOCT) 
      GO TO 190 
  130 CCAP(LX0+LOCT)=AG(1)*QCAP(LX0+LOCT)+AG(2)*QCAP(LX1+LOCT)
     1              +AG(3)*QCAP(LX2+LOCT)+AG(4)*QCAP(LX3+LOCT)
      GO TO 190 
  140 CCAP(LX0+LOCT)=AG(1)*QCAP(LX0+LOCT)+AG(2)*QCAP(LX1+LOCT)
     1              +AG(3)*QCAP(LX2+LOCT)+AG(4)*QCAP(LX3+LOCT)
     2              +AG(5)*QCAP(LX4+LOCT) 
      GO TO 190 
  150 CCAP(LX0+LOCT)=AG(1)*QCAP(LX0+LOCT)+AG(2)*QCAP(LX1+LOCT)
     1              +AG(3)*QCAP(LX2+LOCT)+AG(4)*QCAP(LX3+LOCT)
     2              +AG(5)*QCAP(LX4+LOCT)+AG(6)*QCAP(LX5+LOCT)
      GO TO 190 
  160 CCAP(LX0+LOCT)=AG(1)*QCAP(LX0+LOCT)+AG(2)*QCAP(LX1+LOCT)
     1              +AG(3)*QCAP(LX2+LOCT)+AG(4)*QCAP(LX3+LOCT)
     2              +AG(5)*QCAP(LX4+LOCT)+AG(6)*QCAP(LX5+LOCT)
     3              +AG(7)*QCAP(LX6+LOCT) 
C 
C... CEQ IS THE EQUIVALENT CURRENT APPLICABLE TO LINEAR CAPACITANCE 
C    (INDUCTANCE) ONLY, I.E. Q=C*V
C 
  190 CEQ=CCAP(LX0+LOCT)-AG(1)*QCAP(LX0+LOCT) 
      GEQ=AG(1)*CAPVAL
      RETURN
      END 
$EMA /BLANK/
      SUBROUTINE ITER8( ITLIM ) 
     + ,92091-16005 REV.5000 <890418.1140>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE DRIVES THE NEWTON-RAPHSON ITERATION TECHNIQUE USED TO
C SOLVE THE SET OF NONLINEAR CIRCUIT EQUATIONS. 
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /MISCEL/ ATIME,APROG(3),ADATE,ATITLE(10),DEFL,DEFW,DEFAD,
     1  DEFAS,RSTATS(50),IWIDTH,LWIDTH,NOPAGE 
      COMMON /FLAGS/ IPRNTA,IPRNTL,IPRNTM,IPRNTN,IPRNTO,LIMTIM,LIMPTS,
     1   LVLCOD,LVLTIM,ITL1,ITL2,ITL3,ITL4,ITL5,IGOOF,NOGO,KEOF 
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
$PAGE 
      IGOOF=0 
      ITERNO=0
      NDRFLO=0
      NONCON=0
      IPASS=0 
C 
C  CONSTRUCT LINEAR EQUATIONS AND CHECK CONVERGENCE 
C 
   10 IVMFLG=0
      CALL LOAD 
   15 IF ((MODE.EQ.1).AND.(MODEDC.EQ.2).AND.(NOSOLV.NE.0)) GO TO 300
      ITERNO=ITERNO+1 
      GO TO (20,30,40,60,50,60),INITF 
   20 IF (MODE.NE.1) GO TO 22 
      CALL SIZMEM(NSNOD,NIC)
      IF (NIC.EQ.0) GO TO 22
      IF (IPASS.NE.0) NONCON=IPASS
      IPASS=0 
   22 IF (NONCON.EQ.0) GO TO 300
      GO TO 100 
   30 INITF=3 
      IF (LVLCOD.EQ.3) LVLCOD=2 
      IPIV=1
   40 IF (NONCON.EQ.0) INITF=1
      IPASS=1 
      GO TO 100 
   50 IF (ITERNO.GT.1) GO TO 60 
      IPIV=1
      IF (LVLCOD.EQ.3) LVLCOD=2 
   60 INITF=1 
$PAGE 
C 
C  SOLVE EQUATIONS FOR NEXT ITERATION 
C 
  100 IF (ITERNO.GE.ITLIM) GO TO 200
      IF (LVLCOD.EQ.3) GO TO 110
  102 CALL DCDCMP 
      IF (IGOOF.NE.0) GO TO 400 
      IF (LVLCOD.EQ.1) GO TO 105
      CALL CODGEN 
      CALL SIZMEM(MACINS,MACPTS)
      RSTATS(43)=MACPTS 
  105 CALL DCSOL
      GO TO 120 
  110 CALL SECOND(T1) 
      CALL CODEXC(NODPLC(MACINS+1),VALUE(LVN))
      CALL SECOND(T2) 
      RSTATS(44)=RSTATS(44)+T2-T1 
  120 IF (IGOOF.EQ.0) GO TO 130 
      IGOOF=0 
      IF (LVLCOD.NE.1) LVLCOD=2 
      IPIV=1
      CALL LOAD 
      GO TO 102 
  130 VALUE(LVN+1)=0.0D0
      DO 135 I=1,NSTOP
      J=NODPLC(ICSWPR+I)
      K=NODPLC(IRSWPF+J)
      VALUE(LVNTMP+K)=VALUE(LVNIM1+I) 
  135 CONTINUE
      CALL COPY8(VALUE(LVNTMP+1),VALUE(LVNIM1+1),NSTOP) 
      NTEMP=NONCON
      NONCON=0
      IF (NTEMP.GT.0) GO TO 150 
      IF (ITERNO.EQ.1) GO TO 150
      DO 140 I=2,NUMNOD 
      VOLD=VALUE(LVNIM1+I)
      VNEW=VALUE(LVN+I) 
      TOL=RELTOL*DMAX1(DABS(VOLD),DABS(VNEW))+VNTOL 
      IF (DABS(VOLD-VNEW).LE.TOL) GO TO 140 
      NONCON=NONCON+1 
  140 CONTINUE
  150 DO 160 I=1,NSTOP
      J=NODPLC(ICSWPR+I)
      K=NODPLC(IRSWPF+J)
      VALUE(LVNIM1+I)=VALUE(LVN+K)
  160 CONTINUE
C     WRITE(6,151) (VALUE(LVN+K),K=1,NSTOP) 
C 151 FORMAT(' SOLUTION: '/1P12D10.3) 
      GO TO 10
$PAGE 
C 
C  NO CONVERGENCE 
C 
  200 IGOOF=1 
  300 IF (NDRFLO.EQ.0) GO TO 400
      WRITE (6,301) NDRFLO
  301 FORMAT('0WARNING:  UNDERFLOW OCCURRED ',I4,' TIME(S)')
C 
C  FINISHED 
C 
  400 RETURN
      END 
$EMA /BLANK/                            
      SUBROUTINE JFET 
     + ,92091-16005 REV.5000 <890418.1140>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE PROCESSES JFETS FOR DC AND TRANSIENT ANALYSES. 
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      DIMENSION VGSO(1),VGDO(1),CGO(1),CDO(1),CGDO(1),GMO(1),GDSO(1), 
     1   GGSO(1),GGDO(1),QGS(1),CQGS(1),QGD(1),CQGD(1)
      EQUIVALENCE (VGSO(1),VALUE( 1)),(VGDO(1),VALUE( 2)),
     1            (CGO (1),VALUE( 3)),(CDO (1),VALUE( 4)),
     2            (CGDO(1),VALUE( 5)),(GMO (1),VALUE( 6)),
     3            (GDSO(1),VALUE( 7)),(GGSO(1),VALUE( 8)),
     4            (GGDO(1),VALUE( 9)),(QGS (1),VALUE(10)),
     5            (CQGS(1),VALUE(11)),(QGD (1),VALUE(12)),
     6            (CQGD(1),VALUE(13)) 
C 
      LOC=LOCATE(13)
   10 IF (LOC.EQ.0) RETURN
      LOCV=NODPLC(LOC+1)
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NODE3=NODPLC(LOC+4) 
      NODE4=NODPLC(LOC+5) 
      NODE5=NODPLC(LOC+6) 
      LOCM=NODPLC(LOC+7)
      IOFF=NODPLC(LOC+8)
      TYPE=NODPLC(LOCM+2) 
      LOCM=NODPLC(LOCM+1) 
      LOCT=NODPLC(LOC+19) 
$PAGE 
C 
C  DC MODEL PARAMETERS
C 
      AREA=VALUE(LOCV+1)
      VTO=VALUE(LOCM+1) 
      BETA=VALUE(LOCM+2)*AREA 
      XLAMB=VALUE(LOCM+3) 
      GDPR=VALUE(LOCM+4)*AREA 
      GSPR=VALUE(LOCM+5)*AREA 
      CSAT=VALUE(LOCM+9)*AREA 
      VCRIT=VALUE(LOCM+16)
C 
C  INITIALIZATION 
C 
      ICHECK=1
      GO TO (100,20,30,50,60,70), INITF 
   20 IF (MODE.NE.1.OR.MODEDC.NE.2.OR.NOSOLV.EQ.0) GO TO 25 
      VDS=TYPE*VALUE(LOCV+2)
      VGS=TYPE*VALUE(LOCV+3)
      VGD=VGS-VDS 
      GO TO 300 
   25 IF (IOFF.NE.0) GO TO 40 
      VGS=-1.0D0
      VGD=-1.0D0
      GO TO 300 
   30 IF (IOFF.EQ.0) GO TO 100
   40 VGS=0.0D0 
      VGD=0.0D0 
      GO TO 300 
   50 VGS=VGSO(LX0+LOCT)
      VGD=VGDO(LX0+LOCT)
      GO TO 300 
   60 VGS=VGSO(LX1+LOCT)
      VGD=VGDO(LX1+LOCT)
      GO TO 300 
   70 XFACT=DELTA/DELOLD(2) 
      VGSO(LX0+LOCT)=VGSO(LX1+LOCT) 
      VGS=(1.0D0+XFACT)*VGSO(LX1+LOCT)-XFACT*VGSO(LX2+LOCT) 
      VGDO(LX0+LOCT)=VGDO(LX1+LOCT) 
      VGD=(1.0D0+XFACT)*VGDO(LX1+LOCT)-XFACT*VGDO(LX2+LOCT) 
      CGO(LX0+LOCT)=CGO(LX1+LOCT) 
      CDO(LX0+LOCT)=CDO(LX1+LOCT) 
      CGDO(LX0+LOCT)=CGDO(LX1+LOCT) 
      GMO(LX0+LOCT)=GMO(LX1+LOCT) 
      GDSO(LX0+LOCT)=GDSO(LX1+LOCT) 
      GGSO(LX0+LOCT)=GGSO(LX1+LOCT) 
      GGDO(LX0+LOCT)=GGDO(LX1+LOCT) 
      GO TO 110 
$PAGE 
C 
C  COMPUTE NEW NONLINEAR BRANCH VOLTAGES
C 
  100 VGS=TYPE*(VALUE(LVNIM1+NODE2)-VALUE(LVNIM1+NODE5))
      VGD=TYPE*(VALUE(LVNIM1+NODE2)-VALUE(LVNIM1+NODE4))
  110 DELVGS=VGS-VGSO(LX0+LOCT) 
      DELVGD=VGD-VGDO(LX0+LOCT) 
      DELVDS=DELVGS-DELVGD
      CGHAT=CGO(LX0+LOCT)+GGDO(LX0+LOCT)*DELVGD+GGSO(LX0+LOCT)*DELVGS 
      CDHAT=CDO(LX0+LOCT)+GMO(LX0+LOCT)*DELVGS+GDSO(LX0+LOCT)*DELVDS
     1   -GGDO(LX0+LOCT)*DELVGD 
C 
C  BYPASS IF SOLUTION HAS NOT CHANGED 
C 
      IF (INITF.EQ.6) GO TO 200 
      TOL=RELTOL*DMAX1(DABS(VGS),DABS(VGSO(LX0+LOCT)))+VNTOL
      IF (DABS(DELVGS).GE.TOL) GO TO 200
      TOL=RELTOL*DMAX1(DABS(VGD),DABS(VGDO(LX0+LOCT)))+VNTOL
      IF (DABS(DELVGD).GE.TOL) GO TO 200
      TOL=RELTOL*DMAX1(DABS(CGHAT),DABS(CGO(LX0+LOCT)))+ABSTOL
      IF (DABS(CGHAT-CGO(LX0+LOCT)).GE.TOL) GO TO 200 
      TOL=RELTOL*DMAX1(DABS(CDHAT),DABS(CDO(LX0+LOCT)))+ABSTOL
      IF (DABS(CDHAT-CDO(LX0+LOCT)).GE.TOL) GO TO 200 
      VGS=VGSO(LX0+LOCT)
      VGD=VGDO(LX0+LOCT)
      VDS=VGS-VGD 
      CG=CGO(LX0+LOCT)
      CD=CDO(LX0+LOCT)
      CGD=CGDO(LX0+LOCT)
      GM=GMO(LX0+LOCT)
      GDS=GDSO(LX0+LOCT)
      GGS=GGSO(LX0+LOCT)
      GGD=GGDO(LX0+LOCT)
      GO TO 900 
C 
C  LIMIT NONLINEAR BRANCH VOLTAGES
C 
  200 ICHK1=1 
      CALL PNJLIM(VGS,VGSO(LX0+LOCT),VT,VCRIT,ICHECK) 
      CALL PNJLIM(VGD,VGDO(LX0+LOCT),VT,VCRIT,ICHK1)
      IF (ICHK1.EQ.1) ICHECK=1
      CALL FETLIM(VGS,VGSO(LX0+LOCT),VTO) 
      CALL FETLIM(VGD,VGDO(LX0+LOCT),VTO) 
$PAGE 
C 
C  DETERMINE DC CURRENT AND DERIVATIVES 
C 
  300 VDS=VGS-VGD 
      IF (VGS.GT.-5.0D0*VT) GO TO 310 
      GGS=-CSAT/VGS+GMIN
      CG=GGS*VGS
      GO TO 320 
  310 EVGS=DEXP(VGS/VT) 
      GGS=CSAT*EVGS/VT+GMIN 
      CG=CSAT*(EVGS-1.0D0)+GMIN*VGS 
  320 IF (VGD.GT.-5.0D0*VT) GO TO 330 
      GGD=-CSAT/VGD+GMIN
      CGD=GGD*VGD 
      GO TO 340 
  330 EVGD=DEXP(VGD/VT) 
      GGD=CSAT*EVGD/VT+GMIN 
      CGD=CSAT*(EVGD-1.0D0)+GMIN*VGD
  340 CG=CG+CGD 
C 
C  COMPUTE DRAIN CURRENT AND DERIVITIVES FOR NORMAL MODE
C 
  400 IF (VDS.LT.0.0D0) GO TO 450 
      VGST=VGS-VTO
C 
C  NORMAL MODE, CUTOFF REGION 
C 
      IF (VGST.GT.0.0D0) GO TO 410
      CDRAIN=0.0D0
      GM=0.0D0
      GDS=0.0D0 
      GO TO 490 
C 
C  NORMAL MODE, SATURATION REGION 
C 
  410 BETAP=BETA*(1.0D0+XLAMB*VDS)
      TWOB=BETAP+BETAP
      IF (VGST.GT.VDS) GO TO 420
      CDRAIN=BETAP*VGST*VGST
      GM=TWOB*VGST
      GDS=XLAMB*BETA*VGST*VGST
      GO TO 490 
C 
C  NORMAL MODE, LINEAR REGION 
C 
  420 CDRAIN=BETAP*VDS*(VGST+VGST-VDS)
      GM=TWOB*VDS 
      GDS=TWOB*(VGST-VDS)+XLAMB*BETA*VDS*(VGST+VGST-VDS)
      GO TO 490 
$PAGE 
C 
C  COMPUTE DRAIN CURRENT AND DERIVITIVES FOR INVERSE MODE 
C 
  450 VGDT=VGD-VTO
C 
C  INVERSE MODE, CUTOFF REGION
C 
      IF (VGDT.GT.0.0D0) GO TO 460
      CDRAIN=0.0D0
      GM=0.0D0
      GDS=0.0D0 
      GO TO 490 
C 
C  INVERSE MODE, SATURATION REGION
C 
  460 BETAP=BETA*(1.0D0-XLAMB*VDS)
      TWOB=BETAP+BETAP
      IF (VGDT.GT.-VDS) GO TO 470 
      CDRAIN=-BETAP*VGDT*VGDT 
      GM=-TWOB*VGDT 
      GDS=XLAMB*BETA*VGDT*VGDT-GM 
      GO TO 490 
C 
C  INVERSE MODE, LINEAR REGION
C 
  470 CDRAIN=BETAP*VDS*(VGDT+VGDT+VDS)
      GM=TWOB*VDS 
      GDS=TWOB*VGDT-XLAMB*BETA*VDS*(VGDT+VGDT+VDS)
C 
C  COMPUTE EQUIVALENT DRAIN CURRENT SOURCE
C 
  490 CD=CDRAIN-CGD 
      IF (MODE.NE.1) GO TO 500
      IF ((MODEDC.EQ.2).AND.(NOSOLV.NE.0)) GO TO 500
      IF (INITF.EQ.4) GO TO 500 
      GO TO 700 
C 
C  CHARGE STORAGE ELEMENTS
C 
  500 CZGS=VALUE(LOCM+6)*AREA 
      CZGD=VALUE(LOCM+7)*AREA 
      PHIB=VALUE(LOCM+8)
      TWOP=PHIB+PHIB
      FCPB=VALUE(LOCM+12) 
      FCPB2=FCPB*FCPB 
      F1=VALUE(LOCM+13) 
      F2=VALUE(LOCM+14) 
      F3=VALUE(LOCM+15) 
      CZGSF2=CZGS/F2
      CZGDF2=CZGD/F2
      IF (VGS.GE.FCPB) GO TO 510
      SARG=DSQRT(1.0D0-VGS/PHIB)
      QGS(LX0+LOCT)=TWOP*CZGS*(1.0D0-SARG)
      CAPGS=CZGS/SARG 
      GO TO 520 
  510 QGS(LX0+LOCT)=CZGS*F1+CZGSF2*(F3*(VGS-FCPB) 
     1   +(VGS*VGS-FCPB2)/(TWOP+TWOP))
      CAPGS=CZGSF2*(F3+VGS/TWOP)
  520 IF (VGD.GE.FCPB) GO TO 530
      SARG=DSQRT(1.0D0-VGD/PHIB)
      QGD(LX0+LOCT)=TWOP*CZGD*(1.0D0-SARG)
      CAPGD=CZGD/SARG 
      GO TO 560 
  530 QGD(LX0+LOCT)=CZGD*F1+CZGDF2*(F3*(VGD-FCPB) 
     1   +(VGD*VGD-FCPB2)/(TWOP+TWOP))
      CAPGD=CZGDF2*(F3+VGD/TWOP)
C 
C  STORE SMALL-SIGNAL PARAMETERS
C 
  560 IF ((MODE.EQ.1).AND.(MODEDC.EQ.2).AND.(NOSOLV.NE.0)) GO TO 700
      IF (INITF.NE.4) GO TO 600 
      VALUE(LX0+LOCT+9)=CAPGS 
      VALUE(LX0+LOCT+11)=CAPGD
      GO TO 1000
C 
C  TRANSIENT ANALYSIS 
C 
  600 IF (INITF.NE.5) GO TO 610 
      QGS(LX1+LOCT)=QGS(LX0+LOCT) 
      QGD(LX1+LOCT)=QGD(LX0+LOCT) 
  610 CALL INTGR8(GEQ,CEQ,CAPGS,LOCT+9) 
      GGS=GGS+GEQ 
      CG=CG+CQGS(LX0+LOCT)
      CALL INTGR8(GEQ,CEQ,CAPGD,LOCT+11)
      GGD=GGD+GEQ 
      CG=CG+CQGD(LX0+LOCT)
      CD=CD-CQGD(LX0+LOCT)
      CGD=CGD+CQGD(LX0+LOCT)
      IF (INITF.NE.5) GO TO 700 
      CQGS(LX1+LOCT)=CQGS(LX0+LOCT) 
      CQGD(LX1+LOCT)=CQGD(LX0+LOCT) 
C 
C  CHECK CONVERGENCE
C 
  700 IF (INITF.NE.3) GO TO 710 
      IF (IOFF.EQ.0) GO TO 710
      GO TO 750 
  710 IF (ICHECK.EQ.1) GO TO 720
      TOL=RELTOL*DMAX1(DABS(CGHAT),DABS(CG))+ABSTOL 
      IF (DABS(CGHAT-CG).GE.TOL) GO TO 720
      TOL=RELTOL*DMAX1(DABS(CDHAT),DABS(CD))+ABSTOL 
      IF (DABS(CDHAT-CD).LE.TOL) GO TO 750
  720 NONCON=NONCON+1 
  750 VGSO(LX0+LOCT)=VGS
      VGDO(LX0+LOCT)=VGD
      CGO(LX0+LOCT)=CG
      CDO(LX0+LOCT)=CD
      CGDO(LX0+LOCT)=CGD
      GMO(LX0+LOCT)=GM
      GDSO(LX0+LOCT)=GDS
      GGSO(LX0+LOCT)=GGS
      GGDO(LX0+LOCT)=GGD
C 
C  LOAD CURRENT VECTOR
C 
  900 CEQGD=TYPE*(CGD-GGD*VGD)
      CEQGS=TYPE*((CG-CGD)-GGS*VGS) 
      CDREQ=TYPE*((CD+CGD)-GDS*VDS-GM*VGS)
      VALUE(LVN+NODE2)=VALUE(LVN+NODE2)-CEQGS-CEQGD 
      VALUE(LVN+NODE4)=VALUE(LVN+NODE4)-CDREQ+CEQGD 
      VALUE(LVN+NODE5)=VALUE(LVN+NODE5)+CDREQ+CEQGS 
C 
C  LOAD Y MATRIX
C 
      LOCY=LVN+NODPLC(LOC+20) 
      VALUE(LOCY)=VALUE(LOCY)+GDPR
      LOCY=LVN+NODPLC(LOC+21) 
      VALUE(LOCY)=VALUE(LOCY)+GGD+GGS 
      LOCY=LVN+NODPLC(LOC+22) 
      VALUE(LOCY)=VALUE(LOCY)+GSPR
      LOCY=LVN+NODPLC(LOC+23) 
      VALUE(LOCY)=VALUE(LOCY)+GDPR+GDS+GGD
      LOCY=LVN+NODPLC(LOC+24) 
      VALUE(LOCY)=VALUE(LOCY)+GSPR+GDS+GM+GGS 
      LOCY=LVN+NODPLC(LOC+9)
      VALUE(LOCY)=VALUE(LOCY)-GDPR
      LOCY=LVN+NODPLC(LOC+10) 
      VALUE(LOCY)=VALUE(LOCY)-GGD 
      LOCY=LVN+NODPLC(LOC+11) 
      VALUE(LOCY)=VALUE(LOCY)-GGS 
      LOCY=LVN+NODPLC(LOC+12) 
      VALUE(LOCY)=VALUE(LOCY)-GSPR
      LOCY=LVN+NODPLC(LOC+13) 
      VALUE(LOCY)=VALUE(LOCY)-GDPR
      LOCY=LVN+NODPLC(LOC+14) 
      VALUE(LOCY)=VALUE(LOCY)+GM-GGD
      LOCY=LVN+NODPLC(LOC+15) 
      VALUE(LOCY)=VALUE(LOCY)-GDS-GM
      LOCY=LVN+NODPLC(LOC+16) 
      VALUE(LOCY)=VALUE(LOCY)-GGS-GM
      LOCY=LVN+NODPLC(LOC+17) 
      VALUE(LOCY)=VALUE(LOCY)-GSPR
      LOCY=LVN+NODPLC(LOC+18) 
      VALUE(LOCY)=VALUE(LOCY)-GDS 
 1000 LOC=NODPLC(LOC) 
      GO TO 10
      END 
$EMA /BLANK/                            
      SUBROUTINE LIMVDS( VNEW , VOLD )
     + ,92091-16005 REV.5000 <890418.1140>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE LIMITS THE PER-ITERATION CHANGE OF FET VDS.
C 
      IF (VOLD.LT.3.5D0) GO TO 200
C 
      IF (VNEW.LE.VOLD) GO TO 100 
      VNEW=DMIN1(VNEW,3.0D0*VOLD+2.0D0) 
      GO TO 500 
  100 IF (VNEW.LT.3.5D0) VNEW=DMAX1(VNEW,2.0D0) 
      GO TO 500 
C 
  200 IF (VNEW.LE.VOLD) GO TO 300 
      VNEW=DMIN1(VNEW,4.0D0)
      GO TO 500 
  300 VNEW=DMAX1(VNEW,-0.5D0) 
C 
  500 RETURN
      END 
$EMA /BLANK/                            
$ALIAS MOSFET='MOSFET_LINK' 
      SUBROUTINE LOAD 
     + ,92091-16005 REV.5000 <890418.1140>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE ZEROES-OUT AND THEN LOADS THE COEFFICIENT MATRIX.
C THE ACTIVE DEVICES AND THE CONTROLLED SOURCES ARE LOADED BY SEPARATE
C SUBROUTINES.
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /MISCEL/ ATIME,APROG(3),ADATE,ATITLE(10),DEFL,DEFW,DEFAD,
     1  DEFAS,RSTATS(50),IWIDTH,LWIDTH,NOPAGE 
      COMMON /FLAGS/ IPRNTA,IPRNTL,IPRNTM,IPRNTN,IPRNTO,LIMTIM,LIMPTS,
     1   LVLCOD,LVLTIM,ITL1,ITL2,ITL3,ITL4,ITL5,IGOOF,NOGO,KEOF 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      DIMENSION QCAP(1),CCAP(1) 
      EQUIVALENCE (QCAP(1),VALUE(1)),(CCAP(1),VALUE(2)) 
      DIMENSION FIND(1),VIND(1) 
      EQUIVALENCE (FIND(1),VALUE(1)),(VIND(1),VALUE(2)) 
C 
      CALL SECOND(T1) 
C 
C  ZERO Y MATRIX AND CURRENT VECTOR 
C 
      CALL ZERO8(VALUE(LVN+1),NSTOP+NTTBR)
$PAGE 
C 
C  RESISTORS
C 
      LOC=LOCATE(1) 
   20 IF (LOC.EQ.0) GO TO 30
      LOCV=NODPLC(LOC+1)
      VAL=VALUE(LOCV+1) 
      LOCY=LVN+NODPLC(LOC+6)
      VALUE(LOCY)=VALUE(LOCY)+VAL 
      LOCY=LVN+NODPLC(LOC+7)
      VALUE(LOCY)=VALUE(LOCY)+VAL 
      LOCY=LVN+NODPLC(LOC+4)
      VALUE(LOCY)=VALUE(LOCY)-VAL 
      LOCY=LVN+NODPLC(LOC+5)
      VALUE(LOCY)=VALUE(LOCY)-VAL 
      LOC=NODPLC(LOC) 
      GO TO 20
C 
C  CAPACITORS 
C 
   30 LOC=LOCATE(2) 
      IF ((MODE.EQ.1).AND.(MODEDC.NE.2)) GO TO 100
   40 IF (LOC.EQ.0) GO TO 100 
      LOCV=NODPLC(LOC+1)
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      LOCT=NODPLC(LOC+8)
      IPOLY=NODPLC(LOC+4) 
      IF (IPOLY.EQ.1) GO TO 43
      LCOEF=NODPLC(LOC+7) 
      CALL SIZMEM(NODPLC(LOC+7),NCOEF)
   43 VCAP=VALUE(LOCV+2)
      IF ((MODE.EQ.1).AND.(INITF.EQ.2)) GO TO 45
      IF ((NOSOLV.NE.0).AND.(INITF.EQ.5)) GO TO 45
      VCAP=VALUE(LVNIM1+NODE1)-VALUE(LVNIM1+NODE2)
   45 VALUE(LOCV+3)=VCAP
      IF (MODE.EQ.1) GO TO 60 
   47 IF (INITF.NE.6) GO TO 50
      QCAP(LX0+LOCT)=QCAP(LX1+LOCT) 
      GO TO 60
   50 IF (IPOLY.EQ.0) GO TO 53
      QCAP(LX0+LOCT)=VALUE(LOCV+1)*VCAP 
      IF (INITF.NE.5) GO TO 60
      IF (NOSOLV.NE.0) QCAP(LX0+LOCT)=VALUE(LOCV+1)*VALUE(LOCV+2) 
      QCAP(LX1+LOCT)=QCAP(LX0+LOCT) 
      GO TO 60
   53 CALL EVPOLY(QCAP(LX0+LOCT),-1,LCOEF,NCOEF,LOCV+2,1,LOC+8) 
      IF (INITF.NE.5) GO TO 60
      IF (NOSOLV.EQ.0) GO TO 55 
      VCAP=VALUE(LOCV+2)
      VALUE(LOCV+3)=VCAP
      CALL EVPOLY(QCAP(LX0+LOCT),-1,LCOEF,NCOEF,LOCV+2,1,LOC+8) 
   55 QCAP(LX1+LOCT)=QCAP(LX0+LOCT) 
   60 IF (IPOLY.EQ.1) GO TO 62
      CALL EVPOLY(VALUE(LOCV+1),0,LCOEF,NCOEF,LOCV+2,1,LOC+8) 
   62 IF (MODE.EQ.1) GO TO 90 
      CALL INTGR8(GEQ,CEQ,VALUE(LOCV+1),LOCT) 
      IF (IPOLY.EQ.1) GO TO 65
      CEQ=CEQ-GEQ*VCAP+AG(1)*QCAP(LX0+LOCT) 
   65 IF (INITF.NE.5) GO TO 70  
      CCAP(LX1+LOCT)=CCAP(LX0+LOCT) 
   70 LOCY=LVN+NODPLC(LOC+10) 
      VALUE(LOCY)=VALUE(LOCY)+GEQ 
      LOCY=LVN+NODPLC(LOC+11) 
      VALUE(LOCY)=VALUE(LOCY)+GEQ 
      LOCY=LVN+NODPLC(LOC+5)
      VALUE(LOCY)=VALUE(LOCY)-GEQ 
      LOCY=LVN+NODPLC(LOC+6)
      VALUE(LOCY)=VALUE(LOCY)-GEQ 
      VALUE(LVN+NODE1)=VALUE(LVN+NODE1)-CEQ 
      VALUE(LVN+NODE2)=VALUE(LVN+NODE2)+CEQ 
   90 LOC=NODPLC(LOC) 
      GO TO 40
C 
C  INDUCTORS
C 
  100 IF (JELCNT(3).EQ.0) GO TO 400 
      IF (MODE.EQ.1) GO TO 150
      IF (INITF.EQ.6) GO TO 150 
      LOC=LOCATE(3) 
  110 IF (LOC.EQ.0) GO TO 120 
      LOCV=NODPLC(LOC+1)
      IPTR=NODPLC(LOC+5)
      LOCT=NODPLC(LOC+11) 
      IPOLY=NODPLC(LOC+4) 
      IF (IPOLY.EQ.0) GO TO 115 
      FIND(LX0+LOCT)=VALUE(LOCV+1)*VALUE(LVNIM1+IPTR) 
      IF ((INITF.EQ.5).AND.(NOSOLV.NE.0)) 
     1   FIND(LX0+LOCT)=VALUE(LOCV+1)*VALUE(LOCV+2) 
      GO TO 118 
  115 LCOEF=NODPLC(LOC+10)
      CALL SIZMEM(NODPLC(LOC+10),NCOEF) 
      CIND=VALUE(LVNIM1+IPTR) 
      IF ((INITF.EQ.5).AND.(NOSOLV.NE.0)) CIND=VALUE(LOCV+2)
      VALUE(LOCV+3)=CIND
      CALL EVPOLY(FIND(LX0+LOCT),-1,LCOEF,NCOEF,LOCV+2,1,LOC+11)
  118 LOC=NODPLC(LOC) 
      GO TO 110 
  120 LOC=LOCATE(4) 
  130 IF (LOC.EQ.0) GO TO 150 
      LOCV=NODPLC(LOC+1)
      NL1=NODPLC(LOC+2) 
      NL2=NODPLC(LOC+3) 
      IPTR1=NODPLC(NL1+5) 
      IPTR2=NODPLC(NL2+5) 
      LOCT1=NODPLC(NL1+11)
      LOCT2=NODPLC(NL2+11)
      FIND(LX0+LOCT1)=FIND(LX0+LOCT1)+VALUE(LOCV+1)*VALUE(LVNIM1+IPTR2) 
      FIND(LX0+LOCT2)=FIND(LX0+LOCT2)+VALUE(LOCV+1)*VALUE(LVNIM1+IPTR1) 
      LOC=NODPLC(LOC) 
      GO TO 130 
  150 LOC=LOCATE(3) 
  160 IF (LOC.EQ.0) GO TO 300 
      LOCV=NODPLC(LOC+1)
      IPTR=NODPLC(LOC+5)
      LOCT=NODPLC(LOC+11) 
      IPOLY=NODPLC(LOC+4) 
      IF (IPOLY.EQ.1) GO TO 170 
      LCOEF=NODPLC(LOC+10)
      CALL SIZMEM(NODPLC(LOC+10),NCOEF) 
  170 CIND=VALUE(LVNIM1+IPTR) 
      IF ((NOSOLV.NE.0).AND.(INITF.EQ.5)) CIND=VALUE(LOCV+2)
      VALUE(LOCV+3)=CIND
  180 IF (MODE.NE.1) GO TO 200
      VEQ=0.0D0 
      REQ=0.0D0 
      GO TO 210 
  200 IF (INITF.NE.6) GO TO 205 
      FIND(LX0+LOCT)=FIND(LX1+LOCT) 
      GO TO 210 
  205 IF (INITF.NE.5) GO TO 210 
      FIND(LX1+LOCT)=FIND(LX0+LOCT) 
  210 IF (IPOLY.EQ.1) GO TO 220 
      CALL EVPOLY(VALUE(LOCV+1),0,LCOEF,NCOEF,LOCV+2,1,LOC+11)
  220 IF (MODE.EQ.1) GO TO 250
      CALL INTGR8(REQ,VEQ,VALUE(LOCV+1),LOCT) 
      IF (IPOLY.EQ.1) GO TO 250 
      VEQ=VEQ-REQ*CIND+AG(1)*FIND(LX0+LOCT) 
  250 VALUE(LVN+IPTR)=VEQ 
      IF (INITF.NE.5) GO TO 260 
      VIND(LX1+LOCT)=VIND(LX0+LOCT) 
  260 LOCY=LVN+NODPLC(LOC+13) 
      VALUE(LOCY)=-REQ
      LOCY=LVN+NODPLC(LOC+6)
      VALUE(LOCY)=1.0D0 
      LOCY=LVN+NODPLC(LOC+7)
      VALUE(LOCY)=-1.0D0
      LOCY=LVN+NODPLC(LOC+8)
      VALUE(LOCY)=1.0D0 
      LOCY=LVN+NODPLC(LOC+9)
      VALUE(LOCY)=-1.0D0
      LOC=NODPLC(LOC) 
      GO TO 160 
$PAGE 
C 
C  MUTUAL INDUCTANCES 
C 
  300 LOC=LOCATE(4) 
  310 IF (LOC.EQ.0) GO TO 400 
      LOCV=NODPLC(LOC+1)
      REQ=AG(1)*VALUE(LOCV+1) 
      LOCY=LVN+NODPLC(LOC+4)
      VALUE(LOCY)=-REQ
      LOCY=LVN+NODPLC(LOC+5)
      VALUE(LOCY)=-REQ
      LOC=NODPLC(LOC) 
      GO TO 310 
C 
C  NONLINEAR CONTROLLED SOURCES 
C 
  400 CALL NLCSRC 
C 
C  VOLTAGE SOURCES
C 
      LOC=LOCATE(9) 
  610 IF (LOC.EQ.0) GO TO 700 
      LOCV=NODPLC(LOC+1)
      IPTR=NODPLC(LOC+6)
      VALUE(LVN+IPTR)=VALUE(LOCV+1) 
      LOCY=LVN+NODPLC(LOC+7)
      VALUE(LOCY)=VALUE(LOCY)+1.0D0 
      LOCY=LVN+NODPLC(LOC+8)
      VALUE(LOCY)=VALUE(LOCY)-1.0D0 
      LOCY=LVN+NODPLC(LOC+9)
      VALUE(LOCY)=VALUE(LOCY)+1.0D0 
      LOCY=LVN+NODPLC(LOC+10) 
      VALUE(LOCY)=VALUE(LOCY)-1.0D0 
      LOC=NODPLC(LOC) 
      GO TO 610 
C 
C  CURRENT SOURCES
C 
  700 LOC=LOCATE(10)
  710 IF (LOC.EQ.0) GO TO 800 
      LOCV=NODPLC(LOC+1)
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      VAL=VALUE(LOCV+1) 
      VALUE(LVN+NODE1)=VALUE(LVN+NODE1)-VAL 
      VALUE(LVN+NODE2)=VALUE(LVN+NODE2)+VAL 
      LOC=NODPLC(LOC) 
      GO TO 710 
$PAGE 
C 
C  CALL DEVICE MODEL ROUTINES 
C 
  800 CALL DIODE
      CALL BJT
      CALL JFET 
      CALL MOSFET 
C 
C  TRANSMISSION LINES 
C 
      LOC=LOCATE(17)
  910 IF (LOC.EQ.0) GO TO 980 
      LOCV=NODPLC(LOC+1)
      Z0=VALUE(LOCV+1)
      Y0=1.0D0/Z0 
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NODE3=NODPLC(LOC+4) 
      NODE4=NODPLC(LOC+5) 
      IBR1=NODPLC(LOC+8)
      IBR2=NODPLC(LOC+9)
      LOCY=LVN+NODPLC(LOC+10) 
      VALUE(LOCY)=VALUE(LOCY)+Y0
      LOCY=LVN+NODPLC(LOC+11) 
      VALUE(LOCY)=-Y0 
      LOCY=LVN+NODPLC(LOC+12) 
      VALUE(LOCY)=-1.0D0
      LOCY=LVN+NODPLC(LOC+13) 
      VALUE(LOCY)=VALUE(LOCY)+Y0
      LOCY=LVN+NODPLC(LOC+14) 
      VALUE(LOCY)=-1.0D0
      LOCY=LVN+NODPLC(LOC+15) 
      VALUE(LOCY)=-Y0 
      LOCY=LVN+NODPLC(LOC+16) 
      VALUE(LOCY)=+Y0 
      LOCY=LVN+NODPLC(LOC+17) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LOC+18) 
      VALUE(LOCY)=+Y0 
      LOCY=LVN+NODPLC(LOC+19) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LOC+20) 
      VALUE(LOCY)=-1.0D0
      LOCY=LVN+NODPLC(LOC+23) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LOC+27) 
      VALUE(LOCY)=-1.0D0
      LOCY=LVN+NODPLC(LOC+28) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LOC+31) 
      VALUE(LOCY)=-Y0 
      LOCY=LVN+NODPLC(LOC+32) 
      VALUE(LOCY)=-Y0 
      IF (MODE.NE.1) GO TO 920
      LOCY=LVN+NODPLC(LOC+21) 
      VALUE(LOCY)=-1.0D0
      LOCY=LVN+NODPLC(LOC+22) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LOC+24) 
      VALUE(LOCY)=-(1.0D0-GMIN)*Z0
      LOCY=LVN+NODPLC(LOC+25) 
      VALUE(LOCY)=-1.0D0
      LOCY=LVN+NODPLC(LOC+26) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LOC+29) 
      VALUE(LOCY)=-(1.0D0-GMIN)*Z0
      GO TO 950 
  920 IF (INITF.NE.5) GO TO 930 
      IF (NOSOLV.NE.0) GO TO 925
      VALUE(LOCV+3)=VALUE(LVNIM1+NODE3)-VALUE(LVNIM1+NODE4) 
     1   +VALUE(LVNIM1+IBR2)*Z0 
      VALUE(LOCV+4)=VALUE(LVNIM1+NODE1)-VALUE(LVNIM1+NODE2) 
     1   +VALUE(LVNIM1+IBR1)*Z0 
      GO TO 930 
  925 VALUE(LOCV+3)=VALUE(LOCV+7)+VALUE(LOCV+8)*Z0
      VALUE(LOCV+4)=VALUE(LOCV+5)+VALUE(LOCV+6)*Z0
  930 VALUE(LVN+IBR1)=VALUE(LOCV+3) 
      VALUE(LVN+IBR2)=VALUE(LOCV+4) 
  950 LOC=NODPLC(LOC) 
      GO TO 910 
C 
C  INITIALIZE NODES 
C 
  980 IF (MODE.NE.1) GO TO 995  
      IF (INITF.NE.3.AND.INITF.NE.2) GO TO 1000 
      CALL SIZMEM(NSNOD,NIC)
      IF (NIC.EQ.0) GO TO 995 
      CALL SIZMEM(ICNOD,NTEST)
      IF (MODEDC.EQ.2.AND.NTEST.NE.0) GO TO 995 
      G=1.0D0 
      DO 990 I=1,NIC
      LOCY=LVN+NODPLC(NSMAT+I)
      VALUE(LOCY)=VALUE(LOCY)+G 
      NODE=NODPLC(NSNOD+I)
      VALUE(LVN+NODE)=VALUE(LVN+NODE)+VALUE(NSVAL+I)*G
  990 CONTINUE
$PAGE 
C 
C  TRANSIENT INITIAL CONDITIONS (UIC NOT SPECIFIED) 
C 
  995 IF (MODE.NE.1) GO TO 1000 
      IF (MODEDC.NE.2) GO TO 1000 
      IF (NOSOLV.NE.0) GO TO 1000 
      CALL SIZMEM(ICNOD,NIC)
      IF (NIC.EQ.0) GO TO 1000
      G=1.0D0 
      DO 996 I=1,NIC
      LOCY=LVN+NODPLC(ICMAT+I)
      VALUE(LOCY)=VALUE(LOCY)+G 
      NODE=NODPLC(ICNOD+I)
      VALUE(LVN+NODE)=VALUE(LVN+NODE)+VALUE(ICVAL+I)*G
  996 CONTINUE
C 
C  FINISHED 
C 
 1000 CALL SECOND(T2) 
      RSTATS(45)=RSTATS(45)+T2-T1 
      RETURN
      END 
