/*****************************************************************/
/* Berkeley BSIM3v3.2.0 & BSIM3v3.2.4 (default) Verilog-A model  */
/*****************************************************************/

//      UPDATED March 19 2004 
//      Contributed By:
//      Geoffrey Coram, Ph.D Senior CAD Engineer Analog Devices, Inc.


`define VOLTAGE_MAXDELTA 0.3

`include "discipline.h"

// Following line must be uncomment for using NQS charge model (NQSMOD=1)
//`define NQSMOD

//****** Physical constants ******//
`define EPSOX         3.453133e-11
`define KboQ          8.617087e-5 
`define EPSSI         1.03594e-10
`define Charge_q      1.60219e-19
`define CONSTvt0      0.02586419
`define CONSTroot2    1.41421356

//****** Mathematical constants and constants of limitation ******//
`define PI            3.141592654
`define EXP_THRESHOLD 34.0
`define MIN_EXP       1.713908431e-15
`define MAX_EXP       5.834617425e14

// //****** Constants for the model ******//
`define DELTA_1 0.02
`define DELTA_3 0.02
`define DELTA_4 0.02

module bsim3mos (drain, gate, source, bulk);
   inout drain, gate, source, bulk;
   electrical drain, gate, source, bulk;
   electrical drainp, sourcep; // internal nodes
`ifdef NQSMOD
   electrical q;               // NQS charge model node
