icepack_atmo: compute Jordan stability function in psi_stable_jordan

The computation of the Jordan stability function, used for the stable
case, is currently done in subroutine compute_stability_function,
whereas the computations for the unstable case are done in their own
functions, psi_scalar_unstable and psi_momentum_unstable.

Move the computation of the (stable) Jordan stability function to its
own function, psi_stable_jordan, and call it in
compute_stability_function. This will allow us to call psi_stable_jordan
directly in a following commit.

columnphysics/icepack_atmo.F90

@@ -1054,9 +1054,7 @@ subroutine compute_stability_function(qty, hol, stable, psi)
 
       stable = p5 + sign(p5 , hol)
 
-      psi_stable = -(0.7_dbl_kind*hol &
-                 + 0.75_dbl_kind*(hol-14.3_dbl_kind) &
-                 * exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind)
+      psi_stable = psi_stable_jordan(hol)
 
       if(trim(qty) == 'momentum') then
          psi_unstable = psi_momentum_unstable(hol)
@@ -1104,6 +1102,18 @@ real(kind=dbl_kind) function psi_scalar_unstable(hol)
 
       end function psi_scalar_unstable
 
+!=======================================================================
+
+      real(kind=dbl_kind) function psi_stable_jordan(hol) result(psi_stable)
+
+      real(kind=dbl_kind), intent(in) :: hol
+
+      psi_stable = -(0.7_dbl_kind*hol &
+                 + 0.75_dbl_kind*(hol-14.3_dbl_kind) &
+                 * exp(-0.35_dbl_kind*hol) + 10.7_dbl_kind)
+
+      end function psi_stable_jordan
+
 !=======================================================================
 
       real(kind=dbl_kind) function capital_X(hol)