/* moseq2.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 moseq2(vds,vbs,vgs,gm,gds,gmbs, >*/
/*<      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb) >*/
/* Subroutine */ int moseq2_(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 sig1[4] = { 1.,-1.,1.,-1. };
    static doublereal sig2[4] = { 1.,1.,-1.,-1. };

    /* System generated locals */
    integer i_1;
    doublereal d_1;

    /* Builtin functions */
    double sqrt(), log(), exp(), atan(), cos();

    /* Local variables */
    static doublereal barg, argd, argg, ueff, sarg, vbin, args, body, sphi, 
	    argv;
    static integer iknt, jknt;
    static doublereal expg, vpof, vgst, vgsx, sarg3, argv1, sphi3, d2bdb2, a, 
	    b, c;
    extern /* Subroutine */ int mosq2_();
    static doublereal poly4[8];
    static integer i, j;
    static doublereal d2sdb2, p, r, s, gdbdv, dfact, bsarg, argsd, ufact, 
	    xleff, argxd, cdson, dbxwd, ddxwd, gdson, argsq, gbson, argss, 
	    tsarg, daddb2, a1, argxs, a4[4], b4[4], dgddb2, dbxws, b1, c1, 
	    dasdb2, d1, r3, s2, gammd2, v1, v2, p0, p2, x4[8], y3, p3, p4, a3,
	     b3, delta4, bodys, sargv, vdson, vdsat1, fi, dbrgdb, dbargd, 
	    dgddvb;
    static integer icharg;
    static doublereal gamasd, gammad, ro, dsrgdb, dbargs, dgdvds, gamass, 
	    factor, dgsdvb, cdonco, xn, sbiarg, dgdvbs, dodvbs, dxndvb, 
	    dxndvd, dodvds, udenom, dudvgs, dudvds, dudvbs, dsdvgs, dsdvbs, 
	    xv, xvalid, dbsrdb, gdbdvs, xlfact, dldsat, vqchan, dqdsat, vl, 
	    dfunds, dfundg, dfundb, dldvgs, dldvds, dldvbs, clfact, deltal, 
	    didvds;
    extern /* Subroutine */ int mqspof_();
    static doublereal eta, arg, cfs, gmw, xwb, vth, xwd, xdv, xls, xlv, xws, 
	    arg1, arg2;

/*<       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 */

/* 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 >*/

/*<       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/ >*/

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

/*<       icharg=1 >*/
    icharg = 1;
/*<       if (mode.ne.1.and.xqco.le.0.5d0) go to 100 >*/
    if (status_1.mode != 1 && mosarg_1.xqco <= .5) {
	goto L100;
    }
/*<       icharg=0 >*/
    icharg = 0;
/*<       if (xqco.gt.0.5d0) go to 100 >*/
    if (mosarg_1.xqco > .5) {
	goto L100;
    }
/*<       if (modedc.eq.2.and.nosolv.ne.0) icharg=1 >*/
    if (status_1.modedc == 2 && status_1.nosolv != 0) {
	icharg = 1;
    }
/*<       if (initf.eq.4) icharg=1 >*/
    if (status_1.initf == 4) {
	icharg = 1;
    }

/*  compute some useful quantities */

/*<   100 if (vbs.gt.0.0d0) go to 110 >*/
L100:
    if (*vbs > 0.) {
	goto L110;
    }
/*<       sarg=dsqrt(phi-vbs) >*/
    sarg = sqrt(mosarg_1.phi - *vbs);
/*<       tsarg=sarg+sarg >*/
    tsarg = sarg + sarg;
/*<       dsrgdb=-0.5d0/sarg >*/
    dsrgdb = -.5 / sarg;
/*<       d2sdb2=+0.5d0*dsrgdb/(phi-vbs) >*/
    d2sdb2 = dsrgdb * .5 / (mosarg_1.phi - *vbs);
/*<       go to 120 >*/
    goto L120;
/*<   110 sphi=dsqrt(phi) >*/
L110:
    sphi = sqrt(mosarg_1.phi);
/*<       sphi3=phi*sphi >*/
    sphi3 = mosarg_1.phi * sphi;
/*<       sarg=sphi/(1.0d0+0.5d0*vbs/phi) >*/
    sarg = sphi / (*vbs * .5 / mosarg_1.phi + 1.);
/*<       tsarg=sarg+sarg >*/
    tsarg = sarg + sarg;
/*<       dsrgdb=-0.5d0*sarg*sarg/sphi3 >*/
    dsrgdb = sarg * -.5 * sarg / sphi3;
/*<       d2sdb2=-dsrgdb*sarg/sphi3 >*/
    d2sdb2 = -dsrgdb * sarg / sphi3;
/*<   120 if ((vds-vbs).lt.0.0d0) go to 130 >*/
L120:
    if (*vds - *vbs < 0.) {
	goto L130;
    }
/*<       barg=dsqrt(phi+vds-vbs) >*/
    barg = sqrt(mosarg_1.phi + *vds - *vbs);
/*<       dbrgdb=-0.5d0/barg >*/
    dbrgdb = -.5 / barg;
/*<       d2bdb2=+0.5d0*dbrgdb/(phi+vds-vbs) >*/
    d2bdb2 = dbrgdb * .5 / (mosarg_1.phi + *vds - *vbs);
/*<       go to 200 >*/
    goto L200;
/*<   130 barg=sphi/(1.0d0+0.5d0*(vbs-vds)/phi) >*/
L130:
    barg = sphi / ((*vbs - *vds) * .5 / mosarg_1.phi + 1.);
/*<       dbrgdb=-0.5d0*barg*barg/sphi3 >*/
    dbrgdb = barg * -.5 * barg / sphi3;
/*<       d2bdb2=-dbrgdb*barg/sphi3 >*/
    d2bdb2 = -dbrgdb * barg / sphi3;

/*  calculate threshold voltage (von) */
/*     narrow-channel effect */

/*<   200 factor=0.125d0*fnarrw*twopi*epssil/cox*xl >*/
L200:
    factor = mosarg_1.fnarrw * .125 * knstnt_1.twopi * knstnt_1.epssil / 
	    mosarg_1.cox * mosarg_1.xl;
/*<       eta=1.0d0+factor >*/
    eta = factor + 1.;
/*<       vbin=vbi+factor*(phi-vbs) >*/
    vbin = mosarg_1.vbi + factor * (mosarg_1.phi - *vbs);
/*<       if (gamma.le.0.0d0) go to 215 >*/
    if (mosarg_1.gamma <= 0.) {
	goto L215;
    }
/*<       if (xnsub.le.0.0d0) go to 215 >*/
    if (mosarg_1.xnsub <= 0.) {
	goto L215;
    }
/*<       xwd=xd*barg >*/
    xwd = mosarg_1.xd * barg;
/*<       xws=xd*sarg >*/
    xws = mosarg_1.xd * sarg;

/*     short-channel effect with vds .ne. 0.0d0 */

/*<       argss=0.0d0 >*/
    argss = 0.;
/*<       argsd=0.0d0 >*/
    argsd = 0.;
/*<       dbargs=0.0d0 >*/
    dbargs = 0.;
/*<       dbargd=0.0d0 >*/
    dbargd = 0.;
/*<       dgdvds=0.0d0 >*/
    dgdvds = 0.;
/*<       dgddb2=0.0d0 >*/
    dgddb2 = 0.;
/*<       if (xj.le.0.0d0) go to 205 >*/
    if (mosarg_1.xj <= 0.) {
	goto L205;
    }
/*<       argxs=1.0d0+2.0d0*xws/xj >*/
    argxs = xws * 2. / mosarg_1.xj + 1.;
/*<       args=dsqrt(argxs) >*/
    args = sqrt(argxs);
/*<       argss=0.5d0*xj/xl*(args-1.0d0) >*/
    argss = mosarg_1.xj * .5 / mosarg_1.xl * (args - 1.);
/*<       argxd=1.0d0+2.0d0*xwd/xj >*/
    argxd = xwd * 2. / mosarg_1.xj + 1.;
/*<       argd=dsqrt(argxd) >*/
    argd = sqrt(argxd);
/*<       argsd=0.5d0*xj/xl*(argd-1.0d0) >*/
    argsd = mosarg_1.xj * .5 / mosarg_1.xl * (argd - 1.);
/*<   205 gamasd=gamma*(1.0d0-argss-argsd) >*/
L205:
    gamasd = mosarg_1.gamma * (1. - argss - argsd);
/*<       gamass=gamma*(1.0d0-2.0d0*argss) >*/
    gamass = mosarg_1.gamma * (1. - argss * 2.);
/*<       dbxwd=xd*dbrgdb >*/
    dbxwd = mosarg_1.xd * dbrgdb;
/*<       dbxws=xd*dsrgdb >*/
    dbxws = mosarg_1.xd * dsrgdb;
/*<       if (xj.le.0.0d0) go to 210 >*/
    if (mosarg_1.xj <= 0.) {
	goto L210;
    }
/*<       dbargs=0.5d0/xl*dbxws/args >*/
    dbargs = .5 / mosarg_1.xl * dbxws / args;
/*<       dbargd=0.5d0/xl*dbxwd/argd >*/
    dbargd = .5 / mosarg_1.xl * dbxwd / argd;
/*<       dasdb2=-xd*( d2sdb2+dsrgdb*dsrgdb*xd/(xj*argxs) )/(xl*args) >*/
    dasdb2 = -mosarg_1.xd * (d2sdb2 + dsrgdb * dsrgdb * mosarg_1.xd / (
	    mosarg_1.xj * argxs)) / (mosarg_1.xl * args);
/*<       daddb2=-xd*( d2bdb2+dbrgdb*dbrgdb*xd/(xj*argxd) )/(xl*argd) >*/
    daddb2 = -mosarg_1.xd * (d2bdb2 + dbrgdb * dbrgdb * mosarg_1.xd / (
	    mosarg_1.xj * argxd)) / (mosarg_1.xl * argd);
/*<       dgddb2=-0.5d0*gamma*(dasdb2+daddb2) >*/
    dgddb2 = mosarg_1.gamma * -.5 * (dasdb2 + daddb2);
/*<   210 dgddvb=-gamma*(dbargs+dbargd) >*/
L210:
    dgddvb = -mosarg_1.gamma * (dbargs + dbargd);
/*<       dgsdvb=-2.0d0*gamma*dbargs >*/
    dgsdvb = mosarg_1.gamma * -2. * dbargs;
/*<       if (xj.le.0.0d0) go to 220 >*/
    if (mosarg_1.xj <= 0.) {
	goto L220;
    }
/*<       ddxwd=-dbxwd >*/
    ddxwd = -dbxwd;
/*<       dgdvds=-gamma*0.5d0/xl*ddxwd/argd >*/
    dgdvds = -mosarg_1.gamma * .5 / mosarg_1.xl * ddxwd / argd;
/*<       go to 220 >*/
    goto L220;
/*<   215 gamasd=gamma >*/
L215:
    gamasd = mosarg_1.gamma;
/*<       gamass=gamma >*/
    gamass = mosarg_1.gamma;
/*<       gammad=gamma >*/
    gammad = mosarg_1.gamma;
/*<       dgddvb=0.0d0 >*/
    dgddvb = 0.;
/*<       dgsdvb=0.0d0 >*/
    dgsdvb = 0.;
/*<       dgdvds=0.0d0 >*/
    dgdvds = 0.;
/*<       dgddb2=0.0d0 >*/
    dgddb2 = 0.;
/*<   220 von=vbin+gamasd*sarg >*/
L220:
    mosarg_1.von = vbin + gamasd * sarg;
/*     write(iofile,221) von,vbin,vbi,gamasd,argss,argsd,xj */
/*<   221 format ('0msg1:'/1p7d10.2) >*/
/* L221: */
/*<       vth=von >*/
    vth = mosarg_1.von;
/*<       vdsat=0.0d0 >*/
    mosarg_1.vdsat = 0.;
/*<   225 if (xnfs.eq.0.0d0.or.cox.eq.0.0d0) go to 230 >*/
/* L225: */
    if (mosarg_1.xnfs == 0. || mosarg_1.cox == 0.) {
	goto L230;
    }
/*<       cfs=charge*xnfs >*/
    cfs = knstnt_1.charge * mosarg_1.xnfs;
/*<       cdonco=-(gamasd*dsrgdb+dgddvb*sarg)+factor >*/
    cdonco = -(gamasd * dsrgdb + dgddvb * sarg) + factor;
/*<       xn=1.0d0+cfs/cox*xw*xl+cdonco >*/
    xn = cfs / mosarg_1.cox * mosarg_1.xw * mosarg_1.xl + 1. + cdonco;
/*<       von=von+vt*xn >*/
    mosarg_1.von += status_1.vt * xn;
/*     write (iofile,226) von,cdonco,xn,cfs,xd */
/*<   226 format(' msg2:'/1p6d10.2) >*/
/* L226: */
/*<       argg=1.0d0/(vt*xn) >*/
    argg = 1. / (status_1.vt * xn);
/*<       vgst=vgs-von >*/
    vgst = *vgs - mosarg_1.von;
/*<       go to 300 >*/
    goto L300;
/*<   230 vgst=vgs-von >*/
L230:
    vgst = *vgs - mosarg_1.von;
/*<       if (vgs.gt.von) go to 300 >*/
    if (*vgs > mosarg_1.von) {
	goto L300;
    }

/*  cutoff region */

/*<       gds=0.0d0 >*/
    *gds = 0.;
/*<       go to 1050 >*/
    goto L1050;

/*  compute some more useful quantities */

/*<   300 sarg3=sarg*sarg*sarg >*/
L300:
    sarg3 = sarg * sarg * sarg;
/*<       sbiarg=dsqrt(phib) >*/
    sbiarg = sqrt(mosarg_1.phib);
/*<       gammad=gamasd >*/
    gammad = gamasd;
/*<       dgdvbs=dgddvb >*/
    dgdvbs = dgddvb;
/*<       body=barg*barg*barg-sarg3 >*/
    body = barg * barg * barg - sarg3;
/*<       gdbdv=2.0d0*gammad*(barg*barg*dbrgdb-sarg*sarg*dsrgdb) >*/
    gdbdv = gammad * 2. * (barg * barg * dbrgdb - sarg * sarg * dsrgdb);
/*<       dodvbs=-factor+dgdvbs*sarg+gammad*dsrgdb >*/
    dodvbs = -factor + dgdvbs * sarg + gammad * dsrgdb;
/*<       if (xnfs.eq.0.0d0) go to 400 >*/
    if (mosarg_1.xnfs == 0.) {
	goto L400;
    }
/*<       if (cox.eq.0.0d0) go to 410 >*/
    if (mosarg_1.cox == 0.) {
	goto L410;
    }
/*<       dxndvb=2.0d0*dgdvbs*dsrgdb+gammad*d2sdb2+dgddb2*sarg >*/
    dxndvb = dgdvbs * 2. * dsrgdb + gammad * d2sdb2 + dgddb2 * sarg;
/*<       dodvbs=dodvbs+vt*dxndvb >*/
    dodvbs += status_1.vt * dxndvb;
/*<       dxndvd=dgdvds*dsrgdb >*/
    dxndvd = dgdvds * dsrgdb;
/*<       dodvds=dgdvds*sarg+vt*dxndvd >*/
    dodvds = dgdvds * sarg + status_1.vt * dxndvd;

/*  evaluate effective mobility and its derivatives */

/*<   400 if (cox.le.0.0d0) go to 410 >*/
L400:
    if (mosarg_1.cox <= 0.) {
	goto L410;
    }
/*<       udenom=vgst >*/
    udenom = vgst;
/*<       if (udenom.le.vbp) go to 410 >*/
    if (udenom <= mosarg_1.vbp) {
	goto L410;
    }
/*<       ufact=dexp(uexp*dlog(vbp/udenom)) >*/
    ufact = exp(mosarg_1.uexp * log(mosarg_1.vbp / udenom));
/*<       ueff=uo*ufact >*/
    ueff = mosarg_1.uo * ufact;
/*<       dudvgs=-ufact*uexp/udenom >*/
    dudvgs = -ufact * mosarg_1.uexp / udenom;
/*<       dudvds=0.0d0 >*/
    dudvds = 0.;
/*<       dudvbs=uexp*ufact*dodvbs/vgst >*/
    dudvbs = mosarg_1.uexp * ufact * dodvbs / vgst;
/*<       go to 500 >*/
    goto L500;
/*<   410 ufact=1.0d0 >*/
L410:
    ufact = 1.;
/*<       ueff=uo >*/
    ueff = mosarg_1.uo;
/*<       dudvgs=0.0d0 >*/
    dudvgs = 0.;
/*<       dudvds=0.0d0 >*/
    dudvds = 0.;
/*<       dudvbs=0.0d0 >*/
    dudvbs = 0.;

/*     evaluate saturation voltage and its derivatives according to */
/*     grove-frohman equation */

/*<   500 vgsx=vgs >*/
L500:
    vgsx = *vgs;
/*<       gammad=gamasd/eta >*/
    gammad = gamasd / eta;
/*<       dgdvbs=dgddvb >*/
    dgdvbs = dgddvb;
/*<       if (xnfs.ne.0.0d0.and.cox.ne.0.0d0) >*/
/*<      1   vgsx=dmax1(vgs,von) >*/
    if (mosarg_1.xnfs != 0. && mosarg_1.cox != 0.) {
	vgsx = max(*vgs,mosarg_1.von);
    }
/*<   505 if (gammad.le.0.0d0) go to 535 >*/
/* L505: */
    if (gammad <= 0.) {
	goto L535;
    }
/*<       gammd2=gammad*gammad >*/
    gammd2 = gammad * gammad;
/*<       argv=(vgsx-vbin)/eta+phi-vbs >*/
    argv = (vgsx - vbin) / eta + mosarg_1.phi - *vbs;
/*<       if (argv.le.0.0d0) go to 540 >*/
    if (argv <= 0.) {
	goto L540;
    }
/*<       arg=dsqrt(1.0d0+4.0d0*argv/gammd2) >*/
    arg = sqrt(argv * 4. / gammd2 + 1.);
/*<       vdsat=(vgsx-vbin)/eta+gammd2*(1.0d0-arg)/2.0d0 >*/
    mosarg_1.vdsat = (vgsx - vbin) / eta + gammd2 * (1. - arg) / 2.;
/*<       vdsat=dmax1(vdsat,0.0d0) >*/
    mosarg_1.vdsat = max(mosarg_1.vdsat,0.);
/*<   510 if (icharg.eq.0) go to 530 >*/
/* L510: */
    if (icharg == 0) {
	goto L530;
    }
/*<       arg1=gammd2/(eta*eta) >*/
    arg1 = gammd2 / (eta * eta);
/*<       arg2=vds-0.5d0*arg1 >*/
    arg2 = *vds - arg1 * .5;
/*<       argsq=(arg2+0.5d0*arg1+phi-vbs)*arg1 >*/
    argsq = (arg2 + arg1 * .5 + mosarg_1.phi - *vbs) * arg1;
/*<       if (argsq.ge.0.0d0) go to 515 >*/
    if (argsq >= 0.) {
	goto L515;
    }
/*<       vpof=vth >*/
    vpof = vth;
/*<       go to 520 >*/
    goto L520;
/*<   515 vpof=vbin+eta*(arg2+0.5d0*arg1+dsqrt(argsq)) >*/
L515:
    vpof = vbin + eta * (arg2 + arg1 * .5 + sqrt(argsq));
/*<   520 argv1=(vpof-vbin)/eta+phi-vbs >*/
L520:
    argv1 = (vpof - vbin) / eta + mosarg_1.phi - *vbs;
/*<       if (argv1.gt.0.0d0) go to 525 >*/
    if (argv1 > 0.) {
	goto L525;
    }
/*<       vdsat1=0.0d0 >*/
    vdsat1 = 0.;
/*<       go to 530 >*/
    goto L530;
/*<   525 arg1=dsqrt(1.0d0+4.0d0*argv1/gammd2) >*/
L525:
    arg1 = sqrt(argv1 * 4. / gammd2 + 1.);
/*<       vdsat1=(vpof-vbin)/eta+gammd2*(1.0d0-arg1)/2.0d0 >*/
    vdsat1 = (vpof - vbin) / eta + gammd2 * (1. - arg1) / 2.;
/*<       vdsat1=dmax1(vdsat1,0.0d0) >*/
    vdsat1 = max(vdsat1,0.);
/*<   530 dsdvgs=(1.0d0-1.0d0/arg)/eta >*/
L530:
    dsdvgs = (1. - 1. / arg) / eta;
/*<       dsdvbs=(gammad*(1.0d0-arg)+2.0d0*argv/(gammad*arg))/eta*dgdvbs+ >*/
/*<      1       1.0d0/arg+factor*dsdvgs >*/
    dsdvbs = (gammad * (1. - arg) + argv * 2. / (gammad * arg)) / eta * 
	    dgdvbs + 1. / arg + factor * dsdvgs;
/*<       go to 545 >*/
    goto L545;
/*<   535 vdsat=dmax1((vgsx-vbin)/eta,0.0d0) >*/
L535:
/* Computing MAX */
    d_1 = (vgsx - vbin) / eta;
    mosarg_1.vdsat = max(0.,d_1);
/*<       vpof=dmax1((eta*vds+vbin),0.0d0) >*/
/* Computing MAX */
    d_1 = eta * *vds + vbin;
    vpof = max(0.,d_1);
/*<       vdsat1=dmax1((vpof-vbin)/eta,0.0d0) >*/
/* Computing MAX */
    d_1 = (vpof - vbin) / eta;
    vdsat1 = max(0.,d_1);
/*<       dsdvgs=1.0d0 >*/
    dsdvgs = 1.;
/*<       dsdvbs=0.0d0 >*/
    dsdvbs = 0.;
/*<       go to 545 >*/
    goto L545;
/*<   540 vdsat=0.0d0 >*/
L540:
    mosarg_1.vdsat = 0.;
/*<       vpof=vth >*/
    vpof = vth;
/*<       vdsat1=0.0d0 >*/
    vdsat1 = 0.;
/*<       dsdvgs=0.0d0 >*/
    dsdvgs = 0.;
/*<       dsdvbs=0.0d0 >*/
    dsdvbs = 0.;

/*     store vdsat as above in vpof (pinch-off) */

/*<   545 if (vmax.le.0.0d0) go to 600 >*/
L545:
    if (mosarg_1.vmax <= 0.) {
	goto L600;
    }

/*     evaluate saturation voltage and its derivatives according to */
/*     baum's theory of scattering velocity saturation */

/*<       gammd2=gammad*gammad >*/
    gammd2 = gammad * gammad;
/*<       v1=(vgsx-vbin)/eta+phi-vbs >*/
    v1 = (vgsx - vbin) / eta + mosarg_1.phi - *vbs;
/*<       v2=phi-vbs >*/
    v2 = mosarg_1.phi - *vbs;
/*<       xv=vmax*xl/ueff >*/
    xv = mosarg_1.vmax * mosarg_1.xl / ueff;
/*<       a1=gammad/0.75d0 >*/
    a1 = gammad / .75;
/*<       b1=-2.0d0*(v1+xv) >*/
    b1 = (v1 + xv) * -2.;
/*<       c1=-2.0d0*gammad*xv >*/
    c1 = gammad * -2. * xv;
/*<       d1=2.0d0*v1*(v2+xv)-v2*v2-4.0d0/3.0d0*gammad*sarg3 >*/
    d1 = v1 * 2. * (v2 + xv) - v2 * v2 - gammad * 1.3333333333333333 * sarg3;
/*<       a=-b1 >*/
    a = -b1;
/*<       b=a1*c1-4.0d0*d1 >*/
    b = a1 * c1 - d1 * 4.;
/*<       c=-d1*(a1*a1-4.0d0*b1)-c1*c1 >*/
    c = -d1 * (a1 * a1 - b1 * 4.) - c1 * c1;
/*<       r=-a*a/3.0d0+b >*/
    r = -a * a / 3. + b;
/*<       s=2.0d0*a*a*a/27.0d0-a*b/3.0d0+c >*/
    s = a * 2. * a * a / 27. - a * b / 3. + c;
/*<       r3=r*r*r >*/
    r3 = r * r * r;
/*<       s2=s*s >*/
    s2 = s * s;
/*<       p=s2/4.0d0+r3/27.0d0 >*/
    p = s2 / 4. + r3 / 27.;
/*<       p0=dabs(p) >*/
    p0 = abs(p);
/*<       p2=dsqrt(p0) >*/
    p2 = sqrt(p0);
/*<       if (p.ge.0.0d0) go to 550 >*/
    if (p >= 0.) {
	goto L550;
    }
/*<       ro=dsqrt(s2/4.0d0+p0) >*/
    ro = sqrt(s2 / 4. + p0);
/*<       ro=dlog(ro)/3.0d0 >*/
    ro = log(ro) / 3.;
/*<       ro=dexp(ro) >*/
    ro = exp(ro);
/*<       fi=datan(-2.0d0*p2/s) >*/
    fi = atan(p2 * -2. / s);
/*<       y3=2.0d0*ro*dcos(fi/3.0d0)-a/3.0d0 >*/
    y3 = ro * 2. * cos(fi / 3.) - a / 3.;
/*<       go to 560 >*/
    goto L560;
/*<   550 p3=dexp(dlog(dabs(-s/2.0d0+p2))/3.0d0) >*/
L550:
    p3 = exp(log((d_1 = -s / 2. + p2, abs(d_1))) / 3.);
/*<       p4=dexp(dlog(dabs(-s/2.0d0-p2))/3.0d0) >*/
    p4 = exp(log((d_1 = -s / 2. - p2, abs(d_1))) / 3.);
/*<       y3=p3+p4-a/3.0d0 >*/
    y3 = p3 + p4 - a / 3.;
/*<   560 iknt=0 >*/
L560:
    iknt = 0;
/*<       a3=dsqrt(a1*a1/4.0d0-b1+y3) >*/
    a3 = sqrt(a1 * a1 / 4. - b1 + y3);
/*<       b3=dsqrt(y3*y3/4.0d0-d1) >*/
    b3 = sqrt(y3 * y3 / 4. - d1);
/*<       do 570 i=1,4 >*/
    for (i = 1; i <= 4; ++i) {
/*<       a4(i)=a1/2.0d0+sig1(i)*a3 >*/
	a4[i - 1] = a1 / 2. + sig1[i - 1] * a3;
/*<       b4(i)=y3/2.0d0+sig2(i)*b3 >*/
	b4[i - 1] = y3 / 2. + sig2[i - 1] * b3;
/*<       delta4=a4(i)*a4(i)/4.0d0-b4(i) >*/
	delta4 = a4[i - 1] * a4[i - 1] / 4. - b4[i - 1];
/*<       if (delta4.lt.0.0d0) go to 570 >*/
	if (delta4 < 0.) {
	    goto L570;
	}
/*<       iknt=iknt+1 >*/
	++iknt;
/*<       x4(iknt)=-a4(i)/2.0d0+dsqrt(delta4) >*/
	x4[iknt - 1] = -a4[i - 1] / 2. + sqrt(delta4);
/*<       iknt=iknt+1 >*/
	++iknt;
/*<       x4(iknt)=-a4(i)/2.0d0-dsqrt(delta4) >*/
	x4[iknt - 1] = -a4[i - 1] / 2. - sqrt(delta4);
/*<   570 continue >*/
L570:
    ;}
/*<       jknt=0 >*/
    jknt = 0;
/*<       do 580 j=1,iknt >*/
    i_1 = iknt;
    for (j = 1; j <= i_1; ++j) {
/*<       if (x4(j).le.0.0d0) go to 580 >*/
	if (x4[j - 1] <= 0.) {
	    goto L580;
	}
/*<       poly4(j)=x4(j)*x4(j)*x4(j)*x4(j)+a1*x4(j)*x4(j)*x4(j) >*/
	poly4[j - 1] = x4[j - 1] * x4[j - 1] * x4[j - 1] * x4[j - 1] + a1 * 
		x4[j - 1] * x4[j - 1] * x4[j - 1];
/*<       poly4(j)=poly4(j)+b1*x4(j)*x4(j)+c1*x4(j)+d1 >*/
	poly4[j - 1] = poly4[j - 1] + b1 * x4[j - 1] * x4[j - 1] + c1 * x4[j 
		- 1] + d1;
/*<       if (dabs(poly4(j)).gt.1.0d-6) go to 580 >*/
	if ((d_1 = poly4[j - 1], abs(d_1)) > 1e-6) {
	    goto L580;
	}
/*<       jknt=jknt+1 >*/
	++jknt;
/*<       if (jknt.gt.1) go to 575 >*/
	if (jknt > 1) {
	    goto L575;
	}
/*<       xvalid=x4(j) >*/
	xvalid = x4[j - 1];
/*<   575 if (x4(j).gt.xvalid) go to 580 >*/
L575:
	if (x4[j - 1] > xvalid) {
	    goto L580;
	}
/*<       xvalid=x4(j) >*/
	xvalid = x4[j - 1];
/*<   580 continue >*/
L580:
    ;}
/*<       if (jknt.gt.0) go to 590 >*/
    if (jknt > 0) {
	goto L590;
    }
/*<       ivmflg=ivmflg+1 >*/
    ++status_1.ivmflg;
/*<       go to 600 >*/
    goto L600;
/*<   590 vdsat=xvalid*xvalid+vbs-phi >*/
L590:
    mosarg_1.vdsat = xvalid * xvalid + *vbs - mosarg_1.phi;

/*  evaluate effective channel length and its derivatives */

/*<   600 if (vds.eq.0.0d0) go to 610 >*/
L600:
    if (*vds == 0.) {
	goto L610;
    }
/*<       gammad=gamasd >*/
    gammad = gamasd;
/*<       if ((vbs-vdsat).gt.0.0d0) go to 601 >*/
    if (*vbs - mosarg_1.vdsat > 0.) {
	goto L601;
    }
/*<       bsarg=dsqrt(vdsat-vbs+phi) >*/
    bsarg = sqrt(mosarg_1.vdsat - *vbs + mosarg_1.phi);
/*<       dbsrdb=-0.5d0/bsarg >*/
    dbsrdb = -.5 / bsarg;
/*<       go to 602 >*/
    goto L602;
/*<   601 bsarg=sphi/(1.0d0+0.5d0*(vbs-vdsat)/phi) >*/
L601:
    bsarg = sphi / ((*vbs - mosarg_1.vdsat) * .5 / mosarg_1.phi + 1.);
/*<       dbsrdb=-0.5d0*bsarg*bsarg/sphi3 >*/
    dbsrdb = bsarg * -.5 * bsarg / sphi3;
/*<   602 bodys=bsarg*bsarg*bsarg-sarg3 >*/
L602:
    bodys = bsarg * bsarg * bsarg - sarg3;
/*<       gdbdvs=2.0d0*gammad*(bsarg*bsarg*dbsrdb-sarg*sarg*dsrgdb) >*/
    gdbdvs = gammad * 2. * (bsarg * bsarg * dbsrdb - sarg * sarg * dsrgdb);
/*<       if (vmax.gt.0.0d0) go to 603 >*/
    if (mosarg_1.vmax > 0.) {
	goto L603;
    }
/*<       if (xnsub.eq.0.0d0) go to 610 >*/
    if (mosarg_1.xnsub == 0.) {
	goto L610;
    }
/*<       if (xlamda.gt.0.0d0) go to 610 >*/
    if (mosarg_1.xlamda > 0.) {
	goto L610;
    }
/*<       argv=(vds-vdsat)/4.0d0 >*/
    argv = (*vds - mosarg_1.vdsat) / 4.;
/*<       sargv=dsqrt(1.0d0+argv*argv) >*/
    sargv = sqrt(argv * argv + 1.);
/*<       arg=dsqrt(argv+sargv) >*/
    arg = sqrt(argv + sargv);
/*<       xlfact=xd/(xl*vds) >*/
    xlfact = mosarg_1.xd / (mosarg_1.xl * *vds);
/*<       xlamda=xlfact*arg >*/
    mosarg_1.xlamda = xlfact * arg;
/*<       dldsat=vds*xlfact*arg/(8.0d0*sargv) >*/
    dldsat = *vds * xlfact * arg / (sargv * 8.);
/*<       go to 605 >*/
    goto L605;
/*<   603 argv=(vgsx-vbin)/eta-vdsat >*/
L603:
    argv = (vgsx - vbin) / eta - mosarg_1.vdsat;
/*<       xdv=xd/dsqrt(xneff) >*/
    xdv = mosarg_1.xd / sqrt(mosarg_1.xneff);
/*<       xlv=vmax*xdv/(2.0d0*ueff) >*/
    xlv = mosarg_1.vmax * xdv / (ueff * 2.);
/*<       vqchan=argv-gammad*bsarg >*/
    vqchan = argv - gammad * bsarg;
/*<       dqdsat=-1.0d0+gammad*dbsrdb >*/
    dqdsat = gammad * dbsrdb - 1.;
/*<       vl=vmax*xl >*/
    vl = mosarg_1.vmax * mosarg_1.xl;
/*<       dfunds=vl*dqdsat-ueff*vqchan >*/
    dfunds = vl * dqdsat - ueff * vqchan;
/*<       dfundg=(vl-ueff*vdsat)/eta >*/
    dfundg = (vl - ueff * mosarg_1.vdsat) / eta;
/*<       dfundb=-vl*(1.0d0+dqdsat-factor/eta)+ >*/
/*<      1        ueff*(gdbdvs-dgdvbs*bodys/1.5d0)/eta >*/
    dfundb = -vl * (dqdsat + 1. - factor / eta) + ueff * (gdbdvs - dgdvbs * 
	    bodys / 1.5) / eta;
/*<       dsdvgs=-dfundg/dfunds >*/
    dsdvgs = -dfundg / dfunds;
/*<       dsdvbs=-dfundb/dfunds >*/
    dsdvbs = -dfundb / dfunds;
/*<       if (xnsub.eq.0.0d0) go to 610 >*/
    if (mosarg_1.xnsub == 0.) {
	goto L610;
    }
/*<       if (xlamda.gt.0.0d0) go to 610 >*/
    if (mosarg_1.xlamda > 0.) {
	goto L610;
    }
/*<       argv=dmax1(vds-vdsat,0.0d0) >*/
/* Computing MAX */
    d_1 = *vds - mosarg_1.vdsat;
    argv = max(0.,d_1);
/*<       xls=dsqrt(xlv*xlv+argv) >*/
    xls = sqrt(xlv * xlv + argv);
/*<       dldsat=xdv/(2.0d0*xls) >*/
    dldsat = xdv / (xls * 2.);
/*<       xlfact=xdv/(xl*vds) >*/
    xlfact = xdv / (mosarg_1.xl * *vds);
/*<       xlamda=xlfact*(xls-xlv) >*/
    mosarg_1.xlamda = xlfact * (xls - xlv);
/*<       dldsat=dldsat/xl >*/
    dldsat /= mosarg_1.xl;
/*<   605 dldvgs=dldsat*dsdvgs >*/
L605:
    dldvgs = dldsat * dsdvgs;
/*<       dldvds=-xlamda+dldsat >*/
    dldvds = -mosarg_1.xlamda + dldsat;
/*<       dldvbs=dldsat*dsdvbs >*/
    dldvbs = dldsat * dsdvbs;
/*<       go to 620 >*/
    goto L620;
/*<   610 dldvgs=0.0d0 >*/
L610:
    dldvgs = 0.;
/*<       dldvds=0.0d0 >*/
    dldvds = 0.;
/*<       dldvbs=0.0d0 >*/
    dldvbs = 0.;

/*     limit channel shortening at punch-through */

/*<   620 xwb=xd*sbiarg >*/
L620:
    xwb = mosarg_1.xd * sbiarg;
/*<       xld=xl-xwb >*/
    mosarg_1.xld = mosarg_1.xl - xwb;
/*<       clfact=1.0d0-xlamda*vds >*/
    clfact = 1. - mosarg_1.xlamda * *vds;
/*<       dldvds=-xlamda-dldvds >*/
    dldvds = -mosarg_1.xlamda - dldvds;
/*<       xleff=xl*clfact >*/
    xleff = mosarg_1.xl * clfact;
/*<       deltal=xlamda*vds*xl >*/
    deltal = mosarg_1.xlamda * *vds * mosarg_1.xl;
/*<       if (xnsub.eq.0.0d0) xwb=0.25d-6 >*/
    if (mosarg_1.xnsub == 0.) {
	xwb = 2.5e-7;
    }
/*<       if (xleff.ge.xwb) go to 700 >*/
    if (xleff >= xwb) {
	goto L700;
    }
/*<       xleff=xwb/(1.0d0+(deltal-xld)/xwb) >*/
    xleff = xwb / ((deltal - mosarg_1.xld) / xwb + 1.);
/*<       clfact=xleff/xl >*/
    clfact = xleff / mosarg_1.xl;
/*<       dfact=xleff*xleff/(xwb*xwb) >*/
    dfact = xleff * xleff / (xwb * xwb);
/*<       dldvgs=dfact*dldvgs >*/
    dldvgs = dfact * dldvgs;
/*<       dldvds=dfact*dldvds >*/
    dldvds = dfact * dldvds;
/*<       dldvbs=dfact*dldvbs >*/
    dldvbs = dfact * dldvbs;

/*  evaluate effective beta (effective kp) */

/*<   700 beta1=beta*ufact/clfact >*/
L700:
    mosarg_1.beta1 = mosarg_1.beta * ufact / clfact;

/*  test for mode of operation and branch appropriately */

/*<       gammad=gamasd >*/
    gammad = gamasd;
/*<       dgdvbs=dgddvb >*/
    dgdvbs = dgddvb;
/*<       if (vds.gt.1.0d-8) go to 730 >*/
    if (*vds > 1e-8) {
	goto L730;
    }
/*<       if (vgs.gt.von) go to 720 >*/
    if (*vgs > mosarg_1.von) {
	goto L720;
    }
/*<       if ((xnfs.ne.0.0d0).and.(cox.ne.0.0d0)) go to 710 >*/
    if (mosarg_1.xnfs != 0. && mosarg_1.cox != 0.) {
	goto L710;
    }
/*<       gds=0.0d0 >*/
    *gds = 0.;
/*<       go to 1050 >*/
    goto L1050;

/*<   710 gds=beta1*(von-vbin-gammad*sarg)*dexp(argg*(vgs-von)) >*/
L710:
    *gds = mosarg_1.beta1 * (mosarg_1.von - vbin - gammad * sarg) * exp(argg *
	     (*vgs - mosarg_1.von));
/*<       go to 1050 >*/
    goto L1050;


/*<   720 gds=beta1*(vgs-vbin-gammad*sarg) >*/
L720:
    *gds = mosarg_1.beta1 * (*vgs - vbin - gammad * sarg);
/*<       go to 1050 >*/
    goto L1050;

/*<   730 if (vgs.gt.von) go to 900 >*/
L730:
    if (*vgs > mosarg_1.von) {
	goto L900;
    }

/*  subthreshold region */

/*<       if (vdsat.gt.0.0d0) go to 830 >*/
    if (mosarg_1.vdsat > 0.) {
	goto L830;
    }
/*<       gds=0.0d0 >*/
    *gds = 0.;
/*<       if (vgs.gt.vth) go to 1020 >*/
    if (*vgs > vth) {
	goto L1020;
    }
/*<       go to 1050 >*/
    goto L1050;
/*<   830 vdson=dmin1(vdsat,vds) >*/
L830:
    vdson = min(mosarg_1.vdsat,*vds);
/*<       if (vds.le.vdsat) go to 850 >*/
    if (*vds <= mosarg_1.vdsat) {
	goto L850;
    }
/*<       barg=bsarg >*/
    barg = bsarg;
/*<       dbrgdb=dbsrdb >*/
    dbrgdb = dbsrdb;
/*<       body=bodys >*/
    body = bodys;
/*<       gdbdv=gdbdvs >*/
    gdbdv = gdbdvs;
/*<   850 cdson=beta1*((von-vbin-eta*vdson*0.5d0)*vdson-gammad*body/1.5d0) >*/
L850:
    cdson = mosarg_1.beta1 * ((mosarg_1.von - vbin - eta * vdson * .5) * 
	    vdson - gammad * body / 1.5);
/*<       didvds=beta1*(von-vbin-eta*vdson-gammad*barg) >*/
    didvds = mosarg_1.beta1 * (mosarg_1.von - vbin - eta * vdson - gammad * 
	    barg);
/*<       gdson=-cdson*dldvds/clfact-beta1*dgdvds*body/1.5d0 >*/
    gdson = -cdson * dldvds / clfact - mosarg_1.beta1 * dgdvds * body / 1.5;
/*<       if (vds.lt.vdsat) gdson=gdson+didvds >*/
    if (*vds < mosarg_1.vdsat) {
	gdson += didvds;
    }
/*<       gbson=-cdson*dldvbs/clfact >*/
/*<      1      +beta1*(dodvbs*vdson+factor*vdson-dgdvbs*body/1.5d0-gdbdv) >*/
    gbson = -cdson * dldvbs / clfact + mosarg_1.beta1 * (dodvbs * vdson + 
	    factor * vdson - dgdvbs * body / 1.5 - gdbdv);
/*<       if (vds.gt.vdsat) gbson=gbson+didvds*dsdvbs >*/
    if (*vds > mosarg_1.vdsat) {
	gbson += didvds * dsdvbs;
    }
/*<       expg=dexp(argg*(vgs-von)) >*/
    expg = exp(argg * (*vgs - mosarg_1.von));
/*<       cdrain=cdson*expg >*/
    mosarg_1.cdrain = cdson * expg;
/*<       gmw=cdrain*argg >*/
    gmw = mosarg_1.cdrain * argg;
/*<       gm=gmw >*/
    *gm = gmw;
/*<       if (vds.gt.vdsat) gm=gmw+didvds*dsdvgs*expg >*/
    if (*vds > mosarg_1.vdsat) {
	*gm = gmw + didvds * dsdvgs * expg;
    }
/*<       gds=gdson*expg-gm*dodvds-gmw*(vgs-von)*dxndvd/xn >*/
    *gds = gdson * expg - *gm * dodvds - gmw * (*vgs - mosarg_1.von) * dxndvd 
	    / xn;
/*<       gmbs=gbson*expg-gm*dodvbs-gmw*(vgs-von)*dxndvb/xn >*/
    *gmbs = gbson * expg - *gm * dodvbs - gmw * (*vgs - mosarg_1.von) * 
	    dxndvb / xn;
/*<       go to 1020 >*/
    goto L1020;


/*<   900 if (vds.gt.vdsat) go to 1000 >*/
L900:
    if (*vds > mosarg_1.vdsat) {
	goto L1000;
    }

/*  linear region */

/*<       cdrain=beta1*((vgs-vbin-eta*vds/2.0d0)*vds-gammad*body/1.5d0) >*/
    mosarg_1.cdrain = mosarg_1.beta1 * ((*vgs - vbin - eta * *vds / 2.) * *
	    vds - gammad * body / 1.5);
/*<       arg=cdrain*(dudvgs/ufact-dldvgs/clfact) >*/
    arg = mosarg_1.cdrain * (dudvgs / ufact - dldvgs / clfact);
/*<       gm=arg+beta1*vds >*/
    *gm = arg + mosarg_1.beta1 * *vds;
/*<       arg=cdrain*(dudvds/ufact-dldvds/clfact) >*/
    arg = mosarg_1.cdrain * (dudvds / ufact - dldvds / clfact);
/*<       gds=arg+beta1*(vgs-vbin-eta*vds- >*/
/*<      1   gammad*barg-dgdvds*body/1.5d0) >*/
    *gds = arg + mosarg_1.beta1 * (*vgs - vbin - eta * *vds - gammad * barg - 
	    dgdvds * body / 1.5);
/*<       arg=cdrain*(dudvbs/ufact-dldvbs/clfact) >*/
    arg = mosarg_1.cdrain * (dudvbs / ufact - dldvbs / clfact);
/*<       gmbs=arg-beta1*(gdbdv+dgdvbs*body/1.5d0-factor*vds) >*/
    *gmbs = arg - mosarg_1.beta1 * (gdbdv + dgdvbs * body / 1.5 - factor * *
	    vds);
/*<       go to 1020 >*/
    goto L1020;

/*  saturation region */

/*<  1000 cdrain=beta1*((vgs-vbin-eta*vdsat/2.0d0)*vdsat-gammad*bodys/1.5d0) >*/
L1000:
    mosarg_1.cdrain = mosarg_1.beta1 * ((*vgs - vbin - eta * mosarg_1.vdsat / 
	    2.) * mosarg_1.vdsat - gammad * bodys / 1.5);
/*<       arg=cdrain*(dudvgs/ufact-dldvgs/clfact) >*/
    arg = mosarg_1.cdrain * (dudvgs / ufact - dldvgs / clfact);
/*<       gm=arg+beta1*vdsat+ >*/
/*<      1   beta1*(vgs-vbin-eta*vdsat-gammad*bsarg)*dsdvgs >*/
    *gm = arg + mosarg_1.beta1 * mosarg_1.vdsat + mosarg_1.beta1 * (*vgs - 
	    vbin - eta * mosarg_1.vdsat - gammad * bsarg) * dsdvgs;
/*<       gds=-cdrain*dldvds/clfact-beta1*dgdvds*bodys/1.5d0 >*/
    *gds = -mosarg_1.cdrain * dldvds / clfact - mosarg_1.beta1 * dgdvds * 
	    bodys / 1.5;
/*<       arg=cdrain*(dudvbs/ufact-dldvbs/clfact) >*/
    arg = mosarg_1.cdrain * (dudvbs / ufact - dldvbs / clfact);
/*<       gmbs=arg-beta1*(gdbdvs+dgdvbs*bodys/1.5d0-factor*vdsat)+ >*/
/*<      1     beta1*(vgs-vbin-eta*vdsat-gammad*bsarg)*dsdvbs >*/
    *gmbs = arg - mosarg_1.beta1 * (gdbdvs + dgdvbs * bodys / 1.5 - factor * 
	    mosarg_1.vdsat) + mosarg_1.beta1 * (*vgs - vbin - eta * 
	    mosarg_1.vdsat - gammad * bsarg) * dsdvbs;

/*     compute charges for "on" region */

/*<  1020 if (icharg.eq.0) go to 1500 >*/
L1020:
    if (icharg == 0) {
	goto L1500;
    }
/*<       if (vgs.le.vth) go to 1070 >*/
    if (*vgs <= vth) {
	goto L1070;
    }
/*<       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;

/*  finish special cases */

/*<  1050 cdrain=0.0d0 >*/
L1050:
    mosarg_1.cdrain = 0.;
/*<       gm=0.0d0 >*/
    *gm = 0.;
/*<       gmbs=0.0d0 >*/
    *gmbs = 0.;
/*<  1070 xqc=xqco >*/
L1070:
    mosarg_1.xqc = mosarg_1.xqco;
/*<       if (icharg.eq.0) go to 1500 >*/
    if (icharg == 0) {
	goto L1500;
    }
/*<       call mosq2(vds,vbs,vgs,vdsat,vth,vbin,gamasd,cox,phi, >*/
/*<      1   qg,qc,qb,cggb,cgdb,cgsb,cbgb,cbdb,cbsb) >*/
    mosq2_(vds, vbs, vgs, &mosarg_1.vdsat, &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;

/*  finished */

/*<  1500 qg=0.0d0 >*/
L1500:
    *qg = 0.;
/*<       qb=0.0d0 >*/
    *qb = 0.;
/*<       qc=0.0d0 >*/
    *qc = 0.;
/*<       qspof=0.0d0 >*/
    mosarg_1.qspof = 0.;
/*<  2000 return >*/
L2000:
    return 0;
/*<       end >*/
} /* moseq2_ */

