Four flexible water models based an the harmonic and a modified anharmonic form of Toukan and Rahman intramolecular potential models grafted onto 3-site and 4-site intermolecular models were proposed, and the influence of their flexibility and intermolecular potential on the phase behavior were examined. A global optimization of the potential parameters was performed according to the physical properties and IR spectra of water at the ambient state. The critical properties of the four models were estimated by using direct molecular dynamics (MD) simulation of liquid-vapor coexisting phases and fitting the coexisting densities to the law of rectilinear diameter and the scaling law for density with a scaling exponent of beta(c) = 0.325. Based on a similar intermolecular model (the simple point charge model, SPC) we found that the anharmonic model gave a critical point (T-c = 370.3 +/- 3.6 degrees C; rho(c) = 0.32 +/- 0.006 g/cm(3)) which agreed well with the experimental properties (T-c = 374.0 degrees C; rho(c) = 0.322 g/cm(3)), while the harmonic model gave low values (T-c = 337.3 +/- 8.5 degrees C; rho(c) = 0.290 +/- 0.007 g/cm(3)). On the other hand, for anharmonic flexible water models adopting the experimental value, the model grafted onto a 3-site TIPS model (the transferable intermolecular potential functions, TIP3P) did not provide satisfactory reproduction of the critical point (T-c = 320.1 +/- 5.2 degrees C; rho(c) = 0.288 +/- 0.005 g/cm(3)), while the model grafted onto a 4-site TIPS model (modified TIP4P model with M-site on the center of mass, referred to as cm4P) provided good representation of critical properties (T-c = 368.4 +/- 6.2 degrees C; rho(c) = 0.307 +/- 0.005 g/cm(3)). The flexible models with anharmonic form of flexibility grafted the SPC or cm4P model are suitable model for study of water and aqueous solutions in state points ranging from liquid to near- and supercritical region.