`endif

   //****** Device Parameters ******//
   parameter real L      = 5.0e-6;
   parameter real W      = 5.0e-6;
   parameter real PS     = 0.0;
   parameter real PD     = 0.0;
   parameter real AS     = 0.0;
   parameter real AD     = 0.0;
   parameter real NRS    = 1.0;     // Number of source diffusion squares
   parameter real NRD    = 1.0;     // Number of drain diffusion squares
   parameter real NQSMOD = 0;       // Non-quasi-static model selector
   
   parameter real GMIN = 1e-12;
   parameter real TYPE = 1.0;

   // Versions can be 3.20 or 3.24
   // (BSIM3v3.2.0 or BSIM3v3.2.4)
   parameter real VERSION  = 3.24;
   parameter real PARAMCHK = 0;
   
   //****** Model Selectors/Controllers ******//
   parameter real MOBMOD   = 1;     // Mobility model selector
   parameter real CAPMOD   = 3;     // Capacitance model selector
   parameter real NOIMOD   = 1;
   parameter real BINUNIT  = 1;     // Bin  unit  selector
   
   parameter real TOX  = 150.0e-10; // Gate oxide thickness in meters
   parameter real TOXM = TOX;       // Gate oxide thickness used in extraction
 
   parameter real CDSC    = 2.4e-4;  // Drain/Source and channel coupling capacitance
   parameter real CDSCB   = 0.0;     // Body-bias dependence of cdsc
   parameter real CDSCD   = 0.0;     // Drain-bias dependence of cdsc
   parameter real CIT     = 0.0;     // Interface state capacitance
   parameter real NFACTOR = 1;       // Subthreshold swing Coefficient
   parameter real XJ      = 0.15e-6; // Junction depth in meters
   parameter real VSAT    = 8.0e4;   // Saturation velocity at tnom
   parameter real AT      = 3.3e4;   // Temperature coefficient of vsat
   parameter real A0      = 1.0;     // Non-uniform depletion width effect coefficient.
   parameter real AGS     = 0.0;     // Gate bias  coefficient of Abulk.
   parameter real A1      = 0.0;     // Non-saturation effect coefficient
   parameter real A2      = 1.0;     // Non-saturation effect coefficient
   parameter real KETA    = -0.047;  // Body-bias coefficient of non-uniform depletion width effect.
   parameter real NSUB    = -99.0;   // Substrate doping concentration
   parameter real NCH     = -99.0;   // Channel doping concentration
   parameter real NGATE   = 0;       // Poly-gate doping concentration
   parameter real GAMMA1  = -99.0;   // Vth body coefficient
   parameter real GAMMA2  = -99.0;   // Vth body coefficient
   parameter real VBX     = -99.0;   // Vth transition body Voltage
   parameter real VBM     = -3.0;    // Maximum body voltage

   parameter real XT    = -99.0;   // Doping depth
   parameter real K1    = -99.0;   // Bulk effect coefficient 1
   parameter real KT1   = -0.11;   // Temperature coefficient of Vth
   parameter real KT1L  = 0.0;     // Temperature coefficient of Vth
   parameter real KT2   = 0.022;   // Body-coefficient of kt1
   parameter real K2    = -99.0;   // Bulk effect coefficient 2
   parameter real K3    = 80.0;    // Narrow width effect coefficient
   parameter real K3B   = 0.0;     // Body effect coefficient of k3
   parameter real W0    = 2.5e-6;  // Narrow width effect parameter
   parameter real NLX   = 1.74e-7; // Lateral non-uniform doping effect
   parameter real DVT0  = 2.2;     // Short channel effect coeff. 0
   parameter real DVT1  = 0.53;    // Short channel effect coeff. 1
   parameter real DVT2  = -0.032;  // Short channel effect coeff. 2
   parameter real DVT0W = 0.0;     // Narrow Width coeff. 0
   parameter real DVT1W = 5.3e6;   // Narrow Width effect coeff. 1
   parameter real DVT2W = -0.032;  // Narrow Width effect coeff. 2
   parameter real DROUT = 0.56;    // DIBL coefficient of output resistance
   parameter real DSUB  = DROUT;   // DIBL coefficient in the subthreshold region
   parameter real VTHO  = -99.0;   // Threshold voltage
   parameter real VTH0  = -99.0;   // Threshold voltage
   
   parameter real UA    = 2.25e-9;   // Linear gate dependence of mobility
   parameter real UA1   = 4.31e-9;   // Temperature coefficient of ua
   parameter real UB    = 5.87e-19;  // Quadratic gate dependence of mobility
   parameter real UB1   = -7.61e-18; // Temperature coefficient of ub
   parameter real UC    = -99.0;     // Body-bias dependence of mobility
   parameter real UC1   = -99.0;     // Temperature coefficient of uc
   parameter real U0    =-99.0;     // -99.0 Low-field mobility at Tnom
   parameter real UTE   = -1.5;      // Temperature coefficient of mobility
   parameter real VOFF  = -0.08;     // Threshold voltage offset
   parameter real TNOM  = -99.0;     // Parameter measurement temperature in C degree
   parameter real CGSO  = -99.0;     // Gate-source overlap capacitance per width
   parameter real CGDO  = -99.0;     // Gate-drain overlap capacitance per width
   parameter real CGBO  = -99.0;     // Gate-bulk overlap capacitance per length
   parameter real XPART = 0.0;       // Channel charge partitioning
   parameter real ELM   = 5.0;       // Non-quasi-static Elmore Constant Parameter
   parameter real DELTA = 0.01;      // Effective Vds parameter
   parameter real RSH   = 0.0;       // Source-drain sheet resistance
   parameter real RDSW  = 0;         // Source-drain resistance per width
   
   parameter real PRWG     = 0.0;      // Gate-bias effect on parasitic resistance     
   parameter real PRWB     = 0.0;      // Body-effect on parasitic resistance     
   parameter real PRT      = 0.0;      // Temperature coefficient of parasitic resistance     
   parameter real ETA0     = 0.08;     // Subthreshold region DIBL coefficient
   parameter real ETAB     = -0.07;    // Subthreshold region DIBL coefficient
   parameter real PCLM     = 1.3;      // Channel length modulation Coefficient
   parameter real PDIBLC1  = 0.39;     // Drain-induced barrier lowering coefficient   
   parameter real PDIBLC2  = 0.0086;   // Drain-induced barrier lowering coefficient   
   parameter real PDIBLCB  = 0.0;      // Body-effect on drain-induced barrier lowering   
   parameter real PSCBE1   = 4.24e8;   // Substrate current body-effect coefficient   
   parameter real PSCBE2   = 1.0e-5;   // Substrate current body-effect coefficient   
   parameter real PVAG     = 0.0;      // Gate dependence of output resistance parameter
   parameter real JS       = 1.0E-4;   // Source/drain junction reverse saturation current density
   parameter real JSW      = 0.0;      // Sidewall junction reverse saturation current density
   parameter real PB       = 1.0;      // Source/drain junction built-in potential
   parameter real NJ       = 1.0;      // Source/drain junction emission coefficient
   parameter real XTI      = 3.0;      // Junction current temperature exponent
   parameter real MJ       = 0.5;      // Source/drain bottom junction capacitance grading coefficient
   parameter real PBSW     = 1.0;      // Source/drain sidewall junction capacitance built in potential
   parameter real MJSW     = 0.33;     // Source/drain sidewall junction capacitance grading coefficient
   parameter real PBSWG    = PBSW;     // Source/drain (gate side) sidewall junction capacitance built in potential
   parameter real MJSWG    = MJSW;     // Source/drain (gate side) sidewall junction capacitance grading coefficient
   parameter real CJ       = 5.0E-4;   // Source/drain bottom junction capacitance per unit area
   parameter real VFBCV    = -1.0;     // Flat Band Voltage parameter for capmod=0 only
   parameter real VFB      = -99.0;    // Flat Band Voltage
   parameter real CJSW     = 5.0E-10;  // Source/drain sidewall junction capacitance per unit periphery
   parameter real CJSWG    = CJSW;     // Source/drain (gate side) sidewall junction capacitance per unit width
   parameter real TPB      = 0.0;      // Temperature coefficient of pb
   parameter real TCJ      = 0.0;      // Temperature coefficient of cj
   parameter real TPBSW    = 0.0;      // Temperature coefficient of pbsw
   parameter real TCJSW    = 0.0;      // Temperature coefficient of cjsw
   parameter real TPBSWG   = 0.0;      // Temperature coefficient of pbswg
   parameter real TCJSWG   = 0.0;      // Temperature coefficient of cjswg
   parameter real ACDE     = 1.0;      // Exponential coefficient for finite charge thickness
   parameter real MOIN     = 15.0;     // Coefficient for gate-bias dependent surface potential
   parameter real NOFF     = 1.0;      // C-V turn-on/off parameter
   parameter real VOFFCV   = 0.0;      // C-V lateral-shift parameter
   parameter real LINT     = 0.0;      // Length reduction parameter
   parameter real LL       = 0.0;      // Length reduction parameter
   parameter real LLC      = LL;       // Length reduction parameter for CV
   parameter real LLN      = 1.0;      // Length reduction parameter
   parameter real LW       = 0.0;      // Length reduction parameter
   parameter real LWC      = LW;       // Length reduction parameter for CV
   parameter real LWN      = 1.0;      // Length reduction parameter
   parameter real LWL      = 0.0;      // Length reduction parameter
   parameter real LWLC     = LWL;      // Length reduction parameter for CV
   parameter real LMIN     = 0.0;      // Minimum length for the model
   parameter real LMAX     = 1.0;      // Maximum length for the model
   parameter real WR       = 1.0;      // Width dependence of rds
   parameter real WINT     = 0.0;      // Width reduction parameter
   parameter real DWG      = 0.0;      // Width reduction parameter
   parameter real DWB      = 0.0;      // Width reduction parameter
   parameter real WL       = 0.0;      // Width reduction parameter
   parameter real WLC      = WL;       // Width reduction parameter for CV
   parameter real WLN      = 1.0;      // Width reduction parameter
   parameter real WW       = 0.0;      // Width reduction parameter
   parameter real WWC      = WW;       // Width reduction parameter for CV
   parameter real WWN      = 1.0;      // Width reduction parameter
   parameter real WWL      = 0.0;      // Width reduction parameter
   parameter real WWLC     = WWL;      // Width reduction parameter for CV
   parameter real WMIN     = 0.0;      // Minimum width for the model
   parameter real WMAX     = 1.0;      // Maximum width for the model
   parameter real B0       = 0.0;      // Abulk narrow width parameter
   parameter real B1       = 0.0;      // Abulk narrow width parameter
   parameter real CGSL     = 0.0;      // New C-V model parameter
   parameter real CGDL     = 0.0;      // New C-V model parameter
   parameter real CKAPPA   = 0.6;      // New C-V model parameter
   parameter real CF       = -99.0;    // Fringe capacitance parameter
   parameter real CLC      = 0.1e-6;   // Vdsat parameter for C-V model
   parameter real CLE      = 0.6;      // Vdsat parameter for C-V model
   parameter real DWC      = WINT;     // Delta W for C-V model
   parameter real DLC      = -99.0;    // Delta L for C-V model
   parameter real ALPHA0   = 0.0;      // substrate current model parameter
   parameter real ALPHA1   = 0.0;      // substrate current model parameter
   parameter real BETA0    = 30.0;     // substrate current model parameter
   parameter real IJTH     = 0.1;      // Diode limiting current

   /*** Length dependance model parameters ***/
   parameter real LCDSC    = 0.0;   // Length dependence of cdsc
   parameter real LCDSCB   = 0.0;   // Length dependence of cdscb
   parameter real LCDSCD   = 0.0;   // Length dependence of cdscd
   parameter real LCIT     = 0.0;   // Length dependence of cit
   parameter real LNFACTOR = 0.0;   // Length dependence of nfactor
   parameter real LXJ      = 0.0;   // Length dependence of xj
   parameter real LVSAT    = 0.0;   // Length dependence of vsat
   parameter real LAT      = 0.0;   // Length dependence of at
   parameter real LA0      = 0.0;   // Length dependence of a0 
   parameter real LAGS     = 0.0;   // Length dependence of ags 
   parameter real LA1      = 0.0;   // Length dependence of a1
   parameter real LA2      = 0.0;   // Length dependence of a2
   parameter real LKETA    = 0.0;   // Length dependence of keta
   parameter real LNSUB    = 0.0;   // Length dependence of nsub
   parameter real LNCH     = 0.0;   // Length dependence of nch
   parameter real LNGATE   = 0.0;   // Length dependence of ngate
   parameter real LGAMMA1  = -99.0; // Length dependence of gamma1
   parameter real LGAMMA2  = -99.0; // Length dependence of gamma2
   parameter real LVBX     = -99.0; // Length dependence of vbx
   parameter real LVBM     = 0.0;   // Length dependence of vbm
   parameter real LXT      = 0.0;   // Length dependence of xt
   parameter real LK1      = -99.0; // Length dependence of k1
   parameter real LKT1     = 0.0;   // Length dependence of kt1
   parameter real LKT1L    = 0.0;   // Length dependence of kt1l
   parameter real LKT2     = 0.0;   // Length dependence of kt2
   parameter real LK2      = -99.0; // Length dependence of k2
   parameter real LK3      = 0.0;   // Length dependence of k3
   parameter real LK3B     = 0.0;   // Length dependence of k3b
   parameter real LW0      = 0.0;   // Length dependence of w0
   parameter real LNLX     = 0.0;   // Length dependence of nlx
   parameter real LDVT0    = 0.0;   // Length dependence of dvt0
   parameter real LDVT1    = 0.0;   // Length dependence of dvt1
   parameter real LDVT2    = 0.0;   // Length dependence of dvt2
   parameter real LDVT0W   = 0.0;   // Length dependence of dvt0w
   parameter real LDVT1W   = 0.0;   // Length dependence of dvt1w
   parameter real LDVT2W   = 0.0;   // Length dependence of dvt2w
   parameter real LDROUT   = 0.0;   // Length dependence of drout
   parameter real LDSUB    = 0.0;   // Length dependence of dsub
   parameter real LVTH0    = 0.0;   // Length dependence of vto
   parameter real LVTHO    = 0.0;   // Length dependence of vto
   parameter real LUA      = 0.0;   // Length dependence of ua
   parameter real LUA1     = 0.0;   // Length dependence of ua1
   parameter real LUB      = 0.0;   // Length dependence of ub
   parameter real LUB1     = 0.0;   // Length dependence of ub1
   parameter real LUC      = 0.0;   // Length dependence of uc
   parameter real LUC1     = 0.0;   // Length dependence of uc1
   parameter real LU0      = 0.0;   // Length dependence of u0
   parameter real LUTE     = 0.0;   // Length dependence of ute
   parameter real LVOFF    = 0.0;   // Length dependence of voff
   parameter real LELM     = 0.0;   // Length dependence of elm
   parameter real LDELTA   = 0.0;   // Length dependence of delta
   parameter real LRDSW    = 0.0;   // Length dependence of rdsw     
   parameter real LPRWG    = 0.0;   // Length dependence of prwg     
   parameter real LPRWB    = 0.0;   // Length dependence of prwb 
   parameter real LPRT     = 0.0;   // Length dependence of prt     
   parameter real LETA0    = 0.0;   // Length dependence of eta0   
   parameter real LETAB    = 0.0;   // Length dependence of etab   
   parameter real LPCLM    = 0.0;   // Length dependence of pclm   
   parameter real LPDIBLC1 = 0.0;   // Length dependence of pdiblc1   
   parameter real LPDIBLC2 = 0.0;   // Length dependence of pdiblc2   
   parameter real LPDIBLCB = 0.0;   // Length dependence of pdiblcb   
   parameter real LPSCBE1  = 0.0;   // Length dependence of pscbe1   
   parameter real LPSCBE2  = 0.0;   // Length dependence of pscbe2   
   parameter real LPVAG    = 0.0;   // Length dependence of pvag   
   parameter real LWR      = 0.0;   // Length dependence of wr
   parameter real LDWG     = 0.0;   // Length dependence of dwg
   parameter real LDWB     = 0.0;   // Length dependence of dwb
   parameter real LB0      = 0.0;   // Length dependence of b0
   parameter real LB1      = 0.0;   // Length dependence of b1
   parameter real LCGSL    = 0.0;   // Length dependence of cgsl
   parameter real LCGDL    = 0.0;   // Length dependence of cgdl
   parameter real LCKAPPA  = 0.0;   // Length dependence of ckappa
   parameter real LCF      = 0.0;   // Length dependence of cf
   parameter real LCLC     = 0.0;   // Length dependence of clc
   parameter real LCLE     = 0.0;   // Length dependence of cle
   parameter real LALPHA0  = 0.0;   // Length dependence of alpha0
   parameter real LALPHA1  = 0.0;   // Length dependence of alpha1
   parameter real LBETA0   = 0.0;   // Length dependence of beta0
   parameter real LVFBCV   = 0.0;   // Length dependence of vfbcv
   parameter real LVFB     = 0.0;   // Length dependence of vfb
   parameter real LACDE    = 0.0;   // Length dependence of acde
   parameter real LMOIN    = 0.0;   // Length dependence of moin
   parameter real LNOFF    = 0.0;   // Length dependence of noff
   parameter real LVOFFCV  = 0.0;   // Length dependence of voffcv

   parameter real WCDSC    = 0.0;   // Width dependence of cdsc
   parameter real WCDSCB   = 0.0;   // Width dependence of cdscb  
   parameter real WCDSCD   = 0.0;   // Width dependence of cdscd  
   parameter real WCIT     = 0.0;   // Width dependence of cit
   parameter real WNFACTOR = 0.0;   // Width dependence of nfactor
   parameter real WXJ      = 0.0;   // Width dependence of xj
   parameter real WVSAT    = 0.0;   // Width dependence of vsat
   parameter real WAT      = 0.0;   // Width dependence of at
   parameter real WA0      = 0.0;   // Width dependence of a0 
   parameter real WAGS     = 0.0;   // Width dependence of ags 
   parameter real WA1      = 0.0;   // Width dependence of a1
   parameter real WA2      = 0.0;   // Width dependence of a2
   parameter real WKETA    = 0.0;   // Width dependence of keta
   parameter real WNSUB    = 0.0;   // Width dependence of nsub
   parameter real WNCH     = 0.0;   // Width dependence of nch
   parameter real WNGATE   = 0.0;   // Width dependence of ngate
   parameter real WGAMMA1  = -99.0; // Width dependence of gamma1
   parameter real WGAMMA2  = -99.0; // Width dependence of gamma2
   parameter real WVBX     = -99.0; // Width dependence of vbx
   parameter real WVBM     = 0.0;   // Width dependence of vbm
   parameter real WXT      = 0.0;   // Width dependence of xt
   parameter real WK1      = -99.0; // Width dependence of k1
   parameter real WKT1     = 0.0;   // Width dependence of kt1
   parameter real WKT1L    = 0.0;   // Width dependence of kt1l
   parameter real WKT2     = 0.0;   // Width dependence of kt2
   parameter real WK2      = -99.0; // Width dependence of k2
   parameter real WK3      = 0.0;   // Width dependence of k3
   parameter real WK3B     = 0.0;   // Width dependence of k3b
   parameter real WW0      = 0.0;   // Width dependence of w0
   parameter real WNLX     = 0.0;   // Width dependence of nlx
   parameter real WDVT0    = 0.0;   // Width dependence of dvt0
   parameter real WDVT1    = 0.0;   // Width dependence of dvt1
   parameter real WDVT2    = 0.0;   // Width dependence of dvt2
   parameter real WDVT0W   = 0.0;   // Width dependence of dvt0w
   parameter real WDVT1W   = 0.0;   // Width dependence of dvt1w
   parameter real WDVT2W   = 0.0;   // Width dependence of dvt2w
   parameter real WDROUT   = 0.0;   // Width dependence of drout
   parameter real WDSUB    = 0.0;   // Width dependence of dsub
   parameter real WVTH0    = 0.0;   // Width dependence of vto
   parameter real WVTHO    = 0.0;   // Width dependence of vto
   parameter real WUA      = 0.0;   // Width dependence of ua
   parameter real WUA1     = 0.0;   // Width dependence of ua1
   parameter real WUB      = 0.0;   // Width dependence of ub
   parameter real WUB1     = 0.0;   // Width dependence of ub1
   parameter real WUC      = 0.0;   // Width dependence of uc
   parameter real WUC1     = 0.0;   // Width dependence of uc1
   parameter real WU0      = 0.0;   // Width dependence of u0
   parameter real WUTE     = 0.0;   // Width dependence of ute
   parameter real WVOFF    = 0.0;   // Width dependence of voff
   parameter real WELM     = 0.0;   // Width dependence of elm
   parameter real WDELTA   = 0.0;   // Width dependence of delta
   parameter real WRDSW    = 0.0;   // Width dependence of rdsw 
   parameter real WPRWG    = 0.0;   // Width dependence of prwg 
   parameter real WPRWB    = 0.0;   // Width dependence of prwb 
   parameter real WPRT     = 0.0;   // Width dependence of prt
   parameter real WETA0    = 0.0;   // Width dependence of eta0   
   parameter real WETAB    = 0.0;   // Width dependence of etab   
   parameter real WPCLM    = 0.0;   // Width dependence of pclm   
   parameter real WPDIBLC1 = 0.0;   // Width dependence of pdiblc1   
   parameter real WPDIBLC2 = 0.0;   // Width dependence of pdiblc2   
   parameter real WPDIBLCB = 0.0;   // Width dependence of pdiblcb   
   parameter real WPSCBE1  = 0.0;   // Width dependence of pscbe1   
   parameter real WPSCBE2  = 0.0;   // Width dependence of pscbe2   
   parameter real WPVAG    = 0.0;   // Width dependence of pvag   
   parameter real WWR      = 0.0;   // Width dependence of wr
   parameter real WDWG     = 0.0;   // Width dependence of dwg
   parameter real WDWB     = 0.0;   // Width dependence of dwb
   parameter real WB0      = 0.0;   // Width dependence of b0
   parameter real WB1      = 0.0;   // Width dependence of b1
   parameter real WCGSL    = 0.0;   // Width dependence of cgsl
   parameter real WCGDL    = 0.0;   // Width dependence of cgdl
   parameter real WCKAPPA  = 0.0;   // Width dependence of ckappa
   parameter real WCF      = 0.0;   // Width dependence of cf
   parameter real WCLC     = 0.0;   // Width dependence of clc
   parameter real WCLE     = 0.0;   // Width dependence of cle
   parameter real WALPHA0  = 0.0;   // Width dependence of alpha0
   parameter real WALPHA1  = 0.0;   // Width dependence of alpha1
   parameter real WBETA0   = 0.0;   // Width dependence of beta0
   parameter real WVFBCV   = 0.0;   // Width dependence of vfbcv
   parameter real WVFB     = 0.0;   // Width dependence of vfb
   parameter real WACDE    = 0.0;   // Width dependence of acde
   parameter real WMOIN    = 0.0;   // Width dependence of moin
   parameter real WNOFF    = 0.0;   // Width dependence of noff
   parameter real WVOFFCV  = 0.0;   // Width dependence of voffcv
   
   parameter real PCDSC    = 0.0;   // Cross-term dependence of cdsc
   parameter real PCDSCB   = 0.0;   // Cross-term dependence of cdscb 
   parameter real PCDSCD   = 0.0;   // Cross-term dependence of cdscd
   parameter real PCIT     = 0.0;   // Cross-term dependence of cit
   parameter real PNFACTOR = 0.0;   // Cross-term dependence of nfactor
   parameter real PXJ      = 0.0;   // Cross-term dependence of xj
   parameter real PVSAT    = 0.0;   // Cross-term dependence of vsat
   parameter real PAT      = 0.0;   // Cross-term dependence of at
   parameter real PA0      = 0.0;   // Cross-term dependence of a0 
   parameter real PAGS     = 0.0;   // Cross-term dependence of ags
   parameter real PA1      = 0.0;   // Cross-term dependence of a1
   parameter real PA2      = 0.0;   // Cross-term dependence of a2
   parameter real PKETA    = 0.0;   // Cross-term dependence of keta
   parameter real PNSUB    = 0.0;   // Cross-term dependence of nsub
   parameter real PNCH     = 0.0;   // Cross-term dependence of nch
   parameter real PNGATE   = 0.0;   // Cross-term dependence of ngate
   parameter real PGAMMA1  = -99.0; // Cross-term dependence of gamma1
   parameter real PGAMMA2  = -99.0; // Cross-term dependence of gamma2
   parameter real PVBX     = -99.0; // Cross-term dependence of vbx
   parameter real PVBM     = 0.0;   // Cross-term dependence of vbm
   parameter real PXT      = 0.0;   // Cross-term dependence of xt
   parameter real PK1      = -99.0; // Cross-term dependence of k1
   parameter real PKT1     = 0.0;   // Cross-term dependence of kt1
   parameter real PKT1L    = 0.0;   // Cross-term dependence of kt1l
   parameter real PKT2     = 0.0;   // Cross-term dependence of kt2
   parameter real PK2      = -99.0; // Cross-term dependence of k2
   parameter real PK3      = 0.0;   // Cross-term dependence of k3
   parameter real PK3B     = 0.0;   // Cross-term dependence of k3b
   parameter real PW0      = 0.0;   // Cross-term dependence of w0
   parameter real PNLX     = 0.0;   // Cross-term dependence of nlx
   parameter real PDVT0    = 0.0;   // Cross-term dependence of dvt0
   parameter real PDVT1    = 0.0;   // Cross-term dependence of dvt1
   parameter real PDVT2    = 0.0;   // Cross-term dependence of dvt2
   parameter real PDVT0W   = 0.0;   // Cross-term dependence of dvt0w
   parameter real PDVT1W   = 0.0;   // Cross-term dependence of dvt1w
   parameter real PDVT2W   = 0.0;   // Cross-term dependence of dvt2w
   parameter real PDROUT   = 0.0;   // Cross-term dependence of drout
   parameter real PDSUB    = 0.0;   // Cross-term dependence of dsub
   parameter real PVTH0    = 0.0;   // Cross-term dependence of vto
   parameter real PVTHO    = 0.0;   // Cross-term dependence of vto
   parameter real PUA      = 0.0;   // Cross-term dependence of ua
   parameter real PUA1     = 0.0;   // Cross-term dependence of ua1
   parameter real PUB      = 0.0;   // Cross-term dependence of ub
   parameter real PUB1     = 0.0;   // Cross-term dependence of ub1
   parameter real PUC      = 0.0;   // Cross-term dependence of uc
   parameter real PUC1     = 0.0;   // Cross-term dependence of uc1
   parameter real PU0      = 0.0;   // Cross-term dependence of u0
   parameter real PUTE     = 0.0;   // Cross-term dependence of ute
   parameter real PVOFF    = 0.0;   // Cross-term dependence of voff
   parameter real PELM     = 0.0;   // Cross-term dependence of elm
   parameter real PDELTA   = 0.0;   // Cross-term dependence of delta
   parameter real PRDSW    = 0.0;   // Cross-term dependence of rdsw     
   parameter real PPRWG    = 0.0;   // Cross-term dependence of prwg     
   parameter real PPRWB    = 0.0;   // Cross-term dependence of prwb     
   parameter real PPRT     = 0.0;   // Cross-term dependence of prt 
   parameter real PETA0    = 0.0;   // Cross-term dependence of eta0
   parameter real PETAB    = 0.0;   // Cross-term dependence of etab
   parameter real PPCLM    = 0.0;   // Cross-term dependence of pclm
   parameter real PPDIBLC1 = 0.0;   // Cross-term dependence of pdiblc1
   parameter real PPDIBLC2 = 0.0;   // Cross-term dependence of pdiblc2
   parameter real PPDIBLCB = 0.0;   // Cross-term dependence of pdiblcb
   parameter real PPSCBE1  = 0.0;   // Cross-term dependence of pscbe1
   parameter real PPSCBE2  = 0.0;   // Cross-term dependence of pscbe2
   parameter real PPVAG    = 0.0;   // Cross-term dependence of pvag   
   parameter real PWR      = 0.0;   // Cross-term dependence of wr
   parameter real PDWG     = 0.0;   // Cross-term dependence of dwg
   parameter real PDWB     = 0.0;   // Cross-term dependence of dwb
   parameter real PB0      = 0.0;   // Cross-term dependence of b0
   parameter real PB1      = 0.0;   // Cross-term dependence of b1
   parameter real PCGSL    = 0.0;   // Cross-term dependence of cgsl
   parameter real PCGDL    = 0.0;   // Cross-term dependence of cgdl
   parameter real PCKAPPA  = 0.0;   // Cross-term dependence of ckappa
   parameter real PCF      = 0.0;   // Cross-term dependence of cf
   parameter real PCLC     = 0.0;   // Cross-term dependence of clc
   parameter real PCLE     = 0.0;   // Cross-term dependence of cle
   parameter real PALPHA0  = 0.0;   // Cross-term dependence of alpha0
   parameter real PALPHA1  = 0.0;   // Cross-term dependence of alpha1
   parameter real PBETA0   = 0.0;   // Cross-term dependence of beta0
   parameter real PVFBCV   = 0.0;   // Cross-term dependence of vfbcv
   parameter real PVFB     = 0.0;   // Cross-term dependence of vfb
   parameter real PACDE    = 0.0;   // Cross-term dependence of acde
   parameter real PMOIN    = 0.0;   // Cross-term dependence of moin
   parameter real PNOFF    = 0.0;   // Cross-term dependence of noff
   parameter real PVOFFCV  = 0.0;   // Cross-term dependence of voffcv

   integer   Fatal_Flag;
	  
   real      tox, cox, vth0, nch, ngate, u0, uc, uc1, tnom;
   real      dlcGiven, cgso, cgdo, cgbo;
   real      pdibl1, pdibl2, pdiblb;
   
   real      jctSatCurDensity;
   real      jctSidewallSatCurDensity;
   real      bulkJctPotential;
   real      jctEmissionCoeff;
   real      jctTempExponent;
   real      bulkJctBotGradingCoeff;
   real      sidewallJctPotential;
   real      bulkJctSideGradingCoeff;
   real      GatesidewallJctPotential;
   real      bulkJctGateSideGradingCoeff;
   real      unitAreaJctCap;
   real      unitAreaTempJctCap;
   real      unitLengthSidewallJctCap;
   real      unitLengthSidewallTempJctCap;
   real      unitLengthGateSidewallJctCap;
   real      unitLengthGateSidewallTempJctCap;
   
   
   real      cf;
   real      lnpeak, lngate, lvth0;
   real      lpdibl1, lpdibl2, lpdiblb;
   real      wnpeak, wngate, wvth0;
   real      wpdibl1, wpdibl2, wpdiblb;
   real      pnpeak, pngate, pvth0;
   real      ppdibl1, ppdibl2, ppdiblb;

   real      Temp_, Tnom, TRatio, delTemp;
   real      factor1, Vtm0, Eg0, ni, vtm, Eg;
   real      T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10;
   real      jctTempSatCurDensity, jctSidewallTempSatCurDensity;
   real      PhiB, PhiBSW, PhiBSWG;
   
   real      Ldrn, Wdrn;
   real      tmp, tmp1, tmp2, tmp3, tmp4;
   real      dl, dlc, dw, dwc, leff, weff;
   real      leffCV, weffCV;     
   real      Inv_L, Inv_W, Inv_LW;
   
   real      cdsc, cdscb, cdscd, cit, nfactor, xj;
   real      vsat, at, a0, ags, a1, a2, keta, nsub;
   real      npeak;
   real      gamma1, lgamma1, wgamma1, pgamma1, gamma1Given;
   real      gamma2, lgamma2, wgamma2, pgamma2, gamma2Given;
   real      vbx, lvbx, wvbx, pvbx, vbxGiven;
   real      vfb, vfbGiven;
   real      k1, lk1, wk1, pk1, k1Given;
   real      k2, lk2, wk2, pk2, k2Given;
   
   real      vbm, xt, kt1, kt1l, kt2, k3;
   real      k3b, w0, nlx, dvt0, dvt1, dvt2, dvt0w, dvt1w, dvt2w, drout;
   real      dsub, ua, ua1, ub, ub1, ute;
   real      voff, delta, rdsw, prwg, prwb, prt, eta0, etab, pclm;
   real      pscbe1, pscbe2, pvag, wr, dwg, dwb, b0, b1;
   real      alpha0, alpha1, beta0;
   real      elm, cgsl, cgdl, ckappa, clc, cle, vfbcv, acde, moin, noff_param, Noff2, voffcv_param;

   real      abulkCVfactor, u0temp, vsattemp, rds0, tconst;
   real      phi, sqrtPhi, phis3, Xdep0, litl, vbi, cdep0, ldeb;
   real      nsubGiven, xtGiven, npeakGiven, vth0Given;
   real      vbsc, k1ox, k2ox, theta0vb0, thetaRout, vfbzb;
   real      sheetResistance;
   
   real      sourceArea,        drainArea;
   real      sourceSquares,     drainSquares;
   real      sourcePerimeter,   drainPerimeter;
   real      SourceSatCurrent,  DrainSatCurrent;
   real      sourceConductance, drainConductance;
   
   real      Nvtm, vjsm, vjdm, IsEvjsm, IsEvjdm;

   real      vbs, vgs, vds;
   real      vbd, vgd, vgb;

   // diode equations variables
   real      gbs, cbs, evbs;
   real      gbd, cbd, evbd;

   real      mode, Vds, Vgs, Vbs;
   real      Vbseff, Phis, sqrtPhis, Xdep, V0, lt1, ltw;
   real      Leff, Vtm, Theta0, thetavth, Delt_vth;
   real      temp, TempRatio, dDIBL_Sft_dVd, DIBL_Sft, Vth;

   real      n, Vgs_eff, Vgst, VgstNVt, ExpArg;
   real      Vgsteff, ExpVgst;
   real      Weff, Rds, Abulk0, Abulk, Denomi, ueff;
   real      WVCox, WVCoxRds, Esat, EsatL, Lambda, Vgst2Vtm, Vdsat, vdsat;
   real      Vdseff, Vasat, diffVds, VACLM, VADIBL, Va, VASCBE;
   real      CoxWovL, beta, fgche1, fgche2, gche, Idl;
   real      Idsa, Ids, Isub, cdrain, csub;

   // Depletion capacitance related variables
   real      czbd, czbs;
   real      czbdswg, czbdsw;
   real      czbsswg, czbssw;
   real      mj, mjsw, mjswg;
   real      qbs, capbs;
   real      qbd, capbd;
   real      arg, sarg;
   
   
   // Charge model related variables
   real      qgate, qdrn, qsrc, qbulk;
   real      Qgate, Qdrn, Qsrc, Qbulk;
   
   real      Vfb, CoxWL, Arg1, qinv, Qac0, Vfbeff;
   real      One_Third_CoxWL, Two_Third_CoxWL, AbulkCV;
   real      VdsatCV, Alphaz, VbseffCV, Qsub0, qinoi, VdseffCV;
   real      Cox, V3, Tox, Tox2, LINK, V4, Coxeff, CoxWLcen, Tcen, Ccen;
   real      DeltaPhi;
   
   // extrinsic capacitance related variables
   real      cgso_param, cgdo_param, cgbo_param;
   real      qgso, qgdo;
   real      qgs, qgd, qgb;

   // NQS model
`ifdef NQSMOD
   real sxpart, dxpart;
   real qcheq , qcdump, qdef, gtau_drift, gtau_diff, ScalingFactor, gtau;
   real cqcheq, cqdef;
