/* diode.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 diode >*/
/* Subroutine */ int diode_()
{
    /* System generated locals */
    integer i_1;
    doublereal d_1, d_2, d_3;

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

    /* Local variables */
    static doublereal capd, area, fcpb, cdeq;
    static integer ioff, locm;
    static doublereal csat, sarg;
    static integer locv, loct;
    static doublereal vlim, gspr;
    static integer locy, node1, node2, node3;
    static doublereal czof2, cdhat, delvd, xfact, vcrit, evrev, czero, f1, f2,
	     f3, cd, gd;
    extern /* Subroutine */ int intgr8_();
    static doublereal pb;
#define qd ((doublereal *)&blank_1 + 3)
    static integer icheck;
    static doublereal bv, vd, xm;
#define nodplc ((integer *)&blank_1)
#define cvalue ((complex *)&blank_1)
    extern /* Subroutine */ int pnjlim_();
    static doublereal vdtemp;
#define cdo ((doublereal *)&blank_1 + 1)
#define cqd ((doublereal *)&blank_1 + 4)
#define gdo ((doublereal *)&blank_1 + 2)
    static doublereal arg, ceq;
    static integer loc;
    static doublereal evd, geq;
#define vdo ((doublereal *)&blank_1)
    static doublereal tau, vte, tol;

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

/*     this routine processes diodes 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 vdo(1),cdo(1),gdo(1),qd(1),cqd(1) >*/
/*<       equivalence (vdo(1),value(1)),(cdo(1),value(2)), >*/
/*<      1   (gdo(1),value(3)),(qd(1),value(4)),(cqd(1),value(5)) >*/


/*<       loc=locate(11) >*/
    loc = cirdat_1.locate[10];
/*<    10 if ((loc.eq.0).or.(nodplc(loc+16).ne.0)) return >*/
L10:
    if (loc == 0 || nodplc[loc + 15] != 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];
/*<       locm=nodplc(loc+5) >*/
    locm = nodplc[loc + 4];
/*<       ioff=nodplc(loc+6) >*/
    ioff = nodplc[loc + 5];
/*<       locm=nodplc(locm+1) >*/
    locm = nodplc[locm];
/*<       loct=nodplc(loc+11) >*/
    loct = nodplc[loc + 10];

/*  dc model parameters */

/*<       area=value(locv+1) >*/
    area = blank_1.value[locv];
/*<       csat=value(locm+1)*area >*/
    csat = blank_1.value[locm] * area;
/*<       gspr=value(locm+2)*area >*/
    gspr = blank_1.value[locm + 1] * area;
/*<       vte=value(locm+3)*vt >*/
    vte = blank_1.value[locm + 2] * status_1.vt;
/*<       bv=value(locm+13) >*/
    bv = blank_1.value[locm + 12];
/*<       vcrit=value(locm+18) >*/
    vcrit = blank_1.value[locm + 17];

/*  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;
    }
/*<       vd=value(locv+2) >*/
    vd = blank_1.value[locv + 1];
/*<       go to 300 >*/
    goto L300;
/*<    25 if(ioff.ne.0) go to 40 >*/
L25:
    if (ioff != 0) {
	goto L40;
    }
/*<       vd=vcrit >*/
    vd = vcrit;
/*<       go to 300 >*/
    goto L300;
/*<    30 if (ioff.eq.0) go to 100 >*/
L30:
    if (ioff == 0) {
	goto L100;
    }
/*<    40 vd=0.0d0 >*/
L40:
    vd = 0.;
/*<       go to 300 >*/
    goto L300;
/*<    50 vd=vdo(lx0+loct) >*/
L50:
    vd = vdo[tabinf_1.lx0 + loct - 1];
/*<       go to 300 >*/
    goto L300;
/*<    60 vd=vdo(lx1+loct) >*/
L60:
    vd = vdo[tabinf_1.lx1 + loct - 1];
/*<       go to 300 >*/
    goto L300;
/*<    70 xfact=delta/delold(2) >*/
L70:
    xfact = status_1.delta / status_1.delold[1];
/*<       vdo(lx0+loct)=vdo(lx1+loct) >*/
    vdo[tabinf_1.lx0 + loct - 1] = vdo[tabinf_1.lx1 + loct - 1];
/*<       vd=(1.0d0+xfact)*vdo(lx1+loct)-xfact*vdo(lx2+loct) >*/
    vd = (xfact + 1.) * vdo[tabinf_1.lx1 + loct - 1] - xfact * vdo[
	    tabinf_1.lx2 + loct - 1];
/*<       cdo(lx0+loct)=cdo(lx1+loct) >*/
    cdo[tabinf_1.lx0 + loct - 1] = cdo[tabinf_1.lx1 + loct - 1];
/*<       gdo(lx0+loct)=gdo(lx1+loct) >*/
    gdo[tabinf_1.lx0 + loct - 1] = gdo[tabinf_1.lx1 + loct - 1];
/*<       go to 110 >*/
    goto L110;

/*  compute new nonlinear branch voltage */

/*<   100 vd=value(lvnim1+node3)-value(lvnim1+node2) >*/
L100:
    vd = blank_1.value[tabinf_1.lvnim1 + node3 - 1] - blank_1.value[
	    tabinf_1.lvnim1 + node2 - 1];
/*<   110 delvd=vd-vdo(lx0+loct) >*/
L110:
    delvd = vd - vdo[tabinf_1.lx0 + loct - 1];
/*<       cdhat=cdo(lx0+loct)+gdo(lx0+loct)*delvd >*/
    cdhat = cdo[tabinf_1.lx0 + loct - 1] + gdo[tabinf_1.lx0 + loct - 1] * 
	    delvd;

/*  bypass if solution has not changed */

/*<       if (initf.eq.6) go to 200 >*/
    if (status_1.initf == 6) {
	goto L200;
    }
/*<       tol=reltol*dmax1(dabs(vd),dabs(vdo(lx0+loct)))+vntol >*/
/* Computing MAX */
    d_2 = abs(vd), d_3 = (d_1 = vdo[tabinf_1.lx0 + loct - 1], abs(d_1));
    tol = knstnt_1.reltol * max(d_3,d_2) + knstnt_1.vntol;
/*<       if (dabs(delvd).ge.tol) go to 200 >*/
    if (abs(delvd) >= 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;
    }
/*<       vd=vdo(lx0+loct) >*/
    vd = vdo[tabinf_1.lx0 + loct - 1];
/*<       cd=cdo(lx0+loct) >*/
    cd = cdo[tabinf_1.lx0 + loct - 1];
/*<       gd=gdo(lx0+loct) >*/
    gd = gdo[tabinf_1.lx0 + loct - 1];
/*<       go to 800 >*/
    goto L800;

/*  limit new junction voltage */

/*<   200 vlim=vte+vte >*/
L200:
    vlim = vte + vte;
/*<       if(bv.eq.0.0d0) go to 205 >*/
    if (bv == 0.) {
	goto L205;
    }
/*<       if (vd.lt.dmin1(0.0d0,-bv+10.0d0*vte)) go to 210 >*/
/* Computing MAX */
    d_1 = 0., d_2 = -bv + vte * 10.;
    if (vd < min(d_2,d_1)) {
	goto L210;
    }
/*<   205 call pnjlim(vd,vdo(lx0+loct),vte,vcrit,icheck) >*/
L205:
    pnjlim_(&vd, &vdo[tabinf_1.lx0 + loct - 1], &vte, &vcrit, &icheck);
/*<       go to 300 >*/
    goto L300;
/*<   210 vdtemp=-(vd+bv) >*/
L210:
    vdtemp = -(vd + bv);
/*<       call pnjlim(vdtemp,-(vdo(lx0+loct)+bv),vte,vcrit,icheck) >*/
    d_1 = -(vdo[tabinf_1.lx0 + loct - 1] + bv);
    pnjlim_(&vdtemp, &d_1, &vte, &vcrit, &icheck);
/*<       vd=-(vdtemp+bv) >*/
    vd = -(vdtemp + bv);

/*  compute dc current and derivitives */

/*<   300 if (vd.lt.-5.0d0*vte) go to 310 >*/
L300:
    if (vd < vte * -5.) {
	goto L310;
    }
/*<       evd=dexp(vd/vte) >*/
    evd = exp(vd / vte);
/*<       cd=csat*(evd-1.0d0)+gmin*vd >*/
    cd = csat * (evd - 1.) + knstnt_1.gmin * vd;
/*<       gd=csat*evd/vte+gmin >*/
    gd = csat * evd / vte + knstnt_1.gmin;
/*<       go to 330 >*/
    goto L330;
/*<   310 if(bv.eq.0.0d0) go to 315 >*/
L310:
    if (bv == 0.) {
	goto L315;
    }
/*<       if(vd.lt.-bv) go to 320 >*/
    if (vd < -bv) {
	goto L320;
    }
/*<   315 gd=-csat/vd+gmin >*/
L315:
    gd = -csat / vd + knstnt_1.gmin;
/*<       cd=gd*vd >*/
    cd = gd * vd;
/*<       go to 330 >*/
    goto L330;
/*<   320 evrev=dexp(-(bv+vd)/vt) >*/
L320:
    evrev = exp(-(bv + vd) / status_1.vt);
/*<       cd=-csat*(evrev-1.0d0+bv/vt) >*/
    cd = -csat * (evrev - 1. + bv / status_1.vt);
/*<       gd=csat*evrev/vt >*/
    gd = csat * evrev / status_1.vt;
/*<   330 if (mode.ne.1) go to 500 >*/
L330:
    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 tau=value(locm+4) >*/
L500:
    tau = blank_1.value[locm + 3];
/*<       czero=value(locm+5)*area >*/
    czero = blank_1.value[locm + 4] * area;
/*<       pb=value(locm+6) >*/
    pb = blank_1.value[locm + 5];
/*<       xm=value(locm+7) >*/
    xm = blank_1.value[locm + 6];
/*<       fcpb=value(locm+12) >*/
    fcpb = blank_1.value[locm + 11];
/*<       if (vd.ge.fcpb) go to 510 >*/
    if (vd >= fcpb) {
	goto L510;
    }
/*<       arg=1.0d0-vd/pb >*/
    arg = 1. - vd / pb;
/*<       sarg=dexp(-xm*dlog(arg)) >*/
    sarg = exp(-xm * log(arg));
/*<       qd(lx0+loct)=tau*cd+pb*czero*(1.0d0-arg*sarg)/(1.0d0-xm) >*/
    qd[tabinf_1.lx0 + loct - 1] = tau * cd + pb * czero * (1. - arg * sarg) / 
	    (1. - xm);
/*<       capd=tau*gd+czero*sarg >*/
    capd = tau * gd + czero * sarg;
/*<       go to 520 >*/
    goto L520;
/*<   510 f1=value(locm+15) >*/
L510:
    f1 = blank_1.value[locm + 14];
/*<       f2=value(locm+16) >*/
    f2 = blank_1.value[locm + 15];
/*<       f3=value(locm+17) >*/
    f3 = blank_1.value[locm + 16];
/*<       czof2=czero/f2 >*/
    czof2 = czero / f2;
/*<       qd(lx0+loct)=tau*cd+czero*f1+czof2*(f3*(vd-fcpb) >*/
/*<      1   +(xm/(pb+pb))*(vd*vd-fcpb*fcpb)) >*/
    qd[tabinf_1.lx0 + loct - 1] = tau * cd + czero * f1 + czof2 * (f3 * (vd - 
	    fcpb) + xm / (pb + pb) * (vd * vd - fcpb * fcpb));
/*<       capd=tau*gd+czof2*(f3+xm*vd/pb) >*/
    capd = tau * gd + czof2 * (f3 + xm * vd / pb);

/*  store small-signal parameters */

/*<   520 if ((mode.eq.1).and.(modedc.eq.2).and.(nosolv.ne.0)) go to 700 >*/
L520:
    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+4)=capd >*/
    blank_1.value[tabinf_1.lx0 + loct + 3] = capd;
