FTN66,Y,J,E 
*************************************************************** 
* (C) COPYRIGHT HEWLETT-PACKARD COMPANY 1980.  ALL RIGHTS     * 
* RESERVED.  NO PART OF THIS PROGRAM MAY BE PHOTOCOPIED, RE-  * 
* REPRODUCED, OR TRANSLATED TO ANOTHER PROGRAM LANGUAGE WITH- * 
* OUT THE PRIOR WRITTEN CONSENT OF HEWLETT-PACKARD COMPANY.   * 
*************************************************************** 
*  NAME:   DCTR2                      * 
*  SOURCE: 92091-18022                * 
*  RELOC:  92091-16022                * 
*  PGMR:   C. GEBER                   * 
*************************************** 
$PAGE 
$EMA /BLANK/
$CDS ON 
      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 ) 
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (A-H,O-Z) 
      COMMON /MOSARG/ VTO,BETA,GAMMA,PHI,PHIB,COX,XNSUB,XNFS,XD,XJ,XLD, 
     1   XLAMDA,UO,UEXP,VBP,UTRA,VMAX,XNEFF,XL,XW,VBI,VON,VDSAT,QSPOF,
     2   BETA0,BETA1,CDRAIN,XQCO,XQC,FNARRW,FSHORT,LEV
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
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)
$PAGE 
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 )
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (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 
      COMMON /MOSARG/ VTO,BETA,GAMMA,PHI,PHIB,COX,XNSUB,XNFS,XD,XJ,XLD, 
     1   XLAMDA,UO,UEXP,VBP,UTRA,VMAX,XNEFF,XL,XW,VBI,VON,VDSAT,QSPOF,
     2   BETA0,BETA1,CDRAIN,XQCO,XQC,FNARRW,FSHORT,LEV
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      VBD=VBS-VDS 
      VGB=VGS-VBS 
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
$PAGE 
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 )
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (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 
      COMMON /MOSARG/ VTO,BETA,GAMMA,PHI,PHIB,COX,XNSUB,XNFS,XD,XJ,XLD, 
     1   XLAMDA,UO,UEXP,VBP,UTRA,VMAX,XNEFF,XL,XW,VBI,VON,VDSAT,QSPOF,
     2   BETA0,BETA1,CDRAIN,XQCO,XQC,FNARRW,FSHORT,LEV
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
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 
$PAGE 
  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
$PAGE 
C     WRITE(6,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 (6,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
$PAGE 
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
$PAGE 
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 
  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 
  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
$PAGE 
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 )
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (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 
      COMMON /MOSARG/ VTO,BETA,GAMMA,PHI,PHIB,COX,XNSUB,XNFS,XD,XJ,XLD, 
     1   XLAMDA,UO,UEXP,VBP,UTRA,VMAX,XNEFF,XL,XW,VBI,VON,VDSAT,QSPOF,
     2   BETA0,BETA1,CDRAIN,XQCO,XQC,FNARRW,FSHORT,LEV
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
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)
$PAGE 
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 
$PAGE 
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 
$PAGE 
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
$PAGE 
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 
$PAGE 
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 
$PAGE 
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 
$EMA /BLANK/                            
      SUBROUTINE MOSFET 
     + ,92091-16022 REV.5000 <890418.1141>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE PROCESSES MOSFETS FOR DC AND TRANSIENT ANALYSES. 
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /MOSARG/ VTO,BETA,GAMMA,PHI,PHIB,COX,XNSUB,XNFS,XD,XJ,XLD, 
     1   XLAMDA,UO,UEXP,VBP,UTRA,VMAX,XNEFF,XL,XW,VBI,VON,VDSAT,QSPOF,
     2   BETA0,BETA1,CDRAIN,XQCO,XQC,FNARRW,FSHORT,LEV
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      COMMON/DEBUG/ IDEBUG(20)
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
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)
$PAGE 
      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 
      LOC=LOCATE(14)
   10 IF (LOC.EQ.0) RETURN
      LOCV=NODPLC(LOC+1)
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NODE3=NODPLC(LOC+4) 
      NODE4=NODPLC(LOC+5) 
      NODE5=NODPLC(LOC+6) 
      NODE6=NODPLC(LOC+7) 
      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 
$PAGE 
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)
$PAGE 
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 
$PAGE 
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) 
  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 (IXXOR(DEVMOD,-1.0D0).EQ.0) 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 
$PAGE 
  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 
C.. INTEGRATE QB
C 
  810 CALL INTGR8(GEQ,CEQ,0.0D0,LOCT+12)
C 
C.. INTEGRATE QG
C 
      CALL INTGR8(GEQ,CEQ,0.0D0,LOCT+14)
C 
C.. INTEGRATE QD
C 
      CALL INTGR8(GEQ,CEQ,0.0D0,LOCT+16)
      GO TO 860 
$PAGE 
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
$PAGE 
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 , 
     1    COX , PHI , QG , QC , QB , CGGB , CGDB , CGSB , CBGB , CBDB , 
     2    CBSB )
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (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 , 
     1   COX , PHI , QG , QC , QB , CGGB , CGDB , CGSB , CBGB , CBDB ,
     2   CBSB ) 
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (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
$PAGE 
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
$PAGE 
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 ) 
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
      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 
