diff --git a/examples/cider/surfmob/Tbicmpd1.cir b/examples/cider/surfmob/Tbicmpd1.cir
new file mode 100644
index 000000000..2da793fbf
--- /dev/null
+++ b/examples/cider/surfmob/Tbicmpd1.cir
@@ -0,0 +1,111 @@
+BiCMOS Pulldown Circuit
+
+VSS 2 0 0v
+
+VIN 3 2 0v (PULSE 0.0v 4.2v 0ns 1ns 1ns 9ns 20ns)
+
+M1  8 3 5 11 M_NMOS_1 W=4u L=1u
+VD  4 8 0v
+VBK 11 2 0v
+
+Q1  10 7 9   M_NPN AREA=8
+VC  4 10 0v
+VB  5 7 0v
+VE  9 2 0v
+
+CL  4 6 1pF
+VL  6 2 0v
+
+.IC V(10)=5.0v V(7)=0.0v 
+*.TRAN 0.1ns 5ns 0ns 0.1ns
+.TRAN 0.1ns 0.3ns 0ns 0.1ns
+.PLOT TRAN I(VIN)
+
+*.include bicmos.lib
+.MODEL M_NPN nbjt level=2
++ title TWO-DIMENSIONAL NUMERICAL POLYSILICON EMITTER BIPOLAR TRANSISTOR
++ * Since, we are only simulating half of a device, we double the unit width
++ * 1.0 um emitter length
++ options defw=2.0u
++ 
++ *x.mesh w=2.5 n=5
++ x.mesh w=2.0  h.e=0.05 h.m=0.2 r=1.5
++ x.mesh w=0.5  h.s=0.05 h.m=0.1 r=1.5
++ 
++ y.mesh l=-0.2 n=1
++ y.mesh l= 0.0 n=5
++ y.mesh w=0.10 h.e=0.002 h.m=0.01  r=1.5
++ y.mesh w=0.15 h.s=0.002 h.m=0.01  r=1.5
++ y.mesh w=0.35 h.s=0.01  h.m=0.2   r=1.5
++ y.mesh w=0.40 h.e=0.05  h.m=0.2   r=1.5
++ y.mesh w=0.30 h.s=0.05  h.m=0.1   r=1.5
++
++ domain num=1 material=1 x.l=2.0 y.h=0.0
++ domain num=2 material=2 x.h=2.0 y.h=0.0
++ domain num=3 material=3 y.l=0.0
++ material num=1 polysilicon
++ material num=2 oxide
++ material num=3 silicon
++
++ elec num=1 x.l=0.0  x.h=0.0  y.l=1.1  y.h=1.3
++ elec num=2 x.l=0.0  x.h=0.5  y.l=0.0  y.h=0.0
++ elec num=3 x.l=2.0  x.h=3.0  y.l=-0.2 y.h=-0.2
++
++ doping gauss n.type conc=3e20 x.l=2.0 x.h=3.0 y.l=-0.2 y.h=0.0
++ + char.l=0.047 lat.rotate
++ doping gauss p.type conc=1e19 x.l=0.0 x.h=5.0 y.l=-0.2 y.h=0.0
++ + char.l=0.094 lat.rotate
++ doping unif  n.type conc=1e16 x.l=0.0 x.h=5.0 y.l=0.0 y.h=1.3
++ doping gauss n.type conc=5e19 x.l=0.0 x.h=5.0 y.l=1.3 y.h=1.3
++ + char.l=0.100 lat.rotate
++
++ method ac=direct itlim=10
++ models bgn srh auger conctau concmob fieldmob
+
+.MODEL M_NMOS_1 numos
++ title 1.0um NMOS Device
++ 
++ x.mesh w=0.9 h.e=0.020 h.m=0.2 r=2.0
++ x.mesh w=0.2 h.e=0.005 h.m=0.02 r=2.0
++ x.mesh w=0.4 h.s=0.005 h.m=0.1 r=2.0
++ x.mesh w=0.4 h.e=0.005 h.m=0.1 r=2.0
++ x.mesh w=0.2 h.e=0.005 h.m=0.02 r=2.0
++ x.mesh w=0.9 h.s=0.020 h.m=0.2 r=2.0
++
++ y.mesh l=-.0200 n=1
++ y.mesh l=0.0 n=6
++ y.mesh w=0.15 h.s=0.0001 h.max=.02 r=2.0
++ y.mesh w=0.45 h.s=0.02 h.max=0.2 r=2.0
++ y.mesh w=1.40 h.s=0.20 h.max=0.4 r=2.0
++
++ region num=1 material=1 y.h=0.0
++ region num=2 material=2 y.l=0.0
++ interface dom=2 nei=1 x.l=1.0 x.h=2.0 layer.width=0.0
++ material num=1 oxide
++ material num=2 silicon
++
++ elec num=1 x.l=2.5 x.h=3.1	y.l=0.0 y.h=0.0
++ elec num=2 x.l=1.0  x.h=2.0	iy.l=1  iy.h=1
++ elec num=3 x.l=-0.1 x.h=0.5	y.l=0.0 y.h=0.0
++ elec num=4 x.l=-0.1 x.h=3.1   y.l=2.0 y.h=2.0
++
++ doping gauss p.type conc=1.0e17 x.l=-0.1 x.h=3.1 y.l=0.0
++ + char.l=0.30
++ doping unif p.type conc=5.0e15 x.l=-0.1 x.h=3.1 y.l=0.0 y.h=2.1
++ doping gauss n.type conc=4e17  x.l=-0.1 x.h=1.0 y.l=0.0 y.h=0.0
++ + char.l=0.16 lat.rotate ratio=0.65
++ doping gauss n.type conc=1e20  x.l=-0.1 x.h=0.95 y.l=0.0 y.h=0.08
++ + char.l=0.03 lat.rotate ratio=0.65
++ doping gauss n.type conc=4e17  x.l=2.0 x.h=3.1 y.l=0.0 y.h=0.0
++ + char.l=0.16 lat.rotate ratio=0.65
++ doping gauss n.type conc=1e20  x.l=2.05 x.h=3.1 y.l=0.0 y.h=0.08
++ + char.l=0.03 lat.rotate ratio=0.65
++
++ contact num=2 workf=4.10
++ models concmob fieldmob surfmob srh auger conctau bgn ^aval
++ method ac=direct itlim=10 onec
+
+
+* .OPTIONS ACCT BYPASS=1 filetype=ascii
+* .OPTIONS filetype=ascii
+.END
diff --git a/examples/cider/surfmob/Tbicmpd1xy.cir b/examples/cider/surfmob/Tbicmpd1xy.cir
new file mode 100644
index 000000000..12fc7541b
--- /dev/null
+++ b/examples/cider/surfmob/Tbicmpd1xy.cir
@@ -0,0 +1,111 @@
+BiCMOS Pulldown Circuit
+
+VSS 2 0 0v
+
+VIN 3 2 0v (PULSE 0.0v 4.2v 0ns 1ns 1ns 9ns 20ns)
+
+M1  8 3 5 11 M_NMOS_1 W=4u L=1u
+VD  4 8 0v
+VBK 11 2 0v
+
+Q1  10 7 9   M_NPN AREA=8
+VC  4 10 0v
+VB  5 7 0v
+VE  9 2 0v
+
+CL  4 6 1pF
+VL  6 2 0v
+
+.IC V(10)=5.0v V(7)=0.0v 
+*.TRAN 0.1ns 5ns 0ns 0.1ns
+.TRAN 0.1ns 0.3ns 0ns 0.1ns
+.PLOT TRAN I(VIN)
+
+*.include bicmos.lib
+.MODEL M_NPN nbjt level=2
++ title TWO-DIMENSIONAL NUMERICAL POLYSILICON EMITTER BIPOLAR TRANSISTOR
++ * Since, we are only simulating half of a device, we double the unit width
++ * 1.0 um emitter length
++ options defw=2.0u
++ 
++ *y.mesh w=2.5 n=5
++ y.mesh w=2.0  h.e=0.05 h.m=0.2 r=1.5
++ y.mesh w=0.5  h.s=0.05 h.m=0.1 r=1.5
++ 
++ x.mesh l=-0.2 n=1
++ x.mesh l= 0.0 n=5
++ x.mesh w=0.10 h.e=0.002 h.m=0.01  r=1.5
++ x.mesh w=0.15 h.s=0.002 h.m=0.01  r=1.5
++ x.mesh w=0.35 h.s=0.01  h.m=0.2   r=1.5
++ x.mesh w=0.40 h.e=0.05  h.m=0.2   r=1.5
++ x.mesh w=0.30 h.s=0.05  h.m=0.1   r=1.5
++
++ domain num=1 material=1 y.l=2.0 x.h=0.0
++ domain num=2 material=2 y.h=2.0 x.h=0.0
++ domain num=3 material=3 x.l=0.0
++ material num=1 polysilicon
++ material num=2 oxide
++ material num=3 silicon
++
++ elec num=1 y.l=0.0  y.h=0.0  x.l=1.1  x.h=1.3
++ elec num=2 y.l=0.0  y.h=0.5  x.l=0.0  x.h=0.0
++ elec num=3 y.l=2.0  y.h=3.0  x.l=-0.2 x.h=-0.2
++
++ doping gauss n.type conc=3e20 y.l=2.0 y.h=3.0 x.l=-0.2 x.h=0.0
++ + char.l=0.047 lat.rotate
++ doping gauss p.type conc=1e19 y.l=0.0 y.h=5.0 x.l=-0.2 x.h=0.0
++ + char.l=0.094 lat.rotate
++ doping unif  n.type conc=1e16 y.l=0.0 y.h=5.0 x.l=0.0 x.h=1.3
++ doping gauss n.type conc=5e19 y.l=0.0 y.h=5.0 x.l=1.3 x.h=1.3
++ + char.l=0.100 lat.rotate
++
++ method ac=direct itlim=10
++ models bgn srh auger conctau concmob fieldmob
+
+.MODEL M_NMOS_1 numos
++ title 1.0um NMOS Device
++ 
++ y.mesh w=0.9 h.e=0.020 h.m=0.2 r=2.0
++ y.mesh w=0.2 h.e=0.005 h.m=0.02 r=2.0
++ y.mesh w=0.4 h.s=0.005 h.m=0.1 r=2.0
++ y.mesh w=0.4 h.e=0.005 h.m=0.1 r=2.0
++ y.mesh w=0.2 h.e=0.005 h.m=0.02 r=2.0
++ y.mesh w=0.9 h.s=0.020 h.m=0.2 r=2.0
++
++ x.mesh l=-.0200 n=1
++ x.mesh l=0.0 n=6
++ x.mesh w=0.15 h.s=0.0001 h.max=.02 r=2.0
++ x.mesh w=0.45 h.s=0.02 h.max=0.2 r=2.0
++ x.mesh w=1.40 h.s=0.20 h.max=0.4 r=2.0
++
++ region num=1 material=1 x.h=0.0
++ region num=2 material=2 x.l=0.0
++ interface dom=2 nei=1 y.l=1.0 y.h=2.0 layer.width=0.0
++ material num=1 oxide
++ material num=2 silicon
++
++ elec num=1 y.l=2.5 y.h=3.1	x.l=0.0 x.h=0.0
++ elec num=2 y.l=1.0  y.h=2.0	ix.l=1  ix.h=1
++ elec num=3 y.l=-0.1 y.h=0.5	x.l=0.0 x.h=0.0
++ elec num=4 y.l=-0.1 y.h=3.1   x.l=2.0 x.h=2.0
++
++ doping gauss p.type conc=1.0e17 y.l=-0.1 y.h=3.1 x.l=0.0
++ + char.l=0.30
++ doping unif p.type conc=5.0e15 y.l=-0.1 y.h=3.1 x.l=0.0 x.h=2.1
++ doping gauss n.type conc=4e17  y.l=-0.1 y.h=1.0 x.l=0.0 x.h=0.0
++ + char.l=0.16 lat.rotate ratio=0.65
++ doping gauss n.type conc=1e20  y.l=-0.1 y.h=0.95 x.l=0.0 x.h=0.08
++ + char.l=0.03 lat.rotate ratio=0.65
++ doping gauss n.type conc=4e17  y.l=2.0 y.h=3.1 x.l=0.0 x.h=0.0
++ + char.l=0.16 lat.rotate ratio=0.65
++ doping gauss n.type conc=1e20  y.l=2.05 y.h=3.1 x.l=0.0 x.h=0.08
++ + char.l=0.03 lat.rotate ratio=0.65
++
++ contact num=2 workf=4.10
++ models concmob fieldmob surfmob srh auger conctau bgn ^aval
++ method ac=direct itlim=10 onec
+
+
+* .OPTIONS ACCT BYPASS=1 filetype=ascii
+* .OPTIONS filetype=ascii
+.END
diff --git a/examples/cider/surfmob/Tbicmpu1.cir b/examples/cider/surfmob/Tbicmpu1.cir
new file mode 100644
index 000000000..d1e381786
--- /dev/null
+++ b/examples/cider/surfmob/Tbicmpu1.cir
@@ -0,0 +1,105 @@
+BICMOS INVERTER PULLUP CIRCUIT
+
+VDD 1 0 5.0V
+VSS 2 0 0.0V
+
+VIN 3 0 0.75V
+
+VC  1 11 0.0V
+VB  5 15 0.0V
+
+Q1  11 15 4 M_NPNS   AREA=8
+M1  5 3 1 1 M_PMOS_1 W=10U L=1U
+
+CL  4 0 5.0PF
+
+.IC V(4)=0.75V V(5)=0.0V
+
+.MODEL M_NPNS nbjt level=2
++ title TWO-DIMENSIONAL NUMERICAL POLYSILICON EMITTER BIPOLAR TRANSISTOR
++ * Since half the device is simulated, double the unit width to get
++ * 1.0 um emitter.  Use a small mesh for this model.
++ options defw=2.0u
++ 
++ x.mesh w=2.0 h.e=0.02 h.m=0.5 r=2.0
++ x.mesh w=0.5 h.s=0.02 h.m=0.2 r=2.0
++ 
++ y.mesh l=-0.2 n=1
++ y.mesh l= 0.0 n=5
++ y.mesh w=0.10 h.e=0.004 h.m=0.05  r=2.5
++ y.mesh w=0.15 h.s=0.004 h.m=0.02  r=2.5
++ y.mesh w=1.05 h.s=0.02  h.m=0.1   r=2.5
++
++ domain num=1 material=1 x.l=2.0 y.h=0.0
++ domain num=2 material=2 x.h=2.0 y.h=0.0
++ domain num=3 material=3 y.l=0.0
++ material num=1 polysilicon
++ material num=2 oxide
++ material num=3 silicon
++
++ elec num=1 x.l=0.0  x.h=0.0  y.l=1.1  y.h=1.3
++ elec num=2 x.l=0.0  x.h=0.5  y.l=0.0  y.h=0.0
++ elec num=3 x.l=2.0  x.h=3.0  y.l=-0.2 y.h=-0.2
++
++ doping gauss n.type conc=3e20 x.l=2.0 x.h=3.0 y.l=-0.2 y.h=0.0
++ + char.l=0.047 lat.rotate
++ doping gauss p.type conc=5e18 x.l=0.0 x.h=5.0 y.l=-0.2 y.h=0.0
++ + char.l=0.100 lat.rotate
++ doping gauss p.type conc=1e20 x.l=0.0 x.h=0.5 y.l=-0.2 y.h=0.0
++ + char.l=0.100 lat.rotate ratio=0.7
++ doping unif  n.type conc=1e16 x.l=0.0 x.h=5.0 y.l=0.0 y.h=1.3
++ doping gauss n.type conc=5e19 x.l=0.0 x.h=5.0 y.l=1.3 y.h=1.3
++ + char.l=0.100 lat.rotate
++
++ method ac=direct itlim=10
++ models bgn srh auger conctau concmob fieldmob
+
+.MODEL M_PMOS_1 numos
++ 
++ x.mesh w=0.9 h.e=0.020 h.m=0.2 r=2.0
++ x.mesh w=0.2 h.e=0.005 h.m=0.02 r=2.0
++ x.mesh w=0.4 h.s=0.005 h.m=0.1 r=2.0
++ x.mesh w=0.4 h.e=0.005 h.m=0.1 r=2.0
++ x.mesh w=0.2 h.e=0.005 h.m=0.02 r=2.0
++ x.mesh w=0.9 h.s=0.020 h.m=0.2 r=2.0
++
++ y.mesh l=-.0200 n=1
++ y.mesh l=0.0 n=6
++ y.mesh w=0.15 h.s=0.0001 h.max=.02 r=2.0
++ y.mesh w=0.45 h.s=0.02 h.max=0.2 r=2.0
++ y.mesh w=1.40 h.s=0.20 h.max=0.4 r=2.0
++
++ region num=1 material=1 y.h=0.0
++ region num=2 material=2 y.l=0.0
++ interface dom=2 nei=1 x.l=1 x.h=2 layer.width=0.0
++ material num=1 oxide
++ material num=2 silicon
++
++ elec num=1 x.l=2.5  x.h=3.1  y.l=0.0 y.h=0.0
++ elec num=2 x.l=1 x.h=2  iy.l=1  iy.h=1
++ elec num=3 x.l=-0.1  x.h=0.5  y.l=0.0 y.h=0.0
++ elec num=4 x.l=-0.1  x.h=3.1 y.l=2.0 y.h=2.0
++
++ doping gauss n.type conc=1.0e17 x.l=-0.1 x.h=3.1 y.l=0.0
++ + char.l=0.30
++ doping unif n.type conc=5.0e15 x.l=-0.1 x.h=3.1 y.l=0.0 y.h=2.1
++ doping gauss p.type conc=4e17  x.l=-0.1 x.h=1 y.l=0.0 y.h=0.0
++ + char.l=0.16 lat.rotate ratio=0.65
++ doping gauss p.type conc=1e20  x.l=-0.1 x.h=0.95 y.l=0.0 y.h=0.08
++ + char.l=0.03 lat.rotate ratio=0.65
++ doping gauss p.type conc=4e17  x.l=2 x.h=3.1 y.l=0.0 y.h=0.0
++ + char.l=0.16 lat.rotate ratio=0.65
++ doping gauss p.type conc=1e20  x.l=2.05 x.h=3.1 y.l=0.0 y.h=0.08
++ + char.l=0.03 lat.rotate ratio=0.65
++
++ contact num=2 workf=5.29
++ models concmob surfmob transmob fieldmob srh auger conctau bgn 
++ method ac=direct itlim=10 onec
+
+
+.TRAN 0.5NS 3.0NS
+.PRINT TRAN V(3) V(4)
+.PLOT TRAN V(3) V(4)
+
+* .OPTION ACCT BYPASS=1
+.END
diff --git a/examples/cider/surfmob/Tbicmpu1xy.cir b/examples/cider/surfmob/Tbicmpu1xy.cir
new file mode 100644
index 000000000..229ea31db
--- /dev/null
+++ b/examples/cider/surfmob/Tbicmpu1xy.cir
@@ -0,0 +1,105 @@
+BICMOS INVERTER PULLUP CIRCUIT
+
+VDD 1 0 5.0V
+VSS 2 0 0.0V
+
+VIN 3 0 0.75V
+
+VC  1 11 0.0V
+VB  5 15 0.0V
+
+Q1  11 15 4 M_NPNS   AREA=8
+M1  5 3 1 1 M_PMOS_1 W=10U L=1U
+
+CL  4 0 5.0PF
+
+.IC V(4)=0.75V V(5)=0.0V
+
+.MODEL M_NPNS nbjt level=2
++ title TWO-DIMENSIONAL NUMERICAL POLYSILICON EMITTER BIPOLAR TRANSISTOR
++ * Since half the device is simulated, double the unit width to get
++ * 1.0 um emitter.  Use a small mesh for this model.
++ options defw=2.0u
++ 
++ x.mesh w=2.0 h.e=0.02 h.m=0.5 r=2.0
++ x.mesh w=0.5 h.s=0.02 h.m=0.2 r=2.0
++ 
++ y.mesh l=-0.2 n=1
++ y.mesh l= 0.0 n=5
++ y.mesh w=0.10 h.e=0.004 h.m=0.05  r=2.5
++ y.mesh w=0.15 h.s=0.004 h.m=0.02  r=2.5
++ y.mesh w=1.05 h.s=0.02  h.m=0.1   r=2.5
++
++ domain num=1 material=1 x.l=2.0 y.h=0.0
++ domain num=2 material=2 x.h=2.0 y.h=0.0
++ domain num=3 material=3 y.l=0.0
++ material num=1 polysilicon
++ material num=2 oxide
++ material num=3 silicon
++
++ elec num=1 x.l=0.0  x.h=0.0  y.l=1.1  y.h=1.3
++ elec num=2 x.l=0.0  x.h=0.5  y.l=0.0  y.h=0.0
++ elec num=3 x.l=2.0  x.h=3.0  y.l=-0.2 y.h=-0.2
++
++ doping gauss n.type conc=3e20 x.l=2.0 x.h=3.0 y.l=-0.2 y.h=0.0
++ + char.l=0.047 lat.rotate
++ doping gauss p.type conc=5e18 x.l=0.0 x.h=5.0 y.l=-0.2 y.h=0.0
++ + char.l=0.100 lat.rotate
++ doping gauss p.type conc=1e20 x.l=0.0 x.h=0.5 y.l=-0.2 y.h=0.0
++ + char.l=0.100 lat.rotate ratio=0.7
++ doping unif  n.type conc=1e16 x.l=0.0 x.h=5.0 y.l=0.0 y.h=1.3
++ doping gauss n.type conc=5e19 x.l=0.0 x.h=5.0 y.l=1.3 y.h=1.3
++ + char.l=0.100 lat.rotate
++
++ method ac=direct itlim=10
++ models bgn srh auger conctau concmob fieldmob
+
+.MODEL M_PMOS_1 numos
++ 
++ y.mesh w=0.9 h.e=0.020 h.m=0.2 r=2.0
++ y.mesh w=0.2 h.e=0.005 h.m=0.02 r=2.0
++ y.mesh w=0.4 h.s=0.005 h.m=0.1 r=2.0
++ y.mesh w=0.4 h.e=0.005 h.m=0.1 r=2.0
++ y.mesh w=0.2 h.e=0.005 h.m=0.02 r=2.0
++ y.mesh w=0.9 h.s=0.020 h.m=0.2 r=2.0
++
++ x.mesh l=-.0200 n=1
++ x.mesh l=0.0 n=6
++ x.mesh w=0.15 h.s=0.0001 h.max=.02 r=2.0
++ x.mesh w=0.45 h.s=0.02 h.max=0.2 r=2.0
++ x.mesh w=1.40 h.s=0.20 h.max=0.4 r=2.0
++
++ region num=1 material=1 x.h=0.0
++ region num=2 material=2 x.l=0.0
++ interface dom=2 nei=1 y.l=1 y.h=2 layer.width=0.0
++ material num=1 oxide
++ material num=2 silicon
++
++ elec num=1 y.l=2.5   y.h=3.1  x.l=0.0 x.h=0.0
++ elec num=2 y.l=1     y.h=2    ix.l=1  ix.h=1
++ elec num=3 y.l=-0.1  y.h=0.5  x.l=0.0 x.h=0.0
++ elec num=4 y.l=-0.1  y.h=3.1  x.l=2.0 x.h=2.0
++
++ doping gauss n.type conc=1.0e17 y.l=-0.1 y.h=3.1 x.l=0.0
++ + char.l=0.30
++ doping unif n.type conc=5.0e15 y.l=-0.1 y.h=3.1 x.l=0.0 x.h=2.1
++ doping gauss p.type conc=4e17  y.l=-0.1 y.h=1 x.l=0.0 x.h=0.0
++ + char.l=0.16 lat.rotate ratio=0.65
++ doping gauss p.type conc=1e20  y.l=-0.1 y.h=0.95 x.l=0.0 x.h=0.08
++ + char.l=0.03 lat.rotate ratio=0.65
++ doping gauss p.type conc=4e17  y.l=2 y.h=3.1 x.l=0.0 x.h=0.0
++ + char.l=0.16 lat.rotate ratio=0.65
++ doping gauss p.type conc=1e20  y.l=2.05 y.h=3.1 x.l=0.0 x.h=0.08
++ + char.l=0.03 lat.rotate ratio=0.65
++
++ contact num=2 workf=5.29
++ models concmob surfmob transmob fieldmob srh auger conctau bgn 
++ method ac=direct itlim=10 onec
+
+
+.TRAN 0.5NS 3.0NS
+.PRINT TRAN V(3) V(4)
+.PLOT TRAN V(3) V(4)
+
+* .OPTION ACCT BYPASS=1
+.END
diff --git a/examples/cider/surfmob/cd4007n_200402idvg.net b/examples/cider/surfmob/cd4007n_200402idvg.net
new file mode 100644
index 000000000..6d8c45e1a
--- /dev/null
+++ b/examples/cider/surfmob/cd4007n_200402idvg.net
@@ -0,0 +1,159 @@
+$ CD4007, version: cd4007n_200402idvg 
+$ Generated using Matlab automation script. 
+
+
+$ BEGIN CIDER 
+
+.model CD4007UBENUMOSN numos
+$ -- x Grid 
+$ ---- GOX 
++ x.mesh (location = -64.0000e-3) (node = 1)
++ x.mesh (location = -48.0000e-3) (node = 2)
++ x.mesh (location = -32.0000e-3) (node = 3)
++ x.mesh (location = -16.0000e-3) (node = 4)
++ x.mesh (location = 0.0000e-3) (node = 5)
+$ ---- Channel and S/D 
++ x.mesh (location = 20.0000e-3) (node = 6)
++ x.mesh (location = 48.4511e-3) (node = 7)
++ x.mesh (location = 86.9893e-3) (node = 8)
++ x.mesh (location = 139.1911e-3) (node = 9)
++ x.mesh (location = 209.9008e-3) (node = 10)
++ x.mesh (location = 305.6803e-3) (node = 11)
++ x.mesh (location = 435.4180e-3) (node = 12)
++ x.mesh (location = 611.1537e-3) (node = 13)
++ x.mesh (location = 849.1956e-3) (node = 14)
++ x.mesh (location = 1171.6341e-3) (node = 15)
++ x.mesh (location = 1608.3918e-3) (node = 16)
++ x.mesh (location = 2200.0000e-3) (node = 17)
++ x.mesh (location = 3000.0000e-3) (node = 18)
+$ ---- Mid substrate 
++ x.mesh (location = 3800.0000e-3) (node = 19)
++ x.mesh (location = 4911.8170e-3) (node = 20)
++ x.mesh (location = 6290.4701e-3) (node = 21)
++ x.mesh (location = 8000.0000e-3) (node = 22)
++ x.mesh (location = 10000.0000e-3) (node = 23)
++ 
+$ -- y Grid 
+$ ---- Body well 
++ y.mesh (location = -35000.0000e-3) (node = 1)
++ y.mesh (location = -30000.0000e-3) (node = 2)
++ y.mesh (location = -25000.0000e-3) (node = 3)
++ y.mesh (location = -20000.0000e-3) (node = 4)
+$ ---- Body-Source 
++ y.mesh (location = -17500.0000e-3) (node = 5)
++ y.mesh (location = -15000.0000e-3) (node = 6)
++ y.mesh (location = -12500.0000e-3) (node = 7)
++ y.mesh (location = -10000.0000e-3) (node = 8)
+$ ---- Mid source 
++ y.mesh (location = -7500.0000e-3) (node = 9)
++ y.mesh (location = -4980.0000e-3) (node = 10)
++ y.mesh (location = -3000.0000e-3) (node = 11)
++ y.mesh (location = -2000.0000e-3) (node = 12)
+$ ---- Near source 
++ y.mesh (location = -1000.0000e-3) (node = 13)
++ y.mesh (location = -200.0000e-3) (node = 14)
++ y.mesh (location = 0.0000e-3) (node = 15)
+$ ---- Channel near source 
++ y.mesh (location = 200.0000e-3) (node = 16)
++ y.mesh (location = 442.1759e-3) (node = 17)
++ y.mesh (location = 732.7869e-3) (node = 18)
++ y.mesh (location = 1081.5201e-3) (node = 19)
++ y.mesh (location = 1500.0000e-3) (node = 20)
++ y.mesh (location = 2000.0000e-3) (node = 21)
+$ ---- Channel center 
++ y.mesh (location = 2500.0000e-3) (node = 22)
++ y.mesh (location = 3000.0000e-3) (node = 23)
++ y.mesh (location = 3500.0000e-3) (node = 24)
++ y.mesh (location = 4000.0000e-3) (node = 25)
++ y.mesh (location = 4500.0000e-3) (node = 26)
++ y.mesh (location = 5000.0000e-3) (node = 27)
+$ ---- Channel near drain 
++ y.mesh (location = 5500.0000e-3) (node = 28)
++ y.mesh (location = 5961.8362e-3) (node = 29)
++ y.mesh (location = 6376.3046e-3) (node = 30)
++ y.mesh (location = 6748.2634e-3) (node = 31)
++ y.mesh (location = 7082.0725e-3) (node = 32)
++ y.mesh (location = 7381.6449e-3) (node = 33)
++ y.mesh (location = 7650.4919e-3) (node = 34)
++ y.mesh (location = 7891.7648e-3) (node = 35)
++ y.mesh (location = 8108.2918e-3) (node = 36)
++ y.mesh (location = 8302.6109e-3) (node = 37)
++ y.mesh (location = 8476.9998e-3) (node = 38)
++ y.mesh (location = 8633.5027e-3) (node = 39)
++ y.mesh (location = 8773.9540e-3) (node = 40)
++ y.mesh (location = 8900.0000e-3) (node = 41)
++ y.mesh (location = 9000.0000e-3) (node = 42)
+$ ---- Near drain 
++ y.mesh (location = 9200.0000e-3) (node = 43)
++ y.mesh (location = 9612.8440e-3) (node = 44)
++ y.mesh (location = 10190.8257e-3) (node = 45)
++ y.mesh (location = 11000.0000e-3) (node = 46)
++ y.mesh (location = 12000.0000e-3) (node = 47)
+$ ---- Mid drain 
++ y.mesh (location = 13000.0000e-3) (node = 48)
++ y.mesh (location = 15500.0000e-3) (node = 49)
++ y.mesh (location = 18000.0000e-3) (node = 50)
++ 
+$ -- Regions 
+$ ---- substrate 
++ region (num=1 material=1) x.l = 0.0000e-3
+$ ---- GOX 
++ region (num=2 material=2) x.h = 0.0000e-3
++ 
+$ -- Materials 
++ material num=1 silicon
++ mobility material=1 concmod=sg fieldmod=sg
++ Mobility Material=1 elec major MuMax=1200.0 MuMin=200.0 Vsat=3.2e7
++ Mobility Material=1 elec MuS=600.0 EC.A=3.0e5
++ mobility material=1 elec minor MuMax=1200.0 MuMin=200.0
++ mobility material=1 hole major
++ mobility material=1 hole minor
++ material num=2 oxide
++ 
+$ -- Contacts and workfunction 
+$ ---- Drain 
++ elec num=1 (x.l = 0.0000e-3 x.h = 0.0000e-3) (y.l = 12000.0000e-3 y.h = 16000.0000e-3)
+$ ---- Gate contact 
++ elec num=2 (x.l = -64.0000e-3 x.h = -64.0000e-3) (y.l = 0.0000e-3 y.h = 9000.0000e-3)
+$ ---- Source 
++ elec num=3 (x.l = 0.0000e-3 x.h = 0.0000e-3) (y.l = -6000.0000e-3 y.h = -3000.0000e-3)
+$ ---- Body 
++ elec num=4 (x.l = 0.0000e-3 x.h = 0.0000e-3) (y.l = -30000.0000e-3 y.h = -25000.0000e-3)
++ contact num=1 workf=4.0500
++ contact num=2 workf=4.0500
++ contact num=3 workf=4.0500
++ contact num=4 workf=5.1700
++ 
+$ -- Doping profiles 
+$ ---- p subs 
++ doping Domains = 1 Uniform P.Type Conc = 1e16
+$ ---- Threshold Adjustment Diffuse 
++ doping Domains = 1 Gaussian Lat.Erfc X.Axis P.Type x.l=0.0 x.h=0.5 y.l=-0 y.h=10 Peak.Conc=2e16 Location=0.1 Char.Length=0.8 Ratio.Lat=0.8
+$ ---- Body Diffuse 
++ doping Domains = 1 Gaussian Lat.Erfc X.Axis P.Type x.l=0.0 x.h=3.0 y.l=-30 y.h=-20 Peak.Conc=2e20 Location=0 Char.Length=0.6 Ratio.Lat=0.8
+$ ---- Source Diffuse 
++ doping Domains = 1 Gaussian Lat.Erfc X.Axis N.Type x.l=0.0 x.h=3.0 y.l=-10 y.h=0 Peak.Conc=2e20 Location=0 Char.Length=0.6 Ratio.Lat=0.8
+$ ---- Drain Diffuse 
++ doping Domains = 1 Gaussian Lat.Erfc X.Axis N.Type x.l=0.0 x.h=3.0 y.l=9.0 y.h=18.0 Peak.Conc=2e20 Location=0 Char.Length=0.6 Ratio.Lat=0.8
++ 
++ models concmob fieldmob
++ method ac=direct ^onec
++ 
++ interface domain=1 neighbor=2 y.l=0 y.h=8.0 layer.width=0.5
++ models surfmob transmob ^matchmob
+
+$ Total grid points X23 x Y50 = 1150 
+
+$ END OF CIDER 
+
+
+$ PSPICE Schematics Netlist 
+VId VPP VDD 0
+MN1 VDD VGG 0 0 CD4007UBENUMOSN w=298u
+* .save i(VId) *
+VVDS VPP 0   DC  0
+VVGS VGG 0   DC  0
+.dc VVGS 0.00 5.00 0.25 VVDS 0.05 0.05 0.50
+$ VGS x VDS sweep points: 51 x 1 = 51  
+.print dc VPP VGG VPP
+.end
diff --git a/examples/cider/surfmob/cdxy.net b/examples/cider/surfmob/cdxy.net
new file mode 100644
index 000000000..f420c0c80
--- /dev/null
+++ b/examples/cider/surfmob/cdxy.net
@@ -0,0 +1,159 @@
+$ CD4007, version: cd4007n_200402idvg 
+$ Generated using Matlab automation script. 
+
+
+$ BEGIN CIDER 
+
+.model CD4007UBENUMOSN numos
+$ -- y Grid 
+$ ---- GOX 
++ y.mesh (location = -64.0000e-3) (node = 1)
++ y.mesh (location = -48.0000e-3) (node = 2)
++ y.mesh (location = -32.0000e-3) (node = 3)
++ y.mesh (location = -16.0000e-3) (node = 4)
++ y.mesh (location = 0.0000e-3) (node = 5)
+$ ---- Channel and S/D 
++ y.mesh (location = 20.0000e-3) (node = 6)
++ y.mesh (location = 48.4511e-3) (node = 7)
++ y.mesh (location = 86.9893e-3) (node = 8)
++ y.mesh (location = 139.1911e-3) (node = 9)
++ y.mesh (location = 209.9008e-3) (node = 10)
++ y.mesh (location = 305.6803e-3) (node = 11)
++ y.mesh (location = 435.4180e-3) (node = 12)
++ y.mesh (location = 611.1537e-3) (node = 13)
++ y.mesh (location = 849.1956e-3) (node = 14)
++ y.mesh (location = 1171.6341e-3) (node = 15)
++ y.mesh (location = 1608.3918e-3) (node = 16)
++ y.mesh (location = 2200.0000e-3) (node = 17)
++ y.mesh (location = 3000.0000e-3) (node = 18)
+$ ---- Mid substrate 
++ y.mesh (location = 3800.0000e-3) (node = 19)
++ y.mesh (location = 4911.8170e-3) (node = 20)
++ y.mesh (location = 6290.4701e-3) (node = 21)
++ y.mesh (location = 8000.0000e-3) (node = 22)
++ y.mesh (location = 10000.0000e-3) (node = 23)
++ 
+$ -- x Grid 
+$ ---- Body well 
++ x.mesh (location = -35000.0000e-3) (node = 1)
++ x.mesh (location = -30000.0000e-3) (node = 2)
++ x.mesh (location = -25000.0000e-3) (node = 3)
++ x.mesh (location = -20000.0000e-3) (node = 4)
+$ ---- Body-Source 
++ x.mesh (location = -17500.0000e-3) (node = 5)
++ x.mesh (location = -15000.0000e-3) (node = 6)
++ x.mesh (location = -12500.0000e-3) (node = 7)
++ x.mesh (location = -10000.0000e-3) (node = 8)
+$ ---- Mid source 
++ x.mesh (location = -7500.0000e-3) (node = 9)
++ x.mesh (location = -4980.0000e-3) (node = 10)
++ x.mesh (location = -3000.0000e-3) (node = 11)
++ x.mesh (location = -2000.0000e-3) (node = 12)
+$ ---- Near source 
++ x.mesh (location = -1000.0000e-3) (node = 13)
++ x.mesh (location = -200.0000e-3) (node = 14)
++ x.mesh (location = 0.0000e-3) (node = 15)
+$ ---- Channel near source 
++ x.mesh (location = 200.0000e-3) (node = 16)
++ x.mesh (location = 442.1759e-3) (node = 17)
++ x.mesh (location = 732.7869e-3) (node = 18)
++ x.mesh (location = 1081.5201e-3) (node = 19)
++ x.mesh (location = 1500.0000e-3) (node = 20)
++ x.mesh (location = 2000.0000e-3) (node = 21)
+$ ---- Channel center 
++ x.mesh (location = 2500.0000e-3) (node = 22)
++ x.mesh (location = 3000.0000e-3) (node = 23)
++ x.mesh (location = 3500.0000e-3) (node = 24)
++ x.mesh (location = 4000.0000e-3) (node = 25)
++ x.mesh (location = 4500.0000e-3) (node = 26)
++ x.mesh (location = 5000.0000e-3) (node = 27)
+$ ---- Channel near drain 
++ x.mesh (location = 5500.0000e-3) (node = 28)
++ x.mesh (location = 5961.8362e-3) (node = 29)
++ x.mesh (location = 6376.3046e-3) (node = 30)
++ x.mesh (location = 6748.2634e-3) (node = 31)
++ x.mesh (location = 7082.0725e-3) (node = 32)
++ x.mesh (location = 7381.6449e-3) (node = 33)
++ x.mesh (location = 7650.4919e-3) (node = 34)
++ x.mesh (location = 7891.7648e-3) (node = 35)
++ x.mesh (location = 8108.2918e-3) (node = 36)
++ x.mesh (location = 8302.6109e-3) (node = 37)
++ x.mesh (location = 8476.9998e-3) (node = 38)
++ x.mesh (location = 8633.5027e-3) (node = 39)
++ x.mesh (location = 8773.9540e-3) (node = 40)
++ x.mesh (location = 8900.0000e-3) (node = 41)
++ x.mesh (location = 9000.0000e-3) (node = 42)
+$ ---- Near drain 
++ x.mesh (location = 9200.0000e-3) (node = 43)
++ x.mesh (location = 9612.8440e-3) (node = 44)
++ x.mesh (location = 10190.8257e-3) (node = 45)
++ x.mesh (location = 11000.0000e-3) (node = 46)
++ x.mesh (location = 12000.0000e-3) (node = 47)
+$ ---- Mid drain 
++ x.mesh (location = 13000.0000e-3) (node = 48)
++ x.mesh (location = 15500.0000e-3) (node = 49)
++ x.mesh (location = 18000.0000e-3) (node = 50)
++ 
+$ -- Regions 
+$ ---- substrate 
++ region (num=1 material=1) y.l = 0.0000e-3
+$ ---- GOX 
++ region (num=2 material=2) y.h = 0.0000e-3
++ 
+$ -- Materials 
++ material num=1 silicon
++ mobility material=1 concmod=sg fieldmod=sg
++ Mobility Material=1 elec major MuMax=1200.0 MuMin=200.0 Vsat=3.2e7
++ Mobility Material=1 elec MuS=600.0 EC.A=3.0e5
++ mobility material=1 elec minor MuMax=1200.0 MuMin=200.0
++ mobility material=1 hole major
++ mobility material=1 hole minor
++ material num=2 oxide
++ 
+$ -- Contacts and workfunction 
+$ ---- Drain 
++ elec num=1 (y.l = 0.0000e-3 y.h = 0.0000e-3) (x.l = 12000.0000e-3 x.h = 16000.0000e-3)
+$ ---- Gate contact 
++ elec num=2 (y.l = -64.0000e-3 y.h = -64.0000e-3) (x.l = 0.0000e-3 x.h = 9000.0000e-3)
+$ ---- Source 
++ elec num=3 (y.l = 0.0000e-3 y.h = 0.0000e-3) (x.l = -6000.0000e-3 x.h = -3000.0000e-3)
+$ ---- Body 
++ elec num=4 (y.l = 0.0000e-3 y.h = 0.0000e-3) (x.l = -30000.0000e-3 x.h = -25000.0000e-3)
++ contact num=1 workf=4.0500
++ contact num=2 workf=4.0500
++ contact num=3 workf=4.0500
++ contact num=4 workf=5.1700
++ 
+$ -- Doping profiles 
+$ ---- p subs 
++ doping Domains = 1 Uniform P.Type Conc = 1e16
+$ ---- Threshold Adjustment Diffuse 
++ doping Domains = 1 Gaussian Lat.Erfc Y.Axis P.Type y.l=0.0 y.h=0.5 x.l=-0 x.h=10 Peak.Conc=2e16 Location=0.1 Char.Length=0.8 Ratio.Lat=0.8
+$ ---- Body Diffuse 
++ doping Domains = 1 Gaussian Lat.Erfc Y.Axis P.Type y.l=0.0 y.h=3.0 x.l=-30 x.h=-20 Peak.Conc=2e20 Location=0 Char.Length=0.6 Ratio.Lat=0.8
+$ ---- Source Diffuse 
++ doping Domains = 1 Gaussian Lat.Erfc Y.Axis N.Type y.l=0.0 y.h=3.0 x.l=-10 x.h=0 Peak.Conc=2e20 Location=0 Char.Length=0.6 Ratio.Lat=0.8
+$ ---- Drain Diffuse 
++ doping Domains = 1 Gaussian Lat.Erfc Y.Axis N.Type y.l=0.0 y.h=3.0 x.l=9.0 x.h=18.0 Peak.Conc=2e20 Location=0 Char.Length=0.6 Ratio.Lat=0.8
++ 
++ models concmob fieldmob
++ method ac=direct ^onec
++ 
++ interface domain=1 neighbor=2 x.l=0 x.h=8.0 layer.width=0.5
++ models surfmob transmob ^matchmob
+
+$ Total grid points X23 x Y50 = 1150 
+
+$ END OF CIDER 
+
+
+$ PSPICE Schematics Netlist 
+VId VPP VDD 0
+MN1 VDD VGG 0 0 CD4007UBENUMOSN w=298u
+* .save i(VId) *
+VVDS VPP 0   DC  0
+VVGS VGG 0   DC  0
+.dc VVGS 0.00 5.00 0.25 VVDS 0.05 0.05 0.50
+$ VGS x VDS sweep points: 51 x 1 = 51  
+.print dc VPP VGG VPP
+.end
diff --git a/examples/cider/surfmob/runtests.sh b/examples/cider/surfmob/runtests.sh
new file mode 100755
index 000000000..b137641c7
--- /dev/null
+++ b/examples/cider/surfmob/runtests.sh
@@ -0,0 +1,9 @@
+#!/bin/sh
+set -xv
+
+ngspice -b cd4007n_200402idvg.net
+ngspice -b cdxy.net
+ngspice -b Tbicmpd1.cir
+ngspice -b Tbicmpd1xy.cir
+ngspice -b Tbicmpu1.cir
+ngspice -b Tbicmpu1xy.cir
diff --git a/examples/cider/surfmob/spiceinit b/examples/cider/surfmob/spiceinit
new file mode 100644
index 000000000..d08d85362
--- /dev/null
+++ b/examples/cider/surfmob/spiceinit
@@ -0,0 +1,2 @@
+option klu
+set numdgt=3
diff --git a/src/ciderlib/oned/oneadmit.c b/src/ciderlib/oned/oneadmit.c
index 686e96a98..38c251003 100644
--- a/src/ciderlib/oned/oneadmit.c
+++ b/src/ciderlib/oned/oneadmit.c
@@ -14,6 +14,8 @@ Author:	1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
 #include "ngspice/complex.h"
 #include "ngspice/spmatrix.h"
 #include "ngspice/ifsim.h"
