We investigate the numerical approximation of solutions to some variational inequalities modeling the humid atmosphere when the saturation of water vapor in the air is taken into account. In order to overcome the difficulties caused by the constraints on the humidity $q$ ($0 \leq q \leq 1$) and the discontinuity in the variational inequalities, we construct a penalized and regularized implicit Euler method. We manage to show that the approximation functions associated with the numerical scheme converge to the solutions of the variational inequalities through deriving various delicate a priori estimates and by using compactness arguments.