diff --git a/diag/pbldia_new_sfx.f90 b/diag/pbldia_new_sfx.f90
index 3773d0674761774a1e662d0c95ba75b66c0f6ca5..963b30515efe3fb57a72a54dc2e1072043554c7e 100755
--- a/diag/pbldia_new_sfx.f90
+++ b/diag/pbldia_new_sfx.f90
@@ -88,12 +88,15 @@ module pbldia_new_sfx
             TSTAR = AR2(9) * ARDIN(2)
             QSTAR = AR2(9) * ARDIN(3)
       
-           if(zeta.ne.0)then
-      
-      
-           if(zeta.gt.zetalim) L_lim = HIN/zetalim
+           if (zeta.ne.0) then
+               L = HIN/zeta
+               if(zeta.gt.zetalim) then
+                    L_lim = HIN/zetalim
+               else
+                   L_lim = L
+               end if
       
-            L = HIN/zeta
+
       
            call get_psi_sheba(psi_m_hs,psi_h_hs,HIN/L_lim)
            call get_psi_sheba(psi_m,psi_h,HWIND/L_lim)
@@ -144,9 +147,13 @@ module pbldia_new_sfx
             QSTAR = AR2(9) * ARDIN(3)
       
            if(zeta.ne.0)then
-      
-           if(zeta.gt.zetalim) L_lim = HIN/zetalim
-           L = HIN/zeta
+
+               L = HIN/zeta
+               if(zeta.gt.zetalim) then
+                   L_lim = HIN/zetalim
+               else
+                   L_lim = L
+               end if
       
            if(zeta.gt.0)then ! stable stratification
       
@@ -208,9 +215,13 @@ module pbldia_new_sfx
             QSTAR = AR2(9) * ARDIN(3)
       
            if(zeta.ne.0)then
-      
-           if(zeta.gt.zetalim) L_lim=HIN/zetalim
-            L = HIN/zeta
+
+               L = HIN/zeta
+               if(zeta.gt.zetalim) then
+                   L_lim = HIN/zetalim
+               else
+                   L_lim = L
+               end if
       
            if(zeta.gt.0)then ! stable stratification
       
@@ -275,9 +286,13 @@ module pbldia_new_sfx
             QSTAR = AR2(9) * ARDIN(3)
       
            if(zeta.ne.0)then
-      
-           if(zeta.gt.zetalim) L_lim=HIN/zetalim
-            L = HIN/zeta
+
+               L = HIN/zeta
+               if(zeta.gt.zetalim) then
+                   L_lim = HIN/zetalim
+               else
+                   L_lim = L
+               end if
       
            call get_psi_most(psi_m_hs,psi_h_hs,HIN/L_lim)
            call get_psi_most(psi_m,psi_h,HWIND/L_lim)
@@ -323,9 +338,12 @@ module pbldia_new_sfx
             QSTAR = AR2(9) * ARDIN(3)
       
            if(zeta.ne.0)then
-      
-           if(zeta.gt.zetalim) L_lim=HIN/zetalim
-            L = HIN/zeta
+               L = HIN/zeta
+               if(zeta.gt.zetalim) then
+                   L_lim = HIN/zetalim
+               else
+                   L_lim = L
+               end if
       
            call get_psi_esm1(psi_m,psi_h,HIN,HWIND,L_lim)
            UFWIND = psi_m