      subroutine magphs(cvar,xmag,xphs)
      implicit double precision (a-h,o-z)
c
c     this routine computes the magnitude and phase of its complex arg-
c ument cvar, storing the results in xmag and xphs.
c
c spice version 2g.6  sccsid=knstnt 3/15/83
      common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
     1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox,
     2   pivtol,pivrel
      complex cvar
c
c
      xreal=dble(real(cvar))
      ximag=dble(aimag(cvar))
      xmag=dsqrt(xreal*xreal+ximag*ximag)
      if (xmag.ge.1.0d-20) go to 10
      xmag=1.0d-20
      xphs=0.0d0
      return
   10 xphs=rad*datan2(ximag,xreal)
      return
      end
      subroutine matloc
      implicit double precision (a-h,o-z)
c
c     this routine stores the locations of the various matrix terms to
c which the different circuit elements contribute.
c
c spice version 2g.6  sccsid=tabinf 3/15/83
      common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
     1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
     2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
     3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
     4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
     5   imynl,imvn,lcvn,nsnod,nsmat,nsval,icnod,icmat,icval,
     6   loutpt,lpol,lzer,irswpf,irswpr,icswpf,icswpr,irpt,jcpt,
     7   irowno,jcolno,nttbr,nttar,lvntmp
c spice version 2g.6  sccsid=cirdat 3/15/83
      common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
     1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs,numalt,numcyc
c spice version 2g.6  sccsid=status 3/15/83
      common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
     1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
     2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
c spice version 2g.6  sccsid=ac 3/15/83
      common /ac/ fstart,fstop,fincr,skw2,refprl,spw2,jacflg,idfreq,
     1   inoise,nosprt,nosout,nosin,idist,idprt
c spice version 2g.6  sccsid=blank 3/15/83
      common /blank/ value(50000)
      integer nodplc(64)
      complex cvalue(32)
      equivalence (value(1),nodplc(1),cvalue(1))
c
c  resistors
c
      loc=locate(1)
  690 if ((loc.eq.0).or.(nodplc(loc+8).ne.0)) go to 700
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      nodplc(loc+4)=indxx(node1,node2)
      nodplc(loc+5)=indxx(node2,node1)
      nodplc(loc+6)=indxx(node1,node1)
      nodplc(loc+7)=indxx(node2,node2)
      loc=nodplc(loc)
      go to 690
c
c  capacitors
c
  700 loc=locate(2)
  710 if ((loc.eq.0).or.(nodplc(loc+12).ne.0)) go to 720
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      nodplc(loc+5)=indxx(node1,node2)
      nodplc(loc+6)=indxx(node2,node1)
      nodplc(loc+10)=indxx(node1,node1)
      nodplc(loc+11)=indxx(node2,node2)
      loc=nodplc(loc)
      go to 710
c
c  inductors
c
  720 loc=locate(3)
  730 if ((loc.eq.0).or.(nodplc(loc+14).ne.0)) go to 740
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ibr=nodplc(loc+5)
      nodplc(loc+6)=indxx(node1,ibr)
      nodplc(loc+7)=indxx(node2,ibr)
      nodplc(loc+8)=indxx(ibr,node1)
      nodplc(loc+9)=indxx(ibr,node2)
      nodplc(loc+13)=indxx(ibr,ibr)
      loc=nodplc(loc)
      go to 730
c
c  mutual inductances
c
  740 loc=locate(4)
  750 if ((loc.eq.0).or.(nodplc(loc+6).ne.0)) go to 760
      nl1=nodplc(loc+2)
      nl2=nodplc(loc+3)
      ibr1=nodplc(nl1+5)
      ibr2=nodplc(nl2+5)
      nodplc(loc+4)=indxx(ibr1,ibr2)
      nodplc(loc+5)=indxx(ibr2,ibr1)
      loc=nodplc(loc)
      go to 750
c
c  nonlinear voltage controlled current sources
c
  760 loc=locate(5)
  762 if ((loc.eq.0).or.(nodplc(loc+13).ne.0)) go to 764
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ndim=nodplc(loc+4)
      lnod=nodplc(loc+6)
      lmat=nodplc(loc+7)
      do 763 i=1,ndim
      node3=nodplc(lnod+1)
      node4=nodplc(lnod+2)
      lnod=lnod+2
      nodplc(lmat+1)=indxx(node1,node3)
      nodplc(lmat+2)=indxx(node1,node4)
      nodplc(lmat+3)=indxx(node2,node3)
      nodplc(lmat+4)=indxx(node2,node4)
      lmat=lmat+4
  763 continue
      loc=nodplc(loc)
      go to 762
c
c  nonlinear voltage controlled voltage sources
c
  764 loc=locate(6)
  766 if ((loc.eq.0).or.(nodplc(loc+14).ne.0)) go to 768
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ndim=nodplc(loc+4)
      ibr=nodplc(loc+6)
      lnod=nodplc(loc+7)
      lmat=nodplc(loc+8)
      nodplc(lmat+1)=indxx(node1,ibr)
      nodplc(lmat+2)=indxx(node2,ibr)
      nodplc(lmat+3)=indxx(ibr,node1)
      nodplc(lmat+4)=indxx(ibr,node2)
      lmat=lmat+4
      do 767 i=1,ndim
      node3=nodplc(lnod+1)
      node4=nodplc(lnod+2)
      lnod=lnod+2
      nodplc(lmat+1)=indxx(ibr,node3)
      nodplc(lmat+2)=indxx(ibr,node4)
      lmat=lmat+2
  767 continue
      loc=nodplc(loc)
      go to 766
c
c  nonlinear current controlled current sources
c
  768 loc=locate(7)
  770 if ((loc.eq.0).or.(nodplc(loc+13).ne.0)) go to 772
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ndim=nodplc(loc+4)
      locvs=nodplc(loc+6)
      lmat=nodplc(loc+7)
      do 771 i=1,ndim
      locvst=nodplc(locvs+i)
      ibr=nodplc(locvst+6)
      nodplc(lmat+1)=indxx(node1,ibr)
      nodplc(lmat+2)=indxx(node2,ibr)
      lmat=lmat+2
  771 continue
      loc=nodplc(loc)
      go to 770
c
c  nonlinear current controlled voltage sources
c
  772 loc=locate(8)
  774 if ((loc.eq.0).or.(nodplc(loc+14).ne.0)) go to 780
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ndim=nodplc(loc+4)
      ibr=nodplc(loc+6)
      locvs=nodplc(loc+7)
      lmat=nodplc(loc+8)
      nodplc(lmat+1)=indxx(node1,ibr)
      nodplc(lmat+2)=indxx(node2,ibr)
      nodplc(lmat+3)=indxx(ibr,node1)
      nodplc(lmat+4)=indxx(ibr,node2)
      lmat=lmat+4
      do 775 i=1,ndim
      locvst=nodplc(locvs+i)
      kbr=nodplc(locvst+6)
      nodplc(lmat+i)=indxx(ibr,kbr)
  775 continue
      loc=nodplc(loc)
      go to 774
c
c  voltage sources
c
  780 loc=locate(9)
  790 if ((loc.eq.0).or.(nodplc(loc+11).ne.0)) go to 800
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      iptr=nodplc(loc+6)
      nodplc(loc+7)=indxx(node1,iptr)
      nodplc(loc+8)=indxx(node2,iptr)
      nodplc(loc+9)=indxx(iptr,node1)
      nodplc(loc+10)=indxx(iptr,node2)
      loc=nodplc(loc)
      go to 790
c
c  diodes
c
  800 loc=locate(11)
  810 if ((loc.eq.0).or.(nodplc(loc+16).ne.0)) go to 820
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      nodplc(loc+7)=indxx(node1,node3)
      nodplc(loc+8)=indxx(node2,node3)
      nodplc(loc+9)=indxx(node3,node1)
      nodplc(loc+10)=indxx(node3,node2)
      nodplc(loc+13)=indxx(node1,node1)
      nodplc(loc+14)=indxx(node2,node2)
      nodplc(loc+15)=indxx(node3,node3)
      loc=nodplc(loc)
      go to 810
c
c  transistors
c
  820 loc=locate(12)
  830 if ((loc.eq.0).or.(nodplc(loc+36).ne.0)) go to 840
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      node4=nodplc(loc+5)
      node5=nodplc(loc+6)
      node6=nodplc(loc+7)
      node7=nodplc(loc+30)
      nodplc(loc+10)=indxx(node1,node4)
      nodplc(loc+11)=indxx(node2,node5)
      nodplc(loc+12)=indxx(node3,node6)
      nodplc(loc+13)=indxx(node4,node1)
      nodplc(loc+14)=indxx(node4,node5)
      nodplc(loc+15)=indxx(node4,node6)
      nodplc(loc+16)=indxx(node5,node2)
      nodplc(loc+17)=indxx(node5,node4)
      nodplc(loc+18)=indxx(node5,node6)
      nodplc(loc+19)=indxx(node6,node3)
      nodplc(loc+20)=indxx(node6,node4)
      nodplc(loc+21)=indxx(node6,node5)
      nodplc(loc+24)=indxx(node1,node1)
      nodplc(loc+25)=indxx(node2,node2)
      nodplc(loc+26)=indxx(node3,node3)
      nodplc(loc+27)=indxx(node4,node4)
      nodplc(loc+28)=indxx(node5,node5)
      nodplc(loc+29)=indxx(node6,node6)
      nodplc(loc+31)=indxx(node7,node7)
      nodplc(loc+32)=indxx(node4,node7)
      nodplc(loc+33)=indxx(node7,node4)
      nodplc(loc+34)=indxx(node2,node4)
      nodplc(loc+35)=indxx(node4,node2)
      loc=nodplc(loc)
      go to 830
c
c  jfets
c
  840 loc=locate(13)
  850 if ((loc.eq.0).or.(nodplc(loc+25).ne.0)) go to 860
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      node4=nodplc(loc+5)
      node5=nodplc(loc+6)
      nodplc(loc+9)=indxx(node1,node4)
      nodplc(loc+10)=indxx(node2,node4)
      nodplc(loc+11)=indxx(node2,node5)
      nodplc(loc+12)=indxx(node3,node5)
      nodplc(loc+13)=indxx(node4,node1)
      nodplc(loc+14)=indxx(node4,node2)
      nodplc(loc+15)=indxx(node4,node5)
      nodplc(loc+16)=indxx(node5,node2)
      nodplc(loc+17)=indxx(node5,node3)
      nodplc(loc+18)=indxx(node5,node4)
      nodplc(loc+20)=indxx(node1,node1)
      nodplc(loc+21)=indxx(node2,node2)
      nodplc(loc+22)=indxx(node3,node3)
      nodplc(loc+23)=indxx(node4,node4)
      nodplc(loc+24)=indxx(node5,node5)
      loc=nodplc(loc)
      go to 850
c
c  mosfets
c
  860 loc=locate(14)
  870 if ((loc.eq.0).or.(nodplc(loc+33).ne.0)) go to 900
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      node4=nodplc(loc+5)
      node5=nodplc(loc+6)
      node6=nodplc(loc+7)
      nodplc(loc+10)=indxx(node1,node5)
      nodplc(loc+11)=indxx(node2,node4)
      nodplc(loc+12)=indxx(node2,node5)
      nodplc(loc+13)=indxx(node2,node6)
      nodplc(loc+14)=indxx(node3,node6)
      nodplc(loc+15)=indxx(node4,node2)
      nodplc(loc+16)=indxx(node4,node5)
      nodplc(loc+17)=indxx(node4,node6)
      nodplc(loc+18)=indxx(node5,node1)
      nodplc(loc+19)=indxx(node5,node2)
      nodplc(loc+20)=indxx(node5,node4)
      nodplc(loc+21)=indxx(node5,node6)
      nodplc(loc+22)=indxx(node6,node2)
      nodplc(loc+23)=indxx(node6,node3)
      nodplc(loc+24)=indxx(node6,node4)
      nodplc(loc+25)=indxx(node6,node5)
      nodplc(loc+27)=indxx(node1,node1)
      nodplc(loc+28)=indxx(node2,node2)
      nodplc(loc+29)=indxx(node3,node3)
      nodplc(loc+30)=indxx(node4,node4)
      nodplc(loc+31)=indxx(node5,node5)
      nodplc(loc+32)=indxx(node6,node6)
      loc=nodplc(loc)
      go to 870
c
c  transmission lines
c
  900 loc=locate(17)
  910 if ((loc.eq.0).or.(nodplc(loc+33).ne.0)) go to 1000
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      node4=nodplc(loc+5)
      ni1=nodplc(loc+6)
      ni2=nodplc(loc+7)
      ibr1=nodplc(loc+8)
      ibr2=nodplc(loc+9)
      nodplc(loc+10)=indxx(node1,node1)
      nodplc(loc+11)=indxx(node1,ni1)
      nodplc(loc+12)=indxx(node2,ibr1)
      nodplc(loc+13)=indxx(node3,node3)
      nodplc(loc+14)=indxx(node4,ibr2)
      nodplc(loc+15)=indxx(ni1,node1)
      nodplc(loc+16)=indxx(ni1,ni1)
      nodplc(loc+17)=indxx(ni1,ibr1)
      nodplc(loc+18)=indxx(ni2,ni2)
      nodplc(loc+19)=indxx(ni2,ibr2)
      nodplc(loc+20)=indxx(ibr1,node2)
      nodplc(loc+21)=indxx(ibr1,node3)
      nodplc(loc+22)=indxx(ibr1,node4)
      nodplc(loc+23)=indxx(ibr1,ni1)
      nodplc(loc+24)=indxx(ibr1,ibr2)
      nodplc(loc+25)=indxx(ibr2,node1)
      nodplc(loc+26)=indxx(ibr2,node2)
      nodplc(loc+27)=indxx(ibr2,node4)
      nodplc(loc+28)=indxx(ibr2,ni2)
      nodplc(loc+29)=indxx(ibr2,ibr1)
      nodplc(loc+31)=indxx(node3,ni2)
      nodplc(loc+32)=indxx(ni2,node3)
      loc=nodplc(loc)
      go to 910
c
c  finished
c
 1000 return
      end
      subroutine matptr
      implicit double precision (a-h,o-z)
c
c     this routine (by calls to the routine reserve) establishes the
c nonzero-element structure of the circuit equation coefficient matrix.
c
c spice version 2g.6  sccsid=tabinf 3/15/83
      common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
     1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
     2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
     3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
     4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
     5   imynl,imvn,lcvn,nsnod,nsmat,nsval,icnod,icmat,icval,
     6   loutpt,lpol,lzer,irswpf,irswpr,icswpf,icswpr,irpt,jcpt,
     7   irowno,jcolno,nttbr,nttar,lvntmp
c spice version 2g.6  sccsid=cirdat 3/15/83
      common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
     1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs,numalt,numcyc
c spice version 2g.6  sccsid=blank 3/15/83
      common /blank/ value(50000)
      integer nodplc(64)
      complex cvalue(32)
      equivalence (value(1),nodplc(1),cvalue(1))
c
c  allocate and initialize storage
c
      call getm4(isr,nstop+1)
      numvs=jelcnt(3)+jelcnt(6)+jelcnt(8)+jelcnt(9)+2*jelcnt(17)
      call getm4(iseq,numvs)
      call getm4(iseq1,numvs)
      call getm4(neqn,numvs)
      call getm4(nodevs,numnod)
      call getm4(ndiag,nstop)
      call getm4(nmoffc,nstop)
      call getm4(numoff,nstop)
      call getm4(irpt,nstop)
      call getm4(jcpt,nstop)
      call getm4(irowno,nstop)
      call getm4(jcolno,nstop)
      call slpmem(irpt,nstop)
      call slpmem(jcpt,nstop)
      call slpmem(irowno,nstop)
      call slpmem(jcolno,nstop)
      call crunch
c
      call zero4(nodplc(irpt+1),nstop)
      call zero4(nodplc(jcpt+1),nstop)
      call zero4(nodplc(irowno+1),nstop)
      call zero4(nodplc(jcolno+1),nstop)
      call zero4(nodplc(iseq1+1),numvs)
      call zero4(nodplc(nodevs+1),numnod)
      call zero4(nodplc(ndiag+1),nstop)
      call zero4(nodplc(nmoffc+1),nstop)
      call zero4(nodplc(numoff+1),nstop)
c
      numvs=0
      nxtrm=0
      ndist=0
      ntlin=1
      ibr=numnod
c
c  resistors
c
      loc=locate(1)
  110 if ((loc.eq.0).or.(nodplc(loc+8).ne.0)) go to 120
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      call reserv(node1,node1)
      call reserv(node1,node2)
      call reserv(node2,node1)
      call reserv(node2,node2)
      loc=nodplc(loc)
      go to 110
c
c  capacitors
c
  120 loc=locate(2)
  130 if ((loc.eq.0).or.(nodplc(loc+12).ne.0)) go to 400
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      call reserv(node1,node2)
      call reserv(node2,node1)
      ntemp=nodplc(ndiag+node1)
      call reserv(node1,node1)
      nodplc(ndiag+node1)=ntemp
      ntemp=nodplc(ndiag+node2)
      call reserv(node2,node2)
      nodplc(ndiag+node2)=ntemp
      nodplc(loc+8)=nxtrm+1
      nxtrm=nxtrm+2
      loc=nodplc(loc)
      go to 130
c
c  inductors
c
  400 loc=locate(3)
  430 if ((loc.eq.0).or.(nodplc(loc+14).ne.0)) go to 440
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ibr=ibr+1
      nodplc(loc+5)=ibr
      call reserv(node1,ibr)
      call reserv(node2,ibr)
      call reserv(ibr,node1)
      call reserv(ibr,node2)
      ntemp=nodplc(ndiag+ibr)
      call reserv(ibr,ibr)
      nodplc(ndiag+ibr)=ntemp
      numvs=numvs+1
      nodplc(iseq+numvs)=loc
      nodplc(neqn+numvs)=ibr
      nodplc(nodevs+node1)=nodplc(nodevs+node1)+1
      nodplc(nodevs+node2)=nodplc(nodevs+node2)+1
      nodplc(loc+11)=nxtrm+1
      nxtrm=nxtrm+2
      loc=nodplc(loc)
      go to 430
