diff --git a/columnphysics/icepack_atmo.F90 b/columnphysics/icepack_atmo.F90
index c5d7712c..f1a661ac 100644
--- a/columnphysics/icepack_atmo.F90
+++ b/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)