C     VDSAT1=DMAX1(VDS,VDSAT1)+1.0D-3 
C 
      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 
C      write(6,*) "VGS,VDS,QC,CGGB,CGDB,CGSB,CBGB,CBDB,CBSB =", 
C    1   vgs,vds,qc,cggb,cgdb,cgsb,cbgb,cbdb,cbsb 
C      write(6,*) "VPOF,VDSAT,VDSAT1,QCPOF1,QCPOF2,QS,CSGB1 =", 
C    1   vpof,vdsat,vdsat1,qcpof1,qcpof2,qs,csgb1 
C 
      QSPOF2=(1.0D0-XQCO)*QCPOF2
      IF (DABS(QS) .LT. DABS(QSPOF2)) QS=QSPOF2 
      IF ( DABS( QS ) .GE. 0.5D0 * DABS( QC ) ) GOTO 200
C 
C     CSDB=-0.25D0*(CGDB+CBDB)
C     QS=QS+CSDB*(VDSAT-VDS)
C     XQC=DMIN1(0.5D0,(QC-QS)/QC) 
C 
      XQC=0.5D0 
C 
C      write(6,*) "QS,XQC =", 
C    1   qs,xqc 
C 
      GOTO 1000 
  200 QD = QC - QS
      XQC = QD / QC 
C 
C     write(6,*) "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 
C 
 1000 RETURN
      END 
$EMA /BLANK/
      SUBROUTINE NLCSRC 
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE LOADS THE NONLINEAR CONTROLLED SOURCES INTO THE
C COEFFICIENT MATRIX. 
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /FLAGS/ IPRNTA,IPRNTL,IPRNTM,IPRNTN,IPRNTO,LIMTIM,LIMPTS,
     1   LVLCOD,LVLTIM,ITL1,ITL2,ITL3,ITL4,ITL5,IGOOF,NOGO,KEOF 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
C  NONLINEAR VOLTAGE-CONTROLLED CURRENT SOURCES 
C 
      LOC=LOCATE(5) 
   10 IF (LOC.EQ.0) GO TO 100 
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NDIM=NODPLC(LOC+4)
      LNOD=NODPLC(LOC+6)
      LMAT=NODPLC(LOC+7)
      LCOEF=NODPLC(LOC+8) 
      CALL SIZMEM(NODPLC(LOC+8),NCOEF)
      LARG=NODPLC(LOC+9)
      LEXP=NODPLC(LOC+10) 
      LIC=NODPLC(LOC+11)
      LOCT=NODPLC(LOC+12)+1 
      ICHECK=0
      DO 20 I=1,NDIM
      CALL UPDATE(VALUE(LIC+I),LOCT,NODPLC(LNOD+1),NODPLC(LNOD+2),2,
     1   ICHECK)
      VALUE(LARG+I)=VALUE(LX0+LOCT) 
      LOCT=LOCT+2 
      LNOD=LNOD+2 
   20 CONTINUE
$PAGE 
      CALL EVPOLY(COLD,0,LCOEF,NCOEF,LARG,NDIM,LEXP)
      LOCT=NODPLC(LOC+12) 
      IF (ICHECK.EQ.1) GO TO 30 
      IF (INITF.EQ.6) GO TO 30
      TOL=RELTOL*DMAX1(DABS(COLD),DABS(VALUE(LX0+LOCT)))+ABSTOL 
      IF (DABS(COLD-VALUE(LX0+LOCT)).LT.TOL) GO TO 40 
   30 NONCON=NONCON+1 
   40 VALUE(LX0+LOCT)=COLD
      CEQ=COLD
      DO 50 I=1,NDIM
      CALL EVPOLY(GEQ,I,LCOEF,NCOEF,LARG,NDIM,LEXP) 
      LOCT=LOCT+2 
      VALUE(LX0+LOCT)=GEQ 
      CEQ=CEQ-GEQ*VALUE(LARG+I) 
      LOCY=LVN+NODPLC(LMAT+1) 
      VALUE(LOCY)=VALUE(LOCY)+GEQ 
      LOCY=LVN+NODPLC(LMAT+2) 
      VALUE(LOCY)=VALUE(LOCY)-GEQ 
      LOCY=LVN+NODPLC(LMAT+3) 
      VALUE(LOCY)=VALUE(LOCY)-GEQ 
      LOCY=LVN+NODPLC(LMAT+4) 
      VALUE(LOCY)=VALUE(LOCY)+GEQ 
      LMAT=LMAT+4 
   50 CONTINUE
      VALUE(LVN+NODE1)=VALUE(LVN+NODE1)-CEQ 
      VALUE(LVN+NODE2)=VALUE(LVN+NODE2)+CEQ 
      LOC=NODPLC(LOC) 
      GO TO 10
