Added the option of including the atmospheric or ice pressure and sea surface
height displacements from the global reference height in the pressures used in
the density calculations for Boussinesq pressure gradient calculations.  Note
that the full pressures were already being used everywhere apart from the
calculation of the equation of state.  This capability is controlled by the new
runtime parameter SSH_IN_EOS_PRESSURE_FOR_PGF.  This commit changes the Z_0p
argument to int_density_dz and 8 other routines (int_density_dz_generic_pcm,
int_density_dz_generic_plm, int_density_dz_generic_ppm, analytic_int_density_dz,
int_density_dz_wright, int_density_dz_wright_full, int_density_dz_wright_red
and int_density_dz_linear) from scalars into 2-d arrays, as were the internal
z0pres arrays in most of these routines.  By default, all answers are bitwise
identical, but there is a new runtime parameter in the MOM_parameter_doc files
for Boussinesq cases.

  Rearranged the calculation of the layer average pressure used in the density
integrals to have one less multiply and be clearer to read.  This could change
answers at roundoff if REFERENCE_HEIGHT is nonzero (it is 0 by default and in
all known cases) or if SSH_IN_EOS_PRESSURE_FOR_PGF is set to true.  Because the
former is mostly used for approximate self-consistency testing and the latter
option has just been added, it is unlikely that answers will change for any
production runs.