/* jfet.f -- translated by f2c (version of 3 February 1990  3:36:42).
   You must link the resulting object file with the libraries:
	-lF77 -lI77 -lm -lc   (in that order)
*/

#include "f2c.h"

/* Common Block Declarations */

struct {
    integer ielmnt, isbckt, nsbckt, iunsat, nunsat, itemps, numtem, isens, 
	    nsens, ifour, nfour, ifield, icode, idelim, icolum, insize, 
	    junode, lsbkpt, numbkp, iorder, jmnode, iur, iuc, ilc, ilr, 
	    numoff, isr, nmoffc, iseq, iseq1, neqn, nodevs, ndiag, iswap, 
	    iequa, macins, lvnim1, lx0, lvn, lynl, lyu, lyl, lx1, lx2, lx3, 
	    lx4, lx5, lx6, lx7, ld0, ld1, ltd, imynl, imvn, lcvn, nsnod, 
	    nsmat, nsval, icnod, icmat, icval, loutpt, lpol, lzer, irswpf, 
	    irswpr, icswpf, icswpr, irpt, jcpt, irowno, jcolno, nttbr, nttar, 
	    lvntmp;
} tabinf_;

#define tabinf_1 tabinf_

struct {
    integer locate[50], jelcnt[50], nunods, ncnods, numnod, nstop, nut, nlt, 
	    nxtrm, ndist, ntlin, ibr, numvs, numalt, numcyc;
} cirdat_;

#define cirdat_1 cirdat_

struct {
    doublereal omega, time, delta, delold[7], ag[7], vt, xni, egfet, xmu, 
	    sfactr;
    integer mode, modedc, icalc, initf, method, iord, maxord, noncon, iterno, 
	    itemno, nosolv, modac, ipiv, ivmflg, ipostp, iscrch, iofile;
} status_;

#define status_1 status_

struct {
    doublereal twopi, xlog2, xlog10, root2, rad, boltz, charge, ctok, gmin, 
	    reltol, abstol, vntol, trtol, chgtol, eps0, epssil, epsox, pivtol,
	     pivrel;
} knstnt_;

#define knstnt_1 knstnt_

struct {
    doublereal value[200000];
} blank_;

#define blank_1 blank_