C 
C  NONLINEAR VOLTAGE CONTROLLED VOLTAGE SOURCES 
C 
  100 LOC=LOCATE(6) 
  110 IF (LOC.EQ.0) GO TO 200 
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NDIM=NODPLC(LOC+4)
      IPTR=NODPLC(LOC+6)
      LNOD=NODPLC(LOC+7)
      LMAT=NODPLC(LOC+8)
      LCOEF=NODPLC(LOC+9) 
      CALL SIZMEM(NODPLC(LOC+9),NCOEF)
      LARG=NODPLC(LOC+10) 
      LEXP=NODPLC(LOC+11) 
      LIC=NODPLC(LOC+12)
      LOCT=NODPLC(LOC+13)+2 
      ICHECK=0
      DO 120 I=1,NDIM 
      CALL UPDATE(VALUE(LIC+I),LOCT,NODPLC(LNOD+1),NODPLC(LNOD+2),2,
     1   ICHECK)
      VALUE(LARG+I)=VALUE(LX0+LOCT) 
      LOCT=LOCT+2 
      LNOD=LNOD+2 
  120 CONTINUE
      CALL EVPOLY(VOLT,0,LCOEF,NCOEF,LARG,NDIM,LEXP)
      LOCT=NODPLC(LOC+13) 
      IF (ICHECK.EQ.1) GO TO 130
      IF (INITF.EQ.6) GO TO 130 
      TOL=RELTOL*DMAX1(DABS(VOLT),DABS(VALUE(LX0+LOCT)))+VNTOL
      IF (DABS(VOLT-VALUE(LX0+LOCT)).LT.TOL) GO TO 140
  130 NONCON=NONCON+1 
  140 VALUE(LX0+LOCT)=VOLT
      VALUE(LX0+LOCT+1)=VALUE(LVNIM1+IPTR)
      VEQ=VOLT
      LOCY=LVN+NODPLC(LMAT+1) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LMAT+2) 
      VALUE(LOCY)=-1.0D0
      LOCY=LVN+NODPLC(LMAT+3) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LMAT+4) 
      VALUE(LOCY)=-1.0D0
      LMAT=LMAT+4 
      LOCT=LOCT+1 
      DO 150 I=1,NDIM 
      CALL EVPOLY(VGAIN,I,LCOEF,NCOEF,LARG,NDIM,LEXP) 
      LOCT=LOCT+2 
      VALUE(LX0+LOCT)=VGAIN 
      VEQ=VEQ-VGAIN*VALUE(LARG+I) 
      LOCY=LVN+NODPLC(LMAT+1) 
      VALUE(LOCY)=VALUE(LOCY)-VGAIN 
      LOCY=LVN+NODPLC(LMAT+2) 
      VALUE(LOCY)=VALUE(LOCY)+VGAIN 
      LMAT=LMAT+2 
  150 CONTINUE
      VALUE(LVN+IPTR)=VEQ 
      LOC=NODPLC(LOC) 
      GO TO 110 
C 
C  NONLINEAR CURRENT-CONTROLLED CURRENT SOURCES 
C 
  200 LOC=LOCATE(7) 
  210 IF (LOC.EQ.0) GO TO 300 
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NDIM=NODPLC(LOC+4)
      LVS=NODPLC(LOC+6) 
      LMAT=NODPLC(LOC+7)
      LCOEF=NODPLC(LOC+8) 
      CALL SIZMEM(NODPLC(LOC+8),NCOEF)
      LARG=NODPLC(LOC+9)
      LEXP=NODPLC(LOC+10) 
      LIC=NODPLC(LOC+11)
      LOCT=NODPLC(LOC+12)+1 
      ICHECK=0
      DO 220 I=1,NDIM 
      IPTR=NODPLC(LVS+I)
      IPTR=NODPLC(IPTR+6) 
      CALL UPDATE(VALUE(LIC+I),LOCT,IPTR,1,2,ICHECK)
      VALUE(LARG+I)=VALUE(LX0+LOCT) 
      LOCT=LOCT+2 
  220 CONTINUE
      CALL EVPOLY(CSRC,0,LCOEF,NCOEF,LARG,NDIM,LEXP)
      LOCT=NODPLC(LOC+12) 
      IF (ICHECK.EQ.1) GO TO 230
      IF (INITF.EQ.6) GO TO 230 
      TOL=RELTOL*DMAX1(DABS(CSRC),DABS(VALUE(LX0+LOCT)))+ABSTOL 
      IF (DABS(CSRC-VALUE(LX0+LOCT)).LT.TOL) GO TO 240
  230 NONCON=NONCON+1 
  240 VALUE(LX0+LOCT)=CSRC
      CEQ=CSRC
      DO 250 I=1,NDIM 
      CALL EVPOLY(CGAIN,I,LCOEF,NCOEF,LARG,NDIM,LEXP) 
      LOCT=LOCT+2 
      VALUE(LX0+LOCT)=CGAIN 
      CEQ=CEQ-CGAIN*VALUE(LARG+I) 
      LOCY=LVN+NODPLC(LMAT+1) 
      VALUE(LOCY)=VALUE(LOCY)+CGAIN 
      LOCY=LVN+NODPLC(LMAT+2) 
      VALUE(LOCY)=VALUE(LOCY)-CGAIN 
      LMAT=LMAT+2 
  250 CONTINUE
      VALUE(LVN+NODE1)=VALUE(LVN+NODE1)-CEQ 
      VALUE(LVN+NODE2)=VALUE(LVN+NODE2)+CEQ 
      LOC=NODPLC(LOC) 
      GO TO 210 