c
c  mutual inductors
c
  440 loc=locate(4)
  450 if ((loc.eq.0).or.(nodplc(loc+6).ne.0)) go to 460
      nl1=nodplc(loc+2)
      nl2=nodplc(loc+3)
      nl1=nodplc(nl1+5)
      nl2=nodplc(nl2+5)
      call reserv(nl1,nl2)
      call reserv(nl2,nl1)
      loc=nodplc(loc)
      go to 450
c
c  nonlinear voltage-controlled current sources
c
  460 loc=locate(5)
  462 if ((loc.eq.0).or.(nodplc(loc+13).ne.0)) go to 464
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ndim=nodplc(loc+4)
      ndim2=ndim+ndim
      locn=nodplc(loc+6)
      do 463 i=1,ndim2
      node=nodplc(locn+i)
      call reserv(node1,node)
      call reserv(node2,node)
  463 continue
      nodplc(loc+12)=nxtrm+1
      nxtrm=nxtrm+1+ndim2
      loc=nodplc(loc)
      go to 462
c
c  nonlinear voltage controlled voltage sources
c
  464 loc=locate(6)
  466 if ((loc.eq.0).or.(nodplc(loc+14).ne.0)) go to 468
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ibr=ibr+1
      nodplc(loc+6)=ibr
      call reserv(node1,ibr)
      call reserv(node2,ibr)
      call reserv(ibr,node1)
      call reserv(ibr,node2)
      numvs=numvs+1
      nodplc(iseq+numvs)=loc
      nodplc(neqn+numvs)=ibr
      nodplc(nodevs+node1)=nodplc(nodevs+node1)+1
      nodplc(nodevs+node2)=nodplc(nodevs+node2)+1
      ndim=nodplc(loc+4)
      ndim2=ndim+ndim
      locn=nodplc(loc+7)
      do 467 i=1,ndim2
      node=nodplc(locn+i)
      call reserv(ibr,node)
  467 continue
      nodplc(loc+13)=nxtrm+1
      nxtrm=nxtrm+2+ndim2
      loc=nodplc(loc)
      go to 466
c
c  voltage sources
c
  468 loc=locate(9)
  470 if ((loc.eq.0).or.(nodplc(loc+11).ne.0)) go to 472
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ibr=ibr+1
      nodplc(loc+6)=ibr
      call reserv(node1,ibr)
      call reserv(node2,ibr)
      call reserv(ibr,node1)
      call reserv(ibr,node2)
      numvs=numvs+1
      nodplc(iseq+numvs)=loc
      nodplc(neqn+numvs)=ibr
      nodplc(nodevs+node1)=nodplc(nodevs+node1)+1
      nodplc(nodevs+node2)=nodplc(nodevs+node2)+1
      loc=nodplc(loc)
      go to 470
c
c  nonlinear current controlled current sources
c
  472 loc=locate(7)
  474 if ((loc.eq.0).or.(nodplc(loc+13).ne.0)) go to 476
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ndim=nodplc(loc+4)
      locvs=nodplc(loc+6)
      do 475 i=1,ndim
      locvst=nodplc(locvs+i)
      kbr=nodplc(locvst+6)
      call reserv(node1,kbr)
      call reserv(node2,kbr)
  475 continue
      nodplc(loc+12)=nxtrm+1
      nxtrm=nxtrm+1+ndim+ndim
      loc=nodplc(loc)
      go to 474
c
c  nonlinear current controlled voltage sources
c
  476 loc=locate(8)
  478 if ((loc.eq.0).or.(nodplc(loc+14).ne.0)) go to 500
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      ibr=ibr+1
      nodplc(loc+6)=ibr
      call reserv(node1,ibr)
      call reserv(node2,ibr)
      call reserv(ibr,node1)
      call reserv(ibr,node2)
      numvs=numvs+1
      nodplc(iseq+numvs)=loc
      nodplc(neqn+numvs)=ibr
      nodplc(nodevs+node1)=nodplc(nodevs+node1)+1
      nodplc(nodevs+node2)=nodplc(nodevs+node2)+1
      ndim=nodplc(loc+4)
      locvs=nodplc(loc+7)
      do 479 i=1,ndim
      locvst=nodplc(locvs+i)
      kbr=nodplc(locvst+6)
      call reserv(ibr,kbr)
  479 continue
      nodplc(loc+13)=nxtrm+1
      nxtrm=nxtrm+2+ndim+ndim
      loc=nodplc(loc)
      go to 478
c
c  diodes
c
  500 loc=locate(11)
  510 if ((loc.eq.0).or.(nodplc(loc+16).ne.0)) go to 520
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      call reserv(node1,node1)
      call reserv(node2,node2)
      call reserv(node3,node3)
      call reserv(node1,node3)
      call reserv(node2,node3)
      call reserv(node3,node1)
      call reserv(node3,node2)
      nodplc(loc+11)=nxtrm+1
      nxtrm=nxtrm+5
      nodplc(loc+12)=ndist+1
      ndist=ndist+7
      loc=nodplc(loc)
      go to 510
c
c  transistors
c
  520 loc=locate(12)
  530 if ((loc.eq.0).or.(nodplc(loc+36).ne.0)) go to 540
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      node4=nodplc(loc+5)
      node5=nodplc(loc+6)
      node6=nodplc(loc+7)
      node7=nodplc(loc+30)
      locm=nodplc(loc+8)
      locm=nodplc(locm+1)
      cdis=value(locm+32)
      call reserv(node1,node1)
      call reserv(node2,node2)
      call reserv(node3,node3)
      call reserv(node4,node4)
      call reserv(node5,node5)
      call reserv(node6,node6)
      call reserv(node1,node4)
      call reserv(node2,node5)
      call reserv(node3,node6)
      call reserv(node4,node5)
      call reserv(node4,node6)
      call reserv(node5,node6)
      call reserv(node4,node1)
      call reserv(node5,node2)
      call reserv(node6,node3)
      call reserv(node5,node4)
      call reserv(node6,node4)
      call reserv(node6,node5)
      call reserv(node7,node7)
      call reserv(node4,node7)
      call reserv(node7,node4)
      if (cdis.lt.1.0d0) call reserv(node2,node4)
      if (cdis.lt.1.0d0) call reserv(node4,node2)
      nodplc(loc+22)=nxtrm+1
      nxtrm=nxtrm+19
      nodplc(loc+23)=ndist+1
      ndist=ndist+21
      loc=nodplc(loc)
      go to 530
c
c  jfets
c
  540 loc=locate(13)
  550 if ((loc.eq.0).or.(nodplc(loc+25).ne.0)) go to 560
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      node4=nodplc(loc+5)
      node5=nodplc(loc+6)
      call reserv(node1,node1)
      call reserv(node2,node2)
      call reserv(node3,node3)
      call reserv(node4,node4)
      call reserv(node5,node5)
      call reserv(node1,node4)
      call reserv(node2,node4)
      call reserv(node2,node5)
      call reserv(node3,node5)
      call reserv(node4,node5)
      call reserv(node4,node1)
      call reserv(node4,node2)
      call reserv(node5,node2)
      call reserv(node5,node3)
      call reserv(node5,node4)
      nodplc(loc+19)=nxtrm+1
      nxtrm=nxtrm+13
      loc=nodplc(loc)
      go to 550
c
c  mosfets
c
  560 loc=locate(14)
  570 if ((loc.eq.0).or.(nodplc(loc+33).ne.0)) go to 600
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      node4=nodplc(loc+5)
      node5=nodplc(loc+6)
      node6=nodplc(loc+7)
      call reserv(node1,node1)
      call reserv(node2,node2)
      call reserv(node3,node3)
      call reserv(node4,node4)
      call reserv(node5,node5)
      call reserv(node6,node6)
      call reserv(node1,node5)
      call reserv(node2,node4)
      call reserv(node2,node5)
      call reserv(node2,node6)
      call reserv(node3,node6)
      call reserv(node4,node5)
      call reserv(node4,node6)
      call reserv(node5,node6)
      call reserv(node5,node1)
      call reserv(node4,node2)
      call reserv(node5,node2)
      call reserv(node6,node2)
      call reserv(node6,node3)
      call reserv(node5,node4)
      call reserv(node6,node4)
      call reserv(node6,node5)
      nodplc(loc+26)=nxtrm+1
      nxtrm=nxtrm+28
      loc=nodplc(loc)
      go to 570
c
c  transmission lines
c
  600 loc=locate(17)
  610 if ((loc.eq.0).or.(nodplc(loc+33).ne.0)) go to 1000
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      node4=nodplc(loc+5)
      ni1=nodplc(loc+6)
      ni2=nodplc(loc+7)
      ibr1=ibr+1
      ibr2=ibr+2
      ibr=ibr+2
      nodplc(loc+8)=ibr1
      nodplc(loc+9)=ibr2
      call reserv(node1,node1)
      call reserv(node1,ni1)
      call reserv(node2,ibr1)
      call reserv(node3,node3)
      call reserv(node4,ibr2)
      call reserv(ni1,node1)
      call reserv(ni1,ni1)
      call reserv(ni1,ibr1)
      call reserv(ni2,ni2)
      call reserv(ni2,ibr2)
      call reserv(ibr1,node2)
      call reserv(ibr1,node3)
      call reserv(ibr1,node4)
      call reserv(ibr1,ni1)
      call reserv(ibr1,ibr2)
      call reserv(ibr2,node1)
      call reserv(ibr2,node2)
      call reserv(ibr2,node4)
      call reserv(ibr2,ni2)
      call reserv(ibr2,ibr1)
      call reserv(node3,ni2)
      call reserv(ni2,node3)
      numvs=numvs+1
      nodplc(iseq+numvs)=loc
      nodplc(iseq1+numvs)=1
      nodplc(neqn+numvs)=ibr1
      nodplc(nodevs+ni1)=nodplc(nodevs+ni1)+1
      nodplc(nodevs+node2)=nodplc(nodevs+node2)+1
      numvs=numvs+1
      nodplc(iseq+numvs)=loc
      nodplc(iseq1+numvs)=2
      nodplc(neqn+numvs)=ibr2
      nodplc(nodevs+ni2)=nodplc(nodevs+ni2)+1
      nodplc(nodevs+node4)=nodplc(nodevs+node4)+1
      nodplc(loc+30)=ntlin+1
      ntlin=ntlin+2
      loc=nodplc(loc)
      go to 610
c
c  finished
c
 1000 return
      end
      subroutine memadj
      implicit double precision (a-h,o-z)
c spice version 2g.6  sccsid=memmgr 3/15/83
      common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
     1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
     2   nwd8,nwd16
   50 maxuse=max0(maxuse,(ldval-memavl-ifwa))
      memdec=2*nxtmem(1)
      if (memavl.lt.memdec) return
c...  compress current allocations of memory
      call comprs(0,ldval)
c...  adjust memory size
      memdel=0
   60 icore=icore-memdec
      memdel=memdel+memdec
      memavl=memavl-memdec
      if (memavl.ge.memdec) go to 60
      ltab1=ldval-ntab
      istack(ltab1+2)=istack(ltab1+2)-memdel
c...  relocate block entry table
      nwords=numblk*ntab
      cpyknt=cpyknt+dble(nwords)
      call copy4(istack(loctab+1),istack(loctab-memdel+1),nwords)
      loctab=loctab-memdel
      ldval=ldval-memdel
      call memory
      return
      end
      subroutine memory
      implicit double precision (a-h,o-z)
c spice version 2g.6  sccsid=memmgr 3/15/83
      common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
     1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
     2   nwd8,nwd16
c spice version 2g.6  sccsid=flags 3/15/83
      common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
     1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,itl6,igoof,nogo,keof
      if(icore.le.maxmem) go to 10
    5 memerr=4
      return
   10 continue
c
c  set job field length to 'icore+locf(istack(ifwa))'
c
      return
      end
      logical function memptr(ipntr)
      implicit double precision (a-h,o-z)
c
c      this routine checks whether *ipntr* is a valid block pointer.
c if it is valid, *ltab* is set to point to the corresponding entry in
c the block table.
c
c... ipntr is an array to avoid 'call by value' problems (see setmem)
      dimension ipntr(1)
c spice version 2g.6  sccsid=memmgr 3/15/83
      common /memmgr/ cpyknt,istack(1),lorg,icore,maxcor,maxuse,memavl,
     1   ldval,numblk,loctab,ltab,ifwa,nwoff,ntab,maxmem,memerr,nwd4,
     2   nwd8,nwd16
      memptr=.false.
      ltab=loctab
      locpnt=locf(ipntr(1))
      do 20 i=1,numblk
      if (locpnt.ne.istack(ltab+4)) go to 10
      if (ipntr(1)*istack(ltab+5).ne.istack(ltab+1)) go to 10
      memptr=.true.
      go to 30
   10 ltab=ltab+ntab
   20 continue
   30 return
      end
      subroutine modchk
      implicit double precision (a-h,o-z)
c
c     this routine performs one-time processing of device model para-
c meters and prints out a device model summary.  it also reserves the
c additional nodes required by nonzero device extrinsic resistances.
c
c spice version 2g.6  sccsid=tabinf 3/15/83
      common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
     1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
     2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
     3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
     4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
     5   imynl,imvn,lcvn,nsnod,nsmat,nsval,icnod,icmat,icval,
     6   loutpt,lpol,lzer,irswpf,irswpr,icswpf,icswpr,irpt,jcpt,
     7   irowno,jcolno,nttbr,nttar,lvntmp
c spice version 2g.6  sccsid=miscel 3/15/83
      common /miscel/ atime,aprog(3),adate,atitle(10),defl,defw,defad,
     1  defas,rstats(50),iwidth,lwidth,nopage
c spice version 2g.6  sccsid=cirdat 3/15/83
      common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
     1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs,numalt,numcyc
c spice version 2g.6  sccsid=status 3/15/83
      common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
     1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
     2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
c spice version 2g.6  sccsid=flags 3/15/83
      common /flags/ iprnta,iprntl,iprntm,iprntn,iprnto,limtim,limpts,
     1   lvlcod,lvltim,itl1,itl2,itl3,itl4,itl5,itl6,igoof,nogo,keof
c spice version 2g.6  sccsid=knstnt 3/15/83
      common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
     1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox,
     2   pivtol,pivrel
c spice version 2g.6  sccsid=blank 3/15/83
      common /blank/ value(50000)
      integer nodplc(64)
      complex cvalue(32)
      equivalence (value(1),nodplc(1),cvalue(1))
c
c
      dimension itab(50),atable(12)
      dimension cpar(2),btable(12)
      dimension antype(4),aptype(4)
      dimension ipar(5),ampar(115),defval(115),ifmt(115),ivchk(115)
      dimension titled(4),titleb(4),titlej(4),titlem(4)
      data titled / 8hdiode mo, 8hdel para, 8hmeters  , 8h         /
      data titleb / 8hbjt mode, 8hl parame, 8hters    , 8h         /
      data titlej / 8hjfet mod, 8hel param, 8heters   , 8h         /
      data titlem / 8hmosfet m, 8hodel par, 8hameters , 8h         /
      data antype /1h ,3hnpn,3hnjf,4hnmos/
      data aptype /1h ,3hpnp,3hpjf,4hpmos/
      data ipar /0,14,60,72,114/
      data cpar / 3hc2 ,3hc4 /
      data aundef /2h.u/
      data ampar /
     1   6his    ,6hrs    ,6hn     ,6htt    ,6hcjo   ,6hvj    ,6hm     ,
     2   6heg    ,6hxti   ,6hkf    ,6haf    ,6hfc    ,6hbv    ,6hibv   ,
     1   6his    ,6hbf    ,6hnf    ,6hvaf   ,6hikf   ,6hise   ,6hne    ,
     2   6hbr    ,6hnr    ,6hvar   ,6hikr   ,6hisc   ,6hnc    ,6h0     ,
     3   6h0     ,6hrb    ,6hirb   ,6hrbm   ,6hre    ,6hrc    ,6hcje   ,
     4   6hvje   ,6hmje   ,6htf    ,6hxtf   ,6hvtf   ,6hitf   ,6hptf   ,
     5   6hcjc   ,6hvjc   ,6hmjc   ,6hxcjc  ,6htr    ,6h0     ,6h0     ,
     6   6h0     ,6h0     ,6hcjs   ,6hvjs   ,6hmjs   ,6hxtb   ,6heg    ,
     7   6hxti   ,6hkf    ,6haf    ,6hfc    ,
     1   6hvto   ,6hbeta  ,6hlambda,6hrd    ,6hrs    ,6hcgs   ,6hcgd   ,
     2   6hpb    ,6his    ,6hkf    ,6haf    ,6hfc    ,
     1   6hlevel ,6hvto   ,6hkp    ,6hgamma ,6hphi   ,6hlambda,6hrd    ,
     2   6hrs    ,6hcbd   ,6hcbs   ,6his    ,6hpb    ,6hcgso  ,6hcgdo  ,
     3   6hcgbo  ,6hrsh   ,6hcj    ,6hmj    ,6hcjsw  ,6hmjsw  ,6hjs    ,
     4   6htox   ,6hnsub  ,6hnss   ,6hnfs   ,6htpg   ,6hxj    ,6hld    ,
     5   6huo    ,6hucrit ,6huexp  ,6hutra  ,6hvmax  ,6hneff  ,6hxqc   ,
     6   6hkf    ,6haf    ,6hfc    ,6hdelta ,6htheta ,6heta   ,6hkappa ,
     7   0.0d0   /
      data defval /
     1   1.0d-14,  0.0d0,  1.0d0,2*0.0d0,  1.0d0,  0.5d0, 1.11d0,
     2     3.0d0,  0.0d0,  1.0d0,  0.5d0,  0.0d0, 1.0d-3,
     1   1.0d-16,100.0d0,  1.0d0,3*0.0d0,  1.5d0,2*1.0d0,3*0.0d0,
     2     2.0d0,  0.0d0,  1.0d0,6*0.0d0, 0.75d0, 0.33d0,2*0.0d0,
     3   4*0.0d0, 0.75d0, 0.33d0,  1.0d0,6*0.0d0, 0.75d0,2*0.0d0,
     4    1.11d0,  3.0d0,  0.0d0,  1.0d0,  0.5d0,
     1    -2.0d0, 1.0d-4,5*0.0d0,  1.0d0,1.0d-14,  0.0d0,  1.0d0,
     2     0.5d0,
     1     1.0d0,  0.0d0, 2.0d-5,  0.0d0,  0.6d0,5*0.0d0,1.0d-14,
     2     0.8d0,5*0.0d0,  0.5d0,  0.0d0, 0.33d0,5*0.0d0,  1.0d0,
     3   2*0.0d0,600.0d0, 1.0d+4,3*0.0d0,  1.0d0,  1.0d0,  0.0d0,
     4     1.0d0,  0.5d0,3*0.0d0,  0.2d0,
     5     0.0d0/
      data ifmt /
     1   4,1,1,2,2,1,1,1,1,2,1,1,2,2,
     2   4,3,3,2,2,2,1,3,3,2,2,2,1,0,0,1,2,1,1,1,2,1,1,2,2,2,2,1,2,1,
     2   1,1,2,0,0,0,0,2,1,1,2,1,1,2,2,2,
     3   3,4,1,1,1,2,2,1,2,2,1,1,
     4   3,3,4,1,1,2,1,1,2,2,2,1,2,2,2,1,2,1,2,1,2,2,2,2,2,1,2,2,
     4   1,2,1,1,2,1,1,2,1,1,1,1,1,1,
     5   0/
      data ivchk /
     1   0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     2   0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     2   0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,0,0,
     3   -1,0,0,0,0,0,0,0,0,0,0,0,
     4   0,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,-1,0,-1,
     4   0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
     5   0/
