A computational algorithm for tracking liquid- vapor interfaces moving with curvature-dependent velocity, based on, the Gibbs energy formulation. and the level-set approach is described. The algorithm is used for finding equilibrium distributions of condensate in several samples of 3-D reconstructed porous media under varying relative saturation of the surrounding atmosphere. Primary and higher-order hysteresis cycles are calculated and their dependence on various physico-chemical parameters (temperature, wetting angle) as well as on geometrical properties of the porous medium is investigated. An explanation of the influence of the morphology of a porous medium on the character of corresponding hysteresis loops is provided.