C 
C  NONLINEAR CURRENT CONTROLLED VOLTAGE SOURCES 
C 
  300 LOC=LOCATE(8) 
  310 IF (LOC.EQ.0) GO TO 1000
      NODE1=NODPLC(LOC+2) 
      NODE2=NODPLC(LOC+3) 
      NDIM=NODPLC(LOC+4)
      IBR=NODPLC(LOC+6) 
      LVS=NODPLC(LOC+7) 
      LMAT=NODPLC(LOC+8)
      LCOEF=NODPLC(LOC+9) 
      CALL SIZMEM(NODPLC(LOC+9),NCOEF)
      LARG=NODPLC(LOC+10) 
      LEXP=NODPLC(LOC+11) 
      LIC=NODPLC(LOC+12)
      LOCT=NODPLC(LOC+13)+2 
      ICHECK=0
      DO 320 I=1,NDIM 
      IPTR=NODPLC(LVS+I)
      IPTR=NODPLC(IPTR+6) 
      CALL UPDATE(VALUE(LIC+I),LOCT,IPTR,1,2,ICHECK)
      VALUE(LARG+I)=VALUE(LX0+LOCT) 
      LOCT=LOCT+2 
  320 CONTINUE
      CALL EVPOLY(VOLT,0,LCOEF,NCOEF,LARG,NDIM,LEXP)
      LOCT=NODPLC(LOC+13) 
      IF (ICHECK.EQ.1) GO TO 330
      IF (INITF.EQ.6) GO TO 330 
      TOL=RELTOL*DMAX1(DABS(VOLT),DABS(VALUE(LX0+LOCT)))+VNTOL
      IF (DABS(VOLT-VALUE(LX0+LOCT)).LT.TOL) GO TO 340
  330 NONCON=NONCON+1 
  340 VALUE(LX0+LOCT)=VOLT
      VALUE(LX0+LOCT+1)=VALUE(LVNIM1+IBR) 
      VEQ=VOLT
      LOCY=LVN+NODPLC(LMAT+1) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LMAT+2) 
      VALUE(LOCY)=-1.0D0
      LOCY=LVN+NODPLC(LMAT+3) 
      VALUE(LOCY)=+1.0D0
      LOCY=LVN+NODPLC(LMAT+4) 
      VALUE(LOCY)=-1.0D0
      LMAT=LMAT+4 
      LOCT=LOCT+1 
      DO 350 I=1,NDIM 
      CALL EVPOLY(TRANSR,I,LCOEF,NCOEF,LARG,NDIM,LEXP)
      LOCT=LOCT+2 
      VALUE(LX0+LOCT)=TRANSR
      VEQ=VEQ-TRANSR*VALUE(LARG+I)
      LOCY=LVN+NODPLC(LMAT+I) 
      VALUE(LOCY)=VALUE(LOCY)-TRANSR
  350 CONTINUE
      VALUE(LVN+IBR)=VEQ
      LOC=NODPLC(LOC) 
      GO TO 310 
C 
C  FINISHED 
C 
 1000 RETURN
      END 
      SUBROUTINE NXTPWR( PWRSEQ , PDIM )
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE DETERMINES THE 'NEXT' SET OF EXPONENTS FOR THE 
C DIFFERENT DIMENSIONS OF A POLYNOMIAL. 
C 
      INTEGER PWRSEQ(1),PDIM,PSUM 
C 
      IF (PDIM.EQ.1) GO TO 80 
      K=PDIM
   10 IF (PWRSEQ(K).NE.0) GO TO 20
      K=K-1 
      IF (K.NE.0) GO TO 10
      GO TO 80
   20 IF (K.EQ.PDIM) GO TO 30 
      PWRSEQ(K)=PWRSEQ(K)-1 
      PWRSEQ(K+1)=PWRSEQ(K+1)+1 
      GO TO 100 
   30 KM1=K-1 
      DO 40 I=1,KM1 
      IF (PWRSEQ(I).NE.0) GO TO 50
   40 CONTINUE
      PWRSEQ(1)=PWRSEQ(PDIM)+1
      PWRSEQ(PDIM)=0
      GO TO 100 
   50 PSUM=1
      K=PDIM
   60 IF (PWRSEQ(K-1).GE.1) GO TO 70
      PSUM=PSUM+PWRSEQ(K) 
      PWRSEQ(K)=0 
      K=K-1 
      GO TO 60
   70 PWRSEQ(K)=PWRSEQ(K)+PSUM
      PWRSEQ(K-1)=PWRSEQ(K-1)-1 
      GO TO 100 
   80 PWRSEQ(1)=PWRSEQ(1)+1 
