function [Dne,Vpe,Vpep,Sne,...
    dDne_dx,dDne_dxdot,dDne_du,...
    dVpe_dx,dVpe_dxdot,dVpe_du,...
    dVpep_dx,dVpep_dxdot,dVpep_du,...
    dSne_dx,dSne_dxdot,dSne_du] = density_transport(model,geop,mode)

zz = zeros(model.rgrid.nrhogauss,1);
dzz_dx = zeros(model.rgrid.nrhogauss,model.dims.nx);
dzz_du = zeros(model.rgrid.nrhogauss,model.dims.nu);

Dne = 1*ones(size(zz));
Vpe = -5*ones(size(zz));
Vpep = zz;
dDne_dx = dzz_dx;
dDne_dxdot = dzz_dx;
dDne_du = dzz_du;
dVpe_dx = dzz_dx;
dVpe_dxdot = dzz_dx;
dVpe_du = dzz_du;
dVpep_dx = dzz_dx;
dVpep_dxdot = dzz_dx;
dVpep_du = dzz_du;

Sne = zz;
dSne_dx = dzz_dx;
dSne_dxdot = dzz_dx;
dSne_du = dzz_du;

switch mode
    case 'dummy'
        % return defaults
       return 
    case 'benchmark'
        % From Y.S. Na ITPA-IOS benchmark effort
        rhogauss = model.rgrid.rhogauss;
        R0 = 6.20;
        PhiN = rhogauss.^2;
        PhiNped = 0.88;
        icore = (PhiN<PhiNped);
        
        D0 = 0.5;
        D1 = 1.0;
        D2 = 0.11;
        V0 = 1.385;

        Dne(icore) = D0 + D1*(PhiN(icore)).^2;
        Dne(~icore) = D2;
         
        v0 = V0*rhogauss.*Dne/R0;

        Vpe = -v0.*(geop.g1./(geop.g0.*geop.Vp)); % = v0*<gradrho^2>/<gradrho> = <gradV^2>/<gradV>/dVdrho
        Vpep = [0;diff(Vpe)./diff(rhogauss)];
        
        S0 = 7.5e20;
        C = 0.25e24;
        d = 0.225;
        Sne_edge = 0.5*S0*exp(15*(PhiN - 1));
        Sne_core = 0.5*C*d^2* PhiN.^(6.5).*(1-PhiN).^(8.5)./(d^2+(PhiN-0.5).^2) ;
      
        Sne = (Sne_core+Sne_edge);
    otherwise
        error('not supported')
end
return
