function [g0] = eval_g0p(~,g,~,model,gauss)
% function [g0] = eval_g0p(x,g,v,model,whichgrid)
% dg0/drho = d/drho <|nabla V|>
% #codegen

% g0p = gg0p / H - gg0/H^2* Hp

H = eval_H([],g,[],model,gauss); % (dpsi_dV)^2

if ~gauss
    gg0 = model.geom.LamG*g(model.geom.ind_gg0,:);
    gg0p = model.geom.LamGp*g(model.geom.ind_gg0,:);
    Hp   = model.geom.LamGp*g(model.geom.ind_H,:);
else
    gg0 = (model.geom.LamGgauss*g(model.geom.ind_gg0,:));
    gg0p = (model.geom.LamGpgauss*g(model.geom.ind_gg0,:));
    Hp   = model.geom.LamGpgauss*g(model.geom.ind_H,:);
end

diverted = (H(end,:)==0);
g0 = gg0p./H - gg0./(H.*H).*Hp;
g0(end,diverted) = inf;

return