C 
C  FINISHED 
C 
  100 RETURN
      END 
      SUBROUTINE PNJLIM( VNEW , VOLD , VT , VCRIT , ICHECK )
     + ,92091-16022 REV.5000 <890418.1141>
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE LIMITS THE CHANGE-PER-ITERATION OF DEVICE PN-JUNCTION
C VOLTAGES. 
C 
      IF (VNEW.LE.VCRIT) GO TO 30 
      VLIM=VT+VT
      DELV=VNEW-VOLD
      IF (DABS(DELV).LE.VLIM) GO TO 30
      IF (VOLD.LE.0.0D0) GO TO 20 
      ARG=1.0D0+DELV/VT 
      IF (ARG.LE.0.0D0) GO TO 10
      VNEW=VOLD+VT*DLOG(ARG)
      GO TO 100 
   10 VNEW=VCRIT
      GO TO 100 
   20 VNEW=VT*DLOG(VNEW/VT) 
      GO TO 100 
   30 ICHECK=0
C 
C  FINISHED 
C 
  100 RETURN
      END 
$EMA /BLANK/                            
      SUBROUTINE SORUPD 
     + ,92091-16022 REV.5000 <890418.1141>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE UPDATES THE INDEPENDENT VOLTAGE AND CURRENT SOURCES
C USED IN THE CIRCUIT.  IT ALSO UPDATES THE LTD TABLE (WHICH CONTAINS 
C PREVIOUS (DELAYED) VALUES OF THE SOURCES USED TO MODEL TRANSMISSION 
C LINES). 
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /FLAGS/ IPRNTA,IPRNTL,IPRNTM,IPRNTN,IPRNTO,LIMTIM,LIMPTS,
     1   LVLCOD,LVLTIM,ITL1,ITL2,ITL3,ITL4,ITL5,IGOOF,NOGO,KEOF 
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      DO 500 ID=9,10
      LOC=LOCATE(ID)
   10 IF (LOC.EQ.0) GO TO 500 
      LOCV=NODPLC(LOC+1)
      LOCP=NODPLC(LOC+5)
      ITYPE=NODPLC(LOC+4)+1 
      GO TO (490,100,200,300,400,450), ITYPE
$PAGE 
C 
C  PULSE SOURCE 
C 
  100 V1=VALUE(LOCP+1)
      V2=VALUE(LOCP+2)
      T1=VALUE(LOCP+3)
      T2=VALUE(LOCP+4)
      T3=VALUE(LOCP+5)
      T4=VALUE(LOCP+6)
      PERIOD=VALUE(LOCP+7)
      TIME1=TIME
      IF (TIME1.LE.0.0D0) GO TO 160 
  110 IF (TIME1.LT.T1+PERIOD) GO TO 120 
      TIME1=TIME1-PERIOD
      GO TO 110 
  120 IF (TIME1.LT.T4) GO TO 130
      VALUE(LOCV+1)=V1
      GO TO 490 
  130 IF (TIME1.LT.T3) GO TO 140
      VALUE(LOCV+1)=V2+(TIME1-T3)*(V1-V2)/(T4-T3) 
      GO TO 490 
  140 IF (TIME1.LT.T2) GO TO 150
      VALUE(LOCV+1)=V2
      GO TO 490 
  150 IF (TIME1.LT.T1) GO TO 160
      VALUE(LOCV+1)=V1+(TIME1-T1)*(V2-V1)/(T2-T1) 
      GO TO 490 
  160 VALUE(LOCV+1)=V1
      GO TO 490 
C 
C  SINUSOIDAL SOURCE
C 
  200 V1=VALUE(LOCP+1)
      V2=VALUE(LOCP+2)
      OMEG=VALUE(LOCP+3)
      T1=VALUE(LOCP+4)
      THETA=VALUE(LOCP+5) 
      TIME1=TIME-T1 
      IF (TIME1.GT.0.0D0) GO TO 210 
      VALUE(LOCV+1)=V1
      GO TO 490 
  210 IF (THETA.NE.0.0D0) GO TO 220 
      VALUE(LOCV+1)=V1+V2*DSIN(OMEG*TIME1)
      GO TO 490 
  220 VALUE(LOCV+1)=V1+V2*DSIN(OMEG*TIME1)*DEXP(-TIME1*THETA) 
      GO TO 490 
$PAGE 
C 
C  EXPONENTIAL SOURCE 
C 
  300 V1=VALUE(LOCP+1)
      V2=VALUE(LOCP+2)
      T1=VALUE(LOCP+3)
      TAU1=VALUE(LOCP+4)
      T2=VALUE(LOCP+5)
      TAU2=VALUE(LOCP+6)
      TIME1=TIME
      IF (TIME1.GT.T1) GO TO 310
      VALUE(LOCV+1)=V1
      GO TO 490 
  310 IF (TIME1.GT.T2) GO TO 320
      VALUE(LOCV+1)=V1+(V2-V1)*(1.0D0-DEXP((T1-TIME1)/TAU1))
      GO TO 490 
  320 VALUE(LOCV+1)=V1+(V2-V1)*(1.0D0-DEXP((T1-TIME1)/TAU1))
     1   +(V1-V2)*(1.0D0-DEXP((T2-TIME1)/TAU2)) 
      GO TO 490 
