We introduce a new effective strategy to assign group and cluster membership probabilities Pmem 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 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.25L * within r200 of each simulated halo and for a statistical photometric redshift accuracy σ((zs −zp)/(1+zs)) = 0.03.The mean values p = 56% and c = 93% are consistent with the median and have negligible sub-percent uncertainties. Accurate photometric redshifts (σ((zs − zp)/(1 + zs)) 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 1.5 halos, faint (∼ 0.25 L * ) galaxies, and those at the outskirts of the halo (at cluster-centric projected distances ∼ r200) for which the purity is decreased, ∆p 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 (Ntrue) vs. estimated richness (λ = Pmem) yields on average to unbiased results, Log(λ/Ntrue) = −0.0051 ± 0.15. The scatter around the mean of the logarithmic difference between λ and the halo mass is 0.10 dex for massive halos 10 14.5 M . 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.