/* moseq3.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 {
    doublereal vto, beta, gamma, phi, phib, cox, xnsub, xnfs, xd, xj, xld, 
	    xlamda, uo, uexp, vbp, utra, vmax, xneff, xl, xw, vbi, von, vdsat,
	     qspof, beta0, beta1, cdrain, xqco, xqc, fnarrw, fshort;
    integer lev;
} mosarg_;

#define mosarg_1 mosarg_

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_

/*<       subroutine moseq3(vds,vbs,vgs,gm,gds,gmbs, >*/
/*<      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb) >*/
/* Subroutine */ int moseq3_(vds, vbs, vgs, gm, gds, gmbs, qg, qc, qb, cggb, 
	cgdb, cgsb, cbgb, cbdb, cbsb)
doublereal *vds, *vbs, *vgs, *gm, *gds, *gmbs, *qg, *qc, *qb, *cggb, *cgdb, *
	cgsb, *cbgb, *cbdb, *cbsb;
{
    /* Initialized data */

    static doublereal coeff0 = .0631353;
    static doublereal coeff1 = .8013292;
    static doublereal coeff2 = -.01110777;

    /* System generated locals */
    doublereal d_1, d_2;

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

    /* Local variables */
    static doublereal arga, argc, argb, onfg, emax, vdsc, vbin, vbix, vpof, 
	    onxj, onxl, onxn, vdsx, vgsx;
    extern /* Subroutine */ int mosq3_();
    static doublereal dadvb, dbdvb, diddl;
#define alpha ((doublereal *)&mosarg_1 + 11)
    static doublereal fgate, dgdvb, dgdvg, dgdvd, dldvd, cdsat, dldem, gdsat, 
	    fbody, wfact;
#define theta ((doublereal *)&mosarg_1 + 14)
    static doublereal phibs, delxl, ondvt, vdsat1, sqphs3, dfbdvb, dfddvb;
    static integer icharg;
    static doublereal dfgdvd, dfgdvb, dfgdvg, dfsdvb, gammas, dqbdvb, cdonco, 
	    xn, dsqdvb;
#define xkappa ((doublereal *)&mosarg_1 + 15)
    static doublereal dwpdvb, fbodys, onfbdy, djonxj, qbonco, sqphbs, dvtdvd, 
	    sqphis, wconxj, dvtdvb, csonco, dxndvb, dvodvd, dvodvb, us, 
	    xjonxl, dvsdvg, dvsdvd, wponxj, dvsdvb, onvdsc, dvsdga, dcodvg, 
	    dcodvd, dcodvb, cdnorm, cd1, fdrain, fd2, dfddvg, dfddvd, gdoncd, 
	    gdonfd, gdonfg, emoncd, emongd, demdvg, demdvd, demdvb, ddldvg, 
	    ddldvd, ddldvb, dlonxl, xlfact;
    extern /* Subroutine */ int mqspof_();
    static doublereal gamasd, cdo;
#define eta ((doublereal *)&mosarg_1 + 13)
    static doublereal vgb, vfb, wcs, gms, vth, gmw, wps, gds0, wcs2;

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

/*     this routine evaluates the drain current, its derivatives and */
/*     the charges associated with the gate, channel and bulk */
/*     for mosfets based on semi-empirical equations */

/* spice version 2g.6  sccsid=mosarg 3/15/83 */
/*<       common /mosarg/ vto,beta,gamma,phi,phib,cox,xnsub,xnfs,xd,xj,xld, >*/
/*<      1   xlamda,uo,uexp,vbp,utra,vmax,xneff,xl,xw,vbi,von,vdsat,qspof, >*/
/*<      2   beta0,beta1,cdrain,xqco,xqc,fnarrw,fshort,lev >*/
/* 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 >*/

/*<       equivalence (xlamda,alpha),(vbp,theta),(uexp,eta),(utra,xkappa) >*/
/*<       data coeff0/0.0631353d0/,coeff1/0.8013292d0/,coeff2/-0.01110777d0/ >*/

/*     icharg=1 causes charges to be computed */
/*     icharg=0 bypasses the computation of charges */

/*     icharg=1 */
/*     if (mode.ne.1) go to 10 */
/*     icharg=0 */
/*     if (modedc.eq.2.and.nosolv.ne.0) icharg=1 */
/*     if (initf.eq.4) icharg=1 */

/*     reference cdrain equations to source and */
/*     charge equations to bulk */

/*< 10    continue >*/
/* L10: */
/*<       icharg=0 >*/
    icharg = 0;
/*<       vgb=vgs-vbs >*/
    vgb = *vgs - *vbs;
/*<       vfb=vbi-phi >*/
    vfb = mosarg_1.vbi - mosarg_1.phi;
/*<       vdsat=0.0d0 >*/
    mosarg_1.vdsat = 0.;
/*<       qg=0.0d0 >*/
    *qg = 0.;
/*<       qb=0.0d0 >*/
    *qb = 0.;
/*<       qc=0.0d0 >*/
    *qc = 0.;
/*<       cgdb=0.0d0 >*/
    *cgdb = 0.;
/*<       cbdb=0.0d0 >*/
    *cbdb = 0.;
/*<       onxl=1.0d0/xl >*/
    onxl = 1. / mosarg_1.xl;
/*<       eta=eta/(xl*xl*xl) >*/
    *eta /= mosarg_1.xl * mosarg_1.xl * mosarg_1.xl;

/* .....square root term */

/*<       if ( vbs.gt.0.0d0 ) go to 120 >*/
    if (*vbs > 0.) {
	goto L120;
    }
/*<       phibs=phi-vbs >*/
    phibs = mosarg_1.phi - *vbs;
/*<       sqphbs=dsqrt(phibs) >*/
    sqphbs = sqrt(phibs);
/*<       dsqdvb=-0.5d0/sqphbs >*/
    dsqdvb = -.5 / sqphbs;
/*<       go to 200 >*/
    goto L200;
/*< 120   continue >*/
L120:
/*<       sqphis=dsqrt(phi) >*/
    sqphis = sqrt(mosarg_1.phi);
/*<       sqphs3=phi*sqphis >*/
    sqphs3 = mosarg_1.phi * sqphis;
/*<       sqphbs=sqphis/(1.0d0+vbs/(phi+phi)) >*/
    sqphbs = sqphis / (*vbs / (mosarg_1.phi + mosarg_1.phi) + 1.);
/*<       phibs=sqphbs*sqphbs >*/
    phibs = sqphbs * sqphbs;
/*<       dsqdvb=-phibs/(sqphs3+sqphs3) >*/
    dsqdvb = -phibs / (sqphs3 + sqphs3);

/* .....short channel effect factor */

/*< 200   continue >*/
L200:
/*<       if ( (xj.eq.0.0d0).or.(xd.eq.0.0d0) ) go to 210 >*/
    if (mosarg_1.xj == 0. || mosarg_1.xd == 0.) {
	goto L210;
    }
/*<       wps=xd*sqphbs >*/
    wps = mosarg_1.xd * sqphbs;
/*<       onxj=1.0d0/xj >*/
    onxj = 1. / mosarg_1.xj;
/*<       xjonxl=xj*onxl >*/
    xjonxl = mosarg_1.xj * onxl;
/*<       djonxj=xld*onxj >*/
    djonxj = mosarg_1.xld * onxj;
/*<       wponxj=wps*onxj >*/
    wponxj = wps * onxj;
/*<       wconxj=coeff0+coeff1*wponxj+coeff2*wponxj*wponxj >*/
    wconxj = coeff0 + coeff1 * wponxj + coeff2 * wponxj * wponxj;
/*<       wcs=wconxj*xj >*/
    wcs = wconxj * mosarg_1.xj;
/*<       arga=wconxj+djonxj >*/
    arga = wconxj + djonxj;
/*<       argc=wponxj/(1.0d0+wponxj) >*/
    argc = wponxj / (wponxj + 1.);
/*<       argb=dsqrt(1.0d0-argc*argc) >*/
    argb = sqrt(1. - argc * argc);
/*<       fshort=1.0d0-xjonxl*(arga*argb-djonxj) >*/
    mosarg_1.fshort = 1. - xjonxl * (arga * argb - djonxj);
/*<       dwpdvb=xd*dsqdvb >*/
    dwpdvb = mosarg_1.xd * dsqdvb;
/*<       dadvb=(coeff1+coeff2*(wponxj+wponxj))*dwpdvb*onxj >*/
    dadvb = (coeff1 + coeff2 * (wponxj + wponxj)) * dwpdvb * onxj;
/*<       dbdvb=-argc*argc*(1.0d0-argc)*dwpdvb/(argb*wps) >*/
    dbdvb = -argc * argc * (1. - argc) * dwpdvb / (argb * wps);
/*<       dfsdvb=-xjonxl*(dadvb*argb+arga*dbdvb) >*/
    dfsdvb = -xjonxl * (dadvb * argb + arga * dbdvb);
/*<       go to 220 >*/
    goto L220;
/*< 210   continue >*/
L210:
/*<       fshort=1.0d0 >*/
    mosarg_1.fshort = 1.;
/*<       dfsdvb=0.0d0 >*/
    dfsdvb = 0.;
/*<       wcs=0.05d-6 >*/
    wcs = 5e-8;

/* .....body effect */

/*< 220   continue >*/
L220:
/*<       gammas=gamma*fshort >*/
    gammas = mosarg_1.gamma * mosarg_1.fshort;
/*<       fbodys=0.5d0*gammas/(sqphbs+sqphbs) >*/
    fbodys = gammas * .5 / (sqphbs + sqphbs);
/*<       fbody=fbodys+fnarrw >*/
    fbody = fbodys + mosarg_1.fnarrw;
/*<       onfbdy=1.0d0/(1.0d0+fbody) >*/
    onfbdy = 1. / (fbody + 1.);
/*<       dfbdvb=-fbodys*dsqdvb/sqphbs+fbodys*dfsdvb/fshort >*/
    dfbdvb = -fbodys * dsqdvb / sqphbs + fbodys * dfsdvb / mosarg_1.fshort;
/*<       qbonco=gammas*sqphbs+fnarrw*phibs >*/
    qbonco = gammas * sqphbs + mosarg_1.fnarrw * phibs;
/*<       dqbdvb=gammas*dsqdvb+gamma*dfsdvb*sqphbs-fnarrw >*/
    dqbdvb = gammas * dsqdvb + mosarg_1.gamma * dfsdvb * sqphbs - 
	    mosarg_1.fnarrw;

/* .....static feedback effect */

/*<       vbix=vbi-eta*vds >*/
    vbix = mosarg_1.vbi - *eta * *vds;

/* .....threshold voltage */

/*<       vth=vbix+qbonco >*/
    vth = vbix + qbonco;
/*<       dvtdvd=-eta >*/
    dvtdvd = -(*eta);
/*<       dvtdvb=dqbdvb >*/
    dvtdvb = dqbdvb;

/* .....joint weak inversion and strong inversion */

/*<       von=vth >*/
    mosarg_1.von = vth;
/*<       if ( xnfs.eq.0.0d0 ) go to 250 >*/
    if (mosarg_1.xnfs == 0.) {
	goto L250;
    }
/*<            csonco=charge*xnfs*xl*xw/cox >*/
    csonco = knstnt_1.charge * mosarg_1.xnfs * mosarg_1.xl * mosarg_1.xw / 
	    mosarg_1.cox;
/*<            cdonco=qbonco/(phibs+phibs) >*/
    cdonco = qbonco / (phibs + phibs);
/*<            xn=1.0d0+csonco+cdonco >*/
    xn = csonco + 1. + cdonco;
/*<            von=vth+vt*xn >*/
    mosarg_1.von = vth + status_1.vt * xn;
/*<            dxndvb=dqbdvb/(phibs+phibs)-qbonco*dsqdvb/(phibs*sqphbs) >*/
    dxndvb = dqbdvb / (phibs + phibs) - qbonco * dsqdvb / (phibs * sqphbs);
/*<            dvodvd=dvtdvd >*/
    dvodvd = dvtdvd;
/*<            dvodvb=dvtdvb+vt*dxndvb >*/
    dvodvb = dvtdvb + status_1.vt * dxndvb;
/*<            go to 300 >*/
    goto L300;

/* .....cutoff region */

/*< 250   continue >*/
L250:
/*<       if ( vgs.gt.von ) go to 300 >*/
    if (*vgs > mosarg_1.von) {
	goto L300;
    }
/*<       cdrain=0.0d0 >*/
    mosarg_1.cdrain = 0.;
/*<       gm=0.0d0 >*/
    *gm = 0.;
/*<       gds=0.0d0 >*/
    *gds = 0.;
/*<       gmbs=0.0d0 >*/
    *gmbs = 0.;
/*<       if ( icharg.ne.0 ) go to 800 >*/
    if (icharg != 0) {
	goto L800;
    }
/*<       go to 1000 >*/
    goto L1000;

/* .....device is on */

/*< 300   continue >*/
L300:
/*<       vgsx=dmax1(vgs,von) >*/
    vgsx = max(*vgs,mosarg_1.von);

/* .....mobility modulation by gate voltage */

/*<       onfg=1.0d0+theta*(vgsx-vth) >*/
    onfg = *theta * (vgsx - vth) + 1.;
/*<       fgate=1.0d0/onfg >*/
    fgate = 1. / onfg;
/*<       us=uo*fgate >*/
    us = mosarg_1.uo * fgate;
/*<       dfgdvg=-theta*fgate*fgate >*/
    dfgdvg = -(*theta) * fgate * fgate;
/*<       dfgdvd=-dfgdvg*dvtdvd >*/
    dfgdvd = -dfgdvg * dvtdvd;
/*<       dfgdvb=-dfgdvg*dvtdvb >*/
    dfgdvb = -dfgdvg * dvtdvb;

/* .....saturation voltage */

/*<       vdsat=(vgsx-vth)*onfbdy >*/
    mosarg_1.vdsat = (vgsx - vth) * onfbdy;
/*<       vpof=vdsat >*/
    vpof = mosarg_1.vdsat;
/*<       if ( vmax.gt.0.0d0 ) go to 310 >*/
    if (mosarg_1.vmax > 0.) {
	goto L310;
    }
/*<       dvsdvg=onfbdy >*/
    dvsdvg = onfbdy;
/*<       dvsdvd=-dvsdvg*dvtdvd >*/
    dvsdvd = -dvsdvg * dvtdvd;
/*<       dvsdvb=-dvsdvg*dvtdvb-vdsat*dfbdvb*onfbdy >*/
    dvsdvb = -dvsdvg * dvtdvb - mosarg_1.vdsat * dfbdvb * onfbdy;
/*<       go to 400 >*/
    goto L400;
/*<   310 vdsc=xl*vmax/us >*/
L310:
    vdsc = mosarg_1.xl * mosarg_1.vmax / us;
/*<       onvdsc=1.0d0/vdsc >*/
    onvdsc = 1. / vdsc;
/*<       arga=(vgsx-vth)*onfbdy >*/
    arga = (vgsx - vth) * onfbdy;
/*<       argb=dsqrt(arga*arga+vdsc*vdsc) >*/
    argb = sqrt(arga * arga + vdsc * vdsc);
/*<       vdsat=arga+vdsc-argb >*/
    mosarg_1.vdsat = arga + vdsc - argb;
/*<       dvsdga=(1.0d0-arga/argb)*onfbdy >*/
    dvsdga = (1. - arga / argb) * onfbdy;
/*<       dvsdvg=dvsdga-(1.0d0-vdsc/argb)*vdsc*dfgdvg*onfg >*/
    dvsdvg = dvsdga - (1. - vdsc / argb) * vdsc * dfgdvg * onfg;
/*<       dvsdvd=-dvsdvg*dvtdvd >*/
    dvsdvd = -dvsdvg * dvtdvd;
/*<       dvsdvb=-dvsdvg*dvtdvb-arga*dvsdga*dfbdvb >*/
    dvsdvb = -dvsdvg * dvtdvb - arga * dvsdga * dfbdvb;

/* .....current factors in linear region */

/*< 400   continue >*/
L400:
/*<       vdsx=dmin1(vds,vdsat) >*/
    vdsx = min(*vds,mosarg_1.vdsat);
/*<       if ( vdsx.eq.0.0d0 ) go to 900 >*/
    if (vdsx == 0.) {
	goto L900;
    }
/*<       cdo=vgsx-vth-0.5d0*(1.0d0+fbody)*vdsx >*/
    cdo = vgsx - vth - (fbody + 1.) * .5 * vdsx;
/*<       dcodvg=1.0d0 >*/
    dcodvg = 1.;
/*<       if (vds.lt.vdsat) dcodvd=-dvtdvd-0.5d0*(1.0d0+fbody) >*/
    if (*vds < mosarg_1.vdsat) {
	dcodvd = -dvtdvd - (fbody + 1.) * .5;
    }
/*<       dcodvb=-dvtdvb-0.5d0*dfbdvb*vdsx >*/
    dcodvb = -dvtdvb - dfbdvb * .5 * vdsx;

/* .....normalized drain current */

/*< 410   continue >*/
/* L410: */
/*<       cdnorm=cdo*vdsx >*/
    cdnorm = cdo * vdsx;
/*<       gm=vdsx >*/
    *gm = vdsx;
/*<       gds=vgsx-vth-(1.0d0+fbody+dvtdvd)*vdsx >*/
    *gds = vgsx - vth - (fbody + 1. + dvtdvd) * vdsx;
/*<       gmbs=dcodvb*vdsx >*/
    *gmbs = dcodvb * vdsx;

/* .....drain current without velocity saturation effect */

/*<       cd1=beta*cdnorm >*/
    cd1 = mosarg_1.beta * cdnorm;
/*<       beta=beta*fgate >*/
    mosarg_1.beta *= fgate;
/*<       cdrain=beta*cdnorm >*/
    mosarg_1.cdrain = mosarg_1.beta * cdnorm;
/*<       gm=beta*gm+dfgdvg*cd1 >*/
    *gm = mosarg_1.beta * *gm + dfgdvg * cd1;
/*<       gds=beta*gds+dfgdvd*cd1 >*/
    *gds = mosarg_1.beta * *gds + dfgdvd * cd1;
/*<       gmbs=beta*gmbs >*/
    *gmbs = mosarg_1.beta * *gmbs;

/* .....velocity saturation factor */

/*<       if ( vmax.eq.0.0d0 ) go to 500 >*/
    if (mosarg_1.vmax == 0.) {
	goto L500;
    }
/*<       fdrain=1.0d0/(1.0d0+vdsx*onvdsc) >*/
    fdrain = 1. / (vdsx * onvdsc + 1.);
/*<       fd2=fdrain*fdrain >*/
    fd2 = fdrain * fdrain;
/*<       arga=fd2*vdsx*onvdsc*onfg >*/
    arga = fd2 * vdsx * onvdsc * onfg;
/*<       dfddvg=-dfgdvg*arga >*/
    dfddvg = -dfgdvg * arga;
/*<       dfddvd=-dfgdvd*arga-fd2*onvdsc >*/
    dfddvd = -dfgdvd * arga - fd2 * onvdsc;
/*<       dfddvb=-dfgdvb*arga >*/
    dfddvb = -dfgdvb * arga;

/* .....drain current */

/*<       gm=fdrain*gm+dfddvg*cdrain >*/
    *gm = fdrain * *gm + dfddvg * mosarg_1.cdrain;
/*<       gds=fdrain*gds+dfddvd*cdrain >*/
    *gds = fdrain * *gds + dfddvd * mosarg_1.cdrain;
/*<       gmbs=fdrain*gmbs+dfddvb*cdrain >*/
    *gmbs = fdrain * *gmbs + dfddvb * mosarg_1.cdrain;
/*<       cdrain=fdrain*cdrain >*/
    mosarg_1.cdrain = fdrain * mosarg_1.cdrain;
/*<       beta=beta*fdrain >*/
    mosarg_1.beta *= fdrain;

/* .....channel length modulation */

/*< 500   continue >*/
L500:
/*<       if ( vds.le.vdsat ) go to 700 >*/
    if (*vds <= mosarg_1.vdsat) {
	goto L700;
    }
/*<       if ( vmax.eq.0.0d0 ) go to 510 >*/
    if (mosarg_1.vmax == 0.) {
	goto L510;
    }
/*<       if (alpha.eq.0.0d0) go to 700 >*/
    if (*alpha == 0.) {
	goto L700;
    }
/*<       cdsat=cdrain >*/
    cdsat = mosarg_1.cdrain;
/*<       gdsat=cdsat*(1.0d0-fdrain)*onvdsc >*/
    gdsat = cdsat * (1. - fdrain) * onvdsc;
/*<       gdsat=dmax1(1.0d-12,gdsat) >*/
    gdsat = max(1e-12,gdsat);
/*<       gdoncd=gdsat/cdsat >*/
    gdoncd = gdsat / cdsat;
/*<       gdonfd=gdsat/(1.0d0-fdrain) >*/
    gdonfd = gdsat / (1. - fdrain);
/*<       gdonfg=gdsat*onfg >*/
    gdonfg = gdsat * onfg;
/*<       dgdvg=gdoncd*gm-gdonfd*dfddvg+gdonfg*dfgdvg >*/
    dgdvg = gdoncd * *gm - gdonfd * dfddvg + gdonfg * dfgdvg;
/*<       dgdvd=gdoncd*gds-gdonfd*dfddvd+gdonfg*dfgdvd >*/
    dgdvd = gdoncd * *gds - gdonfd * dfddvd + gdonfg * dfgdvd;
/*<       dgdvb=gdoncd*gmbs-gdonfd*dfddvb+gdonfg*dfgdvb >*/
    dgdvb = gdoncd * *gmbs - gdonfd * dfddvb + gdonfg * dfgdvb;

/*<       emax=cdsat*onxl/gdsat >*/
    emax = cdsat * onxl / gdsat;
/*<       emoncd=emax/cdsat >*/
    emoncd = emax / cdsat;
/*<       emongd=emax/gdsat >*/
    emongd = emax / gdsat;
/*<       demdvg=emoncd*gm-emongd*dgdvg >*/
    demdvg = emoncd * *gm - emongd * dgdvg;
/*<       demdvd=emoncd*gds-emongd*dgdvd >*/
    demdvd = emoncd * *gds - emongd * dgdvd;
/*<       demdvb=emoncd*gmbs-emongd*dgdvb >*/
    demdvb = emoncd * *gmbs - emongd * dgdvb;

/*<       arga=0.5d0*emax*alpha >*/
    arga = emax * .5 * *alpha;
/*<       argc=xkappa*alpha >*/
    argc = *xkappa * *alpha;
/*<       argb=dsqrt(arga*arga+argc*(vds-vdsat)) >*/
    argb = sqrt(arga * arga + argc * (*vds - mosarg_1.vdsat));
/*<       delxl=argb-arga >*/
    delxl = argb - arga;
/*<       dldvd=argc/(argb+argb) >*/
    dldvd = argc / (argb + argb);
/*<       dldem=0.5d0*(arga/argb-1.0d0)*alpha >*/
    dldem = (arga / argb - 1.) * .5 * *alpha;
/*<       ddldvg=dldem*demdvg >*/
    ddldvg = dldem * demdvg;
/*<       ddldvd=dldem*demdvd-dldvd >*/
    ddldvd = dldem * demdvd - dldvd;
/*<       ddldvb=dldem*demdvb >*/
    ddldvb = dldem * demdvb;
/*<       go to 520 >*/
    goto L520;
/*< 510   continue >*/
L510:
/*<       delxl=dsqrt(xkappa*(vds-vdsat)*alpha) >*/
    delxl = sqrt(*xkappa * (*vds - mosarg_1.vdsat) * *alpha);
/*<       dldvd=0.5d0*delxl/(vds-vdsat) >*/
    dldvd = delxl * .5 / (*vds - mosarg_1.vdsat);
/*<       ddldvg=0.0d0 >*/
    ddldvg = 0.;
/*<       ddldvd=-dldvd >*/
    ddldvd = -dldvd;
/*<       ddldvb=0.0d0 >*/
    ddldvb = 0.;

/* .....punch through approximation */

/*< 520   continue >*/
L520:
/*<       if ( delxl.le.(0.5d0*xl) ) go to 600 >*/
    if (delxl <= mosarg_1.xl * .5) {
	goto L600;
    }
/*<       wcs2=wcs*wcs >*/
    wcs2 = wcs * wcs;
/*<       delxl=xl-(xl**2/(4.0d0*delxl)) >*/
/* Computing 2nd power */
    d_1 = mosarg_1.xl;
    delxl = mosarg_1.xl - d_1 * d_1 / (delxl * 4.);
/*<       arga=4.0d0*(xl-delxl)**2/xl**2 >*/
/* Computing 2nd power */
    d_1 = mosarg_1.xl - delxl;
/* Computing 2nd power */
    d_2 = mosarg_1.xl;
    arga = d_1 * d_1 * 4. / (d_2 * d_2);
/*<       ddldvg=ddldvg*arga >*/
    ddldvg *= arga;
/*<       ddldvd=ddldvd*arga >*/
    ddldvd *= arga;
/*<       ddldvb=ddldvb*arga >*/
    ddldvb *= arga;
/*<        dldvd= dldvd*arga >*/
    dldvd *= arga;

/* .....saturation region */

/*< 600   continue >*/
L600:
/*<       dlonxl=delxl*onxl >*/
    dlonxl = delxl * onxl;
/*<       xlfact=1.0d0/(1.0d0-dlonxl) >*/
    xlfact = 1. / (1. - dlonxl);
/*<       cdrain=cdrain*xlfact >*/
    mosarg_1.cdrain *= xlfact;
/*<       diddl=cdrain/(xl-delxl) >*/
    diddl = mosarg_1.cdrain / (mosarg_1.xl - delxl);
/*<       gm=gm*xlfact+diddl*ddldvg >*/
    *gm = *gm * xlfact + diddl * ddldvg;
/*<       gds0=gds*xlfact+diddl*ddldvd >*/
    gds0 = *gds * xlfact + diddl * ddldvd;
/*<       gmbs=gmbs*xlfact+diddl*ddldvb >*/
    *gmbs = *gmbs * xlfact + diddl * ddldvb;
/*<       gm=gm+gds0*dvsdvg >*/
    *gm += gds0 * dvsdvg;
/*<       gmbs=gmbs+gds0*dvsdvb >*/
    *gmbs += gds0 * dvsdvb;
/*<       gds=gds0*dvsdvd+diddl*dldvd >*/
    *gds = gds0 * dvsdvd + diddl * dldvd;

/* .....finish strong inversion case */

/*< 700   continue >*/
L700:
/*<       if ( vgs.ge.von ) go to 750 >*/
    if (*vgs >= mosarg_1.von) {
	goto L750;
    }

/* .....weak inversion */

/*<                 onxn=1.0d0/xn >*/
    onxn = 1. / xn;
/*<                 ondvt=onxn/vt >*/
    ondvt = onxn / status_1.vt;
/*<                 wfact=dexp( (vgs-von)*ondvt ) >*/
    wfact = exp((*vgs - mosarg_1.von) * ondvt);
/*<                 cdrain=cdrain*wfact >*/
    mosarg_1.cdrain *= wfact;
/*<                 gms=gm*wfact >*/
    gms = *gm * wfact;
/*<                 gmw=cdrain*ondvt >*/
    gmw = mosarg_1.cdrain * ondvt;
/*<                 gm=gmw >*/
    *gm = gmw;
/*<                 if (vds.gt.vdsat) gm=gm+gds0*dvsdvg*wfact >*/
    if (*vds > mosarg_1.vdsat) {
	*gm += gds0 * dvsdvg * wfact;
    }
/*<                 gds=gds*wfact+(gms-gmw)*dvodvd >*/
    *gds = *gds * wfact + (gms - gmw) * dvodvd;
/*<                 gmbs=gmbs*wfact+(gms-gmw)*dvodvb >*/
/*<      1                         -gmw*(vgs-von)*onxn*dxndvb >*/
    *gmbs = *gmbs * wfact + (gms - gmw) * dvodvb - gmw * (*vgs - mosarg_1.von)
	     * onxn * dxndvb;

/* .....charge computation */

/*<   750 continue >*/
L750:
/*<       if (icharg.eq.0) go to 1000 >*/
    if (icharg == 0) {
	goto L1000;
    }
/*<       if (vgs.le.vth) go to 800 >*/
    if (*vgs <= vth) {
	goto L800;
    }
/*<       call mqspof(vds,vbs,vgs,vpof,vdsat1,vth,vbin,gamasd, >*/
/*<      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb) >*/
    mqspof_(vds, vbs, vgs, &vpof, &vdsat1, &vth, &vbin, &gamasd, qg, qc, qb, 
	    cggb, cgdb, cgsb, cbgb, cbdb, cbsb);
/*<       go to 2000 >*/
    goto L2000;

/* .....charge computation for vgs<vth */

/*< 800   continue >*/
L800:
/*<       xqc=xqco >*/
    mosarg_1.xqc = mosarg_1.xqco;
/*<       call mosq3(vds,vbs,vpof,vdsat1,vth,vbin,gamasd,cox,phi, >*/
/*<      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb) >*/
    mosq3_(vds, vbs, &vpof, &vdsat1, &vth, &vbin, &gamasd, &mosarg_1.cox, &
	    mosarg_1.phi, qg, qc, qb, cggb, cgdb, cgsb, cbgb, cbdb, cbsb);
/*<       qspof=0.0d0 >*/
    mosarg_1.qspof = 0.;
/*<       go to 2000 >*/
    goto L2000;

/* .....special case of vds=0.0d0 */

/*< 900   continue >*/
L900:
/*<       beta=beta*fgate >*/
    mosarg_1.beta *= fgate;
/*<       cdrain=0.0d0 >*/
    mosarg_1.cdrain = 0.;
/*<       gm=0.0d0 >*/
    *gm = 0.;
/*<       gds=beta*(vgsx-vth) >*/
    *gds = mosarg_1.beta * (vgsx - vth);
/*<       gmbs=0.0d0 >*/
    *gmbs = 0.;
/*<            if ( (xnfs.ne.0.0d0).and.(vgs.lt.von) ) >*/
/*<      1          gds=gds*dexp((vgs-von)/(vt*xn)) >*/
    if (mosarg_1.xnfs != 0. && *vgs < mosarg_1.von) {
	*gds *= exp((*vgs - mosarg_1.von) / (status_1.vt * xn));
    }
/*<       if (icharg.eq.0) go to 1000 >*/
    if (icharg == 0) {
	goto L1000;
    }
/*<       call mosq3(vds,vbs,vpof,vdsat1,vth,vbin,gamasd,cox,phi, >*/
/*<      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb) >*/
    mosq3_(vds, vbs, &vpof, &vdsat1, &vth, &vbin, &gamasd, &mosarg_1.cox, &
	    mosarg_1.phi, qg, qc, qb, cggb, cgdb, cgsb, cbgb, cbdb, cbsb);
/*< 1000  qspof=0.0d0 >*/
L1000:
    mosarg_1.qspof = 0.;

/* .....done */

/*<  2000 return >*/
L2000:
    return 0;
/*<       end >*/
} /* moseq3_ */

#undef eta
#undef xkappa
#undef theta
#undef alpha