+#include "ngspice/cktdefs.h"
+#include "ngspice/ftedefs.h"
 
 #include "onedext.h"
 #include "oneddefs.h"
@@ -29,7 +31,17 @@ extern IFfrontEnd *SPfrontEnd;
  */
  SPcomplex yAc;
 
-
+#ifdef KLU
+static void small_signal_check(char *where)
+{
+  if (ft_curckt->ci_ckt->CKTmode == MODEAC
+   || ft_curckt->ci_ckt->CKTmode == MODEACNOISE) {
+    fprintf(stderr, "Error: CIDER %s small signal simulation is not (yet) supported with 'option klu'.\n", where);
+    fprintf(stderr, "    Use 'option sparse' instead.\n");
+    controlled_exit(1);
+  }
+}
+#endif
 
 int
 NUMDadmittance(ONEdevice *pDevice, double omega, SPcomplex *yd)
@@ -45,6 +57,11 @@ NUMDadmittance(ONEdevice *pDevice, double omega, SPcomplex *yd)
   bool SORFailed;
   double startTime;
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NUMDadmittance");
+  }
+#endif
 
   /* Each time we call this counts as one AC iteration. */
   pDevice->pStats->numIters[STAT_AC] += 1;
@@ -119,6 +136,7 @@ NUMDadmittance(ONEdevice *pDevice, double omega, SPcomplex *yd)
 
 #ifdef KLU
     if (pDevice->matrix->CKTkluMODE) {
+      // Not implemented
       pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUMatrixComplex ;
     } else {
 #endif
@@ -196,6 +214,11 @@ NBJTadmittance(ONEdevice *pDevice, double omega, SPcomplex *yIeVce,
   SPcomplex cOmega, pIeVce, pIcVce, pIeVbe, pIcVbe;
   double startTime;
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NBJTadmittance");
+  }
+#endif
   /* Each time we call this counts as one AC iteration. */
   pDevice->pStats->numIters[STAT_AC] += 1;
 
@@ -315,6 +338,7 @@ NBJTadmittance(ONEdevice *pDevice, double omega, SPcomplex *yIeVce,
 
 #ifdef KLU
     if (pDevice->matrix->CKTkluMODE) {
+      // Not implemented
       pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUMatrixComplex ;
     } else {
 #endif
@@ -552,6 +576,11 @@ NUMDys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yd)
   SPcomplex *y;
 
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NUMDys");
+  }
+#endif
   /*
    * change context names of solution vectors for ac analysis dcDeltaSolution
    * stores the real part and copiedSolution stores the imaginary part of the
@@ -580,6 +609,7 @@ NUMDys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yd)
 
 #ifdef KLU
   if (pDevice->matrix->CKTkluMODE) {
+    // Not implemented
     pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUMatrixComplex ;
   } else {
 #endif
@@ -640,6 +670,11 @@ NBJTys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yIeVce,
   SPcomplex temp, cOmega;
   SPcomplex pIeVce, pIcVce, pIeVbe, pIcVbe;
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NBJTys");
+  }
+#endif
   /*
    * change context names of solution vectors for ac analysis dcDeltaSolution
    * stores the real part and copiedSolution stores the imaginary part of the
@@ -667,6 +702,7 @@ NBJTys(ONEdevice *pDevice, SPcomplex *s, SPcomplex *yIeVce,
 
 #ifdef KLU
   if (pDevice->matrix->CKTkluMODE) {
+    // Not implemented
     pDevice->matrix->SMPkluMatrix->KLUmatrixIsComplex = KLUMatrixComplex ;
   } else {
 #endif
diff --git a/src/ciderlib/twod/twoadmit.c b/src/ciderlib/twod/twoadmit.c
index 8f3a4e1b5..787c10805 100644
--- a/src/ciderlib/twod/twoadmit.c
+++ b/src/ciderlib/twod/twoadmit.c
@@ -16,6 +16,9 @@ Author: 1987 Kartikeya Mayaram, U. C. Berkeley CAD Group
 #include "ngspice/bool.h"
 #include "ngspice/macros.h"
 #include "ngspice/ifsim.h"
+#include "ngspice/cktdefs.h"
+#include "ngspice/ftedefs.h"
+
 #include "twoddefs.h"  
 #include "twodext.h"
 #include "ngspice/cidersupt.h"
@@ -28,6 +31,18 @@ extern IFfrontEnd *SPfrontEnd;
  */
  SPcomplex yTotal;
 
+#ifdef KLU
+static void small_signal_check(char *where)
+{
+  if (ft_curckt->ci_ckt->CKTmode == MODEAC
+   || ft_curckt->ci_ckt->CKTmode == MODEACNOISE) {
+    fprintf(stderr, "Error: CIDER %s small signal simulation is not (yet) supported with 'option klu'.\n", where);
+    fprintf(stderr, "    Use 'option sparse' instead.\n");
+    controlled_exit(1);
+  }
+}
+#endif
+
 int
 NUMD2admittance(TWOdevice *pDevice, double omega, SPcomplex *yd)
 {
@@ -42,6 +57,11 @@ NUMD2admittance(TWOdevice *pDevice, double omega, SPcomplex *yd)
   bool SORFailed;
   double startTime;
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NUMD2admittance");
+  }
+#endif
   /* Each time we call this counts as one AC iteration. */
   pDevice->pStats->numIters[STAT_AC] += 1;
 
@@ -192,6 +212,11 @@ NBJT2admittance(TWOdevice *pDevice, double omega, SPcomplex *yIeVce,
   SPcomplex cOmega;
   double startTime;
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NBJT2admittance");
+  }
+#endif
   /* Each time we call this counts as one AC iteration. */
   pDevice->pStats->numIters[STAT_AC] += 1;
 
@@ -416,6 +441,11 @@ NUMOSadmittance(TWOdevice *pDevice, double omega, struct mosAdmittances *yAc)
   SPcomplex *y, cOmega;
   double startTime;
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NUMOSadmittance");
+  }
+#endif
   /* Each time we call this counts as one AC iteration. */
   pDevice->pStats->numIters[STAT_AC] += 1;
 