/*<       subroutine jfet >*/
/* Subroutine */ int jfet_()
{
    /* System generated locals */
    integer i_1;
    doublereal d_1, d_2, d_3;

    /* Builtin functions */
    double exp(), sqrt();

    /* Local variables */
    static doublereal area, fcpb, beta;
#define cgdo ((doublereal *)&blank_1 + 4)
#define cqgd ((doublereal *)&blank_1 + 12)
#define ggdo ((doublereal *)&blank_1 + 8)
    static doublereal phib;
    static integer ioff;
    static doublereal evgd, czgd;
    static integer locm;
#define gdso ((doublereal *)&blank_1 + 6)
#define vgdo ((doublereal *)&blank_1 + 1)
#define ggso ((doublereal *)&blank_1 + 7)
#define cqgs ((doublereal *)&blank_1 + 10)
    static integer locv, loct;
    static doublereal gdpr;
#define vgso ((doublereal *)&blank_1)
    static doublereal type, gspr, csat, evgs, vgst, twob, vgdt, czgs, twop, 
	    sarg, fcpb2;
    static integer locy, ichk1, node1, node2, node3, node4, node5;
    static doublereal capgd, cdhat, cghat, ceqgd, betap, capgs, cdreq, ceqgs, 
	    xlamb, xfact, vcrit, f1, f2, f3, czgdf2, czgsf2, cd, cg;
    extern /* Subroutine */ int intgr8_();
    static doublereal gm;
    static integer icheck;
    static doublereal cdrain, delvgd;
#define nodplc ((integer *)&blank_1)
#define cvalue ((complex *)&blank_1)
    static doublereal delvgs, delvds;
    extern /* Subroutine */ int pnjlim_(), fetlim_();
    static doublereal cgd, ggd;
#define cdo ((doublereal *)&blank_1 + 3)
#define cgo ((doublereal *)&blank_1 + 2)
    static doublereal ceq;
#define qgd ((doublereal *)&blank_1 + 11)
    static integer loc;
    static doublereal gds, geq, vgd;
#define gmo ((doublereal *)&blank_1 + 5)
    static doublereal ggs;
#define qgs ((doublereal *)&blank_1 + 9)
    static doublereal vds, vgs, tol, vto;

/*<       implicit double precision (a-h,o-z) >*/

/*     this routine processes jfets for dc and transient analyses. */

/* 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 >*/
/* 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 >*/
/* 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 >*/
/* 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 >*/
/* spice version 2g.6  sccsid=blank 3/15/83 */
/*<       common /blank/ value(200000) >*/
/*<       integer nodplc(64) >*/
/*<       complex cvalue(32) >*/
/*<       equivalence (value(1),nodplc(1),cvalue(1)) >*/


/*<       dimension vgso(1),vgdo(1),cgo(1),cdo(1),cgdo(1),gmo(1),gdso(1), >*/
/*<      1   ggso(1),ggdo(1),qgs(1),cqgs(1),qgd(1),cqgd(1) >*/
/*<       equivalence (vgso(1),value( 1)),(vgdo(1),value( 2)), >*/
/*<      1            (cgo (1),value( 3)),(cdo (1),value( 4)), >*/
/*<      2            (cgdo(1),value( 5)),(gmo (1),value( 6)), >*/
/*<      3            (gdso(1),value( 7)),(ggso(1),value( 8)), >*/
/*<      4            (ggdo(1),value( 9)),(qgs (1),value(10)), >*/
/*<      5            (cqgs(1),value(11)),(qgd (1),value(12)), >*/
/*<      6            (cqgd(1),value(13)) >*/


/*<       loc=locate(13) >*/
    loc = cirdat_1.locate[12];
/*<    10 if ((loc.eq.0).or.(nodplc(loc+25).ne.0)) return >*/
L10:
    if (loc == 0 || nodplc[loc + 24] != 0) {
	return 0;
    }
/*<       locv=nodplc(loc+1) >*/
    locv = nodplc[loc];
/*<       node1=nodplc(loc+2) >*/
    node1 = nodplc[loc + 1];
/*<       node2=nodplc(loc+3) >*/
    node2 = nodplc[loc + 2];
/*<       node3=nodplc(loc+4) >*/
    node3 = nodplc[loc + 3];
/*<       node4=nodplc(loc+5) >*/
    node4 = nodplc[loc + 4];
/*<       node5=nodplc(loc+6) >*/
    node5 = nodplc[loc + 5];
/*<       locm=nodplc(loc+7) >*/
    locm = nodplc[loc + 6];
/*<       ioff=nodplc(loc+8) >*/
    ioff = nodplc[loc + 7];
/*<       type=nodplc(locm+2) >*/
    type = (doublereal) nodplc[locm + 1];
/*<       locm=nodplc(locm+1) >*/
    locm = nodplc[locm];
/*<       loct=nodplc(loc+19) >*/
    loct = nodplc[loc + 18];

/*  dc model parameters */

/*<       area=value(locv+1) >*/
    area = blank_1.value[locv];
/*<       vto=value(locm+1) >*/
    vto = blank_1.value[locm];
/*<       beta=value(locm+2)*area >*/
    beta = blank_1.value[locm + 1] * area;
/*<       xlamb=value(locm+3) >*/
    xlamb = blank_1.value[locm + 2];
/*<       gdpr=value(locm+4)*area >*/
    gdpr = blank_1.value[locm + 3] * area;
/*<       gspr=value(locm+5)*area >*/
    gspr = blank_1.value[locm + 4] * area;
/*<       csat=value(locm+9)*area >*/
    csat = blank_1.value[locm + 8] * area;
/*<       vcrit=value(locm+16) >*/
    vcrit = blank_1.value[locm + 15];

/*  initialization */

/*<       icheck=1 >*/
    icheck = 1;
/*<       go to (100,20,30,50,60,70), initf >*/
    switch (status_1.initf) {
	case 1:  goto L100;
	case 2:  goto L20;
	case 3:  goto L30;
	case 4:  goto L50;
	case 5:  goto L60;
	case 6:  goto L70;
    }
/*<    20 if(mode.ne.1.or.modedc.ne.2.or.nosolv.eq.0) go to 25 >*/
L20:
    if (status_1.mode != 1 || status_1.modedc != 2 || status_1.nosolv == 0) {
	goto L25;
    }
/*<       vds=type*value(locv+2) >*/
    vds = type * blank_1.value[locv + 1];
/*<       vgs=type*value(locv+3) >*/
    vgs = type * blank_1.value[locv + 2];
/*<       vgd=vgs-vds >*/
    vgd = vgs - vds;
/*<       go to 300 >*/
    goto L300;
/*<    25 if(ioff.ne.0) go to 40 >*/
L25:
    if (ioff != 0) {
	goto L40;
    }
/*<       vgs=-1.0d0 >*/
    vgs = -1.;
/*<       vgd=-1.0d0 >*/
    vgd = -1.;
/*<       go to 300 >*/
    goto L300;
/*<    30 if (ioff.eq.0) go to 100 >*/
L30:
    if (ioff == 0) {
	goto L100;
    }
/*<    40 vgs=0.0d0 >*/
L40:
    vgs = 0.;
/*<       vgd=0.0d0 >*/
    vgd = 0.;
/*<       go to 300 >*/
    goto L300;
/*<    50 vgs=vgso(lx0+loct) >*/
L50:
    vgs = vgso[tabinf_1.lx0 + loct - 1];
/*<       vgd=vgdo(lx0+loct) >*/
    vgd = vgdo[tabinf_1.lx0 + loct - 1];
/*<       go to 300 >*/
    goto L300;
/*<    60 vgs=vgso(lx1+loct) >*/
L60:
    vgs = vgso[tabinf_1.lx1 + loct - 1];
/*<       vgd=vgdo(lx1+loct) >*/
    vgd = vgdo[tabinf_1.lx1 + loct - 1];
/*<       go to 300 >*/
    goto L300;
/*<    70 xfact=delta/delold(2) >*/
L70:
    xfact = status_1.delta / status_1.delold[1];
/*<       vgso(lx0+loct)=vgso(lx1+loct) >*/
    vgso[tabinf_1.lx0 + loct - 1] = vgso[tabinf_1.lx1 + loct - 1];
/*<       vgs=(1.0d0+xfact)*vgso(lx1+loct)-xfact*vgso(lx2+loct) >*/
    vgs = (xfact + 1.) * vgso[tabinf_1.lx1 + loct - 1] - xfact * vgso[
	    tabinf_1.lx2 + loct - 1];
/*<       vgdo(lx0+loct)=vgdo(lx1+loct) >*/
    vgdo[tabinf_1.lx0 + loct - 1] = vgdo[tabinf_1.lx1 + loct - 1];
/*<       vgd=(1.0d0+xfact)*vgdo(lx1+loct)-xfact*vgdo(lx2+loct) >*/
    vgd = (xfact + 1.) * vgdo[tabinf_1.lx1 + loct - 1] - xfact * vgdo[
	    tabinf_1.lx2 + loct - 1];
/*<       cgo(lx0+loct)=cgo(lx1+loct) >*/
    cgo[tabinf_1.lx0 + loct - 1] = cgo[tabinf_1.lx1 + loct - 1];
/*<       cdo(lx0+loct)=cdo(lx1+loct) >*/
    cdo[tabinf_1.lx0 + loct - 1] = cdo[tabinf_1.lx1 + loct - 1];
/*<       cgdo(lx0+loct)=cgdo(lx1+loct) >*/
    cgdo[tabinf_1.lx0 + loct - 1] = cgdo[tabinf_1.lx1 + loct - 1];
/*<       gmo(lx0+loct)=gmo(lx1+loct) >*/
    gmo[tabinf_1.lx0 + loct - 1] = gmo[tabinf_1.lx1 + loct - 1];
/*<       gdso(lx0+loct)=gdso(lx1+loct) >*/
    gdso[tabinf_1.lx0 + loct - 1] = gdso[tabinf_1.lx1 + loct - 1];
/*<       ggso(lx0+loct)=ggso(lx1+loct) >*/
    ggso[tabinf_1.lx0 + loct - 1] = ggso[tabinf_1.lx1 + loct - 1];
/*<       ggdo(lx0+loct)=ggdo(lx1+loct) >*/
    ggdo[tabinf_1.lx0 + loct - 1] = ggdo[tabinf_1.lx1 + loct - 1];
/*<       go to 110 >*/
    goto L110;

/*  compute new nonlinear branch voltages */

/*<   100 vgs=type*(value(lvnim1+node2)-value(lvnim1+node5)) >*/
L100:
    vgs = type * (blank_1.value[tabinf_1.lvnim1 + node2 - 1] - blank_1.value[
	    tabinf_1.lvnim1 + node5 - 1]);
/*<       vgd=type*(value(lvnim1+node2)-value(lvnim1+node4)) >*/
    vgd = type * (blank_1.value[tabinf_1.lvnim1 + node2 - 1] - blank_1.value[
	    tabinf_1.lvnim1 + node4 - 1]);
/*<   110 delvgs=vgs-vgso(lx0+loct) >*/
L110:
    delvgs = vgs - vgso[tabinf_1.lx0 + loct - 1];
/*<       delvgd=vgd-vgdo(lx0+loct) >*/
    delvgd = vgd - vgdo[tabinf_1.lx0 + loct - 1];
/*<       delvds=delvgs-delvgd >*/
    delvds = delvgs - delvgd;
/*<       cghat=cgo(lx0+loct)+ggdo(lx0+loct)*delvgd+ggso(lx0+loct)*delvgs >*/
    cghat = cgo[tabinf_1.lx0 + loct - 1] + ggdo[tabinf_1.lx0 + loct - 1] * 
	    delvgd + ggso[tabinf_1.lx0 + loct - 1] * delvgs;
/*<       cdhat=cdo(lx0+loct)+gmo(lx0+loct)*delvgs+gdso(lx0+loct)*delvds >*/
/*<      1   -ggdo(lx0+loct)*delvgd >*/
    cdhat = cdo[tabinf_1.lx0 + loct - 1] + gmo[tabinf_1.lx0 + loct - 1] * 
	    delvgs + gdso[tabinf_1.lx0 + loct - 1] * delvds - ggdo[
	    tabinf_1.lx0 + loct - 1] * delvgd;

/*  bypass if solution has not changed */

/*<       if (initf.eq.6) go to 200 >*/
    if (status_1.initf == 6) {
	goto L200;
    }
/*<       tol=reltol*dmax1(dabs(vgs),dabs(vgso(lx0+loct)))+vntol >*/
/* Computing MAX */
    d_2 = abs(vgs), d_3 = (d_1 = vgso[tabinf_1.lx0 + loct - 1], abs(d_1));
    tol = knstnt_1.reltol * max(d_3,d_2) + knstnt_1.vntol;
/*<       if (dabs(delvgs).ge.tol) go to 200 >*/
    if (abs(delvgs) >= tol) {
	goto L200;
    }
/*<       tol=reltol*dmax1(dabs(vgd),dabs(vgdo(lx0+loct)))+vntol >*/
/* Computing MAX */
    d_2 = abs(vgd), d_3 = (d_1 = vgdo[tabinf_1.lx0 + loct - 1], abs(d_1));
    tol = knstnt_1.reltol * max(d_3,d_2) + knstnt_1.vntol;
/*<       if (dabs(delvgd).ge.tol) go to 200 >*/
    if (abs(delvgd) >= tol) {
	goto L200;
    }
/*<       tol=reltol*dmax1(dabs(cghat),dabs(cgo(lx0+loct)))+abstol >*/
/* Computing MAX */
    d_2 = abs(cghat), d_3 = (d_1 = cgo[tabinf_1.lx0 + loct - 1], abs(d_1));
    tol = knstnt_1.reltol * max(d_3,d_2) + knstnt_1.abstol;
/*<       if (dabs(cghat-cgo(lx0+loct)).ge.tol) go to 200 >*/
    if ((d_1 = cghat - cgo[tabinf_1.lx0 + loct - 1], abs(d_1)) >= tol) {
	goto L200;
    }
/*<       tol=reltol*dmax1(dabs(cdhat),dabs(cdo(lx0+loct)))+abstol >*/
/* Computing MAX */
    d_2 = abs(cdhat), d_3 = (d_1 = cdo[tabinf_1.lx0 + loct - 1], abs(d_1));
    tol = knstnt_1.reltol * max(d_3,d_2) + knstnt_1.abstol;
/*<       if (dabs(cdhat-cdo(lx0+loct)).ge.tol) go to 200 >*/
    if ((d_1 = cdhat - cdo[tabinf_1.lx0 + loct - 1], abs(d_1)) >= tol) {
	goto L200;
    }
/*<       vgs=vgso(lx0+loct) >*/
    vgs = vgso[tabinf_1.lx0 + loct - 1];
/*<       vgd=vgdo(lx0+loct) >*/
    vgd = vgdo[tabinf_1.lx0 + loct - 1];
/*<       vds=vgs-vgd >*/
    vds = vgs - vgd;
/*<       cg=cgo(lx0+loct) >*/
    cg = cgo[tabinf_1.lx0 + loct - 1];
/*<       cd=cdo(lx0+loct) >*/
    cd = cdo[tabinf_1.lx0 + loct - 1];
/*<       cgd=cgdo(lx0+loct) >*/
    cgd = cgdo[tabinf_1.lx0 + loct - 1];
/*<       gm=gmo(lx0+loct) >*/
    gm = gmo[tabinf_1.lx0 + loct - 1];
/*<       gds=gdso(lx0+loct) >*/
    gds = gdso[tabinf_1.lx0 + loct - 1];
/*<       ggs=ggso(lx0+loct) >*/
    ggs = ggso[tabinf_1.lx0 + loct - 1];
/*<       ggd=ggdo(lx0+loct) >*/
    ggd = ggdo[tabinf_1.lx0 + loct - 1];
/*<       go to 900 >*/
    goto L900;

/*  limit nonlinear branch voltages */

/*<   200 ichk1=1 >*/
L200:
    ichk1 = 1;
/*<       call pnjlim(vgs,vgso(lx0+loct),vt,vcrit,icheck) >*/
    pnjlim_(&vgs, &vgso[tabinf_1.lx0 + loct - 1], &status_1.vt, &vcrit, &
	    icheck);
/*<       call pnjlim(vgd,vgdo(lx0+loct),vt,vcrit,ichk1) >*/
    pnjlim_(&vgd, &vgdo[tabinf_1.lx0 + loct - 1], &status_1.vt, &vcrit, &
	    ichk1);
/*<       if (ichk1.eq.1) icheck=1 >*/
    if (ichk1 == 1) {
	icheck = 1;
    }
/*<       call fetlim(vgs,vgso(lx0+loct),vto) >*/
    fetlim_(&vgs, &vgso[tabinf_1.lx0 + loct - 1], &vto);
/*<       call fetlim(vgd,vgdo(lx0+loct),vto) >*/
    fetlim_(&vgd, &vgdo[tabinf_1.lx0 + loct - 1], &vto);

/*  determine dc current and derivatives */

/*<   300 vds=vgs-vgd >*/
L300:
    vds = vgs - vgd;
/*<       if (vgs.gt.-5.0d0*vt) go to 310 >*/
    if (vgs > status_1.vt * -5.) {
	goto L310;
    }
/*<       ggs=-csat/vgs+gmin >*/
    ggs = -csat / vgs + knstnt_1.gmin;
/*<       cg=ggs*vgs >*/
    cg = ggs * vgs;
/*<       go to 320 >*/
    goto L320;
/*<   310 evgs=dexp(vgs/vt) >*/
L310:
    evgs = exp(vgs / status_1.vt);
/*<       ggs=csat*evgs/vt+gmin >*/
    ggs = csat * evgs / status_1.vt + knstnt_1.gmin;
/*<       cg=csat*(evgs-1.0d0)+gmin*vgs >*/
    cg = csat * (evgs - 1.) + knstnt_1.gmin * vgs;
/*<   320 if (vgd.gt.-5.0d0*vt) go to 330 >*/
L320:
    if (vgd > status_1.vt * -5.) {
	goto L330;
    }
/*<       ggd=-csat/vgd+gmin >*/
    ggd = -csat / vgd + knstnt_1.gmin;
/*<       cgd=ggd*vgd >*/
    cgd = ggd * vgd;
/*<       go to 340 >*/
    goto L340;
/*<   330 evgd=dexp(vgd/vt) >*/
L330:
    evgd = exp(vgd / status_1.vt);
/*<       ggd=csat*evgd/vt+gmin >*/
    ggd = csat * evgd / status_1.vt + knstnt_1.gmin;
/*<       cgd=csat*(evgd-1.0d0)+gmin*vgd >*/
    cgd = csat * (evgd - 1.) + knstnt_1.gmin * vgd;
/*<   340 cg=cg+cgd >*/
L340:
    cg += cgd;

/*  compute drain current and derivitives for normal mode */

/*<   400 if (vds.lt.0.0d0) go to 450 >*/
/* L400: */
    if (vds < 0.) {
	goto L450;
    }
/*<       vgst=vgs-vto >*/
    vgst = vgs - vto;

/*  normal mode, cutoff region */

/*<       if (vgst.gt.0.0d0) go to 410 >*/
    if (vgst > 0.) {
	goto L410;
    }
/*<       cdrain=0.0d0 >*/
    cdrain = 0.;
/*<       gm=0.0d0 >*/
    gm = 0.;
/*<       gds=0.0d0 >*/
    gds = 0.;
/*<       go to 490 >*/
    goto L490;

/*  normal mode, saturation region */

/*<   410 betap=beta*(1.0d0+xlamb*vds) >*/
L410:
    betap = beta * (xlamb * vds + 1.);
/*<       twob=betap+betap >*/
    twob = betap + betap;
/*<       if (vgst.gt.vds) go to 420 >*/
    if (vgst > vds) {
	goto L420;
    }
/*<       cdrain=betap*vgst*vgst >*/
    cdrain = betap * vgst * vgst;
/*<       gm=twob*vgst >*/
    gm = twob * vgst;
/*<       gds=xlamb*beta*vgst*vgst >*/
    gds = xlamb * beta * vgst * vgst;
/*<       go to 490 >*/
    goto L490;

/*  normal mode, linear region */

/*<   420 cdrain=betap*vds*(vgst+vgst-vds) >*/
L420:
    cdrain = betap * vds * (vgst + vgst - vds);
/*<       gm=twob*vds >*/
    gm = twob * vds;
/*<       gds=twob*(vgst-vds)+xlamb*beta*vds*(vgst+vgst-vds) >*/
    gds = twob * (vgst - vds) + xlamb * beta * vds * (vgst + vgst - vds);
/*<       go to 490 >*/
    goto L490;

/*  compute drain current and derivitives for inverse mode */

/*<   450 vgdt=vgd-vto >*/
L450:
    vgdt = vgd - vto;

/*  inverse mode, cutoff region */

/*<       if (vgdt.gt.0.0d0) go to 460 >*/
    if (vgdt > 0.) {
	goto L460;
    }
/*<       cdrain=0.0d0 >*/
    cdrain = 0.;
/*<       gm=0.0d0 >*/
    gm = 0.;
/*<       gds=0.0d0 >*/
    gds = 0.;
/*<       go to 490 >*/
    goto L490;

/*  inverse mode, saturation region */

/*<   460 betap=beta*(1.0d0-xlamb*vds) >*/
L460:
    betap = beta * (1. - xlamb * vds);
/*<       twob=betap+betap >*/
    twob = betap + betap;
/*<       if (vgdt.gt.-vds) go to 470 >*/
    if (vgdt > -vds) {
	goto L470;
    }
/*<       cdrain=-betap*vgdt*vgdt >*/
    cdrain = -betap * vgdt * vgdt;
/*<       gm=-twob*vgdt >*/
    gm = -twob * vgdt;
/*<       gds=xlamb*beta*vgdt*vgdt-gm >*/
    gds = xlamb * beta * vgdt * vgdt - gm;
/*<       go to 490 >*/
    goto L490;

/*  inverse mode, linear region */

/*<   470 cdrain=betap*vds*(vgdt+vgdt+vds) >*/
L470:
    cdrain = betap * vds * (vgdt + vgdt + vds);
/*<       gm=twob*vds >*/
    gm = twob * vds;
/*<       gds=twob*vgdt-xlamb*beta*vds*(vgdt+vgdt+vds) >*/
    gds = twob * vgdt - xlamb * beta * vds * (vgdt + vgdt + vds);

/*  compute equivalent drain current source */

/*<   490 cd=cdrain-cgd >*/
L490:
    cd = cdrain - cgd;
/*<       if (mode.ne.1) go to 500 >*/
    if (status_1.mode != 1) {
	goto L500;
    }
/*<       if ((modedc.eq.2).and.(nosolv.ne.0)) go to 500 >*/
    if (status_1.modedc == 2 && status_1.nosolv != 0) {
	goto L500;
    }
/*<       if (initf.eq.4) go to 500 >*/
    if (status_1.initf == 4) {
	goto L500;
    }
/*<       go to 700 >*/
    goto L700;

/*  charge storage elements */

/*<   500 czgs=value(locm+6)*area >*/
L500:
    czgs = blank_1.value[locm + 5] * area;
/*<       czgd=value(locm+7)*area >*/
    czgd = blank_1.value[locm + 6] * area;
/*<       phib=value(locm+8) >*/
    phib = blank_1.value[locm + 7];
/*<       twop=phib+phib >*/
    twop = phib + phib;
/*<       fcpb=value(locm+12) >*/
    fcpb = blank_1.value[locm + 11];
/*<       fcpb2=fcpb*fcpb >*/
    fcpb2 = fcpb * fcpb;
/*<       f1=value(locm+13) >*/
    f1 = blank_1.value[locm + 12];
/*<       f2=value(locm+14) >*/
    f2 = blank_1.value[locm + 13];
/*<       f3=value(locm+15) >*/
    f3 = blank_1.value[locm + 14];
/*<       czgsf2=czgs/f2 >*/
    czgsf2 = czgs / f2;
/*<       czgdf2=czgd/f2 >*/
    czgdf2 = czgd / f2;
/*<       if (vgs.ge.fcpb) go to 510 >*/
    if (vgs >= fcpb) {
	goto L510;
    }
/*<       sarg=dsqrt(1.0d0-vgs/phib) >*/
    sarg = sqrt(1. - vgs / phib);
/*<       qgs(lx0+loct)=twop*czgs*(1.0d0-sarg) >*/
    qgs[tabinf_1.lx0 + loct - 1] = twop * czgs * (1. - sarg);
/*<       capgs=czgs/sarg >*/
    capgs = czgs / sarg;
/*<       go to 520 >*/
    goto L520;
/*<   510 qgs(lx0+loct)=czgs*f1+czgsf2*(f3*(vgs-fcpb) >*/
/*<      1   +(vgs*vgs-fcpb2)/(twop+twop)) >*/
L510:
    qgs[tabinf_1.lx0 + loct - 1] = czgs * f1 + czgsf2 * (f3 * (vgs - fcpb) + (
	    vgs * vgs - fcpb2) / (twop + twop));
/*<       capgs=czgsf2*(f3+vgs/twop) >*/
    capgs = czgsf2 * (f3 + vgs / twop);
/*<   520 if (vgd.ge.fcpb) go to 530 >*/
L520:
    if (vgd >= fcpb) {
	goto L530;
    }
/*<       sarg=dsqrt(1.0d0-vgd/phib) >*/
    sarg = sqrt(1. - vgd / phib);
/*<       qgd(lx0+loct)=twop*czgd*(1.0d0-sarg) >*/
    qgd[tabinf_1.lx0 + loct - 1] = twop * czgd * (1. - sarg);
/*<       capgd=czgd/sarg >*/
    capgd = czgd / sarg;
/*<       go to 560 >*/
    goto L560;
/*<   530 qgd(lx0+loct)=czgd*f1+czgdf2*(f3*(vgd-fcpb) >*/
/*<      1   +(vgd*vgd-fcpb2)/(twop+twop)) >*/
L530:
    qgd[tabinf_1.lx0 + loct - 1] = czgd * f1 + czgdf2 * (f3 * (vgd - fcpb) + (
	    vgd * vgd - fcpb2) / (twop + twop));
/*<       capgd=czgdf2*(f3+vgd/twop) >*/
    capgd = czgdf2 * (f3 + vgd / twop);

/*  store small-signal parameters */

/*<   560 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700 >*/
L560:
    if (status_1.mode == 1 && status_1.modedc == 2 && status_1.nosolv != 0) {
	goto L700;
    }
/*<       if (initf.ne.4) go to 600 >*/
    if (status_1.initf != 4) {
	goto L600;
    }
/*<       value(lx0+loct+9)=capgs >*/
    blank_1.value[tabinf_1.lx0 + loct + 8] = capgs;
/*<       value(lx0+loct+11)=capgd >*/
    blank_1.value[tabinf_1.lx0 + loct + 10] = capgd;
/*<       go to 1000 >*/
    goto L1000;

/*  transient analysis */

/*<   600 if (initf.ne.5) go to 610 >*/
L600:
    if (status_1.initf != 5) {
	goto L610;
    }
/*<       qgs(lx1+loct)=qgs(lx0+loct) >*/
    qgs[tabinf_1.lx1 + loct - 1] = qgs[tabinf_1.lx0 + loct - 1];
/*<       qgd(lx1+loct)=qgd(lx0+loct) >*/
    qgd[tabinf_1.lx1 + loct - 1] = qgd[tabinf_1.lx0 + loct - 1];
/*<   610 call intgr8(geq,ceq,capgs,loct+9) >*/
L610:
    i_1 = loct + 9;
    intgr8_(&geq, &ceq, &capgs, &i_1);
/*<       ggs=ggs+geq >*/
    ggs += geq;
/*<       cg=cg+cqgs(lx0+loct) >*/
    cg += cqgs[tabinf_1.lx0 + loct - 1];
/*<       call intgr8(geq,ceq,capgd,loct+11) >*/
    i_1 = loct + 11;
    intgr8_(&geq, &ceq, &capgd, &i_1);
/*<       ggd=ggd+geq >*/
    ggd += geq;
/*<       cg=cg+cqgd(lx0+loct) >*/
    cg += cqgd[tabinf_1.lx0 + loct - 1];
/*<       cd=cd-cqgd(lx0+loct) >*/
    cd -= cqgd[tabinf_1.lx0 + loct - 1];
/*<       cgd=cgd+cqgd(lx0+loct) >*/
    cgd += cqgd[tabinf_1.lx0 + loct - 1];
/*<       if (initf.ne.5) go to 700 >*/
    if (status_1.initf != 5) {
	goto L700;
    }
/*<       cqgs(lx1+loct)=cqgs(lx0+loct) >*/
    cqgs[tabinf_1.lx1 + loct - 1] = cqgs[tabinf_1.lx0 + loct - 1];
/*<       cqgd(lx1+loct)=cqgd(lx0+loct) >*/
    cqgd[tabinf_1.lx1 + loct - 1] = cqgd[tabinf_1.lx0 + loct - 1];

/*  check convergence */

/*<   700 if (initf.ne.3) go to 710 >*/
L700:
    if (status_1.initf != 3) {
	goto L710;
    }
/*<       if (ioff.eq.0) go to 710 >*/
    if (ioff == 0) {
	goto L710;
    }
/*<       go to 750 >*/
    goto L750;
/*<   710 if (icheck.eq.1) go to 720 >*/
L710:
    if (icheck == 1) {
	goto L720;
    }
/*<       tol=reltol*dmax1(dabs(cghat),dabs(cg))+abstol >*/
/* Computing MAX */
    d_1 = abs(cghat), d_2 = abs(cg);
    tol = knstnt_1.reltol * max(d_2,d_1) + knstnt_1.abstol;
/*<       if (dabs(cghat-cg).ge.tol) go to 720 >*/
    if ((d_1 = cghat - cg, abs(d_1)) >= tol) {
	goto L720;
    }
/*<       tol=reltol*dmax1(dabs(cdhat),dabs(cd))+abstol >*/
/* Computing MAX */
    d_1 = abs(cdhat), d_2 = abs(cd);
    tol = knstnt_1.reltol * max(d_2,d_1) + knstnt_1.abstol;
/*<       if (dabs(cdhat-cd).le.tol) go to 750 >*/
    if ((d_1 = cdhat - cd, abs(d_1)) <= tol) {
	goto L750;
    }
/*<   720 noncon=noncon+1 >*/
L720:
    ++status_1.noncon;
/*<   750 vgso(lx0+loct)=vgs >*/
L750:
    vgso[tabinf_1.lx0 + loct - 1] = vgs;
/*<       vgdo(lx0+loct)=vgd >*/
    vgdo[tabinf_1.lx0 + loct - 1] = vgd;
/*<       cgo(lx0+loct)=cg >*/
    cgo[tabinf_1.lx0 + loct - 1] = cg;
/*<       cdo(lx0+loct)=cd >*/
    cdo[tabinf_1.lx0 + loct - 1] = cd;
/*<       cgdo(lx0+loct)=cgd >*/
    cgdo[tabinf_1.lx0 + loct - 1] = cgd;
/*<       gmo(lx0+loct)=gm >*/
    gmo[tabinf_1.lx0 + loct - 1] = gm;
/*<       gdso(lx0+loct)=gds >*/
    gdso[tabinf_1.lx0 + loct - 1] = gds;
/*<       ggso(lx0+loct)=ggs >*/
    ggso[tabinf_1.lx0 + loct - 1] = ggs;
/*<       ggdo(lx0+loct)=ggd >*/
    ggdo[tabinf_1.lx0 + loct - 1] = ggd;

/*  load current vector */

/*<   900 ceqgd=type*(cgd-ggd*vgd) >*/
L900:
    ceqgd = type * (cgd - ggd * vgd);
/*<       ceqgs=type*((cg-cgd)-ggs*vgs) >*/
    ceqgs = type * (cg - cgd - ggs * vgs);
/*<       cdreq=type*((cd+cgd)-gds*vds-gm*vgs) >*/
    cdreq = type * (cd + cgd - gds * vds - gm * vgs);
/*<       value(lvn+node2)=value(lvn+node2)-ceqgs-ceqgd >*/
    blank_1.value[tabinf_1.lvn + node2 - 1] = blank_1.value[tabinf_1.lvn + 
	    node2 - 1] - ceqgs - ceqgd;
/*<       value(lvn+node4)=value(lvn+node4)-cdreq+ceqgd >*/
    blank_1.value[tabinf_1.lvn + node4 - 1] = blank_1.value[tabinf_1.lvn + 
	    node4 - 1] - cdreq + ceqgd;
/*<       value(lvn+node5)=value(lvn+node5)+cdreq+ceqgs >*/
    blank_1.value[tabinf_1.lvn + node5 - 1] = blank_1.value[tabinf_1.lvn + 
	    node5 - 1] + cdreq + ceqgs;

/*  load y matrix */

/*<       locy=lvn+nodplc(loc+20) >*/
    locy = tabinf_1.lvn + nodplc[loc + 19];
/*<       value(locy)=value(locy)+gdpr >*/
    blank_1.value[locy - 1] += gdpr;
/*<       locy=lvn+nodplc(loc+21) >*/
    locy = tabinf_1.lvn + nodplc[loc + 20];
/*<       value(locy)=value(locy)+ggd+ggs >*/
    blank_1.value[locy - 1] = blank_1.value[locy - 1] + ggd + ggs;
/*<       locy=lvn+nodplc(loc+22) >*/
    locy = tabinf_1.lvn + nodplc[loc + 21];
/*<       value(locy)=value(locy)+gspr >*/
    blank_1.value[locy - 1] += gspr;
/*<       locy=lvn+nodplc(loc+23) >*/
    locy = tabinf_1.lvn + nodplc[loc + 22];
/*<       value(locy)=value(locy)+gdpr+gds+ggd >*/
    blank_1.value[locy - 1] = blank_1.value[locy - 1] + gdpr + gds + ggd;
/*<       locy=lvn+nodplc(loc+24) >*/
    locy = tabinf_1.lvn + nodplc[loc + 23];
/*<       value(locy)=value(locy)+gspr+gds+gm+ggs >*/
    blank_1.value[locy - 1] = blank_1.value[locy - 1] + gspr + gds + gm + ggs;

/*<       locy=lvn+nodplc(loc+9) >*/
    locy = tabinf_1.lvn + nodplc[loc + 8];
/*<       value(locy)=value(locy)-gdpr >*/
    blank_1.value[locy - 1] -= gdpr;
/*<       locy=lvn+nodplc(loc+10) >*/
    locy = tabinf_1.lvn + nodplc[loc + 9];
/*<       value(locy)=value(locy)-ggd >*/
    blank_1.value[locy - 1] -= ggd;
/*<       locy=lvn+nodplc(loc+11) >*/
    locy = tabinf_1.lvn + nodplc[loc + 10];
/*<       value(locy)=value(locy)-ggs >*/
    blank_1.value[locy - 1] -= ggs;
/*<       locy=lvn+nodplc(loc+12) >*/
    locy = tabinf_1.lvn + nodplc[loc + 11];
/*<       value(locy)=value(locy)-gspr >*/
    blank_1.value[locy - 1] -= gspr;
/*<       locy=lvn+nodplc(loc+13) >*/
    locy = tabinf_1.lvn + nodplc[loc + 12];
/*<       value(locy)=value(locy)-gdpr >*/
    blank_1.value[locy - 1] -= gdpr;
/*<       locy=lvn+nodplc(loc+14) >*/
    locy = tabinf_1.lvn + nodplc[loc + 13];
/*<       value(locy)=value(locy)+gm-ggd >*/
    blank_1.value[locy - 1] = blank_1.value[locy - 1] + gm - ggd;
/*<       locy=lvn+nodplc(loc+15) >*/
    locy = tabinf_1.lvn + nodplc[loc + 14];
/*<       value(locy)=value(locy)-gds-gm >*/
    blank_1.value[locy - 1] = blank_1.value[locy - 1] - gds - gm;
/*<       locy=lvn+nodplc(loc+16) >*/
    locy = tabinf_1.lvn + nodplc[loc + 15];
/*<       value(locy)=value(locy)-ggs-gm >*/
    blank_1.value[locy - 1] = blank_1.value[locy - 1] - ggs - gm;
/*<       locy=lvn+nodplc(loc+17) >*/
    locy = tabinf_1.lvn + nodplc[loc + 16];
/*<       value(locy)=value(locy)-gspr >*/
    blank_1.value[locy - 1] -= gspr;
/*<       locy=lvn+nodplc(loc+18) >*/
    locy = tabinf_1.lvn + nodplc[loc + 17];
/*<       value(locy)=value(locy)-gds >*/
    blank_1.value[locy - 1] -= gds;
/*<  1000 loc=nodplc(loc) >*/
L1000:
    loc = nodplc[loc - 1];
/*<       go to 10 >*/
    goto L10;
/*<       end >*/
} /* jfet_ */

#undef qgs
#undef gmo
#undef qgd
#undef cgo
#undef cdo
#undef cvalue
#undef nodplc
#undef vgso
#undef cqgs
#undef ggso
#undef vgdo
#undef gdso
#undef ggdo
#undef cqgd
#undef cgdo