`endif
   
   real cqgate, cqdrn, cqbulk;
   
   analog
     begin
       	
       	@(initial_step)
          begin
	     tox = TOX;

	     cox = 3.453133e-11 / tox;

`ifdef NQSMOD
	     ScalingFactor = 1.0e-9;
`endif
	     
	     if (NSUB == -99.0)
	       begin
		  nsub      = 6.0e16;
		  nsubGiven = 0;
	       end
	     else
	       begin
		  nsub      = NSUB;
		  nsubGiven = 1;
	       end

	     if (XT == -99.0)
	       begin
		  xt      = 1.55e-7;
		  xtGiven = 0;
	       end
	     else
	       begin
		  xt      = XT;
		  xtGiven = 1;
	       end
	     
	     if (NCH == -99.0)
	       begin
		  npeak      = 1.7e17;
		  npeakGiven = 0;
	       end
	     else
	       begin
		  npeak      = NCH;
		  npeakGiven = 1;
	       end

	     if (GAMMA1 == -99.0)
	       begin
		  gamma1      = 0.0;
		  gamma1Given = 0;
	       end
	     else
	       begin
		  gamma1      = GAMMA1;
		  gamma1Given = 1;
	       end

	     if (LGAMMA1 == -99.0)
	       lgamma1 = 0.0;
	     else
	       lgamma1 = LGAMMA1;

	     if (WGAMMA1 == -99.0)
	       wgamma1      = 0.0;
	     else
	       wgamma1      = WGAMMA1;

	     if (PGAMMA1 == -99.0)
	       pgamma1 = 0.0;
	     else
	       pgamma1 = PGAMMA1;

	     if (GAMMA2 == -99.0)
	       begin
		  gamma2      = 0.0;
		  gamma2Given = 0;
	       end
	     else
	       begin
		  gamma2      = GAMMA2;
		  gamma2Given = 1;
	       end

	     if (LGAMMA2 == -99.0)
	       lgamma2 = 0.0;
	     else
	       lgamma2 = LGAMMA2;

	     if (WGAMMA2 == -99.0)
	       wgamma2      = 0.0;
	     else
	       wgamma2      = WGAMMA2;
	     
	     if (PGAMMA2 == -99.0)
	       pgamma2 = 0.0;
	     else
	       pgamma2 = PGAMMA2;
	     
	     if (VBX == -99.0)
	       begin
		  vbx      = 0.0;
		  vbxGiven = 0.0;
	       end
	     else
	       begin
		  vbx      = VBX;
		  vbxGiven = 1;
	       end

	     if (LVBX == -99.0)
	       lvbx = 0.0;
	     else
	       lvbx = LVBX;

	     if (WVBX == -99.0)
	       wvbx = 0.0;
	     else
	       wvbx = WVBX;

	     if (PVBX == -99.0)
	       pvbx = 0.0;
	     else
	       pvbx = PVBX;
		  
	     if (VFB == -99.0)
	       begin
		  vfb      = 0.0;
		  vfbGiven = 0.0;
	       end
	     else
	       begin
		  vfb      = VFB;
		  vfbGiven = 1.0;
	       end

	     if (K1 == -99.0)
	       begin
		  k1      = 0.0;
		  k1Given = 0.0;
	       end
	     else
	       begin
		  k1      = K1;
		  k1Given = 1.0;
	       end

	     if (LK1 == -99.0)
	       lk1 = 0.0;
	     else
	       lk1 = LK1;

	     if (WK1 == -99.0)
	       wk1 = 0.0;
	     else
	       wk1 = WK1;

	     if (PK1 == -99.0)
	       pk1 = 0.0;
	     else
	       pk1 = PK1;
	     
	     if (K2 == -99.0)
	       begin
		  k2      = 0.0;
		  k2Given = 0.0;
	       end
	     else
	       begin
		  k2      = K2;
		  k2Given = 1.0;
	       end

	     if (LK2 == -99.0)
	       lk2 = 0.0;
	     else
	       lk2 = LK2;

	     if (WK2 == -99.0)
	       wk2 = 0.0;
	     else
	       wk2 = WK2;
	     
	     if (PK2 == -99.0)
	       pk2 = 0.0;
	     else
	       pk2 = PK2;
	     
	     if (NCH > 1.0e20)
	       npeak = npeak * 1.0e-6;
	     else
	       npeak = npeak;
	     
	     if (NGATE > 1.0e23)
	       ngate = NGATE * 1.0e-6;
	     else
	       ngate = NGATE;

	     // calculating ungiven parameters
	     if ( (VTHO == -99.0) && (VTH0 == -99.0) )
	       begin
		  vth0      = (TYPE > 0.0) ? 0.7 : -0.7;
		  vth0Given = 0;
	       end
	     else
	       begin
		  vth0Given = 1;
		  if (VTH0 == -99)
		    vth0 = VTHO;
		  if (VTHO == -99)
		    vth0 = VTH0;
	       end

	     if (UC  == -99.0 )
	       uc  = (MOBMOD == 3) ? -0.0465 : -0.0465e-9;
	     else
	       uc  = UC;
	     
	     if (UC1 == -99.0)
	       uc1 = (MOBMOD == 3) ? -0.056  : -0.056e-9;
	     else
	       uc1 = UC1;

	     if (U0 == -99.0)
	       u0 = (TYPE == 1) ? 0.067 : 0.025;
	     else
	       u0 = U0;

	     if (TNOM == -99.0)
	       tnom = 300.15;
	     else
	       tnom = TNOM + 273.15;

	     if (DLC == -99.0)
	       begin
		  dlcGiven = 0;
		  dlc      = LINT;
	       end
	     else
	       begin
		  dlcGiven = 1;
		  dlc      = DLC;
	       end
	     
	     if (CGSO == -99.0)
	       if ( dlcGiven && (dlc > 0.0) )
	    	 cgso_param = dlc * cox - CGSL;
	       else
		 cgso_param = 0.6 * XJ * cox;
	     else
	       cgso_param = CGSO;
	     
	     if (CGDO == -99.0)
	       if ( dlcGiven && (dlc > 0.0) )
	    	 cgdo_param = dlc * cox - CGDL;
	       else
		 cgdo_param = 0.6 * XJ * cox;
	     else
	       cgdo_param = CGDO;
	     
	     if (CGBO == -99.0)
	       cgbo_param = 2.0 * DWC * cox;
	     else
	       cgbo_param = CGBO;
	     
	     pdibl1 = PDIBLC1;
	     pdibl2 = PDIBLC2;
	     pdiblb = PDIBLCB;

	     jctSatCurDensity             = JS;
	     jctSidewallSatCurDensity     = JSW;
	     bulkJctPotential             = PB;
	     jctEmissionCoeff             = NJ;
	     jctTempExponent              = XTI;
	     bulkJctBotGradingCoeff       = MJ;
	     sidewallJctPotential         = PBSW;
	     bulkJctSideGradingCoeff      = MJSW;
	     GatesidewallJctPotential     = PBSWG;
	     bulkJctGateSideGradingCoeff  = MJSWG;
	     unitAreaJctCap               = CJ;
	     unitLengthSidewallJctCap     = CJSW;
	     unitLengthGateSidewallJctCap = CJSWG;
	     
	     sheetResistance              = RSH;
	     drainSquares                 = NRD;
	     sourceSquares                = NRS;
	     drainArea                    = AD;
	     sourceArea                   = AS;
	     drainPerimeter               = PD;
	     sourcePerimeter              = PS;
	     
	     if (CF == -99.0)
	       cf = 2.0 * `EPSOX / `PI * ln(1.0 + 0.4e-6 / tox);
	     else
	       cf = CF;

	     // Channel length dependance parameters
	     if (LNCH > 1.0e20)
	       lnpeak = LNCH * 1.0e-6;
	     else
	       lnpeak = LNCH;

	     if (LNGATE > 1.0e23)
	       lngate = LNGATE * 1.0e-6;
	     else
	       lngate = LNGATE;

	     if (LVTH0 != 0.0)
	       lvth0 = LVTH0;
	     if (LVTHO != 0.0)
	       lvth0 = LVTHO;
	     
	     lpdibl1 = LPDIBLC1;
	     lpdibl2 = LPDIBLC2;
	     lpdiblb = LPDIBLCB;

	     // Channel width dependance parameters
	     if (WNCH > 1.0e20)
	       wnpeak = WNCH * 1.0e-6;
	     else
	       wnpeak = WNCH;

	     if (WNGATE > 1.0e23)
	       wngate = WNGATE * 1.0e-6;
	     else
	       wngate = WNGATE;

	     if (WVTH0 != 0.0)
	       wvth0 = WVTH0;
	     if (WVTHO != 0.0)
	       wvth0 = WVTHO;
	     
	     wpdibl1 = WPDIBLC1;
	     wpdibl2 = WPDIBLC2;
	     wpdiblb = WPDIBLCB;
	     
	     // Cross-term dependence parameters
	     if (PNCH > 1.0e20)
	       pnpeak = PNCH * 1.0e-6;
	     else
	       pnpeak = PNCH;

	     if (PNGATE > 1.0e20)
	       pngate = PNGATE * 1.0e-6;
	     else
	       pngate = PNGATE;
	     
	     if (PVTH0 != 0.0)
	       pvth0 = PVTH0;
	     if (PVTHO != 0.0)
	       pvth0 = PVTHO;

	     ppdibl1 = PPDIBLC1;
	     ppdibl2 = PPDIBLC2;
	     ppdiblb = PPDIBLCB;
	     
             if (bulkJctPotential < 0.1)
	       begin
		  bulkJctPotential = 0.1;
		  $strobe ("Given pb is less than 0.1. Pb is set to 0.1.");
	       end
	     if (sidewallJctPotential < 0.1)
	       begin
		  sidewallJctPotential = 0.1;
		  $strobe ("Given pbsw is less than 0.1. Pbsw is set to 0.1.");
	       end
	     if (GatesidewallJctPotential < 0.1)
	       begin
		  GatesidewallJctPotential = 0.1;
		  $strobe ("Given pbswg is less than 0.1. Pbswg is set to 0.1.");
	       end

	     Temp_   = $temperature;
	     Tnom   = tnom;
	     TRatio = Temp_ / Tnom;

	     factor1 = sqrt(`EPSSI / `EPSOX * tox);

	     Vtm0    = `KboQ * Tnom;
	     Eg0     = 1.16 - 7.02e-4 * Tnom * Tnom / (Tnom + 1108.0);
             ni      = 1.45e10 * (Tnom / 300.15) * sqrt(Tnom / 300.15) * exp(21.5565981 - Eg0 / (2.0 * Vtm0));

             vtm     = `KboQ * Temp_;
             Eg      = 1.16 - 7.02e-4 * Temp_ * Temp_ / (Temp_ + 1108.0);

	     if (Temp_ != Tnom)
	       begin
		  T0 = Eg0 / Vtm0 - Eg / vtm + jctTempExponent * ln(Temp_ / Tnom);
		  T1 = exp(T0 / jctEmissionCoeff);
		  
		  jctTempSatCurDensity         = jctSatCurDensity         * T1;
		  jctSidewallTempSatCurDensity = jctSidewallSatCurDensity * T1;
	       end
	     else
	       begin
		  jctTempSatCurDensity         = jctSatCurDensity;
		  jctSidewallTempSatCurDensity = jctSidewallSatCurDensity;
	       end

	     if (jctTempSatCurDensity < 0.0)
	       jctTempSatCurDensity = 0.0;
	     if (jctSidewallTempSatCurDensity < 0.0)
	       jctSidewallTempSatCurDensity = 0.0;

	     /* Temperature dependence of D/B and S/B diode capacitance */
	     delTemp = $temperature - tnom;
	     T0      = TCJ * delTemp;

	     if (T0 >= -1.0)
	       begin
		  if (VERSION == 3.24)
		    unitAreaTempJctCap = unitAreaJctCap * (1.0 + T0);
		  else
		    unitAreaJctCap     = unitAreaJctCap * (1.0 + T0);
	       end
	     else if (unitAreaJctCap > 0.0)
	       begin
		  if (VERSION == 3.24)
		    unitAreaTempJctCap = 0.0;
		  else
		    unitAreaJctCap = 0.0;
		  $strobe ("Temperature effect has caused cj to be negative. Cj is clamped to zero.");
	       end

	     T0 = TCJSW * delTemp;
	     if (T0 >= -1.0)
	       begin
		  if (VERSION == 3.24)
		    unitLengthSidewallTempJctCap = unitLengthSidewallJctCap * (1.0 + T0);
		  else
		    unitLengthSidewallJctCap     = unitLengthSidewallJctCap * (1.0 + T0);
	       end
	     else if (unitLengthSidewallJctCap > 0.0)
	       begin
		  if (VERSION == 3.24)
		    unitLengthSidewallTempJctCap = 0.0;
		  else
		    unitLengthSidewallJctCap = 0.0;
		  $strobe ("Temperature effect has caused cjsw to be negative. Cjsw is clamped to zero.");
	       end

	     T0 = TCJSWG * delTemp;
	     if (T0 >= -1.0)
	       begin
		  if (VERSION == 3.24)
		    unitLengthGateSidewallTempJctCap = unitLengthGateSidewallJctCap * (1.0 + T0);
		  else
		    unitLengthGateSidewallJctCap     = unitLengthGateSidewallJctCap * (1.0 + T0);
	       end
	     else if (unitLengthGateSidewallJctCap > 0.0)
	       begin
		  if (VERSION == 3.24)
		    unitLengthGateSidewallTempJctCap = 0.0;
		  else
		    unitLengthGateSidewallJctCap = 0.0;
		  $strobe ("Temperature effect has caused cjswg to be negative. Cjswg is clamped to zero.");
	       end

	     PhiB = bulkJctPotential - TPB * delTemp;
             if (PhiB < 0.01)
	       begin
		  PhiB = 0.01;
		  $strobe ("Temperature effect has caused pb to be less than 0.01. Pb is clamped to 0.01.");
	       end
	     
             PhiBSW = sidewallJctPotential - TPBSW * delTemp;
             if (PhiBSW <= 0.01)
	       begin
		  PhiBSW = 0.01;
		  $strobe ("Temperature effect has caused pbsw to be less than 0.01. Pbsw is clamped to 0.01.");
	       end
	     
	     PhiBSWG = GatesidewallJctPotential - TPBSWG * delTemp;
             if (PhiBSWG <= 0.01)
	       begin
		  PhiBSWG = 0.01;
		  $strobe ("Temperature effect has caused pbswg to be less than 0.01. Pbswg is clamped to 0.01.");
	       end

	     /* End of junction capacitance */

	     /*** Effective channel length and width calculation ***/
	     Ldrn = L;
	     Wdrn = W;
	     
	     T0 = pow(Ldrn, LLN);
             T1 = pow(Wdrn, LWN);

	     tmp1 = LL / T0  + LW / T1  + LWL / (T0 * T1);
             dl   = LINT + tmp1;
	     
             tmp2 = LLC / T0 + LWC / T1 + LWLC / (T0 * T1);
             dlc  = dlc  + tmp2;

	     T2 = pow(Ldrn, WLN);
             T3 = pow(Wdrn, WWN);
	     
	     tmp1 = WL / T2 + WW / T3 + WWL / (T2 * T3);
             dw   = WINT + tmp1;
             tmp2 = WLC / T2 + WWC / T3 + WWLC / (T2 * T3);
             dwc  = DWC  + tmp2;

	     leff = L - 2.0 * dl;
	     
             if (leff <= 0.0)
	       begin
		  $strobe ("BSIM3: device %m: Effective channel length <= 0");
		  // $finish(1);
	       end

	     weff = W - 2.0 * dw;
	     if (leff <= 0.0)
	       begin
		  $strobe ("BSIM3: device %m: Effective channel width <= 0");
		  // $finish(1);
	       end

	     leffCV = L - 2.0 * dlc;
             if (leffCV <= 0.0)
	       begin
		  $strobe ("BSIM3: device %m: Effective channel length for C-V <= 0");
		  // $finish(1);
               end

	     weffCV = W - 2.0 * dwc;
             if (weffCV <= 0.0)
	       begin
		  $strobe ("BSIM3: device %m: Effective channel width for C-V <= 0");
		  // $finish(1);
	       end

	     if (BINUNIT == 1)
	       begin
		  Inv_L  = 1.0e-6  / leff;
		  Inv_W  = 1.0e-6  / weff;
		  Inv_LW = 1.0e-12 / (leff * weff);
	       end
	     else
	       begin
		  Inv_L  = 1.0 / leff;
		  Inv_W  = 1.0 / weff;
		  Inv_LW = 1.0 / (leff * weff);
	       end

	     cdsc  = CDSC  + LCDSC  * Inv_L + WCDSC  * Inv_W + PCDSC  * Inv_LW;
	     cdscb = CDSCB + LCDSCB * Inv_L + WCDSCB * Inv_W + PCDSCB * Inv_LW; 
    	     cdscd = CDSCD + LCDSCD * Inv_L + WCDSCD * Inv_W + PCDSCD * Inv_LW; 
	     
	     cit     = CIT     + LCIT     * Inv_L + WCIT     * Inv_W + PCIT     * Inv_LW;
	     nfactor = NFACTOR + LNFACTOR * Inv_L + WNFACTOR * Inv_W + PNFACTOR * Inv_LW;
	     
	     xj   = XJ   + LXJ   * Inv_L + WXJ   * Inv_W + PXJ   * Inv_LW;
	     vsat = VSAT + LVSAT * Inv_L + WVSAT * Inv_W + PVSAT * Inv_LW;
	     
	     at = AT + LAT * Inv_L + WAT * Inv_W + PAT * Inv_LW;
	     a0 = A0 + LA0 * Inv_L + WA0 * Inv_W + PA0 * Inv_LW; 
	     
	     ags  = AGS  + LAGS  * Inv_L + WAGS  * Inv_W + PAGS  * Inv_LW;
	     a1   = A1   + LA1   * Inv_L + WA1   * Inv_W + PA1   * Inv_LW;
	     a2   = A2   + LA2   * Inv_L + WA2   * Inv_W + PA2   * Inv_LW;
	     keta = KETA + LKETA * Inv_L + WKETA * Inv_W + PKETA * Inv_LW;
	     nsub = nsub + LNSUB * Inv_L + WNSUB * Inv_W + PNSUB * Inv_LW;
	     
	     npeak = npeak + lnpeak * Inv_L + wnpeak * Inv_W + pnpeak * Inv_LW;
	     ngate = ngate + lngate * Inv_L + wngate * Inv_W + pngate * Inv_LW;
	     
	     gamma1 = gamma1 + lgamma1 * Inv_L + wgamma1 * Inv_W + pgamma1 * Inv_LW;
	     gamma2 = gamma2 + lgamma2 * Inv_L + wgamma2 * Inv_W + pgamma2 * Inv_LW;
	     
	     vbx = vbx + lvbx * Inv_L + wvbx * Inv_W + pvbx * Inv_LW;
	     vbm = VBM + LVBM * Inv_L + WVBM * Inv_W + PVBM * Inv_LW;
	     xt  = xt  + LXT  * Inv_L + WXT  * Inv_W + PXT  * Inv_LW;
             vfb = vfb + LVFB * Inv_L + WVFB * Inv_W + PVFB * Inv_LW;
	     
	     k1   = k1   + lk1   * Inv_L + wk1   * Inv_W + pk1   * Inv_LW;
	     kt1  = KT1  + LKT1  * Inv_L + WKT1  * Inv_W + PKT1  * Inv_LW;
	     kt1l = KT1L + LKT1L * Inv_L + WKT1L * Inv_W + PKT1L * Inv_LW;
	     k2   = k2   + lk2   * Inv_L + wk2   * Inv_W + pk2   * Inv_LW;
	     kt2  = KT2  + LKT2  * Inv_L + WKT2  * Inv_W + PKT2  * Inv_LW;
	     k3   = K3   + LK3   * Inv_L + WK3   * Inv_W + PK3   * Inv_LW;
	     k3b  = K3B  + LK3B  * Inv_L + WK3B  * Inv_W + PK3B  * Inv_LW;
	     
	     w0  = W0  + LW0  * Inv_L + WW0  * Inv_W + PW0  * Inv_LW;
	     nlx = NLX + LNLX * Inv_L + WNLX * Inv_W + PNLX * Inv_LW;
	     
	     dvt0 = DVT0 + LDVT0 * Inv_L + WDVT0 * Inv_W + PDVT0 * Inv_LW;
	     dvt1 = DVT1 + LDVT1 * Inv_L + WDVT1 * Inv_W + PDVT1 * Inv_LW;
	     dvt2 = DVT2 + LDVT2 * Inv_L + WDVT2 * Inv_W + PDVT2 * Inv_LW;
	     
	     dvt0w = DVT0W + LDVT0W * Inv_L + WDVT0W * Inv_W + PDVT0W * Inv_LW;
	     dvt1w = DVT1W + LDVT1W * Inv_L + WDVT1W * Inv_W + PDVT1W * Inv_LW;
	     dvt2w = DVT2W + LDVT2W * Inv_L + WDVT2W * Inv_W + PDVT2W * Inv_LW;
	     drout = DROUT + LDROUT * Inv_L + WDROUT * Inv_W + PDROUT * Inv_LW;
	     
	     dsub = DSUB + LDSUB * Inv_L + WDSUB * Inv_W + PDSUB * Inv_LW;
	     vth0 = vth0 + lvth0 * Inv_L + wvth0 * Inv_W + pvth0 * Inv_LW;
	     
	     ua  = UA  + LUA  * Inv_L + WUA  * Inv_W + PUA  * Inv_LW;
	     ua1 = UA1 + LUA1 * Inv_L + WUA1 * Inv_W + PUA1 * Inv_LW;
	     ub  = UB  + LUB  * Inv_L + WUB  * Inv_W + PUB  * Inv_LW;
	     ub1 = UB1 + LUB1 * Inv_L + WUB1 * Inv_W + PUB1 * Inv_LW;
	     uc  = uc  + LUC  * Inv_L + WUC  * Inv_W + PUC  * Inv_LW;
	     uc1 = uc1 + LUC1 * Inv_L + WUC1 * Inv_W + PUC1 * Inv_LW;
	     u0  = u0  + LU0  * Inv_L + WU0  * Inv_W + PU0  * Inv_LW;
	     ute = UTE + LUTE * Inv_L + WUTE * Inv_W + PUTE * Inv_LW;
	     
	     voff  = VOFF  + LVOFF  * Inv_L + WVOFF  * Inv_W + PVOFF  * Inv_LW;
	     delta = DELTA + LDELTA * Inv_L + WDELTA * Inv_W + PDELTA * Inv_LW;
	     
	     rdsw = RDSW + LRDSW * Inv_L + WRDSW * Inv_W + PRDSW * Inv_LW;
	     prwg = PRWG + LPRWG * Inv_L + WPRWG * Inv_W + PPRWG * Inv_LW;
	     prwb = PRWB + LPRWB * Inv_L + WPRWB * Inv_W + PPRWB * Inv_LW;
	     prt  = PRT  + LPRT  * Inv_L + WPRT  * Inv_W + PPRT  * Inv_LW;
	     eta0 = ETA0 + LETA0 * Inv_L + WETA0 * Inv_W + PETA0 * Inv_LW;
	     etab = ETAB + LETAB * Inv_L + WETAB * Inv_W + PETAB * Inv_LW;
	     pclm = PCLM + LPCLM * Inv_L + WPCLM * Inv_W + PPCLM * Inv_LW;
	     
	     pdibl1 = pdibl1 + lpdibl1 * Inv_L + wpdibl1 * Inv_W + ppdibl1 * Inv_LW;
	     pdibl2 = pdibl2 + lpdibl2 * Inv_L + wpdibl2 * Inv_W + ppdibl2 * Inv_LW;
	     pdiblb = pdiblb + lpdiblb * Inv_L + wpdiblb * Inv_W + ppdiblb * Inv_LW;
	     pscbe1 = PSCBE1 + LPSCBE1 * Inv_L + WPSCBE1 * Inv_W + PPSCBE1 * Inv_LW;
	     pscbe2 = PSCBE2 + LPSCBE2 * Inv_L + WPSCBE2 * Inv_W + PPSCBE2 * Inv_LW;
	     pvag   = PVAG   + LPVAG   * Inv_L + WPVAG   * Inv_W + PPVAG   * Inv_LW;
	     
	     wr  = WR  + LWR  * Inv_L + WWR  * Inv_W + PWR  * Inv_LW;
	     dwg = DWG + LDWG * Inv_L + WDWG * Inv_W + PDWG * Inv_LW;
	     dwb = DWB + LDWB * Inv_L + WDWB * Inv_W + PDWB * Inv_LW;
	     b0  = B0  + LB0  * Inv_L + WB0  * Inv_W + PB0  * Inv_LW;
	     b1  = B1  + LB1  * Inv_L + WB1  * Inv_W + PB1  * Inv_LW;
	     
	     alpha0 = ALPHA0 + LALPHA0 * Inv_L + WALPHA0 * Inv_W + PALPHA0 * Inv_LW;
             alpha1 = ALPHA1 + LALPHA1 * Inv_L + WALPHA1 * Inv_W + PALPHA1 * Inv_LW;
	     beta0  = BETA0  + LBETA0  * Inv_L + WBETA0  * Inv_W + PBETA0  * Inv_LW;
	     
	     /* CV model */
	     elm  = ELM  + LELM  * Inv_L + WELM  * Inv_W + PELM  * Inv_LW;
	     cgsl = CGSL + LCGSL * Inv_L + WCGSL * Inv_W + PCGSL * Inv_LW;
	     cgdl = CGDL + LCGDL * Inv_L + WCGDL * Inv_W + PCGDL * Inv_LW;
	     
	     ckappa = CKAPPA + LCKAPPA * Inv_L + WCKAPPA * Inv_W + PCKAPPA * Inv_LW;
	     
	     cf  = cf  + LCF  * Inv_L + WCF  * Inv_W + PCF  * Inv_LW;
	     clc = CLC + LCLC * Inv_L + WCLC * Inv_W + PCLC * Inv_LW;
	     cle = CLE + LCLE * Inv_L + WCLE * Inv_W + PCLE * Inv_LW;
	     
	     vfbcv  = VFBCV  + LVFBCV  * Inv_L + WVFBCV  * Inv_W + PVFBCV  * Inv_LW;
             acde   = ACDE   + LACDE   * Inv_L + WACDE   * Inv_W + PACDE   * Inv_LW;
             moin   = MOIN   + LMOIN   * Inv_L + WMOIN   * Inv_W + PMOIN   * Inv_LW;
             noff_param   = NOFF   + LNOFF   * Inv_L + WNOFF   * Inv_W + PNOFF   * Inv_LW;
             voffcv_param = VOFFCV + LVOFFCV * Inv_L + WVOFFCV * Inv_W + PVOFFCV * Inv_LW;
	     
	     abulkCVfactor = 1.0 + pow((clc / leffCV), cle);

	     T0 = (TRatio - 1.0);

	     ua = ua + ua1 * T0;
	     ub = ub + ub1 * T0;
	     uc = uc + uc1 * T0;

	     if (u0 > 1.0) 
               u0 = u0 / 1.0e4;
	     
             u0temp   = u0 * pow(TRatio, ute); 
             vsattemp = vsat - at * T0;
	     rds0     = (rdsw + prt * T0) / pow(weff * 1E6, wr);

	     // *** check model and instance parameters ***
	     Fatal_Flag = 0;
	     
	     if (nlx < -leff)
	       begin
	          $strobe ("Fatal: Nlx = %g is less than -Leff.", nlx);
		  Fatal_Flag = 1;
               end

	     if (tox <= 0.0)
	       begin
 		  $strobe ("Fatal: Tox = %g is not positive.", tox);
		  Fatal_Flag = 1;
	       end

             if (TOXM <= 0.0)
               begin
		  $strobe ("Fatal: Toxm = %g is not positive.", TOXM);
		  Fatal_Flag = 1;
               end

	     if (npeak <= 0.0)
	       begin
		  $strobe ("Fatal: Nch = %g is not positive.", npeak);
		  Fatal_Flag = 1;
	       end
	     if (nsub <= 0.0)
	       begin
		  $strobe ("Fatal: Nsub = %g is not positive.", nsub);
		  Fatal_Flag = 1;
	       end
	     if (ngate < 0.0)
	       begin
		  $strobe ("Fatal: Ngate = %g Ngate is not positive.", ngate);
		  Fatal_Flag = 1;
	       end
	     if (ngate > 1.0e25)
	       begin
		  $strobe ("Fatal: Ngate = %g Ngate is too high", ngate);
		  Fatal_Flag = 1;
	       end
	     if (xj <= 0.0)
	       begin
		  $strobe ("Fatal: Xj = %g is not positive.", xj);
		  Fatal_Flag = 1;
	       end

	     if (dvt1 < 0.0)
	       begin
		  $strobe ("Fatal: Dvt1 = %g is negative.", dvt1);   
		  Fatal_Flag = 1;
	       end
	     
	     if (dvt1w < 0.0)
	       begin
		  $strobe ("Fatal: Dvt1w = %g is negative.", dvt1w);
		  Fatal_Flag = 1;
	       end
	     
	     if (w0 == -weff)
	       begin
		  $strobe ("Fatal: (W0 + Weff) = 0 causing divided-by-zero.");
		  Fatal_Flag = 1;
               end   

	     if (dsub < 0.0)
	       begin
		  $strobe ("Fatal: Dsub = %g is negative.", dsub);
		  Fatal_Flag = 1;
	       end
	     if (b1 == -weff)
	       begin
		  $strobe ("Fatal: (B1 + Weff) = 0 causing divided-by-zero.");
		  Fatal_Flag = 1;
               end  
             if (u0temp <= 0.0)
	       begin
		  $strobe ("Fatal: u0 at current temperature = %g is not positive.", u0temp);
		  Fatal_Flag = 1;
               end
	     
	     /* Check delta parameter */      
             if (delta < 0.0)
	       begin
		  $strobe ("Fatal: Delta = %g is less than zero.", delta);
		  Fatal_Flag = 1;
               end      

	     if (vsattemp <= 0.0)
	       begin
		  $strobe ("Fatal: Vsat at current temperature = %g is not positive.", vsattemp);
		  Fatal_Flag = 1;
	       end
	     /* Check Rout parameters */
	     if (pclm <= 0.0)
	       begin
		  $strobe ("Fatal: Pclm = %g is not positive.", pclm);
		  Fatal_Flag = 1;
	       end

	     if (drout < 0.0)
	       begin
		  $strobe ("Fatal: Drout = %g is negative.", drout);
		  Fatal_Flag = 1;
	       end

             if (pscbe2 <= 0.0)
               begin
		  $strobe ("Warning: Pscbe2 = %g is not positive.", pscbe2);
               end

	     if (unitLengthSidewallJctCap > 0.0 || unitLengthGateSidewallJctCap > 0.0)
	       begin
		  if (drainPerimeter < weff)
		    begin
		       $strobe ("Warning: Pd = %g is less than W.", drainPerimeter);
		    end
		  if (sourcePerimeter < weff)
		    begin
		       $strobe ("Warning: Ps = %g is less than W.", sourcePerimeter);
		    end
	       end

             if (noff_param < 0.1)
      	       $strobe ("Warning: Noff = %g is too small.", noff_param);
             if (noff_param > 4.0)
               $strobe ("Warning: Noff = %g is too large.", noff_param);
             if (voffcv_param < -0.5)
               $strobe ("Warning: Voffcv = %g is too small.", voffcv_param);
             if (voffcv_param > 0.5)
               $strobe ("Warning: Voffcv = %g is too large.", voffcv_param);
             
             if (IJTH < 0.0)
               begin
		  $strobe ("Fatal: Ijth = %g cannot be negative.", IJTH);
		  Fatal_Flag = 1;
               end

	     /* Check capacitance parameters */
             if (clc < 0.0)
	       begin
		  $strobe ("Fatal: Clc = %g is negative.", clc);
		  Fatal_Flag = 1;
               end      

             if (moin < 5.0)
               $strobe ("Warning: Moin = %g is too small.", moin);
             if (moin > 25.0)
               $strobe ("Warning: Moin = %g is too large.", moin);
             if (((acde < 0.4) && !(VERSION == 3.24)) ||
		 ((acde < 0.4) &&  (VERSION == 3.24) && (CAPMOD == 3.0)))
               $strobe ("Warning: Acde = %g is too small.", acde);
             if (((acde > 1.6) && !(VERSION == 3.24)) ||
		 ((acde > 1.6) &&  (VERSION == 3.24) && (CAPMOD == 3.0)))
               $strobe ("Warning: Acde = %g is too large.", acde);
             // *** end of parameters checking ***

	     if (PARAMCHK ==1)
	       begin
		  /* Check L and W parameters */ 
		  if (leff <= 5.0e-8)
		    $strobe ("Warning: Leff = %g may be too small.", leff);
		  if (leffCV <= 5.0e-8)
		    $strobe ("Warning: Leff for CV = %g may be too small.", leffCV);
		  if (weff <= 1.0e-7)
		    $strobe ("Warning: Weff = %g may be too small.", weff);
		  if (weffCV <= 1.0e-7)
		    $strobe ("Warning: Weff for CV = %g may be too small.", weffCV);
		  
		  /* Check threshold voltage parameters */
		  if (nlx < 0.0)
		    $strobe ("Warning: Nlx = %g is negative.", nlx);
		  if (tox < 1.0e-9)
		    $strobe ("Warning: Tox = %g is less than 10A.", tox);

		  if (npeak <= 1.0e15)
		    $strobe ("Warning: Nch = %g may be too small.", npeak);
		  else if (npeak >= 1.0e21)
		    $strobe ("Warning: Nch = %g may be too large.", npeak);

		  if (nsub <= 1.0e14)
		    $strobe ("Warning: Nsub = %g may be too small.", nsub);
		  else if (nsub >= 1.0e21)
		    $strobe ("Warning: Nsub = %g may be too large.", nsub);
		  
		  if ((ngate > 0.0) && (ngate <= 1.0e18))
		    $strobe ("Warning: Ngate = %g is less than 1.E18cm^-3.", ngate);
		  
		  if (dvt0 < 0.0)
		    $strobe ("Warning: Dvt0 = %g is negative.", dvt0);   
		  if (abs(1.0e-6 / (w0 + weff)) > 10.0)
		    $strobe ("Warning: (W0 + Weff) may be too small.");
		  
		  /* Check subthreshold parameters */
		  if (nfactor < 0.0)
		    $strobe ("Warning: Nfactor = %g is negative.", nfactor);
		  if (cdsc < 0.0)
		    $strobe ("Warning: Cdsc = %g is negative.", cdsc);
		  if (cdscd < 0.0)
		    $strobe ("Warning: Cdscd = %g is negative.", cdscd);
		  /* Check DIBL parameters */
		  if (eta0 < 0.0)
		    $strobe ("Warning: Eta0 = %g is negative.", eta0); 		  
		  /* Check Abulk parameters */	    
		  if (abs(1.0e-6 / (b1 + weff)) > 10.0)
       		    $strobe ("Warning: (B1 + Weff) may be too small.");

		  /* Check Saturation parameters */
     		  if (a2 < 0.01)
		    begin
		       $strobe ("Warning: A2 = %g is too small. Set to 0.01.", a2);
		       a2 = 0.01;
		    end
		  else if (a2 > 1.0)
		    begin
		       $strobe ("Warning: A2 = %g is larger than 1. A2 is set to 1 and A1 is set to 0.", a2);
		       a2 = 1.0;
		       a1 = 0.0;
		    end

		  if (rdsw < 0.0)
		    begin
		       $strobe ("Warning: Rdsw = %g is negative. Set to zero.", rdsw);
		       rdsw = 0.0;
		       rds0 = 0.0;
		    end
		  else if ((rds0 > 0.0) && (rds0 < 0.001))
		    begin
		       $strobe ("Warning: Rds at current temperature = %g is less than 0.001 ohm. Set to zero.", rds0);
		       rds0 = 0.0;
		    end
		  if (vsattemp < 1.0e3)
		    $strobe ("Warning: Vsat at current temperature = %g may be too small.", vsattemp);
		  if (pdibl1 < 0.0)
		    $strobe ("Warning: Pdibl1 = %g is negative.", pdibl1);
		  if (pdibl2 < 0.0)
		    $strobe ("Warning: Pdibl2 = %g is negative.", pdibl2);

		  /* Check overlap capacitance parameters */
		  if (cgdo_param < 0.0)
		    begin
		       $strobe ("Warning: cgdo = %g is negative. Set to zero.", cgdo_param);
		       cgdo_param = 0.0;
		    end      
		  if (cgso_param < 0.0)
		    begin
		       $strobe ("Warning: cgso = %g is negative. Set to zero.", cgso_param);
		       cgso_param = 0.0;
		    end      
		  if (cgbo_param < 0.0)
		    begin
		       $strobe ("Warning: cgbo = %g is negative. Set to zero.", cgbo_param);
		       cgbo_param = 0.0;
		    end

	       end/* loop for the parameter check for warning messages */      
	     
	     if (Fatal_Flag)
	       // $finish(1);
	     	     
	     cgdo_param = (cgdo_param + cf) * weffCV;
             cgso_param = (cgso_param + cf) * weffCV;
             cgbo_param =  cgbo_param       * leffCV;

	     T0     = leffCV * leffCV;
             tconst = u0temp * elm / (cox * weffCV * leffCV * T0);
	     
             if ( !npeakGiven && gamma1Given )
               begin
		  T0 = gamma1 * cox;
                  npeak = 3.021E22 * T0 * T0;
	       end

	     phi     = 2.0 * Vtm0 * ln(npeak / ni);
	     sqrtPhi = sqrt(phi);
	     phis3   = sqrtPhi * phi;

             Xdep0     = sqrt(2.0 * `EPSSI / (`Charge_q * npeak * 1.0e6)) * sqrtPhi;
	     
             litl = sqrt(3.0 * xj * tox);
             vbi  = Vtm0 * ln(1.0e20 * npeak / (ni * ni));
	     
             cdep0 = sqrt(`Charge_q * `EPSSI * npeak * 1.0e6 / 2.0 / phi);

             ldeb = sqrt(`EPSSI * Vtm0 / (`Charge_q * npeak * 1.0e6)) / 3.0;
             acde = acde * pow((npeak / 2.0e16), -0.25);
	     
	     if ( k1Given || k2Given )
	       begin
		  if (!k1Given)
	            begin
		       $strobe ("Warning: k1 should be specified with k2.");
                       k1 = 0.53;
                    end
                  if (!k2Given)
	            begin
		       $strobe ("Warning: k2 should be specified with k1.");
                       k2 = -0.0186;
                    end
		  
                  if (nsubGiven)
                    $strobe ("Warning: nsub is ignored because k1 or k2 is given.");
                  if (xtGiven)
                    $strobe ("Warning: xt is ignored because k1 or k2 is given.");
                  if (vbxGiven)
                    $strobe ("Warning: vbx is ignored because k1 or k2 is given.");
                  if (gamma1Given)
                    $strobe ("Warning: gamma1 is ignored because k1 or k2 is given.");
                  if (gamma2Given)
                    $strobe ("Warning: gamma2 is ignored because k1 or k2 is given.");
               end
             else
	       begin
		  if (!vbxGiven)
                    vbx = phi - 7.7348e-4 * npeak * xt * xt;
	          if (vbx > 0.0)
		    vbx = -vbx;
	          if (vbm > 0.0)
                    vbm = -vbm;
		  
                  if (!gamma1Given)
                    gamma1 = 5.753e-12 * sqrt(npeak) / cox;
                  if (!gamma2Given)
                    gamma2 = 5.753e-12 * sqrt(nsub)  / cox;

                  T0 = gamma1 - gamma2;
                  T1 = sqrt(phi - vbx) - sqrtPhi;
                  T2 = sqrt(phi * (phi - vbm)) - phi;
		  
                  k2 = T0 * T1 / (2.0 * T2 + vbm);
                  k1 = gamma2 - 2.0 * k2 * sqrt(phi - vbm);

	       end
	     
	     if (k2 < 0.0)
	       begin
		  T0   = 0.5 * k1 / k2;
                  vbsc = 0.9 * (phi - T0 * T0);
		  
		  if (vbsc > -3.0)
		    vbsc = -3.0;
		  else if (vbsc < -30.0)
		    vbsc = -30.0;
	       end
	     else
	       vbsc = -30.0;

	     if (vbsc > vbm)
	       vbsc = vbm;

             if (!vfbGiven)
               begin
		  if (vth0Given)
		    vfb = TYPE * vth0 - phi - k1 * sqrtPhi;
                  else
                    vfb = -1.0;
               end
	     
             if (!vth0Given)
               vth0 = TYPE * (vfb + phi + k1 * sqrtPhi);
	     
	     k1ox = k1 * tox / TOXM;
             k2ox = k2 * tox / TOXM;
	     
             T1 = sqrt(`EPSSI / `EPSOX * tox * Xdep0);
             T0 = exp(-0.5 * dsub * leff / T1);
	     
             theta0vb0 = (T0 + 2.0 * T0 * T0);
	     
             T0 = exp(-0.5 * drout * leff / T1);
             T2 = (T0 + 2.0 * T0 * T0);
	     
             thetaRout = pdibl1 * T2 + pdibl2;
	     
	     /* vfbzb for capMod 1, 2 & 3 */
	     tmp  = sqrt(Xdep0);
             tmp1 = vbi - phi;
             tmp2 = factor1 * tmp;
	     
	     T0 = -0.5 * dvt1w * weff * leff / tmp2;
	     
             if (T0 > -`EXP_THRESHOLD)
               begin
		  T1 = exp(T0);
                  T2 = T1 * (1.0 + 2.0 * T1);
               end
             else
               begin
		  T1 = `MIN_EXP;
                  T2 = T1 * (1.0 + 2.0 * T1);
               end
	     
             T0 = dvt0w * T2;
             T2 = T0 * tmp1;
	     
             T0 = -0.5 * dvt1 * leff / tmp2;
	     
             if (T0 > -`EXP_THRESHOLD)
               begin
		  T1 = exp(T0);
                  T3 = T1 * (1.0 + 2.0 * T1);
               end
             else
	       begin
		  T1 = `MIN_EXP;
		  T3 = T1 * (1.0 + 2.0 * T1);
               end

	     T3 = dvt0 * T3 * tmp1;
	     
             T4 = tox * phi / (weff + w0);
	     
             T0 = sqrt(1.0 + nlx / leff);
             T5 = k1ox * (T0 - 1.0) * sqrtPhi + (kt1 + kt1l / leff) * (TRatio - 1.0);
	     
             tmp3  = TYPE * vth0 - T2 - T3 + k3 * T4 + T5;
             vfbzb = tmp3 - phi - k1 * sqrtPhi;

	     // End of vfbzb calculation

	     // process source/drain series resistance
	     drainConductance = sheetResistance  * drainSquares;
	     
             if (drainConductance > 0.0)
               drainConductance = 1.0 / drainConductance;
	     else
               drainConductance = 0.0;
                  
             sourceConductance = sheetResistance * sourceSquares;
	     
             if (sourceConductance > 0.0) 
               sourceConductance = 1.0 / sourceConductance;
	     else
               sourceConductance = 0.0;
	     
	     Nvtm = vtm * jctEmissionCoeff;
	     
             if ((sourceArea <= 0.0) && (sourcePerimeter <= 0.0))
               SourceSatCurrent = 1.0e-14;
             else
               SourceSatCurrent =   sourceArea      * jctTempSatCurDensity
                                  + sourcePerimeter * jctSidewallTempSatCurDensity;
             
             if ((SourceSatCurrent > 0.0) && (IJTH > 0.0))
	       begin
		  vjsm = Nvtm * ln(IJTH / SourceSatCurrent + 1.0);
		  
		  if (VERSION == 3.24)
		    IsEvjsm = SourceSatCurrent * exp(vjsm / Nvtm);
	       end
             
             if ((drainArea <= 0.0) && (drainPerimeter <= 0.0))
               DrainSatCurrent = 1.0e-14;
             else
               DrainSatCurrent =   drainArea      * jctTempSatCurDensity
                                 + drainPerimeter * jctSidewallTempSatCurDensity;
             
             if ((DrainSatCurrent > 0.0) && (IJTH > 0.0))
	       begin
		  vjdm = Nvtm * ln(IJTH / DrainSatCurrent + 1.0);
		  
		  if (VERSION == 3.24)
		    IsEvjdm = DrainSatCurrent * exp(vjdm / Nvtm);
	       end

         end
	
      	//*********************************//
      	//****** End of initial_step ******//
      	//*********************************//
	
	vbs  = TYPE * V(bulk,   sourcep);
	vgs  = TYPE * V(gate,   sourcep);
	vds  = TYPE * V(drainp, sourcep);
`ifdef NQSMOD
	qdef = TYPE * V(q);
