icepack_atmo: use separate stability factors for momentum and scalars

In 37f2a17 (Add support for staggered atmospheric levels (CICE-Consortium#364),
2021-05-28), we introduced the subroutine compute_stability_function,
which computes the stability function during the iteration in
'icepack_amto::atmo_boundary_layer'. This subroutine also returns the
unit step function 'stable' at 'hol', since it is used at the end of the
computation to compute the diagnostic temperature and humidity. Note
that we overwrite 'stable' in the second call to
compute_stability_function, but this is OK since we call it first for
momentum and then for scalars, and it is the value for scalars that we
want to use later for diagnostic temperature and humdity.

In a following commit we will want to compute the diagnostic _velocity_
in a different way than what is done now, and we will need the unit step
function computed at the momentum level (which is currently
overwritten). In preparation for that change, add a new variable
'stablem' for the unit step function at the momentum level, and rename
the existing 'stable' variable to 'stables' (for scalars). Adjust the
computation of 'psix2' accordingly.

columnphysics/icepack_atmo.F90

@@ -166,7 +166,8 @@ subroutine atmo_boundary_layer (sfctype,            &
          cp    , & ! specific heat of moist air
          holm  , & ! H (at zlvl  ) over L
          hols  , & ! H (at zlvs  ) over L (if zlvs present)
-         stable, & ! stability factor
+         stablem, & ! stability factor (momentum)
+         stables, & ! stability factor (scalars)
          cpvir , & ! defined as cp_wv/cp_air - 1.
          psixh     ! stability function at zlvl (at zlvs if present) (heat and water)
 
@@ -297,8 +298,8 @@ subroutine atmo_boundary_layer (sfctype,            &
             hols = holm
          endif
 
-         call compute_stability_function('momentum', holm, stable, psimh)
-         call compute_stability_function('scalar'  , hols, stable, psixh)
+         call compute_stability_function('momentum', holm, stablem, psimh)
+         call compute_stability_function('scalar'  , hols, stables, psixh)
 
          ! shift all coeffs to measurement height and stability
          rd = rdn / (c1+rdn/vonkar*(alzm-psimh))
@@ -381,7 +382,7 @@ subroutine atmo_boundary_layer (sfctype,            &
       else
          hols  = hols*zTrf/zlvl
       endif
-      psix2 = -c5*hols*stable + (c1-stable)*psi_scalar_unstable(hols)
+      psix2 = -c5*hols*stables + (c1-stables)*psi_scalar_unstable(hols)
       fac   = (rh/vonkar) &
             * (alzs + al2 - psixh + psix2)
       Tref  = potT - delt*fac