@@ -1151,6 +1181,11 @@ NUMD2ys(TWOdevice *pDevice, SPcomplex *s, SPcomplex *yIn)
   bool deltaVContact = FALSE;
   SPcomplex temp, cOmega;
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NUMD2ys");
+  }
+#endif
   /*
    * change context names of solution vectors for ac analysis dcDeltaSolution
    * stores the real part and copiedSolution stores the imaginary part of the
@@ -1252,6 +1287,11 @@ NBJT2ys(TWOdevice *pDevice, SPcomplex *s, SPcomplex *yIeVce, SPcomplex *yIcVce,
   SPcomplex pIeVce, pIcVce, pIeVbe, pIcVbe;
   SPcomplex temp, cOmega;
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NBJT2ys");
+  }
+#endif
   pDevice->solverType = SLV_SMSIG;
   rhsReal = pDevice->rhs;
   rhsImag = pDevice->rhsImag;
@@ -1375,6 +1415,11 @@ NUMOSys(TWOdevice *pDevice, SPcomplex *s, struct mosAdmittances *yAc)
   SPcomplex *y;
   SPcomplex temp, cOmega;
 
+#ifdef KLU
+  if (pDevice->matrix->CKTkluMODE) {
+    small_signal_check("NUMOSys");
+  }
+#endif
   pDevice->solverType = SLV_SMSIG;
   rhsReal = pDevice->rhs;
   rhsImag = pDevice->rhsImag;
diff --git a/src/ciderlib/twod/twocont.c b/src/ciderlib/twod/twocont.c
index 4e79d8415..1fce53439 100644
--- a/src/ciderlib/twod/twocont.c
+++ b/src/ciderlib/twod/twocont.c
@@ -1149,7 +1149,8 @@ TWObindCSC (TWOdevice *pDevice)
 
 	    }
 	  } else { /* Horizontal Slice */