c
c
      tnom=value(itemps+1)+ctok
      xkt=boltz*tnom
      vt=xkt/charge
      xni=1.45d16
      egfet=1.16d0-(7.02d-4*tnom*tnom)/(tnom+1108.0d0)
      nummod=jelcnt(21)+jelcnt(22)+jelcnt(23)+jelcnt(24)
      if (nummod.eq.0) go to 1000
c
c  special preprocessing for mosfet models
c
      loc=locate(24)
    5 if (loc.eq.0) go to 35
      locv=nodplc(loc+1)
      type=nodplc(loc+2)
c
c     default parameters for higher level mos models
c
      lev=value(locv+1)
      if (value(locv+1).eq.aundef) lev=1
      if (value(locv+23).ne.aundef) xnsub=value(locv+23)*1.0d6
      if (value(locv+22).eq.aundef.and.lev.gt.1) value(locv+22)=1.0d-7
      if (value(locv+22).eq.aundef) go to 33
      cox=epsox/value(locv+22)
c
c     compute kp, if not input, using default mobility 600 cm**2/v*sec
c
      if (value(locv+3).ne.aundef) go to 10
      if (value(locv+29).eq.aundef) value(locv+29)=600.0d0
      value(locv+3)=value(locv+29)*cox*1.0d-4
   10 if (value(locv+23).eq.aundef) go to 33
      if (xnsub.le.xni) go to 30
c
c     nsub nonzero => process oriented model
c
      if (value(locv+5).eq.aundef) value(locv+5)=
     1   dmax1((2.0d0*vt*dlog(xnsub/xni)),0.1d0)
      fermis=type*0.5d0*value(locv+5)
      wkfng=3.2d0
      if (value(locv+26).eq.aundef) value(locv+26)=1.0d0
      if (value(locv+26).eq.0.0d0) go to 15
c
c  polysilicon gate
c
      fermig=type*value(locv+26)*0.5d0*egfet
      wkfng=3.25d0+0.5d0*egfet-fermig
   15 wkfngs=wkfng-(3.25d0+0.5d0*egfet+fermis)
      if (value(locv+4).eq.aundef)
     1   value(locv+4)=dsqrt(2.0d0*epssil*charge*xnsub)/cox
c
c     computed vto
c
      if (value(locv+2).ne.aundef) go to 20
      if (value(locv+24).eq.aundef) value(locv+24)=0.0d0
      value(locv+44)=wkfngs-value(locv+24)*1.0d4*charge/cox
      value(locv+2)=value(locv+44)
     1   +type*(value(locv+4)*dsqrt(value(locv+5))+value(locv+5))
      go to 25
c
c     measured vto has been input
c
   20 value(locv+44)=value(locv+2)
     1   -type*(value(locv+4)*dsqrt(value(locv+5))+value(locv+5))
   25 value(locv+45)=dsqrt((epssil+epssil)/(charge*xnsub))
      go to 33
   30 value(locv+23)=0.0d0
      write (iofile,31) value(locv)
   31 format('0*error*:  nsub <= ni in mosfet model ',a8,/)
      nogo=1
c
c   special processing for mos3: limit kappa>0,
c   set to zero lambda,ucrit,uexp and utar
c
   33 if (lev.ne.3) go to 34
      if (value(locv+42).eq.aundef) value(locv+42)=0.2d0
      value(locv+6)=0.0d0
      value(locv+30)=0.0d0
      value(locv+31)=0.0d0
      value(locv+32)=0.0d0
   34 loc=nodplc(loc)
      go to 5
c
c     cycle thru devices
c
   35 kntlim=lwidth/11
      do 390 id=1,4
      if (jelcnt(id+20).eq.0) go to 390
      locm=ipar(id)
      nopar=ipar(id+1)-locm
      do 45 i=1,nopar
      if (ifmt(locm+i).ge.3) go to 40
      itab(i)=0
      go to 45
   40 itab(i)=ifmt(locm+i)-2
   45 continue
c
c  assign default values
c
      loc=locate(id+20)
   50 if (loc.eq.0) go to 70
      locv=nodplc(loc+1)
      do 65 i=1,nopar
      if (value(locv+i).eq.aundef) go to 62
      if (ivchk(locm+i).lt.0) go to 55
      if (value(locv+i).lt.0.0d0) go to 62
   55 if (itab(i).ne.0) go to 65
      itab(i)=ifmt(locm+i)
      go to 65
   62 value(locv+i)=defval(locm+i)
   65 continue
      loc=nodplc(loc)
      go to 50
c
c     limit model values
c
   70 go to (80,85,90,95), id
c...  diodes
   80 loc=locate(21)
   82 if (loc.eq.0) go to 130
      locv=nodplc(loc+1)
      value(locv+7)=dmin1(value(locv+7),0.9d0)
      value(locv+8)=dmax1(value(locv+8),0.1d0)
      value(locv+11)=dmax1(value(locv+11),0.1d0)
      value(locv+12)=dmin1(value(locv+12),0.95d0)
      loc=nodplc(loc)
      go to 82
c...  bipolar transistors
   85 loc=locate(22)
   87 if (loc.eq.0) go to 130
      locv=nodplc(loc+1)
      value(locv+23)=dmin1(value(locv+23),0.9d0)
      if (value(locv+24).eq.0.0d0) value(locv+28)=0.0d0
      value(locv+31)=dmin1(value(locv+31),0.9d0)
      value(locv+32)=dmin1(value(locv+32),1.0d0)
      value(locv+40)=dmin1(value(locv+40),0.9d0)
      value(locv+42)=dmax1(value(locv+42),0.1d0)
      value(locv+45)=dmax1(value(locv+45),0.1d0)
      value(locv+46)=dmin1(value(locv+46),0.9999d0)
      loc=nodplc(loc)
      if (value(locv+18).eq.0.0d0) value(locv+18)=value(locv+16)
      if (value(locv+16).ge.value(locv+18)) go to 87
      write(iofile,89) value(locv)
   89 format('0warning:  minimum base resistance (rbm) is less than '
     1       ,'total (rb) for model ',a8,/10x,' rbm set equal to rb',/)
      value(locv+18)=value(locv+16)
      go to 87
c...  jfets
   90 loc=locate(23)
   92 if (loc.eq.0) go to 130
      locv=nodplc(loc+1)
      value(locv+11)=dmax1(value(locv+11),0.1d0)
      value(locv+12)=dmin1(value(locv+12),0.95d0)
      loc=nodplc(loc)
      go to 92
c...  mosfets
   95 loc=locate(24)
   97 if (loc.eq.0) go to 130
      locv=nodplc(loc+1)
c
  100 value(locv+37)=dmax1(value(locv+37),0.1d0)
      value(locv+38)=dmin1(value(locv+38),0.95d0)
      if (value(locv+23).le.0.0d0) go to 120
      cj=dsqrt(epssil*charge*value(locv+23)*1.0d6/
     1   (2.0d0*value(locv+12)))
      if (value(locv+9).le.0.0d0) go to 105
      itab(9)=2
  105 if (value(locv+10).le.0.0d0) go to 110
      itab(10)=2
      go to 115
  110 if (value(locv+17).le.0.0d0) value(locv+17)=cj
      itab(17)=2
  115 if ((value(locv+7).le.0.0d0).and.
     1    (value(locv+8).le.0.0d0)) go to 120
      itab(7)=2
      itab(8)=2
  120 if (value(locv+6).ge.0.2d0) write (iofile,121) value(locv)
  121 format ('0warning:  the value of lambda for mosfet model ',a8,/,
     1   ' is unusually large and might cause nonconvergence',/)
      if (lev.ne.2) value(locv+35)=1.0d0
      if (lev.ne.3) go to 125
      itab(40)=1
      itab(41)=1
      itab(42)=1
      itab(43)=1
  125 loc=nodplc(loc)
      go to 97
