  We introduce a new effective strategy to assign group and cluster membership
probabilities $P_{mem}$ to galaxies using photometric redshift information.
Large dynamical ranges both in halo mass and cosmic time are considered. The
method takes the magnitude distribution of both cluster and field galaxies as
well as the radial distribution of galaxies in clusters into account using a
non-parametric formalism and relies on Bayesian inference to take photometric
redshift uncertainties into account. We successfully test the method against
1,208 galaxy clusters within redshifts $z=0.05-2.58$ and masses
$10^{13.29-14.80}~M_\odot$ drawn from wide field simulated galaxy mock catalogs
developed for the Euclid mission. Median purity $(55^{+17}_{-15})\%$ and
completeness $(95^{+5}_{-10})\%$ are reached for galaxies brighter than
0.25$L_\ast$ within $r_{200}$ of each simulated halo and for a statistical
photometric redshift accuracy $\sigma((z_s-z_p)/(1+z_s))=0.03$. The mean values
$p=56\%$ and $c=93\%$ have sub-percent uncertainties. Accurate photometric
redshifts ($\sigma((z_s-z_p)/(1+z_s))\lesssim0.05$) and robust estimates for
the cluster redshift and the center coordinates are required. The method is
applied to derive accurate richness estimates. A statistical comparison between
the true ($N_{\rm true}$) vs estimated richness ($\lambda=\sum P_{mem}$) yields
on average to unbiased results, $Log(\lambda/N_{\rm true})=-0.0051\pm0.15$. The
scatter around the mean of the logarithmic difference between $\lambda$ and the
halo mass is 0.10~dex, for massive halos $\gtrsim10^{14.5}~M_\odot$. Our
estimates could be useful to calibrate independent cluster mass estimates from
weak lensing, SZ, and X-ray studies. Our method can be applied to any list of
galaxy clusters or groups in both present and forthcoming surveys such as SDSS,
CFHTLS, DES, LSST, and Euclid.
