Return to dctrans.f CVS log

Up to [CSRG BSD Unix] / 3BSD / cmd / spice

Annotation of 3BSD/cmd/spice/dctrans.f, revision 1.1

1.1     ! root        1:       subroutine dctran
        !             2:       implicit double precision (a-h,o-z)
        !             3: c
        !             4: c
        !             5: c     this routine controls the dc transfer curve, dc operating point,
        !             6: c and transient analyses.  the variables mode and modedc (defined below)
        !             7: c determine exactly which analysis is performed.
        !             8: c
        !             9:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !            10:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !            11:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !            12:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !            13:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !            14:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !            15:       common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
        !            16:      1  defas,rstats(50),iwidth,lwidth,nopage
        !            17:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !            18:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !            19:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !            20:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !            21:      2   itemno,nosolv,ipostp,iscrch
        !            22:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
        !            23:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
        !            24:       common /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop,
        !            25:      1   kinel,kidin,kovar,kidout
        !            26:       common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg
        !            27:       common /cje/ maxtim,itime,icost
        !            28:       common /blank/ value(1000)
        !            29:       integer nodplc(64)
        !            30:       complex*16 cvalue(32)
        !            31:       equivalence (value(1),nodplc(1),cvalue(1))
        !            32: c
        !            33: c
        !            34:       dimension subtit(4,2)
        !            35:       dimension avhdr(3),avfrm(4)
        !            36:        data aendor /9.87654321d+27/
        !            37:       data avhdr / 8h( (2x,a4, 8h,3x,a7,3, 5hx)//) /
        !            38:       data avfrm / 8h( (1h ,a, 8h1,i3,1h), 8h,f10.4,3, 4hx)/) /
        !            39:       data anode, avltg / 4hnode, 7hvoltage /
        !            40:       data subtit / 8hsmall si, 8hgnal bia, 8hs soluti, 8hon      ,
        !            41:      1              8hinitial , 8htransien, 8ht soluti, 8hon      /
        !            42:       data lprn /1h(/
        !            43: c
        !            44: c      the variables *mode*, *modedc*, and *initf* are used by spice to
        !            45: c keep track of the state of the analysis.  the values of these flags
        !            46: c (and the corresponding meanings) are as follows:
        !            47: c
        !            48: c        flag    value    meaning
        !            49: c        ----    -----    -------
        !            50: c
        !            51: c        mode      1      dc analysis (subtype defined by *modedc*)
        !            52: c                  2      transient analysis
        !            53: c                  3      ac analysis (small signal)
        !            54: c
        !            55: c        modedc    1      dc operating point
        !            56: c                  2      initial operating point for transient analysis
        !            57: c                  3      dc transfer curve computation
        !            58: c
        !            59: c        initf     1      converge with 'off' devices allowed to float
        !            60: c                  2      initialize junction voltages
        !            61: c                  3      converge with 'off' devices held 'off'
        !            62: c                  4      store small-signal parameters away
        !            63: c                  5      first timepoint in transient analysis
        !            64: c                  6      prediction step
        !            65: c
        !            66: c note:  *modedc* is only significant if *mode* = 1.
        !            67: c
        !            68: c
        !            69: c  initialize
        !            70: c
        !            71:       call second(t1)
        !            72: c.. don't take any chances with lx3, set to large number
        !            73:       lx3=20000000
        !            74:       lx2=20000000
        !            75:       nolx2=0
        !            76:       nolx3=0
        !            77:       loctim=5
        !            78: c.. post-processing initialization
        !            79:       if(ipostp.eq.0) go to 1
        !            80:       numcur=jelcnt(9)
        !            81:       numpos=nunods+numcur
        !            82:       call getm8(ibuff,numpos)
        !            83:       numpos=numpos*4
        !            84:       if(numcur.eq.0) go to 1
        !            85:       loc=locate(9)
        !            86:       loccur=nodplc(loc+6)-1
        !            87: c...  set up format
        !            88:     1 nvprln=4+(lwidth-72)/19
        !            89:       nvprln=min0(nvprln,ncnods-1)
        !            90:       ipos=2
        !            91:       call alfnum(nvprln,avfrm,ipos)
        !            92:       ipos=2
        !            93:       call alfnum(nvprln,avhdr,ipos)
        !            94: c...  allocate storage
        !            95:       if (mode.eq.2) go to 5
        !            96:       need=4*nstop+nut+nlt+nxtrm
        !            97:       call avlm8(navl)
        !            98:       if(need.le.navl) go to 4
        !            99: c...  not enough memory for dc operating point analysis
        !           100:       write(6,3) need,navl
        !           101:     3 format('0insufficient memory available for dc analysis.',/
        !           102:      1' memory required ',i6,', memory available ',i6,'.')
        !           103:       nogo=1
        !           104:       go to 1100
        !           105:     4 call getm8(lvnim1,nstop)
        !           106:       call getm8(ndiag,nstop)
        !           107:       call getm8(lvn,nstop+nstop+nut+nlt)
        !           108:       call getm8(lx0,nxtrm)
        !           109:       if (modedc.ne.3) go to 15
        !           110:     5 call getm8(lx1,nxtrm)
        !           111:       if(nolx2.eq.0) call getm8(lx2,nxtrm)
        !           112:       if (mode.ne.2) go to 12
        !           113:       if(nolx3.eq.0) call getm8(lx3,nxtrm)
        !           114:       call getm8(ltd,0)
        !           115:    12 call getm8(loutpt,0)
        !           116:    15 call crunch
        !           117:       lynl=lvn+nstop
        !           118:       lyu=lynl+nstop
        !           119:       lyl=lyu+nut
        !           120:    20 if (mode.eq.2) go to 500
        !           121:       time=0.0d0
        !           122:       ag(1)=0.0d0
        !           123:       call sorupd
        !           124:       if (modedc.eq.3) go to 300
        !           125: c
        !           126: c  ....  single point dc analysis
        !           127: c
        !           128: c
        !           129: c  compute dc operating point
        !           130: c
        !           131:   100 initf=2
        !           132:       call iter8(itl1)
        !           133:       rstats(6)=rstats(6)+iterno
        !           134:       if (igoof.ne.0) go to 150
        !           135:       if (modedc.ne.1) go to 120
        !           136:       initf=4
        !           137:       call diode
        !           138:       call bjt
        !           139:       call jfet
        !           140:       call mosfet
        !           141: c
        !           142: c  print operating point
        !           143: c
        !           144:   120 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 1000
        !           145:       call title(-1,lwidth,1,subtit(1,modedc))
        !           146:       write (6,avhdr) (anode,avltg,i=1,nvprln)
        !           147:       write (6,avfrm) (lprn,nodplc(junode+i),value(lvnim1+i),i=2,ncnods)
        !           148:       go to 1000
        !           149: c
        !           150: c  no convergence
        !           151: c
        !           152:   150 nogo=1
        !           153:       write (6,151)
        !           154:   151 format('1*error*:  no convergence in dc analysis'/'0last node vol'
        !           155:      1   ,'tages:'/)
        !           156:       write (6,avhdr) (anode,avltg,i=1,nvprln)
        !           157:       write (6,avfrm) (lprn,nodplc(junode+i),value(lvnim1+i),i=2,ncnods)
        !           158:       go to 1000
        !           159: c
        !           160: c  ....  dc transfer curves
        !           161: c
        !           162:   300 numout=jelcnt(41)+1
        !           163:       if(ipostp.ne.0) call pheadr(atitle)
        !           164:       itemp=itcelm(1)
        !           165:       locs=nodplc(itemp+1)
        !           166:       temval=value(locs+1)
        !           167:       icvfl2=1
        !           168:       if(itcelm(2).eq.0) go to 310
        !           169:       itemp=itcelm(2)
        !           170:       locs2=nodplc(itemp+1)
        !           171:       temv2=value(locs2+1)
        !           172:       value(locs2+1)=tcstar(2)
        !           173:       temp=dabs((tcstop(2)-tcstar(2))/tcincr(2))+0.5d0
        !           174:       icvfl2=idint(temp)+1
        !           175:       icvfl2=max0(icvfl2,1)
        !           176:   310 delta=tcincr(1)
        !           177:       do 320 i=1,7
        !           178:       delold(i)=delta
        !           179:   320 continue
        !           180:       icvfl1=icvflg/icvfl2
        !           181:       value(locs+1)=tcstar(1)
        !           182:       icalc=0
        !           183:       ical2=0
        !           184:       loctim=3
        !           185:   340 initf=2
        !           186:       call iter8(itl1)
        !           187:       rstats(4)=rstats(4)+iterno
        !           188:       call copy8(value(lx0+1),value(lx1+1),nxtrm)
        !           189:       if(nolx2.eq.0) call copy8(value(lx0+1),value(lx2+1),nxtrm)
        !           190:       if (igoof.ne.0) go to 450
        !           191:       go to 360
        !           192:   350 call getcje
        !           193:       if ((maxtim-itime).le.limtim) go to 460
        !           194:       initf=6
        !           195:       call iter8(itl2)
        !           196:       rstats(4)=rstats(4)+iterno
        !           197:       if (igoof.ne.0) go to 340
        !           198: c
        !           199: c  store outputs
        !           200: c
        !           201:   360 call extmem(loutpt,numout)
        !           202:       loco=loutpt+icalc*numout
        !           203:       icalc=icalc+1
        !           204:       ical2=ical2+1
        !           205:       value(loco+1)=value(locs+1)
        !           206:       loc=locate(41)
        !           207:   370 if (loc.eq.0) go to 400
        !           208:       if (nodplc(loc+5).ne.0) go to 380
        !           209:       node1=nodplc(loc+2)
        !           210:       node2=nodplc(loc+3)
        !           211:       iseq=nodplc(loc+4)
        !           212:       value(loco+iseq)=value(lvnim1+node1)-value(lvnim1+node2)
        !           213:       loc=nodplc(loc)
        !           214:       go to 370
        !           215:   380 iptr=nodplc(loc+2)
        !           216:       iptr=nodplc(iptr+6)
        !           217:       iseq=nodplc(loc+4)
        !           218:       value(loco+iseq)=value(lvnim1+iptr)
        !           219:       loc=nodplc(loc)
        !           220:       go to 370
        !           221: c
        !           222: c  increment source value
        !           223: c
        !           224:   400 if(ipostp.eq.0) go to 410
        !           225:       value(ibuff+1)=value(locs+1)
        !           226:       call copy8(value(lvnim1+2),value(ibuff+2),nunods-1)
        !           227:       if(numcur.ne.0) call copy8(value(lvnim1+loccur+1),
        !           228:      1  value(ibuff+nunods+1),numcur)
        !           229:       call fwrite(value(ibuff+1),numpos)
        !           230:   410 if (icalc.ge.icvflg) go to 490
        !           231:       if(ical2.ge.icvfl1) go to 480
        !           232:       if(nolx2.ne.0) go to 420
        !           233:       call ptrmem(lx2,itemp)
        !           234:       call ptrmem(lx1,lx2)
        !           235:       go to 430
        !           236:   420 call ptrmem(lx1,itemp)
        !           237:   430 call ptrmem(lx0,lx1)
        !           238:       call ptrmem(itemp,lx0)
        !           239:       value(locs+1)=tcstar(1)+dfloat(ical2)*delta
        !           240:       go to 350
        !           241: c
        !           242: c  no convergence
        !           243: c
        !           244:   450 itemp=itcelm(1)
        !           245:       loce=nodplc(itemp+1)
        !           246:       write (6,451) value(loce),value(locs+1)
        !           247:   451 format('1*error*:  no convergence in dc transfer curves at ',a8,
        !           248:      1   ' = ',1pd10.3/'0last node voltages:'/)
        !           249:       write (6,avhdr) (anode,avltg,i=1,nvprln)
        !           250:       write (6,avfrm) (lprn,nodplc(junode+i),value(lvnim1+i),i=2,ncnods)
        !           251:       go to 470
        !           252:   460 write (6,461)
        !           253:   461 format('0*error*:  cpu time limit exceeded ... analysis stopped'/)
        !           254:   470 nogo=1
        !           255:       go to 490
        !           256: c... reset first sweep variable ... step second
        !           257:   480 ical2=0
        !           258:       value(locs+1)=tcstar(1)
        !           259:       value(locs2+1)=value(locs2+1)+tcincr(2)
        !           260:       go to 340
        !           261: c
        !           262: c  finished with dc transfer curves
        !           263: c
        !           264:   490 value(locs+1)=temval
        !           265:       if(itcelm(2).ne.0) value(locs2+1)=temv2
        !           266:       if(ipostp.eq.0) go to 1000
        !           267:       value(ibuff+1)=aendor
        !           268:       call fwrite(value(ibuff+1),numpos)
        !           269:       go to 1000
        !           270: c
        !           271: c  ....  transient analysis
        !           272: c
        !           273:   500 numout=jelcnt(42)+1
        !           274:       if(ipostp.ne.0) call pheadr(atitle)
        !           275:       numese=jelcnt(2)+jelcnt(3)+jelcnt(11)+jelcnt(12)+jelcnt(13)
        !           276:      1   +jelcnt(14)
        !           277:       if (numese.eq.0) delmax=dmin1(delmax,tstep)
        !           278:       initf=5
        !           279:       iord=1
        !           280:       loctim=9
        !           281:       icalc=0
        !           282:       numtp=0
        !           283:       numrtp=0
        !           284:       numnit=0
        !           285:       time=0.0d0
        !           286:       ibkflg=1
        !           287:       delbkp=delmax
        !           288:       nbkpt=1
        !           289:       delta=delmax
        !           290:       do 510 i=1,7
        !           291:       delold(i)=delta
        !           292:   510 continue
        !           293:       delmin=1.0d-9*delmax
        !           294:       go to 650
        !           295: c
        !           296: c  increment time, update sources, and solve next timepoint
        !           297: c
        !           298:   600 time=time+delta
        !           299:       call sorupd
        !           300:       if (nogo.ne.0) go to 950
        !           301:       call getcje
        !           302:       if ((maxtim-itime).le.limtim) go to 920
        !           303:       if ((itl5.ne.0).and.(numnit.ge.itl5)) go to 905
        !           304:       call comcof
        !           305:       if (initf.ne.5) initf=6
        !           306:       itrlim=itl4
        !           307:       if ((numtp.eq.0).and.(nosolv.ne.0)) itrlim=itl1
        !           308:       call iter8(itrlim)
        !           309:       if(itl5.ne.0) numnit=numnit+iterno
        !           310:       numtp=numtp+1
        !           311:       if (numtp.ne.1) go to 605
        !           312:       if(nolx2.eq.0) call copy8(value(lx1+1),value(lx2+1),nxtrm)
        !           313:       if(nolx3.eq.0) call copy8(value(lx1+1),value(lx3+1),nxtrm)
        !           314: c.. note that time-point is cut when itrlim exceeded regardless
        !           315: c.. of which time-step contol is specified thru 'lvltim'.
        !           316:   605 if (igoof.eq.0) go to 610
        !           317:       jord=iord
        !           318:       iord=1
        !           319:       if (jord.ge.5) call clrmem(lx7)
        !           320:       if (jord.ge.4) call clrmem(lx6)
        !           321:       if (jord.ge.3) call clrmem(lx5)
        !           322:       if ((jord.ge.2).and.(method.ne.1)) call clrmem(lx4)
        !           323:       igoof=0
        !           324:       time=time-delta
        !           325:       delta=delta/8.0d0
        !           326:       go to 620
        !           327:   610 delnew=delta
        !           328:       if (numtp.eq.1) go to 630
        !           329:       call trunc(delnew)
        !           330:       if (delnew.ge.(0.9d0*delta)) go to 630
        !           331:       time=time-delta
        !           332:       delta=delnew
        !           333:   620 numrtp=numrtp+1
        !           334:       ibkflg=0
        !           335:       delold(1)=delta
        !           336:       if (delta.ge.delmin) go to 600
        !           337:       time=time+delta
        !           338:       go to 900
        !           339: c.. time-point accepted
        !           340:   630 call copy8(delold(1),delold(2),6)
        !           341:       delta=delnew
        !           342:       delold(1)=delta
        !           343: c
        !           344: c  determine order of integration method
        !           345: c
        !           346: c...  skip if trapezoidal algorithm used
        !           347:       if ((method.eq.1).and.(iord.eq.2)) go to 650
        !           348:       if (numtp.eq.1) go to 650
        !           349:       ordrat=1.05d0
        !           350:       if (iord.gt.1) go to 635
        !           351:       iord=2
        !           352:       call trunc(delnew)
        !           353:       iord=1
        !           354:       if ((delnew/delta).le.ordrat) go to 650
        !           355:       if (maxord.le.1) go to 650
        !           356:       iord=2
        !           357:       if (method.eq.1) go to 650
        !           358:       call getm8(lx4,nxtrm)
        !           359:       go to 650
        !           360:   635 if (iord.lt.maxord) go to 640
        !           361:       iord=iord-1
        !           362:       call trunc(delnew)
        !           363:       iord=iord+1
        !           364:       if ((delnew/delta).le.ordrat) go to 650
        !           365:       go to 642
        !           366:   640 iord=iord-1
        !           367:       call trunc(delnew)
        !           368:       iord=iord+1
        !           369:       if ((delnew/delta).le.ordrat) go to 645
        !           370:   642 iord=iord-1
        !           371:       if (iord.eq.1) call clrmem(lx4)
        !           372:       if (iord.eq.2) call clrmem(lx5)
        !           373:       if (iord.eq.3) call clrmem(lx6)
        !           374:       if (iord.eq.4) call clrmem(lx7)
        !           375:       go to 650
        !           376:   645 iord=iord+1
        !           377:       call trunc(delnew)
        !           378:       iord=iord-1
        !           379:       if ((delnew/delta).le.ordrat) go to 650
        !           380:       iord=iord+1
        !           381:       if (iord.eq.2) call getm8(lx4,nxtrm)
        !           382:       if (iord.eq.3) call getm8(lx5,nxtrm)
        !           383:       if (iord.eq.4) call getm8(lx6,nxtrm)
        !           384:       if (iord.eq.5) call getm8(lx7,nxtrm)
        !           385: c
        !           386: c  store outputs
        !           387: c
        !           388:   650 if ((time+delta).le.tstart) go to 685
        !           389:       if ((numtp.eq.0).and.(nosolv.ne.0)) go to 685
        !           390:       call extmem(loutpt,numout)
        !           391:       loco=loutpt+icalc*numout
        !           392:       icalc=icalc+1
        !           393:       value(loco+1)=time
        !           394:       loc=locate(42)
        !           395:   670 if (loc.eq.0) go to 682
        !           396:       if (nodplc(loc+5).ne.0) go to 680
        !           397:       node1=nodplc(loc+2)
        !           398:       node2=nodplc(loc+3)
        !           399:       iseq=nodplc(loc+4)
        !           400:       value(loco+iseq)=value(lvnim1+node1)-value(lvnim1+node2)
        !           401:       loc=nodplc(loc)
        !           402:       go to 670
        !           403:   680 iptr=nodplc(loc+2)
        !           404:       iptr=nodplc(iptr+6)
        !           405:       iseq=nodplc(loc+4)
        !           406:       value(loco+iseq)=value(lvnim1+iptr)
        !           407:       loc=nodplc(loc)
        !           408:       go to 670
        !           409:   682 if(ipostp.eq.0) go to 684
        !           410:       value(ibuff+1)=time
        !           411:       call copy8(value(lvnim1+2),value(ibuff+2),nunods-1)
        !           412:       if(numcur.ne.0) call copy8(value(lvnim1+loccur+1),
        !           413:      1  value(ibuff+nunods+1),numcur)
        !           414:       call fwrite(value(ibuff+1),numpos)
        !           415:   684 continue
        !           416:   685 if (jelcnt(17).eq.0) go to 694
        !           417:       call sizmem(ltd,ltdsiz)
        !           418:       numtd=ltdsiz/ntlin
        !           419:       if (numtd.le.3) go to 689
        !           420:       baktim=time-tdmax
        !           421:       if (baktim.lt.0.0d0) go to 689
        !           422:       lcntr=0
        !           423:       ltemp=ltd
        !           424:       do 686 i=1,numtd
        !           425:       if (value(ltemp+1).ge.baktim) go to 687
        !           426:       ltemp=ltemp+ntlin
        !           427:       lcntr=lcntr+1
        !           428:   686 continue
        !           429:       go to 689
        !           430:   687 if (lcntr.le.2) go to 689
        !           431:       lcntr=lcntr-2
        !           432:       nwords=lcntr*ntlin
        !           433:       ltemp=ltemp-ntlin-ntlin
        !           434:       call copy8(value(ltemp+1),value(ltd+1),ltdsiz-nwords)
        !           435:       call relmem(ltd,nwords)
        !           436:       call sizmem(ltd,ltdsiz)
        !           437:   689 call extmem(ltd,ntlin)
        !           438:       ltdptr=ltd+ltdsiz
        !           439:       value(ltdptr+1)=time
        !           440:       loc=locate(17)
        !           441:   690 if (loc.eq.0) go to 693
        !           442:       locv=nodplc(loc+1)
        !           443:       z0=value(locv+1)
        !           444:       node1=nodplc(loc+2)
        !           445:       node2=nodplc(loc+3)
        !           446:       node3=nodplc(loc+4)
        !           447:       node4=nodplc(loc+5)
        !           448:       ibr1=nodplc(loc+8)
        !           449:       ibr2=nodplc(loc+9)
        !           450:       lspot=nodplc(loc+30)+ltdptr
        !           451:       if ((initf.eq.5).and.(nosolv.ne.0)) go to 691
        !           452:       value(lspot)=value(lvnim1+node3)-value(lvnim1+node4)
        !           453:      1   +value(lvnim1+ibr2)*z0
        !           454:       value(lspot+1)=value(lvnim1+node1)-value(lvnim1+node2)
        !           455:      1   +value(lvnim1+ibr1)*z0
        !           456:       go to 692
        !           457:   691 value(lspot)=value(locv+7)+value(locv+8)*z0
        !           458:       value(lspot+1)=value(locv+5)+value(locv+6)*z0
        !           459:   692 loc=nodplc(loc)
        !           460:       go to 690
        !           461: c
        !           462: c  add two *fake* backpoints to ltd for interpolation near time=0.0d0
        !           463: c
        !           464:   693 if (numtd.ne.0) go to 694
        !           465:       call extmem(ltd,ntlin+ntlin)
        !           466:       call copy8(value(ltd+1),value(ltd+ntlin+1),ntlin)
        !           467:       call copy8(value(ltd+1),value(ltd+2*ntlin+1),ntlin)
        !           468:       value(ltd+2*ntlin+1)=time
        !           469:       value(ltd+ntlin+1)=time-delta
        !           470:       value(ltd+1)=time-delta-delta
        !           471: c
        !           472: c  rotate state vector storage
        !           473: c
        !           474:   694 go to (710,706,702,698,696,696), iord
        !           475:   696 call ptrmem(lx7,itemp)
        !           476:       call ptrmem(lx6,lx7)
        !           477:       go to 700
        !           478:   698 call ptrmem(lx6,itemp)
        !           479:   700 call ptrmem(lx5,lx6)
        !           480:       go to 704
        !           481:   702 call ptrmem(lx5,itemp)
        !           482:   704 call ptrmem(lx4,lx5)
        !           483:       go to 708
        !           484:   706 if (method.eq.1) go to 710
        !           485:       call ptrmem(lx4,itemp)
        !           486:   708 call ptrmem(lx3,lx4)
        !           487:       go to 713
        !           488:   710 if(nolx3.eq.0) go to 712
        !           489:       if(nolx2.eq.0) go to 711
        !           490:       call ptrmem(lx1,itemp)
        !           491:       go to 714
        !           492:   711 call ptrmem(lx2,itemp)
        !           493:       call ptrmem(lx1,lx2)
        !           494:       go to 714
        !           495:   712 call ptrmem(lx3,itemp)
        !           496:   713 call ptrmem(lx2,lx3)
        !           497:       call ptrmem(lx1,lx2)
        !           498:   714 call ptrmem(lx0,lx1)
        !           499:       call ptrmem(itemp,lx0)
        !           500: c
        !           501: c  check breakpoints
        !           502: c
        !           503:   750 if (ibkflg.eq.0) go to 760
        !           504: c.. just accepted analysis at breakpoint
        !           505:       jord=iord
        !           506:       iord=1
        !           507:       if (jord.ge.5) call clrmem(lx7)
        !           508:       if (jord.ge.4) call clrmem(lx6)
        !           509:       if (jord.ge.3) call clrmem(lx5)
        !           510:       if ((jord.ge.2).and.(method.ne.1)) call clrmem(lx4)
        !           511:       ibkflg=0
        !           512:       nbkpt=nbkpt+1
        !           513:       if (nbkpt.gt.numbkp) go to 950
        !           514:       temp=dmin1(delbkp,value(lsbkpt+nbkpt)-time)
        !           515:       delta=dmin1(delta,0.1d0*temp,delmax)
        !           516:       if (numtp.eq.0) delta=delta/10.0d0
        !           517:       delold(1)=delta
        !           518:       go to 600
        !           519:   760 del1=value(lsbkpt+nbkpt)-time
        !           520:       if ((1.01d0*delta).le.del1) go to 600
        !           521:       ibkflg=1
        !           522:       delbkp=delta
        !           523:       delta=del1
        !           524:       delold(1)=delta
        !           525:       go to 600
        !           526: c
        !           527: c  transient analysis failed
        !           528: c
        !           529:   900 write (6,901)
        !           530:   901 format('1*error*:  internal timestep too small in transient analys
        !           531:      1is'/)
        !           532:       go to 910
        !           533:   905 write (6,906) itl5
        !           534:   906 format('1*error*:  transient analysis iterations exceed limit of '
        !           535:      1,i5,/'0this limit may be overridden using the itl5 parameter on th
        !           536:      2e .option card')
        !           537:   910 write (6,911) time,delta,numnit
        !           538:   911 format(1h0,10x,'time = ',1pd12.5,';  delta = ',d12.5,';  numnit =
        !           539:      1',i6/)
        !           540:       write (6,916)
        !           541:   916 format(1h0/'0last node voltages:'/)
        !           542:       write (6,avhdr) (anode,avltg,i=1,nvprln)
        !           543:       write (6,avfrm) (lprn,nodplc(junode+i),value(lvnim1+i),i=2,ncnods)
        !           544:       go to 930
        !           545:   920 write (6,921) time
        !           546:   921 format('0*error*:  cpu time limit exceeded in transient analysis '
        !           547:      1   ,'at time = ',1pd13.6/)
        !           548:   930 nogo=1
        !           549: c
        !           550: c  finished with transient analysis
        !           551: c
        !           552:   950 rstats(10)=rstats(10)+numnit
        !           553:       rstats(30)=rstats(30)+numtp
        !           554:       rstats(31)=rstats(31)+numrtp
        !           555:       rstats(32)=rstats(32)+numnit
        !           556:       if(ipostp.eq.0) go to 1000
        !           557:       value(ibuff+1)=aendor
        !           558:       call fwrite(value(ibuff+1),numpos)
        !           559: c
        !           560: c  return unneeded memory
        !           561: c
        !           562:  1000 if (mode.eq.2) go to 1010
        !           563:       if (modedc.ne.3) go to 1100
        !           564:  1010 call clrmem(lvnim1)
        !           565:       call clrmem(lx0)
        !           566:       call clrmem(ndiag)
        !           567:       call clrmem(lvn)
        !           568:       call clrmem(lx1)
        !           569:       if(nolx2.eq.0) call clrmem(lx2)
        !           570:       if ((mode.eq.1).and.(modedc.eq.3)) go to 1020
        !           571:       if(nolx3.eq.0) call clrmem(lx3)
        !           572:       if (mode.eq.1) go to 1020
        !           573:       call clrmem(ltd)
        !           574:       if (iord.eq.1) go to 1020
        !           575:       if (method.eq.1) go to 1020
        !           576:       call clrmem(lx4)
        !           577:       if (iord.eq.2) go to 1020
        !           578:       call clrmem(lx5)
        !           579:       if (iord.eq.3) go to 1020
        !           580:       call clrmem(lx6)
        !           581:       if (iord.eq.4) go to 1020
        !           582:       call clrmem(lx7)
        !           583:  1020 call extmem(loutpt,2*numout)
        !           584:  1100 if(ipostp.ne.0) call clrmem(ibuff)
        !           585:       call second(t2)
        !           586:       rstats(loctim)=rstats(loctim)+t2-t1
        !           587:       return
        !           588:       end
        !           589:       subroutine fwrite(iarray,numwds)
        !           590: c
        !           591: c   write onto 'punch' file numwds 16-bit words starting
        !           592: c   with location iarray(/1/)
        !           593: c
        !           594:       integer iarray(1)
        !           595: c
        !           596: c... data must be written out in 40 word (80 byte) chunks
        !           597: c
        !           598:       integer idata(20)
        !           599:       numwd4=(numwds+1)/2
        !           600:       numblk=(numwd4-1)/20+1
        !           601:       kntr=1
        !           602:       numlft=numwd4
        !           603:       do 10 i=1,numblk
        !           604:       kstop=min0(numlft,20)
        !           605:       call copy4(iarray(kntr),idata(1),kstop)
        !           606:       write(7) idata
        !           607:       kntr=kntr+20
        !           608:       numlft=numlft-20
        !           609:    10 continue
        !           610:       return
        !           611:       end
        !           612:       subroutine pheadr(aheadr)
        !           613:       implicit double precision (a-h,o-z)
        !           614:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !           615:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !           616:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !           617:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !           618:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !           619:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !           620:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !           621:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !           622:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !           623:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !           624:      2   itemno,nosolv,ipostp,iscrch
        !           625:       common /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop,
        !           626:      1   kinel,kidin,kovar,kidout
        !           627:       common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
        !           628:      1  defas,rstats(50),iwidth,lwidth,nopage
        !           629:       common /blank/ value(1000)
        !           630:       integer nodplc(64)
        !           631:       complex*16 cvalue(32)
        !           632:       integer*2 int2,nodpl2(128)
        !           633:       equivalence (value(1),nodpl2(1))
        !           634:       equivalence (value(1),nodplc(1),cvalue(1))
        !           635:       dimension aheadr(10)
        !           636: c
        !           637: c  put out the header records onto the post-processing file
        !           638: c  routine is used for all analysis modes (mode=1,2,3)
        !           639: c
        !           640:       dimension xtype(2)
        !           641:       data xtype /4htime,4hfreq/
        !           642:       data ablnk,aletv,aleti /1h ,1hv,1hi/
        !           643:       call getm8(ibuff,12)
        !           644:       call copy8(aheadr(1),value(ibuff+1),10)
        !           645:       value(ibuff+11)=adate
        !           646:       value(ibuff+12)=atime
        !           647:       call fwrite(value(ibuff+1),48)
        !           648:       numout=nunods+jelcnt(9)
        !           649:       info=3
        !           650:       call getm8(inames,numout)
        !           651:       call getm4(itypes,numout)
        !           652:       call getm4(iseqs,numout)
        !           653:       itype2=itypes*2
        !           654:       iseq2=iseqs*2
        !           655:       iknt=1
        !           656:       nodpl2(iseq2+1)=1
        !           657:       if(mode.ne.1) go to 10
        !           658:       loc=itcelm(1)
        !           659:       locv=nodplc(loc+1)
        !           660:       value(inames+1)=value(locv)
        !           661:       anam=ablnk
        !           662:       call move(anam,1,value(locv),1,1)
        !           663:       ityp=0
        !           664:       if(anam.eq.aletv) ityp=3
        !           665:       if(anam.eq.aleti) ityp=4
        !           666:       nodpl2(itype2+1)=ityp
        !           667:       go to 20
        !           668:    10 value(inames+1)=xtype(mode-1)
        !           669:       nodpl2(itype2+1)=mode-1
        !           670:    20 do 30 i=2,nunods
        !           671:       nodpl2(itype2+i)=3
        !           672:       nodpl2(iseq2+i)=i
        !           673:       value(inames+i)=ablnk
        !           674:       ipos=1
        !           675:       call alfnum(nodplc(junode+i),value(inames+i),ipos)
        !           676:    30 continue
        !           677:       loc=locate(9)
        !           678:       iknt=nunods
        !           679:    40 if(loc.eq.0) go to 50
        !           680:       iknt=iknt+1
        !           681:       nodpl2(itype2+iknt)=4
        !           682:       nodpl2(iseq2+iknt)=iknt
        !           683:       locv=nodplc(loc+1)
        !           684:       value(inames+iknt)=value(locv)
        !           685:       loc=nodplc(loc)
        !           686:       go to 40
        !           687:    50 int2=numout
        !           688:       call fwrite(int2,1)
        !           689:       int2=info
        !           690:       call fwrite(int2,1)
        !           691:       nwds=numout*4
        !           692:       call fwrite(value(inames+1),nwds)
        !           693:       call fwrite(nodpl2(itype2+1),numout)
        !           694:       call fwrite(nodpl2(iseq2+1),numout)
        !           695:       call clrmem(ibuff)
        !           696:       call clrmem(inames)
        !           697:       call clrmem(itypes)
        !           698:       call clrmem(iseqs)
        !           699:       return
        !           700:       end
        !           701:       subroutine comcof
        !           702:       implicit double precision (a-h,o-z)
        !           703: c
        !           704: c     this routine calculates the timestep-dependent terms used in the
        !           705: c numerical integration.
        !           706: c
        !           707:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !           708:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !           709:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !           710:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !           711:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !           712:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !           713:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !           714:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !           715:      2   itemno,nosolv,ipostp,iscrch
        !           716:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
        !           717:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
        !           718:       common /blank/ value(1000)
        !           719:       integer nodplc(64)
        !           720:       complex*16 cvalue(32)
        !           721:       equivalence (value(1),nodplc(1),cvalue(1))
        !           722:       dimension gmat(7,7)
        !           723: c
        !           724: c  compute coefficients for particular integration method
        !           725: c
        !           726:       if (method.ne.1) go to 5
        !           727:       if (iord.eq.1) go to 5
        !           728: c...  trapezoidal method
        !           729:       ag(1)=1.0d0/delta/(1.0d0-xmu)
        !           730:       ag(2)=xmu/(1.0d0-xmu)
        !           731:       go to 200
        !           732: c
        !           733: c  construct gear coefficient matrix
        !           734: c
        !           735:     5 istop=iord+1
        !           736:       call zero8(ag,istop)
        !           737:       ag(2)=-1.0d0
        !           738:       do 10 i=1,istop
        !           739:       gmat(1,i)=1.0d0
        !           740:    10 continue
        !           741:       do 20 i=2,istop
        !           742:       gmat(i,1)=0.0d0
        !           743:    20 continue
        !           744:       arg=0.0d0
        !           745:       do 40 i=2,istop
        !           746:       arg=arg+delold(i-1)
        !           747:       arg1=1.0d0
        !           748:       do 30 j=2,istop
        !           749:       arg1=arg1*arg
        !           750:       gmat(j,i)=arg1
        !           751:    30 continue
        !           752:    40 continue
        !           753: c
        !           754: c  solve for gear coefficients ag(*)
        !           755: c
        !           756: c
        !           757: c  lu decomposition
        !           758: c
        !           759:       do 70 i=2,istop
        !           760:       jstart=i+1
        !           761:       if (jstart.gt.istop) go to 70
        !           762:       do 60 j=jstart,istop
        !           763:       gmat(j,i)=gmat(j,i)/gmat(i,i)
        !           764:       do 50 k=jstart,istop
        !           765:       gmat(j,k)=gmat(j,k)-gmat(j,i)*gmat(i,k)
        !           766:    50 continue
        !           767:    60 continue
        !           768:    70 continue
        !           769: c
        !           770: c  forward substitution
        !           771: c
        !           772:       do 90 i=2,istop
        !           773:       jstart=i+1
        !           774:       if (jstart.gt.istop) go to 90
        !           775:       do 80 j=jstart,istop
        !           776:       ag(j)=ag(j)-gmat(j,i)*ag(i)
        !           777:    80 continue
        !           778:    90 continue
        !           779: c
        !           780: c  backward substitution
        !           781: c
        !           782:       ag(istop)=ag(istop)/gmat(istop,istop)
        !           783:       ir=istop
        !           784:       do 110 i=2,istop
        !           785:       jstart=ir
        !           786:       ir=ir-1
        !           787:       do 100 j=jstart,istop
        !           788:       ag(ir)=ag(ir)-gmat(ir,j)*ag(j)
        !           789:   100 continue
        !           790:       ag(ir)=ag(ir)/gmat(ir,ir)
        !           791:   110 continue
        !           792: c
        !           793: c  finished
        !           794: c
        !           795:   200 return
        !           796:       end
        !           797:       subroutine trunc(delnew)
        !           798:       implicit double precision (a-h,o-z)
        !           799: c
        !           800: c     this routine determines the new transient stepsize by either
        !           801: c calling terr to estimate the local truncation error, or by checking
        !           802: c on the number of iterations needed to converge at the last timepoint.
        !           803: c
        !           804:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !           805:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !           806:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !           807:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !           808:      2   itemno,nosolv,ipostp,iscrch
        !           809:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
        !           810:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
        !           811:       common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg
        !           812:       common /blank/ value(1000)
        !           813:       integer nodplc(64)
        !           814:       complex*16 cvalue(32)
        !           815:       equivalence (value(1),nodplc(1),cvalue(1))
        !           816: c
        !           817: c
        !           818:       if (lvltim.ne.0) go to 5
        !           819:       delnew=dmin1(tstep,delmax)
        !           820:       return
        !           821:     5 if (lvltim.ne.1) go to 10
        !           822:       delnew=delta
        !           823:       if (iterno.gt.itl3) return
        !           824:       delnew=dmin1(2.0d0*delta,tstep,delmax)
        !           825:       return
        !           826: c
        !           827: c  capacitors
        !           828: c
        !           829:    10 delnew=1.0d20
        !           830:       loc=locate(2)
        !           831:    20 if (loc.eq.0) go to 30
        !           832:       loct=nodplc(loc+8)
        !           833:       call terr(loct,delnew)
        !           834:       loc=nodplc(loc)
        !           835:       go to 20
        !           836: c
        !           837: c  inductors
        !           838: c
        !           839:    30 loc=locate(3)
        !           840:    40 if (loc.eq.0) go to 50
        !           841:       loct=nodplc(loc+11)
        !           842:       call terr(loct,delnew)
        !           843:       loc=nodplc(loc)
        !           844:       go to 40
        !           845: c
        !           846: c  diodes
        !           847: c
        !           848:    50 loc=locate(11)
        !           849:    60 if (loc.eq.0) go to 70
        !           850:       loct=nodplc(loc+11)
        !           851:       call terr(loct+3,delnew)
        !           852:       loc=nodplc(loc)
        !           853:       go to 60
        !           854: c
        !           855: c  bjts
        !           856: c
        !           857:    70 loc=locate(12)
        !           858:    80 if (loc.eq.0) go to 90
        !           859:       loct=nodplc(loc+22)
        !           860:       call terr(loct+8,delnew)
        !           861:       call terr(loct+10,delnew)
        !           862:       call terr(loct+12,delnew)
        !           863:       loc=nodplc(loc)
        !           864:       go to 80
        !           865: c
        !           866: c  jfets
        !           867: c
        !           868:    90 loc=locate(13)
        !           869:   100 if (loc.eq.0) go to 110
        !           870:       loct=nodplc(loc+19)
        !           871:       call terr(loct+9,delnew)
        !           872:       call terr(loct+11,delnew)
        !           873:       loc=nodplc(loc)
        !           874:       go to 100
        !           875: c
        !           876: c  mosfets
        !           877: c
        !           878:   110 loc=locate(14)
        !           879:   120 if (loc.eq.0) go to 200
        !           880:       loct=nodplc(loc+26)
        !           881:       call terr(loct+12,delnew)
        !           882:       call terr(loct+14,delnew)
        !           883:       call terr(loct+16,delnew)
        !           884:       call terr(loct+18,delnew)
        !           885:       call terr(loct+20,delnew)
        !           886:       loc=nodplc(loc)
        !           887:       go to 120
        !           888: c
        !           889: c  delta is allowed only to double at each timepoint
        !           890: c
        !           891:   200 delnew=dmin1(2.0d0*delta,delnew,delmax)
        !           892:       return
        !           893:       end
        !           894:       subroutine terr(loct,delnew)
        !           895:       implicit double precision (a-h,o-z)
        !           896: c
        !           897: c     this routine estimates the local truncation error for a particular
        !           898: c circuit element.  it then computes the appropriate stepsize which
        !           899: c should be used.
        !           900: c
        !           901:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !           902:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !           903:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !           904:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !           905:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !           906:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !           907:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !           908:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !           909:      2   itemno,nosolv,ipostp,iscrch
        !           910:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
        !           911:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
        !           912:       common /blank/ value(1000)
        !           913:       integer nodplc(64)
        !           914:       complex*16 cvalue(32)
        !           915:       equivalence (value(1),nodplc(1),cvalue(1))
        !           916: c
        !           917: c
        !           918:       dimension qcap(1),ccap(1),diff(8),deltmp(7),coef(6)
        !           919:       equivalence (qcap(1),value(1)),(ccap(1),value(2))
        !           920:       data coef / 5.000000000d-1, 2.222222222d-1, 1.363636364d-1,
        !           921:      1            9.600000000d-2, 7.299270073d-2, 5.830903790d-2 /
        !           922:       data xtwelv / 8.333333333d-2 /
        !           923: c
        !           924: c
        !           925:       tol=reltol*dmax1(dabs(ccap(lx0+loct)),dabs(ccap(lx1+loct)))+abstol
        !           926:       ctol=reltol*dmax1(dabs(qcap(lx0+loct)),dabs(qcap(lx1+loct)),
        !           927:      1   chgtol)/delta
        !           928:       tol=dmax1(tol,ctol)
        !           929: c
        !           930: c  determine divided differences
        !           931: c
        !           932:       go to (6,5,4,3,2,1), iord
        !           933:     1 diff(8)=qcap(lx7+loct)
        !           934:     2 diff(7)=qcap(lx6+loct)
        !           935:     3 diff(6)=qcap(lx5+loct)
        !           936:     4 diff(5)=qcap(lx4+loct)
        !           937:     5 diff(4)=qcap(lx3+loct)
        !           938:     6 diff(3)=qcap(lx2+loct)
        !           939:       diff(2)=qcap(lx1+loct)
        !           940:       diff(1)=qcap(lx0+loct)
        !           941:       istop=iord+1
        !           942:       do 10 i=1,istop
        !           943:       deltmp(i)=delold(i)
        !           944:    10 continue
        !           945:    20 do 30 i=1,istop
        !           946:       diff(i)=(diff(i)-diff(i+1))/deltmp(i)
        !           947:    30 continue
        !           948:       istop=istop-1
        !           949:       if (istop.eq.0) go to 100
        !           950:       do 40 i=1,istop
        !           951:       deltmp(i)=deltmp(i+1)+delold(i)
        !           952:    40 continue
        !           953:       go to 20
        !           954: c
        !           955: c  diff(1) contains divided difference
        !           956: c
        !           957:   100 const=coef(iord)
        !           958:       if ((method.eq.1).and.(iord.eq.2)) const=xtwelv
        !           959:       del=trtol*tol/dmax1(abstol,const*dabs(diff(1)))
        !           960:       if (iord.eq.1) go to 200
        !           961:       if (iord.ge.3) go to 150
        !           962:       del=dsqrt(del)
        !           963:       go to 200
        !           964:   150 del=dexp(dlog(del)/dfloat(iord))
        !           965:   200 delnew=dmin1(delnew,del)
        !           966:       return
        !           967:       end
        !           968:       subroutine sorupd
        !           969:       implicit double precision (a-h,o-z)
        !           970: c
        !           971: c     this routine updates the independent voltage and current sources
        !           972: c used in the circuit.  it also updates the ltd table (which contains
        !           973: c previous (delayed) values of the sources used to model transmission
        !           974: c lines).
        !           975: c
        !           976:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !           977:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !           978:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !           979:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !           980:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !           981:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !           982:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !           983:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !           984:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !           985:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !           986:      2   itemno,nosolv,ipostp,iscrch
        !           987:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
        !           988:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
        !           989:       common /blank/ value(1000)
        !           990:       integer nodplc(64)
        !           991:       complex*16 cvalue(32)
        !           992:       equivalence (value(1),nodplc(1),cvalue(1))
        !           993: c
        !           994: c
        !           995:       do 500 id=9,10
        !           996:       loc=locate(id)
        !           997:    10 if (loc.eq.0) go to 500
        !           998:       locv=nodplc(loc+1)
        !           999:       locp=nodplc(loc+5)
        !          1000:       itype=nodplc(loc+4)+1
        !          1001:       go to (490,100,200,300,400,450), itype
        !          1002: c
        !          1003: c  pulse source
        !          1004: c
        !          1005:   100 v1=value(locp+1)
        !          1006:       v2=value(locp+2)
        !          1007:       t1=value(locp+3)
        !          1008:       t2=value(locp+4)
        !          1009:       t3=value(locp+5)
        !          1010:       t4=value(locp+6)
        !          1011:       period=value(locp+7)
        !          1012:       time1=time
        !          1013:       if (time1.le.0.0d0) go to 160
        !          1014:   110 if (time1.lt.t1+period) go to 120
        !          1015:       time1=time1-period
        !          1016:       go to 110
        !          1017:   120 if (time1.lt.t4) go to 130
        !          1018:       value(locv+1)=v1
        !          1019:       go to 490
        !          1020:   130 if (time1.lt.t3) go to 140
        !          1021:       value(locv+1)=v2+(time1-t3)*(v1-v2)/(t4-t3)
        !          1022:       go to 490
        !          1023:   140 if (time1.lt.t2) go to 150
        !          1024:       value(locv+1)=v2
        !          1025:       go to 490
        !          1026:   150 if (time1.lt.t1) go to 160
        !          1027:       value(locv+1)=v1+(time1-t1)*(v2-v1)/(t2-t1)
        !          1028:       go to 490
        !          1029:   160 value(locv+1)=v1
        !          1030:       go to 490
        !          1031: c
        !          1032: c  sinusoidal source
        !          1033: c
        !          1034:   200 v1=value(locp+1)
        !          1035:       v2=value(locp+2)
        !          1036:       omeg=value(locp+3)
        !          1037:       t1=value(locp+4)
        !          1038:       theta=value(locp+5)
        !          1039:       time1=time-t1
        !          1040:       if (time1.gt.0.0d0) go to 210
        !          1041:       value(locv+1)=v1
        !          1042:       go to 490
        !          1043:   210 if (theta.ne.0.0d0) go to 220
        !          1044:       value(locv+1)=v1+v2*dsin(omeg*time1)
        !          1045:       go to 490
        !          1046:   220 value(locv+1)=v1+v2*dsin(omeg*time1)*dexp(-time1*theta)
        !          1047:       go to 490
        !          1048: c
        !          1049: c  exponential source
        !          1050: c
        !          1051:   300 v1=value(locp+1)
        !          1052:       v2=value(locp+2)
        !          1053:       t1=value(locp+3)
        !          1054:       tau1=value(locp+4)
        !          1055:       t2=value(locp+5)
        !          1056:       tau2=value(locp+6)
        !          1057:       time1=time
        !          1058:       if (time1.gt.t1) go to 310
        !          1059:       value(locv+1)=v1
        !          1060:       go to 490
        !          1061:   310 if (time1.gt.t2) go to 320
        !          1062:       value(locv+1)=v1+(v2-v1)*(1.0d0-dexp((t1-time1)/tau1))
        !          1063:       go to 490
        !          1064:   320 value(locv+1)=v1+(v2-v1)*(1.0d0-dexp((t1-time1)/tau1))
        !          1065:      1   +(v1-v2)*(1.0d0-dexp((t2-time1)/tau2))
        !          1066:       go to 490
        !          1067: c
        !          1068: c  piecewise-linear source
        !          1069: c
        !          1070:   400 t1=value(locp+1)
        !          1071:       v1=value(locp+2)
        !          1072:       t2=value(locp+3)
        !          1073:       v2=value(locp+4)
        !          1074:       iknt=4
        !          1075:   410 if (time.le.t2) go to 420
        !          1076:       t1=t2
        !          1077:       v1=v2
        !          1078:       t2=value(locp+iknt+1)
        !          1079:       v2=value(locp+iknt+2)
        !          1080:       iknt=iknt+2
        !          1081:       go to 410
        !          1082:   420 value(locv+1)=v1+((time-t1)/(t2-t1))*(v2-v1)
        !          1083:       go to 490
        !          1084: c
        !          1085: c  single-frequency fm
        !          1086: c
        !          1087:   450 v1=value(locp+1)
        !          1088:       v2=value(locp+2)
        !          1089:       omegc=value(locp+3)
        !          1090:       xmod=value(locp+4)
        !          1091:       omegs=value(locp+5)
        !          1092:       value(locv+1)=v1+v2*dsin(omegc*time+xmod*dsin(omegs*time))
        !          1093:   490 loc=nodplc(loc)
        !          1094:       go to 10
        !          1095:   500 continue
        !          1096: c
        !          1097: c  update transmission line sources
        !          1098: c
        !          1099:       if (jelcnt(17).eq.0) go to 1000
        !          1100:       if (mode.ne.2) go to 1000
        !          1101:       call sizmem(ltd,ltdsiz)
        !          1102:       numtd=ltdsiz/ntlin
        !          1103:       if (numtd.lt.3) go to 900
        !          1104:       loc=locate(17)
        !          1105:   610 if (loc.eq.0) go to 1000
        !          1106:       locv=nodplc(loc+1)
        !          1107:       td=value(locv+2)
        !          1108:       baktim=time-td
        !          1109:       if (baktim.lt.0.0d0) go to 640
        !          1110:       ltdptr=nodplc(loc+30)
        !          1111:       icntr=2
        !          1112:       l1=ltd
        !          1113:       l2=l1+ntlin
        !          1114:       l3=l2+ntlin
        !          1115:       t1=value(l1+1)
        !          1116:       t2=value(l2+1)
        !          1117:   620 t3=value(l3+1)
        !          1118:       icntr=icntr+1
        !          1119:       if (baktim.le.t3) go to 630
        !          1120:       if (icntr.eq.numtd) go to 900
        !          1121:       l1=l2
        !          1122:       l2=l3
        !          1123:       l3=l2+ntlin
        !          1124:       t1=t2
        !          1125:       t2=t3
        !          1126:       go to 620
        !          1127:   630 dt1t2=t1-t2
        !          1128:       dt1t3=t1-t3
        !          1129:       dt2t3=t2-t3
        !          1130:       tdnom1=1.0d0/(dt1t2*dt1t3)
        !          1131:       tdnom2=-1.0d0/(dt1t2*dt2t3)
        !          1132:       tdnom3=1.0d0/(dt2t3*dt1t3)
        !          1133:       dtt1=baktim-t1
        !          1134:       dtt2=baktim-t2
        !          1135:       dtt3=baktim-t3
        !          1136:       tfact1=dtt2*dtt3*tdnom1
        !          1137:       tfact2=dtt1*dtt3*tdnom2
        !          1138:       tfact3=dtt1*dtt2*tdnom3
        !          1139:       value(locv+3)=value(l1+ltdptr+0)*tfact1+value(l2+ltdptr+0)*tfact2
        !          1140:      1   +value(l3+ltdptr+0)*tfact3
        !          1141:       value(locv+4)=value(l1+ltdptr+1)*tfact1+value(l2+ltdptr+1)*tfact2
        !          1142:      1   +value(l3+ltdptr+1)*tfact3
        !          1143:   640 loc=nodplc(loc)
        !          1144:       go to 610
        !          1145: c
        !          1146: c  internal logic error:  less than 3 entries in ltd
        !          1147: c
        !          1148:   900 nogo=1
        !          1149:       write (6,901) numtd,icntr
        !          1150:   901 format('0*abort*:  internal spice error:  sorupd:  ',2i5/)
        !          1151: c
        !          1152: c  finished
        !          1153: c
        !          1154:  1000 return
        !          1155:       end
        !          1156:       subroutine iter8(itlim)
        !          1157:       implicit double precision (a-h,o-z)
        !          1158: c
        !          1159: c     this routine drives the newton-raphson iteration technique used to
        !          1160: c solve the set of nonlinear circuit equations.
        !          1161: c
        !          1162:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !          1163:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !          1164:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !          1165:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !          1166:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !          1167:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !          1168:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
        !          1169:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
        !          1170:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !          1171:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !          1172:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          1173:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          1174:      2   itemno,nosolv,ipostp,iscrch
        !          1175:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
        !          1176:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
        !          1177:       common /blank/ value(1000)
        !          1178:       integer nodplc(64)
        !          1179:       complex*16 cvalue(32)
        !          1180:       equivalence (value(1),nodplc(1),cvalue(1))
        !          1181: c
        !          1182: c
        !          1183:       igoof=0
        !          1184:       iterno=0
        !          1185:       ndrflo=0
        !          1186:       noncon=0
        !          1187:       ipass=0
        !          1188: c
        !          1189: c  construct linear equations and check convergence
        !          1190: c
        !          1191:    10 call load
        !          1192:       if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 300
        !          1193:       iterno=iterno+1
        !          1194:       go to (20,30,40,50,50,50),initf
        !          1195:    20 if(mode.ne.1) go to 22
        !          1196:       call sizmem(nsnod,nic)
        !          1197:       if(nic.eq.0) go to 22
        !          1198:       if(ipass.ne.0) noncon=ipass
        !          1199:       ipass=0
        !          1200:    22 if(noncon.eq.0) go to 300
        !          1201:       go to 100
        !          1202:    30 initf=3
        !          1203:    40 if (noncon.eq.0) initf=1
        !          1204:       ipass=1
        !          1205:       go to 100
        !          1206:    50 initf=1
        !          1207: c
        !          1208: c  solve equations for next iteration
        !          1209: c
        !          1210:   100 if (iterno.ge.itlim) go to 200
        !          1211:   110 call dcdcmp
        !          1212:       call dcsol
        !          1213:   120 if (igoof.eq.0) go to 130
        !          1214:       ndrflo=ndrflo+1
        !          1215:       igoof=0
        !          1216:   130 value(lvn+1)=0.0d0
        !          1217:       ntemp=noncon
        !          1218:       noncon=0
        !          1219:       if (ntemp.gt.0) go to 150
        !          1220:       if (iterno.eq.1) go to 150
        !          1221:       do 140 i=2,numnod
        !          1222:       vold=value(lvnim1+i)
        !          1223:       vnew=value(lvn+i)
        !          1224:       tol=reltol*dmax1(dabs(vold),dabs(vnew))+vntol
        !          1225:       if (dabs(vold-vnew).le.tol) go to 140
        !          1226:       noncon=noncon+1
        !          1227:   140 continue
        !          1228:   150 call copy8(value(lvn+1),value(lvnim1+1),nstop)
        !          1229:       go to 10
        !          1230: c
        !          1231: c  no convergence
        !          1232: c
        !          1233:   200 igoof=1
        !          1234:   300 if (ndrflo.eq.0) go to 400
        !          1235:       write (6,301) ndrflo
        !          1236:   301 format('0warning:  underflow occurred ',i4,' time(s)')
        !          1237: c
        !          1238: c  finished
        !          1239: c
        !          1240:   400 return
        !          1241:       end
        !          1242:       subroutine load
        !          1243:       implicit double precision (a-h,o-z)
        !          1244: c
        !          1245: c     this routine zeroes-out and then loads the coefficient matrix.
        !          1246: c the active devices and the controlled sources are loaded by separate
        !          1247: c subroutines.
        !          1248: c
        !          1249:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !          1250:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !          1251:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !          1252:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !          1253:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !          1254:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !          1255:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !          1256:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !          1257:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          1258:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          1259:      2   itemno,nosolv,ipostp,iscrch
        !          1260:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
        !          1261:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
        !          1262:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
        !          1263:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
        !          1264:       common /blank/ value(1000)
        !          1265:       integer nodplc(64)
        !          1266:       complex*16 cvalue(32)
        !          1267:       equivalence (value(1),nodplc(1),cvalue(1))
        !          1268: c
        !          1269: c
        !          1270:       dimension qcap(1),ccap(1)
        !          1271:       equivalence (qcap(1),value(1)),(ccap(1),value(2))
        !          1272:       dimension find(1),vind(1)
        !          1273:       equivalence (find(1),value(1)),(vind(1),value(2))
        !          1274: c
        !          1275: c  zero y matrix and current vector
        !          1276: c
        !          1277:       call zero8(value(lvn+1),nstop+nstop+nut+nlt)
        !          1278: c
        !          1279: c  resistors
        !          1280: c
        !          1281:       loc=locate(1)
        !          1282:    20 if (loc.eq.0) go to 30
        !          1283:       locv=nodplc(loc+1)
        !          1284:       val=value(locv+1)
        !          1285:       locy=lynl+nodplc(loc+6)
        !          1286:       value(locy)=value(locy)+val
        !          1287:       locy=lynl+nodplc(loc+7)
        !          1288:       value(locy)=value(locy)+val
        !          1289:       locy=lynl+nodplc(loc+4)
        !          1290:       value(locy)=value(locy)-val
        !          1291:       locy=lynl+nodplc(loc+5)
        !          1292:       value(locy)=value(locy)-val
        !          1293:       loc=nodplc(loc)
        !          1294:       go to 20
        !          1295: c
        !          1296: c  capacitors
        !          1297: c
        !          1298:    30 loc=locate(2)
        !          1299:    40 if (loc.eq.0) go to 100
        !          1300:       locv=nodplc(loc+1)
        !          1301:       node1=nodplc(loc+2)
        !          1302:       node2=nodplc(loc+3)
        !          1303:       lcoef=nodplc(loc+7)
        !          1304:       call sizmem(nodplc(loc+7),ncoef)
        !          1305:       loct=nodplc(loc+8)
        !          1306:       vcap=value(locv+2)
        !          1307:       if ((mode.eq.1).and.(initf.eq.2)) go to 45
        !          1308:       if ((nosolv.ne.0).and.(initf.eq.5)) go to 45
        !          1309:       vcap=value(lvnim1+node1)-value(lvnim1+node2)
        !          1310:    45 value(locv+3)=vcap
        !          1311:       if (mode.eq.1) go to 60
        !          1312:       if (initf.ne.6) go to 50
        !          1313:       qcap(lx0+loct)=qcap(lx1+loct)
        !          1314:       go to 60
        !          1315:    50 call evpoly(qcap(lx0+loct),-1,lcoef,ncoef,locv+2,1,loc+8)
        !          1316:       if (initf.ne.5) go to 60
        !          1317:       if (nosolv.eq.0) go to 55
        !          1318:       vcap=value(locv+2)
        !          1319:       value(locv+3)=vcap
        !          1320:       call evpoly(qcap(lx0+loct),-1,lcoef,ncoef,locv+2,1,loc+8)
        !          1321:    55 qcap(lx1+loct)=qcap(lx0+loct)
        !          1322:    60 call evpoly(value(locv+1),0,lcoef,ncoef,locv+2,1,loc+8)
        !          1323:       if (mode.eq.1) go to 90
        !          1324:       call intgr8(geq,ceq,value(locv+1),loct)
        !          1325:       ceq=ceq-geq*vcap+ag(1)*qcap(lx0+loct)
        !          1326:       if(initf.ne.5) go to 70
        !          1327:       ccap(lx1+loct)=ccap(lx0+loct)
        !          1328:    70 locy=lynl+nodplc(loc+10)
        !          1329:       value(locy)=value(locy)+geq
        !          1330:       locy=lynl+nodplc(loc+11)
        !          1331:       value(locy)=value(locy)+geq
        !          1332:       locy=lynl+nodplc(loc+5)
        !          1333:       value(locy)=value(locy)-geq
        !          1334:       locy=lynl+nodplc(loc+6)
        !          1335:       value(locy)=value(locy)-geq
        !          1336:       value(lvn+node1)=value(lvn+node1)-ceq
        !          1337:       value(lvn+node2)=value(lvn+node2)+ceq
        !          1338:    90 loc=nodplc(loc)
        !          1339:       go to 40
        !          1340: c
        !          1341: c  inductors
        !          1342: c
        !          1343:   100 if (jelcnt(3).eq.0) go to 400
        !          1344:       if (mode.eq.1) go to 150
        !          1345:       if (initf.eq.6) go to 150
        !          1346:       loc=locate(3)
        !          1347:   110 if (loc.eq.0) go to 120
        !          1348:       locv=nodplc(loc+1)
        !          1349:       iptr=nodplc(loc+5)
        !          1350:       lcoef=nodplc(loc+10)
        !          1351:       call sizmem(nodplc(loc+10),ncoef)
        !          1352:       loct=nodplc(loc+11)
        !          1353:       cind=value(lvnim1+iptr)
        !          1354:       if ((mode.eq.1).and.(initf.eq.2)) cind=value(locv+2)
        !          1355:       if ((initf.eq.5).and.(nosolv.ne.0)) cind=value(locv+2)
        !          1356:       value(locv+3)=cind
        !          1357:       call evpoly(find(lx0+loct),-1,lcoef,ncoef,locv+2,1,loc+11)
        !          1358:       loc=nodplc(loc)
        !          1359:       go to 110
        !          1360:   120 loc=locate(4)
        !          1361:   130 if (loc.eq.0) go to 150
        !          1362:       locv=nodplc(loc+1)
        !          1363:       nl1=nodplc(loc+2)
        !          1364:       nl2=nodplc(loc+3)
        !          1365:       iptr1=nodplc(nl1+5)
        !          1366:       iptr2=nodplc(nl2+5)
        !          1367:       loct1=nodplc(nl1+11)
        !          1368:       loct2=nodplc(nl2+11)
        !          1369:       find(lx0+loct1)=find(lx0+loct1)+value(locv+1)*value(lvnim1+iptr2)
        !          1370:       find(lx0+loct2)=find(lx0+loct2)+value(locv+1)*value(lvnim1+iptr1)
        !          1371:       loc=nodplc(loc)
        !          1372:       go to 130
        !          1373:   150 loc=locate(3)
        !          1374:   160 if (loc.eq.0) go to 300
        !          1375:       locv=nodplc(loc+1)
        !          1376:       iptr=nodplc(loc+5)
        !          1377:       lcoef=nodplc(loc+10)
        !          1378:       call sizmem(nodplc(loc+10),ncoef)
        !          1379:       loct=nodplc(loc+11)
        !          1380:       cind=value(lvnim1+iptr)
        !          1381:       if ((mode.eq.1).and.(initf.eq.2)) cind=value(locv+2)
        !          1382:       if ((nosolv.ne.0).and.(initf.eq.5)) cind=value(locv+2)
        !          1383:       value(locv+3)=cind
        !          1384:       if (mode.ne.1) go to 200
        !          1385:       veq=0.0d0
        !          1386:       req=0.0d0
        !          1387:       go to 210
        !          1388:   200 if (initf.ne.6) go to 205
        !          1389:       find(lx0+loct)=find(lx1+loct)
        !          1390:       go to 210
        !          1391:   205 if (initf.ne.5) go to 210
        !          1392:       find(lx1+loct)=find(lx0+loct)
        !          1393:   210 call evpoly(value(locv+1),0,lcoef,ncoef,locv+2,1,loc+11)
        !          1394:       if (mode.eq.1) go to 250
        !          1395:       call intgr8(req,veq,value(locv+1),loct)
        !          1396:       if (ncoef.eq.1) go to 250
        !          1397:       veq=veq-req*cind+ag(1)*find(lx0+loct)
        !          1398:   250 value(lvn+iptr)=veq
        !          1399:       if(initf.ne.5) go to 260
        !          1400:       vind(lx1+loct)=vind(lx0+loct)
        !          1401:   260 locy=lynl+nodplc(loc+13)
        !          1402:       value(locy)=-req
        !          1403:       locy=lynl+nodplc(loc+6)
        !          1404:       value(locy)=1.0d0
        !          1405:       locy=lynl+nodplc(loc+7)
        !          1406:       value(locy)=-1.0d0
        !          1407:       locy=lynl+nodplc(loc+8)
        !          1408:       value(locy)=1.0d0
        !          1409:       locy=lynl+nodplc(loc+9)
        !          1410:       value(locy)=-1.0d0
        !          1411:       loc=nodplc(loc)
        !          1412:       go to 160
        !          1413: c
        !          1414: c  mutual inductances
        !          1415: c
        !          1416:   300 loc=locate(4)
        !          1417:   310 if (loc.eq.0) go to 400
        !          1418:       locv=nodplc(loc+1)
        !          1419:       req=ag(1)*value(locv+1)
        !          1420:       locy=lynl+nodplc(loc+4)
        !          1421:       value(locy)=-req
        !          1422:       locy=lynl+nodplc(loc+5)
        !          1423:       value(locy)=-req
        !          1424:       loc=nodplc(loc)
        !          1425:       go to 310
        !          1426: c
        !          1427: c  nonlinear controlled sources
        !          1428: c
        !          1429:   400 call nlcsrc
        !          1430: c
        !          1431: c  voltage sources
        !          1432: c
        !          1433:       loc=locate(9)
        !          1434:   610 if (loc.eq.0) go to 700
        !          1435:       locv=nodplc(loc+1)
        !          1436:       iptr=nodplc(loc+6)
        !          1437:       value(lvn+iptr)=value(locv+1)
        !          1438:       locy=lynl+nodplc(loc+7)
        !          1439:       value(locy)=value(locy)+1.0d0
        !          1440:       locy=lynl+nodplc(loc+8)
        !          1441:       value(locy)=value(locy)-1.0d0
        !          1442:       locy=lynl+nodplc(loc+9)
        !          1443:       value(locy)=value(locy)+1.0d0
        !          1444:       locy=lynl+nodplc(loc+10)
        !          1445:       value(locy)=value(locy)-1.0d0
        !          1446:       loc=nodplc(loc)
        !          1447:       go to 610
        !          1448: c
        !          1449: c  current sources
        !          1450: c
        !          1451:   700 loc=locate(10)
        !          1452:   710 if (loc.eq.0) go to 800
        !          1453:       locv=nodplc(loc+1)
        !          1454:       node1=nodplc(loc+2)
        !          1455:       node2=nodplc(loc+3)
        !          1456:       val=value(locv+1)
        !          1457:       value(lvn+node1)=value(lvn+node1)-val
        !          1458:       value(lvn+node2)=value(lvn+node2)+val
        !          1459:       loc=nodplc(loc)
        !          1460:       go to 710
        !          1461: c
        !          1462: c  call device model routines
        !          1463: c
        !          1464:   800 call diode
        !          1465:       call bjt
        !          1466:       call jfet
        !          1467:       call mosfet
        !          1468: c
        !          1469: c  transmission lines
        !          1470: c
        !          1471:       loc=locate(17)
        !          1472:   910 if (loc.eq.0) go to 980
        !          1473:       locv=nodplc(loc+1)
        !          1474:       z0=value(locv+1)
        !          1475:       y0=1.0d0/z0
        !          1476:       node1=nodplc(loc+2)
        !          1477:       node2=nodplc(loc+3)
        !          1478:       node3=nodplc(loc+4)
        !          1479:       node4=nodplc(loc+5)
        !          1480:       ibr1=nodplc(loc+8)
        !          1481:       ibr2=nodplc(loc+9)
        !          1482:       locy=lynl+nodplc(loc+10)
        !          1483:       value(locy)=value(locy)+y0
        !          1484:       locy=lynl+nodplc(loc+11)
        !          1485:       value(locy)=-y0
        !          1486:       locy=lynl+nodplc(loc+12)
        !          1487:       value(locy)=-1.0d0
        !          1488:       locy=lynl+nodplc(loc+13)
        !          1489:       value(locy)=value(locy)+y0
        !          1490:       locy=lynl+nodplc(loc+14)
        !          1491:       value(locy)=-1.0d0
        !          1492:       locy=lynl+nodplc(loc+15)
        !          1493:       value(locy)=-y0
        !          1494:       locy=lynl+nodplc(loc+16)
        !          1495:       value(locy)=+y0
        !          1496:       locy=lynl+nodplc(loc+17)
        !          1497:       value(locy)=+1.0d0
        !          1498:       locy=lynl+nodplc(loc+18)
        !          1499:       value(locy)=+y0
        !          1500:       locy=lynl+nodplc(loc+19)
        !          1501:       value(locy)=+1.0d0
        !          1502:       locy=lynl+nodplc(loc+20)
        !          1503:       value(locy)=-1.0d0
        !          1504:       locy=lynl+nodplc(loc+23)
        !          1505:       value(locy)=+1.0d0
        !          1506:       locy=lynl+nodplc(loc+27)
        !          1507:       value(locy)=-1.0d0
        !          1508:       locy=lynl+nodplc(loc+28)
        !          1509:       value(locy)=+1.0d0
        !          1510:       locy=lynl+nodplc(loc+31)
        !          1511:       value(locy)=-y0
        !          1512:       locy=lynl+nodplc(loc+32)
        !          1513:       value(locy)=-y0
        !          1514:       if (mode.ne.1) go to 920
        !          1515:       locy=lynl+nodplc(loc+21)
        !          1516:       value(locy)=-1.0d0
        !          1517:       locy=lynl+nodplc(loc+22)
        !          1518:       value(locy)=+1.0d0
        !          1519:       locy=lynl+nodplc(loc+24)
        !          1520:       value(locy)=-(1.0d0-gmin)*z0
        !          1521:       locy=lynl+nodplc(loc+25)
        !          1522:       value(locy)=-1.0d0
        !          1523:       locy=lynl+nodplc(loc+26)
        !          1524:       value(locy)=+1.0d0
        !          1525:       locy=lynl+nodplc(loc+29)
        !          1526:       value(locy)=-(1.0d0-gmin)*z0
        !          1527:       go to 950
        !          1528:   920 if (initf.ne.5) go to 930
        !          1529:       if (nosolv.ne.0) go to 925
        !          1530:       value(locv+3)=value(lvnim1+node3)-value(lvnim1+node4)
        !          1531:      1   +value(lvnim1+ibr2)*z0
        !          1532:       value(locv+4)=value(lvnim1+node1)-value(lvnim1+node2)
        !          1533:      1   +value(lvnim1+ibr1)*z0
        !          1534:       go to 930
        !          1535:   925 value(locv+3)=value(locv+7)+value(locv+8)*z0
        !          1536:       value(locv+4)=value(locv+5)+value(locv+6)*z0
        !          1537:   930 value(lvn+ibr1)=value(locv+3)
        !          1538:       value(lvn+ibr2)=value(locv+4)
        !          1539:   950 loc=nodplc(loc)
        !          1540:       go to 910
        !          1541: c
        !          1542: c  initialize nodes
        !          1543: c
        !          1544:   980 if(mode.ne.1) go to 995
        !          1545:       if(initf.ne.3.and.initf.ne.2) go to 1000
        !          1546:       call sizmem(nsnod,nic)
        !          1547:       if(nic.eq.0) go to 995
        !          1548:       call sizmem(icnod,ntest)
        !          1549:       if(modedc.eq.2.and.ntest.ne.0) go to 995
        !          1550:       g=1.0d0
        !          1551:       do 990 i=1,nic
        !          1552:       locy=lynl+nodplc(nsmat+i)
        !          1553:       value(locy)=value(locy)+g
        !          1554:       node=nodplc(nsnod+i)
        !          1555:       value(lvn+node)=value(lvn+node)+value(nsval+i)*g
        !          1556:   990 continue
        !          1557: c
        !          1558: c  transient initial conditions (uic not specified)
        !          1559: c
        !          1560:   995 if(mode.ne.1) go to 1000
        !          1561:       if(modedc.ne.2) go to 1000
        !          1562:       if(nosolv.ne.0) go to 1000
        !          1563:       call sizmem(icnod,nic)
        !          1564:       if(nic.eq.0) go to 1000
        !          1565:       g=1.0d0
        !          1566:       do 996 i=1,nic
        !          1567:       locy=lynl+nodplc(icmat+i)
        !          1568:       value(locy)=value(locy)+g
        !          1569:       node=nodplc(icnod+i)
        !          1570:       value(lvn+node)=value(lvn+node)+value(icval+i)*g
        !          1571:   996 continue
        !          1572: c
        !          1573: c  reorder right-hand side
        !          1574: c
        !          1575:  1000 do 1010 i=2,nstop
        !          1576:       j=nodplc(iswap+i)
        !          1577:       value(ndiag+i)=value(lvn+j)
        !          1578:  1010 continue
        !          1579:       call copy8(value(ndiag+1),value(lvn+1),nstop)
        !          1580: c
        !          1581: c  finished
        !          1582: c
        !          1583:       return
        !          1584:       end
        !          1585:       subroutine nlcsrc
        !          1586:       implicit double precision (a-h,o-z)
        !          1587: c
        !          1588: c     this routine loads the nonlinear controlled sources into the
        !          1589: c coefficient matrix.
        !          1590: c
        !          1591:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !          1592:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !          1593:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !          1594:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !          1595:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !          1596:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !          1597:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !          1598:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !          1599:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          1600:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          1601:      2   itemno,nosolv,ipostp,iscrch
        !          1602:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
        !          1603:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
        !          1604:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
        !          1605:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
        !          1606:       common /blank/ value(1000)
        !          1607:       integer nodplc(64)
        !          1608:       complex*16 cvalue(32)
        !          1609:       equivalence (value(1),nodplc(1),cvalue(1))
        !          1610: c
        !          1611: c  nonlinear voltage-controlled current sources
        !          1612: c
        !          1613:       loc=locate(5)
        !          1614:    10 if (loc.eq.0) go to 100
        !          1615:       node1=nodplc(loc+2)
        !          1616:       node2=nodplc(loc+3)
        !          1617:       ndim=nodplc(loc+4)
        !          1618:       lnod=nodplc(loc+6)
        !          1619:       lmat=nodplc(loc+7)
        !          1620:       lcoef=nodplc(loc+8)
        !          1621:       call sizmem(nodplc(loc+8),ncoef)
        !          1622:       larg=nodplc(loc+9)
        !          1623:       lexp=nodplc(loc+10)
        !          1624:       lic=nodplc(loc+11)
        !          1625:       loct=nodplc(loc+12)+1
        !          1626:       icheck=0
        !          1627:       do 20 i=1,ndim
        !          1628:       call update(value(lic+i),loct,nodplc(lnod+1),nodplc(lnod+2),2,
        !          1629:      1   icheck)
        !          1630:       value(larg+i)=value(lx0+loct)
        !          1631:       loct=loct+2
        !          1632:       lnod=lnod+2
        !          1633:    20 continue
        !          1634:       call evpoly(cold,0,lcoef,ncoef,larg,ndim,lexp)
        !          1635:       loct=nodplc(loc+12)
        !          1636:       if (icheck.eq.1) go to 30
        !          1637:       if (initf.eq.6) go to 30
        !          1638:       tol=reltol*dmax1(dabs(cold),dabs(value(lx0+loct)))+abstol
        !          1639:       if (dabs(cold-value(lx0+loct)).lt.tol) go to 40
        !          1640:    30 noncon=noncon+1
        !          1641:    40 value(lx0+loct)=cold
        !          1642:       ceq=cold
        !          1643:       do 50 i=1,ndim
        !          1644:       call evpoly(geq,i,lcoef,ncoef,larg,ndim,lexp)
        !          1645:       loct=loct+2
        !          1646:       value(lx0+loct)=geq
        !          1647:       ceq=ceq-geq*value(larg+i)
        !          1648:       locy=lynl+nodplc(lmat+1)
        !          1649:       value(locy)=value(locy)+geq
        !          1650:       locy=lynl+nodplc(lmat+2)
        !          1651:       value(locy)=value(locy)-geq
        !          1652:       locy=lynl+nodplc(lmat+3)
        !          1653:       value(locy)=value(locy)-geq
        !          1654:       locy=lynl+nodplc(lmat+4)
        !          1655:       value(locy)=value(locy)+geq
        !          1656:       lmat=lmat+4
        !          1657:    50 continue
        !          1658:       value(lvn+node1)=value(lvn+node1)-ceq
        !          1659:       value(lvn+node2)=value(lvn+node2)+ceq
        !          1660:       loc=nodplc(loc)
        !          1661:       go to 10
        !          1662: c
        !          1663: c  nonlinear voltage controlled voltage sources
        !          1664: c
        !          1665:   100 loc=locate(6)
        !          1666:   110 if (loc.eq.0) go to 200
        !          1667:       node1=nodplc(loc+2)
        !          1668:       node2=nodplc(loc+3)
        !          1669:       ndim=nodplc(loc+4)
        !          1670:       iptr=nodplc(loc+6)
        !          1671:       lnod=nodplc(loc+7)
        !          1672:       lmat=nodplc(loc+8)
        !          1673:       lcoef=nodplc(loc+9)
        !          1674:       call sizmem(nodplc(loc+9),ncoef)
        !          1675:       larg=nodplc(loc+10)
        !          1676:       lexp=nodplc(loc+11)
        !          1677:       lic=nodplc(loc+12)
        !          1678:       loct=nodplc(loc+13)+2
        !          1679:       icheck=0
        !          1680:       do 120 i=1,ndim
        !          1681:       call update(value(lic+i),loct,nodplc(lnod+1),nodplc(lnod+2),2,
        !          1682:      1   icheck)
        !          1683:       value(larg+i)=value(lx0+loct)
        !          1684:       loct=loct+2
        !          1685:       lnod=lnod+2
        !          1686:   120 continue
        !          1687:       call evpoly(volt,0,lcoef,ncoef,larg,ndim,lexp)
        !          1688:       loct=nodplc(loc+13)
        !          1689:       if (icheck.eq.1) go to 130
        !          1690:       if (initf.eq.6) go to 130
        !          1691:       tol=reltol*dmax1(dabs(volt),dabs(value(lx0+loct)))+vntol
        !          1692:       if (dabs(volt-value(lx0+loct)).lt.tol) go to 140
        !          1693:   130 noncon=noncon+1
        !          1694:   140 value(lx0+loct)=volt
        !          1695:       value(lx0+loct+1)=value(lvnim1+iptr)
        !          1696:       veq=volt
        !          1697:       locy=lynl+nodplc(lmat+1)
        !          1698:       value(locy)=+1.0d0
        !          1699:       locy=lynl+nodplc(lmat+2)
        !          1700:       value(locy)=-1.0d0
        !          1701:       locy=lynl+nodplc(lmat+3)
        !          1702:       value(locy)=+1.0d0
        !          1703:       locy=lynl+nodplc(lmat+4)
        !          1704:       value(locy)=-1.0d0
        !          1705:       lmat=lmat+4
        !          1706:       loct=loct+1
        !          1707:       do 150 i=1,ndim
        !          1708:       call evpoly(vgain,i,lcoef,ncoef,larg,ndim,lexp)
        !          1709:       loct=loct+2
        !          1710:       value(lx0+loct)=vgain
        !          1711:       veq=veq-vgain*value(larg+i)
        !          1712:       locy=lynl+nodplc(lmat+1)
        !          1713:       value(locy)=value(locy)-vgain
        !          1714:       locy=lynl+nodplc(lmat+2)
        !          1715:       value(locy)=value(locy)+vgain
        !          1716:       lmat=lmat+2
        !          1717:   150 continue
        !          1718:       value(lvn+iptr)=veq
        !          1719:       loc=nodplc(loc)
        !          1720:       go to 110
        !          1721: c
        !          1722: c  nonlinear current-controlled current sources
        !          1723: c
        !          1724:   200 loc=locate(7)
        !          1725:   210 if (loc.eq.0) go to 300
        !          1726:       node1=nodplc(loc+2)
        !          1727:       node2=nodplc(loc+3)
        !          1728:       ndim=nodplc(loc+4)
        !          1729:       lvs=nodplc(loc+6)
        !          1730:       lmat=nodplc(loc+7)
        !          1731:       lcoef=nodplc(loc+8)
        !          1732:       call sizmem(nodplc(loc+8),ncoef)
        !          1733:       larg=nodplc(loc+9)
        !          1734:       lexp=nodplc(loc+10)
        !          1735:       lic=nodplc(loc+11)
        !          1736:       loct=nodplc(loc+12)+1
        !          1737:       icheck=0
        !          1738:       do 220 i=1,ndim
        !          1739:       iptr=nodplc(lvs+i)
        !          1740:       iptr=nodplc(iptr+6)
        !          1741:       call update(value(lic+i),loct,iptr,1,2,icheck)
        !          1742:       value(larg+i)=value(lx0+loct)
        !          1743:       loct=loct+2
        !          1744:   220 continue
        !          1745:       call evpoly(csrc,0,lcoef,ncoef,larg,ndim,lexp)
        !          1746:       loct=nodplc(loc+12)
        !          1747:       if (icheck.eq.1) go to 230
        !          1748:       if (initf.eq.6) go to 230
        !          1749:       tol=reltol*dmax1(dabs(csrc),dabs(value(lx0+loct)))+abstol
        !          1750:       if (dabs(csrc-value(lx0+loct)).lt.tol) go to 240
        !          1751:   230 noncon=noncon+1
        !          1752:   240 value(lx0+loct)=csrc
        !          1753:       ceq=csrc
        !          1754:       do 250 i=1,ndim
        !          1755:       call evpoly(cgain,i,lcoef,ncoef,larg,ndim,lexp)
        !          1756:       loct=loct+2
        !          1757:       value(lx0+loct)=cgain
        !          1758:       ceq=ceq-cgain*value(larg+i)
        !          1759:       locy=lynl+nodplc(lmat+1)
        !          1760:       value(locy)=value(locy)+cgain
        !          1761:       locy=lynl+nodplc(lmat+2)
        !          1762:       value(locy)=value(locy)-cgain
        !          1763:       lmat=lmat+2
        !          1764:   250 continue
        !          1765:       value(lvn+node1)=value(lvn+node1)-ceq
        !          1766:       value(lvn+node2)=value(lvn+node2)+ceq
        !          1767:       loc=nodplc(loc)
        !          1768:       go to 210
        !          1769: c
        !          1770: c  nonlinear current controlled voltage sources
        !          1771: c
        !          1772:   300 loc=locate(8)
        !          1773:   310 if (loc.eq.0) go to 1000
        !          1774:       node1=nodplc(loc+2)
        !          1775:       node2=nodplc(loc+3)
        !          1776:       ndim=nodplc(loc+4)
        !          1777:       ibr=nodplc(loc+6)
        !          1778:       lvs=nodplc(loc+7)
        !          1779:       lmat=nodplc(loc+8)
        !          1780:       lcoef=nodplc(loc+9)
        !          1781:       call sizmem(nodplc(loc+9),ncoef)
        !          1782:       larg=nodplc(loc+10)
        !          1783:       lexp=nodplc(loc+11)
        !          1784:       lic=nodplc(loc+12)
        !          1785:       loct=nodplc(loc+13)+2
        !          1786:       icheck=0
        !          1787:       do 320 i=1,ndim
        !          1788:       iptr=nodplc(lvs+i)
        !          1789:       iptr=nodplc(iptr+6)
        !          1790:       call update(value(lic+i),loct,iptr,1,2,icheck)
        !          1791:       value(larg+i)=value(lx0+loct)
        !          1792:       loct=loct+2
        !          1793:   320 continue
        !          1794:       call evpoly(volt,0,lcoef,ncoef,larg,ndim,lexp)
        !          1795:       loct=nodplc(loc+13)
        !          1796:       if (icheck.eq.1) go to 330
        !          1797:       if (initf.eq.6) go to 330
        !          1798:       tol=reltol*dmax1(dabs(volt),dabs(value(lx0+loct)))+vntol
        !          1799:       if (dabs(volt-value(lx0+loct)).lt.tol) go to 340
        !          1800:   330 noncon=noncon+1
        !          1801:   340 value(lx0+loct)=volt
        !          1802:       value(lx0+loct+1)=value(lvnim1+ibr)
        !          1803:       veq=volt
        !          1804:       locy=lynl+nodplc(lmat+1)
        !          1805:       value(locy)=+1.0d0
        !          1806:       locy=lynl+nodplc(lmat+2)
        !          1807:       value(locy)=-1.0d0
        !          1808:       locy=lynl+nodplc(lmat+3)
        !          1809:       value(locy)=+1.0d0
        !          1810:       locy=lynl+nodplc(lmat+4)
        !          1811:       value(locy)=-1.0d0
        !          1812:       lmat=lmat+4
        !          1813:       loct=loct+1
        !          1814:       do 350 i=1,ndim
        !          1815:       call evpoly(transr,i,lcoef,ncoef,larg,ndim,lexp)
        !          1816:       loct=loct+2
        !          1817:       value(lx0+loct)=transr
        !          1818:       veq=veq-transr*value(larg+i)
        !          1819:       locy=lynl+nodplc(lmat+i)
        !          1820:       value(locy)=value(locy)-transr
        !          1821:   350 continue
        !          1822:       value(lvn+ibr)=veq
        !          1823:       loc=nodplc(loc)
        !          1824:       go to 310
        !          1825: c
        !          1826: c  finished
        !          1827: c
        !          1828:  1000 return
        !          1829:       end
        !          1830:       subroutine update(vinit,loct,node1,node2,nupda,icheck)
        !          1831:       implicit double precision (a-h,o-z)
        !          1832: c
        !          1833: c     this routine updates and limits the controlling variables for the
        !          1834: c nonlinear controlled sources.
        !          1835: c
        !          1836:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !          1837:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !          1838:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !          1839:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !          1840:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !          1841:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !          1842:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          1843:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          1844:      2   itemno,nosolv,ipostp,iscrch
        !          1845:       common /blank/ value(1000)
        !          1846:       integer nodplc(64)
        !          1847:       complex*16 cvalue(32)
        !          1848:       equivalence (value(1),nodplc(1),cvalue(1))
        !          1849: c
        !          1850: c
        !          1851:       go to (40,10,40,20,30,50), initf
        !          1852:    10 vnew=vinit
        !          1853:       go to 70
        !          1854:    20 vnew=value(lx0+loct)
        !          1855:       go to 70
        !          1856:    30 vnew=value(lx1+loct)
        !          1857:       go to 70
        !          1858:    40 vnew=value(lvnim1+node1)-value(lvnim1+node2)
        !          1859:       go to 60
        !          1860:    50 call copy8(value(lx1+loct),value(lx0+loct),nupda)
        !          1861:       xfact=delta/delold(2)
        !          1862:       vnew=(1.0d0+xfact)*value(lx1+loct)-xfact*value(lx2+loct)
        !          1863:    60 if (dabs(vnew).le.1.0d0) go to 80
        !          1864:       delv=vnew-value(lx0+loct)
        !          1865:       if (dabs(delv).le.0.1d0) go to 80
        !          1866:       vlim=dmax1(dabs(0.1d0*value(lx0+loct)),0.1d0)
        !          1867:       vnew=value(lx0+loct)+dsign(dmin1(dabs(delv),vlim),delv)
        !          1868:       go to 70
        !          1869:    70 icheck=1
        !          1870:    80 value(lx0+loct)=vnew
        !          1871:       return
        !          1872:       end
        !          1873:       subroutine evpoly(result,itype,lcoef,ncoef,larg,
        !          1874:      1  narg,lexp)
        !          1875:       implicit double precision (a-h,o-z)
        !          1876: c
        !          1877: c     this routine evaluates a polynomial.  lcoef points to the coef-
        !          1878: c ficients, and larg points to the values of the polynomial argument(s).
        !          1879: c
        !          1880:       common /blank/ value(1000)
        !          1881:       integer nodplc(64)
        !          1882:       complex*16 cvalue(32)
        !          1883:       equivalence (value(1),nodplc(1),cvalue(1))
        !          1884: c
        !          1885: c
        !          1886:       if (itype) 100,200,300
        !          1887: c
        !          1888: c  integration (polynomial *must* be one-dimensional)
        !          1889: c
        !          1890:   100 result=0.0d0
        !          1891:       arg=1.0d0
        !          1892:       arg1=value(larg+1)
        !          1893:       do 110 i=1,ncoef
        !          1894:       arg=arg*arg1
        !          1895:       result=result+value(lcoef+i)*arg/dfloat(i)
        !          1896:   110 continue
        !          1897:       go to 1000
        !          1898: c
        !          1899: c  evaluation of the polynomial
        !          1900: c
        !          1901:   200 result=value(lcoef+1)
        !          1902:       if (ncoef.eq.1) go to 1000
        !          1903:       call zero4(nodplc(lexp+1),narg)
        !          1904:       do 220 i=2,ncoef
        !          1905:       call nxtpwr(nodplc(lexp+1),narg)
        !          1906:       if (value(lcoef+i).eq.0.0d0) go to 220
        !          1907:       arg=1.0d0
        !          1908:       do 210 j=1,narg
        !          1909:       call evterm(val,value(larg+j),nodplc(lexp+j))
        !          1910:       arg=arg*val
        !          1911:   210 continue
        !          1912:       result=result+value(lcoef+i)*arg
        !          1913:   220 continue
        !          1914:       go to 1000
        !          1915: c
        !          1916: c  partial derivative with respect to the itype*th variable
        !          1917: c
        !          1918:   300 result=0.0d0
        !          1919:       if (ncoef.eq.1) go to 1000
        !          1920:       call zero4(nodplc(lexp+1),narg)
        !          1921:       do 330 i=2,ncoef
        !          1922:       call nxtpwr(nodplc(lexp+1),narg)
        !          1923:       if (nodplc(lexp+itype).eq.0) go to 330
        !          1924:       if (value(lcoef+i).eq.0.0d0) go to 330
        !          1925:       arg=1.0d0
        !          1926:       do 320 j=1,narg
        !          1927:       if (j.eq.itype) go to 310
        !          1928:       call evterm(val,value(larg+j),nodplc(lexp+j))
        !          1929:       arg=arg*val
        !          1930:       go to 320
        !          1931:   310 call evterm(val,value(larg+j),nodplc(lexp+j)-1)
        !          1932:       arg=arg*dfloat(nodplc(lexp+j))*val
        !          1933:   320 continue
        !          1934:       result=result+value(lcoef+i)*arg
        !          1935:   330 continue
        !          1936: c
        !          1937: c  finished
        !          1938: c
        !          1939:  1000 return
        !          1940:       end
        !          1941:       subroutine evterm(val,arg,iexp)
        !          1942:       implicit double precision (a-h,o-z)
        !          1943: c
        !          1944: c     this routine evaluates one term of a polynomial.
        !          1945: c
        !          1946:       jexp=iexp+1
        !          1947:       if (jexp.ge.6) go to 60
        !          1948:       go to (10,20,30,40,50), jexp
        !          1949:    10 val=1.0d0
        !          1950:       go to 100
        !          1951:    20 val=arg
        !          1952:       go to 100
        !          1953:    30 val=arg*arg
        !          1954:       go to 100
        !          1955:    40 val=arg*arg*arg
        !          1956:       go to 100
        !          1957:    50 val=arg*arg
        !          1958:       val=val*val
        !          1959:       go to 100
        !          1960:    60 if (arg.eq.0.0d0) go to 70
        !          1961:       argexp=dfloat(iexp)*dlog(dabs(arg))
        !          1962:       if (argexp.lt.-200.0d0) go to 70
        !          1963:       val=dexp(argexp)
        !          1964:       if((iexp/2)*2.eq.iexp) go to 100
        !          1965:       val=dsign(val,arg)
        !          1966:       go to 100
        !          1967:    70 val=0.0d0
        !          1968: c
        !          1969: c  finished
        !          1970: c
        !          1971:   100 return
        !          1972:       end
        !          1973:       subroutine nxtpwr(pwrseq,pdim)
        !          1974:       implicit double precision (a-h,o-z)
        !          1975: c
        !          1976: c     this routine determines the 'next' set of exponents for the
        !          1977: c different dimensions of a polynomial.
        !          1978: c
        !          1979:       integer pwrseq(1),pdim,psum
        !          1980: c
        !          1981: c
        !          1982:       if (pdim.eq.1) go to 80
        !          1983:       k=pdim
        !          1984:    10 if (pwrseq(k).ne.0) go to 20
        !          1985:       k=k-1
        !          1986:       if (k.ne.0) go to 10
        !          1987:       go to 80
        !          1988:    20 if (k.eq.pdim) go to 30
        !          1989:       pwrseq(k)=pwrseq(k)-1
        !          1990:       pwrseq(k+1)=pwrseq(k+1)+1
        !          1991:       go to 100
        !          1992:    30 km1=k-1
        !          1993:       do 40 i=1,km1
        !          1994:       if (pwrseq(i).ne.0) go to 50
        !          1995:    40 continue
        !          1996:       pwrseq(1)=pwrseq(pdim)+1
        !          1997:       pwrseq(pdim)=0
        !          1998:       go to 100
        !          1999:    50 psum=1
        !          2000:       k=pdim
        !          2001:    60 if (pwrseq(k-1).ge.1) go to 70
        !          2002:       psum=psum+pwrseq(k)
        !          2003:       pwrseq(k)=0
        !          2004:       k=k-1
        !          2005:       go to 60
        !          2006:    70 pwrseq(k)=pwrseq(k)+psum
        !          2007:       pwrseq(k-1)=pwrseq(k-1)-1
        !          2008:       go to 100
        !          2009:    80 pwrseq(1)=pwrseq(1)+1
        !          2010: c
        !          2011: c  finished
        !          2012: c
        !          2013:   100 return
        !          2014:       end
        !          2015:       subroutine intgr8(geq,ceq,capval,loct)
        !          2016:       implicit double precision (a-h,o-z)
        !          2017: c
        !          2018: c     this routine performs the actual numerical integration for each
        !          2019: c circuit element.
        !          2020: c
        !          2021:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !          2022:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !          2023:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !          2024:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !          2025:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !          2026:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !          2027:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          2028:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          2029:      2   itemno,nosolv,ipostp,iscrch
        !          2030:       common /blank/ value(1000)
        !          2031:       integer nodplc(64)
        !          2032:       complex*16 cvalue(32)
        !          2033:       equivalence (value(1),nodplc(1),cvalue(1))
        !          2034: c
        !          2035: c
        !          2036:       dimension qcap(1),ccap(1)
        !          2037:       equivalence (qcap(1),value(1)),(ccap(1),value(2))
        !          2038: c
        !          2039: c
        !          2040:       if (method.eq.2) go to 100
        !          2041: c
        !          2042: c  trapezoidal algorithm
        !          2043: c
        !          2044:       if (iord.eq.1) go to 100
        !          2045:       ccap(lx0+loct)=-ccap(lx1+loct)*ag(2)
        !          2046:      1   +ag(1)*(qcap(lx0+loct)-qcap(lx1+loct))
        !          2047:       go to 190
        !          2048: c
        !          2049: c  gears algorithm
        !          2050: c
        !          2051:   100 go to (110,120,130,140,150,160), iord
        !          2052:   110 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
        !          2053:       go to 190
        !          2054:   120 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
        !          2055:      1              +ag(3)*qcap(lx2+loct)
        !          2056:       go to 190
        !          2057:   130 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
        !          2058:      1              +ag(3)*qcap(lx2+loct)+ag(4)*qcap(lx3+loct)
        !          2059:       go to 190
        !          2060:   140 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
        !          2061:      1              +ag(3)*qcap(lx2+loct)+ag(4)*qcap(lx3+loct)
        !          2062:      2              +ag(5)*qcap(lx4+loct)
        !          2063:       go to 190
        !          2064:   150 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
        !          2065:      1              +ag(3)*qcap(lx2+loct)+ag(4)*qcap(lx3+loct)
        !          2066:      2              +ag(5)*qcap(lx4+loct)+ag(6)*qcap(lx5+loct)
        !          2067:       go to 190
        !          2068:   160 ccap(lx0+loct)=ag(1)*qcap(lx0+loct)+ag(2)*qcap(lx1+loct)
        !          2069:      1              +ag(3)*qcap(lx2+loct)+ag(4)*qcap(lx3+loct)
        !          2070:      2              +ag(5)*qcap(lx4+loct)+ag(6)*qcap(lx5+loct)
        !          2071:      3              +ag(7)*qcap(lx6+loct)
        !          2072: c... ceq is the equivalent current applicable to linear capacitance
        !          2073: c    (inductance) only, i.e. q=c*v
        !          2074:   190 ceq=ccap(lx0+loct)-ag(1)*qcap(lx0+loct)
        !          2075:       geq=ag(1)*capval
        !          2076:       return
        !          2077:       end
        !          2078:       subroutine pnjlim(vnew,vold,vt,vcrit,icheck)
        !          2079:       implicit double precision (a-h,o-z)
        !          2080: c
        !          2081: c     this routine limits the change-per-iteration of device pn-junction
        !          2082: c voltages.
        !          2083: c
        !          2084:       if (vnew.le.vcrit) go to 30
        !          2085:       vlim=vt+vt
        !          2086:       delv=vnew-vold
        !          2087:       if (dabs(delv).le.vlim) go to 30
        !          2088:       if (vold.le.0.0d0) go to 20
        !          2089:       arg=1.0d0+delv/vt
        !          2090:       if (arg.le.0.0d0) go to 10
        !          2091:       vnew=vold+vt*dlog(arg)
        !          2092:       icheck=1
        !          2093:       go to 100
        !          2094:    10 vnew=vcrit
        !          2095:       icheck=1
        !          2096:       go to 100
        !          2097:    20 vnew=vt*dlog(vnew/vt)
        !          2098:       icheck=1
        !          2099:       go to 100
        !          2100:    30 continue
        !          2101: c
        !          2102: c  finished
        !          2103: c
        !          2104:   100 return
        !          2105:       end
        !          2106:       subroutine diode
        !          2107:       implicit double precision (a-h,o-z)
        !          2108: c
        !          2109: c     this routine processes diodes for dc and transient analyses.
        !          2110: c
        !          2111:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !          2112:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !          2113:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !          2114:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !          2115:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !          2116:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !          2117:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !          2118:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !          2119:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          2120:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          2121:      2   itemno,nosolv,ipostp,iscrch
        !          2122:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
        !          2123:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
        !          2124:       common /blank/ value(1000)
        !          2125:       integer nodplc(64)
        !          2126:       complex*16 cvalue(32)
        !          2127:       equivalence (value(1),nodplc(1),cvalue(1))
        !          2128: c
        !          2129: c
        !          2130:       dimension vdo(1),cdo(1),gdo(1),qd(1),cqd(1)
        !          2131:       equivalence (vdo(1),value(1)),(cdo(1),value(2)),
        !          2132:      1   (gdo(1),value(3)),(qd(1),value(4)),(cqd(1),value(5))
        !          2133: c
        !          2134: c
        !          2135:       loc=locate(11)
        !          2136:    10 if (loc.eq.0) return
        !          2137:       locv=nodplc(loc+1)
        !          2138:       node1=nodplc(loc+2)
        !          2139:       node2=nodplc(loc+3)
        !          2140:       node3=nodplc(loc+4)
        !          2141:       locm=nodplc(loc+5)
        !          2142:       ioff=nodplc(loc+6)
        !          2143:       locm=nodplc(locm+1)
        !          2144:       loct=nodplc(loc+11)
        !          2145: c
        !          2146: c  dc model parameters
        !          2147: c
        !          2148:       area=value(locv+1)
        !          2149:       csat=value(locm+1)*area
        !          2150:       gspr=value(locm+2)*area
        !          2151:       vte=value(locm+3)*vt
        !          2152:       bv=value(locm+13)
        !          2153:       vcrit=value(locm+18)
        !          2154: c
        !          2155: c  initialization
        !          2156: c
        !          2157:       icheck=1
        !          2158:       go to (100,20,30,50,60,70),initf
        !          2159:    20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25
        !          2160:       vd=value(locv+2)
        !          2161:       go to 300
        !          2162:    25 if(ioff.ne.0) go to 40
        !          2163:       vd=vcrit
        !          2164:       go to 300
        !          2165:    30 if (ioff.eq.0) go to 100
        !          2166:    40 vd=0.0d0
        !          2167:       go to 300
        !          2168:    50 vd=vdo(lx0+loct)
        !          2169:       go to 300
        !          2170:    60 vd=vdo(lx1+loct)
        !          2171:       go to 300
        !          2172:    70 xfact=delta/delold(2)
        !          2173:       vdo(lx0+loct)=vdo(lx1+loct)
        !          2174:       vd=(1.0d0+xfact)*vdo(lx1+loct)-xfact*vdo(lx2+loct)
        !          2175:       cdo(lx0+loct)=cdo(lx1+loct)
        !          2176:       gdo(lx0+loct)=gdo(lx1+loct)
        !          2177:       go to 110
        !          2178: c
        !          2179: c  compute new nonlinear branch voltage
        !          2180: c
        !          2181:   100 vd=value(lvnim1+node3)-value(lvnim1+node2)
        !          2182:   110 delvd=vd-vdo(lx0+loct)
        !          2183:       cdhat=cdo(lx0+loct)+gdo(lx0+loct)*delvd
        !          2184: c
        !          2185: c  bypass if solution has not changed
        !          2186: c
        !          2187:       if (6    .eq.6) go to 200
        !          2188:       tol=reltol*dmax1(dabs(vd),dabs(vdo(lx0+loct)))+vntol
        !          2189:       if (dabs(delvd).ge.tol) go to 200
        !          2190:       tol=reltol*dmax1(dabs(cdhat),dabs(cdo(lx0+loct)))+abstol
        !          2191:       if (dabs(cdhat-cdo(lx0+loct)).ge.tol) go to 200
        !          2192:       vd=vdo(lx0+loct)
        !          2193:       cd=cdo(lx0+loct)
        !          2194:       gd=gdo(lx0+loct)
        !          2195:       go to 800
        !          2196: c
        !          2197: c  limit new junction voltage
        !          2198: c
        !          2199:   200 vlim=vte+vte
        !          2200:       if(bv.eq.0.0d0) go to 205
        !          2201:       if (vd.lt.dmin1(0.0d0,-bv+10.0d0*vte)) go to 210
        !          2202:   205 icheck=0
        !          2203:       call pnjlim(vd,vdo(lx0+loct),vte,vcrit,icheck)
        !          2204:       go to 300
        !          2205:   210 vdtemp=-(vd+bv)
        !          2206:       call pnjlim(vdtemp,-(vdo(lx0+loct)+bv),vte,vcrit,icheck)
        !          2207:       vd=-(vdtemp+bv)
        !          2208: c
        !          2209: c  compute dc current and derivitives
        !          2210: c
        !          2211:   300 if (vd.lt.-5.0d0*vte) go to 310
        !          2212:       evd=dexp(vd/vte)
        !          2213:       cd=csat*(evd-1.0d0)+gmin*vd
        !          2214:       gd=csat*evd/vte+gmin
        !          2215:       go to 330
        !          2216:   310 if(bv.eq.0.0d0) go to 315
        !          2217:       if(vd.lt.-bv) go to 320
        !          2218:   315 gd=-csat/vd+gmin
        !          2219:       cd=gd*vd
        !          2220:       go to 330
        !          2221:   320 evrev=dexp(-(bv+vd)/vt)
        !          2222:       cd=-csat*(evrev-1.0d0+bv/vt)
        !          2223:       gd=csat*evrev/vt
        !          2224:   330 if (mode.ne.1) go to 500
        !          2225:       if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
        !          2226:       if (initf.eq.4) go to 500
        !          2227:       go to 700
        !          2228: c
        !          2229: c  charge storage elements
        !          2230: c
        !          2231:   500 tau=value(locm+4)
        !          2232:       czero=value(locm+5)*area
        !          2233:       pb=value(locm+6)
        !          2234:       xm=value(locm+7)
        !          2235:       fcpb=value(locm+12)
        !          2236:       if (vd.ge.fcpb) go to 510
        !          2237:       arg=1.0d0-vd/pb
        !          2238:       sarg=dexp(-xm*dlog(arg))
        !          2239:       qd(lx0+loct)=tau*cd+pb*czero*(1.0d0-arg*sarg)/(1.0d0-xm)
        !          2240:       capd=tau*gd+czero*sarg
        !          2241:       go to 520
        !          2242:   510 f1=value(locm+15)
        !          2243:       f2=value(locm+16)
        !          2244:       f3=value(locm+17)
        !          2245:       czof2=czero/f2
        !          2246:       qd(lx0+loct)=tau*cd+czero*f1+czof2*(f3*(vd-fcpb)
        !          2247:      1   +(xm/(pb+pb))*(vd*vd-fcpb*fcpb))
        !          2248:       capd=tau*gd+czof2*(f3+xm*vd/pb)
        !          2249: c
        !          2250: c  store small-signal parameters
        !          2251: c
        !          2252:   520 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700
        !          2253:       if (initf.ne.4) go to 600
        !          2254:       value(lx0+loct+4)=capd
        !          2255:       go to 1000
        !          2256: c
        !          2257: c  transient analysis
        !          2258: c
        !          2259:   600 if (initf.ne.5) go to 610
        !          2260:       qd(lx1+loct)=qd(lx0+loct)
        !          2261:   610 call intgr8(geq,ceq,capd,loct+3)
        !          2262:       gd=gd+geq
        !          2263:       cd=cd+cqd(lx0+loct)
        !          2264:       if (initf.ne.5) go to 700
        !          2265:       cqd(lx1+loct)=cqd(lx0+loct)
        !          2266: c
        !          2267: c  check convergence
        !          2268: c
        !          2269:   700 if (initf.ne.3) go to 710
        !          2270:       if (ioff.eq.0) go to 710
        !          2271:       go to 750
        !          2272:   710 if (icheck.eq.1) go to 720
        !          2273:       tol=reltol*dmax1(dabs(cdhat),dabs(cd))+abstol
        !          2274:       if (dabs(cdhat-cd).le.tol) go to 750
        !          2275:   720 noncon=noncon+1
        !          2276:   750 vdo(lx0+loct)=vd
        !          2277:       cdo(lx0+loct)=cd
        !          2278:       gdo(lx0+loct)=gd
        !          2279: c
        !          2280: c  load current vector
        !          2281: c
        !          2282:   800 cdeq=cd-gd*vd
        !          2283:       value(lvn+node2)=value(lvn+node2)+cdeq
        !          2284:       value(lvn+node3)=value(lvn+node3)-cdeq
        !          2285: c
        !          2286: c  load matrix
        !          2287: c
        !          2288:       locy=lynl+nodplc(loc+13)
        !          2289:       value(locy)=value(locy)+gspr
        !          2290:       locy=lynl+nodplc(loc+14)
        !          2291:       value(locy)=value(locy)+gd
        !          2292:       locy=lynl+nodplc(loc+15)
        !          2293:       value(locy)=value(locy)+gd+gspr
        !          2294:       locy=lynl+nodplc(loc+7)
        !          2295:       value(locy)=value(locy)-gspr
        !          2296:       locy=lynl+nodplc(loc+8)
        !          2297:       value(locy)=value(locy)-gd
        !          2298:       locy=lynl+nodplc(loc+9)
        !          2299:       value(locy)=value(locy)-gspr
        !          2300:       locy=lynl+nodplc(loc+10)
        !          2301:       value(locy)=value(locy)-gd
        !          2302:  1000 loc=nodplc(loc)
        !          2303:       go to 10
        !          2304:       end
        !          2305:       subroutine bjt
        !          2306:       implicit double precision (a-h,o-z)
        !          2307: c
        !          2308: c     this routine processes bjts for dc and transient analyses.
        !          2309: c
        !          2310:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !          2311:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !          2312:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !          2313:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !          2314:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !          2315:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !          2316:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !          2317:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !          2318:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          2319:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          2320:      2   itemno,nosolv,ipostp,iscrch
        !          2321:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
        !          2322:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
        !          2323:       common /blank/ value(1000)
        !          2324:       integer nodplc(64)
        !          2325:       complex*16 cvalue(32)
        !          2326:       equivalence (value(1),nodplc(1),cvalue(1))
        !          2327: c
        !          2328: c
        !          2329:       dimension vbeo(1),vbco(1),cco(1),cbo(1),gpio(1),gmuo(1),gmo(1),
        !          2330:      1   goo(1),qbe(1),cqbe(1),qbc(1),cqbc(1),qcs(1),cqcs(1),qbx(1),
        !          2331:      2   cqbx(1),gxo(1),cexbc(1)
        !          2332:       equivalence (vbeo(1),value(1)),(vbco(1),value(2)),
        !          2333:      1   (cco(1),value(3)),(cbo(1),value(4)),(gpio(1),value(5)),
        !          2334:      2   (gmuo(1),value(6)),(gmo(1),value(7)),(goo(1),value(8)),
        !          2335:      3   (qbe(1),value(9)),(cqbe(1),value(10)),(qbc(1),value(11)),
        !          2336:      4   (cqbc(1),value(12)),(qcs(1),value(13)),(cqcs(1),value(14)),
        !          2337:      5   (qbx(1),value(15)),(cqbx(1),value(16)),(gxo(1),value(17)),
        !          2338:      6   (cexbc(1),value(18))
        !          2339: c
        !          2340: c
        !          2341:       loc=locate(12)
        !          2342:    10 if (loc.eq.0) return
        !          2343:       locv=nodplc(loc+1)
        !          2344:       node1=nodplc(loc+2)
        !          2345:       node2=nodplc(loc+3)
        !          2346:       node3=nodplc(loc+4)
        !          2347:       node4=nodplc(loc+5)
        !          2348:       node5=nodplc(loc+6)
        !          2349:       node6=nodplc(loc+7)
        !          2350:       node7=nodplc(loc+30)
        !          2351:       locm=nodplc(loc+8)
        !          2352:       ioff=nodplc(loc+9)
        !          2353:       type=nodplc(locm+2)
        !          2354:       locm=nodplc(locm+1)
        !          2355:       loct=nodplc(loc+22)
        !          2356:       gccs=0.0d0
        !          2357:       ceqcs=0.0d0
        !          2358:       geqbx=0.0d0
        !          2359:       ceqbx=0.0d0
        !          2360:       geqcb=0.0d0
        !          2361: c
        !          2362: c  dc model paramters
        !          2363: c
        !          2364:       area=value(locv+1)
        !          2365:       bfm=value(locm+2)
        !          2366:       brm=value(locm+8)
        !          2367:       csat=value(locm+1)*area
        !          2368:       rbpr=value(locm+18)/area
        !          2369:       rbpi=value(locm+16)/area-rbpr
        !          2370:       gcpr=value(locm+20)*area
        !          2371:       gepr=value(locm+19)*area
        !          2372:       ova=value(locm+4)
        !          2373:       ovb=value(locm+10)
        !          2374:       oik=value(locm+5)/area
        !          2375:       c2=value(locm+6)*area
        !          2376:       vte=value(locm+7)*vt
        !          2377:       oikr=value(locm+11)/area
        !          2378:       c4=value(locm+12)*area
        !          2379:       vtc=value(locm+13)*vt
        !          2380:       vcrit=value(locm+54)
        !          2381:       td=value(locm+28)
        !          2382:       xjrb=value(locm+17)*area
        !          2383: c
        !          2384: c  initialization
        !          2385: c
        !          2386:       icheck=1
        !          2387:       go to (100,20,30,50,60,70),initf
        !          2388:    20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25
        !          2389:       vbe=type*value(locv+2)
        !          2390:       vce=type*value(locv+3)
        !          2391:       vbc=vbe-vce
        !          2392:       vbx=vbc
        !          2393:       vcs=0.0d0
        !          2394:       go to 300
        !          2395:    25 if(ioff.ne.0) go to 40
        !          2396:       vbe=vcrit
        !          2397:       vbc=0.0d0
        !          2398:       go to 300
        !          2399:    30 if (ioff.eq.0) go to 100
        !          2400:    40 vbe=0.0d0
        !          2401:       vbc=0.0d0
        !          2402:       go to 300
        !          2403:    50 vbe=vbeo(lx0+loct)
        !          2404:       vbc=vbco(lx0+loct)
        !          2405:       vbx=type*(value(lvnim1+node2)-value(lvnim1+node4))
        !          2406:       vcs=type*(value(lvnim1+node7)-value(lvnim1+node4))
        !          2407:       go to 300
        !          2408:    60 vbe=vbeo(lx1+loct)
        !          2409:       vbc=vbco(lx1+loct)
        !          2410:       vbx=type*(value(lvnim1+node2)-value(lvnim1+node4))
        !          2411:       vcs=type*(value(lvnim1+node7)-value(lvnim1+node4))
        !          2412:       if(mode.ne.2.or.nosolv.eq.0) go to 300
        !          2413:       vbx=type*(value(locv+2)-value(locv+3))
        !          2414:       vcs=0.0d0
        !          2415:       go to 300
        !          2416:    70 xfact=delta/delold(2)
        !          2417:       vbeo(lx0+loct)=vbeo(lx1+loct)
        !          2418:       vbe=(1.0d0+xfact)*vbeo(lx1+loct)-xfact*vbeo(lx2+loct)
        !          2419:       vbco(lx0+loct)=vbco(lx1+loct)
        !          2420:       vbc=(1.0d0+xfact)*vbco(lx1+loct)-xfact*vbco(lx2+loct)
        !          2421:       cco(lx0+loct)=cco(lx1+loct)
        !          2422:       cbo(lx0+loct)=cbo(lx1+loct)
        !          2423:       gpio(lx0+loct)=gpio(lx1+loct)
        !          2424:       gmuo(lx0+loct)=gmuo(lx1+loct)
        !          2425:       gmo(lx0+loct)=gmo(lx1+loct)
        !          2426:       goo(lx0+loct)=goo(lx1+loct)
        !          2427:       go to 110
        !          2428: c
        !          2429: c  compute new nonlinear branch voltages
        !          2430: c
        !          2431:   100 vbe=type*(value(lvnim1+node5)-value(lvnim1+node6))
        !          2432:       vbc=type*(value(lvnim1+node5)-value(lvnim1+node4))
        !          2433:   110 delvbe=vbe-vbeo(lx0+loct)
        !          2434:       delvbc=vbc-vbco(lx0+loct)
        !          2435:       vbx=type*(value(lvnim1+node2)-value(lvnim1+node4))
        !          2436:       vcs=type*(value(lvnim1+node7)-value(lvnim1+node4))
        !          2437:       cchat=cco(lx0+loct)+(gmo(lx0+loct)+goo(lx0+loct))*delvbe
        !          2438:      1   -(goo(lx0+loct)+gmuo(lx0+loct))*delvbc
        !          2439:       cbhat=cbo(lx0+loct)+gpio(lx0+loct)*delvbe+gmuo(lx0+loct)*delvbc
        !          2440: c
        !          2441: c  limit nonlinear branch voltages
        !          2442: c
        !          2443:   200 icheck=0
        !          2444:       call pnjlim(vbe,vbeo(lx0+loct),vt,vcrit,icheck)
        !          2445:       call pnjlim(vbc,vbco(lx0+loct),vt,vcrit,icheck)
        !          2446: c
        !          2447: c  determine dc current and derivitives
        !          2448: c
        !          2449:   300 vtn=vt*value(locm+3)
        !          2450:       if(vbe.le.-5.0d0*vtn) go to 320
        !          2451:       evbe=dexp(vbe/vtn)
        !          2452:       cbe=csat*(evbe-1.0d0)+gmin*vbe
        !          2453:       gbe=csat*evbe/vtn+gmin
        !          2454:       if (c2.ne.0.0d0) go to 310
        !          2455:       cben=0.0d0
        !          2456:       gben=0.0d0
        !          2457:       go to 350
        !          2458:   310 evben=dexp(vbe/vte)
        !          2459:       cben=c2*(evben-1.0d0)
        !          2460:       gben=c2*evben/vte
        !          2461:       go to 350
        !          2462:   320 gbe=-csat/vbe+gmin
        !          2463:       cbe=gbe*vbe
        !          2464:       gben=-c2/vbe
        !          2465:       cben=gben*vbe
        !          2466:   350 vtn=vt*value(locm+9)
        !          2467:       if(vbc.le.-5.0d0*vtn) go to 370
        !          2468:       evbc=dexp(vbc/vtn)
        !          2469:       cbc=csat*(evbc-1.0d0)+gmin*vbc
        !          2470:       gbc=csat*evbc/vtn+gmin
        !          2471:       if (c4.ne.0.0d0) go to 360
        !          2472:       cbcn=0.0d0
        !          2473:       gbcn=0.0d0
        !          2474:       go to 400
        !          2475:   360 evbcn=dexp(vbc/vtc)
        !          2476:       cbcn=c4*(evbcn-1.0d0)
        !          2477:       gbcn=c4*evbcn/vtc
        !          2478:       go to 400
        !          2479:   370 gbc=-csat/vbc+gmin
        !          2480:       cbc=gbc*vbc
        !          2481:       gbcn=-c4/vbc
        !          2482:       cbcn=gbcn*vbc
        !          2483: c
        !          2484: c  determine base charge terms
        !          2485: c
        !          2486:   400 q1=1.0d0/(1.0d0-ova*vbc-ovb*vbe)
        !          2487:       if (oik.ne.0.0d0) go to 405
        !          2488:       if (oikr.ne.0.0d0) go to 405
        !          2489:       qb=q1
        !          2490:       dqbdve=q1*qb*ovb
        !          2491:       dqbdvc=q1*qb*ova
        !          2492:       go to 410
        !          2493:   405 q2=oik*cbe+oikr*cbc
        !          2494:       arg=dmax1(0.0d0,1.0d0+4.0d0*q2)
        !          2495:       sqarg=1.0d0
        !          2496:       if(arg.ne.0.0d0) sqarg=dsqrt(arg)
        !          2497:       qb=q1*(1.0d0+sqarg)/2.0d0
        !          2498:       dqbdve=q1*(qb*ovb+oik*gbe/sqarg)
        !          2499:       dqbdvc=q1*(qb*ova+oikr*gbc/sqarg)
        !          2500: c
        !          2501: c  weil's approx. for excess phase applied with backward-
        !          2502: c  euler integration
        !          2503: c
        !          2504:   410 cc=0.0d0
        !          2505:       cex=cbe
        !          2506:       gex=gbe
        !          2507:       if(mode.eq.1) go to 420
        !          2508:       if(td.eq.0.0d0) go to 420
        !          2509:       arg1=delta/td
        !          2510:       arg2=3.0d0*arg1
        !          2511:       arg1=arg2*arg1
        !          2512:       denom=1.0d0+arg1+arg2
        !          2513:       arg3=arg1/denom
        !          2514:       if(initf.ne.5) go to 411
        !          2515:       cexbc(lx1+loct)=cbe/qb
        !          2516:       cexbc(lx2+loct)=cexbc(lx1+loct)
        !          2517:   411 cc=(cexbc(lx1+loct)*(1.0d0+delta/delold(2)+arg2)
        !          2518:      1  -cexbc(lx2+loct)*delta/delold(2))/denom
        !          2519:       cex=cbe*arg3
        !          2520:       gex=gbe*arg3
        !          2521:       cexbc(lx0+loct)=cc+cex/qb
        !          2522: c
        !          2523: c  determine dc incremental conductances
        !          2524: c
        !          2525:   420 cc=cc+(cex-cbc)/qb-cbc/brm-cbcn
        !          2526:       cb=cbe/bfm+cben+cbc/brm+cbcn
        !          2527:       gx=rbpr+rbpi/qb
        !          2528:       if(xjrb.eq.0.0d0) go to 430
        !          2529:       arg1=dmax1(cb/xjrb,1.0d-9)
        !          2530:       arg2=(-1.0d0+dsqrt(1.0d0+14.59025d0*arg1))/2.4317d0/dsqrt(arg1)
        !          2531:       arg1=dtan(arg2)
        !          2532:       gx=rbpr+3.0d0*rbpi*(arg1-arg2)/arg2/arg1/arg1
        !          2533:   430 if(gx.ne.0.0d0) gx=1.0d0/gx
        !          2534:       gpi=gbe/bfm+gben
        !          2535:       gmu=gbc/brm+gbcn
        !          2536:       go=(gbc+(cex-cbc)*dqbdvc/qb)/qb
        !          2537:       gm=(gex-(cex-cbc)*dqbdve/qb)/qb-go
        !          2538:       if (mode.ne.1) go to 500
        !          2539:       if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
        !          2540:       if (initf.eq.4) go to 500
        !          2541:       go to 700
        !          2542: c
        !          2543: c  charge storage elements
        !          2544: c
        !          2545:   500 tf=value(locm+24)
        !          2546:       tr=value(locm+33)
        !          2547:       czbe=value(locm+21)*area
        !          2548:       pe=value(locm+22)
        !          2549:       xme=value(locm+23)
        !          2550:       cdis=value(locm+32)
        !          2551:       ctot=value(locm+29)*area
        !          2552:       czbc=ctot*cdis
        !          2553:       czbx=ctot-czbc
        !          2554:       pc=value(locm+30)
        !          2555:       xmc=value(locm+31)
        !          2556:       fcpe=value(locm+46)
        !          2557:       czcs=value(locm+38)*area
        !          2558:       ps=value(locm+39)
        !          2559:       xms=value(locm+40)
        !          2560:       xtf=value(locm+25)
        !          2561:       ovtf=value(locm+26)
        !          2562:       xjtf=value(locm+27)*area
        !          2563:       if(tf.eq.0.0d0) go to 505
        !          2564:       if(vbe.le.0.0d0) go to 505
        !          2565:       argtf=0.0d0
        !          2566:       arg2=0.0d0
        !          2567:       arg3=0.0d0
        !          2568:       if(xtf.eq.0.0d0) go to 504
        !          2569:       argtf=xtf
        !          2570:       if(ovtf.ne.0.0d0) argtf=argtf*dexp(vbc*ovtf)
        !          2571:       arg2=argtf
        !          2572:       if(xjtf.eq.0.0d0) go to 503
        !          2573:       temp=cbe/(cbe+xjtf)
        !          2574:       argtf=argtf*temp*temp
        !          2575:       arg2=argtf*(3.0d0-temp-temp)
        !          2576:   503 arg3=cbe*argtf*ovtf
        !          2577:   504 cbe=cbe*(1.0d0+argtf)/qb
        !          2578:       gbe=(gbe*(1.0d0+arg2)-cbe*dqbdve)/qb
        !          2579:       geqcb=tf*(arg3-cbe*dqbdvc)/qb
        !          2580:   505 if (vbe.ge.fcpe) go to 510
        !          2581:       arg=1.0d0-vbe/pe
        !          2582:       sarg=dexp(-xme*dlog(arg))
        !          2583:       qbe(lx0+loct)=tf*cbe+pe*czbe*(1.0d0-arg*sarg)/(1.0d0-xme)
        !          2584:       capbe=tf*gbe+czbe*sarg
        !          2585:       go to 520
        !          2586:   510 f1=value(locm+47)
        !          2587:       f2=value(locm+48)
        !          2588:       f3=value(locm+49)
        !          2589:       czbef2=czbe/f2
        !          2590:       qbe(lx0+loct)=tf*cbe+czbe*f1+czbef2*(f3*(vbe-fcpe)
        !          2591:      1   +(xme/(pe+pe))*(vbe*vbe-fcpe*fcpe))
        !          2592:       capbe=tf*gbe+czbef2*(f3+xme*vbe/pe)
        !          2593:   520 fcpc=value(locm+50)
        !          2594:       f1=value(locm+51)
        !          2595:       f2=value(locm+52)
        !          2596:       f3=value(locm+53)
        !          2597:       if (vbc.ge.fcpc) go to 530
        !          2598:       arg=1.0d0-vbc/pc
        !          2599:       sarg=dexp(-xmc*dlog(arg))
        !          2600:       qbc(lx0+loct)=tr*cbc+pc*czbc*(1.0d0-arg*sarg)/(1.0d0-xmc)
        !          2601:       capbc=tr*gbc+czbc*sarg
        !          2602:       go to 540
        !          2603:   530 czbcf2=czbc/f2
        !          2604:       qbc(lx0+loct)=tr*cbc+czbc*f1+czbcf2*(f3*(vbc-fcpc)
        !          2605:      1   +(xmc/(pc+pc))*(vbc*vbc-fcpc*fcpc))
        !          2606:       capbc=tr*gbc+czbcf2*(f3+xmc*vbc/pc)
        !          2607:   540 if(vbx.ge.fcpc) go to 550
        !          2608:       arg=1.0d0-vbx/pc
        !          2609:       sarg=dexp(-xmc*dlog(arg))
        !          2610:       qbx(lx0+loct)=pc*czbx*(1.0d0-arg*sarg)/(1.0d0-xmc)
        !          2611:       capbx=czbx*sarg
        !          2612:       go to 560
        !          2613:   550 czbxf2=czbx/f2
        !          2614:       qbx(lx0+loct)=czbx*f1+czbxf2*(f3*(vbx-fcpc)+(xmc/(pc+pc))*
        !          2615:      1   (vbx*vbx-fcpc*fcpc))
        !          2616:       capbx=czbxf2*(f3+xmc*vbx/pc)
        !          2617:   560 if(vcs.ge.0.0d0) go to 570
        !          2618:       arg=1.0d0-vcs/ps
        !          2619:       sarg=dexp(-xms*dlog(arg))
        !          2620:       qcs(lx0+loct)=ps*czcs*(1.0d0-arg*sarg)/(1.0d0-xms)
        !          2621:       capcs=czcs*sarg
        !          2622:       go to 580
        !          2623:   570 qcs(lx0+loct)=vcs*czcs*(1.0d0+xms*vcs/(2.0d0*ps))
        !          2624:       capcs=czcs*(1.0d0+xms*vcs/ps)
        !          2625: c
        !          2626: c  store small-signal parameters
        !          2627: c
        !          2628:   580 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700
        !          2629:       if (initf.ne.4) go to 600
        !          2630:       value(lx0+loct+9)=capbe
        !          2631:       value(lx0+loct+11)=capbc
        !          2632:       value(lx0+loct+13)=capcs
        !          2633:       value(lx0+loct+15)=capbx
        !          2634:       value(lx0+loct+17)=geqcb
        !          2635:       go to 1000
        !          2636: c
        !          2637: c  transient analysis
        !          2638: c
        !          2639:   600 if (initf.ne.5) go to 610
        !          2640:       qbe(lx1+loct)=qbe(lx0+loct)
        !          2641:       qbc(lx1+loct)=qbc(lx0+loct)
        !          2642:       qbx(lx1+loct)=qbx(lx0+loct)
        !          2643:       qcs(lx1+loct)=qcs(lx0+loct)
        !          2644:   610 call intgr8(geq,ceq,capbe,loct+8)
        !          2645:       geqcb=geqcb*ag(1)
        !          2646:       gpi=gpi+geq
        !          2647:       cb=cb+cqbe(lx0+loct)
        !          2648:       call intgr8(geq,ceq,capbc,loct+10)
        !          2649:       gmu=gmu+geq
        !          2650:       cb=cb+cqbc(lx0+loct)
        !          2651:       cc=cc-cqbc(lx0+loct)
        !          2652:       call intgr8(gccs,ceq,capcs,loct+12)
        !          2653:       ceqcs=type*(cqcs(lx0+loct)-vcs*gccs)
        !          2654:       call intgr8(geqbx,ceq,capbx,loct+14)
        !          2655:       ceqbx=type*(cqbx(lx0+loct)-vbx*geqbx)
        !          2656:       if (initf.ne.5) go to 700
        !          2657:       cqbe(lx1+loct)=cqbe(lx0+loct)
        !          2658:       cqbc(lx1+loct)=cqbc(lx0+loct)
        !          2659:       cqbx(lx1+loct)=cqbx(lx0+loct)
        !          2660:       cqcs(lx1+loct)=cqcs(lx0+loct)
        !          2661: c
        !          2662: c  check convergence
        !          2663: c
        !          2664:   700 if (initf.ne.3) go to 710
        !          2665:       if (ioff.eq.0) go to 710
        !          2666:       go to 750
        !          2667:   710 if (icheck.eq.1) go to 720
        !          2668:       tol=reltol*dmax1(dabs(cchat),dabs(cc))+abstol
        !          2669:       if (dabs(cchat-cc).gt.tol) go to 720
        !          2670:       tol=reltol*dmax1(dabs(cbhat),dabs(cb))+abstol
        !          2671:       if (dabs(cbhat-cb).le.tol) go to 750
        !          2672:   720 noncon=noncon+1
        !          2673:   750 vbeo(lx0+loct)=vbe
        !          2674:       vbco(lx0+loct)=vbc
        !          2675:       cco(lx0+loct)=cc
        !          2676:       cbo(lx0+loct)=cb
        !          2677:       gpio(lx0+loct)=gpi
        !          2678:       gmuo(lx0+loct)=gmu
        !          2679:       gmo(lx0+loct)=gm
        !          2680:       goo(lx0+loct)=go
        !          2681:       gxo(lx0+loct)=gx
        !          2682: c
        !          2683: c  load current excitation vector
        !          2684: c
        !          2685:   900 ceqbe=type*(cc+cb-vbe*(gm+go+gpi)+vbc*(go-geqcb))
        !          2686:       ceqbc=type*(-cc+vbe*(gm+go)-vbc*(gmu+go))
        !          2687:       value(lvn+node2)=value(lvn+node2)-ceqbx
        !          2688:       value(lvn+node4)=value(lvn+node4)+ceqcs+ceqbx+ceqbc
        !          2689:       value(lvn+node5)=value(lvn+node5)-ceqbe-ceqbc
        !          2690:       value(lvn+node6)=value(lvn+node6)+ceqbe
        !          2691:       value(lvn+node7)=value(lvn+node7)-ceqcs
        !          2692: c
        !          2693: c  load y matrix
        !          2694: c
        !          2695:       locy=lynl+nodplc(loc+24)
        !          2696:       value(locy)=value(locy)+gcpr
        !          2697:       locy=lynl+nodplc(loc+25)
        !          2698:       value(locy)=value(locy)+gx+geqbx
        !          2699:       locy=lynl+nodplc(loc+26)
        !          2700:       value(locy)=value(locy)+gepr
        !          2701:       locy=lynl+nodplc(loc+27)
        !          2702:       value(locy)=value(locy)+gmu+go+gcpr+gccs+geqbx
        !          2703:       locy=lynl+nodplc(loc+28)
        !          2704:       value(locy)=value(locy)+gx  +gpi+gmu+geqcb
        !          2705:       locy=lynl+nodplc(loc+29)
        !          2706:       value(locy)=value(locy)+gpi+gepr+gm+go
        !          2707:       locy=lynl+nodplc(loc+10)
        !          2708:       value(locy)=value(locy)-gcpr
        !          2709:       locy=lynl+nodplc(loc+11)
        !          2710:       value(locy)=value(locy)-gx
        !          2711:       locy=lynl+nodplc(loc+12)
        !          2712:       value(locy)=value(locy)-gepr
        !          2713:       locy=lynl+nodplc(loc+13)
        !          2714:       value(locy)=value(locy)-gcpr
        !          2715:       locy=lynl+nodplc(loc+14)
        !          2716:       value(locy)=value(locy)-gmu+gm
        !          2717:       locy=lynl+nodplc(loc+15)
        !          2718:       value(locy)=value(locy)-gm-go
        !          2719:       locy=lynl+nodplc(loc+16)
        !          2720:       value(locy)=value(locy)-gx
        !          2721:       locy=lynl+nodplc(loc+17)
        !          2722:       value(locy)=value(locy)-gmu-geqcb
        !          2723:       locy=lynl+nodplc(loc+18)
        !          2724:       value(locy)=value(locy)-gpi
        !          2725:       locy=lynl+nodplc(loc+19)
        !          2726:       value(locy)=value(locy)-gepr
        !          2727:       locy=lynl+nodplc(loc+20)
        !          2728:       value(locy)=value(locy)-go+geqcb
        !          2729:       locy=lynl+nodplc(loc+21)
        !          2730:       value(locy)=value(locy)-gpi-gm-geqcb
        !          2731:       locy=lynl+nodplc(loc+31)
        !          2732:       value(locy)=value(locy)+gccs
        !          2733:       locy=lynl+nodplc(loc+32)
        !          2734:       value(locy)=value(locy)-gccs
        !          2735:       locy=lynl+nodplc(loc+33)
        !          2736:       value(locy)=value(locy)-gccs
        !          2737:       locy=lynl+nodplc(loc+34)
        !          2738:       value(locy)=value(locy)-geqbx
        !          2739:       locy=lynl+nodplc(loc+35)
        !          2740:       value(locy)=value(locy)-geqbx
        !          2741:  1000 loc=nodplc(loc)
        !          2742:       go to 10
        !          2743:       end
        !          2744:       subroutine fetlim(vnew,vold,vto,icheck)
        !          2745: c     
        !          2746: c     *** fetlim is not used in this version ***
        !          2747: c     *   if problems arrise with the conver-  *
        !          2748: c     *   gence of MOSFET circuit it should be *
        !          2749: c     *   re-installed.   R.Newton.            *
        !          2750: c     ***                                    ***
        !          2751: c
        !          2752:       return
        !          2753:       end
        !          2754:       subroutine jfet
        !          2755:       implicit double precision (a-h,o-z)
        !          2756: c
        !          2757: c     this routine processes jfets for dc and transient analyses.
        !          2758: c
        !          2759:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !          2760:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !          2761:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !          2762:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !          2763:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !          2764:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !          2765:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !          2766:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !          2767:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          2768:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          2769:      2   itemno,nosolv,ipostp,iscrch
        !          2770:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
        !          2771:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
        !          2772:       common /blank/ value(1000)
        !          2773:       integer nodplc(64)
        !          2774:       complex*16 cvalue(32)
        !          2775:       equivalence (value(1),nodplc(1),cvalue(1))
        !          2776: c
        !          2777: c
        !          2778:       dimension vgso(1),vgdo(1),cgo(1),cdo(1),cgdo(1),gmo(1),gdso(1),
        !          2779:      1   ggso(1),ggdo(1),qgs(1),cqgs(1),qgd(1),cqgd(1)
        !          2780:       equivalence (vgso(1),value( 1)),(vgdo(1),value( 2)),
        !          2781:      1            (cgo (1),value( 3)),(cdo (1),value( 4)),
        !          2782:      2            (cgdo(1),value( 5)),(gmo (1),value( 6)),
        !          2783:      3            (gdso(1),value( 7)),(ggso(1),value( 8)),
        !          2784:      4            (ggdo(1),value( 9)),(qgs (1),value(10)),
        !          2785:      5            (cqgs(1),value(11)),(qgd (1),value(12)),
        !          2786:      6            (cqgd(1),value(13))
        !          2787: c
        !          2788: c
        !          2789:       loc=locate(13)
        !          2790:    10 if (loc.eq.0) return
        !          2791:       locv=nodplc(loc+1)
        !          2792:       node1=nodplc(loc+2)
        !          2793:       node2=nodplc(loc+3)
        !          2794:       node3=nodplc(loc+4)
        !          2795:       node4=nodplc(loc+5)
        !          2796:       node5=nodplc(loc+6)
        !          2797:       locm=nodplc(loc+7)
        !          2798:       ioff=nodplc(loc+8)
        !          2799:       type=nodplc(locm+2)
        !          2800:       locm=nodplc(locm+1)
        !          2801:       loct=nodplc(loc+19)
        !          2802: c
        !          2803: c  dc model parameters
        !          2804: c
        !          2805:       area=value(locv+1)
        !          2806:       vto=value(locm+1)
        !          2807:       beta=value(locm+2)*area
        !          2808:       xlamb=value(locm+3)
        !          2809:       gdpr=value(locm+4)*area
        !          2810:       gspr=value(locm+5)*area
        !          2811:       csat=value(locm+9)*area
        !          2812:       vcrit=value(locm+16)
        !          2813: c
        !          2814: c  initialization
        !          2815: c
        !          2816:       icheck=1
        !          2817:       go to (100,20,30,50,60,70), initf
        !          2818:    20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25
        !          2819:       vds=type*value(locv+2)
        !          2820:       vgs=type*value(locv+3)
        !          2821:       vgd=vgs-vds
        !          2822:       go to 300
        !          2823:    25 if(ioff.ne.0) go to 40
        !          2824:       vgs=-1.0d0
        !          2825:       vgd=-1.0d0
        !          2826:       go to 300
        !          2827:    30 if (ioff.eq.0) go to 100
        !          2828:    40 vgs=0.0d0
        !          2829:       vgd=0.0d0
        !          2830:       go to 300
        !          2831:    50 vgs=vgso(lx0+loct)
        !          2832:       vgd=vgdo(lx0+loct)
        !          2833:       go to 300
        !          2834:    60 vgs=vgso(lx1+loct)
        !          2835:       vgd=vgdo(lx1+loct)
        !          2836:       go to 300
        !          2837:    70 xfact=delta/delold(2)
        !          2838:       vgso(lx0+loct)=vgso(lx1+loct)
        !          2839:       vgs=(1.0d0+xfact)*vgso(lx1+loct)-xfact*vgso(lx2+loct)
        !          2840:       vgdo(lx0+loct)=vgdo(lx1+loct)
        !          2841:       vgd=(1.0d0+xfact)*vgdo(lx1+loct)-xfact*vgdo(lx2+loct)
        !          2842:       cgo(lx0+loct)=cgo(lx1+loct)
        !          2843:       cdo(lx0+loct)=cdo(lx1+loct)
        !          2844:       cgdo(lx0+loct)=cgdo(lx1+loct)
        !          2845:       gmo(lx0+loct)=gmo(lx1+loct)
        !          2846:       gdso(lx0+loct)=gdso(lx1+loct)
        !          2847:       ggso(lx0+loct)=ggso(lx1+loct)
        !          2848:       ggdo(lx0+loct)=ggdo(lx1+loct)
        !          2849:       go to 110
        !          2850: c
        !          2851: c  compute new nonlinear branch voltages
        !          2852: c
        !          2853:   100 vgs=type*(value(lvnim1+node2)-value(lvnim1+node5))
        !          2854:       vgd=type*(value(lvnim1+node2)-value(lvnim1+node4))
        !          2855:   110 delvgs=vgs-vgso(lx0+loct)
        !          2856:       delvgd=vgd-vgdo(lx0+loct)
        !          2857:       delvds=delvgs-delvgd
        !          2858:       cghat=cgo(lx0+loct)+ggdo(lx0+loct)*delvgd+ggso(lx0+loct)*delvgs
        !          2859:       cdhat=cdo(lx0+loct)+gmo(lx0+loct)*delvgs+gdso(lx0+loct)*delvds
        !          2860:      1   -ggdo(lx0+loct)*delvgd
        !          2861: c
        !          2862: c  bypass if solution has not changed
        !          2863: c
        !          2864:       if (initf.eq.6) go to 200
        !          2865:       tol=reltol*dmax1(dabs(vgs),dabs(vgso(lx0+loct)))+vntol
        !          2866:       if (dabs(delvgs).ge.tol) go to 200
        !          2867:       tol=reltol*dmax1(dabs(vgd),dabs(vgdo(lx0+loct)))+vntol
        !          2868:       if (dabs(delvgd).ge.tol) go to 200
        !          2869:       tol=reltol*dmax1(dabs(cghat),dabs(cgo(lx0+loct)))+abstol
        !          2870:       if (dabs(cghat-cgo(lx0+loct)).ge.tol) go to 200
        !          2871:       tol=reltol*dmax1(dabs(cdhat),dabs(cdo(lx0+loct)))+abstol
        !          2872:       if (dabs(cdhat-cdo(lx0+loct)).ge.tol) go to 200
        !          2873:       vgs=vgso(lx0+loct)
        !          2874:       vgd=vgdo(lx0+loct)
        !          2875:       vds=vgs-vgd
        !          2876:       cg=cgo(lx0+loct)
        !          2877:       cd=cdo(lx0+loct)
        !          2878:       cgd=cgdo(lx0+loct)
        !          2879:       gm=gmo(lx0+loct)
        !          2880:       gds=gdso(lx0+loct)
        !          2881:       ggs=ggso(lx0+loct)
        !          2882:       ggd=ggdo(lx0+loct)
        !          2883:       go to 900
        !          2884: c
        !          2885: c  limit nonlinear branch voltages
        !          2886: c
        !          2887:   200 icheck=0
        !          2888:       call pnjlim(vgs,vgso(lx0+loct),vt,vcrit,icheck)
        !          2889:       call pnjlim(vgd,vgdo(lx0+loct),vt,vcrit,icheck)
        !          2890:       call fetlim(vgs,vgso(lx0+loct),vto,icheck)
        !          2891:       call fetlim(vgd,vgdo(lx0+loct),vto,icheck)
        !          2892: c
        !          2893: c  determine dc current and derivatives
        !          2894: c
        !          2895:   300 vds=vgs-vgd
        !          2896:       if (vgs.gt.-5.0d0*vt) go to 310
        !          2897:       ggs=-csat/vgs+gmin
        !          2898:       cg=ggs*vgs
        !          2899:       go to 320
        !          2900:   310 evgs=dexp(vgs/vt)
        !          2901:       ggs=csat*evgs/vt+gmin
        !          2902:       cg=csat*(evgs-1.0d0)+gmin*vgs
        !          2903:   320 if (vgd.gt.-5.0d0*vt) go to 330
        !          2904:       ggd=-csat/vgd+gmin
        !          2905:       cgd=ggd*vgd
        !          2906:       go to 340
        !          2907:   330 evgd=dexp(vgd/vt)
        !          2908:       ggd=csat*evgd/vt+gmin
        !          2909:       cgd=csat*(evgd-1.0d0)+gmin*vgd
        !          2910:   340 cg=cg+cgd
        !          2911: c
        !          2912: c  compute drain current and derivitives for normal mode
        !          2913: c
        !          2914:   400 if (vds.lt.0.0d0) go to 450
        !          2915:       vgst=vgs-vto
        !          2916: c
        !          2917: c  normal mode, cutoff region
        !          2918: c
        !          2919:       if (vgst.gt.0.0d0) go to 410
        !          2920:       cdrain=0.0d0
        !          2921:       gm=0.0d0
        !          2922:       gds=0.0d0
        !          2923:       go to 490
        !          2924: c
        !          2925: c  normal mode, saturation region
        !          2926: c
        !          2927:   410 betap=beta*(1.0d0+xlamb*vds)
        !          2928:       twob=betap+betap
        !          2929:       if (vgst.gt.vds) go to 420
        !          2930:       cdrain=betap*vgst*vgst
        !          2931:       gm=twob*vgst
        !          2932:       gds=xlamb*beta*vgst*vgst
        !          2933:       go to 490
        !          2934: c
        !          2935: c  normal mode, linear region
        !          2936: c
        !          2937:   420 cdrain=betap*vds*(vgst+vgst-vds)
        !          2938:       gm=twob*vds
        !          2939:       gds=twob*(vgst-vds)+xlamb*beta*vds*(vgst+vgst-vds)
        !          2940:       go to 490
        !          2941: c
        !          2942: c  compute drain current and derivitives for inverse mode
        !          2943: c
        !          2944:   450 vgdt=vgd-vto
        !          2945: c
        !          2946: c  inverse mode, cutoff region
        !          2947: c
        !          2948:       if (vgdt.gt.0.0d0) go to 460
        !          2949:       cdrain=0.0d0
        !          2950:       gm=0.0d0
        !          2951:       gds=0.0d0
        !          2952:       go to 490
        !          2953: c
        !          2954: c  inverse mode, saturation region
        !          2955: c
        !          2956:   460 betap=beta*(1.0d0-xlamb*vds)
        !          2957:       twob=betap+betap
        !          2958:       if (vgdt.gt.-vds) go to 470
        !          2959:       cdrain=-betap*vgdt*vgdt
        !          2960:       gm=-twob*vgdt
        !          2961:       gds=xlamb*beta*vgdt*vgdt-gm
        !          2962:       go to 490
        !          2963: c
        !          2964: c  inverse mode, linear region
        !          2965: c
        !          2966:   470 cdrain=betap*vds*(vgdt+vgdt+vds)
        !          2967:       gm=twob*vds
        !          2968:       gds=twob*vgdt-xlamb*beta*vds*(vgdt+vgdt+vds)
        !          2969: c
        !          2970: c  compute equivalent drain current source
        !          2971: c
        !          2972:   490 cd=cdrain-cgd
        !          2973:       if (mode.ne.1) go to 500
        !          2974:       if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
        !          2975:       if (initf.eq.4) go to 500
        !          2976:       go to 700
        !          2977: c
        !          2978: c  charge storage elements
        !          2979: c
        !          2980:   500 czgs=value(locm+6)*area
        !          2981:       czgd=value(locm+7)*area
        !          2982:       phib=value(locm+8)
        !          2983:       twop=phib+phib
        !          2984:       fcpb=value(locm+12)
        !          2985:       fcpb2=fcpb*fcpb
        !          2986:       f1=value(locm+13)
        !          2987:       f2=value(locm+14)
        !          2988:       f3=value(locm+15)
        !          2989:       czgsf2=czgs/f2
        !          2990:       czgdf2=czgd/f2
        !          2991:       if (vgs.ge.fcpb) go to 510
        !          2992:       sarg=dsqrt(1.0d0-vgs/phib)
        !          2993:       qgs(lx0+loct)=twop*czgs*(1.0d0-sarg)
        !          2994:       capgs=czgs/sarg
        !          2995:       go to 520
        !          2996:   510 qgs(lx0+loct)=czgs*f1+czgsf2*(f3*(vgs-fcpb)
        !          2997:      1   +(vgs*vgs-fcpb2)/(twop+twop))
        !          2998:       capgs=czgsf2*(f3+vgs/twop)
        !          2999:   520 if (vgd.ge.fcpb) go to 530
        !          3000:       sarg=dsqrt(1.0d0-vgd/phib)
        !          3001:       qgd(lx0+loct)=twop*czgd*(1.0d0-sarg)
        !          3002:       capgd=czgd/sarg
        !          3003:       go to 560
        !          3004:   530 qgd(lx0+loct)=czgd*f1+czgdf2*(f3*(vgd-fcpb)
        !          3005:      1   +(vgd*vgd-fcpb2)/(twop+twop))
        !          3006:       capgd=czgdf2*(f3+vgd/twop)
        !          3007: c
        !          3008: c  store small-signal parameters
        !          3009: c
        !          3010:   560 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700
        !          3011:       if (initf.ne.4) go to 600
        !          3012:       value(lx0+loct+9)=capgs
        !          3013:       value(lx0+loct+11)=capgd
        !          3014:       go to 1000
        !          3015: c
        !          3016: c  transient analysis
        !          3017: c
        !          3018:   600 if (initf.ne.5) go to 610
        !          3019:       qgs(lx1+loct)=qgs(lx0+loct)
        !          3020:       qgd(lx1+loct)=qgd(lx0+loct)
        !          3021:   610 call intgr8(geq,ceq,capgs,loct+9)
        !          3022:       ggs=ggs+geq
        !          3023:       cg=cg+cqgs(lx0+loct)
        !          3024:       call intgr8(geq,ceq,capgd,loct+11)
        !          3025:       ggd=ggd+geq
        !          3026:       cg=cg+cqgd(lx0+loct)
        !          3027:       cd=cd-cqgd(lx0+loct)
        !          3028:       cgd=cgd+cqgd(lx0+loct)
        !          3029:       if (initf.ne.5) go to 700
        !          3030:       cqgs(lx1+loct)=cqgs(lx0+loct)
        !          3031:       cqgd(lx1+loct)=cqgd(lx0+loct)
        !          3032: c
        !          3033: c  check convergence
        !          3034: c
        !          3035:   700 if (initf.ne.3) go to 710
        !          3036:       if (ioff.eq.0) go to 710
        !          3037:       go to 750
        !          3038:   710 if (icheck.eq.1) go to 720
        !          3039:       tol=reltol*dmax1(dabs(cghat),dabs(cg))+abstol
        !          3040:       if (dabs(cghat-cg).ge.tol) go to 720
        !          3041:       tol=reltol*dmax1(dabs(cdhat),dabs(cd))+abstol
        !          3042:       if (dabs(cdhat-cd).le.tol) go to 750
        !          3043:   720 noncon=noncon+1
        !          3044:   750 vgso(lx0+loct)=vgs
        !          3045:       vgdo(lx0+loct)=vgd
        !          3046:       cgo(lx0+loct)=cg
        !          3047:       cdo(lx0+loct)=cd
        !          3048:       cgdo(lx0+loct)=cgd
        !          3049:       gmo(lx0+loct)=gm
        !          3050:       gdso(lx0+loct)=gds
        !          3051:       ggso(lx0+loct)=ggs
        !          3052:       ggdo(lx0+loct)=ggd
        !          3053: c
        !          3054: c  load current vector
        !          3055: c
        !          3056:   900 ceqgd=type*(cgd-ggd*vgd)
        !          3057:       ceqgs=type*((cg-cgd)-ggs*vgs)
        !          3058:       cdreq=type*((cd+cgd)-gds*vds-gm*vgs)
        !          3059:       value(lvn+node2)=value(lvn+node2)-ceqgs-ceqgd
        !          3060:       value(lvn+node4)=value(lvn+node4)-cdreq+ceqgd
        !          3061:       value(lvn+node5)=value(lvn+node5)+cdreq+ceqgs
        !          3062: c
        !          3063: c  load y matrix
        !          3064: c
        !          3065:       locy=lynl+nodplc(loc+20)
        !          3066:       value(locy)=value(locy)+gdpr
        !          3067:       locy=lynl+nodplc(loc+21)
        !          3068:       value(locy)=value(locy)+ggd+ggs
        !          3069:       locy=lynl+nodplc(loc+22)
        !          3070:       value(locy)=value(locy)+gspr
        !          3071:       locy=lynl+nodplc(loc+23)
        !          3072:       value(locy)=value(locy)+gdpr+gds+ggd
        !          3073:       locy=lynl+nodplc(loc+24)
        !          3074:       value(locy)=value(locy)+gspr+gds+gm+ggs
        !          3075:       locy=lynl+nodplc(loc+9)
        !          3076:       value(locy)=value(locy)-gdpr
        !          3077:       locy=lynl+nodplc(loc+10)
        !          3078:       value(locy)=value(locy)-ggd
        !          3079:       locy=lynl+nodplc(loc+11)
        !          3080:       value(locy)=value(locy)-ggs
        !          3081:       locy=lynl+nodplc(loc+12)
        !          3082:       value(locy)=value(locy)-gspr
        !          3083:       locy=lynl+nodplc(loc+13)
        !          3084:       value(locy)=value(locy)-gdpr
        !          3085:       locy=lynl+nodplc(loc+14)
        !          3086:       value(locy)=value(locy)+gm-ggd
        !          3087:       locy=lynl+nodplc(loc+15)
        !          3088:       value(locy)=value(locy)-gds-gm
        !          3089:       locy=lynl+nodplc(loc+16)
        !          3090:       value(locy)=value(locy)-ggs-gm
        !          3091:       locy=lynl+nodplc(loc+17)
        !          3092:       value(locy)=value(locy)-gspr
        !          3093:       locy=lynl+nodplc(loc+18)
        !          3094:       value(locy)=value(locy)-gds
        !          3095:  1000 loc=nodplc(loc)
        !          3096:       go to 10
        !          3097:       end
        !          3098:       subroutine mosfet
        !          3099:       implicit double precision (a-h,o-z)
        !          3100: c
        !          3101: c     this routine processes mosfets for dc and transient analyses.
        !          3102: c
        !          3103:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
        !          3104:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
        !          3105:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
        !          3106:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
        !          3107:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
        !          3108:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
        !          3109:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
        !          3110:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
        !          3111:       common /mosarg/ gamma,beta,vto,phi,cox,vbi,xnfs,xnsub,xd,xj,xl,
        !          3112:      1   xlamda,utra,uexp,vbp,von,vdsat,theta,vcrit,vtra,gleff,cdrain
        !          3113:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          3114:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          3115:      2   itemno,nosolv,ipostp,iscrch
        !          3116:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
        !          3117:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
        !          3118:       common /blank/ value(1000)
        !          3119:       integer nodplc(64)
        !          3120:       complex*16 cvalue(32)
        !          3121:       equivalence (value(1),nodplc(1),cvalue(1))
        !          3122: c
        !          3123: c
        !          3124:       dimension vbdo(1),vbso(1),vgbo(1),gdgo(1),cdo(1),cbo(1),
        !          3125:      1   gddo(1),gdso(1),gbgo(1),gbdo(1),gbso(1),qb(1),cqb(1),qg(1),
        !          3126:      2   cqg(1),qd(1),cqd(1)
        !          3127: c.. note: for direct indexing with 'value', use, e.g. loct+2 to find vgbo
        !          3128:       equivalence (vbdo (1),value( 1)),(vbso (1),value( 2)),
        !          3129:      1            (vgbo (1),value( 3)),
        !          3130:      2            (cdo  (1),value( 5)),(cbo  (1),value( 6)),
        !          3131:      3            (gddo (1),value( 7)),(gdgo (1),value(8)),
        !          3132:      4            (gdso (1),value( 9)),(gbgo (1),value(10)),
        !          3133:      5            (gbdo (1),value(11)),(gbso (1),value(12)),
        !          3134:      6            (qb   (1),value(13)),(cqb  (1),value(14)),
        !          3135:      7            (qg   (1),value(15)),(cqg  (1),value(16)),
        !          3136:      8            (qd   (1),value(17)),(cqd  (1),value(18))
        !          3137: c
        !          3138: c
        !          3139:       loc=locate(14)
        !          3140:    10 if (loc.eq.0) return
        !          3141:       locm=nodplc(loc+8)
        !          3142:       if(nodplc(locm+2).ne.0) go to 15
        !          3143:       call gasfet(loc)
        !          3144:       go to 1000
        !          3145:    15 locv=nodplc(loc+1)
        !          3146:       node1=nodplc(loc+2)
        !          3147:       node2=nodplc(loc+3)
        !          3148:       node3=nodplc(loc+4)
        !          3149:       node4=nodplc(loc+5)
        !          3150:       node5=nodplc(loc+6)
        !          3151:       node6=nodplc(loc+7)
        !          3152:       ioff=nodplc(loc+9)
        !          3153:       type=nodplc(locm+2)
        !          3154:       locm=nodplc(locm+1)
        !          3155:       loct=nodplc(loc+26)
        !          3156: c
        !          3157: c  dc model parameters
        !          3158: c
        !          3159:       xj=value(locm+19)
        !          3160:       xl=value(locv+1)-2.0d0*value(locm+20)
        !          3161:       xw=value(locv+2)-2.0d0*value(locm+36)
        !          3162:       devmod=value(locv+8)
        !          3163:       vto=type*value(locm+1)
        !          3164:       beta=value(locm+2)*xw/xl
        !          3165:       gamma=value(locm+3)
        !          3166:       phi=value(locm+4)
        !          3167:       xlamda=value(locm+5)
        !          3168:       csat=value(locm+15)
        !          3169:       ad=value(locv+3)
        !          3170:       as=value(locv+4)
        !          3171:       cdsat=csat*ad
        !          3172:       cssat=csat*as
        !          3173:       gdpr=value(locv+11)
        !          3174:       gspr=value(locv+12)
        !          3175:       covlgs=value(locm+8)*xw
        !          3176:       covlgd=value(locm+9)*xw
        !          3177:       covlgb=value(locm+10)*xl
        !          3178:       cox=value(locm+13)
        !          3179:       xnsub=value(locm+16)
        !          3180:       xnfs=value(locm+18)
        !          3181:       vbp=value(locm+24)
        !          3182:       uexp=value(locm+25)
        !          3183:       utra=value(locm+26)
        !          3184:       vbi=type*value(locm+34)
        !          3185:       xd=value(locm+35)
        !          3186:       vcrit=value(locm+37)*xl
        !          3187:       vtra=value(locm+38)*xl
        !          3188:       theta=value(locm+39)
        !          3189:       gleff=value(locm+40)
        !          3190: c
        !          3191: c  initialization
        !          3192: c
        !          3193:       icheck=1
        !          3194:       go to (100,20,30,50,60,70), initf
        !          3195:    20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25
        !          3196:       vbs=type*value(locv+7)
        !          3197:       vbd=type*value(locv+5)-vbs
        !          3198:       vgb=type*value(locv+6)-vbs
        !          3199:       go to 300
        !          3200:    25 if(ioff.ne.0) go to 40
        !          3201:       vbs=0.0d0
        !          3202:       vbd=0.0d0
        !          3203:       vgb=vto
        !          3204:       go to 300
        !          3205:    30 if (ioff.eq.0) go to 100
        !          3206:    40 vbs=0.0d0
        !          3207:       vbd=0.0d0
        !          3208:       vgb=0.0d0
        !          3209:       go to 300
        !          3210:    50 vbs=vbso(lx0+loct)
        !          3211:       vbd=vbdo(lx0+loct)
        !          3212:       vgb=vgbo(lx0+loct)
        !          3213:       go to 300
        !          3214:    60 vbs=vbso(lx1+loct)
        !          3215:       vbd=vbdo(lx1+loct)
        !          3216:       vgb=vgbo(lx1+loct)
        !          3217:       go to 300
        !          3218:    70 xfact=delta/delold(2)
        !          3219:       vbso(lx0+loct)=vbso(lx1+loct)
        !          3220:       vbs=(1.0d0+xfact)*vbso(lx1+loct)-xfact*vbso(lx2+loct)
        !          3221:       vbdo(lx0+loct)=vbdo(lx1+loct)
        !          3222:       vbd=(1.0d0+xfact)*vbdo(lx1+loct)-xfact*vbdo(lx2+loct)
        !          3223:       vgbo(lx0+loct)=vgbo(lx1+loct)
        !          3224:       vgb=(1.0d0+xfact)*vgbo(lx1+loct)-xfact*vgbo(lx2+loct)
        !          3225:       cdo(lx0+loct)=cdo(lx1+loct)
        !          3226:       cbo(lx0+loct)=cbo(lx1+loct)
        !          3227:       gdgo(lx0+loct)=gdgo(lx1+loct)
        !          3228:       gddo(lx0+loct)=gddo(lx1+loct)
        !          3229:       gdso(lx0+loct)=gdso(lx1+loct)
        !          3230:       gbgo(lx0+loct)=gbgo(lx1+loct)
        !          3231:       gbdo(lx0+loct)=gbdo(lx1+loct)
        !          3232:       gbso(lx0+loct)=gbso(lx1+loct)
        !          3233:       go to 110
        !          3234: c
        !          3235: c  compute new nonlinear branch voltages
        !          3236: c
        !          3237:   100 vbs=type*(value(lvnim1+node4)-value(lvnim1+node6))
        !          3238:       vbd=type*(value(lvnim1+node4)-value(lvnim1+node5))
        !          3239:       vgb=type*(value(lvnim1+node2)-value(lvnim1+node4))
        !          3240:   110 delvbs=vbs-vbso(lx0+loct)
        !          3241:       delvbd=vbd-vbdo(lx0+loct)
        !          3242:       delvgb=vgb-vgbo(lx0+loct)
        !          3243:       cdhat=cdo(lx0+loct)+gdgo(lx0+loct)*delvgb-gddo(lx0+loct)*delvbd
        !          3244:      1  -gdso(lx0+loct)*delvbs
        !          3245:       cbhat=cbo(lx0+loct)+gbgo(lx0+loct)*delvgb-gbdo(lx0+loct)*delvbd
        !          3246:      1   -gbso(lx0+loct)*delvbs
        !          3247: c
        !          3248: c  bypass if solution has not changed
        !          3249: c
        !          3250: c********** kill bypass for now!!!!!
        !          3251:       if (6    .eq.6) go to 200
        !          3252:       tol=reltol*dmax1(dabs(vbs),dabs(vbso(lx0+loct)))+vntol
        !          3253:       if (dabs(delvbs).ge.tol) go to 200
        !          3254:       tol=reltol*dmax1(dabs(vbd),dabs(vbdo(lx0+loct)))+vntol
        !          3255:       if (dabs(delvbd).ge.tol) go to 200
        !          3256:       tol=reltol*dmax1(dabs(vgb),dabs(vgbo(lx0+loct)))+vntol
        !          3257:       if (dabs(delvgb).ge.tol) go to 200
        !          3258:       tol=reltol*dmax1(dabs(cdhat),dabs(cdo(lx0+loct)))+abstol
        !          3259:       if (dabs(cdhat-cdo(lx0+loct)).ge.tol) go to 200
        !          3260:       tol=reltol*dmax1(dabs(cbhat),dabs(cbo(lx0+loct)))+abstol
        !          3261:       if (dabs(cbhat-(cbo(lx0+loct))).ge.tol) go to 200
        !          3262:       vbd=vbdo(lx0+loct)
        !          3263:       vbs=vbso(lx0+loct)
        !          3264:       vgb=vgbo(lx0+loct)
        !          3265:       cdrain=cdo(lx0+loct)
        !          3266:       cbulk=cbo(lx0+loct)
        !          3267:       gccdg=gdgo(lx0+loct)
        !          3268:       gccdd=gddo(lx0+loct)
        !          3269:       gccds=gdso(lx0+loct)
        !          3270:       gccbg=gbgo(lx0+loct)
        !          3271:       gccbd=gbdo(lx0+loct)
        !          3272:       gccbs=gbso(lx0+loct)
        !          3273:       go to 900
        !          3274: c
        !          3275: c  limit nonlinear branch voltages
        !          3276: c
        !          3277:   200 von=type*value(locv+9)
        !          3278:       icheck=0
        !          3279:       call fetlim(vgb,vgbo(lx0+loct),von,icheck)
        !          3280:       vcornr=0.0d0
        !          3281: c     if(vbs.gt.0.0d0) vcornr=vt*dlog(vt/(root2*cssat))
        !          3282:       call pnjlim(vbs,vbso(lx0+loct),vt,vcornr,icheck)
        !          3283: c     vbs=dmax1(vbso(lx0+loct)-10.0d0,vbs)
        !          3284:       vcornr=0.0d0
        !          3285: c     if(vbd.gt.0.0d0) vcornr=vt*dlog(vt/(root2*cdsat))
        !          3286:       call pnjlim(vbd,vbdo(lx0+loct),vt,vcornr,icheck)
        !          3287: c     vbd=dmax1(vbdo(lx0+loct)-10.0d0,vbd)
        !          3288: c
        !          3289: c  determine bulk-drain and bulk-source diode terms
        !          3290: c
        !          3291:   300 fivevt=-5.0d0*vt
        !          3292:       if (vbs.gt.fivevt) go to 310
        !          3293:       geqbs=-cssat/vbs+gmin
        !          3294:       cbulk=geqbs*vbs
        !          3295:       go to 320
        !          3296:   310 evbs=dexp(vbs/vt)
        !          3297:       geqbs=cssat*evbs/vt+gmin
        !          3298:       cbulk=cssat*(evbs-1.0d0)+gmin*vbs
        !          3299:   320 if (vbd.gt.fivevt) go to 330
        !          3300:       geqbd=-cdsat/vbd+gmin
        !          3301:       cbd=geqbd*vbd
        !          3302:       cbulk=cbulk+cbd
        !          3303:       go to 340
        !          3304:   330 evbd=dexp(vbd/vt)
        !          3305:       geqbd=cdsat*evbd/vt+gmin
        !          3306:       cbd=cdsat*(evbd-1.0d0)+gmin*vbd
        !          3307:       cbulk=cbulk+cbd
        !          3308: c.. cbd must also be subtracted from drain current
        !          3309:   340 continue
        !          3310:       gccdd=geqbd
        !          3311:       gccbd=-geqbd
        !          3312:       gccss=geqbs
        !          3313:       gccbs=-geqbs
        !          3314:       if (mode.ne.1) go to 350
        !          3315: c.. zero out some conductances and cgate
        !          3316:       cgate=0.0d0
        !          3317:       gccgg=0.0d0
        !          3318:       gccgd=0.0d0
        !          3319:       gccgs=0.0d0
        !          3320:       gccbg=0.0d0
        !          3321: c
        !          3322: c  compute drain current and derivatives
        !          3323: c
        !          3324:   350 cox=cox*xl*xw
        !          3325:       if(vbd.gt.vbs) go to 360
        !          3326: c.. normal operation
        !          3327:       devmod=1.0d0
        !          3328:       call calcq(vgb,vbd,vbs,qgate,qchan,qbulk,
        !          3329:      1 ccgg,ccgd,ccgs,ccbg,ccbd,ccbs,didvg,didvd,didvs)
        !          3330:       didvg=beta*didvg
        !          3331:       didvd=beta*didvd
        !          3332:       didvs=beta*didvs
        !          3333:       go to 370
        !          3334: c.. inverted operation
        !          3335:   360 devmod=-1.0d0
        !          3336:       call calcq(vgb,vbs,vbd,qgate,qchan,qbulk,
        !          3337:      1 ccgg,ccgs,ccgd,ccbg,ccbs,ccbd,didvg,didvs,didvd)
        !          3338:       didvg=-beta*didvg
        !          3339:       didvd=-beta*didvd
        !          3340:       didvs=-beta*didvs
        !          3341:       cdrain=-cdrain
        !          3342:   370 cdrain=beta*cdrain-cbd
        !          3343: c     if(mode.ne.1) write(6,6666) qgate,qchan,qbulk,time,delta
        !          3344:  6666 format(' qg qc qb',1p3e11.3,' time delta ',2e11.3)
        !          3345:       value(locv+8)=devmod
        !          3346:       value(locv+9)=type*von
        !          3347:       value(locv+10)=type*vdsat
        !          3348:       if(mode.ne.1) go to 500
        !          3349:       if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
        !          3350:       if (initf.eq.4) go to 500
        !          3351:       gccdg=didvg
        !          3352:       gccdd=gccdd+didvd
        !          3353:       gccds=didvs
        !          3354:       gccsg=-didvg
        !          3355:       gccsd=-didvd
        !          3356:       gccss=gccss-didvs
        !          3357:       go to 700
        !          3358: c
        !          3359: c  charge storage elements
        !          3360: c
        !          3361: c.. bulk-drain and bulk-source depletion capacitances
        !          3362:   500 czbd=value(locm+11)*ad
        !          3363:       czbs=value(locm+12)*as
        !          3364:       phib=value(locm+14)
        !          3365:       twop=phib+phib
        !          3366:       fcpb=value(locm+29)
        !          3367:       if(vbs.lt.fcpb.and.vbd.lt.fcpb) go to 505
        !          3368:       fcpb2=fcpb*fcpb
        !          3369:       f1=value(locm+30)
        !          3370:       f2=value(locm+31)
        !          3371:       f3=value(locm+32)
        !          3372:       czbsf2=czbs/f2
        !          3373:       czbdf2=czbd/f2
        !          3374:       if (vbs.ge.fcpb) go to 510
        !          3375:   505 sarg=dsqrt(1.0d0-vbs/phib)
        !          3376:       qbs=twop*czbs*(1.0d0-sarg)
        !          3377:       capbs=czbs/sarg
        !          3378:       go to 520
        !          3379:   510 qbs=czbs*f1+czbsf2*(f3*(vbs-fcpb)
        !          3380:      1   +(vbs*vbs-fcpb2)/(twop+twop))
        !          3381:       capbs=czbsf2*(f3+vbs/twop)
        !          3382:   520 if (vbd.ge.fcpb) go to 530
        !          3383:       sarg=dsqrt(1.0d0-vbd/phib)
        !          3384:       qbd=twop*czbd*(1.0d0-sarg)
        !          3385:       capbd=czbd/sarg
        !          3386:       go to 540
        !          3387:   530 qbd=czbd*f1+czbdf2*(f3*(vbd-fcpb)
        !          3388:      1   +(vbd*vbd-fcpb2)/(twop+twop))
        !          3389:       capbd=czbdf2*(f3+vbd/twop)
        !          3390: c.. bulk and channel charge (plus overlaps)
        !          3391:   540 qgd=covlgd*(vgb+vbd)
        !          3392:       qgs=covlgs*(vgb+vbs)
        !          3393:       qgb=covlgb*vgb
        !          3394:       qg(lx0+loct)=qgate+qgb+qgd+qgs
        !          3395:       qd(lx0+loct)=qchan*0.5d0-qgd-qbd
        !          3396:       qb(lx0+loct)=qbulk+qbd+qbs-qgb
        !          3397: c
        !          3398: c  store small-signal parameters
        !          3399: c
        !          3400:   590 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 900
        !          3401:       if (initf.ne.4) go to 600
        !          3402:       value(lx0+loct)=didvg
        !          3403:       value(lx0+loct+1)=didvd
        !          3404:       value(lx0+loct+2)=didvs
        !          3405:       value(lx0+loct+3)=geqbd
        !          3406: c.. cdrain is used in printing as well as noise calculation
        !          3407:       value(lx0+loct+4)=cdrain
        !          3408:       value(lx0+loct+5)=geqbs
        !          3409: c..   (loct+6 not used)
        !          3410: c.. this is the 'gm' term used in the noise calculation
        !          3411:       value(lx0+loct+7)=didvg
        !          3412:       value(lx0+loct+8)=ccgg
        !          3413:       value(lx0+loct+9)=ccgd
        !          3414:       value(lx0+loct+10)=ccgs
        !          3415:       value(lx0+loct+11)=ccbg
        !          3416:       value(lx0+loct+12)=ccbd
        !          3417:       value(lx0+loct+13)=ccbs
        !          3418:       value(lx0+loct+14)=capbd
        !          3419:       value(lx0+loct+15)=capbs
        !          3420:       go to 1000
        !          3421: c
        !          3422: c  transient analysis
        !          3423: c
        !          3424:   600 if (initf.ne.5) go to 610
        !          3425:       qb(lx1+loct)=qb(lx0+loct)
        !          3426:       qg(lx1+loct)=qg(lx0+loct)
        !          3427:       qd(lx1+loct)=qd(lx0+loct)
        !          3428: c.. integrate qb
        !          3429:   610 call intgr8(geq,ceq,0.0d0,loct+12)
        !          3430: c.. integrate qg
        !          3431:       call intgr8(geq,ceq,0.0d0,loct+14)
        !          3432: c.. integrate qd
        !          3433:       call intgr8(geq,ceq,0.0d0,loct+16)
        !          3434: c.. divvey up the channel charge 50/50 to source and drain
        !          3435: c.. note that symmetry also precludes need for 'devmod' decisions
        !          3436:       gccg2=-0.5d0*(ccgg+ccbg)*ag(1)
        !          3437:       gccd2=-0.5d0*(ccgd+ccbd)*ag(1)
        !          3438:       gccs2=-0.5d0*(ccgs+ccbs)*ag(1)
        !          3439:       gccdg=gccg2+didvg-covlgd*ag(1)
        !          3440:       gccdd=gccdd+gccd2+didvd+(capbd+covlgd)*ag(1)
        !          3441:       gccds=gccs2+didvs
        !          3442:       gccsg=gccg2-didvg-covlgs*ag(1)
        !          3443:       gccsd=gccd2-didvd
        !          3444:       gccss=gccss+gccs2-didvs+(capbs+covlgs)*ag(1)
        !          3445:       gccgg=(ccgg+covlgd+covlgs+covlgb)*ag(1)
        !          3446:       gccgd=(ccgd-covlgd)*ag(1)
        !          3447:       gccgs=(ccgs-covlgs)*ag(1)
        !          3448:       gccbg=(ccbg-covlgb)*ag(1)
        !          3449:       gccbd=gccbd+(ccbd-capbd)*ag(1)
        !          3450:       gccbs=gccbs+(ccbs-capbs)*ag(1)
        !          3451:       cgate=cqg(lx0+loct)
        !          3452:       cbulk=cbulk+cqb(lx0+loct)
        !          3453:       cdrain=cdrain+cqd(lx0+loct)
        !          3454:       if (initf.ne.5) go to 700
        !          3455:       cqb(lx1+loct)=cqb(lx0+loct)
        !          3456:       cqg(lx1+loct)=cqg(lx0+loct)
        !          3457:       cqd(lx1+loct)=cqd(lx0+loct)
        !          3458: c
        !          3459: c  check convergence
        !          3460: c
        !          3461:   700 if (initf.ne.3) go to 710
        !          3462:       if (ioff.ne.0) go to 750
        !          3463:   710 if (icheck.eq.1) go to 720
        !          3464: c     tol=reltol*dmax1(dabs(cdhat),dabs(cdrain))+abstol
        !          3465: c     if (dabs(cdhat-cdrain).ge.tol) go to 720
        !          3466: c     tol=reltol*dmax1(dabs(cbhat),dabs(cbulk))+abstol
        !          3467: c     if (dabs(cbhat-cbulk).le.tol) go to 750
        !          3468:       tol=reltol*dabs(vgb)+vntol
        !          3469:       if(dabs(delvgb).ge.tol) go to 720
        !          3470:       tol=reltol*dabs(vbd)+vntol
        !          3471:       if(dabs(delvbd).ge.tol) go to 720
        !          3472:       tol=reltol*dabs(vbs)+vntol
        !          3473:       if(dabs(delvbs).lt.tol) go to 750
        !          3474:   720 noncon=noncon+1
        !          3475:   750 vbdo(lx0+loct)=vbd
        !          3476:       vbso(lx0+loct)=vbs
        !          3477:       vgbo(lx0+loct)=vgb
        !          3478:       cdo(lx0+loct)=cdrain
        !          3479:       cbo(lx0+loct)=cbulk
        !          3480:       gdgo(lx0+loct)=gccdg
        !          3481:       gddo(lx0+loct)=gccdd
        !          3482:       gdso(lx0+loct)=gccds
        !          3483:       gbgo(lx0+loct)=gccbg
        !          3484:       gbdo(lx0+loct)=gccbd
        !          3485:       gbso(lx0+loct)=gccbs
        !          3486: c
        !          3487: c  load current vector
        !          3488: c
        !          3489:   900 ceqg=type*(cgate-gccgg*vgb+gccgd*vbd+gccgs*vbs)
        !          3490:       ceqb=type*(cbulk-gccbg*vgb+gccbd*vbd+gccbs*vbs)
        !          3491:       ceqd=type*(cdrain-gccdg*vgb+gccdd*vbd+gccds*vbs)
        !          3492:       value(lvn+node2)=value(lvn+node2)-ceqg
        !          3493:       value(lvn+node4)=value(lvn+node4)-ceqb
        !          3494:       value(lvn+node5)=value(lvn+node5)-ceqd
        !          3495:       value(lvn+node6)=value(lvn+node6)+ceqd+ceqg+ceqb
        !          3496: c
        !          3497: c  load y matrix
        !          3498: c
        !          3499:       locy=lynl+nodplc(loc+27)
        !          3500:       value(locy)=value(locy)+gdpr
        !          3501:       locy=lynl+nodplc(loc+28)
        !          3502:       value(locy)=value(locy)+gccgg
        !          3503:       locy=lynl+nodplc(loc+29)
        !          3504:       value(locy)=value(locy)+gspr
        !          3505:       locy=lynl+nodplc(loc+30)
        !          3506:       value(locy)=value(locy)-gccbg-gccbd-gccbs
        !          3507:       locy=lynl+nodplc(loc+31)
        !          3508:       value(locy)=value(locy)+gdpr+gccdd
        !          3509:       locy=lynl+nodplc(loc+32)
        !          3510:       value(locy)=value(locy)+gspr+gccss
        !          3511:       locy=lynl+nodplc(loc+10)
        !          3512:       value(locy)=value(locy)-gdpr
        !          3513:       locy=lynl+nodplc(loc+11)
        !          3514:       value(locy)=value(locy)-gccgg-gccgd-gccgs
        !          3515:       locy=lynl+nodplc(loc+12)
        !          3516:       value(locy)=value(locy)+gccgd
        !          3517:       locy=lynl+nodplc(loc+13)
        !          3518:       value(locy)=value(locy)+gccgs
        !          3519:       locy=lynl+nodplc(loc+14)
        !          3520:       value(locy)=value(locy)-gspr
        !          3521:       locy=lynl+nodplc(loc+15)
        !          3522:       value(locy)=value(locy)+gccbg
        !          3523:       locy=lynl+nodplc(loc+16)
        !          3524:       value(locy)=value(locy)+gccbd
        !          3525:       locy=lynl+nodplc(loc+17)
        !          3526:       value(locy)=value(locy)+gccbs
        !          3527:       locy=lynl+nodplc(loc+18)
        !          3528:       value(locy)=value(locy)-gdpr
        !          3529:       locy=lynl+nodplc(loc+19)
        !          3530:       value(locy)=value(locy)+gccdg
        !          3531:       locy=lynl+nodplc(loc+20)
        !          3532:       value(locy)=value(locy)-gccdg-gccdd-gccds
        !          3533:       locy=lynl+nodplc(loc+21)
        !          3534:       value(locy)=value(locy)+gccds
        !          3535:       locy=lynl+nodplc(loc+22)
        !          3536:       value(locy)=value(locy)+gccsg
        !          3537:       locy=lynl+nodplc(loc+23)
        !          3538:       value(locy)=value(locy)-gspr
        !          3539:       locy=lynl+nodplc(loc+24)
        !          3540:       value(locy)=value(locy)-gccsg-gccsd-gccss
        !          3541:       locy=lynl+nodplc(loc+25)
        !          3542:       value(locy)=value(locy)+gccsd
        !          3543:  1000 loc=nodplc(loc)
        !          3544:       go to 10
        !          3545:       end
        !          3546:       subroutine calcq(vgb,vbd,vbs,qg,qc,qb,
        !          3547:      1  ccgg,ccgd,ccgs,ccbg,ccbd,ccbs,
        !          3548:      2  didvg,didvd,didvs)
        !          3549:       implicit double precision (a-h,o-z)
        !          3550:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
        !          3551:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
        !          3552:      2   itemno,nosolv,ipostp,iscrch
        !          3553:       common /mosarg/ gamma,beta,vto,phi,cox,vbi,xnfs,xnsub,xd,xj,xl,
        !          3554:      1   xlamda,utra,uexp,vbp,von,vdsat,theta,vcrit,vtra,gleff,cdrain
        !          3555:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
        !          3556:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
        !          3557:       iflag=1
        !          3558:       if(mode.ne.1) go to 5
        !          3559:       iflag=0
        !          3560:       if(modedc.eq.2.and.nosolv.ne.0) iflag=1
        !          3561:       if(initf.eq.4) iflag=1
        !          3562:     5 vd=dmax1(phi-vbd,1.0d-8)
        !          3563:       vg=vgb-vbi+phi
        !          3564:       vs=dmax1(phi-vbs,1.0d-8)
        !          3565:       vsp5=dsqrt(vs)
        !          3566:       gammad=gamma
        !          3567:       if(gamma.eq.0.0d0) go to 7
        !          3568:       if(xj.eq.0.0d0) go to 7
        !          3569:       arg=dsqrt(1.0d0+xd*2.0d0*vsp5/xj)
        !          3570:       gfact=1.0d0-xj/xl*(arg-1.0d0)
        !          3571:       gfact=dmax1(0.5d0,gfact)
        !          3572:       gammad=gamma*gfact
        !          3573:     7 vth=gammad*vsp5+vs
        !          3574: c.. von is referenced to vgb for 'fetlim'
        !          3575: c.. change mosfet to reference to vgs (von=vth+vbi-vs) for
        !          3576: c.. printing
        !          3577:       von=vth+vbi-phi
        !          3578:       vdsat=0.0d0
        !          3579:       if(vg.lt.vth) go to 100
        !          3580: c
        !          3581: c  'on' region (linear and saturated)
        !          3582: c
        !          3583:       gamma2=gammad*0.5d0
        !          3584:       sqarg=dsqrt(gamma2*gamma2+vg)
        !          3585:       vsat=(sqarg-gamma2)**2
        !          3586:       vsatcl=vsat
        !          3587:       vs2=vs*vs
        !          3588:       vs3=vs2*vs
        !          3589:       vs5=vs3*vs2
        !          3590:       vs1p5=vs*vsp5
        !          3591:       vs2p5=vs1p5*vs
        !          3592:       if(vcrit.eq.0.0d0) go to 9
        !          3593: c...... iterate to new vsat
        !          3594:       iter=1
        !          3595:     8 ve2=vsat*vsat
        !          3596: c     write(6,8878) iter,vsat
        !          3597:  8878 format(' iter vsat ',i4,1pd11.2)
        !          3598:       vep5=dsqrt(vsat)
        !          3599:       ve1p5=vsat*vep5
        !          3600:       arg2=0.5d0*(ve2-vs2)
        !          3601:       arg1p5=gammad*(ve1p5-vs1p5)/1.5d0
        !          3602:       cdrain=vg*(vsat-vs)-arg1p5-arg2
        !          3603:       didve=vg-gammad*vep5-vsat
        !          3604:       d2idve=-0.5d0*gammad/dmax1(vep5,1.0d-5)-1.0d0
        !          3605:       if(vtra.eq.0.0d0) go to 88
        !          3606:       trafac=1.0d0/(1.0d0+(vsat-vs)/vtra)
        !          3607:       dtrdve=-trafac*trafac/vtra
        !          3608:       d2idve=d2idve*trafac+(dtrdve+dtrdve)*(didve-cdrain*trafac/vtra)
        !          3609:       didve=didve*trafac+dtrdve*cdrain
        !          3610:       cdrain=cdrain*trafac
        !          3611:    88 delv=(didve*vcrit-cdrain)/dabs(didve-vcrit*d2idve)
        !          3612: c.. limit voltage excursion to 1/2 old vsat
        !          3613:       if(dabs(delv).gt.0.5d0*vsat) delv=vsat*dsign(0.5d0,delv)
        !          3614:       vsat=vsat+delv
        !          3615: c     vsat=dmax1(vsat,1.0d-5)
        !          3616: c     vsat=dmin1(vsat,vsatcl)
        !          3617:       if(dabs(delv).lt.1.0d-6) go to 9
        !          3618:       iter=iter+1
        !          3619:       if(iter.gt.20) write(6,7777) vg,vs,vsat,delv
        !          3620:  7777 format(' iteration count for vsat hit limit of 20'/,
        !          3621:      1  ' vg vs vsat delv ',1p4d11.3)
        !          3622:       if(iter.gt.20) go to 9
        !          3623:       go to 8
        !          3624: c.. end of iteration loop
        !          3625: c.. vdsat is referenced to vds for printing only
        !          3626:     9 vdsat=vsat-vs
        !          3627:       if(vsat.gt.vsatcl) write(6,9989) vsat,vsatcl
        !          3628:  9989 format(' ********error****** vsat is larger than classical vsat',/
        !          3629:      1' vsat ',1pd11.3,' classical vsat ',d11.3)
        !          3630:  9999 format(' vsat ',1pd11.3,' classical vsat ',d11.3)
        !          3631:       if(vd.ge.vsat) go to 10
        !          3632: c.. linear region
        !          3633:       ve=vd
        !          3634:       dvedvd=1.0d0
        !          3635:       dvedvg=0.0d0
        !          3636:       go to 15
        !          3637: c.. saturated region
        !          3638:    10 ve=vsat
        !          3639:       dvedvd=0.0d0
        !          3640: c**************** zero dvedvg!!!
        !          3641: c     dvedvg=1.0d0-gamma2/sqarg
        !          3642:       dvedvg=0.0d0
        !          3643: c
        !          3644:    15 ve2=ve*ve
        !          3645:       ve3=ve2*ve
        !          3646:       ve5=ve3*ve2
        !          3647:       vep5=dsqrt(ve)
        !          3648:       ve1p5=ve*vep5
        !          3649:       ve2p5=ve1p5*ve
        !          3650:       arg2=0.5d0*(ve2-vs2)
        !          3651:       arg1p5=gammad*(ve1p5-vs1p5)/1.5d0
        !          3652:       cdrain=vg*(ve-vs)-arg1p5-arg2
        !          3653:       didve=vg-gammad*vep5-ve
        !          3654:       didvg=ve-vs+didve*dvedvg
        !          3655:       didvs=-vg+gammad*vsp5+vs
        !          3656:       if(iflag.eq.0) go to 30
        !          3657:       if(dabs(cdrain).gt.1.0d-5) go to 20
        !          3658: c .. special case when ve almost equals vs and regular formulas don't work
        !          3659: c     write(6,5475) time,cdrain
        !          3660:  5475 format(' time = ',1pd15.6,' cdrain = ',e11.3)
        !          3661:       qg=cox*(vg-vs)
        !          3662:       ccgg=cox
        !          3663:       ccgd=-0.5d0*cox
        !          3664:       ccgs=ccgd
        !          3665:       qb=-cox*gammad*vsp5
        !          3666:       ccbg=0.0d0
        !          3667:       ccbd=-cox*0.25d0*gammad/dmax1(vsp5,1.0d-2)
        !          3668:       ccbs=ccbd
        !          3669:       go to 30
        !          3670: c
        !          3671:    20 arg2p5=gammad*0.4d0*(ve2p5-vs2p5)
        !          3672:       varg=(vg*arg2-arg2p5-(ve3-vs3)/3.0d0)/cdrain
        !          3673:       qg=cox*(vg-varg)
        !          3674:       dqgdve=cox/cdrain*(varg-ve)*didve
        !          3675:       ccgg=cox*(1.0d0-(arg2-varg*(ve-vs))/cdrain)
        !          3676:      1  +dqgdve*dvedvg
        !          3677:       ccgd=dqgdve*dvedvd
        !          3678:       ccgs=cox/cdrain*(varg-vs)*didvs
        !          3679:       qb=-cox/cdrain*(vg*arg1p5-gammad*gammad*arg2-arg2p5)
        !          3680:       dqbdve=-cox/cdrain*(gammad*vep5+qb/cox)*didve
        !          3681:       ccbd=dqbdve*dvedvd
        !          3682:       ccbs=-cox/cdrain*(gammad*vsp5+qb/cox)*didvs
        !          3683:       ccbg=-cox/cdrain*(arg1p5+qb/cox*(ve-vs))
        !          3684:      1  +dqbdve*dvedvg
        !          3685: c.. mobility factor (a-la bdm)
        !          3686:    30 if(uexp.eq.0.0d0) go to 35
        !          3687:       vdenom=vg-vth-utra*(ve-vs)
        !          3688:       if(vdenom.le.vbp) go to 45
        !          3689:       arg=vbp/vdenom
        !          3690:       ufact=dexp(uexp*dlog(arg))
        !          3691:       dcoef=-uexp*ufact*arg/vbp
        !          3692:       didvg=ufact*didvg+cdrain*dcoef
        !          3693:       didvs=ufact*didvs-cdrain*dcoef*(0.5d0*gammad/vsp5+1.0d0-utra)
        !          3694:       didve=ufact*didve-cdrain*dcoef*utra
        !          3695:       cdrain=cdrain*ufact
        !          3696:       go to 45
        !          3697: c.. lateral field effects 
        !          3698:    35 if(vtra.eq.0.0d0) go to 40
        !          3699:       trafac=1.0d0/(1.0d0+(ve-vs)/vtra)
        !          3700:       dtrdve=-trafac*trafac/vtra
        !          3701:       didve=didve*trafac+dtrdve*cdrain
        !          3702: c.. note that dtrdvs=-dtrdve
        !          3703:       didvs=didvs*trafac-dtrdve*cdrain
        !          3704:       didvg=didvg*trafac
        !          3705:       cdrain=cdrain*trafac
        !          3706: c.. mobility variation a-la sun-daseking
        !          3707:    40 if(theta.eq.0.0d0) go to 45
        !          3708:       ufact=1.0d0/(1.0d0+theta*(vg-vth))
        !          3709:       dufact=-theta*ufact*ufact
        !          3710:       didve=ufact*didve
        !          3711:       didvg=ufact*didvg+cdrain*dufact
        !          3712:       didvs=ufact*didvs-cdrain*dufact*(0.5d0*gammad/vsp5+1.0d0)
        !          3713:       cdrain=cdrain*ufact
        !          3714: c .. done with 've', use it
        !          3715:    45 didvd=didve*dvedvd
        !          3716: c.. channel length modulation
        !          3717:       if(vcrit.ne.0.0d0) go to 80
        !          3718:       if (xlamda.gt.0.0d0) go to 50
        !          3719:       if (xnsub.eq.0.0d0) go to 50
        !          3720: c.. frohman-grove (lousy) formulation modified a-la newton
        !          3721:       arg1=(vd-vsat)/4.0d0
        !          3722:       arg2=dsqrt(1.0d0+arg1*arg1)
        !          3723:       arg3=dsqrt(arg1+arg2)
        !          3724:       clfact=1.0d0/(1.0d0-xd/xl*arg3)
        !          3725:       if(clfact.le.0.0d0) go to 60
        !          3726:       dclfct=0.125d0*clfact*clfact*xd/xl*(1.0d0+arg1/arg2)/arg3
        !          3727:       didvd=clfact*didvd+cdrain*dclfct
        !          3728:       didvg=clfact*didvg-cdrain*dclfct*(1.0d0-gamma2/sqarg)
        !          3729:       didvs=clfact*didvs
        !          3730:       cdrain=cdrain*clfact
        !          3731:       go to 200
        !          3732: c.. simple (1+vds*lambda/l) formulation
        !          3733:    50 xlfact=xlamda/xl
        !          3734:       clfact=1.0d0+xlfact*(vd-vs)
        !          3735:       didvd=clfact*didvd+cdrain*xlfact
        !          3736:       didvs=clfact*didvs-cdrain*xlfact
        !          3737:       didvg=clfact*didvg
        !          3738:       cdrain=cdrain*clfact
        !          3739:       go to 200
        !          3740: c.. device punched thru
        !          3741:    60 clfact=1000.0d0
        !          3742:       if(ipunch.gt.50) go to 200
        !          3743:       ipunch=ipunch+1
        !          3744:       write(6,61)
        !          3745:    61 format('0warning:  channel length reduced to zero in mosfet')
        !          3746:       go to 200
        !          3747: c.. into saturation with vcrit ne 0
        !          3748:    80 if(vd.le.vsat) go to 200
        !          3749:       xk1=vcrit/2.0d0/xl/gleff
        !          3750:       temp=dsqrt(xk1*xk1+(vd-vsat)/gleff)
        !          3751:       clfact=1.0d0+(temp-xk1)/xl
        !          3752:       dclfct=0.5d0/xl/temp/gleff
        !          3753:       didvd=didvd*clfact+cdrain*dclfct
        !          3754:       didvs=didvs*clfact
        !          3755:       didvg=didvg*clfact
        !          3756:       cdrain=cdrain*clfact
        !          3757: c
        !          3758:       go to 200
        !          3759: c
        !          3760: c.. cut-off region (vg<vth)
        !          3761: c
        !          3762:   100 continue
        !          3763:       cdrain=0.0d0
        !          3764:       didvg=0.0d0
        !          3765:       didvd=0.0d0
        !          3766:       didvs=0.0d0
        !          3767:       if(iflag.eq.0) return
        !          3768:       if(vg.gt.0.0d0) go to 120
        !          3769:       qg=cox*vg
        !          3770:       ccgg=cox
        !          3771:       go to 130
        !          3772:   120 gamma2=gammad*0.5d0
        !          3773:       sqarg=dsqrt(gamma2*gamma2+vg)
        !          3774:       qg=gammad*cox*(sqarg-gamma2)
        !          3775:       ccgg=0.5d0*cox*gammad/sqarg
        !          3776:   130 qb=-qg
        !          3777:       ccbg=-ccgg
        !          3778:       ccgd=0.0d0
        !          3779:       ccgs=0.0d0
        !          3780:       ccbd=0.0d0
        !          3781:       ccbs=0.0d0
        !          3782:   200 qc=-(qg+qb)
        !          3783: c     write(6,7657) time,vg,vs,ve,qg,ccgg,ccgd,ccgs,qb,ccbg,ccbd,ccbs
        !          3784:  7657 format(' time = ',1pd12.5,' vg vs ve ',3d11.3,
        !          3785:      1 /,' qg ',4d11.3,' qb ',4d11.3)
        !          3786:       return
        !          3787:       end
        !          3788:       subroutine gasfet(loc)
        !          3789: c     *** a gasfet (or any other) model may ***
        !          3790: c     *   be inserted here... if you happen   *
        !          3791: c     *   to have a good one!                 *
        !          3792: c     ***                                   ***
        !          3793: c
        !          3794:       return
        !          3795:       end

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.