C 
C  PIECEWISE-LINEAR SOURCE
C 
  400 T1=VALUE(LOCP+1)
      V1=VALUE(LOCP+2)
      T2=VALUE(LOCP+3)
      V2=VALUE(LOCP+4)
      IKNT=4
  410 IF (TIME.LE.T2) GO TO 420 
      T1=T2 
      V1=V2 
      T2=VALUE(LOCP+IKNT+1) 
      V2=VALUE(LOCP+IKNT+2) 
      IKNT=IKNT+2 
      GO TO 410 
  420 VALUE(LOCV+1)=V1+((TIME-T1)/(T2-T1))*(V2-V1)
      GO TO 490 
C 
C  SINGLE-FREQUENCY FM
C 
  450 V1=VALUE(LOCP+1)
      V2=VALUE(LOCP+2)
      OMEGC=VALUE(LOCP+3) 
      XMOD=VALUE(LOCP+4)
      OMEGS=VALUE(LOCP+5) 
      VALUE(LOCV+1)=V1+V2*DSIN(OMEGC*TIME+XMOD*DSIN(OMEGS*TIME))
  490 LOC=NODPLC(LOC) 
      GO TO 10
  500 CONTINUE
$PAGE 
C 
C  UPDATE TRANSMISSION LINE SOURCES 
C 
      IF (JELCNT(17).EQ.0) GO TO 1000 
      IF (MODE.NE.2) GO TO 1000 
      CALL SIZMEM(LTD,LTDSIZ) 
      NUMTD=LTDSIZ/NTLIN
      IF (NUMTD.LT.3) GO TO 900 
      LOC=LOCATE(17)
  610 IF (LOC.EQ.0) GO TO 1000
      LOCV=NODPLC(LOC+1)
      TD=VALUE(LOCV+2)
      BAKTIM=TIME-TD
      IF (BAKTIM.LT.0.0D0) GO TO 640
      LTDPTR=NODPLC(LOC+30) 
      ICNTR=2 
      L1=LTD
      L2=L1+NTLIN 
      L3=L2+NTLIN 
      T1=VALUE(L1+1)
      T2=VALUE(L2+1)
  620 T3=VALUE(L3+1)
      ICNTR=ICNTR+1 
      IF (BAKTIM.LE.T3) GO TO 630 
      IF (ICNTR.EQ.NUMTD) GO TO 900 
      L1=L2 
      L2=L3 
      L3=L2+NTLIN 
      T1=T2 
      T2=T3 
      GO TO 620 
  630 DT1T2=T1-T2 
      DT1T3=T1-T3 
      DT2T3=T2-T3 
      TDNOM1=1.0D0/(DT1T2*DT1T3)
      TDNOM2=-1.0D0/(DT1T2*DT2T3) 
      TDNOM3=1.0D0/(DT2T3*DT1T3)
      DTT1=BAKTIM-T1
      DTT2=BAKTIM-T2
      DTT3=BAKTIM-T3
      TFACT1=DTT2*DTT3*TDNOM1 
      TFACT2=DTT1*DTT3*TDNOM2 
      TFACT3=DTT1*DTT2*TDNOM3 
      VALUE(LOCV+3)=VALUE(L1+LTDPTR+0)*TFACT1+VALUE(L2+LTDPTR+0)*TFACT2 
     1   +VALUE(L3+LTDPTR+0)*TFACT3 
      VALUE(LOCV+4)=VALUE(L1+LTDPTR+1)*TFACT1+VALUE(L2+LTDPTR+1)*TFACT2 
     1   +VALUE(L3+LTDPTR+1)*TFACT3 
  640 LOC=NODPLC(LOC) 
      GO TO 610 
$PAGE 
C 
C  INTERNAL LOGIC ERROR:  LESS THAN 3 ENTRIES IN LTD
C 
  900 NOGO=1
      WRITE (6,901) NUMTD,ICNTR 
  901 FORMAT('0*ABORT*:  INTERNAL SPICE ERROR:  SORUPD:  ',2I5/)
C 
C  FINISHED 
C 
 1000 RETURN
      END 
$EMA /BLANK/                            
      SUBROUTINE TERR( LOCT , DELNEW )
     + ,92091-16022 REV.5000 <890418.1141>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE ESTIMATES THE LOCAL TRUNCATION ERROR FOR A PARTICULAR
C CIRCUIT ELEMENT.  IT THEN COMPUTES THE APPROPRIATE STEPSIZE WHICH 
C SHOULD BE USED. 
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /KNSTNT/ TWOPI,XLOG2,XLOG10,ROOT2,RAD,BOLTZ,CHARGE,CTOK, 
     1   GMIN,RELTOL,ABSTOL,VNTOL,TRTOL,CHGTOL,EPS0,EPSSIL,EPSOX, 
     2   PIVTOL,PIVREL
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      DIMENSION QCAP(1),CCAP(1),DIFF(8),DELTMP(7),COEF(6) 
      EQUIVALENCE (QCAP(1),VALUE(1)),(CCAP(1),VALUE(2)) 
      DATA COEF / 5.000000000D-1, 2.222222222D-1, 1.363636364D-1, 
     1            9.600000000D-2, 7.299270073D-2, 5.830903790D-2 /
      DATA XTWELV / 8.333333333D-2 /
