icepack_atmo: use Jordan stability function for [TQ]ref

In the iteration in icepack_atmo::atmo_boundary_layer, we compute the
stability function in the stable case following Jordan et al 1999 (see
compute_stability_function and psi_stable_jordan). However, when
computing the diagnostic temperature and humidity at the reference level
at the end of the subroutine, we do not use the Jordan formulation and
instead use

     \psi_m = -5 \zeta

See for example p.40 of Kauffman & Large [1].

It is unclear if this is an implementation error [2], or if this is
deliberately done to be consistent with the atmospheric model used in
CESM and E3SM [3]. In our in-house version of CICE4, we are using the
Jordan formulation also for the diagnostic variables [4].

To bring CICE6 closer to our in-house CICE4, adjust the computation of
the stability function at the diagnostic level (psix2) to also use the
Jordan stability formulation for the stable case. This results in
relative differences of the order of 0.3 % for Tref and 6.0 % for Qref.

[1] https://github.com/CICE-Consortium/CICE/blob/master/doc/PDF/KL_NCAR2002.pdf
[2] CICE-Consortium#361 (comment)
[3] CICE-Consortium#361 (comment)
[4] https://gitlab.science.gc.ca/concepts/concepts/blob/master/CICECMC/cice4.0_cmc/source/ice_atmo.F90#L449-469

columnphysics/icepack_atmo.F90

@@ -368,7 +368,7 @@ subroutine atmo_boundary_layer (sfctype,            &
       else
          hols  = hols*zTrf/zlvl
       endif
-      psix2 = -c5*hols*stables + (c1-stables)*psi_scalar_unstable(hols)
+      psix2 = psi_stable_jordan(hols)*stables + (c1-stables)*psi_scalar_unstable(hols)
       fac   = (rh/vonkar) &
             * (alzs + al2 - psixh + psix2)
       Tref  = potT - delt*fac