![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dctrans.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [CSRG BSD Unix] / 3BSD / cmd / spice|
Return to dctrans.f CVS log | Up to [CSRG BSD Unix] / 3BSD / cmd / spice |
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.