`endif
	
        vbd = vbs - vds;
        vgd = vgs - vds;
        vgb = vgs - vbs;

	temp = $temperature;

	// Source/drain junction diode DC model begins	
	if (SourceSatCurrent <= 0.0)
	  begin
	     gbs = GMIN;
             cbs = gbs * vbs;
          end
	else
	  begin
	     if (IJTH == 0.0)
	       begin
		  evbs = exp(vbs / Nvtm);
                  gbs  = SourceSatCurrent * evbs / Nvtm  + GMIN;
                  cbs  = SourceSatCurrent * (evbs - 1.0) + GMIN * vbs; 
               end
             else
               begin
		  if (vbs < vjsm)
                    begin
		       evbs = exp(vbs / Nvtm);
                       gbs  = SourceSatCurrent * evbs / Nvtm  + GMIN;
                       cbs  = SourceSatCurrent * (evbs - 1.0) + GMIN * vbs;
                    end
                  else
                    begin
		       if (VERSION == 3.24)
			 begin
			    T0  = IsEvjsm / Nvtm;
			    cbs = IsEvjsm - SourceSatCurrent
				  + GMIN * vbs + T0 * (vbs - vjsm);
			 end
		       else
			 begin
			    T0  = (SourceSatCurrent + IJTH) / Nvtm;
			    cbs = IJTH + GMIN * vbs + T0 * (vbs - vjsm);
			 end
		       
                       gbs = T0   + GMIN;
                    end
               end
          end

	if (DrainSatCurrent <= 0.0)
	  begin
	     gbd = GMIN;
             cbd = gbd * vbd;
          end
	else
          begin
	     if (IJTH == 0.0)
               begin
		  evbd = exp(vbd / Nvtm);
                  gbd  = DrainSatCurrent * evbd / Nvtm  + GMIN;
                  cbd  = DrainSatCurrent * (evbd - 1.0) + GMIN * vbd;
               end
             else
               begin
		  if (vbd < vjdm)
                    begin
		       evbd = exp(vbd / Nvtm);
                       gbd  = DrainSatCurrent * evbd / Nvtm  + GMIN;
                       cbd  = DrainSatCurrent * (evbd - 1.0) + GMIN * vbd;
                    end
                  else
                    begin
		       if (VERSION == 3.24)
			 begin
			    T0  = IsEvjdm / Nvtm;
			    cbd = IsEvjdm - DrainSatCurrent
                                  + GMIN * vbd + T0 * (vbd - vjdm);
			 end
		       else
			 begin
			    T0  = (DrainSatCurrent + IJTH) / Nvtm;
		            cbd = IJTH + GMIN * vbd + T0 * (vbd - vjdm);
			 end
	               gbd = T0 + GMIN;
		       
                    end
               end
          end
	// End of diode DC model	

	if (vds >= 0.0)
	  begin   /* normal mode */
             mode = 1;
             Vds  = vds;
             Vgs  = vgs;
             Vbs  = vbs;
          end
	else
	  begin   /* inverse mode */
             mode = -1;
             Vds  = -vds;
             Vgs  = vgd;
             Vbs  = vbd;
          end

	T0     = Vbs - vbsc - 0.001;
	T1     = sqrt(T0 * T0 - 0.004 * vbsc);
	Vbseff = vbsc + 0.5 * (T0 + T1);

	// Added to avoid the possible numerical problems due to computer accuracy. 
	// (See comments for diffVds)
	if (Vbseff < Vbs)
	  Vbseff = Vbs;

	if (Vbseff > 0.0)
	  begin
	     T0       = phi / (phi + Vbseff);
             Phis     = phi * T0;
             sqrtPhis = phis3 / (phi + 0.5 * Vbseff);
          end
	else
	  begin
	     Phis     = phi - Vbseff;
             sqrtPhis = sqrt(Phis);
          end

	Xdep = Xdep0 * sqrtPhis / sqrtPhi;

	Leff = leff;
	Vtm  = vtm;
	
	/*** Vth Calculation ***/
	T3 = sqrt(Xdep);
        V0 = vbi - phi;
	
        T0 = dvt2 * Vbseff;
	
        if (T0 >= - 0.5)
	  begin
	     T1 = 1.0 + T0;
	     T2 = dvt2;
	  end
	else /* Added to avoid any discontinuity problems caused by dvt2 */ 
	  begin
	     T4 = 1.0 / (3.0 + 8.0 * T0);
	     T1 = (1.0 + 3.0 * T0) * T4; 
	     T2 = dvt2 * T4 * T4;
	  end

	lt1 = factor1 * T3 * T1;
	
	T0 = dvt2w * Vbseff;
	
	if (T0 >= - 0.5)
	  begin
	     T1 = 1.0 + T0;
	     T2 = dvt2w;
	  end
	else /* Added to avoid any discontinuity problems caused by dvt2w */ 
	  begin
	     T4 = 1.0 / (3.0 + 8.0 * T0);
	     T1 = (1.0 + 3.0 * T0) * T4; 
	     T2 = dvt2w * T4 * T4;
	  end
	
        ltw = factor1 * T3 * T1;

	T0 = -0.5 * dvt1 * Leff / lt1;
	
        if (T0 > -`EXP_THRESHOLD)
          begin
	     T1     = exp(T0);
             Theta0 = T1 * (1.0 + 2.0 * T1);
          end
        else
          begin
	     T1     = `MIN_EXP;
             Theta0 = T1 * (1.0 + 2.0 * T1);
          end

	thetavth = dvt0 * Theta0;
        Delt_vth = thetavth * V0;

	T0 = -0.5 * dvt1w * weff * Leff / ltw;
	
        if (T0 > -`EXP_THRESHOLD)
          begin
	     T1 = exp(T0);
             T2 = T1 * (1.0 + 2.0 * T1);
          end
        else
          begin
	     T1 = `MIN_EXP;
             T2 = T1 * (1.0 + 2.0 * T1);
          end
	
        T0 = dvt0w * T2;
        T2 = T0 * V0;

	TempRatio = temp / tnom - 1.0;
        T0 = sqrt(1.0 + nlx / Leff);
        T1 = k1ox * (T0 - 1.0) * sqrtPhi 
	     + (kt1 + kt1l / Leff + kt2 * Vbseff) * TempRatio;
	
        tmp2 = tox *phi / (weff + w0);

	T3 = eta0 + etab * Vbseff;
	if (T3 < 1.0e-4) /* avoid  discontinuity problems caused by etab */ 
	  begin
	     T9 = 1.0 / (3.0 - 2.0e4 * T3);
	     T3 = (2.0e-4 - T3) * T9;
	     T4 = T9 * T9;
	  end
	else
	  T4 = 1.0;
	
	dDIBL_Sft_dVd = T3 * theta0vb0;
        DIBL_Sft      = dDIBL_Sft_dVd * Vds;
	
        Vth = TYPE * vth0 - k1 * sqrtPhi 
	                  + k1ox * sqrtPhis
	                  - k2ox * Vbseff
	                  - Delt_vth - T2
	                  + (k3 + k3b * Vbseff) * tmp2 + T1 - DIBL_Sft;

	/*** end of Vth calculation ***/

	/* Calculate n */
        tmp2 = nfactor * `EPSSI / Xdep;
        tmp3 = cdsc + cdscb * Vbseff + cdscd * Vds;
	tmp4 = (tmp2 + tmp3 * Theta0 + cit) / cox;

	if (tmp4 >= -0.5)
	  n = 1.0 + tmp4;
	else
	  begin /* avoid  discontinuity problems caused by tmp4 */
	     T0 = 1.0 / (3.0 + 8.0 * tmp4);
	     n  = (1.0 + 3.0 * tmp4) * T0;
	  end
	
	/* Poly Gate Si Depletion Effect */
	T0 = vfb + phi;
	
        if ((ngate > 1.0e18) && (ngate < 1.0e25) && (Vgs > T0))
	  begin /* added to avoid the problem caused by ngate */
             T1 = 1.0e6 * `Charge_q * `EPSSI * ngate / (cox * cox);
             T4 = sqrt(1.0 + 2.0 * (Vgs - T0) / T1);
             T2 = T1 * (T4 - 1.0);
             T3 = 0.5 * T2 * T2 / T1; /* T3 = Vpoly */
             T7 = 1.12 - T3 - 0.05;
             T6 = sqrt(T7 * T7 + 0.224);
             T5 = 1.12 - 0.5 * (T7 + T6);
             Vgs_eff = Vgs - T5;
          end
	else
	  Vgs_eff = Vgs;
	
	Vgst = Vgs_eff - Vth;
	
	/* Effective Vgst (Vgsteff) Calculation */
        T10 = 2.0 * n * Vtm;
        VgstNVt = Vgst / T10;
        ExpArg = (2.0 * voff - Vgst) / T10;
	
	/* MCJ: Very small Vgst */
        if (VgstNVt > `EXP_THRESHOLD)
	  begin
	     Vgsteff = Vgst;
	  end
	else if (ExpArg > `EXP_THRESHOLD)
	  begin
	     T0 = (Vgst - voff) / (n * Vtm);
	     ExpVgst = exp(T0);
	     Vgsteff = Vtm * cdep0 / cox * ExpVgst;
	  end
	else
	  begin
	     ExpVgst = exp(VgstNVt);
             T1 = T10 * ln(1.0 + ExpVgst);
             T2 = 1.0 - T10 * (-cox / (Vtm * cdep0) * exp(ExpArg));
             Vgsteff = T1 / T2;
	  end

	/* Calculate Effective Channel Geometry */
        T9   = sqrtPhis - sqrtPhi;
        Weff = weff - 2.0 * (dwg * Vgsteff + dwb * T9); 

        if (Weff < 2.0e-8) /* to avoid the discontinuity problem due to Weff*/
	  begin
	     T0 = 1.0 / (6.0e-8 - 2.0 * Weff);
	     Weff = 2.0e-8 * (4.0e-8 - Weff) * T0;
          end
	
        T0 = prwg * Vgsteff + prwb * T9;
	if (T0 >= -0.9)
	  Rds = rds0 * (1.0 + T0);
	else /* to avoid the discontinuity problem due to prwg and prwb*/
	  begin
	     T1  = 1.0 / (17.0 + 20.0 * T0);
	     Rds = rds0 * (0.8 + T0) * T1;
	  end

	/* Calculate Abulk */
        T1 = 0.5 * k1ox / sqrtPhis;
	
        T9   = sqrt(xj * Xdep);
        tmp1 = Leff + 2.0 * T9;
        T5   = Leff / tmp1; 
        tmp2 = a0 * T5;
        tmp3 = weff + b1; 
        tmp4 = b0 / tmp3;
        T2   = tmp2 + tmp4;
        T6   = T5 * T5;
        T7   = T5 * T6;
	
        Abulk0 = 1.0 + T1 * T2; 
	
        T8 = ags * a0 * T7;
	
        Abulk = Abulk0 + (-T1 * T8) * Vgsteff; 
	
        if (Abulk0 < 0.1) /* added to avoid the problems caused by Abulk0 */
	  begin
	     T9 = 1.0 / (3.0 - 20.0 * Abulk0);
	     Abulk0 = (0.2 - Abulk0) * T9;
	  end

	if (Abulk < 0.1) /* added to avoid the problems caused by Abulk */
	  begin
	     T9 = 1.0 / (3.0 - 20.0 * Abulk);
	     Abulk = (0.2 - Abulk) * T9;
	  end
	
        T2 = keta * Vbseff;
	if (T2 >= -0.9)
	  T0 = 1.0 / (1.0 + T2);
	else /* added to avoid the problems caused by Keta */
	  begin
	     T1 = 1.0 / (0.8 + T2);
	     T0 = (17.0 + 20.0 * T2) * T1;
	  end
	
	Abulk  = T0 * Abulk;
	Abulk0 = T0 * Abulk0;
	
	/* Mobility calculation */
        if (MOBMOD == 1)
	    begin
	       T0 = Vgsteff + Vth + Vth;
               T2 = ua + uc * Vbseff;
               T3 = T0 / tox;
               T5 = T3 * (T2 + ub * T3);
            end
	else if (MOBMOD == 2)
	     T5 = Vgsteff / tox * (ua + uc * Vbseff + ub * Vgsteff / tox);
	else
	  begin
	     T0 = Vgsteff + Vth + Vth;
             T2 = 1.0 + uc * Vbseff;
             T3 = T0 / tox;
             T4 = T3 * (ua + ub * T3);
	     T5 = T4 * T2;
          end
	
	if (T5 >= -0.8)
	  Denomi = 1.0 + T5;
	else /* Added to avoid the discontinuity problem caused by ua and ub*/ 
	  Denomi = (0.6 + T5) * (1.0 / (7.0 + 10.0 * T5));
	
        ueff = u0temp / Denomi;

	/* Saturation Drain Voltage  Vdsat */
        WVCox    = Weff * vsattemp * cox;
        WVCoxRds = WVCox * Rds; 
	
        Esat = 2.0 * vsattemp / ueff;
        EsatL = Esat * Leff;
  	
	/* Sqrt() */
	if (a1 == 0.0)
	  Lambda = a2;
	else if (a1 > 0.0)
	  /* Added to avoid the discontinuity problem caused by a1 and a2 (Lambda) */
	  begin
	     T0 = 1.0 - a2;
	     T1 = T0 - a1 * Vgsteff - 0.0001;
	     T2 = sqrt(T1 * T1 + 0.0004 * T0);
	     Lambda = a2 + T0 - 0.5 * (T1 + T2);
	  end
	else
	  begin
	     T1 = a2 + a1 * Vgsteff - 0.0001;
	     T2 = sqrt(T1 * T1 + 0.0004 * a2);
	     Lambda = 0.5 * (T1 + T2);
	  end
	
        Vgst2Vtm = Vgsteff + 2.0 * Vtm;
	
        if ((Rds == 0.0) && (Lambda == 1.0))
          begin
	     T0 = 1.0 / (Abulk * EsatL + Vgst2Vtm);
             tmp1 = 0.0;
	     T1 = T0 * T0;
	     T2 = Vgst2Vtm * T0;
             T3 = EsatL * Vgst2Vtm;
             Vdsat = T3 * T0;
          end
        else
          begin
             T9 = Abulk * WVCoxRds;
	     T7 = Vgst2Vtm * T9;
             T6 = Vgst2Vtm * WVCoxRds;
             T0 = 2.0 * Abulk * (T9 - 1.0 + 1.0 / Lambda); 
	     
             T1 = Vgst2Vtm * (2.0 / Lambda - 1.0) + Abulk * EsatL + 3.0 * T7;
             T2 = Vgst2Vtm * (EsatL + 2.0 * T6);
             T3 = sqrt(T1 * T1 - 2.0 * T0 * T2);
	     
             Vdsat = (T1 - T3) / T0;
          end
	
        vdsat = Vdsat;

	/* Effective Vds (Vdseff) Calculation */
        T1 = Vdsat - Vds - delta;
        T2 = sqrt(T1 * T1 + 4.0 * delta * Vdsat);
        Vdseff = Vdsat - 0.5 * (T1 + T2);
	
        if ((Vds == 0.0) && (VERSION == 3.24))
            Vdseff = 0.0;
	     
	/* Calculate VAsat */
        tmp4 = 1.0 - 0.5 * Abulk * Vdsat / Vgst2Vtm;
        T9 = WVCoxRds * Vgsteff;
        T0 = EsatL + Vdsat + 2.0 * T9 * tmp4;
        T9 = WVCoxRds * Abulk; 
        T1 = 2.0 / Lambda - 1.0 + T9; 
	
        Vasat = T0 / T1;
	
	if (Vdseff > Vds)
	  Vdseff = Vds; /* This code is added to fixed the problem
			 caused by computer precision when
			 Vds is very close to Vdseff. */
        diffVds = Vds - Vdseff;

	/* Calculate VACLM */
        if ((pclm > 0.0) && (diffVds > 1.0e-10))
	  begin
	     T0 = 1.0 / (pclm * Abulk * litl);
	     T2 = Vgsteff / EsatL;
             T1 = Leff * (Abulk + T2); 
	     T9 = T0 * T1;
             VACLM = T9 * diffVds;
          end
	else
	  VACLM = `MAX_EXP;

	/* Calculate VADIBL */
        if (thetaRout > 0.0)
	  begin
	     T0 = Vgst2Vtm *  Abulk * Vdsat;
             T1 = Vgst2Vtm + (Abulk * Vdsat);
             VADIBL = (Vgst2Vtm - T0 / T1) / thetaRout;
	     
	     T7 = pdiblb * Vbseff;
	     if (T7 >= -0.9)
	       begin
		  T3 = 1.0 / (1.0 + T7);
                  VADIBL = T3 * VADIBL;
	       end
	     else
	       /* Added to avoid the discontinuity problem caused by pdiblcb */
	       begin
		  T4 = 1.0 / (0.8 + T7);
		  T3 = (17.0 + 20.0 * T7) * T4;
                  VADIBL = T3 * VADIBL;
	       end
          end
	else
	  VADIBL = `MAX_EXP;

	/* Calculate VA */
	T9 = pvag / EsatL * Vgsteff;
	
	if (T9 > -0.9)
	  T0 = 1.0 + T9;
	else /* Added to avoid the discontinuity problems caused by pvag */
	  T0 = (0.8 + T9) * 1.0 / (17.0 + 20.0 * T9);
        
        T1   = VACLM * VADIBL / (VACLM + VADIBL);
        Va   = Vasat + T0 * T1;

	/* Calculate VASCBE */
	if (pscbe2 > 0.0)
	  begin
	     if ( diffVds > (pscbe1 * litl / `EXP_THRESHOLD) )
	       VASCBE = Leff * exp(pscbe1 * litl / diffVds) / pscbe2;
	     else
	       VASCBE = `MAX_EXP * Leff / pscbe2;
	  end
	else
	  VASCBE = `MAX_EXP;
	  
	/* Calculate Ids */
        CoxWovL = cox * Weff / Leff;
        beta    = ueff * CoxWovL;
	
        fgche1 = Vgsteff * (1.0 - 0.5 * Abulk * Vdseff / Vgst2Vtm);
        fgche2 = 1.0 + (Vdseff / EsatL);
 	
        gche = beta * fgche1 / fgche2;
        Idl  = gche * Vdseff / (1.0 + gche * Rds);

        Idsa = Idl  * (1.0 + (diffVds / Va));
        Ids  = Idsa * (1.0 + (diffVds / VASCBE));
	
	/* Substrate current begins */
        tmp = alpha0 + alpha1 * Leff;
	
        if ((tmp <= 0.0) || (beta0 <= 0.0))
	  Isub = 0.0;
	else
	  begin
	     T2 = tmp / Leff;
	     
	     if (diffVds > beta0 / `EXP_THRESHOLD)
	       T1 = T2 * diffVds * exp(-beta0 / diffVds);
	     else
	       T1 = T2 * `MIN_EXP * diffVds;
	     
             Isub = T1 * Idsa;
          end

	cdrain = Ids;
	csub   = Isub;

	// End of I-V model


	// C-V model
		
	// *** Depletion capacitance calculation ***
	
	/*  charge storage elements
         *  bulk-drain and bulk-source depletion capacitances
         *  czbd : zero bias drain junction capacitance
         *  czbs : zero bias source junction capacitance
         *  czbdsw: zero bias drain junction sidewall capacitance
	 along field oxide
         *  czbssw: zero bias source junction sidewall capacitance
	 along field oxide
	 *  czbdswg: zero bias drain junction sidewall capacitance
	 along gate side
	 *  czbsswg: zero bias source junction sidewall capacitance
	 along gate side
         */
	if (VERSION == 3.24)
	  begin
	     czbd = unitAreaTempJctCap * drainArea;
             czbs = unitAreaTempJctCap * sourceArea;
	  end
	else
	  begin
             czbd = unitAreaJctCap * drainArea;
             czbs = unitAreaJctCap * sourceArea;
	  end
	
        if (drainPerimeter < weff)
          begin
	     if (VERSION == 3.24)
	       czbdswg = unitLengthGateSidewallTempJctCap * drainPerimeter;
	     else
               czbdswg = unitLengthGateSidewallJctCap * drainPerimeter;
	     czbdsw  = 0.0;
          end
        else
          begin
	     if (VERSION == 3.24)
	       begin
		  czbdsw  = unitLengthSidewallTempJctCap * (drainPerimeter - weff);
		  czbdswg = unitLengthGateSidewallTempJctCap * weff;
	       end
	     else
	       begin
		  czbdsw  = unitLengthSidewallJctCap * (drainPerimeter - weff);
		  czbdswg = unitLengthGateSidewallJctCap *  weff;
	       end
          end
	
        if (sourcePerimeter < weff)
          begin
             czbssw  = 0.0;
	     if (VERSION == 3.24)
	       czbsswg = unitLengthGateSidewallTempJctCap * sourcePerimeter;
	     else
	       czbsswg = unitLengthGateSidewallJctCap * sourcePerimeter;
          end
        else
          begin
	     if (VERSION == 3.24)
	       begin
		  czbssw  = unitLengthSidewallTempJctCap * (sourcePerimeter - weff);
		  czbsswg = unitLengthGateSidewallTempJctCap * weff;
	       end
	     else
	       begin
		  czbssw  = unitLengthSidewallJctCap * (sourcePerimeter - weff);
		  czbsswg = unitLengthGateSidewallJctCap * weff;
	       end
	  end
	
        mj    = bulkJctBotGradingCoeff;
        mjsw  = bulkJctSideGradingCoeff;
	mjswg = bulkJctGateSideGradingCoeff;


	qbs = 0.0;
	qbd = 0.0;
	
	/* Source Bulk Junction */
	if (vbs == 0.0)
	  begin
	     qbs   = 0.0;
             capbs = czbs + czbssw + czbsswg;
	  end
	else if (vbs < 0.0)
	  begin
	     if (czbs > 0.0)
	       begin
		  arg = 1.0 - vbs / PhiB;
		  
		  if (mj == 0.5)
                    sarg = 1.0 / sqrt(arg);
		  else
                    sarg = exp(-mj * ln(arg));
		  
                  qbs   = PhiB * czbs * (1.0 - arg * sarg) / (1.0 - mj);
		  capbs = czbs * sarg;
	       end
	     else
	       begin
		  qbs   = 0.0;
		  capbs = 0.0;
	       end // else: !if(czbs > 0.0)
	     
	     if (czbssw > 0.0)
	       begin
		  arg = 1.0 - vbs / PhiBSW;
		  
		  if (mjsw == 0.5)
                    sarg = 1.0 / sqrt(arg);
		  else
                    sarg = exp(-mjsw * ln(arg));
		  
                  qbs   = qbs   + PhiBSW * czbssw * (1.0 - arg * sarg) / (1.0 - mjsw);
                  capbs = capbs + czbssw * sarg;
	       end // if (czbssw > 0.0)
	     
	     if (czbsswg > 0.0)
	       begin
		  arg = 1.0 - vbs / PhiBSWG;
		  
		  if (mjswg == 0.5)
                    sarg = 1.0 / sqrt(arg);
		  else
                    sarg = exp(-mjswg * ln(arg));
		  
                  qbs   = qbs   + PhiBSWG * czbsswg * (1.0 - arg * sarg) / (1.0 - mjswg);
                  capbs = capbs + czbsswg * sarg;
	       end
          end
	else
	  begin
	     T0 = czbs + czbssw + czbsswg;
             T1 = vbs * (czbs * mj / PhiB + czbssw * mjsw / PhiBSW + czbsswg * mjswg / PhiBSWG);
	     
             qbs   = vbs * (T0 + 0.5 * T1);
             capbs = T0 + T1;
	  end
	
	/* Drain Bulk Junction */
	if (vbd == 0.0)
	  begin
	     qbd   = 0.0;
             capbd = czbd + czbdsw + czbdswg;
	  end
	else if (vbd < 0.0)
	  begin
	     if (czbd > 0.0)
	       begin
		  arg = 1.0 - vbd / PhiB;
		  
		  if (mj == 0.5)
                    sarg = 1.0 / sqrt(arg);
		  else
                    sarg = exp(-mj * ln(arg));
		  
                  qbd   = PhiB * czbd * (1.0 - arg * sarg) / (1.0 - mj);
                  capbd = czbd * sarg;
	       end
	     else
	       begin
		  qbd   = 0.0;
                  capbd = 0.0;
	       end // else: !if(czbd > 0.0)
	     
	     if (czbdsw > 0.0)
	       begin
		  arg = 1.0 - vbd / PhiBSW;
		  
		  if (mjsw == 0.5)
                    sarg = 1.0 / sqrt(arg);
		  else
                    sarg = exp(-mjsw * ln(arg));
		  
                  qbd   = qbd   + PhiBSW * czbdsw * (1.0 - arg * sarg) / (1.0 - mjsw);
                  capbd = capbd + czbdsw * sarg;
	       end // if (czbdsw > 0.0)
	     
	     if (czbdswg > 0.0)
	       begin
		  arg = 1.0 - vbd / PhiBSWG;
		  
		  if (mjswg == 0.5)
                    sarg = 1.0 / sqrt(arg);
		  else
                    sarg = exp(-mjswg * ln(arg));
		  
                  qbd   = qbd   + PhiBSWG * czbdswg * (1.0 - arg * sarg) / (1.0 - mjswg);
                  capbd = capbd + czbdswg * sarg;
	       end
          end
	else
	  begin
	     T0 = czbd + czbdsw + czbdswg;
             T1 = vbd * (czbd * mj / PhiB + czbdsw * mjsw / PhiBSW + czbdswg * mjswg / PhiBSWG);
	     
             qbd   = vbd * (T0 + 0.5 * T1);
             capbd = T0 + T1; 
          end

	// *** Intrinsic charge calculation ***
	// qdrn, qgate, qsrc, qbulk
	
	qgate = 0.0;
	qdrn  = 0.0;
	qsrc  = 0.0;
	qbulk = 0.0;
	
        if (XPART < 0)
	  begin
	     qgate = 0.0;
	     qdrn  = 0.0;
	     qsrc  = 0.0;
	     qbulk = 0.0;
`ifdef NQSMOD
             gtau  = 0.0;
`endif
          end
	else if (CAPMOD == 0)
	  begin
             if (Vbseff < 0.0)
	       Vbseff = Vbs;
	     else
	       Vbseff = phi - Phis;
	     
             Vfb  = vfbcv;
             Vth  = Vfb + phi + k1ox * sqrtPhis;
	     
             Vgst = Vgs_eff - Vth;
	     
             CoxWL = cox * weffCV * leffCV;
             Arg1  = Vgs_eff - Vbseff - Vfb;
	     
	     if (Arg1 <= 0.0)
	       begin
		  qgate = CoxWL * Arg1;
                  qbulk = -qgate;
                  qdrn = 0.0;
		  
                  qinv = 0.0;
               end
	     else if (Vgst <= 0.0)
	       begin
		  T1 = 0.5 * k1ox;
	          T2 = sqrt(T1 * T1 + Arg1);
	          qgate = CoxWL * k1ox * (T2 - T1);
                  qbulk = -qgate;
                  qdrn = 0.0;
		  
                  qinv = 0.0;
               end
	     else
	       begin
		  One_Third_CoxWL = CoxWL / 3.0;
                  Two_Third_CoxWL = 2.0 * One_Third_CoxWL;
		  
                  AbulkCV = Abulk0 * abulkCVfactor;
	          Vdsat = Vgst / AbulkCV;

                  if (XPART > 0.5)
		    begin
		       /* 0/100 Charge partition model */
		       if (Vdsat <= Vds)
			 begin   /* saturation region */
	                    T1 = Vdsat / 3.0;
	                    qgate = CoxWL * (Vgs_eff - Vfb - phi - T1);
	                    T2 = -Two_Third_CoxWL * Vgst;
	                    qbulk = -(qgate + T2);
	                    qdrn = 0.0;
			    
                            qinv = -(qgate + qbulk);
			 end
		       else
			 begin
			    /* linear region */
			    Alphaz = Vgst / Vdsat;
	                    T1 = 2.0 * Vdsat - Vds;
	                    T2 = Vds / (3.0 * T1);
	                    T3 = T2 * Vds;
	                    T9 = 0.25 * CoxWL;
	                    T4 = T9 * Alphaz;
	                    T7 = 2.0 * Vds - T1 - 3.0 * T3;
	                    T8 = T3 - T1 - 2.0 * Vds;
	                    qgate = CoxWL * (Vgs_eff - Vfb - phi - 0.5 * (Vds - T3));
	                    T10 = T4 * T8;
	                    qdrn = T4 * T7;
	                    qbulk = -(qgate + qdrn + T10);
			    
                            qinv = -(qgate + qbulk);
			  end
		    end // if (XPART > 0.5)
		  
		  else if (XPART < 0.5)
		    
		    begin   /* 40/60 Charge partition model */
		       
		       if (Vds >= Vdsat)
			 begin  /* saturation region */
	                    T1 = Vdsat / 3.0;
	                    qgate = CoxWL * (Vgs_eff - Vfb - phi - T1);
	                    T2 = -Two_Third_CoxWL * Vgst;
	                    qbulk = -(qgate + T2);
	                    qdrn = 0.4 * T2;
			    
                            qinv = -(qgate + qbulk);
			 end
		       else
			 begin   /* linear region  */
			    Alphaz = Vgst / Vdsat;
			    T1 = 2.0 * Vdsat - Vds;
			    T2 = Vds / (3.0 * T1);
			    T3 = T2 * Vds;
			    T9 = 0.25 * CoxWL;
			    T4 = T9 * Alphaz;
			    qgate = CoxWL * (Vgs_eff - Vfb - phi - 0.5 * (Vds - T3));
			    			    
			    T6 = 8.0 * Vdsat * Vdsat - 6.0 * Vdsat * Vds
				 + 1.2 * Vds * Vds;
			    T8 = T2 / T1;
			    T7 = Vds - T1 - T8 * T6;
			    qdrn = T4 * T7;

			    T7 = 2.0 * (T1 + T3);
			    qbulk = -(qgate - T4 * T7);

                            qinv = -(qgate + qbulk);
			 end // else: !if(Vds >= Vdsat)     
		    end // if (XPART < 0.5)
		  else
		    begin   /* 50/50 partitioning */
		       if (Vds >= Vdsat)
			 begin   /* saturation region */
	                    T1 = Vdsat / 3.0;
	                    qgate = CoxWL * (Vgs_eff - Vfb - phi - T1);
	                    T2 = -Two_Third_CoxWL * Vgst;
	                    qbulk = -(qgate + T2);
	                    qdrn = 0.5 * T2;
			    
                            qinv = -(qgate + qbulk);
			 end
		       else
		      	 begin   /* linear region */
			    Alphaz = Vgst / Vdsat;
			    T1 = 2.0 * Vdsat - Vds;
			    T2 = Vds / (3.0 * T1);
			    T3 = T2 * Vds;
			    T9 = 0.25 * CoxWL;
			    T4 = T9 * Alphaz;
			    qgate = CoxWL * (Vgs_eff - Vfb - phi - 0.5 * (Vds - T3));
			    
			    T7 = T1 + T3;
			    qdrn = -T4 * T7;
			    qbulk = - (qgate + qdrn + qdrn);

                            qinv = -(qgate + qbulk);
			 end
	            end
	       end
	  end
	else
	  begin
	     if (Vbseff < 0.0)
	       VbseffCV = Vbseff;
	     else
	       VbseffCV = phi - Phis;
	     
             CoxWL = cox * weffCV * leffCV;
	     
             /* Seperate VgsteffCV with noff and voffcv */
             Noff2   = n   * noff_param;
             T0      = Vtm * Noff2;
             VgstNVt = (Vgst - voffcv_param) / T0;
	     
	     if (VgstNVt > `EXP_THRESHOLD)
	       Vgsteff = Vgst - voffcv_param;
             else if (VgstNVt < -`EXP_THRESHOLD)
	       Vgsteff = T0 * ln(1.0 + `MIN_EXP);
             else
               begin
		  ExpVgst = exp(VgstNVt);
                  Vgsteff = T0 * ln(1.0 + ExpVgst);
	       end /* End of VgsteffCV - Weidong 5/1998 */

	     if (CAPMOD == 1)
	       begin
		  if (VERSION < 3.2)
                    Vfb = Vth - phi - k1ox * sqrtPhis;
                  else
                    Vfb = vfbzb;
		  
                  Arg1 = Vgs_eff - VbseffCV - Vfb - Vgsteff;
		  
                  if (Arg1 <= 0.0)
		    qgate = CoxWL * Arg1;
	          else
	            begin
		       T0 = 0.5 * k1ox;
		       T1 = sqrt(T0 * T0 + Arg1);
                       qgate = CoxWL * k1ox * (T1 - T0);
                    end
		  
	          qbulk = -qgate;

		  One_Third_CoxWL = CoxWL / 3.0;
                  Two_Third_CoxWL = 2.0 * One_Third_CoxWL;
		  
                  AbulkCV = Abulk0 * abulkCVfactor;
	          VdsatCV = Vgsteff / AbulkCV;
		  
                  if (VdsatCV < Vds)
	            begin
		       T0    = Vgsteff - VdsatCV / 3.0;
                       qgate = qgate + CoxWL * T0;
		       
		       T0    = VdsatCV - Vgsteff;
                       qbulk = qbulk + One_Third_CoxWL * T0;
		       
                       if (XPART > 0.5)
		         T0 = -Two_Third_CoxWL;
                       else if (XPART < 0.5)
		         T0 = -0.4 * CoxWL;
                       else
		         T0 = -One_Third_CoxWL;
		       
                       qsrc = T0 * Vgsteff;
                    end
		  else
	            begin
		       T0 = AbulkCV * Vds;
                       T1 = 12.0 * (Vgsteff - 0.5 * T0 + 1.0e-20);
		       T2 = Vds / T1;
		       T3 = T0 * T2;
		       
                       qgate = qgate + CoxWL * (Vgsteff - 0.5 * Vds + T3);
                       qbulk = qbulk + CoxWL * (1.0 - AbulkCV) * (0.5 * Vds - T3);
		       
                       if (XPART > 0.5)
			 begin   /* 0/100 Charge petition model */
                            qsrc = -CoxWL * (0.5 * Vgsteff + 0.25 * T0 - T0 * T0 / (T1 + T1));
			 end
		       else if (XPART < 0.5)
			 begin   /* 40/60 Charge petition model */
		            T2 = 0.5 * CoxWL / (T1 / 12.0 * T1 / 12.0);
		            T3 = Vgsteff * (2.0 * T0 * T0 / 3.0 + Vgsteff * (Vgsteff - 4.0 * T0 / 3.0))
			       - 2.0 * T0 * T0 * T0 / 15.0;
		            qsrc = -T2 * T3;
			 end
		       else
			 begin   /* 50/50 Charge petition model */
                            qsrc = -0.5 * (qgate + qbulk);
			 end
                    end
		  
		  qdrn = -(qgate + qbulk + qsrc);
                  qinv = -(qgate + qbulk);
	       end // if (CAPMOD == 1)
	     
	     else if (CAPMOD == 2)
	       begin
		  if (VERSION < 3.2)
		    Vfb = Vth - phi - k1ox * sqrtPhis;
                  else
                    Vfb = vfbzb;

                  V3 = Vfb - Vgs_eff + VbseffCV - `DELTA_3;
		  if (Vfb <= 0.0)
		    begin
		       T0 = sqrt(V3 * V3 - 4.0 * `DELTA_3 * Vfb);
		       T2 = -`DELTA_3 / T0;
		    end
		  else
		    begin
		       T0 = sqrt(V3 * V3 + 4.0 * `DELTA_3 * Vfb);
		       T2 = `DELTA_3 / T0;
		    end
		  
		  T1 = 0.5 * (1.0 + V3 / T0);
		  Vfbeff = Vfb - 0.5 * (V3 + T0);
		  Qac0 = CoxWL * (Vfbeff - Vfb);

                  T0 = 0.5 * k1ox;
		  T3 = Vgs_eff - Vfbeff - VbseffCV - Vgsteff;
		  
                  if (k1ox == 0.0)
                    begin
		       T1 = 0.0;
                       T2 = 0.0;
                    end
		  else if (T3 < 0.0)
		    begin
		       T1 = T0 + T3 / k1ox;
                       T2 = CoxWL;
		    end
		  else
		    begin
		       T1 = sqrt(T0 * T0 + T3);
                       T2 = CoxWL * T0 / T1;
		    end
		  
		  Qsub0 = CoxWL * k1ox * (T1 - T0);
		  
                  AbulkCV = Abulk0 * abulkCVfactor;
	          VdsatCV = Vgsteff / AbulkCV;

		  V4 = VdsatCV - Vds - `DELTA_4;
		  T0 = sqrt(V4 * V4 + 4.0 * `DELTA_4 * VdsatCV);
		  VdseffCV = VdsatCV - 0.5 * (V4 + T0);

		  /* Added to eliminate non-zero VdseffCV at Vds=0.0 */
                  if ((Vds == 0.0) && (VERSION == 3.24))
                    VdseffCV = 0.0;
		  
	          T0 = AbulkCV * VdseffCV;
                  T1 = 12.0 * (Vgsteff - 0.5 * T0 + 1e-20);
		  T2 = VdseffCV / T1;
		  T3 = T0 * T2;

                  T4 = (1.0 - 12.0 * T2 * T2 * AbulkCV);
                  T5 = (6.0 * T0 * (4.0 * Vgsteff - T0) / (T1 * T1) - 0.5);
                  T6 = 12.0 * T2 * T2 * Vgsteff;

                  qinoi = -CoxWL * (Vgsteff - 0.5 * T0 + AbulkCV * T3);
                  qgate =  CoxWL * (Vgsteff - 0.5 * VdseffCV + T3);

		  T7 = 1.0 - AbulkCV;
                  qbulk = CoxWL * T7 * (0.5 * VdseffCV - T3);

		  if (XPART > 0.5)
		    begin   /* 0/100 Charge petition model */
                       qsrc = -CoxWL * (0.5 * Vgsteff + 0.25 * T0 - T0 * T0 / (T1 + T1));
                    end
		  else if (XPART < 0.5)
		    begin   /* 40/60 Charge petition model */
		       T2 = 0.5 * CoxWL / (T1 / 12.0 * T1 / 12.0);
		       T3 = Vgsteff * (2.0 * T0 * T0 / 3.0 + Vgsteff * (Vgsteff - 4.0 * T0 / 3.0))
		          - 2.0 * T0 * T0 * T0 / 15.0;
		       qsrc = -T2 * T3;
                    end
		  else
		    begin   /* 50/50 Charge petition model */
                       qsrc = -0.5 * (qgate + qbulk);
		    end

		  qgate = qgate + Qac0 + Qsub0;
		  qbulk = qbulk - (Qac0 + Qsub0);
                  qdrn  = -(qgate + qbulk + qsrc);

                  qinv = qinoi;
	       end 
	     /* New Charge-Thickness capMod (CTM) begins - Weidong 7/1997 */
	      else if (CAPMOD == 3)
	      	begin
		   V3 = vfbzb - Vgs_eff + VbseffCV - `DELTA_3;
		   
		   if (vfbzb <= 0.0)
		     begin
			T0 = sqrt(V3 * V3 - 4.0 * `DELTA_3 * vfbzb);
		      	T2 = -`DELTA_3 / T0;
		     end
		   else
		     begin
			T0 = sqrt(V3 * V3 + 4.0 * `DELTA_3 * vfbzb);
		      	T2 = `DELTA_3 / T0;
		     end
		   
		   T1 = 0.5 * (1.0 + V3 / T0);
		   Vfbeff = vfbzb - 0.5 * (V3 + T0);
		   
                   Cox = cox;
                   Tox = 1.0e8 * tox;
                   T0 = (Vgs_eff - VbseffCV - vfbzb) / Tox;
		   
                   tmp = T0 * acde;
		   
                   if ((-`EXP_THRESHOLD < tmp) && (tmp < `EXP_THRESHOLD))
		     Tcen = ldeb * exp(tmp);
                   else if (tmp <= -`EXP_THRESHOLD)
		     Tcen = ldeb * `MIN_EXP;
                   else
                     Tcen = ldeb * `MAX_EXP;
		   
		   LINK = 1.0e-3 * tox;
                   V3   = ldeb - Tcen - LINK;
                   V4   = sqrt(V3 * V3 + 4.0 * LINK * ldeb);
                   Tcen = ldeb - 0.5 * (V3 + V4);
                   T1   = 0.5 * (1.0 + V3 / V4);
		   
                   Ccen     = `EPSSI / Tcen;
                   T2       = Cox / (Cox + Ccen);
                   Coxeff   = T2 * Ccen;
                   T3       = -Ccen / Tcen;
                   CoxWLcen = CoxWL * Coxeff / Cox;
		   
                   Qac0   = CoxWLcen * (Vfbeff - vfbzb);
		   
                   T0 = 0.5 * k1ox;
                   T3 = Vgs_eff - Vfbeff - VbseffCV - Vgsteff;
		   
                   if (k1ox == 0.0)
                     begin
			T1 = 0.0;
                      	T2 = 0.0;
                     end
                   else if (T3 < 0.0)
                     begin
			T1 = T0 + T3 / k1ox;
                      	T2 = CoxWLcen;
                     end
                   else
                     begin
			T1 = sqrt(T0 * T0 + T3);
                      	T2 = CoxWLcen * T0 / T1;
                     end

		   Qsub0 = CoxWLcen * k1ox * (T1 - T0);
		   
		   /* Gate-bias dependent delta Phis begins */
		   if (k1ox <= 0.0)
		     begin
			Denomi = 0.25 * moin * Vtm;
                      	T0     = 0.5 * sqrtPhi;
		     end
		   else
		     begin
			Denomi = moin * Vtm * k1ox * k1ox;
                      	T0     = k1ox * sqrtPhi;
		     end
		   
                   T1 = 2.0 * T0 + Vgsteff;
		   
		   DeltaPhi = Vtm * ln(1.0 + T1 * Vgsteff / Denomi);
		   /* End of delta Phis */
		   
                   T3  = 4.0 * (Vth - vfbzb - phi);
                   Tox2 = Tox + Tox;
		   
                   if (T3 >= 0.0)
                     T0 = (Vgsteff + T3) / Tox2;
                   else
                     T0 = (Vgsteff + 1.0e-20) / Tox2;
		   
                   tmp  = exp(0.7 * ln(T0));
                   T1   = 1.0 + tmp;
                   T2   = 0.7 * tmp / (T0 * Tox2);
                   Tcen = 1.9e-9 / T1;
		   
		   Ccen     = `EPSSI / Tcen;
		   T0       = Cox / (Cox + Ccen);
		   Coxeff   = T0 * Ccen;
		   T1       = -Ccen / Tcen;
		   CoxWLcen = CoxWL * Coxeff / Cox;

		   AbulkCV = Abulk0 * abulkCVfactor;
		   VdsatCV = (Vgsteff - DeltaPhi) / AbulkCV;
		   V4 = VdsatCV - Vds - `DELTA_4;
		   T0 = sqrt(V4 * V4 + 4.0 * `DELTA_4 * VdsatCV);
		   VdseffCV = VdsatCV - 0.5 * (V4 + T0);
		   T1 = 0.5 * (1.0 + V4 / T0);
		   T2 = `DELTA_4 / T0;
		   T3 = (1.0 - T1 - T2) / AbulkCV;

		   if ((Vds == 0.0) && (VERSION == 3.24))
                     VdseffCV = 0.0;
		     
                   T0 = AbulkCV * VdseffCV;
		   T1 = Vgsteff - DeltaPhi;
                   T2 = 12.0 * (T1 - 0.5 * T0 + 1.0e-20);
                   T3 = T0 / T2;
                   T4 = 1.0 - 12.0 * T3 * T3;
                   T5 = AbulkCV * (6.0 * T0 * (4.0 * T1 - T0) / (T2 * T2) - 0.5);
		   T6 = T5 * VdseffCV / AbulkCV;
		   
                   qinoi = CoxWLcen * (T1 - T0 * (0.5 - T3));
		   qgate = qinoi;
		   
		   T7 = 1.0 - AbulkCV;

		   qbulk  = CoxWLcen * T7 * (0.5 * VdseffCV - T0 * VdseffCV / T2);

                   if (XPART > 0.5)
		     begin  /* 0/100 partition */
		      	qsrc   = -CoxWLcen * (T1 / 2.0 + T0 / 4.0 - 0.5 * T0 * T0 / T2);
		     end
		   else if (XPART < 0.5)
		     begin   /* 40/60 partition */
		       	T2 = T2 / 12.0;
		       	T3 = 0.5 * CoxWLcen / (T2 * T2);
                       	T4 = T1 * (2.0 * T0 * T0 / 3.0 + T1 
				   * (T1 - 4.0 * T0 / 3.0)) - 2.0 * T0 * T0 * T0 / 15.0;
		       	qsrc   = -T3 * T4;
                     end
		   else
		     begin   /* 50/50 partition */
                       	qsrc = -0.5 * qgate;
		     end
		   
		   qgate = qgate + Qac0 + Qsub0 - qbulk;
		   qbulk = qbulk - (Qac0 + Qsub0);
                   qdrn  = -(qgate + qbulk + qsrc);
		   
                   qinv  = -qinoi;
	      	end /* End of CTM */
          end
	// *** end of intrinsic charge calculation ***

`ifdef NQSMOD
        /* NQS (Mansun 11/1993) modified by Weidong & Min-Chie 1997-1998 */
        if (NQSMOD)
          begin
	     qcheq = -(qbulk + qgate);
	     
             gtau_drift = abs(tconst * qcheq) * ScalingFactor;
             gtau_diff  = 16.0 * u0temp * vtm / (leffCV * leffCV) * ScalingFactor;
	     
             gtau = gtau_drift + gtau_diff;
          end
`endif
		
 	qgdo = 0.0;
	qgso = 0.0;
	
	// *** overlap capacitance charge calculation ***
	if (CAPMOD == 0.0)
	  begin
	     qgdo = cgdo_param * vgd;
	     qgso = cgso_param * vgs;
	  end
	else if (CAPMOD == 1.0)
	  begin
	     if (vgd < 0.0)
	       begin
		  T1 = sqrt(1.0 - 4.0 * vgd / ckappa);
	          qgdo = cgdo_param * vgd - weffCV * 0.5 * cgdl * ckappa * (T1 - 1.0);
	       end
	     else
	       begin
	          qgdo = (weffCV * cgdl + cgdo_param) * vgd;
	       end

	     if (vgs < 0.0)
	       begin
		  T1 = sqrt(1.0 - 4.0 * vgs / ckappa);
	          qgso = cgso_param * vgs - weffCV * 0.5 * cgsl * ckappa * (T1 - 1.0);
	       end
	     else
	       begin
	          qgso = (weffCV * cgsl + cgso_param) * vgs;
	       end
	  end
	else
	  begin
	     T0 = vgd + `DELTA_1;
	     T1 = sqrt(T0 * T0 + 4.0 * `DELTA_1);
	     T2 = 0.5 * (T0 - T1);
	     
	     T3 = weffCV * cgdl;
	     T4 = sqrt(1.0 - 4.0 * T2 / ckappa);
	     qgdo = (cgdo_param + T3) * vgd - T3 * (T2 + 0.5 * ckappa * (T4 - 1.0));

	     T0 = vgs + `DELTA_1;
	     T1 = sqrt(T0 * T0 + 4.0 * `DELTA_1);
	     T2 = 0.5 * (T0 - T1);
	     T3 = weffCV * cgsl;
	     T4 = sqrt(1.0 - 4.0 * T2 / ckappa);
	     qgso = (cgso_param + T3) * vgs - T3 * (T2 + 0.5 * ckappa * (T4 - 1.0));
	  end

	// Add Overlap capacitance charges contribution to total node charge
	// according to mode and NQS model
	if (mode > 0)
          begin
	     if (NQSMOD == 0)
               begin
                  qgd = qgdo;
                  qgs = qgso;
                  qgb = cgbo_param * vgb;
		  
                  qgate = qgate + qgd + qgs + qgb;
                  qbulk = qbulk - qgb;
                  qdrn  = qdrn  - qgd;
                  qsrc  = -(qgate + qbulk + qdrn);

`ifdef NQSMOD
		  sxpart = 0.6;
                  dxpart = 0.4;
`endif
               end // if (NQSMOD == 0)
	     else
	       begin
		  CoxWL = cox * weffCV * leffCV;

`ifdef NQSMOD
		  if ( abs(qcheq) <= 1.0e-5 * CoxWL )
		    begin
		       if (XPART < 0.5)
		         dxpart = 0.4;
		       else if (XPART > 0.5)
		         dxpart = 0.0;
		       else
		         dxpart = 0.5;
		    end
		  else
		    dxpart = qdrn / qcheq;
		  
		  sxpart = 1.0 - dxpart;
`endif
		  
		  qgd   = qgdo;
                  qgs   = qgso;
                  qgb   = cgbo_param * vgb;
		  
                  qgate = qgd + qgs + qgb;
                  qbulk = -qgb;
                  qdrn  = -qgd;
                  qsrc  = -(qgate + qbulk + qdrn);
	       end
	  end // if (mode > 0)
	else
	  begin
	     if (NQSMOD == 0)
               begin
                  qgd = qgdo;
                  qgs = qgso;
                  qgb = cgbo_param * vgb;
		  
                  qgate = qgate + qgd + qgs + qgb;
                  qbulk = qbulk - qgb;
                  qsrc  = qdrn  - qgs;
                  qdrn  = -(qgate + qbulk + qsrc);

`ifdef NQSMOD		  
		  sxpart = 0.4;
                  dxpart = 0.6;
`endif
	       end // if (NQSMOD == 0)
	     else
	       begin
		  CoxWL = cox * weffCV * leffCV;

`ifdef NQSMOD		  
		  if ( abs(qcheq) <= 1.0e-5 * CoxWL )
		    begin
		       if (XPART < 0.5)
		         sxpart = 0.4;
		       else if (XPART > 0.5)
		         sxpart = 0.0;
		       else
		         sxpart = 0.5;
		    end
		  else
		    sxpart = qdrn / qcheq;
		  
		  dxpart = 1.0 - sxpart;
`endif
  
                  qgd = qgdo;
                  qgs = qgso;
                  qgb = cgbo_param * vgb;
		  
                  qgate = qgd + qgs + qgb;
                  qbulk = -qgb;
                  qsrc  = -qgs;
                  qdrn  = -(qgate + qbulk + qsrc);
	       end
	  end

`ifdef NQSMOD
	if (NQSMOD)
          begin
	     qcdump  = qdef * ScalingFactor;

	     cqdef  = ddt(qcdump);
	     cqcheq = ddt(qcheq);
          end
	
	if (analysis("static"))
	  begin
	     dxpart = (mode > 0) ? 0.4 : 0.6;
	     sxpart = 1.0 - dxpart;
	
	     if (NQSMOD)
               gtau = 16.0 * u0temp * vtm / leffCV / leffCV * ScalingFactor;
	     else
               gtau = 0.0;
	  end
`endif
	
	// Add depletion capacitance charge contribution
	Qdrn  = qdrn  - qbd;
	Qsrc  = qsrc  - qbs;
	Qbulk = qbulk + qbd + qbs;
	Qgate = qgate;

	if (mode > 0)
	  begin
	     I(drainp, sourcep) <+ TYPE * cdrain;
	     I(bulk,   drainp)  <+ TYPE * (cbd - csub);
	     I(bulk,   sourcep) <+ TYPE * cbs;
	  end
	else
	  begin
	     I(drainp, sourcep) <+ TYPE * (-cdrain);
	     I(bulk,   drainp)  <+ TYPE * cbd;
	     I(bulk,   sourcep) <+ TYPE * (cbs - csub);
	  end

	// Process drain/source resistance
	if ( drainConductance > 0.0 )
	  I(drain, drainp) <+ drainConductance * V(drain, drainp);
	else
	  V(drain, drainp) <+ 0.0;

	if ( sourceConductance > 0.0 )
	  I(source, sourcep) <+ sourceConductance * V(source, sourcep);
	else
	  V(source, sourcep) <+ 0.0;
	
	// Charge current including overlap and depletion capacitance contribution
	cqgate = TYPE * ddt(Qgate);
	cqdrn  = TYPE * ddt(Qdrn);
	cqbulk = TYPE * ddt(Qbulk);
	
	I(gate)    <+ cqgate;
	I(drainp)  <+ cqdrn;
	I(bulk)    <+ cqbulk;
	I(sourcep) <+ -(cqgate + cqdrn + cqbulk);
	
`ifdef NQSMOD
	if (NQSMOD)
	  begin
             I(gate)    <+ TYPE * (-1)   * qdef * gtau;
             I(drainp)  <+ TYPE * dxpart * qdef * gtau;
             I(sourcep) <+ TYPE * sxpart * qdef * gtau;
	     
             I(q) <+ -TYPE * ( cqdef - cqcheq );
             I(q) <+ -V(q) * gtau;
	  end
	else
	  begin
	     I(q) <+ GMIN * qdef;
	  end
`endif
	
     end // analog begin

endmodule