-	    if ( nIndex <= 1 ) { /* Top Side */
+	    //if ( nIndex <= 1 ) { /* Top Side */
+ 	    if ( nIndex == 0 || nIndex == 3 ) { /* Left (Top?) Side : bug 483 */
 
               CREATE_KLU_BINDING_TABLE_CIDER(fNPsiIn, fNPsiInBinding, nEqn, psiEqnInM) ;
               CREATE_KLU_BINDING_TABLE_CIDER(fNPsiInP1, fNPsiInP1Binding, nEqn, psiEqnInP) ;
diff --git a/src/ciderlib/twod/twoncont.c b/src/ciderlib/twod/twoncont.c
index 194339dfe..58c984d00 100644
--- a/src/ciderlib/twod/twoncont.c
+++ b/src/ciderlib/twod/twoncont.c
@@ -786,7 +786,8 @@ TWONbindCSC (TWOdevice *pDevice)
 
 	    }
 	  } else { /* Horizontal Slice */
-	    if ( nIndex <= 1 ) { /* Top Side */
+	    //if ( nIndex <= 1 ) { /* Top Side */
+	    if ( nIndex == 0 || nIndex == 3 ) { /* Left (Top?) Side : bug 483 */
 
               CREATE_KLU_BINDING_TABLE_CIDER(fNPsiIn, fNPsiInBinding, nEqn, psiEqnInM) ;
               CREATE_KLU_BINDING_TABLE_CIDER(fNPsiInP1, fNPsiInP1Binding, nEqn, psiEqnInP) ;
diff --git a/src/ciderlib/twod/twopcont.c b/src/ciderlib/twod/twopcont.c
index a3b1e6dd5..41e8811f8 100644
--- a/src/ciderlib/twod/twopcont.c
+++ b/src/ciderlib/twod/twopcont.c
@@ -786,7 +786,8 @@ TWOPbindCSC (TWOdevice *pDevice)
 
 	    }
 	  } else { /* Horizontal Slice */
-	    if ( nIndex <= 1 ) { /* Top Side */
+	    //if ( nIndex <= 1 ) { /* Top Side */
+ 	    if ( nIndex == 0 || nIndex == 3 ) { /* Left (Top?) Side : bug 483 */
 
               CREATE_KLU_BINDING_TABLE_CIDER(fPPsiIn, fPPsiInBinding, pEqn, psiEqnInM) ;
               CREATE_KLU_BINDING_TABLE_CIDER(fPPsiInP1, fPPsiInP1Binding, pEqn, psiEqnInP) ;
diff --git a/src/spicelib/devices/nbjt/nbjtset.c b/src/spicelib/devices/nbjt/nbjtset.c
index 34f62293d..9656ef443 100644
--- a/src/spicelib/devices/nbjt/nbjtset.c
+++ b/src/spicelib/devices/nbjt/nbjtset.c
@@ -44,7 +44,7 @@ NBJTsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *states)
   ONEmaterial *pM, *pMaterial = NULL, *materialList = NULL;
   double startTime;
 
-#ifdef KLU
+#if defined(KLU) && defined(NOT_WITH_CIDER)
     if (ckt->CKTkluMODE) {
         fprintf(stderr, "Error: CIDER simulation is not (yet) supported with 'option klu'.\n");
         fprintf(stderr, "    Use 'option sparse' instead.\n");
diff --git a/src/spicelib/devices/nbjt2/nbt2set.c b/src/spicelib/devices/nbjt2/nbt2set.c
index 3164d329f..d700fa2b0 100644
--- a/src/spicelib/devices/nbjt2/nbt2set.c
+++ b/src/spicelib/devices/nbjt2/nbt2set.c
@@ -47,7 +47,7 @@ NBJT2setup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *states)
   DOPtable *dopTableList = NULL;
   double startTime;
 
-#ifdef KLU
+#if defined(KLU) && defined(NOT_WITH_CIDER)
     if (ckt->CKTkluMODE) {
         fprintf(stderr, "Error: CIDER simulation is not (yet) supported with 'option klu'.\n");
         fprintf(stderr, "    Use 'option sparse' instead.\n");
diff --git a/src/spicelib/devices/numd/numdset.c b/src/spicelib/devices/numd/numdset.c
index 365d290d4..426629062 100644
--- a/src/spicelib/devices/numd/numdset.c
+++ b/src/spicelib/devices/numd/numdset.c
@@ -46,7 +46,7 @@ NUMDsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *states)
   double startTime;
 
 
-#ifdef KLU
+#if defined(KLU) && defined(NOT_WITH_CIDER)
     if (ckt->CKTkluMODE) {
         fprintf(stderr, "Error: CIDER simulation is not (yet) supported with 'option klu'.\n");
         fprintf(stderr, "    Use 'option sparse' instead.\n");
diff --git a/src/spicelib/devices/numd2/nud2set.c b/src/spicelib/devices/numd2/nud2set.c
index a35a3132f..79f0e3ddd 100644
--- a/src/spicelib/devices/numd2/nud2set.c
+++ b/src/spicelib/devices/numd2/nud2set.c
@@ -47,7 +47,7 @@ NUMD2setup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *states)
   DOPtable *dopTableList = NULL;
   double startTime;
 
-#ifdef KLU
+#if defined(KLU) && defined(NOT_WITH_CIDER)
     if (ckt->CKTkluMODE) {
         fprintf(stderr, "Error: CIDER simulation is not (yet) supported with 'option klu'.\n");
         fprintf(stderr, "    Use 'option sparse' instead.\n");
diff --git a/src/spicelib/devices/numos/nummset.c b/src/spicelib/devices/numos/nummset.c
index 7e9838056..bba47089b 100644
--- a/src/spicelib/devices/numos/nummset.c
+++ b/src/spicelib/devices/numos/nummset.c
@@ -47,7 +47,7 @@ NUMOSsetup(SMPmatrix *matrix, GENmodel *inModel, CKTcircuit *ckt, int *states)
   DOPtable *dopTableList = NULL;
   double startTime;
 
-#ifdef KLU
+#if defined(KLU) && defined(NOT_WITH_CIDER)
     if (ckt->CKTkluMODE) {
         fprintf(stderr, "Error: CIDER simulation is not (yet) supported with 'option klu'.\n");
         fprintf(stderr, "    Use 'option sparse' instead.\n");