Non-smoothness arises at cloud edge because, in moist thermodynamics, the thermodynamic properties of the atmosphere are different inside a cloud versus in clear air. In particular, inside a cloud, the vapor pressure of water is constrained by the saturation vapor pressure, which acts as a threshold. Due to this threshold, while the water vapor mixing ratio may vary continuously across cloud edge, its derivatives are not necessarily continuous at cloud edge. Similarly, non-smoothness also arises for buoyancy and other variables. Consequently, this non-smoothness in buoyancy and other variables can cause a degraded accuracy in computational simulations. Here we consider special treatment of numerical methods for the interface that arises from phase changes and cloud edges, in order to enhance the accuracy and potentially achieve second-order accuracy. Numerical solutions are computed for the moist non-precipitating Boussinesq equations as an idealized cloud-resolving model (CRM) with phase changes of water, i.e., with cloud formation. Convergence tests, both spatial and temporal, are conducted to measure the numerical error as the grid spacing and time step are refined. While approximately second-order accuracy is seen in root-mean-square error, the accuracy is degraded in the maximum error.