Mimicking the configuration for those machines in CICE, port Icepack to
'daley' and 'banting'.

* icepack_atmo: add stability functions with more explicit names

The two functions 'psimhu' and 'psixhu' in icepack_atmo compute the
stability functions in the unstable case for momentum and scalars (heat
and water), respectively. They are both function of X := (1 - 16
\zeta)^(1/4) (see [1], section 8.1.1). In the code X is computed first
('xqq') and passed to 'psimhu' and 'psixhu'.

Add two new functions with more explicit names, 'psi_momentum_unstable'
and 'psi_scalar_unstable', that take \zeta ('hol') as argument,
compute X internally and return the evaluated stability function.

A subsequent commit will refactor the code to use them.

[1] B.G. Kauffman and W.G. Large. The CCSM coupler, version 5.0.1. 2002.
URL: https://github.com/CICE-Consortium/CICE/blob/master/doc/PDF/KL_NCAR2002.pdf.

* icepack_atmo: use 'psi_{momentum,scalar}_unstable'

Use the new functions 'psi_momentum_unstable' and 'psi_scalar_unstable'
introduced in the preceding commit to compute the stability function
evaluated at 'hol'. Remove duplicated computation  of 'xqq'.

Also remove the now unused functions 'psimhu' and 'psixhu'.

* icepack_atmo: refactor 'atmo_boundary_layer'

In a subsequent commit, we will support receiving the momentum and
scalar atmospheric variables at different levels. In order to minimize
code duplication, factor out the computation of the stability parameter
'hol' and the stability function 'psi[mx]h' in function
'compute_stability_parameter' and subroutine
'compute_stability_function'.

'compute_stability_function' also makes available the unit step function
evaluated at 'hol', 'stable', since it is used at the end of
'atmo_boundary_layer' to compute the diagnostic temperature and
humidity.

* icepack_atmo: add support for staggered atmospheric levels

In some applications, we might want the atmospheric input variables
(wind, temperature and humidity) to be given at different vertical
levels. For example, we might want to receive the winds at some level
and the scalar variables (temperature and humidity) at another.

To support this use case, add an optional argument 'zlvs' to subroutine
'icepack_atm_boundary', and pass it down to 'atmo_boundary_layer' also
as an optional argument, so that existing applications do not have to be
changed if they do not want to use this feature.

The name 'zlvs' is chosen to mirror the existing variable 'zlvl', but
with the suffix 's' for 'scalars', mimicking the existing convention for
several variables in 'atmo_boundary_layer'.

In 'atmo_boundary_layer', add new local variables 'alzs' and 'hols' to
hold quantities computed at zlvs.  Also, for clarity, rename 'alz' to
'alzm' and 'hol' to 'holm', where the 'm' suffix stands for 'momentum',
again mimicking existing convention in the subroutine.

If 'zlvs' is not passed to 'atmo_boundary_layer', it is assumed that all
variables are given at the same level (zlvl), so in that case there is
no change in the results.

* icepack_therm_vertical: add support for staggered atmospheric level

The previous commit added an optional 'zlvs' argument to
'icepack_atm_boundary' to allow atmospheric input variables to be given
at different vertical levels.

Also add 'zlvs' as an optional argument to 'icepack_step_therm1', to
make the functionality available via the Icepack interfaces.

* doc: document staggered atmospheric levels

Add an index entry for 'zlvs' and mention it in the 'Atmopshere' section
of the science guide.

* code will write warning message but not abort if
hicen_init(n+1) <= hicen_init(n)