function [neGR] = eval_neGR(x,g,v,model)
% function [Pthr] = eval_neGR(x,g,v,mode
% Calculation of the Greenwald density limit

% Plasma current
Ip = eval_Ipl(x,g,v,model,false);
% Inverse aspect ratio
epsilon = g(model.geom.ind_epsilon);
% Minor radius
a = epsilon(end)*model.equi.R0;
% Greenwald density limit
neGR = 1e+20*Ip(end,:)/(1e+6*pi*a^2);

return