/*<       go to 1000 >*/
    goto L1000;

/*  transient analysis */

/*<   600 if (initf.ne.5) go to 610 >*/
L600:
    if (status_1.initf != 5) {
	goto L610;
    }
/*<       qd(lx1+loct)=qd(lx0+loct) >*/
    qd[tabinf_1.lx1 + loct - 1] = qd[tabinf_1.lx0 + loct - 1];
/*<   610 call intgr8(geq,ceq,capd,loct+3) >*/
L610:
    i_1 = loct + 3;
    intgr8_(&geq, &ceq, &capd, &i_1);
/*<       gd=gd+geq >*/
    gd += geq;
/*<       cd=cd+cqd(lx0+loct) >*/
    cd += cqd[tabinf_1.lx0 + loct - 1];
/*<       if (initf.ne.5) go to 700 >*/
    if (status_1.initf != 5) {
	goto L700;
    }
/*<       cqd(lx1+loct)=cqd(lx0+loct) >*/
    cqd[tabinf_1.lx1 + loct - 1] = cqd[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(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 vdo(lx0+loct)=vd >*/
L750:
    vdo[tabinf_1.lx0 + loct - 1] = vd;
/*<       cdo(lx0+loct)=cd >*/
    cdo[tabinf_1.lx0 + loct - 1] = cd;
/*<       gdo(lx0+loct)=gd >*/
    gdo[tabinf_1.lx0 + loct - 1] = gd;

/*  load current vector */

/*<   800 cdeq=cd-gd*vd >*/
L800:
    cdeq = cd - gd * vd;
/*<       value(lvn+node2)=value(lvn+node2)+cdeq >*/
    blank_1.value[tabinf_1.lvn + node2 - 1] += cdeq;
/*<       value(lvn+node3)=value(lvn+node3)-cdeq >*/
    blank_1.value[tabinf_1.lvn + node3 - 1] -= cdeq;

/*  load matrix */

/*<       locy=lvn+nodplc(loc+13) >*/
    locy = tabinf_1.lvn + nodplc[loc + 12];
/*<       value(locy)=value(locy)+gspr >*/
    blank_1.value[locy - 1] += gspr;
/*<       locy=lvn+nodplc(loc+14) >*/
    locy = tabinf_1.lvn + nodplc[loc + 13];
/*<       value(locy)=value(locy)+gd >*/
    blank_1.value[locy - 1] += gd;
/*<       locy=lvn+nodplc(loc+15) >*/
    locy = tabinf_1.lvn + nodplc[loc + 14];
/*<       value(locy)=value(locy)+gd+gspr >*/
    blank_1.value[locy - 1] = blank_1.value[locy - 1] + gd + gspr;
/*<       locy=lvn+nodplc(loc+7) >*/
    locy = tabinf_1.lvn + nodplc[loc + 6];
/*<       value(locy)=value(locy)-gspr >*/
    blank_1.value[locy - 1] -= gspr;
/*<       locy=lvn+nodplc(loc+8) >*/
    locy = tabinf_1.lvn + nodplc[loc + 7];
/*<       value(locy)=value(locy)-gd >*/
    blank_1.value[locy - 1] -= gd;
/*<       locy=lvn+nodplc(loc+9) >*/
    locy = tabinf_1.lvn + nodplc[loc + 8];
/*<       value(locy)=value(locy)-gspr >*/
    blank_1.value[locy - 1] -= gspr;
/*<       locy=lvn+nodplc(loc+10) >*/
    locy = tabinf_1.lvn + nodplc[loc + 9];
/*<       value(locy)=value(locy)-gd >*/
    blank_1.value[locy - 1] -= gd;
/*<  1000 loc=nodplc(loc) >*/
L1000:
    loc = nodplc[loc - 1];
/*<       go to 10 >*/
    goto L10;
/*<       end >*/
} /* diode_ */

#undef vdo
#undef gdo
#undef cqd
#undef cdo
#undef cvalue
#undef nodplc
#undef qd


