fixed zetalim bug in pbldia_new_sfx.f90

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