C 
      TOL=RELTOL*DMAX1(DABS(CCAP(LX0+LOCT)),DABS(CCAP(LX1+LOCT)))+ABSTOL
      CTOL=RELTOL*DMAX1(DABS(QCAP(LX0+LOCT)),DABS(QCAP(LX1+LOCT)),
     1   CHGTOL)/DELTA
      TOL=DMAX1(TOL,CTOL) 
$PAGE 
C 
C  DETERMINE DIVIDED DIFFERENCES
C 
      GO TO (6,5,4,3,2,1), IORD 
    1 DIFF(8)=QCAP(LX7+LOCT)
    2 DIFF(7)=QCAP(LX6+LOCT)
    3 DIFF(6)=QCAP(LX5+LOCT)
    4 DIFF(5)=QCAP(LX4+LOCT)
    5 DIFF(4)=QCAP(LX3+LOCT)
    6 DIFF(3)=QCAP(LX2+LOCT)
      DIFF(2)=QCAP(LX1+LOCT)
      DIFF(1)=QCAP(LX0+LOCT)
      ISTOP=IORD+1
      DO 10 I=1,ISTOP 
      DELTMP(I)=DELOLD(I) 
   10 CONTINUE
   20 DO 30 I=1,ISTOP 
      DIFF(I)=(DIFF(I)-DIFF(I+1))/DELTMP(I) 
   30 CONTINUE
      ISTOP=ISTOP-1 
      IF (ISTOP.EQ.0) GO TO 100 
      DO 40 I=1,ISTOP 
      DELTMP(I)=DELTMP(I+1)+DELOLD(I) 
   40 CONTINUE
      GO TO 20
C 
C  DIFF(1) CONTAINS DIVIDED DIFFERENCE
C 
  100 CONST=COEF(IORD)
      IF ((METHOD.EQ.1).AND.(IORD.EQ.2)) CONST=XTWELV 
      DEL=TRTOL*TOL/DMAX1(ABSTOL,CONST*DABS(DIFF(1))) 
      IF (IORD.EQ.1) GO TO 200
      IF (IORD.GE.3) GO TO 150
      DEL=DSQRT(DEL)
      GO TO 200 
  150 DEL=DEXP(DLOG(DEL)/DFLOAT(IORD))
  200 DELNEW=DMIN1(DELNEW,DEL)
      RETURN
      END 
$EMA /BLANK/                            
      SUBROUTINE TRUNC( DELNEW )
     + ,92091-16022 REV.5000 <890418.1141>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE DETERMINES THE NEW TRANSIENT STEPSIZE BY EITHER
C CALLING TERR TO ESTIMATE THE LOCAL TRUNCATION ERROR, OR BY CHECKING 
C ON THE NUMBER OF ITERATIONS NEEDED TO CONVERGE AT THE LAST TIMEPOINT. 
C 
      COMMON /CIRDAT/ LOCATE(50),JELCNT(50),NUNODS,NCNODS,NUMNOD,NSTOP, 
     1   NUT,NLT,NXTRM,NDIST,NTLIN,IBR,NUMVS
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /FLAGS/ IPRNTA,IPRNTL,IPRNTM,IPRNTN,IPRNTO,LIMTIM,LIMPTS,
     1   LVLCOD,LVLTIM,ITL1,ITL2,ITL3,ITL4,ITL5,IGOOF,NOGO,KEOF 
      COMMON /TRAN/ TSTEP,TSTOP,TSTART,DELMAX,TDMAX,FORFRE,JTRFLG 
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      IF (LVLTIM.NE.0) GO TO 5
      DELNEW=DMIN1(TSTEP,DELMAX)
      RETURN
    5 IF (LVLTIM.NE.1) GO TO 10 
      DELNEW=DELTA
      IF (ITERNO.GT.ITL3) RETURN
      DELNEW=DMIN1(2.0D0*DELTA,TSTEP,DELMAX)
      RETURN
C 
C  CAPACITORS 
C 
   10 DELNEW=1.0D20 
      LOC=LOCATE(2) 
   20 IF (LOC.EQ.0) GO TO 30
      LOCT=NODPLC(LOC+8)
      CALL TERR(LOCT,DELNEW)
      LOC=NODPLC(LOC) 
      GO TO 20
C 
C  INDUCTORS
C 
   30 LOC=LOCATE(3) 
   40 IF (LOC.EQ.0) GO TO 50
      LOCT=NODPLC(LOC+11) 
      CALL TERR(LOCT,DELNEW)
      LOC=NODPLC(LOC) 
      GO TO 40
$PAGE 
C 
C  DIODES 
C 
   50 LOC=LOCATE(11)
   60 IF (LOC.EQ.0) GO TO 70
      LOCT=NODPLC(LOC+11) 
      CALL TERR(LOCT+3,DELNEW)
      LOC=NODPLC(LOC) 
      GO TO 60
