We consider supersymmetric extensions of the standard model with two pairs of Higgs doublets. We study the possibility of spontaneous $CP$ violation in these scenarios and present a model where the origin of $CP$ violation is soft, with all the complex phases in the Lagrangian derived from complex masses and vacuum expectation values (VEVs) of the Higgs fields. The main ingredient of the model is an approximate global symmetry, which determines the order of magnitude of Yukawa couplings and scalar VEVs. We assume that the terms violating this symmetry are suppressed by powers of the small parameter $��_{PQ}=O(m_b/m_t)$. The tree-level flavor changing interactions are small due to a combination of this global symmetry and a flavor symmetry, but they can be the dominant source of $CP$ violation. All $CP$-violating effects occur at order $��_{PQ}^2$ as the result of exchange of {\it almost}-decoupled extra Higgs bosons and/or through the usual mechanisms with an {\it almost}-real CKM matrix. On dimensional grounds, the model gives $��_K\approx ��_{PQ}^2$ and predicts for the neutron electric dipole moment (and possibly also for $��'_K$) a suppression of order $��_{PQ}^2$ with respect to the values obtained in standard and minimal supersymmetric scenarios. The predicted $CP$ asymmetries in $B$ decays are generically too small to be seen in the near future. The mass of the lightest neutral scalar, the strong $CP$ problem, and possible contributions to the $Z$ decay into $b$ quarks (the $R_b$ puzzle) are also briefly addressed in the framework of this model.

Latex, 31 pages including 6 figures which are a part of the Latex file