Enable q_sfc to be calculated by each surface model

[LenkaNovak]

Necessary to be able to use Bucket's q_sfc, which takes into account its internal values of snow and bare soil are fractions. Here enable this, and apply the same procedure for the other surfaces.

Major changes

• coupler can calculate rho_sfc
• all surface models will be able to hold their own rho_sfc and calculate their own q_sfc (including the Stubs)
• correct specification of Ice() when calculating q_vap_saturation_generic (before only Liquid() was enabled)
• initial conditions are re-specified to be more consistent with atmos, enabling more stable runs

To do

• enable coupler to calculate rho_sfc based on surface and atmos states
• add calculation of q_sfc to the remaining component models
• set the new q_sfc in the compute_atmos_turbulent_fluxes! function
• resolve grid imprinting problem when calling q_sfc from land - due to atmos init conditions being too abrupt
  • modified surface initial conditions to be more consistent with the new atmos conditions
  • ensure and document the correct order of field exchange initialization
• resolve why nondeterministic ci runs - atmos now has a random initial perturbation
• add tests for new src functions

QA

• ensure tests are passing (physical, conservation, performance)