c
c     print model parameters
c
  130 if (iprntm.eq.0) go to 360
      locs=locate(id+20)
  140 kntr=0
      loc=locs
      go to (150,160,170,180), id
  150 call title(0,lwidth,1,titled)
      go to 200
  160 call title(0,lwidth,1,titleb)
      go to 200
  170 call title(0,lwidth,1,titlej)
      go to 200
  180 call title(0,lwidth,1,titlem)
  200 if (loc.eq.0) go to 210
      if (kntr.lt.kntlim) go to 220
  210 locn=loc
      go to 240
  220 kntr=kntr+1
      locv=nodplc(loc+1)
      atable(kntr)=value(locv)
  230 loc=nodplc(loc)
      go to 200
  240 write (iofile,241) (atable(k),k=1,kntr)
  241 format(//11x,12(2x,a8))
      if (id.eq.1) go to 300
      kntr=0
      loc=locs
  250 if (loc.eq.0) go to 260
      if (kntr.ge.kntlim) go to 260
      kntr=kntr+1
      atable(kntr)=antype(id)
      if (nodplc(loc+2).eq.-1) atable(kntr)=aptype(id)
      loc=nodplc(loc)
      go to 250
  260 write (iofile,261) (atable(k),k=1,kntr)
  261 format('0type',4x,12(4x,a6))
  300 do 340 i=1,nopar
      if (itab(i).eq.0) go to 340
      kntr=0
      iccflg=0
      loc=locs
  310 if (loc.eq.0) go to 320
      if (kntr.ge.kntlim) go to 320
      locv=nodplc(loc+1)
      kntr=kntr+1
      if (iccflg.ne.0) go to 313
      if (id.ne.2) go to 315
      if ((i.ne.6).and.(i.ne.12)) go to 315
      if (value(locv+i).le.1.0d0) go to 315
      iccflg=i/6
  313 btable(kntr)=value(locv+i)
      value(locv+i)=value(locv+i)*value(locv+1)
  315 atable(kntr)=value(locv+i)
      loc=nodplc(loc)
      go to 310
  320 if (itab(i).eq.2) go to 330
      write (iofile,321) ampar(locm+i),(atable(k),k=1,kntr)
  321 format(1h0,a8,12f10.3)
      go to 340
  330 write (iofile,331) ampar(locm+i),(atable(k),k=1,kntr)
  331 format(1h0,a8,1p12d10.2)
      if (iccflg.eq.0) go to 340
      write (iofile,321) cpar(iccflg),(btable(k),k=1,kntr)
  340 continue
      if (locn.eq.0) go to 390
      locs=locn
      go to 140
c
c  special  treatment for c2 & c4 in the bjt model
c  when no model parameter print
c
  360 if (id.ne.2) go to 390
      loc=locate(id+20)
  370 if (loc.eq.0) go to 390
      locv=nodplc(loc+1)
      if (value(locv+6).ge.1.0d0)
     1   value(locv+6)=value(locv+6)*value(locv+1)
      if (value(locv+12).ge.1.0d0)
     1   value(locv+12)=value(locv+12)*value(locv+1)
      loc=nodplc(loc)
      go to 370
  390 continue
c
c  process model parameters
c
c  diodes
c
  400 loc=locate(21)
  410 if (loc.eq.0) go to 420
      locv=nodplc(loc+1)
      if (value(locv+2).ne.0.0d0) value(locv+2)=1.0d0/value(locv+2)
      pb=value(locv+6)
      xm=value(locv+7)
      fc=value(locv+12)
      value(locv+12)=fc*pb
      xfc=dlog(1.0d0-fc)
      value(locv+15)=pb*(1.0d0-dexp((1.0d0-xm)*xfc))/(1.0d0-xm)
      value(locv+16)=dexp((1.0d0+xm)*xfc)
      value(locv+17)=1.0d0-fc*(1.0d0+xm)
      csat=value(locv+1)
      vte=value(locv+3)*vt
      value(locv+18)=vte*dlog(vte/(root2*csat))
      bv=value(locv+13)
      if (bv.eq.0) go to 418
      cbv=value(locv+14)
      if (cbv.ge.csat*bv/vt) go to 412
      cbv=csat*bv/vt
      write (iofile,411) value(locv),cbv
  411 format('0warning:  in diode model ',a8,' ibv increased to ',1pe10.
     1   3,11x,'to resolve incompatibility with specified is'/)
      xbv=bv
      go to 416
  412 tol=reltol*cbv
      xbv=bv-vt*dlog(1.0d0+cbv/csat)
      iter=0
  413 xbv=bv-vt*dlog(cbv/csat+1.0d0-xbv/vt)
      xcbv=csat*(dexp((bv-xbv)/vt)-1.0d0+xbv/vt)
      if (dabs(xcbv-cbv).le.tol) go to 416
      iter=iter+1
      if (iter.lt.25) go to 413
      write (iofile,415) xbv,xcbv
  415 format('0warning:  unable to match forward and reverse diode regio
     1ns',/,11x,'bv = ',1pd10.3,' and ibv = ',d10.3,/)
  416 value(locv+13)=xbv
  418 loc=nodplc(loc)
      go to 410
c
c  bipolar transistor models
c
  420 loc=locate(22)
  430 if (loc.eq.0) go to 440
      locv=nodplc(loc+1)
      if (value(locv+4).ne.0.0d0) value(locv+4)=1.0d0/value(locv+4)
      if (value(locv+5).ne.0.0d0) value(locv+5)=1.0d0/value(locv+5)
      if (value(locv+10).ne.0.0d0) value(locv+10)=1.0d0/value(locv+10)
      if (value(locv+11).ne.0.0d0) value(locv+11)=1.0d0/value(locv+11)
      if (value(locv+19).ne.0.0d0) value(locv+19)=1.0d0/value(locv+19)
      if (value(locv+20).ne.0.0d0) value(locv+20)=1.0d0/value(locv+20)
      if (value(locv+26).ne.0.0d0) value(locv+26)=1.0d0/value(locv+26)
     1   /1.44d0
      value(locv+28)=value(locv+28)/rad*value(locv+24)
      if (value(locv+35).ne.0.0d0) value(locv+35)=1.0d0/value(locv+35)
     1   /1.44d0
      pe=value(locv+22)
      xme=value(locv+23)
      pc=value(locv+30)
      xmc=value(locv+31)
      fc=value(locv+46)
      value(locv+46)=fc*pe
      xfc=dlog(1.0d0-fc)
      value(locv+47)=pe*(1.0d0-dexp((1.0d0-xme)*xfc))/(1.0d0-xme)
      value(locv+48)=dexp((1.0d0+xme)*xfc)
      value(locv+49)=1.0d0-fc*(1.0d0+xme)
      value(locv+50)=fc*pc
      value(locv+51)=pc*(1.0d0-dexp((1.0d0-xmc)*xfc))/(1.0d0-xmc)
      value(locv+52)=dexp((1.0d0+xmc)*xfc)
      value(locv+53)=1.0d0-fc*(1.0d0+xmc)
      csat=value(locv+1)
      value(locv+54)=vt*dlog(vt/(root2*csat))
      loc=nodplc(loc)
      go to 430
c
c  jfet models
c
  440 loc=locate(23)
  450 if (loc.eq.0) go to 460
      locv=nodplc(loc+1)
      if (value(locv+4).ne.0.0d0) value(locv+4)=1.0d0/value(locv+4)
      if (value(locv+5).ne.0.0d0) value(locv+5)=1.0d0/value(locv+5)
      pb=value(locv+8)
      xm=0.5d0
      fc=value(locv+12)
      value(locv+12)=fc*pb
      xfc=dlog(1.0d0-fc)
      value(locv+13)=pb*(1.0d0-dexp((1.0d0-xm)*xfc))/(1.0d0-xm)
      value(locv+14)=dexp((1.0d0+xm)*xfc)
      value(locv+15)=1.0d0-fc*(1.0d0+xm)
      csat=value(locv+9)
      value(locv+16)=vt*dlog(vt/(root2*csat))
      loc=nodplc(loc)
      go to 450
c
c  mosfet models
c
  460 loc=locate(24)
  470 if (loc.eq.0) go to 600
      locv=nodplc(loc+1)
      type=nodplc(loc+2)
      if (value(locv+7).ne.0.0d0) value(locv+7)=1.0d0/value(locv+7)
      if (value(locv+8).ne.0.0d0) value(locv+8)=1.0d0/value(locv+8)
      if (value(locv+16).ne.0.0d0) value(locv+16)=1.0d0/value(locv+16)
      value(locv+23)=value(locv+23)*1.0d6
      value(locv+24)=value(locv+24)*1.0d4
      value(locv+25)=value(locv+25)*1.0d4
      if (value(locv+22).ne.0.0d0) value(locv+22)=epsox/value(locv+22)
      value(locv+29)=value(locv+29)*1.0d-4
      if (lev.eq.3) go to 472
      value(locv+30)=value(locv+30)*1.0d2
      go to 473
c
c   move mos3 parameters : theta from locations locv+40 to locv+30
c                          eta                       41         31
c                          kappa                     42         32
c   and replace locv+6 by (xd)**2
c
  472 value(locv+39)=value(locv+39)
     1   *0.25d0*twopi*epssil/value(locv+22)
      value(locv+30)=value(locv+40)
      value(locv+31)=value(locv+41)*8.15d-22/value(locv+22)
      value(locv+32)=value(locv+42)
      if (value(locv+23).gt.0.0d0)
     1    value(locv+6)=(epssil+epssil)/(charge*value(locv+23))
c
c   noise parameters
c
  473 pb=value(locv+12)
      xm=0.5d0
      fc=value(locv+38)
      value(locv+38)=fc*pb
      xfc=dlog(1.0d0-fc)
      value(locv+40)=pb*(1.0d0-dexp((1.0d0-xm)*xfc))/(1.0d0-xm)
      value(locv+41)=dexp((1.0d0+xm)*xfc)
      value(locv+42)=1.0d0-fc*(1.0d0+xm)
      value(locv+43)=-1.0d0
      value(locv+44)=value(locv+2)-
     1   type*value(locv+4)*dsqrt(value(locv+5))
  475 if (value(locv+22).ne.0.0d0.and.lev.ne.3)
     1   value(locv+30)=value(locv+30)*epssil/value(locv+22)
      loc=nodplc(loc)
      go to 470
c
c  reserve additional nodes
c
c  diodes
c
  600 loc=locate(11)
  610 if ((loc.eq.0).or.(nodplc(loc+16).ne.0)) go to 700
      locm=nodplc(loc+5)
      locm=nodplc(locm+1)
      if (value(locm+2).eq.0.0d0) go to 620
      numnod=numnod+1
      nodplc(loc+4)=numnod
      go to 630
  620 nodplc(loc+4)=nodplc(loc+2)
  630 loc=nodplc(loc)
      go to 610
c
c  transistors
c
  700 loc=locate(12)
  710 if ((loc.eq.0).or.(nodplc(loc+36).ne.0)) go to 800
      nodplc(loc+30)=nodplc(loc+5)
      locm=nodplc(loc+8)
      locm=nodplc(locm+1)
      if (value(locm+16).eq.0.0d0) go to 720
      numnod=numnod+1
      nodplc(loc+6)=numnod
      go to 730
  720 nodplc(loc+6)=nodplc(loc+3)
  730 if (value(locm+20).eq.0.0d0) go to 740
      numnod=numnod+1
      nodplc(loc+5)=numnod
      go to 750
  740 nodplc(loc+5)=nodplc(loc+2)
  750 if (value(locm+19).eq.0.0d0) go to 760
      numnod=numnod+1
      nodplc(loc+7)=numnod
      go to 770
  760 nodplc(loc+7)=nodplc(loc+4)
  770 loc=nodplc(loc)
      go to 710
c
c  jfets
c
  800 loc=locate(13)
  810 if ((loc.eq.0).or.(nodplc(loc+25).ne.0)) go to 900
      locm=nodplc(loc+7)
      locm=nodplc(locm+1)
      if (value(locm+4).eq.0.0d0) go to 820
      numnod=numnod+1
      nodplc(loc+5)=numnod
      go to 830
  820 nodplc(loc+5)=nodplc(loc+2)
  830 if (value(locm+5).eq.0.0d0) go to 840
      numnod=numnod+1
      nodplc(loc+6)=numnod
      go to 850
  840 nodplc(loc+6)=nodplc(loc+4)
  850 loc=nodplc(loc)
      go to 810
c
c  mosfets
c
  900 loc=locate(14)
  910 if (loc.eq.0) go to 1000
      locm=nodplc(loc+8)
      locm=nodplc(locm+1)
      locv=nodplc(loc+1)
      xleff=value(locv+1)-2.0d0*value(locm+28)
      if (xleff.gt.0.0d0) go to 915
      write(iofile,911) value(locv),value(locm)
  911 format('0*error*:  effective channel length of ',a8,' less than ',
     1   'zero.',/' check value of ld for model ',a8)
      if (nodplc(loc+33).ne.0) go to 960
  915 if ((value(locm+7).eq.0.0d0).and.
     1    (value(locm+16).eq.0.0d0)) go to 920
      numnod=numnod+1
      nodplc(loc+6)=numnod
      go to 930
  920 nodplc(loc+6)=nodplc(loc+2)
  930 if ((value(locm+8).eq.0.0d0).and.
     1    (value(locm+16).eq.0.0d0)) go to 940
      numnod=numnod+1
      nodplc(loc+7)=numnod
      go to 950
  940 nodplc(loc+7)=nodplc(loc+4)
  950 ad=value(locv+3)
      as=value(locv+4)
      if ((ad.le.0.0d0).or.(as.le.0.0d0)
     1   .and.value(locm+11).le.0.0d0)
     2   value(locm+11)=1.0d-14
  960 loc=nodplc(loc)
      go to 910
c
c  transmission lines
c
 1000 loc=locate(17)
 1010 if ((loc.eq.0).or.(nodplc(loc+33).ne.0)) go to 2000
      numnod=numnod+1
      nodplc(loc+6)=numnod
      numnod=numnod+1
      nodplc(loc+7)=numnod
      loc=nodplc(loc)
      go to 1010
c
c  finished
c
 2000 return
      end
      subroutine moscap(vgd,vgs,vgb,covlgd,covlgs,covlgb,
     1   capbd,capbs,cggb,cgdb,cgsb,cbgb,cbdb,cbsb,
     2   gcggb,gcgdb,gcgsb,gcbgb,gcbdb,gcbsb,
     3   gcdgb,gcddb,gcdsb,gcsgb,gcsdb,gcssb,
     4   qgate,qchan,qbulk,qdrn,qsrc)
      implicit double precision (a-h,o-z)
c spice version 2g.6  sccsid=mosarg 3/15/83
      common /mosarg/ vto,beta,gamma,phi,phib,cox,xnsub,xnfs,xd,xj,xld,
     1   xlamda,uo,uexp,vbp,utra,vmax,xneff,xl,xw,vbi,von,vdsat,qspof,
     2   beta0,beta1,cdrain,xqco,xqc,fnarrw,fshort,lev
c spice version 2g.6  sccsid=status 3/15/83
      common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
     1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
     2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
c spice version 2g.6  sccsid=blank 3/15/83
      common /blank/ value(50000)
      integer nodplc(64)
      complex cvalue(32)
      equivalence (value(1),nodplc(1),cvalue(1))
c
c     compute equivalent conductances
c     divide up the channel charge (1-xqc)/xqc to source and drain
c
      gcg=(cggb+cbgb)*ag(1)
      gcd=(cgdb+cbdb)*ag(1)
      gcs=(cgsb+cbsb)*ag(1)
      gcgxd=-xqc*gcg
      gcgxs=-(1.0d0-xqc)*gcg
      gcdxd=-xqc*gcd
      gcdxs=-(1.0d0-xqc)*gcd
      gcsxd=-xqc*gcs
      gcsxs=-(1.0d0-xqc)*gcs
      gcdgb=gcgxd-covlgd*ag(1)
      gcddb=gcdxd+(capbd+covlgd)*ag(1)
      gcdsb=gcsxd
      gcsgb=gcgxs-covlgs*ag(1)
      gcsdb=gcdxs
      gcssb=gcsxs+(capbs+covlgs)*ag(1)
      gcggb=(cggb+covlgd+covlgs+covlgb)*ag(1)
      gcgdb=(cgdb-covlgd)*ag(1)
      gcgsb=(cgsb-covlgs)*ag(1)
      gcbgb=(cbgb-covlgb)*ag(1)
      gcbdb=(cbdb-capbd)*ag(1)
      gcbsb=(cbsb-capbs)*ag(1)
c
c     compute total terminal charges
c
      qgd=covlgd*vgd
      qgs=covlgs*vgs
      qgb=covlgb*vgb
      qgate=qgate+qgd+qgs+qgb
      qbulk=qbulk-qgb
      qdrn=xqc*qchan-qgd
      qsrc=(1.0d0-xqc)*qchan-qgs
c
c     finished
c
      return
      end
      subroutine moseq1(vds,vbs,vgs,gm,gds,gmbs)
      implicit double precision (a-h,o-z)
c
c     this routine evaluates the drain current and its derivatives
c     using the shichman-hodges model and the charges associated
c     with the gate, channel and bulk for mosfets
c
c spice version 2g.6  sccsid=mosarg 3/15/83
      common /mosarg/ vto,beta,gamma,phi,phib,cox,xnsub,xnfs,xd,xj,xld,
     1   xlamda,uo,uexp,vbp,utra,vmax,xneff,xl,xw,vbi,von,vdsat,qspof,
     2   beta0,beta1,cdrain,xqco,xqc,fnarrw,fshort,lev
c spice version 2g.6  sccsid=status 3/15/83
      common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
     1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
     2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
c spice version 2g.6  sccsid=knstnt 3/15/83
      common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
     1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox,
     2   pivtol,pivrel
      vbd=vbs-vds
      vgb=vgs-vbs
c
c
      if (vbs.gt.0.0d0) go to 102
      sarg=dsqrt(phi-vbs)
      go to 104
  102 sarg=dsqrt(phi)
      sarg=sarg-vbs/(sarg+sarg)
      sarg=dmax1(0.0d0,sarg)
  104 von=vbi+gamma*sarg
      vgst=vgs-von
      vdsat=dmax1(vgst,0.0d0)
      if (sarg.gt.0.0d0) go to 105
      arg=0.0d0
      go to 108
  105 arg=gamma/(sarg+sarg)
  108 if (vgst.gt.0.0d0) go to 110
c
c     cutoff region
c
      cdrain=0.0d0
      gm=0.0d0
      gds=0.0d0
      gmbs=0.0d0
      go to 1000
c
c     saturation region
c
  110 betap=beta*(1.0d0+xlamda*vds)
      if (vgst.gt.vds) go to 120
      cdrain=betap*vgst*vgst*0.5d0
      gm=betap*vgst
      gds=xlamda*beta*vgst*vgst*0.5d0
      gmbs=gm*arg
      go to 1000
c
c     linear region
c
  120 cdrain=betap*vds*(vgst-0.5d0*vds)
      gm=betap*vds
      gds=betap*(vgst-vds)+xlamda*beta*vds*(vgst-0.5d0*vds)
      gmbs=gm*arg
c
c     finished
c
 1000 return
      end
      subroutine moseq2(vds,vbs,vgs,gm,gds,gmbs,
     1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      implicit double precision (a-h,o-z)
c
c     this routine evaluates the drain current, its derivatives and
c     the charges associated with the gate, channel and bulk
c     for mosfets
c
c spice version 2g.6  sccsid=mosarg 3/15/83
      common /mosarg/ vto,beta,gamma,phi,phib,cox,xnsub,xnfs,xd,xj,xld,
     1   xlamda,uo,uexp,vbp,utra,vmax,xneff,xl,xw,vbi,von,vdsat,qspof,
     2   beta0,beta1,cdrain,xqco,xqc,fnarrw,fshort,lev
c spice version 2g.6  sccsid=status 3/15/83
      common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
     1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
     2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
c spice version 2g.6  sccsid=knstnt 3/15/83
      common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
     1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox,
     2   pivtol,pivrel
c
      dimension a4(4),b4(4),x4(8),poly4(8),sig1(4),sig2(4)
      data sig1 / 1.0d0, -1.0d0, 1.0d0, -1.0d0/,
     1     sig2 / 1.0d0,  1.0d0,-1.0d0, -1.0d0/
c
c     icharg=1 causes charges to be computed
c     icharg=0 bypasses the computation of charges
c
      icharg=1
      if (mode.ne.1.and.xqco.le.0.5d0) go to 100
      icharg=0
      if (xqco.gt.0.5d0) go to 100
      if (modedc.eq.2.and.nosolv.ne.0) icharg=1
      if (initf.eq.4) icharg=1
c
c  compute some useful quantities
c
  100 if (vbs.gt.0.0d0) go to 110
      sarg=dsqrt(phi-vbs)
      tsarg=sarg+sarg
      dsrgdb=-0.5d0/sarg
      d2sdb2=+0.5d0*dsrgdb/(phi-vbs)
      go to 120
  110 sphi=dsqrt(phi)
      sphi3=phi*sphi
      sarg=sphi/(1.0d0+0.5d0*vbs/phi)
      tsarg=sarg+sarg
      dsrgdb=-0.5d0*sarg*sarg/sphi3
      d2sdb2=-dsrgdb*sarg/sphi3
  120 if ((vds-vbs).lt.0.0d0) go to 130
      barg=dsqrt(phi+vds-vbs)
      dbrgdb=-0.5d0/barg
      d2bdb2=+0.5d0*dbrgdb/(phi+vds-vbs)
      go to 200
  130 barg=sphi/(1.0d0+0.5d0*(vbs-vds)/phi)
      dbrgdb=-0.5d0*barg*barg/sphi3
      d2bdb2=-dbrgdb*barg/sphi3
c
c  calculate threshold voltage (von)
c     narrow-channel effect
c
  200 factor=0.125d0*fnarrw*twopi*epssil/cox*xl
      eta=1.0d0+factor
      vbin=vbi+factor*(phi-vbs)
      if (gamma.le.0.0d0) go to 215
      if (xnsub.le.0.0d0) go to 215
      xwd=xd*barg
      xws=xd*sarg
c
c     short-channel effect with vds .ne. 0.0d0
c
      argss=0.0d0
      argsd=0.0d0
      dbargs=0.0d0
      dbargd=0.0d0
      dgdvds=0.0d0
      dgddb2=0.0d0
      if (xj.le.0.0d0) go to 205
      argxs=1.0d0+2.0d0*xws/xj
      args=dsqrt(argxs)
      argss=0.5d0*xj/xl*(args-1.0d0)
      argxd=1.0d0+2.0d0*xwd/xj
      argd=dsqrt(argxd)
      argsd=0.5d0*xj/xl*(argd-1.0d0)
  205 gamasd=gamma*(1.0d0-argss-argsd)
      gamass=gamma*(1.0d0-2.0d0*argss)
      dbxwd=xd*dbrgdb
      dbxws=xd*dsrgdb
      if (xj.le.0.0d0) go to 210
      dbargs=0.5d0/xl*dbxws/args
      dbargd=0.5d0/xl*dbxwd/argd
      dasdb2=-xd*( d2sdb2+dsrgdb*dsrgdb*xd/(xj*argxs) )/(xl*args)
      daddb2=-xd*( d2bdb2+dbrgdb*dbrgdb*xd/(xj*argxd) )/(xl*argd)
      dgddb2=-0.5d0*gamma*(dasdb2+daddb2)
  210 dgddvb=-gamma*(dbargs+dbargd)
      dgsdvb=-2.0d0*gamma*dbargs
      if (xj.le.0.0d0) go to 220
      ddxwd=-dbxwd
      dgdvds=-gamma*0.5d0/xl*ddxwd/argd
      go to 220
  215 gamasd=gamma
      gamass=gamma
      gammad=gamma
      dgddvb=0.0d0
      dgsdvb=0.0d0
      dgdvds=0.0d0
      dgddb2=0.0d0
  220 von=vbin+gamasd*sarg
c     write(iofile,221) von,vbin,vbi,gamasd,argss,argsd,xj
  221 format ('0msg1:'/1p7d10.2)
      vth=von
      vdsat=0.0d0
  225 if (xnfs.eq.0.0d0.or.cox.eq.0.0d0) go to 230
      cfs=charge*xnfs
      cdonco=-(gamasd*dsrgdb+dgddvb*sarg)+factor
      xn=1.0d0+cfs/cox*xw*xl+cdonco
      von=von+vt*xn
c     write (iofile,226) von,cdonco,xn,cfs,xd
  226 format(' msg2:'/1p6d10.2)
      argg=1.0d0/(vt*xn)
      vgst=vgs-von
      go to 300
  230 vgst=vgs-von
      if (vgs.gt.von) go to 300
c
c  cutoff region
c
      gds=0.0d0
      go to 1050
c
c  compute some more useful quantities
c
  300 sarg3=sarg*sarg*sarg
      sbiarg=dsqrt(phib)
      gammad=gamasd
      dgdvbs=dgddvb
      body=barg*barg*barg-sarg3
      gdbdv=2.0d0*gammad*(barg*barg*dbrgdb-sarg*sarg*dsrgdb)
      dodvbs=-factor+dgdvbs*sarg+gammad*dsrgdb
      if (xnfs.eq.0.0d0) go to 400
      if (cox.eq.0.0d0) go to 410
      dxndvb=2.0d0*dgdvbs*dsrgdb+gammad*d2sdb2+dgddb2*sarg
      dodvbs=dodvbs+vt*dxndvb
      dxndvd=dgdvds*dsrgdb
      dodvds=dgdvds*sarg+vt*dxndvd
c
c  evaluate effective mobility and its derivatives
c
  400 if (cox.le.0.0d0) go to 410
      udenom=vgst
      if (udenom.le.vbp) go to 410
      ufact=dexp(uexp*dlog(vbp/udenom))
      ueff=uo*ufact
      dudvgs=-ufact*uexp/udenom
      dudvds=0.0d0
      dudvbs=uexp*ufact*dodvbs/vgst
      go to 500
  410 ufact=1.0d0
      ueff=uo
      dudvgs=0.0d0
      dudvds=0.0d0
      dudvbs=0.0d0
c
c     evaluate saturation voltage and its derivatives according to
c     grove-frohman equation
c
  500 vgsx=vgs
      gammad=gamasd/eta
      dgdvbs=dgddvb
      if (xnfs.ne.0.0d0.and.cox.ne.0.0d0)
     1   vgsx=dmax1(vgs,von)
  505 if (gammad.le.0.0d0) go to 535
      gammd2=gammad*gammad
      argv=(vgsx-vbin)/eta+phi-vbs
      if (argv.le.0.0d0) go to 540
      arg=dsqrt(1.0d0+4.0d0*argv/gammd2)
      vdsat=(vgsx-vbin)/eta+gammd2*(1.0d0-arg)/2.0d0
      vdsat=dmax1(vdsat,0.0d0)
  510 if (icharg.eq.0) go to 530
      arg1=gammd2/(eta*eta)
      arg2=vds-0.5d0*arg1
      argsq=(arg2+0.5d0*arg1+phi-vbs)*arg1
      if (argsq.ge.0.0d0) go to 515
      vpof=vth
      go to 520
  515 vpof=vbin+eta*(arg2+0.5d0*arg1+dsqrt(argsq))
  520 argv1=(vpof-vbin)/eta+phi-vbs
      if (argv1.gt.0.0d0) go to 525
      vdsat1=0.0d0
      go to 530
  525 arg1=dsqrt(1.0d0+4.0d0*argv1/gammd2)
      vdsat1=(vpof-vbin)/eta+gammd2*(1.0d0-arg1)/2.0d0
      vdsat1=dmax1(vdsat1,0.0d0)
  530 dsdvgs=(1.0d0-1.0d0/arg)/eta
      dsdvbs=(gammad*(1.0d0-arg)+2.0d0*argv/(gammad*arg))/eta*dgdvbs+
     1       1.0d0/arg+factor*dsdvgs
      go to 545
  535 vdsat=dmax1((vgsx-vbin)/eta,0.0d0)
      vpof=dmax1((eta*vds+vbin),0.0d0)
      vdsat1=dmax1((vpof-vbin)/eta,0.0d0)
      dsdvgs=1.0d0
      dsdvbs=0.0d0
      go to 545
  540 vdsat=0.0d0
      vpof=vth
      vdsat1=0.0d0
      dsdvgs=0.0d0
      dsdvbs=0.0d0
c
c     store vdsat as above in vpof (pinch-off)
c
  545 if (vmax.le.0.0d0) go to 600
c
c     evaluate saturation voltage and its derivatives according to
c     baum's theory of scattering velocity saturation
c
      gammd2=gammad*gammad
      v1=(vgsx-vbin)/eta+phi-vbs
      v2=phi-vbs
      xv=vmax*xl/ueff
      a1=gammad/0.75d0
      b1=-2.0d0*(v1+xv)
      c1=-2.0d0*gammad*xv
      d1=2.0d0*v1*(v2+xv)-v2*v2-4.0d0/3.0d0*gammad*sarg3
      a=-b1
      b=a1*c1-4.0d0*d1
      c=-d1*(a1*a1-4.0d0*b1)-c1*c1
      r=-a*a/3.0d0+b
      s=2.0d0*a*a*a/27.0d0-a*b/3.0d0+c
      r3=r*r*r
      s2=s*s
      p=s2/4.0d0+r3/27.0d0
      p0=dabs(p)
      p2=dsqrt(p0)
      if (p.ge.0.0d0) go to 550
      ro=dsqrt(s2/4.0d0+p0)
      ro=dlog(ro)/3.0d0
      ro=dexp(ro)
      fi=datan(-2.0d0*p2/s)
      y3=2.0d0*ro*dcos(fi/3.0d0)-a/3.0d0
      go to 560
  550 p3=dexp(dlog(dabs(-s/2.0d0+p2))/3.0d0)
      p4=dexp(dlog(dabs(-s/2.0d0-p2))/3.0d0)
      y3=p3+p4-a/3.0d0
  560 iknt=0
      a3=dsqrt(a1*a1/4.0d0-b1+y3)
      b3=dsqrt(y3*y3/4.0d0-d1)
      do 570 i=1,4
      a4(i)=a1/2.0d0+sig1(i)*a3
      b4(i)=y3/2.0d0+sig2(i)*b3
      delta4=a4(i)*a4(i)/4.0d0-b4(i)
      if (delta4.lt.0.0d0) go to 570
      iknt=iknt+1
      x4(iknt)=-a4(i)/2.0d0+dsqrt(delta4)
      iknt=iknt+1
      x4(iknt)=-a4(i)/2.0d0-dsqrt(delta4)
  570 continue
      jknt=0
      do 580 j=1,iknt
      if (x4(j).le.0.0d0) go to 580
      poly4(j)=x4(j)*x4(j)*x4(j)*x4(j)+a1*x4(j)*x4(j)*x4(j)
      poly4(j)=poly4(j)+b1*x4(j)*x4(j)+c1*x4(j)+d1
      if (dabs(poly4(j)).gt.1.0d-6) go to 580
      jknt=jknt+1
      if (jknt.gt.1) go to 575
      xvalid=x4(j)
  575 if (x4(j).gt.xvalid) go to 580
      xvalid=x4(j)
  580 continue
      if (jknt.gt.0) go to 590
      ivmflg=ivmflg+1
      go to 600
  590 vdsat=xvalid*xvalid+vbs-phi
c
c  evaluate effective channel length and its derivatives
c
  600 if (vds.eq.0.0d0) go to 610
      gammad=gamasd
      if ((vbs-vdsat).gt.0.0d0) go to 601
      bsarg=dsqrt(vdsat-vbs+phi)
      dbsrdb=-0.5d0/bsarg
      go to 602
  601 bsarg=sphi/(1.0d0+0.5d0*(vbs-vdsat)/phi)
      dbsrdb=-0.5d0*bsarg*bsarg/sphi3
  602 bodys=bsarg*bsarg*bsarg-sarg3
      gdbdvs=2.0d0*gammad*(bsarg*bsarg*dbsrdb-sarg*sarg*dsrgdb)
      if (vmax.gt.0.0d0) go to 603
      if (xnsub.eq.0.0d0) go to 610
      if (xlamda.gt.0.0d0) go to 610
      argv=(vds-vdsat)/4.0d0
      sargv=dsqrt(1.0d0+argv*argv)
      arg=dsqrt(argv+sargv)
      xlfact=xd/(xl*vds)
      xlamda=xlfact*arg
      dldsat=vds*xlfact*arg/(8.0d0*sargv)
      go to 605
  603 argv=(vgsx-vbin)/eta-vdsat
      xdv=xd/dsqrt(xneff)
      xlv=vmax*xdv/(2.0d0*ueff)
      vqchan=argv-gammad*bsarg
      dqdsat=-1.0d0+gammad*dbsrdb
      vl=vmax*xl
      dfunds=vl*dqdsat-ueff*vqchan
      dfundg=(vl-ueff*vdsat)/eta
      dfundb=-vl*(1.0d0+dqdsat-factor/eta)+
     1        ueff*(gdbdvs-dgdvbs*bodys/1.5d0)/eta
      dsdvgs=-dfundg/dfunds
      dsdvbs=-dfundb/dfunds
      if (xnsub.eq.0.0d0) go to 610
      if (xlamda.gt.0.0d0) go to 610
      argv=dmax1(vds-vdsat,0.0d0)
      xls=dsqrt(xlv*xlv+argv)
      dldsat=xdv/(2.0d0*xls)
      xlfact=xdv/(xl*vds)
      xlamda=xlfact*(xls-xlv)
      dldsat=dldsat/xl
  605 dldvgs=dldsat*dsdvgs
      dldvds=-xlamda+dldsat
      dldvbs=dldsat*dsdvbs
      go to 620
  610 dldvgs=0.0d0
      dldvds=0.0d0
      dldvbs=0.0d0
c
c     limit channel shortening at punch-through
c
  620 xwb=xd*sbiarg
      xld=xl-xwb
      clfact=1.0d0-xlamda*vds
      dldvds=-xlamda-dldvds
      xleff=xl*clfact
      deltal=xlamda*vds*xl
      if (xnsub.eq.0.0d0) xwb=0.25d-6
      if (xleff.ge.xwb) go to 700
      xleff=xwb/(1.0d0+(deltal-xld)/xwb)
      clfact=xleff/xl
      dfact=xleff*xleff/(xwb*xwb)
      dldvgs=dfact*dldvgs
      dldvds=dfact*dldvds
      dldvbs=dfact*dldvbs
c
c  evaluate effective beta (effective kp)
c
  700 beta1=beta*ufact/clfact
c
c  test for mode of operation and branch appropriately
c
      gammad=gamasd
      dgdvbs=dgddvb
      if (vds.gt.1.0d-8) go to 730
      if (vgs.gt.von) go to 720
      if ((xnfs.ne.0.0d0).and.(cox.ne.0.0d0)) go to 710
      gds=0.0d0
      go to 1050
c
  710 gds=beta1*(von-vbin-gammad*sarg)*dexp(argg*(vgs-von))
      go to 1050
c
c
  720 gds=beta1*(vgs-vbin-gammad*sarg)
      go to 1050
c
  730 if (vgs.gt.von) go to 900
c
c  subthreshold region
c
      if (vdsat.gt.0.0d0) go to 830
      gds=0.0d0
      if (vgs.gt.vth) go to 1020
      go to 1050
  830 vdson=dmin1(vdsat,vds)
      if (vds.le.vdsat) go to 850
      barg=bsarg
      dbrgdb=dbsrdb
      body=bodys
      gdbdv=gdbdvs
  850 cdson=beta1*((von-vbin-eta*vdson*0.5d0)*vdson-gammad*body/1.5d0)
      didvds=beta1*(von-vbin-eta*vdson-gammad*barg)
      gdson=-cdson*dldvds/clfact-beta1*dgdvds*body/1.5d0
      if (vds.lt.vdsat) gdson=gdson+didvds
      gbson=-cdson*dldvbs/clfact
     1      +beta1*(dodvbs*vdson+factor*vdson-dgdvbs*body/1.5d0-gdbdv)
      if (vds.gt.vdsat) gbson=gbson+didvds*dsdvbs
      expg=dexp(argg*(vgs-von))
      cdrain=cdson*expg
      gmw=cdrain*argg
      gm=gmw
      if (vds.gt.vdsat) gm=gmw+didvds*dsdvgs*expg
      gds=gdson*expg-gm*dodvds-gmw*(vgs-von)*dxndvd/xn
      gmbs=gbson*expg-gm*dodvbs-gmw*(vgs-von)*dxndvb/xn
      go to 1020
c
c
  900 if (vds.gt.vdsat) go to 1000
c
c  linear region
c
      cdrain=beta1*((vgs-vbin-eta*vds/2.0d0)*vds-gammad*body/1.5d0)
      arg=cdrain*(dudvgs/ufact-dldvgs/clfact)
      gm=arg+beta1*vds
      arg=cdrain*(dudvds/ufact-dldvds/clfact)
      gds=arg+beta1*(vgs-vbin-eta*vds-
     1   gammad*barg-dgdvds*body/1.5d0)
      arg=cdrain*(dudvbs/ufact-dldvbs/clfact)
      gmbs=arg-beta1*(gdbdv+dgdvbs*body/1.5d0-factor*vds)
      go to 1020
c
c  saturation region
c
 1000 cdrain=beta1*((vgs-vbin-eta*vdsat/2.0d0)*vdsat-gammad*bodys/1.5d0)
      arg=cdrain*(dudvgs/ufact-dldvgs/clfact)
      gm=arg+beta1*vdsat+
     1   beta1*(vgs-vbin-eta*vdsat-gammad*bsarg)*dsdvgs
      gds=-cdrain*dldvds/clfact-beta1*dgdvds*bodys/1.5d0
      arg=cdrain*(dudvbs/ufact-dldvbs/clfact)
      gmbs=arg-beta1*(gdbdvs+dgdvbs*bodys/1.5d0-factor*vdsat)+
     1     beta1*(vgs-vbin-eta*vdsat-gammad*bsarg)*dsdvbs
c
c     compute charges for "on" region
c
 1020 if (icharg.eq.0) go to 1500
      if (vgs.le.vth) go to 1070
      call mqspof(vds,vbs,vgs,vpof,vdsat1,vth,vbin,gamasd,
     1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      go to 2000
c
c  finish special cases
c
 1050 cdrain=0.0d0
      gm=0.0d0
      gmbs=0.0d0
 1070 xqc=xqco
      if (icharg.eq.0) go to 1500
      call mosq2(vds,vbs,vgs,vdsat,vth,vbin,gamasd,cox,phi,
     1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      qspof=0.0d0
      go to 2000
c
c  finished
c
 1500 qg=0.0d0
      qb=0.0d0
      qc=0.0d0
      qspof=0.0d0
 2000 return
      end
      subroutine moseq3(vds,vbs,vgs,gm,gds,gmbs,
     1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      implicit double precision (a-h,o-z)
c
c     this routine evaluates the drain current, its derivatives and
c     the charges associated with the gate, channel and bulk
c     for mosfets based on semi-empirical equations
c
c spice version 2g.6  sccsid=mosarg 3/15/83
      common /mosarg/ vto,beta,gamma,phi,phib,cox,xnsub,xnfs,xd,xj,xld,
     1   xlamda,uo,uexp,vbp,utra,vmax,xneff,xl,xw,vbi,von,vdsat,qspof,
     2   beta0,beta1,cdrain,xqco,xqc,fnarrw,fshort,lev
c spice version 2g.6  sccsid=status 3/15/83
      common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
     1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
     2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
c spice version 2g.6  sccsid=knstnt 3/15/83
      common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
     1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox,
     2   pivtol,pivrel
c
      equivalence (xlamda,alpha),(vbp,theta),(uexp,eta),(utra,xkappa)
      data coeff0/0.0631353d0/,coeff1/0.8013292d0/,coeff2/-0.01110777d0/
c
c     icharg=1 causes charges to be computed
c     icharg=0 bypasses the computation of charges
c
c     icharg=1
c     if (mode.ne.1) go to 10
c     icharg=0
c     if (modedc.eq.2.and.nosolv.ne.0) icharg=1
c     if (initf.eq.4) icharg=1
c
c     reference cdrain equations to source and
c     charge equations to bulk
c
10    continue
      icharg=0
      vgb=vgs-vbs
      vfb=vbi-phi
      vdsat=0.0d0
      qg=0.0d0
      qb=0.0d0
      qc=0.0d0
      cgdb=0.0d0
      cbdb=0.0d0
      onxl=1.0d0/xl
      eta=eta/(xl*xl*xl)
c
c.....square root term
c
      if ( vbs.gt.0.0d0 ) go to 120
      phibs=phi-vbs
      sqphbs=dsqrt(phibs)
      dsqdvb=-0.5d0/sqphbs
      go to 200
120   continue
      sqphis=dsqrt(phi)
      sqphs3=phi*sqphis
      sqphbs=sqphis/(1.0d0+vbs/(phi+phi))
      phibs=sqphbs*sqphbs
      dsqdvb=-phibs/(sqphs3+sqphs3)
c
c.....short channel effect factor
c
200   continue
      if ( (xj.eq.0.0d0).or.(xd.eq.0.0d0) ) go to 210
      wps=xd*sqphbs
      onxj=1.0d0/xj
      xjonxl=xj*onxl
      djonxj=xld*onxj
      wponxj=wps*onxj
      wconxj=coeff0+coeff1*wponxj+coeff2*wponxj*wponxj
      wcs=wconxj*xj
      arga=wconxj+djonxj
      argc=wponxj/(1.0d0+wponxj)
      argb=dsqrt(1.0d0-argc*argc)
      fshort=1.0d0-xjonxl*(arga*argb-djonxj)
      dwpdvb=xd*dsqdvb
      dadvb=(coeff1+coeff2*(wponxj+wponxj))*dwpdvb*onxj
      dbdvb=-argc*argc*(1.0d0-argc)*dwpdvb/(argb*wps)
      dfsdvb=-xjonxl*(dadvb*argb+arga*dbdvb)
      go to 220
210   continue
      fshort=1.0d0
      dfsdvb=0.0d0
      wcs=0.05d-6
c
c.....body effect
c
220   continue
      gammas=gamma*fshort
      fbodys=0.5d0*gammas/(sqphbs+sqphbs)
      fbody=fbodys+fnarrw
      onfbdy=1.0d0/(1.0d0+fbody)
      dfbdvb=-fbodys*dsqdvb/sqphbs+fbodys*dfsdvb/fshort
      qbonco=gammas*sqphbs+fnarrw*phibs
      dqbdvb=gammas*dsqdvb+gamma*dfsdvb*sqphbs-fnarrw
c
c.....static feedback effect
c
      vbix=vbi-eta*vds
c
c.....threshold voltage
c
      vth=vbix+qbonco
      dvtdvd=-eta
      dvtdvb=dqbdvb
c
c.....joint weak inversion and strong inversion
c
      von=vth
      if ( xnfs.eq.0.0d0 ) go to 250
           csonco=charge*xnfs*xl*xw/cox
           cdonco=qbonco/(phibs+phibs)
           xn=1.0d0+csonco+cdonco
           von=vth+vt*xn
           dxndvb=dqbdvb/(phibs+phibs)-qbonco*dsqdvb/(phibs*sqphbs)
           dvodvd=dvtdvd
           dvodvb=dvtdvb+vt*dxndvb
           go to 300
c
c.....cutoff region
c
250   continue
      if ( vgs.gt.von ) go to 300
      cdrain=0.0d0
      gm=0.0d0
      gds=0.0d0
      gmbs=0.0d0
      if ( icharg.ne.0 ) go to 800
      go to 1000
c
c.....device is on
c
300   continue
      vgsx=dmax1(vgs,von)
c
c.....mobility modulation by gate voltage
c
      onfg=1.0d0+theta*(vgsx-vth)
      fgate=1.0d0/onfg
      us=uo*fgate
      dfgdvg=-theta*fgate*fgate
      dfgdvd=-dfgdvg*dvtdvd
      dfgdvb=-dfgdvg*dvtdvb
c
c.....saturation voltage
c
      vdsat=(vgsx-vth)*onfbdy
      vpof=vdsat
      if ( vmax.gt.0.0d0 ) go to 310
      dvsdvg=onfbdy
      dvsdvd=-dvsdvg*dvtdvd
      dvsdvb=-dvsdvg*dvtdvb-vdsat*dfbdvb*onfbdy
      go to 400
  310 vdsc=xl*vmax/us
      onvdsc=1.0d0/vdsc
      arga=(vgsx-vth)*onfbdy
      argb=dsqrt(arga*arga+vdsc*vdsc)
      vdsat=arga+vdsc-argb
      dvsdga=(1.0d0-arga/argb)*onfbdy
      dvsdvg=dvsdga-(1.0d0-vdsc/argb)*vdsc*dfgdvg*onfg
      dvsdvd=-dvsdvg*dvtdvd
      dvsdvb=-dvsdvg*dvtdvb-arga*dvsdga*dfbdvb
c
c.....current factors in linear region
c
400   continue
      vdsx=dmin1(vds,vdsat)
      if ( vdsx.eq.0.0d0 ) go to 900
      cdo=vgsx-vth-0.5d0*(1.0d0+fbody)*vdsx
      dcodvg=1.0d0
      if (vds.lt.vdsat) dcodvd=-dvtdvd-0.5d0*(1.0d0+fbody)
      dcodvb=-dvtdvb-0.5d0*dfbdvb*vdsx
c
c.....normalized drain current
c
410   continue
      cdnorm=cdo*vdsx
      gm=vdsx
      gds=vgsx-vth-(1.0d0+fbody+dvtdvd)*vdsx
      gmbs=dcodvb*vdsx
c
c.....drain current without velocity saturation effect
c
      cd1=beta*cdnorm
      beta=beta*fgate
      cdrain=beta*cdnorm
      gm=beta*gm+dfgdvg*cd1
      gds=beta*gds+dfgdvd*cd1
      gmbs=beta*gmbs
c
c.....velocity saturation factor
c
      if ( vmax.eq.0.0d0 ) go to 500
      fdrain=1.0d0/(1.0d0+vdsx*onvdsc)
      fd2=fdrain*fdrain
      arga=fd2*vdsx*onvdsc*onfg
      dfddvg=-dfgdvg*arga
      dfddvd=-dfgdvd*arga-fd2*onvdsc
      dfddvb=-dfgdvb*arga
c
c.....drain current
c
      gm=fdrain*gm+dfddvg*cdrain
      gds=fdrain*gds+dfddvd*cdrain
      gmbs=fdrain*gmbs+dfddvb*cdrain
      cdrain=fdrain*cdrain
      beta=beta*fdrain
c
c.....channel length modulation
c
500   continue
      if ( vds.le.vdsat ) go to 700
      if ( vmax.eq.0.0d0 ) go to 510
      if (alpha.eq.0.0d0) go to 700
      cdsat=cdrain
      gdsat=cdsat*(1.0d0-fdrain)*onvdsc
      gdsat=dmax1(1.0d-12,gdsat)
      gdoncd=gdsat/cdsat
      gdonfd=gdsat/(1.0d0-fdrain)
      gdonfg=gdsat*onfg
      dgdvg=gdoncd*gm-gdonfd*dfddvg+gdonfg*dfgdvg
      dgdvd=gdoncd*gds-gdonfd*dfddvd+gdonfg*dfgdvd
      dgdvb=gdoncd*gmbs-gdonfd*dfddvb+gdonfg*dfgdvb
c
      emax=cdsat*onxl/gdsat
      emoncd=emax/cdsat
      emongd=emax/gdsat
      demdvg=emoncd*gm-emongd*dgdvg
      demdvd=emoncd*gds-emongd*dgdvd
      demdvb=emoncd*gmbs-emongd*dgdvb
c
      arga=0.5d0*emax*alpha
      argc=xkappa*alpha
      argb=dsqrt(arga*arga+argc*(vds-vdsat))
      delxl=argb-arga
      dldvd=argc/(argb+argb)
      dldem=0.5d0*(arga/argb-1.0d0)*alpha
      ddldvg=dldem*demdvg
      ddldvd=dldem*demdvd-dldvd
      ddldvb=dldem*demdvb
      go to 520
510   continue
      delxl=dsqrt(xkappa*(vds-vdsat)*alpha)
      dldvd=0.5d0*delxl/(vds-vdsat)
      ddldvg=0.0d0
      ddldvd=-dldvd
      ddldvb=0.0d0
c
c.....punch through approximation
c
520   continue
      if ( delxl.le.(0.5d0*xl) ) go to 600
      wcs2=wcs*wcs
      delxl=xl-(xl**2/(4.0d0*delxl))
      arga=4.0d0*(xl-delxl)**2/xl**2
      ddldvg=ddldvg*arga
      ddldvd=ddldvd*arga
      ddldvb=ddldvb*arga
       dldvd= dldvd*arga
c
c.....saturation region
c
600   continue
      dlonxl=delxl*onxl
      xlfact=1.0d0/(1.0d0-dlonxl)
      cdrain=cdrain*xlfact
      diddl=cdrain/(xl-delxl)
      gm=gm*xlfact+diddl*ddldvg
      gds0=gds*xlfact+diddl*ddldvd
      gmbs=gmbs*xlfact+diddl*ddldvb
      gm=gm+gds0*dvsdvg
      gmbs=gmbs+gds0*dvsdvb
      gds=gds0*dvsdvd+diddl*dldvd
c
c.....finish strong inversion case
c
700   continue
      if ( vgs.ge.von ) go to 750
c
c.....weak inversion
c
                onxn=1.0d0/xn
                ondvt=onxn/vt
                wfact=dexp( (vgs-von)*ondvt )
                cdrain=cdrain*wfact
                gms=gm*wfact
                gmw=cdrain*ondvt
                gm=gmw
                if (vds.gt.vdsat) gm=gm+gds0*dvsdvg*wfact
                gds=gds*wfact+(gms-gmw)*dvodvd
                gmbs=gmbs*wfact+(gms-gmw)*dvodvb
     1                         -gmw*(vgs-von)*onxn*dxndvb
c
c.....charge computation
c
  750 continue
      if (icharg.eq.0) go to 1000
      if (vgs.le.vth) go to 800
      call mqspof(vds,vbs,vgs,vpof,vdsat1,vth,vbin,gamasd,
     1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      go to 2000
c
c.....charge computation for vgs<vth
c
800   continue
      xqc=xqco
      call mosq3(vds,vbs,vpof,vdsat1,vth,vbin,gamasd,cox,phi,
     1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      qspof=0.0d0
      go to 2000
c
c.....special case of vds=0.0d0
c
900   continue
      beta=beta*fgate
      cdrain=0.0d0
      gm=0.0d0
      gds=beta*(vgsx-vth)
      gmbs=0.0d0
           if ( (xnfs.ne.0.0d0).and.(vgs.lt.von) )
     1          gds=gds*dexp((vgs-von)/(vt*xn))
      if (icharg.eq.0) go to 1000
      call mosq3(vds,vbs,vpof,vdsat1,vth,vbin,gamasd,cox,phi,
     1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
1000  qspof=0.0d0
c
c.....done
c
 2000 return
      end
      subroutine mosfet
      implicit double precision (a-h,o-z)
c
c     this routine processes mosfets for dc and transient analyses.
c
c spice version 2g.6  sccsid=tabinf 3/15/83
      common /tabinf/ ielmnt,isbckt,nsbckt,iunsat,nunsat,itemps,numtem,
     1   isens,nsens,ifour,nfour,ifield,icode,idelim,icolum,insize,
     2   junode,lsbkpt,numbkp,iorder,jmnode,iur,iuc,ilc,ilr,numoff,isr,
     3   nmoffc,iseq,iseq1,neqn,nodevs,ndiag,iswap,iequa,macins,lvnim1,
     4   lx0,lvn,lynl,lyu,lyl,lx1,lx2,lx3,lx4,lx5,lx6,lx7,ld0,ld1,ltd,
     5   imynl,imvn,lcvn,nsnod,nsmat,nsval,icnod,icmat,icval,
     6   loutpt,lpol,lzer,irswpf,irswpr,icswpf,icswpr,irpt,jcpt,
     7   irowno,jcolno,nttbr,nttar,lvntmp
c spice version 2g.6  sccsid=cirdat 3/15/83
      common /cirdat/ locate(50),jelcnt(50),nunods,ncnods,numnod,nstop,
     1   nut,nlt,nxtrm,ndist,ntlin,ibr,numvs,numalt,numcyc
c spice version 2g.6  sccsid=mosarg 3/15/83
      common /mosarg/ vto,beta,gamma,phi,phib,cox,xnsub,xnfs,xd,xj,xld,
     1   xlamda,uo,uexp,vbp,utra,vmax,xneff,xl,xw,vbi,von,vdsat,qspof,
     2   beta0,beta1,cdrain,xqco,xqc,fnarrw,fshort,lev
c spice version 2g.6  sccsid=status 3/15/83
      common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
     1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
     2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
c spice version 2g.6  sccsid=knstnt 3/15/83
      common /knstnt/ twopi,xlog2,xlog10,root2,rad,boltz,charge,ctok,
     1   gmin,reltol,abstol,vntol,trtol,chgtol,eps0,epssil,epsox,
     2   pivtol,pivrel
c spice version 2g.6  sccsid=debug 3/15/83
      common/debug/ idebug(20)
c spice version 2g.6  sccsid=blank 3/15/83
      common /blank/ value(50000)
      integer nodplc(64)
      complex cvalue(32)
      equivalence (value(1),nodplc(1),cvalue(1))
c
c
      dimension vbdo(1),vbso(1),vgso(1),vdso(1),cdo(1),cbso(1),cbdo(1),
     1   gmo(1),gdso(1),gmbso(1),gbdo(1),gbso(1),
     2   qb(1),cqb(1),qg(1),cqg(1),qd(1),cqd(1),
     3   cggbo(1),cgdbo(1),cgsbo(1),cbgbo(1),cbdbo(1),cbsbo(1),
     4   cgbo(1),cgdo(1),cgso(1),vono(1),vdsato(1)
      dimension qbd(1),cqbd(1),qbs(1),cqbs(1),
     1          qgs(1),cqgs(1),qgd(1),cqgd(1),qgb(1),cqgb(1)
      equivalence (vbdo (1),value( 1)),(vbso (1),value( 2)),
     1            (vgso (1),value( 3)),(vdso (1),value( 4)),
     2            (cdo  (1),value( 5)),(cbso (1),value( 6)),
     3            (cbdo (1),value( 7)),(gmo  (1),value( 8)),
     4            (gdso (1),value( 9)),(gmbso(1),value(10)),
     5            (gbdo (1),value(11)),(gbso (1),value(12)),
     6            (qb   (1),qgs  ( 1), value(13)),
     7            (cqb  (1),cqgs ( 1), value(14)),
     8            (qg   (1),qgd  ( 1), value(15)),
     9            (cqg  (1),cqgd ( 1), value(16)),
     a            (qd   (1),qgb  ( 1), value(17)),
     b            (cqd  (1),cqgb ( 1), value(18)),
     c            (cggbo(1),cgbo  (1), value(19)),
     d            (cgdbo(1),cgdo  (1), value(20)),
     e            (cgsbo(1),cgso  (1), value(21)),
     f            (cbgbo(1),vono  (1), value(22)),
     g            (cbdbo(1),vdsato(1), value(23)),
     h            (cbsbo(1),           value(24))
      equivalence (qbd  (1),value(25)),(cqbd (1),value(26)),
     1            (qbs  (1),value(27)),(cqbs (1),value(28))
c
c
      loc=locate(14)
   10 if ((loc.eq.0).or.(nodplc(loc+33).ne.0)) return
      locv=nodplc(loc+1)
      node1=nodplc(loc+2)
      node2=nodplc(loc+3)
      node3=nodplc(loc+4)
      node4=nodplc(loc+5)
      node5=nodplc(loc+6)
      node6=nodplc(loc+7)
      locm=nodplc(loc+8)
      ioff=nodplc(loc+9)
      type=nodplc(locm+2)
      locm=nodplc(locm+1)
      loct=nodplc(loc+26)
c
c  dc model parameters
c
      xj=value(locm+27)
      xld=value(locm+28)
      xl=value(locv+1)-2.0d0*xld
      xw=value(locv+2)
      devmod=value(locv+8)
      vto=type*value(locm+2)
      vdsat=type*value(locv+10)
      vinit=value(locm+43)
      ad=value(locv+3)
      as=value(locv+4)
      pd=value(locv+11)
      ps=value(locv+12)
      if (value(locm+21).eq.0.0d0.
     1   or.ad.eq.0.0d0.or.as.eq.0.0d0) go to 12
      cdsat=value(locm+21)*ad
      cssat=value(locm+21)*as
      go to 15
   12 cdsat=value(locm+11)
      cssat=value(locm+11)
   15 if ((value(locm+7).le.0.0d0).and.
     1              (value(locm+8).le.0.0d0)) go to 17
      gdpr=value(locm+7)
      gspr=value(locm+8)
      go to 19
   17 gdpr=value(locm+16)/value(locv+13)
      gspr=value(locm+16)/value(locv+14)
   19 covlgs=value(locm+13)*xw
      covlgd=value(locm+14)*xw
      covlgb=value(locm+15)*xl
      lev=value(locm+1)
c
c     mos1, mos2 and mos3 model parameters
c
      beta=value(locm+3)*xw/xl
      gamma=value(locm+4)
      phi=value(locm+5)
      xlamda=value(locm+6)
      phib=value(locm+12)
      cox=value(locm+22)*xw*xl
      xnsub=value(locm+23)
      xnfs=value(locm+25)
      uo=value(locm+29)
      vbp=value(locm+30)
      uexp=value(locm+31)
      utra=value(locm+32)
      vbi=type*value(locm+44)
      xd=value(locm+45)
      vmax=value(locm+33)
      xneff=value(locm+34)
      xqco=value(locm+35)
      fnarrw=value(locm+39)
      if (lev.eq.3) fnarrw=fnarrw/xw
c
c     initialization
c
      icheck=1
      ibypas=0
      go to (100,20,30,50,60,70), initf
   20 if (ioff.ne.0) go to 40
      vds=type*value(locv+5)
      vgs=type*value(locv+6)
      vbs=type*value(locv+7)
      if (vds.ne.0.0d0) go to 300
      if (vgs.ne.0.0d0) go to 300
      if (vbs.ne.0.0d0) go to 300
      if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 300
      vbs=vinit
      vgs=vto
      vds=0.0d0
      go to 300
   30 if (ioff.eq.0) go to 100
   40 vbs=0.0d0
      vgs=0.0d0
      vds=0.0d0
      go to 300
   50 vbs=vbso(lx0+loct)
      vgs=vgso(lx0+loct)
      vds=vdso(lx0+loct)
      go to 300
   60 vbs=vbso(lx1+loct)
      vgs=vgso(lx1+loct)
      vds=vdso(lx1+loct)
      go to 300
   70 xfact=delta/delold(2)
      vbso(lx0+loct)=vbso(lx1+loct)
      vbs=(1.0d0+xfact)*vbso(lx1+loct)-xfact*vbso(lx2+loct)
      vgso(lx0+loct)=vgso(lx1+loct)
      vgs=(1.0d0+xfact)*vgso(lx1+loct)-xfact*vgso(lx2+loct)
      vdso(lx0+loct)=vdso(lx1+loct)
      vds=(1.0d0+xfact)*vdso(lx1+loct)-xfact*vdso(lx2+loct)
      vbdo(lx0+loct)=vbso(lx0+loct)-vdso(lx0+loct)
      cdo(lx0+loct)=cdo(lx1+loct)
      cbso(lx0+loct)=cbso(lx1+loct)
      cbdo(lx0+loct)=cbdo(lx1+loct)
      gmo(lx0+loct)=gmo(lx1+loct)
      gdso(lx0+loct)=gdso(lx1+loct)
      gmbso(lx0+loct)=gmbso(lx1+loct)
      gbdo(lx0+loct)=gbdo(lx1+loct)
      gbso(lx0+loct)=gbso(lx1+loct)
      cggbo(lx0+loct)=cggbo(lx1+loct)
      cgdbo(lx0+loct)=cgdbo(lx1+loct)
      cgsbo(lx0+loct)=cgsbo(lx1+loct)
      cbgbo(lx0+loct)=cbgbo(lx1+loct)
      cbdbo(lx0+loct)=cbdbo(lx1+loct)
      cbsbo(lx0+loct)=cbsbo(lx1+loct)
      go to 110
c
c  compute new nonlinear branch voltages
c
  100 vbs=type*(value(lvnim1+node4)-value(lvnim1+node6))
      vgs=type*(value(lvnim1+node2)-value(lvnim1+node6))
      vds=type*(value(lvnim1+node5)-value(lvnim1+node6))
  110 vbd=vbs-vds
      vgd=vgs-vds
      vgdo=vgso(lx0+loct)-vdso(lx0+loct)
      delvbs=vbs-vbso(lx0+loct)
      delvbd=vbd-vbdo(lx0+loct)
      delvgs=vgs-vgso(lx0+loct)
      delvds=vds-vdso(lx0+loct)
      delvgd=vgd-vgdo
      if (devmod.lt.0.0d0) go to 120
      cdhat=cdo(lx0+loct)-gbdo(lx0+loct)*delvbd+gmbso(lx0+loct)*delvbs
     1   +gmo(lx0+loct)*delvgs+gdso(lx0+loct)*delvds
      go to 130
  120 cdhat=cdo(lx0+loct)-(gbdo(lx0+loct)-gmbso(lx0+loct))*delvbd
     1   -gmo(lx0+loct)*delvgd+gdso(lx0+loct)*delvds
  130 cbhat=cbso(lx0+loct)+cbdo(lx0+loct)+gbdo(lx0+loct)*delvbd
     1   +gbso(lx0+loct)*delvbs
c
c  bypass if solution has not changed
c
      if (initf.eq.6) go to 200
      tol=reltol*dmax1(dabs(vbs),dabs(vbso(lx0+loct)))+vntol
      if (dabs(delvbs).ge.tol) go to 200
      tol=reltol*dmax1(dabs(vbd),dabs(vbdo(lx0+loct)))+vntol
      if (dabs(delvbd).ge.tol) go to 200
      tol=reltol*dmax1(dabs(vgs),dabs(vgso(lx0+loct)))+vntol
      if (dabs(delvgs).ge.tol) go to 200
      tol=reltol*dmax1(dabs(vds),dabs(vdso(lx0+loct)))+vntol
      if (dabs(delvds).ge.tol) go to 200
      tol=reltol*dmax1(dabs(cdhat),dabs(cdo(lx0+loct)))+abstol
      if (dabs(cdhat-cdo(lx0+loct)).ge.tol) go to 200
      tol=reltol*dmax1(dabs(cbhat),dabs(cbso(lx0+loct)+cbdo(lx0+loct)))
     1   +abstol
      if (dabs(cbhat-(cbso(lx0+loct)+cbdo(lx0+loct))).ge.tol) go to 200
      vbs=vbso(lx0+loct)
      vbd=vbdo(lx0+loct)
      vgs=vgso(lx0+loct)
      vds=vdso(lx0+loct)
      vgd=vgs-vds
      vgb=vgs-vbs
      cd=cdo(lx0+loct)
      cbs=cbso(lx0+loct)
      cbd=cbdo(lx0+loct)
      cdrain=devmod*(cd+cbd)
      gm=gmo(lx0+loct)
      gds=gdso(lx0+loct)
      gmbs=gmbso(lx0+loct)
      gbd=gbdo(lx0+loct)
      gbs=gbso(lx0+loct)
      devmod=value(locv+8)
      if (mode.ne.1) go to 135
      if ((modedc.eq.2).and.(nosolv.ne.0)) go to 135
      if (xqco.gt.0.5d0) go to 742
      go to 850
  135 if (xqco.le.0.5d0) go to 140
      cgb=cgbo(lx0+loct)
      cgd=cgdo(lx0+loct)
      cgs=cgso(lx0+loct)
      vgs1=vgso(lx1+loct)
      vgb1=vgs1-vbso(lx1+loct)
      vgd1=vgs1-vdso(lx1+loct)
      go to 735
  140 cggb=cggbo(lx0+loct)
      cgdb=cgdbo(lx0+loct)
      cgsb=cgsbo(lx0+loct)
      cbgb=cbgbo(lx0+loct)
      cbdb=cbdbo(lx0+loct)
      cbsb=cbsbo(lx0+loct)
      xqc=value(locv+15)
      ibypas=1
      go to 755
c
c  limit nonlinear branch voltages
c
  200 von=type*value(locv+9)
      if (vdso(lx0+loct).lt.0.0d0) go to 205
      call fetlim(vgs,vgso(lx0+loct),von)
      vds=vgs-vgd
      call limvds(vds,vdso(lx0+loct))
      vgd=vgs-vds
      go to 210
  205 call fetlim(vgd,vgdo,von)
      vds=vgs-vgd
      call limvds(-vds,-vdso(lx0+loct))
      vgs=vgd+vds
  210 if (vds.lt.0.0d0) go to 220
      vcrit=vt*dlog(vt/(root2*cssat))
      call pnjlim(vbs,vbso(lx0+loct),vt,vcrit,icheck)
      vbd=vbs-vds
      go to 300
  220 vcrit=vt*dlog(vt/(root2*cdsat))
      call pnjlim(vbd,vbdo(lx0+loct),vt,vcrit,icheck)
      vbs=vbd+vds
c
c  determine dc current and derivatives
c
  300 vbd=vbs-vds
      vgd=vgs-vds
      vgb=vgs-vbs
      if (vbs.gt.0.0d0) go to 310
      gbs=cssat/vt
      cbs=gbs*vbs
      gbs=gbs+gmin
      go to 320
  310 evbs=dexp(vbs/vt)
      gbs=cssat*evbs/vt+gmin
      cbs=cssat*(evbs-1.0d0)
  320 if (vbd.gt.0.0d0) go to 330
      gbd=cdsat/vt
      cbd=gbd*vbd
      gbd=gbd+gmin
      go to 400
  330 evbd=dexp(vbd/vt)
      gbd=cdsat*evbd/vt+gmin
      cbd=cdsat*(evbd-1.0d0)
c
c  compute drain current and derivatives
c
  400 if (vds.lt.0.0d0) go to 450
c
c  normal mode
c
      devmod=1.0d0
      value(locv+8)=devmod
      go to (405,410,415), lev
  405 call moseq1(vds,vbs,vgs,gm,gds,gmbs)
      go to 460
  410 call moseq2(vds,vbs,vgs,gm,gds,gmbs,
     1   qgate,qchan,qbulk,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      go to 460
  415 call moseq3(vds,vbs,vgs,gm,gds,gmbs,
     1   qgate,qchan,qbulk,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      go to 460
c
c  inverse mode
c
  450 devmod=-1.0d0
      value(locv+8)=devmod
      go to (452,453,454), lev
  452 call moseq1(-vds,vbd,vgd,gm,gds,gmbs)
      go to 460
  453 call moseq2(-vds,vbd,vgd,gm,gds,gmbs,
     1   qgate,qchan,qbulk,cggb,cgsb,cgdb,cbgb,cbsb,cbdb)
      go to 460
  454 call moseq3(-vds,vbd,vgd,gm,gds,gmbs,
     1   qgate,qchan,qbulk,cggb,cgsb,cgdb,cbgb,cbsb,cbdb)
  460 value(locv+9)=type*von
      value(locv+10)=type*vdsat
      if (xqco.le.0.5d0) value(locv+15)=xqc
c
c  compute equivalent drain current source
c
  490 cd=devmod*cdrain-cbd
      if (mode.ne.1) go to 500
      if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500
      if (initf.eq.4) go to 500
      go to 650
c
c  charge storage elements
c
c.. bulk-drain and bulk-source depletion capacitances
c
  500 czbd=0.0d0
      czbs=0.0d0
      czbdsw=0.0d0
      czbssw=0.0d0
      if ((value(locm+9).eq.0.0d0).or.(value(locm+10).eq.0.0d0))
     1   go to 505
      czbd=value(locm+9)
      czbs=value(locm+10)
      go to 510
  505 if (value(locm+17).eq.0.0d0) go to 510
      czbd=value(locm+17)*ad
      czbs=value(locm+17)*as
  510 if (value(locm+19).eq.0.0d0) go to 515
      czbdsw=value(locm+19)*pd
      czbssw=value(locm+19)*ps
  515 phib=value(locm+12)
      xmj=value(locm+18)
      xmjsw=value(locm+20)
      twop=phib+phib
      fcpb=value(locm+38)
      fcpb2=fcpb*fcpb
      f1=value(locm+40)
      f2=value(locm+41)
      f3=value(locm+42)
      czsf2=czbs/f2
      czswf2=czbssw/f2
      czdf2=czbd/f2
      czdwf2=czbdsw/f2
      if (vbs.ge.fcpb) go to 520
      arg=1.0d0-vbs/phib
      sarg=dexp(-xmj*dlog(arg))
      sargsw=dexp(-xmjsw*dlog(arg))
      qbs(lx0+loct)=phib*(czbs*(1.0d0-arg*sarg)/(1.0d0-xmj)
     1                +czbssw*(1.0d0-arg*sargsw)/(1.0d0-xmjsw))
      capbs=czbs*sarg+czbssw*sargsw
      go to 525
  520 qbs(lx0+loct)=f1*(czbs+czbssw)+f3*(vbs-fcpb)*(czsf2+czswf2)
     1    +(vbs*vbs-fcpb*fcpb)*(czsf2*xmj+czswf2*xmjsw)
      capbs=f3*(czsf2+czswf2)+vbs/phib*(czsf2*xmj+czswf2*xmjsw)
  525 if (vbd.ge.fcpb) go to 530
      arg=1.0d0-vbd/phib
      sarg=dexp(-xmj*dlog(arg))
      sargsw=dexp(-xmjsw*dlog(arg))
      qbd(lx0+loct)=phib*(czbd*(1.0d0-arg*sarg)/(1.0d0-xmj)
     1              +czbdsw*(1.0d0-arg*sargsw)/(1.0d0-xmjsw))
      capbd=czbd*sarg+czbdsw*sargsw
      go to 560
  530 qbd(lx0+loct)=f1*(czbd+czbdsw)+f3*(vbd-fcpb)*(czdf2+czdwf2)
     1    +(vbd*vbd-fcpb*fcpb)*(czdf2*xmj+czdwf2*xmjsw)
      capbd=f3*(czdf2+czdwf2)+vbd/phib*(czdf2*xmj+czdwf2*xmjsw)
c
  560 if (xqco.le.0.5d0) go to 650
      if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 650
      if (initf.ne.4) go to 600
      go to 705
c
cc    calculate equivalent conductances and currents for
cc    depletion capacitors
c
  600 if (initf.ne.5) go to 610
      qbd(lx1+loct)=qbd(lx0+loct)
      qbs(lx1+loct)=qbs(lx0+loct)
  610 call intgr8(geq,ceq,capbd,loct+24)
      gbd=gbd+geq
      cbd=cbd+cqbd(lx0+loct)
      cd=cd-cqbd(lx0+loct)
      call intgr8(geq,ceq,capbs,loct+26)
      gbs=gbs+geq
      cbs=cbs+cqbs(lx0+loct)
      if (initf.ne.5) go to 650
      cqbd(lx1+loct)=cqbd(lx0+loct)
      cqbs(lx1+loct)=cqbs(lx0+loct)
c
c  check convergence
c
  650 if (initf.ne.3) go to 660
      if (ioff.ne.0) go to 680
  660 if (icheck.eq.1) go to 670
      tol=reltol*dmax1(dabs(cdhat),dabs(cd))+abstol
      if (dabs(cdhat-cd).ge.tol) go to 670
      tol=reltol*dmax1(dabs(cbhat),dabs(cbs+cbd))+abstol
      if (dabs(cbhat-(cbs+cbd)).le.tol) go to 680
  670 noncon=noncon+1
  680 vbso(lx0+loct)=vbs
      vbdo(lx0+loct)=vbd
      vgso(lx0+loct)=vgs
      vdso(lx0+loct)=vds
      cdo(lx0+loct)=cd
      cbso(lx0+loct)=cbs
      cbdo(lx0+loct)=cbd
      gmo(lx0+loct)=gm
      gdso(lx0+loct)=gds
      gmbso(lx0+loct)=gmbs
      gbdo(lx0+loct)=gbd
      gbso(lx0+loct)=gbs
      if (xqco.le.0.5d0) go to 690
      vono(lx0+loct)=von
      vdsato(lx0+loct)=vdsat
      go to 700
  690 cggbo(lx0+loct)=cggb
      cgdbo(lx0+loct)=cgdb
      cgsbo(lx0+loct)=cgsb
      cbgbo(lx0+loct)=cbgb
      cbdbo(lx0+loct)=cbdb
      cbsbo(lx0+loct)=cbsb
      go to 750
c
c     xqco > 0.5d0 use meyer"s capacitor model
c
  700 if (mode.ne.1) go to 705
      if ((modedc.eq.2).and.(nosolv.ne.0)) go to 705
      if (initf.eq.4) go to 705
      go to 742
c
c     calculate meyer's capacitors
c
  705 von1=von
      vgs1=vgs
      vgd1=vgd
      vgb1=vgs-vbs
      vdsat1=vdsat
      if ((mode.ne.2).or.(initf.eq.5)) go to 710
      von1=vono(lx1+loct)
      vgs1=vgso(lx1+loct)
      vgd1=vgs1-vdso(lx1+loct)
      vgb1=vgs1-vbso(lx1+loct)
      vdsat1=vdsato(lx1+loct)
  710 if (devmod.lt.0.0d0) go to 715
      call cmeyer (vgs1,vgd1,vgb1,von1,vdsat1,vgs,vgd,vgb,
     1   covlgs,covlgd,covlgb,cgs1,cgd1,cgb1,cgs,cgd,cgb)
      go to 720
  715 call cmeyer (vgd1,vgs1,vgb1,von1,vdsat1,vgd,vgs,vgb,
     1   covlgd,covlgs,covlgb,cgd1,cgs1,cgb1,cgd,cgs,cgb)
  720 cgs=0.5d0*(cgs+cgs1)
      cgd=0.5d0*(cgd+cgd1)
      cgb=0.5d0*(cgb+cgb1)
c
c     store small-signal parameters (for meyer"s model)
c
      if (mode.ne.1) go to 730
      if (initf.ne.4) go to 730
      value(lx0+loct+24)=capbd
      value(lx0+loct+26)=capbs
      value(lx0+loct+12)=cgs-covlgs
      value(lx0+loct+14)=cgd-covlgd
      value(lx0+loct+16)=cgb-covlgb
      go to 1000
cc
  730 if (initf.ne.6) go to 735
      qgs(lx0+loct)=(1.0d0+xfact)*qgs(lx1+loct)-xfact*qgs(lx2+loct)
      qgd(lx0+loct)=(1.0d0+xfact)*qgd(lx1+loct)-xfact*qgd(lx2+loct)
      qgb(lx0+loct)=(1.0d0+xfact)*qgb(lx1+loct)-xfact*qgb(lx2+loct)
      go to 745
  735 qgs(lx0+loct)=(vgs-vgs1)*cgs
      qgd(lx0+loct)=(vgd-vgd1)*cgd
      qgb(lx0+loct)=(vgb-vgb1)*cgb
      if((mode.ne.2).or.(initf.eq.5))go to 740
      qgs(lx0+loct)=qgs(lx0+loct)+qgs(lx1+loct)
      qgd(lx0+loct)=qgd(lx0+loct)+qgd(lx1+loct)
      qgb(lx0+loct)=qgb(lx0+loct)+qgb(lx1+loct)
  740 if((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 742
      if (initf.ne.5) go to 745
      qgs(lx0+loct)=cgs*vgs
      qgd(lx0+loct)=cgd*vgd
      qgb(lx0+loct)=cgb*vgb
      qgs(lx1+loct)=qgs(lx0+loct)
      qgd(lx1+loct)=qgd(lx0+loct)
      qgb(lx1+loct)=qgb(lx0+loct)
c
c     initialize to zero charge conductances and current
c
  742 gcgs=0.0d0
      ceqgs=0.0d0
      gcgd=0.0d0
      ceqgd=0.0d0
      gcgb=0.0d0
      ceqgb=0.0d0
      go to  870
cc
 745  if(cgs.eq.0.0d0) value(lx0+loct+13)=0.0d0
      if(cgd.eq.0.0d0) value(lx0+loct+15)=0.0d0
      if(cgb.eq.0.0d0) value(lx0+loct+17)=0.0d0
cc
cc    calculate equivalent conductances and currents for
cc    meyer"s capacitors
cc
      call intgr8(gcgs,ceqgs,cgs,loct+12)
      call intgr8(gcgd,ceqgd,cgd,loct+14)
      call intgr8(gcgb,ceqgb,cgb,loct+16)
      ceqgs=ceqgs-gcgs*vgs+ag(1)*qgs(lx0+loct)
      ceqgd=ceqgd-gcgd*vgd+ag(1)*qgd(lx0+loct)
      ceqgb=ceqgb-gcgb*vgb+ag(1)*qgb(lx0+loct)
      if (initf.ne.5) go to 870
      cqgs(lx1+loct)=cqgs(lx0+loct)
      cqgd(lx1+loct)=cqgd(lx0+loct)
      cqgb(lx1+loct)=cqgb(lx0+loct)
      go to 870
c
c.. bulk and channel charge (plus overlaps)
c
  750 if (mode.ne.1) go to 755
      if ((modedc.eq.2).and.(nosolv.ne.0)) go to 755
      if (initf.eq.4) go to 755
      go to 850
  755 if (devmod.eq.-1.0d0) go to 760
      call moscap(vgd,vgs,vgb,covlgd,covlgs,covlgb,
     1   capbd,capbs,cggb,cgdb,cgsb,cbgb,cbdb,cbsb,
     2   gcggb,gcgdb,gcgsb,gcbgb,gcbdb,gcbsb,
     3   gcdgb,gcddb,gcdsb,gcsgb,gcsdb,gcssb,
     4   qgate,qchan,qbulk,qdrn,qsrc)
      go to 780
  760 call moscap(vgs,vgd,vgb,covlgs,covlgd,covlgb,
     1   capbs,capbd,cggb,cgsb,cgdb,cbgb,cbsb,cbdb,
     2   gcggb,gcgsb,gcgdb,gcbgb,gcbsb,gcbdb,
     3   gcsgb,gcssb,gcsdb,gcdgb,gcdsb,gcddb,
     4   qgate,qchan,qbulk,qsrc,qdrn)
  780 if (ibypas.eq.1) go to 860
      qg(lx0+loct)=qgate
      qd(lx0+loct)=qdrn-qbd(lx0+loct)
      qb(lx0+loct)=qbulk+qbd(lx0+loct)+qbs(lx0+loct)
c
c  store small-signal parameters
c
  790 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 850
      if (initf.ne.4) go to 800
      value(lx0+loct+18)=cggb
      value(lx0+loct+19)=cgdb
      value(lx0+loct+20)=cgsb
      value(lx0+loct+21)=cbgb
      value(lx0+loct+22)=cbdb
      value(lx0+loct+23)=cbsb
      value(lx0+loct+24)=capbd
      value(lx0+loct+26)=capbs
      go to 1000
c
c  transient analysis
c
  800 if (initf.ne.5) go to 810
      qb(lx1+loct)=qb(lx0+loct)
      qg(lx1+loct)=qg(lx0+loct)
      qd(lx1+loct)=qd(lx0+loct)
c.. integrate qb
  810 call intgr8(geq,ceq,0.0d0,loct+12)
c.. integrate qg
      call intgr8(geq,ceq,0.0d0,loct+14)
c.. integrate qd
      call intgr8(geq,ceq,0.0d0,loct+16)
      go to 860
c
c     initialize to zero charge conductances and current
c
  850 ceqqg=0.0d0
      ceqqb=0.0d0
      ceqqd=0.0d0
      gcdgb=0.0d0
      gcddb=0.0d0
      gcdsb=0.0d0
      gcsgb=0.0d0
      gcsdb=0.0d0
      gcssb=0.0d0
      gcggb=0.0d0
      gcgdb=0.0d0
      gcgsb=0.0d0
      gcbgb=0.0d0
      gcbdb=0.0d0
      gcbsb=0.0d0
      go to 900
c
c     evaluate equivalent charge currents
c
  860 cgate=cqg(lx0+loct)
      cqbulk=cqb(lx0+loct)
      cqdrn=cqd(lx0+loct)
      ceqqg=cgate-gcggb*vgb+gcgdb*vbd+gcgsb*vbs
      ceqqb=cqbulk-gcbgb*vgb+gcbdb*vbd+gcbsb*vbs
      ceqqd=cqdrn-gcdgb*vgb+gcddb*vbd+gcdsb*vbs
      if (initf.ne.5) go to 900
      cqb(lx1+loct)=cqb(lx0+loct)
      cqg(lx1+loct)=cqg(lx0+loct)
      cqd(lx1+loct)=cqd(lx0+loct)
      go to 900
c
cc   do the mapping from meyer"s capacitor model into the charge
cc   oriented model
cc
  870 ceqqg=ceqgs+ceqgb+ceqgd
      ceqqb=-ceqgb
      ceqqd=-ceqgd
      gcbdb=0.0d0
      gcbsb=0.0d0
      gcdsb=0.0d0
      gcsdb=0.0d0
      gcgdb=-gcgd
      gcgsb=-gcgs
      gcbgb=-gcgb
      gcdgb=-gcgd
      gcsgb=-gcgs
      gcssb=gcgs
      gcddb=gcgd
      gcggb=gcgd+gcgs+gcgb
c
c     store charge storage info for meyer's cap in lx table
c
      cgbo(lx0+loct)=cgb
      cgso(lx0+loct)=cgs
      cgdo(lx0+loct)=cgd
c
c  load current vector
c
  900 ceqbs=type*(cbs-(gbs-gmin)*vbs)
      ceqbd=type*(cbd-(gbd-gmin)*vbd)
      ceqqg=type*ceqqg
      ceqqb=type*ceqqb
      ceqqd=type*ceqqd
      xnrm=1.0d0
      xrev=0.0d0
      if (devmod.lt.0.0d0) go to 910
      cdreq=type*(cdrain-gds*vds-gm*vgs-gmbs*vbs)
      go to 920
  910 xnrm=0.0d0
      xrev=1.0d0
      cdreq=-type*(cdrain-gds*(-vds)-gm*vgd-gmbs*vbd)
  920 value(lvn+node2)=value(lvn+node2)-ceqqg
      value(lvn+node4)=value(lvn+node4)-ceqbs-ceqbd-ceqqb
      value(lvn+node5)=value(lvn+node5)-cdreq+ceqbd-ceqqd
      value(lvn+node6)=value(lvn+node6)+cdreq+ceqbs
     1   +ceqqg+ceqqb+ceqqd
c
c  load y matrix
c
      locy=lvn+nodplc(loc+27)
      value(locy)=value(locy)+gdpr
      locy=lvn+nodplc(loc+28)
      value(locy)=value(locy)+gcggb
      locy=lvn+nodplc(loc+29)
      value(locy)=value(locy)+gspr
      locy=lvn+nodplc(loc+30)
      value(locy)=value(locy)+gbd+gbs-gcbgb-gcbdb-gcbsb
      locy=lvn+nodplc(loc+31)
      value(locy)=value(locy)+gdpr+gds+gbd+xrev*(gm+gmbs)+gcddb
      locy=lvn+nodplc(loc+32)
      value(locy)=value(locy)+gspr+gds+gbs+xnrm*(gm+gmbs)+gcssb
      locy=lvn+nodplc(loc+10)
      value(locy)=value(locy)-gdpr
      locy=lvn+nodplc(loc+11)
      value(locy)=value(locy)-gcggb-gcgdb-gcgsb
      locy=lvn+nodplc(loc+12)
      value(locy)=value(locy)+gcgdb
      locy=lvn+nodplc(loc+13)
      value(locy)=value(locy)+gcgsb
      locy=lvn+nodplc(loc+14)
      value(locy)=value(locy)-gspr
      locy=lvn+nodplc(loc+15)
      value(locy)=value(locy)+gcbgb
      locy=lvn+nodplc(loc+16)
      value(locy)=value(locy)-gbd+gcbdb
      locy=lvn+nodplc(loc+17)
      value(locy)=value(locy)-gbs+gcbsb
      locy=lvn+nodplc(loc+18)
      value(locy)=value(locy)-gdpr
      locy=lvn+nodplc(loc+19)
      value(locy)=value(locy)+(xnrm-xrev)*gm+gcdgb
      locy=lvn+nodplc(loc+20)
      value(locy)=value(locy)-gbd+(xnrm-xrev)*gmbs-
     1   gcdgb-gcddb-gcdsb
      locy=lvn+nodplc(loc+21)
      value(locy)=value(locy)-gds-xnrm*(gm+gmbs)+gcdsb
      locy=lvn+nodplc(loc+22)
      value(locy)=value(locy)-(xnrm-xrev)*gm+gcsgb
      locy=lvn+nodplc(loc+23)
      value(locy)=value(locy)-gspr
      locy=lvn+nodplc(loc+24)
      value(locy)=value(locy)-gbs-(xnrm-xrev)*gmbs-
     1   gcsgb-gcsdb-gcssb
      locy=lvn+nodplc(loc+25)
      value(locy)=value(locy)-gds-xrev*(gm+gmbs)+gcsdb
 1000 loc=nodplc(loc)
      go to 10
      end
      subroutine mosq2(vds,vbs,vgs,vdsat,vth,vbin,gamasd,cox,phi,
     1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      implicit double precision (a-h,o-z)
c
c     initialize charges;
c     change reference voltages for charge computation
c
      qg=0.0d0
      qb=0.0d0
      vbd=vbs-vds
      vgb=vgs-vbs
      vd=dmax1(phi-vbd,1.0d-8)
      vs=dmax1(phi-vbs,1.0d-8)
      vg=vgb-vbin+phi
      vsp5=dsqrt(vs)
c
c     determine operating region
c
      if (vgs.le.vth) go to 1100
c
c     compute charges for "on" region
c
 1020 vsat=vdsat+vs
      vs2=vs*vs
      vs3=vs2*vs
      vs5=vs3*vs2
      vs1p5=vs*vsp5
      vs2p5=vs1p5*vs
 1025 if (vd.ge.vsat) go to 1035
      ve=vd
 1030 dvedvd=1.0d0
      dvedvg=0.0d0
      go to 1040
 1035 ve=vsat
      dvedvd=0.0d0
      dvedvg=0.0d0
 1040 ve2=ve*ve
      ve3=ve2*ve
      ve5=ve2*ve3
      vep5=dsqrt(ve)
      ve1p5=ve*vep5
      ve2p5=ve1p5*ve
      term0=ve+vs
      term1=vep5+vsp5
      term2=vep5*vsp5
      term3=ve2+vs2
      term4=ve*vs
      term5=term0*term1
      term6=(term3+term4)+term2*term0
      term7=(term3+term4)*term1
      term10=vep5+0.5d0*vsp5
      term11=1.5d0*ve+vsp5*term10
      term12=2.0d0*ve1p5+vsp5*term11
      term20=0.5d0*vep5+vsp5
      term21=1.5d0*vs+vep5*term20
      term22=2.0d0*vs1p5+vep5*term21
      argn=0.5d0*vg*term5-0.4d0*gamasd*term6-term7/3.0d0
      argd=vg*term1-gamasd*(term0+term2)/1.5d0-0.5d0*term1*term0
      argd2=argd*argd
      qg=cox*(vg-argn/argd)
      dgndve=0.5d0*vg*term11-0.4d0*gamasd*term12-
     1   (2.5d0*ve2+vsp5*term12)/3.0d0
      dddve=0.5d0*vg-gamasd*term10/1.5d0-0.5d0*term11
      dqgdve=-cox/argd*(dgndve-(vg-qg/cox)*dddve)
      dgndvs=0.5d0*vg*term21-0.4d0*gamasd*term22-
     1   (2.5d0*vs2+vep5*term22)/3.0d0
      dddvs=0.5d0*vg-gamasd*term20/1.5d0-0.5d0*term21
      cgdb=-cox/(argd*vep5)*(dgndve-(vg-qg/cox)*dddve)*dvedvd
      cgsb=-cox/(argd*vsp5)*(dgndvs-(vg-qg/cox)*dddvs)
      cggb=cox*(1.0d0-term1/argd*(0.5d0*term0-vg+qg/cox))
      argn=vg*(term0+term2)/1.5d0-0.5d0*gamasd*term5-0.4d0*term6
      dgndve=vg*term10/1.5d0-0.5d0*gamasd*term11-0.4d0*term12
      dgndvs=vg*term20/1.5d0-0.5d0*gamasd*term21-0.4d0*term22
      qb=-gamasd*cox*argn/argd
      cbdb=-cox/(vep5*argd)*(qb/cox*dddve+gamasd*dgndve)*dvedvd
      cbsb=-cox/(vsp5*argd)*(qb/cox*dddvs+gamasd*dgndvs)
      cbgb=-cox/argd*(gamasd*(term0+term2)/1.5d0+qb/cox*term1)
      go to 2000
c
c  finish special cases
c
 1100 if (vg.gt.0.0d0) go to 1110
      qg=cox*vg
      cggb=cox
      go to 1120
 1110 gamma2=gamasd*0.5d0
      sqarg=dsqrt(gamma2*gamma2+vg)
      qg=gamasd*cox*(sqarg-gamma2)
      cggb=0.5d0*cox*gamasd/sqarg
 1120 qb=-qg
      cbgb=-cggb
      cgdb=0.0d0
      cgsb=0.0d0
      cbdb=0.0d0
      cbsb=0.0d0
c
c  finished
c
 2000 qc=-(qg+qb)
 2050 return
      end
      subroutine mosq3(vds,vbs,vgs,vdsat,vth,vbin,gamasd,cox,phi,
     1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      implicit double precision (a-h,o-z)
c
      equivalence (xlamda,alpha),(vbp,theta),(uexp,eta),(utra,xkappa)
c
c     charge equations are referenced to bulk
c
      vgb=vgs-vbs
      vfb=vbi-phi
      onxl=1.0d0/xl
      phibs=sqphbs*sqphbs
c
c     body effect
c
      gammas=gamma*fshort
      fbodys=gammas/(sqphbs+sqphbs)*0.5d0
      fbody=fbodys+fnarrw
      onfbdy=1.0d0/(1.0d0+fbody)
      dfbdvb=-fbodys*dsqdvb/sqphbs+fbodys*dfsdvb/fshort
      qbonco=gammas*sqphbs+fnarrw*phibs
      dqbdvb=gammas*dsqdvb+gamma*dfsdvb*sqphbs-fnarrw
c
c.....static feedback effect
c
      vbix=vbi-eta*vds
c
c.....threshold voltage
c
      vth=vbix+qbonco
      dvtdvd=-eta
      dvtdvb=dqbdvb
c
c     branch according to region of operation
c
      if (vgs.le.vth) go to 800
      vgsx=dmax1(vgs,von)
c
c     branch on vds=0.0d0
c
      vdsx=dmin1(vds,vdsat)
      if ( vdsx.eq.0.0d0 ) go to 900
      cdo=vgsx-vth-0.5d0*(1.0d0+fbody)*vdsx
      dcodvg=1.0d0
      if (vds.lt.vdsat) dcodvd=-dvtdvd-0.5d0*(1.0d0+fbody)
      dcodvb=-dvtdvb-0.5d0*dfbdvb*vdsx
c
c.....charge terms
c
420   continue
      arga=(1.0d0+fbody)*vdsx*vdsx/(12.0d0*cdo)
      dadco=-arga/cdo
      if (vds.lt.vdsat) dadvd=arga/vdsx
      dadfb=arga*onfbdy
c
c.....gate charge
c
      qg=cox*(vgs-vbix-0.5d0*vdsx+arga)
      cggb=cox*(1.0d0+dadco*dcodvg)
      if (vds.lt.vdsat) cgdb=cox*(-dvtdvd-0.5d0+dadvd+dadco*dcodvd)
      cgsb=-cggb-cgdb-cox*(dadco*dcodvb+dadfb*dfbdvb)
c
c.....bulk charge
c
      arga=arga*fbody
      dadco=dadco*fbody
      if (vds.lt.vdsat) dadvd=dadvd*fbody
      dadfb=dadfb*(1.0d0+fbody+fbody)
c
      qb=-cox*(qbonco+0.5d0*fbody*vdsx-arga)
      cbgb=cox*dadco*dcodvg
      if (vds.lt.vdsat) cbdb=-cox*(0.5d0*fbody-dadvd-dadco*dcodvd)
      cbsb=-cbgb-cbdb
     1          +cox*(dqbdvb+(0.5d0*vdsx-dadfb)*dfbdvb-dadco*dcodvb)
      go to 1000
c
c.....charge terms of vgs<vth
c
800   continue
      if ( vgb.gt.vfb ) go to 810
      qg=cox*(vgb-vfb)
      cggb=cox
      go to  820
810   continue
      gamma2=gammas*0.5d0
      arga=dsqrt(gamma2*gamma2+(vgb-vfb))
      qg=gammas*cox*(arga-gamma2)
      cggb=0.5d0*cox*gammas/arga
820   continue
      qb=-qg
      cbgb=-cggb
      cgdb=0.0d0
      cgsb=0.0d0
      cbdb=0.0d0
      cbsb=0.0d0
      go to 1000
c
c     special case vds=0.0d0
c
  900 qg=cox*(vgs-vbi)
      qb=-cox*qbonco
      cggb=cox
      cgdb=-cox*(0.5d0+dvtdvd)
      cgsb=-cox*(0.5d0-dvtdvb)
      cbgb=0.0d0
      cbdb=-0.5d0*cox*fbody
      cbsb=cox*(dqbdvb+0.5d0*fbody)
c
c     done
c
 1000 qc=-(qg+qb)
      return
      end
      subroutine mqspof(vds,vbs,vgs,vpof,vdsat1,vth,vbin,gamasd,
     $qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      implicit double precision (a-h,o-z)
c
c spice version 2g.6  sccsid=mosarg 3/15/83
      common /mosarg/ vto,beta,gamma,phi,phib,cox,xnsub,xnfs,xd,xj,xld,
     1   xlamda,uo,uexp,vbp,utra,vmax,xneff,xl,xw,vbi,von,vdsat,qspof,
     2   beta0,beta1,cdrain,xqco,xqc,fnarrw,fshort,lev
c spice version 2g.6  sccsid=status 3/15/83
      common /status/ omega,time,delta,delold(7),ag(7),vt,xni,egfet,
     1   xmu,sfactr,mode,modedc,icalc,initf,method,iord,maxord,noncon,
     2   iterno,itemno,nosolv,modac,ipiv,ivmflg,ipostp,iscrch,iofile
c
c     vdsat1=dmax1(vds,vdsat1)+1.0d-3
      if( lev .eq. 3 ) goto 50
      if( lev .ne. 2 ) goto 1000
      call mosq2(vds,vbs,vgs,vdsat,vth,vbin,gamasd,cox,phi,
     $qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
      if (vds.ge.vdsat) go to 80
      call mosq2(vds,vbs,vpof,vdsat1,vth,vbin, gamasd,cox,phi,
     $qg1,qcpof1,qb1,cggb1,cgdb1,cgsb1,cbgb1,cbdb1,cbsb1)
      call mosq2(vdsat,vbs,vgs,vdsat,vth,vbin, gamasd,cox,phi,
     $qg2,qcpof2,qb2,cggb2,cgdb2,cgsb2,cbgb2,cbdb2,cbsb2)
      goto 75
   50 call mosq3(vds,vbs,vpof,vdsat1,vth,vbin, gamasd,cox,phi,
     $qg,qcpof,qb,cggb1,cgdb1,cgsb1,cbgb1,cbdb1,cbsb1)
      call mosq3(vds,vbs,vgs,vdsat,vth,vbin,gamasd,cox,phi,
     $qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb)
   75 if(vgs.gt.vpof. or .vds.lt.vdsat) goto 100
   80 xqc = xqco
      goto 1000
c
c     tangential limiting of qs
c
  100 csgb1=-(1.0d0-xqco)*(cggb1+cbgb1)
      qs=csgb1*(vgs-vpof)
     1   +(1.0d0-xqco)*qcpof1
c      write(iofile,*) "vgs,vds,qc,cggb,cgdb,cgsb,cbgb,cbdb,cbsb =",
c    1   vgs,vds,qc,cggb,cgdb,cgsb,cbgb,cbdb,cbsb
c      write(iofile,*) "vpof,vdsat,vdsat1,qcpof1,qcpof2,qs,csgb1 =",
c    1   vpof,vdsat,vdsat1,qcpof1,qcpof2,qs,csgb1
      qspof2=(1.0d0-xqco)*qcpof2
      if (dabs(qs) .lt. dabs(qspof2)) qs=qspof2
      if( dabs( qs ) .ge. 0.5d0 * dabs( qc ) ) goto 200
c     csdb=-0.25d0*(cgdb+cbdb)
c     qs=qs+csdb*(vdsat-vds)
c     xqc=dmin1(0.5d0,(qc-qs)/qc)
      xqc=0.5d0
c      write(iofile,*) "qs,xqc =",
c    1   qs,xqc
      goto 1000
  200 qd = qc - qs
      xqc = qd / qc
c     write(iofile,*) "200,qs,qd,xqc =",
c    1   qs,qd,xqc
c
c     constant limiting of qs
c
c 100 qdpof = qcpof * xqco
c     qspof = qcpof - qdpof
c     if( dabs( qspof ) .gt. 0.5d0 * dabs( qc ) ) goto 200
c     xqc = 0.5d0
c     goto 1000
c 200 qd = qc - qspof
c     qs = qspof
c     xqc = qd / qc
 1000 return
      end
