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 into account the magnitude distribution of both cluster and field galaxies as well as the radial distribution of galaxies in clusters 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 mainly developed for the forthcoming Euclid mission . Median purity and completeness values of { ( 55 ^ { +17 } _ { -15 } ) \% } and { ( 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 \overline { \mathsf { p } } = { 56 } \% and \overline { \mathsf { c } } = { 93 } \% are consistent with the median and have negligible sub-percent uncertainties . Accurate photometric redshifts ( \sigma ( ( z _ { s } - z _ { p } ) / ( 1 + z _ { s } ) ) \lesssim 0.05 ) and robust estimates for the cluster redshift and cluster center coordinates are required . The dependence of the assignments on photometric redshift accuracy , galaxy magnitude and distance from the halo center , and halo properties such as mass , richness , and redshift are investigated . Variations in the mean values of both purity and completeness are globally limited to a few percent . The largest departures from the mean values are found for galaxies associated with distant z \gtrsim 1.5 halos , faint ( \sim 0.25 L _ { \ast } ) galaxies , and those at the outskirts of the halo ( at cluster-centric projected distances \sim r _ { 200 } ) for which the purity is decreased , \Delta \mathsf { p } \simeq 20 \% at most , with respect to the mean value . The proposed method is applied to derive accurate richness estimates . A statistical comparison between the true ( N _ { true } ) vs. estimated richness ( \lambda = \sum P _ { mem } ) yields on average to unbiased results , Log ( \lambda / N _ { true } ) = -0.0051 \pm 0.15 . The scatter around the mean of the logarithmic difference between \lambda and the halo mass is 0.10 dex for massive halos \gtrsim 10 ^ { 14.5 } ~ { } M _ { \odot } . Our estimates could therefore be useful to constrain the cluster mass function and to calibrate independent cluster mass estimates such as those obtained from weak lensing , Sunyaev-Zel’dovich , 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 , Pan-STARRS , DES , LSST , and Euclid .