A model of charged hole-pair bosons with long-range Coulomb interactions and very weak interlayer coupling is used to calculate the order parameter ⌽ of underdoped cuprates. Model parameters are extracted from experimental superfluid densities and plasma frequencies. The temperature dependence ⌽͑T͒ is characterized by a "trapezoidal" shape. At low temperatures, it declines slowly due to harmonic phase fluctuations which are suppressed by anisotropic plasma gaps. Above the single layer Berezinski-Kosterlitz-Thouless temperature, ⌽͑T͒ falls rapidly toward the three-dimensional transition temperature. The theoretical curves are compared to c-axis superfluid density data by Kitano et al. ͓J. Low Temp. Phys. 117, 1241 ͑1999͔͒ and to the transverse nodal velocity measured by angular resolved photoemission spectra on BSCCO samples by Lee et al. ͓Nature ͑London͒ 450, 81 ͑2007͔͒ and by Kanigel et al. ͓Phys. Rev. Lett. 99, 157001 ͑2007͔͒.