C 
C  BJTS 
C 
   70 LOC=LOCATE(12)
   80 IF (LOC.EQ.0) GO TO 90
      LOCT=NODPLC(LOC+22) 
      CALL TERR(LOCT+8,DELNEW)
      CALL TERR(LOCT+10,DELNEW) 
      CALL TERR(LOCT+12,DELNEW) 
      LOC=NODPLC(LOC) 
      GO TO 80
C 
C  JFETS
C 
   90 LOC=LOCATE(13)
  100 IF (LOC.EQ.0) GO TO 110 
      LOCT=NODPLC(LOC+19) 
      CALL TERR(LOCT+9,DELNEW)
      CALL TERR(LOCT+11,DELNEW) 
      LOC=NODPLC(LOC) 
      GO TO 100 
C 
C  MOSFETS
C 
  110 LOC=LOCATE(14)
  120 IF (LOC.EQ.0) GO TO 200 
      LOCT=NODPLC(LOC+26) 
      CALL TERR(LOCT+12,DELNEW) 
      CALL TERR(LOCT+14,DELNEW) 
      CALL TERR(LOCT+16,DELNEW) 
      LOC=NODPLC(LOC) 
      GO TO 120 
C 
C  DELTA IS ALLOWED ONLY TO DOUBLE AT EACH TIMEPOINT
C 
  200 DELNEW=DMIN1(2.0D0*DELTA,DELNEW,DELMAX) 
      RETURN
      END 
$EMA /BLANK/                            
      SUBROUTINE UPDATE( VINIT , LOCT , NODE1 , NODE2 , NUPDA , ICHECK )
     + ,92091-16022 REV.5000 <890418.1141>  
      IMPLICIT REAL*8 (A-H,O-Z) 
C 
C     THIS ROUTINE UPDATES AND LIMITS THE CONTROLLING VARIABLES FOR THE 
C NONLINEAR CONTROLLED SOURCES. 
C 
      COMMON /TABINF/ IELMNT,ISBCKT,NSBCKT,IUNSAT,NUNSAT,ITEMPS,NUMTEM, 
     1   ISENS,NSENS,IFOUR,NFOUR,IFIELD,ICODE,IDELIM,ICOLUM,INSIZE, 
     2   JUNODE,LSBKPT,NUMBKP,IORDER,JMNODE,IUR,IUC,ILC,ILR,NUMOFF,ISR, 
     3   NMOFFC,ISEQ,ISEQ1,NEQN,NODEVS,NDIAG,ISWAP,IEQUA,MACINS,LVNIM1, 
     4   LX0,LVN,LYNL,LYU,LYL,LX1,LX2,LX3,LX4,LX5,LX6,LX7,LD0,LD1,LTD,
     5   IMYNL,IMVN,LCVN,NSNOD,NSMAT,NSVAL,ICNOD,ICMAT,ICVAL, 
     6   LOUTPT,LPOL,LZER,IRSWPF,IRSWPR,ICSWPF,ICSWPR,IRPT,JCPT,
     7   IROWNO,JCOLNO,NTTBR,NTTAR,LVNTMP 
      COMMON /STATUS/ OMEGA,TIME,DELTA,DELOLD(7),AG(7),VT,XNI,EGFET,
     1   XMU,MODE,MODEDC,ICALC,INITF,METHOD,IORD,MAXORD,NONCON,ITERNO,
     2   ITEMNO,NOSOLV,MODAC,IPIV,IVMFLG,IPOSTP,ISCRCH,IOFILE 
      COMMON /BLANK/ VALUE(64)
      INTEGER NODPLC(64)
      COMPLEX CVALUE(32)
      EQUIVALENCE (VALUE(1),NODPLC(1),CVALUE(1))
C 
      GO TO (40,10,40,20,30,50), INITF
C 
   10 VNEW=VINIT
      GO TO 70
   20 VNEW=VALUE(LX0+LOCT)
      GO TO 70
   30 VNEW=VALUE(LX1+LOCT)
      GO TO 70
   40 VNEW=VALUE(LVNIM1+NODE1)-VALUE(LVNIM1+NODE2)
      GO TO 60
   50 CALL COPY8(VALUE(LX1+LOCT),VALUE(LX0+LOCT),NUPDA) 
      XFACT=DELTA/DELOLD(2) 
      VNEW=(1.0D0+XFACT)*VALUE(LX1+LOCT)-XFACT*VALUE(LX2+LOCT)
   60 IF (DABS(VNEW).LE.1.0D0) GO TO 80 
      DELV=VNEW-VALUE(LX0+LOCT) 
      IF (DABS(DELV).LE.0.1D0) GO TO 80 
      VLIM=DMAX1(DABS(0.1D0*VALUE(LX0+LOCT)),0.1D0) 
      VNEW=VALUE(LX0+LOCT)+DSIGN(DMIN1(DABS(DELV),VLIM),DELV) 
      GO TO 70
   70 ICHECK=1
   80 VALUE(LX0+LOCT)=VNEW
      RETURN
      END 
