I'm curious why I don't see a gate current contribution, at least for charge. I see these contributions:

    I(D)<+ TYPE * (Ids + ddt(2.0*Qc_on+Qext_d));
    I(S)<+ TYPE * (-Ids + ddt(2.0*Qc_ambi+Qext_s));
 

but it seems this will lead to KCL violations: if I change the drain voltage, Qext_d will change per:

    Qext_d = -(Vg-Vd)*Cext*wch;   // External charges in channel (drain)

but there will not be a corresponding change to the source charge, meaning it will look like the device is sourcing or sinking current to some global reference node (ground).