Current and future weak lensing surveys contain significant information about our universe. However, their optimal cosmological analysis is challenging, with traditional analyses often resulting in information loss due to reliance on summary statistics like two-point functions. While deep learning methods offer promise in capturing the complex non-linear features of these cosmological fields, they often suffer from issues such as inadequate uncertainty quantification, susceptibility to distribution shifts, and interpretability limitations, which hinder their scientific applicability. In this talk, I propose a novel approach leveraging generative probabilistic modeling with Normalizing Flows to learn the data likelihood function at the field level, facilitating more effective cosmological information extraction. This framework not only enables anomaly detection of distribution shifts (e.g., noise miscalibration and baryonic effect) to improve the robustness of the analysis, but also fostering interpretability via generated samples. I will also discuss incorporating physical prior knowledge, such as symmetries and multiscale structure, into the model architectures to improve their generalization capabilities. Finally, I will talk about our ongoing work on applying this model to the field-level cosmic shear analysis for Hyper Suprime-Cam (HSC).