Annotation of 3BSD/cmd/spice/roots.f, revision 1.1.1.1

1.1       root        1:        program spice
                      2:       implicit double precision (a-h,o-z)
                      3: c
                      4: c
                      5: c
                      6: c         *** version VAX UNIX 2X.x (19aug79)
                      7: c           developed from
                      8: c         *** version 2d.2 (26sep76) ***
                      9: c         *** version hp 2.0 (6dec77) ***
                     10: c
                     11: c   by       dick dowell  - hewlett packard company
                     12: c            richard newton - uc berkeley
                     13: c
                     14: c     spice is an electronic circuit simulation program that was deve-
                     15: c loped by the integrated circuits group of the electronics research
                     16: c laboratory and the department of electrical engineering and computer
                     17: c sciences at the university of california, berkeley, california.  the
                     18: c program spice is available free of charge to any interested party.
                     19: c the sale, resale, or use of this program for profit without the
                     20: c express written consent of the department of electrical engineering
                     21: c and computer sciences, university of california, berkeley, california,
                     22: c is forbidden.
                     23: c
                     24: c
                     25: c implementation notes:
                     26: c
                     27: c     subroutines mclock and mdate return the time (as hh:mm:ss) and
                     28: c the date (as dd mmm yy), respectively.  subroutine getcje returns in
                     29: c common block /cje/ various attributes of the current job environment.
                     30: c spice expects getcje to set /cje/ variables maxtim, itime, and icost.
                     31: c maxtim is the maximum cpu time in seconds, itime is the elapsed cpu
                     32: c time in seconds, and icost is the job cost in cents.
                     33: c subroutine memory is used to change the number of memory words
                     34: c allocated to spice.  if the amount of memory allocated to a jobstep
                     35: c is fixed, subroutine memory need not be changed.
                     36: c     ifamwa (set in a data statement below) should be set to the
                     37: c address of the first available word of memory (following overlays, if
                     38: c any).  the proper value should be easily obtainable from any load map.
                     39: c     with the exception of most flags, all data in spice is stored in
                     40: c the form of managed tables allocated in the /blank/ array value().
                     41: c     spice is particularly well-suited to being run using a one-level
                     42: c overlay structure beginning with routines spice (the overlay root),
                     43: c readin, errchk, setup, dctran, dcop, acan, and ovtpvt.  the order of
                     44: c the routines in this listing corresponds to that structure.  note
                     45: c that if cdc-style overlay is to be used, an overlay directive card
                     46: c must be inserted before the first line of each of the just-named
                     47: c routines.
                     48: c
                     49: c
                     50:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
                     51:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
                     52:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
                     53:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
                     54:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
                     55:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
                     56:       common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
                     57:      1  defas,rstats(50),iwidth,lwidth,nopage
                     58:       common /line/ achar,afield(15),oldlin(15),kntrc,kntlim
                     59:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
                     60:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
                     61:       common /mosarg/ gamma,beta,vto,phi,cox,vbi,xnfs,xnsub,xd,xj,xl,
                     62:      1   xlamda,utra,uexp,vbp,von,vdsat,theta,vcrit,vtra,gleff,cdrain
                     63:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
                     64:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
                     65:      2   itemno,nosolv,ipostp,iscrch
                     66:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
                     67:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
                     68:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
                     69:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
                     70:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                     71:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                     72:      2   nwd8,nwd16
                     73:       common /dc/ tcstar(2),tcstop(2),tcincr(2),icvflg,itcelm(2),kssop,
                     74:      1   kinel,kidin,kovar,kidout
                     75:       common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq,
                     76:      1   inoise,nosprt,nosout,nosin,idist,idprt
                     77:       common /tran/ tstep,tstop,tstart,delmax,tdmax,forfre,jtrflg
                     78:       common /outinf/ xincr,string(15),xstart,yvar(8),itab(8),itype(8),
                     79:      1   ilogy(8),npoint,numout,kntr,numdgt
                     80:       common /cje/ maxtim,itime,icost
                     81: c.. common for putwds
                     82:       common /blank/ value(50000)
                     83:       integer nodplc(64)
                     84:       complex*16 cvalue(32)
                     85:       equivalence (value(1),nodplc(1),cvalue(1))
                     86:       integer*2 istats(200)
                     87:       real*4 r4stat(100)
                     88:       equivalence (rstats(1),istats(1),r4stat(1))
                     89:       integer*2 inknt
                     90: c
                     91:       dimension acctit(4)
                     92:       dimension remain(4)
                     93:       data amrje /4hmrje/
                     94:       data ablnk /1h  /
                     95:       data acctit / 8hjob stat, 8histics s, 8hummary  , 8h         /
                     96:       data ahdr1, ahdr2, ahdr3 / 8h spice 2,8h.Xx (vax,8h11/780)    /
                     97: c
                     98: c
                     99:       ipostp=0
                    100:        maxtim=1e8
                    101:        maxlin=50000
                    102:        icost=0
                    103:        ilines=0
                    104: c
                    105: c  initialization
                    106: c
                    107:       aprog(1)=ahdr1
                    108:       aprog(2)=ahdr2
                    109:       aprog(3)=ahdr3
                    110:       achar=ablnk
                    111:       keof=0
                    112:       call mclock(atime)
                    113:       call mdate(adate)
                    114:       boltz=1.3806226d-23
                    115:       charge=1.6021918d-19
                    116:       ctok=273.15d0
                    117:       eps0=8.854214871d-14
                    118:       epssil=11.7d0*eps0
                    119:       epsox=3.9d0*eps0
                    120:       twopi=8.0d0*datan2(1.0d0,1.0d0)
                    121:       rad=360.0d0/twopi
                    122:       xlog2=dlog(2.0d0)
                    123:       xlog10=dlog(10.0d0)
                    124:       root2=dsqrt(2.0d0)
                    125:       nodata=1
                    126: c
                    127: c  begin job
                    128: c
                    129:    10 if (keof.eq.1) go to 1000
                    130:       call getcje
                    131:       call second(time1)
                    132:       icost1=icost
                    133:       igoof=0
                    134:       mode=0
                    135:       nogo=0
                    136:       maxmem=100000
                    137:       call setmem(nodplc(1),maxmem)
                    138:       if (nogo.ne.0) go to 1000
                    139:       call zero8(rstats,50)
                    140: c
                    141: c  read remainder of data deck and check for input errors
                    142: c
                    143:       call readin
                    144:       if (nogo.ne.0) go to 300
                    145:       if (keof.eq.1) go to 1000
                    146:       nodata=0
                    147:       call errchk
                    148:       if (nogo.ne.0) go to 300
                    149:       call setup
                    150:       if (nogo.ne.0) go to 300
                    151: c
                    152: c  cycle through temperatures
                    153: c
                    154:       itemno=1
                    155:       if (numtem.eq.1) go to 110
                    156:   100 if (itemno.eq.numtem) go to 320
                    157:       itemno=itemno+1
                    158:       call tmpupd
                    159: c
                    160: c  dc transfer curves
                    161: c
                    162:   110 if (icvflg.eq.0) go to 150
                    163: c...  see routine *dctran* for explanation of *mode*, etc.
                    164:       mode=1
                    165:       modedc=3
                    166:       call dctran
                    167:       call ovtpvt
                    168:       if (nogo.ne.0) go to 300
                    169: c
                    170: c  small signal operating point
                    171: c
                    172:   150 if (kssop.gt.0) go to 170
                    173:       if (jacflg.ne.0) go to 170
                    174:       if ((icvflg+jtrflg).gt.0) go to 250
                    175:   170 mode=1
                    176:       modedc=1
                    177:       call dctran
                    178:       if (nogo.ne.0) go to 300
                    179:       call dcop
                    180:       if (nogo.ne.0) go to 300
                    181: c
                    182: c  ac small signal analysis
                    183: c
                    184:   200 if (jacflg.eq.0) go to 250
                    185:       mode=3
                    186:       call acan
                    187:       call ovtpvt
                    188:       if (nogo.ne.0) go to 300
                    189: c
                    190: c  transient analysis
                    191: c
                    192:   250 if (jtrflg.eq.0) go to 100
                    193:       mode=1
                    194:       modedc=2
                    195:       call dctran
                    196:       if (nogo.ne.0) go to 300
                    197:       call dcop
                    198:       if (nogo.ne.0) go to 300
                    199:       mode=2
                    200:       call dctran
                    201:       call ovtpvt
                    202:       if (nogo.ne.0) go to 300
                    203:       go to 100
                    204: c
                    205: c  job concluded
                    206: c
                    207:   300 write (6,301)
                    208:   301 format(1h0,9x,'***** job aborted')
                    209:       nodata=0
                    210: c
                    211: c  job accounting
                    212: c
                    213:   320 continue
                    214:       numel=0
                    215:       do 360 i=1,18
                    216:   360 numel=numel+jelcnt(i)
                    217:       numtem=max0(numtem-1,1)
                    218:       idist=min0(idist,1)
                    219:       if (iprnta.eq.0) go to 800
                    220:       call title(-1,lwidth,1,acctit)
                    221:       write (6,361) nunods,ncnods,numnod,numel,(jelcnt(i),i=11,14)
                    222:   361 format('   nunods ncnods numnod numel  diodes  bjts  jfets  mfets'
                    223:      1   //,i9,2i7,i6,i8,i6,2i7)
                    224:       write (6,371) numtem,icvflg,jtrflg,jacflg,inoise,idist,nogo
                    225:   371 format(/'0  numtem icvflg jtrflg jacflg inoise  idist   nogo'/,
                    226:      1   2h0 ,7i7)
                    227:       write (6,381) rstats(20),rstats(21),rstats(22),rstats(23),
                    228:      1   rstats(26),rstats(27)
                    229:   381 format(/'0  nstop   nttbr   nttar   ifill    iops    perspa'//,
                    230:      1   1x,5f8.0,f9.3)
                    231:       write (6,391) rstats(30),rstats(31),rstats(32),maxmem,maxuse,
                    232:      1   cpyknt
                    233:   391 format(/'0  numttp  numrtp  numnit  maxmem  memuse  copyknt',//,
                    234:      1   2x,3f8.0,2x,i6,2x,i6,2x,f8.0)
                    235:       write (6,401) (rstats(i),i=1,11)
                    236:   401 format(/,
                    237:      1   1h0,9x,'readin  ',12x,f10.2/,
                    238:      2   1h0,9x,'setup   ',12x,f10.2/,
                    239:      3   1h0,9x,'trcurv  ',12x,f10.2,10x,f6.0/,
                    240:      4   1h0,9x,'dcan    ',12x,f10.2,10x,f6.0/,
                    241:      5   1h0,9x,'acan    ',12x,f10.2,10x,f6.0/,
                    242:      6   1h0,9x,'tranan  ',12x,f10.2,10x,f6.0/,
                    243:      7   1h0,9x,'output  ',12x,f10.2)
                    244:   800 call getcje
                    245:       call second(time2)
                    246:       et=time2-time1
                    247:       tcost=dfloat(icost-icost1)/100.0d0
                    248:       if (iprnta.eq.0) go to 810
                    249:       ohead=et-(rstats(1)+rstats(2)+rstats(3)+rstats(5)+rstats(7)
                    250:      1   +rstats(9)+rstats(11))
                    251:       write (6,801) ohead
                    252:   801 format(1h0,9x,'overhead',12x,f10.2)
                    253:   810 write (6,811) et
                    254:   811 format(1h0,9x,'total job time      ',f10.2)
                    255:       tcost=tcost*11.5d0/23.0d0
                    256: c      write(6,812) tcost
                    257: c  812 format(1h0,9x,'total job cost       $',f8.2,
                    258: c     1  ' @ $11.50 per cpu minute',
                    259: c     2  /'0this lower rate applies for remainder of fiscal 79')
                    260:       rstats(33)=cpyknt
                    261:       rstats(34)=et
                    262:       rstats(35)=tcost
                    263:       rstats(36)=ohead
                    264: c.. convert dble to sgl - 72/2 is how many to convert
                    265: c     call dblsgl(rstats(1),72)
                    266:       istats(73)=nunods
                    267:       istats(74)=ncnods
                    268:       istats(75)=numnod
                    269:       istats(76)=numel
                    270:       istats(77)=jelcnt(11)
                    271:       istats(78)=jelcnt(12)
                    272:       istats(79)=jelcnt(13)
                    273:       istats(80)=jelcnt(14)
                    274:       istats(81)=numtem
                    275:       istats(82)=icvflg
                    276:       istats(83)=jtrflg
                    277:       istats(84)=jacflg
                    278:       istats(85)=inoise
                    279:       istats(86)=idist
                    280:       istats(87)=nogo
                    281:       istats(88)=maxmem
                    282:       istats(89)=maxuse
                    283: c     do 820 i=1,36
                    284: c 820 r4stat(i)=rlconv(r4stat(i))
                    285: c     call cadend(istats,100)
                    286:   900 if ((maxtim-itime).ge.limtim) go to 10
                    287:       write (6,901)
                    288:   901 format('1warning:  further analysis stopped due to cpu time limit'
                    289:      1/)
                    290:  1000 if(nodata.ne.0) write(6,1001)
                    291:  1001 format(/1x,'input deck (file) contains no data.')
                    292:       stop
                    293:       end
                    294:       subroutine tmpupd
                    295:       implicit double precision (a-h,o-z)
                    296: c
                    297: c     this routine updates the temperature-dependent parameters in the
                    298: c device models.  it also updates the values of temperature-dependent
                    299: c resistors.  the updated values are printed.
                    300: c
                    301:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
                    302:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
                    303:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
                    304:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
                    305:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
                    306:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
                    307:       common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
                    308:      1  defas,rstats(50),iwidth,lwidth,nopage
                    309:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
                    310:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
                    311:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
                    312:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
                    313:      2   itemno,nosolv,ipostp,iscrch
                    314:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
                    315:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
                    316:       common /blank/ value(50000)
                    317:       integer nodplc(64)
                    318:       complex*16 cvalue(32)
                    319:       equivalence (value(1),nodplc(1),cvalue(1))
                    320: c
                    321: c
                    322:       dimension tmptit(4)
                    323:       data tmptit / 8htemperat, 8hure-adju, 8hsted val, 8hues      /
                    324: c
                    325: c
                    326:       tempd=value(itemps+itemno)+ctok
                    327:       xkt=boltz*tempd
                    328:       oldvt=vt
                    329:       vt=xkt/charge
                    330:       ratio=tempd/(value(itemps+itemno-1)+ctok)
                    331:       ratlog=dlog(ratio)
                    332:       ratio1=ratio-1.0d0
                    333:       delt=value(itemps+itemno)-value(itemps+1)
                    334:       deltsq=delt*delt
                    335:       reftmp=27.0d0+ctok
                    336:       oldeg=egfet
                    337:       egfet=1.16d0-(7.02d-4*tempd*tempd)/(tempd+1108.0d0)
                    338:       oldxni=xni
                    339:       arg=-egfet/(xkt+xkt)+1.1151d0/(boltz*(reftmp+reftmp))
                    340:       factor=tempd/reftmp
                    341:       factor=factor*dsqrt(factor)
                    342:       xni=1.45d10*factor*dexp(charge*arg)
                    343:       pbfact=(vt+vt)*dlog(oldxni/xni)
                    344:       call title(0,lwidth,1,tmptit)
                    345: c
                    346: c  resistors
                    347: c
                    348:       loc=locate(1)
                    349:       ititle=0
                    350:    10 if (loc.eq.0) go to 100
                    351:       locv=nodplc(loc+1)
                    352:       tc1=value(locv+3)
                    353:       tc2=value(locv+4)
                    354:       if (tc1.ne.0.0d0) go to 20
                    355:       if (tc2.eq.0.0d0) go to 40
                    356:    20 if (ititle.ne.0) go to 30
                    357:       write (6,21)
                    358:    21 format(//'0**** resistors',/,'0name',8x,'value',//)
                    359:       ititle=1
                    360:    30 rnew=value(locv+2)*(1.0d0+tc1*delt+tc2*deltsq)
                    361:       value(locv+1)=1.0d0/rnew
                    362:       write (6,31) value(locv),rnew
                    363:    31 format(1x,a8,1p6d11.2)
                    364:    40 loc=nodplc(loc)
                    365:       go to 10
                    366: c
                    367: c  diode model
                    368: c
                    369:   100 loc=locate(21)
                    370:       if (loc.eq.0) go to 200
                    371:       write (6,101)
                    372:   101 format(//'0**** diode model parameters',/,'0name',9x,'is',9x,'pb',
                    373:      1//)
                    374:   110 if (loc.eq.0) go to 200
                    375:       locv=nodplc(loc+1)
                    376: c...  is(t2)=is(t1)*dexp(eg/(n*vt)*(t2/t1-1))*(t2/t1)^(pt/n)
                    377:       xn=value(locv+3)
                    378:       factor=ratio1*value(locv+8)/(xn*vt)+value(locv+9)/xn*ratlog
                    379:       factor=dexp(factor)
                    380:       value(locv+1)=value(locv+1)*factor
                    381:       oldpb=value(locv+6)
                    382:       value(locv+6)=ratio*oldpb+pbfact
                    383:       value(locv+12)=value(locv+12)*value(locv+6)/oldpb
                    384:       value(locv+15)=value(locv+15)*value(locv+6)/oldpb
                    385:       vte=value(locv+3)*vt
                    386:       value(locv+18)=vte*dlog(vte/(root2*value(locv+1)))
                    387:       write (6,31) value(locv),value(locv+1),value(locv+6)
                    388:       loc=nodplc(loc)
                    389:       go to 110
                    390: c
                    391: c  bipolar transistor model
                    392: c
                    393:   200 loc=locate(22)
                    394:       if (loc.eq.0) go to 300
                    395:       write (6,201)
                    396:   201 format(//'0**** bjt model parameters',/,'0name',9x,'js',8x,'bf ',
                    397:      1   7x,'jle',7x,'br ',7x,'jlc',7x,'vje',7x,'vjc',//)
                    398:   210 if (loc.eq.0) go to 300
                    399:       locv=nodplc(loc+1)
                    400: c...  is(t2)=is(t1)*dexp(eg/vt*(t2/t1-1))*(t2/t1)^pt
                    401:       factln=ratio1*value(locv+42)/vt+value(locv+43)*ratlog
                    402:       factor=dexp(factln)
                    403:       value(locv+1)=value(locv+1)*factor
                    404:       tb=value(locv+41)
                    405:       bfactr=dexp(tb*ratlog)
                    406:       value(locv+2)=value(locv+2)*bfactr
                    407:       value(locv+8)=value(locv+8)*bfactr
                    408:       value(locv+6)=value(locv+6)*dexp(factln/value(locv+7))/bfactr
                    409:       value(locv+12)=value(locv+12)*dexp(factln/value(locv+13))
                    410:      1               /bfactr
                    411:       oldpb=value(locv+22)
                    412:       value(locv+22)=ratio*oldpb+pbfact
                    413:       value(locv+46)=value(locv+46)*value(locv+22)/oldpb
                    414:       value(locv+47)=value(locv+47)*value(locv+22)/oldpb
                    415:       oldpb=value(locv+30)
                    416:       value(locv+30)=ratio*oldpb+pbfact
                    417:       value(locv+50)=value(locv+50)*value(locv+30)/oldpb
                    418:       value(locv+51)=value(locv+51)*value(locv+30)/oldpb
                    419:       value(locv+54)=vt*dlog(vt/(root2*value(locv+1)))
                    420:       write (6,211) value(locv),value(locv+1),value(locv+2),
                    421:      1   value(locv+6),value(locv+8),value(locv+12),value(locv+22),
                    422:      2   value(locv+30)
                    423:   211 format(1x,a8,1p7d10.2)
                    424:       loc=nodplc(loc)
                    425:       go to 210
                    426: c
                    427: c  jfet model
                    428: c
                    429:   300 loc=locate(23)
                    430:       if (loc.eq.0) go to 400
                    431:       write (6,301)
                    432:   301 format(//'0**** jfet model parameters',/,'0name',9x,'is',9x,'pb',
                    433:      1//)
                    434:   310 if (loc.eq.0) go to 400
                    435:       locv=nodplc(loc+1)
                    436:       value(locv+9)=value(locv+9)*dexp(ratio1*1.11d0/vt)
                    437:       oldpb=value(locv+8)
                    438:       value(locv+8)=ratio*oldpb+pbfact
                    439:       value(locv+12)=value(locv+12)*value(locv+8)/oldpb
                    440:       value(locv+13)=value(locv+13)*value(locv+8)/oldpb
                    441:       value(locv+16)=vt*dlog(vt/(root2*value(locv+9)))
                    442:       write (6,31) value(locv),value(locv+9),value(locv+8)
                    443:       loc=nodplc(loc)
                    444:       go to 310
                    445: c
                    446: c  mosfet model
                    447: c
                    448:   400 loc=locate(24)
                    449:       if (loc.eq.0) go to 1000
                    450:       iprnt=1
                    451:   410 if (loc.eq.0) go to 1000
                    452: c.. no temperature effects have been coded for ga-as fets
                    453:       if(nodplc(loc+2).eq.0) go to 430
                    454:       locv=nodplc(loc+1)
                    455:       if(iprnt.ne.0) write (6,401)
                    456:   401 format(//'0**** mosfet model parameters',/,'0name',8x,'vto',8x,
                    457:      1   'phi',9x,'pb',9x,'js',7x,'kp',//)
                    458:       iprnt=0
                    459:       ratio4=ratio*dsqrt(ratio)
                    460:       value(locv+2)=value(locv+2)/ratio4
                    461:       value(locv+23)=value(locv+23)/ratio4
                    462:       oldphi=value(locv+4)
                    463:       value(locv+4)=ratio*oldphi+pbfact
                    464:       phi=value(locv+4)
                    465:       type=nodplc(loc+2)
                    466:       tps=value(locv+22)
                    467:       vfb=value(locv+34)-type*oldphi
                    468:       vstrip=vfb+0.5d0*type*oldphi
                    469:       if(value(locv+21).ne.0.0d0) go to 415
                    470:       vstrip=vstrip+0.5d0*(oldeg-egfet)
                    471:       go to 420
                    472:   415 oldgat=oldvt*dlog(value(locv+21)/oldxni)
                    473:       gatnew=vt*dlog(value(locv+21)/xni)
                    474:       vstrip=vstrip+type*tps*(oldgat-gatnew)
                    475:   420 vfb=vstrip-0.5d0*type*phi
                    476:       value(locv+34)=vfb+type*phi
                    477:       value(locv+1)=value(locv+34)+type*value(locv+3)*dsqrt(phi)
                    478:       value(locv+15)=value(locv+15)*dexp(-egfet/vt+oldeg/oldvt)
                    479:       oldpb=value(locv+14)
                    480:       value(locv+14)=ratio*oldpb+pbfact
                    481:       pb=value(locv+14)
                    482:       ratio2=oldpb/pb
                    483:       ratio3=dsqrt(ratio2)
                    484:       value(locv+11)=value(locv+11)*ratio3
                    485:       value(locv+12)=value(locv+12)*ratio3
                    486:       pbrat=1.0d0/ratio2
                    487:       value(locv+29)=value(locv+29)*pbrat
                    488:       value(locv+30)=value(locv+30)*pbrat
                    489:       write (6,31) value(locv),value(locv+1),value(locv+4),
                    490:      1   value(locv+14),value(locv+15),value(locv+2)
                    491:   430 loc=nodplc(loc)
                    492:       go to 410
                    493: c
                    494: c  finished
                    495: c
                    496:  1000 return
                    497:       end
                    498:       subroutine find(aname,id,loc,iforce)
                    499:       implicit double precision (a-h,o-z)
                    500: c
                    501: c     this routine searches the list with number 'id' for an element
                    502: c with name 'aname'.  loc is set to point to the element.  if iforce is
                    503: c nonzero, then find expects to have to add the element to the list, and
                    504: c reports a fatal error if the element is found.  if subcircuit defini-
                    505: c tion is in progress (nonzero value for nsbckt), then find searches the
                    506: c current subcircuit definition list rather than the nominal element
                    507: c list.
                    508: c
                    509:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
                    510:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
                    511:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
                    512:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
                    513:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
                    514:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
                    515:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
                    516:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
                    517:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
                    518:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
                    519:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                    520:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                    521:      2   nwd8,nwd16
                    522:       common /blank/ value(50000)
                    523:       integer nodplc(64)
                    524:       complex*16 cvalue(32)
                    525:       equivalence (value(1),nodplc(1),cvalue(1))
                    526: c
                    527: c  index to the contents of the various lists:
                    528: c
                    529: c        list      contents
                    530: c        ----      --------
                    531: c
                    532: c          1       resistors
                    533: c          2       nonlinear capacitors
                    534: c          3       nonlinear inductors
                    535: c          4       mutual inductors
                    536: c          5       nonlinear voltage controlled current sources
                    537: c          6       nonlinear voltage controlled voltage sources
                    538: c          7       nonlinear current controlled current sources
                    539: c          8       nonlinear current controlled voltage sources
                    540: c          9       independent voltage sources
                    541: c         10       independent current sources
                    542: c         11       diodes
                    543: c         12       bipolar junction transistors
                    544: c         13       junction field-effect transistors (jfets)
                    545: c         14       metal-oxide-semiconductor junction fets (mosfets)
                    546: c         15       s-parameter 2-port network
                    547: c         16       y-parameter 2-port network
                    548: c         17       transmission lines
                    549: c         18       <unused>
                    550: c         19       subcircuit calls
                    551: c         20       subcircuit definitions
                    552: c         21       diode model
                    553: c         22       bjt model
                    554: c         23       jfet model
                    555: c         24       mosfet model
                    556: c      25-30       <unused>
                    557: c         31       .print dc
                    558: c         32       .print tran
                    559: c         33       .print ac
                    560: c         34       .print noise
                    561: c         35       .print distortion
                    562: c         36       .plot dc
                    563: c         37       .plot tr
                    564: c         38       .plot ac
                    565: c         39       .plot noise
                    566: c         40       .plot distortion
                    567: c         41       outputs for dc
                    568: c         42       outputs for transient
                    569: c         43       outputs for ac
                    570: c         44       outputs for noise
                    571: c         45       outputs for distortion
                    572: c      46-50       <unused>
                    573: c
                    574:       integer xxor
                    575:       dimension lnod(50),lval(50)
                    576:       data lnod / 9,13,15, 7,14,15,14,15,12, 7,
                    577:      1           17,37,26,34, 7, 7,34, 0, 5, 5,
                    578:      2            4, 4, 4, 4, 0, 0, 0, 0, 0, 0,
                    579:      3           21,21,21,21,21,21,21,21,21,21,
                    580:      4            8, 8, 8, 8, 8, 0, 0, 0, 0, 0 /
                    581:       data lval / 5, 4, 4, 2, 1, 1, 1, 1, 4, 4,
                    582:      1            3, 4, 4,13, 1, 1, 9, 0, 1, 1,
                    583:      2           19,55,17,41, 0, 0, 0, 0, 0, 0,
                    584:      3            1, 1, 1, 1, 1,17,17,17,17,17,
                    585:      4            1, 1, 1, 1, 1, 0, 0, 0, 0, 0 /
                    586:       data ndefin /2h.u/
                    587: c
                    588: c
                    589:       anam=aname
                    590:       call sizmem(ielmnt,isize)
                    591:       locn=ielmnt+isize+2
                    592:       if (nsbckt.eq.0) go to 10
                    593:       loct=nodplc(isbckt+nsbckt)
                    594:       loc=nodplc(loct+3)
                    595:       if (loc.ne.0) go to 20
                    596:       nodplc(loct+3)=locn
                    597:       go to 60
                    598:    10 loc=locate(id)
                    599:       if (loc.ne.0) go to 20
                    600:       locate(id)=locn
                    601:       go to 50
                    602: c
                    603: c  search list for a name match
                    604: c
                    605:    20 locv=nodplc(loc+1)
                    606:       if (xxor(anam,value(locv)).ne.0) go to 30
                    607:       if (nsbckt.eq.0) go to 25
                    608:       if (nodplc(loc-1).ne.id) go to 30
                    609:    25 if (nodplc(loc+2).eq.ndefin) go to 200
                    610:       if (iforce.eq.0) go to 200
                    611:       write (6,26) anam
                    612:    26 format('0*error*:  above line attempts to redefine ',a8/)
                    613:       nogo=1
                    614:    30 if (nodplc(loc).eq.0) go to 40
                    615:       loc=nodplc(loc)
                    616:       go to 20
                    617: c
                    618: c  reserve space for this element
                    619: c
                    620:    40 nodplc(loc)=locn
                    621:       if (nsbckt.ne.0) go to 60
                    622:    50 jelcnt(id)=jelcnt(id)+1
                    623:    60 loc=locn
                    624:       itemp=loc+lnod(id)*nwd4-1
                    625:       locv=nxtevn(itemp-1)+1
                    626:       itemp=locv-itemp
                    627:       ktmp=lnod(id)*nwd4+lval(id)*nwd8+itemp
                    628:       call extmem(ielmnt,ktmp)
                    629:       locv=(locv-1)/nwd8+1
                    630:       iptr=0
                    631:       if (nsbckt.eq.0) go to 80
                    632:       iptr=id
                    633:    80 nodplc(loc-1)=iptr
                    634:       nodplc(loc)=0
                    635:       nodplc(loc+1)=locv
                    636:       value(locv)=anam
                    637: c
                    638: c  background storage
                    639: c
                    640:   100 nodplc(loc+2)=ndefin
                    641:       nword=lnod(id)-4
                    642:       if (nword.lt.1) go to 120
                    643:       call zero4(nodplc(loc+3),nword)
                    644:   120 nword=lval(id)-1
                    645:       if (nword.lt.1) go to 200
                    646:       call zero8(value(locv+1),nword)
                    647: c
                    648: c  exit
                    649: c
                    650:   200 return
                    651:       end
                    652:       subroutine title(ifold,len,icom,coment)
                    653:       implicit double precision (a-h,o-z)
                    654: c
                    655: c     this routine writes a title on the output file.  ifold indicates
                    656: c whether the page eject should be to the next concave, convex, or any
                    657: c page fold depending on whether its value is <0, >0, or =0.  the page
                    658: c eject is suppressed (as is much of the heading) if the variable nopage
                    659: c is nonzero.
                    660: c
                    661:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
                    662:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
                    663:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
                    664:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
                    665:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
                    666:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
                    667:       common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
                    668:      1  defas,rstats(50),iwidth,lwidth,nopage
                    669:       common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
                    670:      1   xmu,mode,modedc,icalc,initf,method,iord,maxord,noncon,iterno,
                    671:      2   itemno,nosolv,ipostp,iscrch
                    672:       common /blank/ value(50000)
                    673:       integer nodplc(64)
                    674:       complex*16 cvalue(32)
                    675:       equivalence (value(1),nodplc(1),cvalue(1))
                    676: c
                    677: c
                    678:       dimension coment(4)
                    679: c
                    680: c
                    681:       if(nopage.eq.1) go to 150
                    682: c
                    683:    30 if (len.le.80) go to 100
                    684:       write (6,31) adate,aprog,atime,(atitle(i),i=1,10)
                    685:    31 format(1h1,9(2h* ),a10,1x,11(2h* ),3a8,11(2h* ),a10,9(2h *),//1h0,
                    686:      1   15a8/)
                    687:       if (icom.eq.0) go to 40
                    688:       write (6,36) coment,value(itemps+itemno)
                    689:    36 format(5h0****,17x,4a8,21x,'temperature =',f9.3,' deg c'/)
                    690:    40 write (6,41)
                    691:    41 format(1h0,63(2h* )//)
                    692:       go to 200
                    693: c
                    694: c
                    695:   100 write (6,101) adate,aprog,atime,(atitle(i),i=1,10)
                    696:   101 format(1h1,5(1h*),a10,1x,8(1h*),3a8,8(1h*),a10,5(1h*)//1h0,10a8/)
                    697:       if (icom.eq.0) go to 110
                    698:       write (6,106) coment,value(itemps+itemno)
                    699:   106 format(10h0****     ,4a8,' temperature =',f9.3,' deg c'/)
                    700:   110 write (6,111)
                    701:   111 format(1h0,71(1h*)//)
                    702:       go to 200
                    703: c
                    704: c
                    705:   150 if (icom.eq.0) go to 160
                    706:       write (6,106) coment,value(itemps+itemno)
                    707:       go to 200
                    708:   160 write (6,161) aprog
                    709:   161 format(1h0,3a8,/)
                    710: c
                    711: c  finished
                    712: c
                    713:   200 return
                    714:       end
                    715:       subroutine dcdcmp
                    716:       implicit double precision (a-h,o-z)
                    717: c
                    718: c     this routine performs an in-place lu factorization of the coef-
                    719: c ficient matrix.
                    720: c
                    721:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
                    722:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
                    723:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
                    724:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
                    725:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
                    726:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
                    727:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
                    728:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
                    729:       common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
                    730:      1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,igoof,nogo,keof
                    731:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
                    732:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
                    733:       common /blank/ value(50000)
                    734:       integer nodplc(64)
                    735:       complex*16 cvalue(32)
                    736:       equivalence (value(1),nodplc(1),cvalue(1))
                    737:       data ikount /0/
                    738: c
                    739: c
                    740:       do 100 i=2,nstop
                    741:       io=nodplc(iorder+i)
                    742:       if (dabs(value(lynl+io)).ge.gmin) go to 10
                    743:       value(lynl+io)=gmin
                    744:       igoof=igoof+1
                    745:       if(ikount.gt.20) go to 10
                    746:       ikount=ikount+1
                    747:       if(io.le.nunods) write(6,9) nodplc(junode+io)
                    748:     9 format(' at node ',i5)
                    749:    10 jstart=nodplc(ilc+i)
                    750:       jstop=nodplc(ilc+i+1)-1
                    751:       if (jstart.gt.jstop) go to 100
                    752:       do 90 j=jstart,jstop
                    753:       value(lyl+j)=value(lyl+j)/value(lynl+io)
                    754:       icol=nodplc(ilr+j)
                    755:       kstart=nodplc(iur+i)
                    756:       kstop=nodplc(iur+i+1)-1
                    757:       if (kstart.gt.kstop) go to 90
                    758:       do 80 k=kstart,kstop
                    759:       irow=nodplc(iuc+k)
                    760:       if (icol-irow) 20,60,40
                    761: c
                    762: c  find (icol,irow) matrix term (upper triangle)
                    763: c
                    764:    20 l=nodplc(iur+icol+1)
                    765:    30 l=l-1
                    766:       if (nodplc(iuc+l).ne.irow) go to 30
                    767:       ispot=lyu+l
                    768:       go to 70
                    769: c
                    770: c  find (icol,irow) matrix term (lower triangle)
                    771: c
                    772:    40 l=nodplc(ilc+irow+1)
                    773:    50 l=l-1
                    774:       if (nodplc(ilr+l).ne.icol) go to 50
                    775:       ispot=lyl+l
                    776:       go to 70
                    777: c
                    778: c  find (icol,irow) matrix term (diagonal)
                    779: c
                    780:    60 ispot=lynl+nodplc(iorder+irow)
                    781: c
                    782:    70 value(ispot)=value(ispot)-value(lyl+j)*value(lyu+k)
                    783:    80 continue
                    784:    90 continue
                    785:   100 continue
                    786:       return
                    787:       end
                    788:       subroutine dcsol
                    789:       implicit double precision (a-h,o-z)
                    790: c
                    791: c     this routine solves the system of circuit equations by performing
                    792: c a forward and backward substitution step using the previously-computed
                    793: c lu factors.
                    794: c
                    795:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
                    796:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
                    797:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
                    798:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
                    799:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
                    800:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
                    801:       common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
                    802:      1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs
                    803:       common /blank/ value(50000)
                    804:       integer nodplc(64)
                    805:       complex*16 cvalue(32)
                    806:       equivalence (value(1),nodplc(1),cvalue(1))
                    807: c
                    808: c  forward substitution
                    809: c
                    810:       do 20 i=2,nstop
                    811:       jstart=nodplc(ilc+i)
                    812:       jstop=nodplc(ilc+i+1)-1
                    813:       if (jstart.gt.jstop) go to 20
                    814:       io=nodplc(iorder+i)
                    815:       if (value(lvn+io).eq.0.0d0) go to 20
                    816:       do 10 j=jstart,jstop
                    817:       jo=nodplc(ilr+j)
                    818:       jo=nodplc(iorder+jo)
                    819:       value(lvn+jo)=value(lvn+jo)-value(lyl+j)*value(lvn+io)
                    820:    10 continue
                    821:    20 continue
                    822: c
                    823: c  back substitution
                    824: c
                    825:       k=nstop+1
                    826:       do 50 i=2,nstop
                    827:       k=k-1
                    828:       io=nodplc(iorder+k)
                    829:       jstart=nodplc(iur+k)
                    830:       jstop=nodplc(iur+k+1)-1
                    831:       if (jstart.gt.jstop) go to 40
                    832:       do 30 j=jstart,jstop
                    833:       jo=nodplc(iuc+j)
                    834:       jo=nodplc(iorder+jo)
                    835:       value(lvn+io)=value(lvn+io)-value(lyu+j)*value(lvn+jo)
                    836:    30 continue
                    837:    40 value(lvn+io)=value(lvn+io)/value(lynl+io)
                    838:    50 continue
                    839:       return
                    840:       end
                    841:        subroutine setmem(ipntr,ksize)
                    842:       implicit double precision (a-h,o-z)
                    843: c
                    844: c     this routine performs dynamic memory management.  it is used in
                    845: c     spice2, and useable in any program.
                    846: c
                    847: c     memory is managed within an array selected by the calling program.
                    848: c     one may either dimension this array to the 'maxmem' size, or more
                    849: c     desirably, find the address of the first available word of memory
                    850: c     above your program, and dimension your array to '1'.  passing the
                    851: c     address of the first data word available permits the manager to
                    852: c     use 'illegal' indices into the data area.
                    853: c
                    854: c     this routine must have access to an integer function called 'locf'
                    855: c     which returns the address of its argument.  addresses as used by this
                    856: c     program refer to 'integer' addresses, not byte addresses.
                    857: c
                    858: c entry points:
                    859: c      setmem - set initial memory
                    860: c      getm4  - get block for table of integers
                    861: c      getm8  - get block for table of floating point variables
                    862: c      getm16 - get block for table of complex variables
                    863: c      relmem - release part of block
                    864: c      extmem - extend size of existing block
                    865: c      sizmem - determine size of existing block
                    866: c      clrmem - release block
                    867: c      ptrmem - reset memory pointer
                    868: c      crunch - force memory compaction
                    869: c      avlm4  - amount of space available (integers)
                    870: c      avlm8  - amount of space available (real)
                    871: c      avlm16 - amount of space available (complex)
                    872: c
                    873: c calling sequences:
                    874: c      call setmem(imem(1),maxmem)
                    875: c      call setmem(imem(1),maxmem,kfamwa)  non 3000 machines kfamwa is
                    876: c                                          address of first available word
                    877: c                                          of data
                    878: c      call getm4 (ipntr,blksiz)  where blksize is the number of entries
                    879: c      call getm8 (ipntr,blksiz)
                    880: c      call getm16(ipntr,blksiz)
                    881: c      call relmem(ipntr,relsiz)
                    882: c      call extmem(ipntr,extsiz)  extsiz is the number of entries to be added
                    883: c      call sizmem(ipntr,blksiz)
                    884: c      call clrmem(ipntr)
                    885: c      call ptrmem(ipntr1,ipntr2)
                    886: c      call avlm4(ispace)
                    887: c      call avlm8(ispace)
                    888: c      call avlm16(ispace)
                    889: c      call crunch
                    890: c
                    891: c
                    892: c general comments:
                    893: c      for each block which is allocated, a 5-word entry is maintained
                    894: c in a table kept in high memory, of the form
                    895: c
                    896: c        word      contents
                    897: c        ----      --------
                    898: c
                    899: c          1       index of imem(.) into origin of block
                    900: c                  i.e. contents of pointer (used for error check)
                    901: c          2       block size (in words)
                    902: c          3       number of words in use
                    903: c          4       address of variable containing block origin
                    904: c          5       number of words used per table entry
                    905: c
                    906: c      all allocated blocks are an 'even' (nxtevn) number of words in length,
                    907: c where a 'word' is the storage unit required for an 'integer' variable.
                    908: c      since block repositioning may be necessary, the convention that
                    909: c only one variable contain a block origin should be observed.
                    910: c      for *getmem*, *ipntr* is set such that *array(ipntr+1)* is the
                    911: c first word of the allocated block.  'ipntr' is set to address the first
                    912: c entry of the table when used with the appropriate variable type, i.e.,
                    913: c nodplc(ipntr+1), value(ipntr+1), or cvalue(ipntr+1).
                    914: c      for *clrmem*, *ipntr* is set to 'invalid' to enable rapid detection
                    915: c of an attempt to use a cleared block.
                    916: c      if any fatal errors are found, a message is printed and a flag
                    917: c set inhibiting further action until *setmem* is called.  (in this
                    918: c context, insufficient memory is considered a fatal error.)
                    919: c      throughout this routine, *ldval* always contains the subscript of
                    920: c the last addressable word of memory, *memavl* always contains the
                    921: c number of available words of memory, *numblk* always contains the
                    922: c number of allocated blocks, and istack(*loctab* +1) always contains
                    923: c the first word of the block table.
                    924: c
                    925:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                    926:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                    927:      2   nwd8,nwd16
                    928: c
                    929: c.. arguments to memory manager are set up as arrays, even though
                    930: c.. the calling programs usually use simple variables for arguments.
                    931: c.. this is necessary if we are to guarantee that the parameters are
                    932: c.. passed by 'address' and not by 'value'.  we must insure that locf(arg)
                    933: c.. returns the address of the argument, and not the address of a local
                    934: c.. copy of the argument.  as currently configured, this subroutine should
                    935: c.. work on any ansi fortran compiler, provided the function 'locf' can
                    936: c.. be provided.
                    937:       dimension ipntr(1)
                    938: c
                    939:       logical memptr
                    940: c
                    941: c...  approximate time required to copy *nwords* integer values
                    942:       nwd4=1
                    943:       nwd8=2
                    944:       nwd16=4
                    945:       memerr=0
                    946:       nevn=nxtevn(1)
                    947:       icheck=mod(nevn,nwd4)+mod(nevn,nwd8)+mod(nevn,nwd16)+
                    948:      1  mod(nxtmem(1),nevn)
                    949:       if(icheck.eq.0) go to 2
                    950:       memerr=1
                    951:       call errmem(6,memerr,ipntr(1))
                    952:     2 cpyknt=0.0d0
                    953:        ifamwa=locf(ipntr(1))
                    954:       maxmem=ksize
                    955:       ntab=nxtevn(5)
                    956: c... add 'lorg' to an address and you get the 'istack' index to that word
                    957:       lorg=1-locf(istack(1))
                    958:       ifwa=ifamwa+lorg-1
                    959:       nwoff=locf(ipntr(1))+lorg-1
                    960:       icore=nxtmem(1)
                    961: c... don't take chances, back off from 'end of memory' by nxtevn(1)
                    962:       ldval=ifwa+nxtmem(1)-nxtevn(1)
                    963:       memavl=ldval-ntab-ifwa
                    964:       maxcor=0
                    965:       maxuse=0
                    966:       call memory
                    967:       if(memerr.ne.0) call errmem(6,memerr,ipntr(1))
                    968:       numblk=1
                    969:       loctab=ldval-ntab
                    970:        istack(loctab+1)=0
                    971:       istack(loctab+2)=memavl
                    972:       istack(loctab+3)=0
                    973:       istack(loctab+4)=-1
                    974:       istack(loctab+5)=1
                    975:       return
                    976:       end
                    977:       subroutine getm4(ipntr,ksize)
                    978:       implicit double precision (a-h,o-z)
                    979:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                    980:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                    981:      2   nwd8,nwd16
                    982:       dimension ipntr(1)
                    983:       iwsize=nwd4
                    984:       call getmx(ipntr(1),ksize,iwsize)
                    985:       return
                    986:       end
                    987:       subroutine getm8(ipntr,ksize)
                    988:       implicit double precision (a-h,o-z)
                    989:       dimension ipntr(1)
                    990:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                    991:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                    992:      2   nwd8,nwd16
                    993:       iwsize=nwd8
                    994:       call getmx(ipntr(1),ksize,iwsize)
                    995:       return
                    996:       end
                    997:       subroutine getm16(ipntr,ksize)
                    998:       implicit double precision (a-h,o-z)
                    999:       dimension ipntr(1)
                   1000:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1001:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1002:      2   nwd8,nwd16
                   1003:       iwsize=nwd16
                   1004:       call getmx(ipntr(1),ksize,iwsize)
                   1005:       return
                   1006:       end
                   1007:       subroutine getmx(ipntr,ksize,iwsize)
                   1008:       implicit double precision (a-h,o-z)
                   1009:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1010:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1011:      2   nwd8,nwd16
                   1012:       logical memptr
                   1013:       dimension ipntr(1)
                   1014: c
                   1015: c***  getmem - get block
                   1016: c
                   1017: c
                   1018:       isize=ksize*iwsize
                   1019: c...  check for valid size
                   1020:       if (isize.ge.0) go to 5
                   1021:       memerr=2
                   1022:       call errmem(3,memerr,ipntr(1))
                   1023: c...  check for attempt to reallocate existing block
                   1024:     5 if (.not.memptr(ipntr(1))) go to 8
                   1025:       memerr=3
                   1026:       call errmem(3,memerr,ipntr(1))
                   1027:     8 jsize=nxtevn(isize)
                   1028:       call comprs(0,ldval)
                   1029: c...  check if enough space already there
                   1030:       need=jsize+ntab-memavl
                   1031:       if (need.le.0) go to 10
                   1032: c...  insufficient space -- bump memory size
                   1033:       need=nxtmem(need)
                   1034:       icore=icore+need
                   1035:       call memory
                   1036:       if(memerr.ne.0) call errmem(3,memerr,ipntr(1))
                   1037:       ltab1=ldval-ntab
                   1038:       istack(ltab1+2)=istack(ltab1+2)+need
                   1039: c...  relocate block entry table
                   1040:       nwords=numblk*ntab
                   1041:       cpyknt=cpyknt+dfloat(nwords)
                   1042:       call copy4(istack(loctab+1),istack(loctab+need+1),nwords)
                   1043:       loctab=loctab+need
                   1044:       ldval=ldval+need
                   1045:       memavl=memavl+need
                   1046: c...  a block large enough now exists -- allocate it
                   1047:    10 ltab1=ldval-ntab
                   1048:       morg=istack(ltab1+1)
                   1049:       msiz=istack(ltab1+2)
                   1050:       muse=istack(ltab1+3)
                   1051:       muse=nxtevn(muse)
                   1052:       madr=istack(ltab1+4)
                   1053: c...  construct new table entry
                   1054:    15 istack(ltab1+2)=muse
                   1055:       loctab=loctab-ntab
                   1056:       nwords=numblk*ntab
                   1057:       cpyknt=cpyknt+dfloat(nwords)
                   1058:       call copy4(istack(loctab+ntab+1),istack(loctab+1),nwords)
                   1059:       numblk=numblk+1
                   1060:       memavl=memavl-ntab
                   1061:       istack(ltab1+1)=morg+muse
                   1062:       istack(ltab1+2)=msiz-muse-ntab
                   1063: c...  set user size into table entry for this block
                   1064:    20 istack(ltab1+3)=isize
                   1065:       istack(ltab1+4)=locf(ipntr(1))
                   1066:       istack(ltab1+5)=iwsize
                   1067:       memavl=memavl-jsize
                   1068:       ipntr(1)=istack(ltab1+1)/iwsize
                   1069:       call memadj
                   1070:       return
                   1071:       end
                   1072:       subroutine avlm4(iavl)
                   1073:       implicit double precision (a-h,o-z)
                   1074:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1075:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1076:      2   nwd8,nwd16
                   1077: c
                   1078: c***  avlmem - how much space is available ?
                   1079: c
                   1080:       iavl=((maxmem-icore)/nxtmem(1))*nxtmem(1)-ntab+memavl
                   1081:       iavl=iavl/nwd4
                   1082:       return
                   1083:       end
                   1084:       subroutine avlm8(iavl)
                   1085:       implicit double precision (a-h,o-z)
                   1086:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1087:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1088:      2   nwd8,nwd16
                   1089:       iavl=((maxmem-icore)/nxtmem(1))*nxtmem(1)-ntab+memavl
                   1090:       iavl=iavl/nwd8
                   1091:       return
                   1092:       end
                   1093:       subroutine avlm16(iavl)
                   1094:       implicit double precision (a-h,o-z)
                   1095:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1096:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1097:      2   nwd8,nwd16
                   1098:       iavl=((maxmem-icore)/nxtmem(1))*nxtmem(1)-ntab+memavl
                   1099:       iavl=iavl/nwd16
                   1100:       return
                   1101:       end
                   1102:       subroutine relmem(ipntr,ksize)
                   1103:       implicit double precision (a-h,o-z)
                   1104:       dimension ipntr(1)
                   1105:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1106:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1107:      2   nwd8,nwd16
                   1108:       logical memptr
                   1109: c
                   1110: c***  relmem - release part of block
                   1111: c
                   1112: c
                   1113: c...  check for valid pointer
                   1114:       if (memptr(ipntr(1))) go to 10
                   1115:       memerr=5
                   1116:       call errmem(5,memerr,ipntr(1))
                   1117:    10 isize=ksize*istack(ltab+5)
                   1118: c...  check for valid size
                   1119:       if (isize.ge.0) go to 20
                   1120:       memerr=2
                   1121:       call errmem(5,memerr,ipntr(1))
                   1122:    20 jsize=istack(ltab+3)
                   1123:       if (isize.le.jsize) go to 30
                   1124:       memerr=6
                   1125:       call errmem(5,memerr,ipntr(1))
                   1126:    30 istack(ltab+3)=istack(ltab+3)-isize
                   1127:       memavl=memavl+(nxtevn(jsize)-nxtevn(istack(ltab+3)))
                   1128:       call memadj
                   1129:       return
                   1130:       end
                   1131:       subroutine extmem(ipntr,ksize)
                   1132:       implicit double precision (a-h,o-z)
                   1133:       dimension ipntr(1)
                   1134:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1135:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1136:      2   nwd8,nwd16
                   1137:       logical memptr
                   1138: c
                   1139: c***  extmem - extend size of existing block
                   1140: c
                   1141: c
                   1142: c...  check for valid pointer
                   1143:       if (memptr(ipntr(1))) go to 10
                   1144:       memerr=5
                   1145:       call errmem(2,memerr,ipntr(1))
                   1146:    10 isize=ksize*istack(ltab+5)
                   1147: c...  check for valid size
                   1148:       if (isize.ge.0) go to 20
                   1149:       memerr=2
                   1150:       call errmem(2,memerr,ipntr(1))
                   1151: c...  check if enough space already there
                   1152:    20 if ((istack(ltab+2)-istack(ltab+3)).ge.isize) go to 40
                   1153:       need=nxtevn(isize)-memavl
                   1154:       if (need.le.0) go to 30
                   1155: c...  insufficient space -- bump memory size
                   1156:       need=nxtmem(need)
                   1157:       icore=icore+need
                   1158:       call memory
                   1159:       if(memerr.ne.0) call errmem(2,memerr,ipntr(1))
                   1160:       ltab1=ldval-ntab
                   1161:       istack(ltab1+2)=istack(ltab1+2)+need
                   1162: c...  relocate block entry table
                   1163:       nwords=numblk*ntab
                   1164:       cpyknt=cpyknt+dfloat(nwords)
                   1165:       call copy4(istack(loctab+1),istack(loctab+need+1),nwords)
                   1166:       loctab=loctab+need
                   1167:       ldval=ldval+need
                   1168:       memavl=memavl+need
                   1169:       ltab=ltab+need
                   1170: c...  move blocks to make space
                   1171:    30 continue
                   1172:       call comprs(0,ltab)
                   1173:       call comprs(1,ltab)
                   1174:    40 jsize=istack(ltab+3)
                   1175:       istack(ltab+3)=istack(ltab+3)+isize
                   1176:       memavl=memavl-(nxtevn(istack(ltab+3))-nxtevn(jsize))
                   1177:       call memadj
                   1178:       return
                   1179:       end
                   1180:       subroutine sizmem(ipntr,ksize)
                   1181:       implicit double precision (a-h,o-z)
                   1182:       dimension ipntr(1)
                   1183:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1184:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1185:      2   nwd8,nwd16
                   1186:       logical memptr
                   1187: c
                   1188: c***  sizmem - determine size of existing block
                   1189: c
                   1190: c
                   1191: c...  check for valid pointer
                   1192:       if (memptr(ipntr(1))) go to 10
                   1193:       memerr=5
                   1194:       call errmem(7,memerr,ipntr(1))
                   1195:    10 ksize=istack(ltab+3)/istack(ltab+5)
                   1196:       return
                   1197:       end
                   1198:       subroutine clrmem(ipntr)
                   1199:       implicit double precision (a-h,o-z)
                   1200:       dimension ipntr(1)
                   1201:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1202:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1203:      2   nwd8,nwd16
                   1204:       logical memptr
                   1205: c
                   1206: c***  clrmem - release block
                   1207: c
                   1208: c
                   1209: c...  check that pointer is valid
                   1210:       if (memptr(ipntr(1))) go to 10
                   1211:       memerr=5
                   1212:       call errmem(1,memerr,ipntr(1))
                   1213:    10 msiz=istack(ltab+2)
                   1214:       muse=istack(ltab+3)
                   1215:       memavl=memavl+nxtevn(muse)
                   1216: c...  assumption:  first allocated block is never cleared.
                   1217:       ltab1=ltab-ntab
                   1218:       istack(ltab1+2)=istack(ltab1+2)+msiz
                   1219: c...  reposition the block table
                   1220:       nwords=ltab-loctab
                   1221:       cpyknt=cpyknt+dfloat(nwords)
                   1222:       call copy4(istack(loctab+1),istack(loctab+ntab+1),nwords)
                   1223:       numblk=numblk-1
                   1224:       loctab=loctab+ntab
                   1225:       memavl=memavl+ntab
                   1226:       ltab1=ldval-ntab
                   1227:       istack(ltab1+2)=istack(ltab1+2)+ntab
                   1228:       ipntr(1)=2**31-1
                   1229:       call memadj
                   1230:       return
                   1231:       end
                   1232:       subroutine ptrmem(ipntr,ipntr2)
                   1233:       implicit double precision (a-h,o-z)
                   1234:       dimension ipntr(1),ipntr2(1)
                   1235:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1236:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1237:      2   nwd8,nwd16
                   1238:       logical memptr
                   1239: c
                   1240: c***  ptrmem - reset memory pointer
                   1241: c
                   1242: c...  verify that pointer is valid
                   1243:       if (memptr(ipntr(1))) go to 10
                   1244:       memerr=5
                   1245:       call errmem(4,memerr,ipntr(1))
                   1246: c...  reset block pointer to be *ipntr2*
                   1247:    10 ipntr2(1)=ipntr(1)
                   1248:       istack(ltab+4)=locf(ipntr2(1))
                   1249:       call memadj
                   1250:       return
                   1251:       end
                   1252:       subroutine crunch
                   1253:       implicit double precision (a-h,o-z)
                   1254:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1255:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1256:      2   nwd8,nwd16
                   1257: c
                   1258: c***  crunch - force memory compaction
                   1259: c
                   1260:       call comprs(0,ldval)
                   1261:       call memadj
                   1262:       return
                   1263:       end
                   1264:       subroutine errmem(inam,ierror,ipntr)
                   1265:       implicit double precision (a-h,o-z)
                   1266:       dimension ipntr(1)
                   1267:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1268:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1269:      2   nwd8,nwd16
                   1270:       dimension errnam(7)
                   1271:       data errnam /6hclrmem,6hextmem,6hgetmem,6hptrmem,6hrelmem,
                   1272:      1   6hsetmem,6hsizmem/
                   1273: c
                   1274:       go to (200,410,420,300,510,530),ierror
                   1275: c
                   1276: c*** error(s) found ***
                   1277: c
                   1278: c.. nxtevn and/or nxtmem incompatible with nwd4, nwd8, and nwd16
                   1279: c
                   1280:   200 write(6,201)
                   1281:   201 format('0memory manager variables nwd4-8-16 incompatible with nxte
                   1282:      1vn and nxtmem')
                   1283:       go to 900
                   1284: c
                   1285: c...  memory needs exceed maximum available space
                   1286:   300 write (6,301) maxmem
                   1287:   301 format('0*error*:  memory needs exceed',i6,/,
                   1288:      1  '0probable remedy, replace your "// exec spice" card with',/
                   1289:      2  '0// exec spice,region=2000k')
                   1290:       go to 900
                   1291: c...    *isize* < 0
                   1292:   410 write(6,411)
                   1293:   411 format('0size parameter negative')
                   1294:       go to 900
                   1295: c...  getmem:  attempt to reallocate existing block
                   1296:   420 write(6,421)
                   1297:   421 format('0attempt to reallocate existing table')
                   1298:       go to 900
                   1299: c...    *ipntr* invalid
                   1300:   510 write(6,511)
                   1301:   511 format('0table pointer invalid')
                   1302:       go to 900
                   1303: c...  relmem:  *isize* larger than indicated block
                   1304:   530 write(6,531)
                   1305:   531 format('0attempt to release more than total table')
                   1306: c...  issue error message
                   1307:   900 write (6,901) errnam(inam)
                   1308:   901 format('0*abort*:  internal memory manager error at entry ',
                   1309:      1  a7)
                   1310:   950 call dmpmem(ipntr(1))
                   1311:  1000 stop
                   1312:       end
                   1313:       subroutine memadj
                   1314:       implicit double precision (a-h,o-z)
                   1315:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1316:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1317:      2   nwd8,nwd16
                   1318: c
                   1319: c*** adjust memory downward ***
                   1320: c
                   1321:    50 maxuse=max0(maxuse,ldval-memavl-ifwa)
                   1322:       memdec=2*nxtmem(1)
                   1323:       if (memavl.lt.memdec) return
                   1324: c...  compress current allocations of memory
                   1325:       call comprs(0,ldval)
                   1326: c...  adjust memory size
                   1327:       memdel=0
                   1328:    60 icore=icore-memdec
                   1329:       memdel=memdel+memdec
                   1330:       memavl=memavl-memdec
                   1331:       if (memavl.ge.memdec) go to 60
                   1332:       ltab1=ldval-ntab
                   1333:       istack(ltab1+2)=istack(ltab1+2)-memdel
                   1334: c...  relocate block entry table
                   1335:       nwords=numblk*ntab
                   1336:       cpyknt=cpyknt+dfloat(nwords)
                   1337:       call copy4(istack(loctab+1),istack(loctab-memdel+1),nwords)
                   1338:       loctab=loctab-memdel
                   1339:       ldval=ldval-memdel
                   1340:       call memory
                   1341:       return
                   1342:       end
                   1343:       integer function nxtevn(n)
                   1344: c
                   1345: c.. function returns the smallest value nxtevn greater than or equal to
                   1346: c.. n which is evenly divisible by 'nwd4, nwd8, and nwd16' as defined
                   1347: c.. in setmem
                   1348: c
                   1349:       nxtevn=((n+3)/4)*4
                   1350:       return
                   1351:       end
                   1352:       integer function nxtmem(memwds)
                   1353: c
                   1354: c.. function returns the in nxtmem the next available memory size
                   1355: c.. (which must be evenly divisible by 'nwd4, nwd8, and nwd16' as
                   1356: c.. defined in setmem
                   1357: c
                   1358:       nxtmem=((memwds+1999)/2000)*2000
                   1359:       return
                   1360:       end
                   1361:       subroutine comprs(icode,limit)
                   1362:       implicit double precision (a-h,o-z)
                   1363: c
                   1364: c      this routine compresses all available memory into a single block.
                   1365: c if *icode* is zero, compression of memory from word 1 to *limit* is
                   1366: c done;  otherwise, compression from *ldval* down to *limit* is done.
                   1367: c
                   1368:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1369:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1370:      2   nwd8,nwd16
                   1371: c
                   1372: c...  approximate time required to copy *nwords* real values
                   1373:       if (icode.ne.0) go to 100
                   1374:       nblk=numblk
                   1375:       ltab2=loctab
                   1376:    10 ltab1=ltab2
                   1377:       if (ltab1.ge.limit) go to 200
                   1378:       if (nblk.eq.1) go to 200
                   1379:       nblk=nblk-1
                   1380:       ltab2=ltab1+ntab
                   1381:       morg=istack(ltab1+1)
                   1382:       msiz=istack(ltab1+2)
                   1383:       muse=istack(ltab1+3)
                   1384:       muse=nxtevn(muse)
                   1385:       if (msiz.eq.muse) go to 10
                   1386: c...  move succeeding block down
                   1387:       morg2=istack(ltab2+1)
                   1388:       muse2=istack(ltab2+3)
                   1389:       madr2=istack(ltab2+4)
                   1390:       iwsize=istack(ltab2+5)
                   1391:       if (madr2.ne.0) go to 15
                   1392:       if (muse2.eq.0) go to 20
                   1393:    15 cpyknt=cpyknt+dfloat(muse2)
                   1394:       call copy4(istack(nwoff+morg2+1),istack(nwoff+morg+muse+1),muse2)
                   1395:       istack(lorg+madr2)=(morg+muse)/iwsize
                   1396:    20 istack(ltab1+2)=muse
                   1397:       istack(ltab2+1)=morg+muse
                   1398:       istack(ltab2+2)=istack(ltab2+2)+(msiz-muse)
                   1399:       go to 10
                   1400: c
                   1401: c
                   1402:   100 nblk=numblk
                   1403:       ltab2=ldval-ntab
                   1404:   110 ltab1=ltab2
                   1405:       if (ltab1.le.limit) go to 200
                   1406:       if (nblk.eq.1) go to 200
                   1407:       nblk=nblk-1
                   1408:       ltab2=ltab1-ntab
                   1409:       morg=istack(ltab1+1)
                   1410:       msiz=istack(ltab1+2)
                   1411:       muse=istack(ltab1+3)
                   1412:       muse=nxtevn(muse)
                   1413:       madr=istack(ltab1+4)
                   1414:       iwsize=istack(ltab1+5)
                   1415:       mspc=msiz-muse
                   1416:       if (mspc.eq.0) go to 110
                   1417:       cpyknt=cpyknt+dfloat(muse)
                   1418:       call copy4(istack(nwoff+morg+1),istack(nwoff+morg+mspc+1),muse)
                   1419:       istack(ltab1+1)=morg+mspc
                   1420:       istack(ltab1+2)=muse
                   1421:       istack(ltab2+2)=istack(ltab2+2)+mspc
                   1422:       if (madr.eq.0) go to 110
                   1423:       istack(lorg+madr)=(morg+mspc)/iwsize
                   1424:       go to 110
                   1425: c...  all done
                   1426:   200 return
                   1427:       end
                   1428:       logical function memptr(ipntr)
                   1429:       implicit double precision (a-h,o-z)
                   1430: c
                   1431: c      this routine checks whether *ipntr* is a valid block pointer.
                   1432: c if it is valid, *ltab* is set to point to the corresponding entry in
                   1433: c the block table.
                   1434: c
                   1435: c... ipntr is an array to avoid 'call by value' problems (see setmem)
                   1436:       dimension ipntr(1)
                   1437:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1438:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1439:      2   nwd8,nwd16
                   1440: c
                   1441:       memptr=.false.
                   1442:       ltab=loctab
                   1443:       locpnt=locf(ipntr(1))
                   1444:       do 20 i=1,numblk
                   1445:       if (locpnt.ne.istack(ltab+4)) go to 10
                   1446:       if (ipntr(1)*istack(ltab+5).ne.istack(ltab+1)) go to 10
                   1447:       memptr=.true.
                   1448:       go to 30
                   1449:    10 ltab=ltab+ntab
                   1450:    20 continue
                   1451:    30 return
                   1452:       end
                   1453:       subroutine dmpmem(ipntr)
                   1454:       implicit double precision (a-h,o-z)
                   1455: c
                   1456: c      this routine prints out the current memory allocation map.
                   1457: c *ipntr* is the table pointer of the current memory manager call
                   1458: c
                   1459:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
                   1460:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
                   1461:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
                   1462:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
                   1463:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
                   1464:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
                   1465:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1466:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1467:      2   nwd8,nwd16
                   1468: c... ipntr is an array to avoid 'call by value' problems
                   1469:       dimension ipntr(1)
                   1470:       dimension aptr(61)
                   1471:       data aptr /6hielmnt,6hisbckt,6hnsbckt,6hiunsat,6hnunsat,6hitemps,
                   1472:      1 6hnumtem,6hisens ,6hnsens ,6hifour ,6hnfour ,6hifield,
                   1473:      2 6hicode ,6hidelim,6hicolum,6hinsize,
                   1474:      3 6hjunode,6hlsbkpt,6hnumbkp,6hiorder,6hjmnode,
                   1475:      4 6hiur   ,6hiuc   ,6hilc   ,6hilr   ,6hnumoff,6hisr   ,
                   1476:      5 6hnmoffc,6hiseq  ,6hiseq1  ,6hneqn  ,6hnodevs,
                   1477:      6 6hndiag ,6hiswap ,6hiequa ,6hmacins,6hlvnim1,
                   1478:      7 6hlx0   ,6hlvn   ,6hlynl  ,6hlyu   ,6hlyl   ,
                   1479:      8 6hlx1   ,6hlx2   ,6hlx3   ,6hlx4   ,6hlx5   ,6hlx6   ,
                   1480:      9 6hlx7   ,6hld0   ,6hld1   ,6hltd   ,6himynl ,6himvn  ,6hloutpt,
                   1481:      * 6hnsnod ,6hnsmat ,6hnsval ,6hicnod ,6hicmat ,6hicval /
                   1482:       data ablnk /1h /
                   1483:       iaddr=locf(ielmnt)-1
                   1484:       itemp=locf(ipntr(1))-iaddr
                   1485:       anam=ablnk
                   1486:       if(itemp.gt.0.and.itemp.le.61) anam=aptr(itemp)
                   1487:       iadr=locf(ipntr(1))
                   1488:       write (6,5) anam,iadr,icore,maxmem,memavl,ldval
                   1489:     5 format('0current pointer 'a6,'@ = z',z6,/' corsiz=',i7,
                   1490:      1  /' maxmem=',i7,/' avlspc=',i7,/' ldval=',i7,
                   1491:      2  /1h0,24x,'memory allocation map'/14x,'blknum memorg memsiz',
                   1492:      3  '  memuse usrptr  addr    name')
                   1493:       ltab1=loctab
                   1494:       do 20 i=1,numblk
                   1495:       morg=istack(ltab1+1)
                   1496:       msiz=istack(ltab1+2)
                   1497:       muse=istack(ltab1+3)
                   1498:       madr=istack(ltab1+4)
                   1499:       anam=ablnk
                   1500:       ndex=madr-iaddr
                   1501:       if(ndex.gt.0.and.ndex.le.61) anam=aptr(ndex)
                   1502:       jptr=0
                   1503:       if (madr.gt.0) jptr=istack(lorg+madr)
                   1504:       write (6,11) i,morg,msiz,muse,jptr,madr,anam
                   1505:    11 format(13x,5i7,3x,z7,'z',1x,a6)
                   1506:       ltab1=ltab1+ntab
                   1507:    20 continue
                   1508:       write (6,21)
                   1509:    21 format(1h0,24x,'end of allocation map'/)
                   1510:       return
                   1511:       end
                   1512:       subroutine memory
                   1513:       implicit double precision (a-h,o-z)
                   1514:       common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
                   1515:      1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
                   1516:      2   nwd8,nwd16
                   1517:       if(icore.le.maxmem) go to 10
                   1518:       memerr=4
                   1519:       return
                   1520:    10 continue
                   1521:       return
                   1522:       end
                   1523:       subroutine magphs(cvar,xmag,xphs)
                   1524:       implicit double precision (a-h,o-z)
                   1525: c
                   1526: c     this routine computes the magnitude and phase of its complex arg-
                   1527: c ument cvar, storing the results in xmag and xphs.
                   1528: c
                   1529:       common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
                   1530:      1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox
                   1531:       complex*16 cvar
                   1532: c
                   1533: c
                   1534:       xreal=dreal(cvar)
                   1535:       ximag=dimag(cvar)
                   1536:       xmag=dsqrt(xreal*xreal+ximag*ximag)
                   1537:       if (xmag.ge.1.0d-20) go to 10
                   1538:       xmag=1.0d-20
                   1539:       xphs=0.0d0
                   1540:       return
                   1541:    10 xphs=rad*datan2(ximag,xreal)
                   1542:       return
                   1543:       end
                   1544:       integer function xxor(a,b)
                   1545:       implicit double precision (a-h,o-z)
                   1546: c
                   1547: c     this routine computes a single-precision integer result which is
                   1548: c the result of exclusive-or*ing the two real-valued arguments a and b
                   1549: c together.
                   1550: c
                   1551:       xxor=1
                   1552:       if(a.eq.b) xxor=0
                   1553:       return
                   1554:       end
                   1555:       subroutine outnam(loc,ktype,string,ipos)
                   1556:       implicit double precision (a-h,o-z)
                   1557: c
                   1558: c     this routine constructs the 'name' for the output variable indi-
                   1559: c cated by loc, adding the characters to the character array 'string',
                   1560: c beginning with the position marked by ipos.
                   1561: c
                   1562:       common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
                   1563:      1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
                   1564:      2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
                   1565:      3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
                   1566:      4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
                   1567:      5   imynl,imvn,lcvn,loutpt,nsnod,nsmat,nsval,icnod,icmat,icval
                   1568:       common /blank/ value(50000)
                   1569:       integer nodplc(64)
                   1570:       complex*16 cvalue(32)
                   1571:       equivalence (value(1),nodplc(1),cvalue(1))
                   1572: c
                   1573:       dimension string(1)
                   1574:       dimension aout(19),lenout(19),aopt(5),lenopt(5)
                   1575:       data aout / 6hv     , 6hvm    , 6hvr    , 6hvi    , 6hvp    ,
                   1576:      1            6hvdb   , 6hi     , 6him    , 6hir    , 6hii    ,
                   1577:      2            6hip    , 6hidb   , 6honoise, 6hinoise, 6hhd2   ,
                   1578:      1            6hhd3   , 6hdim2  , 6hsim2  , 6hdim3   /
                   1579:       data lenout / 1,2,2,2,2,3,1,2,2,2,2,3,6,6,3,3,4,4,4 /
                   1580:       data aopt / 5hmag  , 5hreal , 5himag , 5hphase, 5hdb    /
                   1581:       data lenopt / 3,4,4,5,2 /
                   1582:       data alprn, acomma, arprn, ablnk / 1h(, 1h,, 1h), 1h  /
                   1583: c
                   1584: c
                   1585:       ioutyp=nodplc(loc+5)
                   1586:       if (ioutyp.ge.2) go to 10
                   1587:       lout=ktype+ioutyp*6
                   1588:       go to 20
                   1589:    10 lout=ioutyp+11
                   1590:    20 call move(string,ipos,aout(lout),1,lenout(lout))
                   1591:       ipos=ipos+lenout(lout)
                   1592:       if (ioutyp.ge.2) go to 200
                   1593:       call move(string,ipos,alprn,1,1)
                   1594:       ipos=ipos+1
                   1595:       if (ioutyp.ne.0) go to 100
                   1596:       node1=nodplc(loc+2)
                   1597:       call alfnum(nodplc(junode+node1),string,ipos)
                   1598:       node2=nodplc(loc+3)
                   1599:       if (node2.eq.1) go to 30
                   1600:       call move(string,ipos,acomma,1,1)
                   1601:       ipos=ipos+1
                   1602:       call alfnum(nodplc(junode+node2),string,ipos)
                   1603:    30 call move(string,ipos,arprn,1,1)
                   1604:       ipos=ipos+1
                   1605:       go to 1000
                   1606: c
                   1607:   100 locv=nodplc(loc+1)
                   1608:       anam=value(locv)
                   1609:       achar=ablnk
                   1610:       do 110 i=1,8
                   1611:       call move(achar,1,anam,i,1)
                   1612:       if (achar.eq.ablnk) go to 120
                   1613:       call move(string,ipos,achar,1,1)
                   1614:       ipos=ipos+1
                   1615:   110 continue
                   1616:   120 call move(string,ipos,arprn,1,1)
                   1617:       ipos=ipos+1
                   1618:       go to 1000
                   1619: c
                   1620:   200 if (ktype.eq.1) go to 1000
                   1621:       call move(string,ipos,alprn,1,1)
                   1622:       ipos=ipos+1
                   1623:       call move(string,ipos,aopt(ktype-1),1,lenopt(ktype-1))
                   1624:       ipos=ipos+lenopt(ktype-1)
                   1625:       call move(string,ipos,arprn,1,1)
                   1626:       ipos=ipos+1
                   1627: c
                   1628: c  finished
                   1629: c
                   1630:  1000 return
                   1631:       end
                   1632:       subroutine alfnum(number,string,ipos)
                   1633:       implicit double precision (a-h,o-z)
                   1634: c
                   1635: c     this routine converts number into character form, storing the
                   1636: c characters in the character array string, beginning with the position
                   1637: c indicated by ipos.
                   1638: c
                   1639: c **** note that the 'ipos' variable is changed to indicate the position
                   1640: c      of the next unwritten character.  this could clobber constants if
                   1641: c      ipos is not a variable in the calling program
                   1642: c
                   1643:       dimension string(1)
                   1644:       dimension adigit(10)
                   1645:       data adigit / 1h0,1h1,1h2,1h3,1h4,1h5,1h6,1h7,1h8,1h9 /
                   1646:       data aminus / 1h- /
                   1647: c
                   1648: c
                   1649:       num=number
                   1650: c
                   1651: c  check for number < 0
                   1652: c
                   1653:       if (num.ge.0) go to 10
                   1654:       num=-num
                   1655: c...  negative number:  insert minus sign
                   1656:       call move(string,ipos,aminus,1,1)
                   1657:       ipos=ipos+1
                   1658: c
                   1659: c  convert number one digit at a time, in reverse order
                   1660: c
                   1661:    10 istart=ipos
                   1662:    20 numtmp=num/10
                   1663:       idigit=num-numtmp*10
                   1664:       call move(string,ipos,adigit(idigit+1),1,1)
                   1665:       ipos=ipos+1
                   1666:       num=numtmp
                   1667:       if (num.ne.0) go to 20
                   1668:       istop=ipos-1
                   1669: c
                   1670: c  now reverse the order of the digits
                   1671: c
                   1672:    30 if (istop.le.istart) go to 40
                   1673:       call move(tmpdgt,1,string,istart,1)
                   1674:       call move(string,istart,string,istop,1)
                   1675:       call move(string,istop,tmpdgt,1,1)
                   1676:       istart=istart+1
                   1677:       istop=istop-1
                   1678:       go to 30
                   1679: c
                   1680: c  conversion complete
                   1681: c
                   1682:    40 return
                   1683:       end
                   1684:       subroutine getcje
                   1685:       implicit double precision (a-h,o-z)
                   1686:       common /cje/ maxtim,itime,icost
                   1687:       call second(xtime)
                   1688:       itime=xtime
                   1689:       icost=xtime*38.3333
                   1690:       return
                   1691:       end
                   1692:       subroutine move(a,iposa,b,iposb,nchar)
                   1693:       character a(1),b(1)
                   1694:       do 10 i=1,nchar
                   1695:       a(iposa+i-1)=b(iposb+i-1)
                   1696:    10 continue
                   1697:       return
                   1698:       end
                   1699:       subroutine copy4(ifrom,ito,nwords)
                   1700:       implicit double precision (a-h,o-z)
                   1701: c
                   1702:       dimension ifrom(1),ito(1)
                   1703: c     this routine copies a block of #nwords# words (of the appropriate
                   1704: c type) from the array #from# to the array #to#.  it determines from
                   1705: c which end of the block to transfer first, to prevent over-stores which
                   1706: c might over-write the data.
                   1707: c
                   1708:       if (nwords.eq.0) return
                   1709:       if (locf(ifrom(1)).lt.locf(ito(1))) go to 20
                   1710: c...  locf() returns as its value the address of its argument
                   1711:       do 10 i=1,nwords
                   1712:       ito(i)=ifrom(i)
                   1713:    10 continue
                   1714:       return
                   1715: c
                   1716:    20 i=nwords
                   1717:    30 ito(i)=ifrom(i)
                   1718:       i=i-1
                   1719:       if (i.ne.0) go to 30
                   1720:       return
                   1721: c
                   1722: c
                   1723:       end
                   1724:       subroutine copy8(rfrom,rto,nwords)
                   1725:       implicit double precision (a-h,o-z)
                   1726: c
                   1727:       dimension rfrom(1),rto(1)
                   1728:       if (nwords.eq.0) return
                   1729:       if (locf(rfrom(1)).lt.locf(rto(1))) go to 120
                   1730:       do 110 i=1,nwords
                   1731:       rto(i)=rfrom(i)
                   1732:   110 continue
                   1733:       return
                   1734: c
                   1735:   120 i=nwords
                   1736:   130 rto(i)=rfrom(i)
                   1737:       i=i-1
                   1738:       if (i.ne.0) go to 130
                   1739:       return
                   1740: c
                   1741: c
                   1742:       end
                   1743:       subroutine copy16(cfrom,cto,nwords)
                   1744:       implicit double precision (a-h,o-z)
                   1745: c
                   1746:       complex*16 cfrom(1),cto(1)
                   1747:       if (nwords.eq.0) return
                   1748:       if (locf(cfrom(1)).lt.locf(cto(1))) go to 220
                   1749:       do 210 i=1,nwords
                   1750:       cto(i)=cfrom(i)
                   1751:   210 continue
                   1752:       return
                   1753: c
                   1754:   220 i=nwords
                   1755:   230 cto(i)=cfrom(i)
                   1756:       i=i-1
                   1757:       if (i.ne.0) go to 230
                   1758:       return
                   1759:       end
                   1760:       subroutine zero4(iarray,length)
                   1761:       implicit double precision (a-h,o-z)
                   1762: c
                   1763:       dimension iarray(1)
                   1764: c     this routine zeroes the memory locations indicated by array(1)
                   1765: c through array(length).
                   1766: c
                   1767:       if (length.eq.0) return
                   1768:       do 10 i=1,length
                   1769:       iarray(i)=0
                   1770:    10 continue
                   1771:       return
                   1772:       end
                   1773:       subroutine zero8(array,length)
                   1774:       implicit double precision (a-h,o-z)
                   1775: c
                   1776:       dimension array(1)
                   1777: c     this routine zeroes the memory locations indicated by array(1)
                   1778: c through array(length).
                   1779: c
                   1780:       if (length.eq.0) return
                   1781:       do 10 i=1,length
                   1782:       array(i)=0.0d0
                   1783:    10 continue
                   1784:       return
                   1785:       end
                   1786:       subroutine zero16(carray,length)
                   1787:       implicit double precision (a-h,o-z)
                   1788:       complex*16 carray(1)
                   1789: c
                   1790: c     this routine zeroes the memory locations indicated by array(1)
                   1791: c through array(length).
                   1792: c
                   1793:       if (length.eq.0) return
                   1794:       do 10 i=1,length
                   1795:       carray(i)=dcmplx(0.0d0,0.0d0)
                   1796:    10 continue
                   1797:       return
                   1798: c
                   1799: c
                   1800: c
                   1801:       end
                   1802:        integer function locf(ivar)
                   1803:        iabsa=loc(ivar)
                   1804:       locf=iabsa/4
                   1805:        if(iabsa.eq.locf*4) return
                   1806:       write(6,100) iabsa
                   1807:   100 format('0*error*:  system 370 error..address ',t10,
                   1808:      1  ' is not on a 4-byte boundary')
                   1809:       stop
                   1810:       end
                   1811:       subroutine mdate(anam)
                   1812:       implicit double precision (a-h,o-z)
                   1813:       call date(anam)
                   1814:       return
                   1815:       end
                   1816:       subroutine mclock(anam)
                   1817:       implicit double precision (a-h,o-z)
                   1818:       call todalf(anam)
                   1819:   100 return
                   1820:       end
                   1821:       subroutine second(t1)
                   1822:       implicit double precision (a-h,o-z)
                   1823:       dimension ibuff(4)
                   1824:       real*8 t1
                   1825:       call times (ibuff)
                   1826:       t1 = dfloat (ibuff(1)) / 60.d0
                   1827:       return
                   1828:       end
                   1829:       subroutine todalf(anam)
                   1830:       double precision anam
                   1831:       anam=0.0d0
                   1832:       return
                   1833:       end
                   1834:       double precision function cpusec(time)
                   1835:       cpusec=0.0d0
                   1836:       return
                   1837:       end
                   1838:       subroutine date(anam)
                   1839:       double precision anam
                   1840:       anam=1.0d0
                   1841:       return
                   1842:       end
                   1843:       FUNCTION DREAL ( X )
                   1844:       COMPLEX*16 X
                   1845:       DREAL = REAL (X)
                   1846:       RETURN
                   1847:       END
                   1848: cunix FUNCTION DIMAG ( X )
                   1849: cunix COMPLEX*16 X
                   1850: cunix DIMAG = AIMAG (X)
                   1851: cunix RETURN
                   1852: cunix END
                   1853: cunix COMPLEX FUNCTION DCMPLX ( X , Y )
                   1854: cunix DCMPLX = CMPLX ( X , Y )
                   1855: cunix RETURN
                   1856: cunix END
                   1857: cunix COMPLEX FUNCTION DCONJG ( X )
                   1858: cunix COMPLEX*16 X
                   1859: cunix DCONJG = CONJG ( X )
                   1860: cunix RETURN
                   1861: cunix END
                   1862:       FUNCTION IFIXD ( X )
                   1863:       IFIXD = IFIX ( X )
                   1864:       RETURN
                   1865:       END

unix.superglobalmegacorp.com

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