This paper is aimed at a numerical simulation of the main phase of two-ribbon flares by solving viscous MHDequations in terms of a multistep implicit scheme. The numerical results clearly show the transition from linear reconnection of tearing modes to quasisteady re-connection in the neutral sheet, the formation of post-flare loops, and the ejection of plas-moids.Two cases are considered respectively: one includes and the other does not include Joule heating in the energy equation. It is shown that including Joule heating leads to a substantial increase of temperature (with a factor of 2-3), but a negligible variation of velocity for the plasma within the neutral sheet. The plasma motions are always subsonic for both cases and therefore, no fast shocks are formed.The numerical results of including heating show two new features. One is a reduction and a tendency to saturation of the reconnection rate which is caused by the expansion of the high-density plasma in the neutral sheet, leading to an increase in the effective thickness of the neutral sheet. The other is the formation of ascending and descending plasmoids, and the coalesces of the latter with the post-flare loops.