iterations foe snow Ri

srcF/sfx_sheba_noniterative.f90

@@ -280,7 +280,12 @@ contains
         ! --- define relative height [thermal]
         h0_t = h / z0_t
       
-       
+        ! --- define Ri-bulk
+        Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
+         
+        ! --- define free convection transition zeta = z/L value 
+        call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
+                h0_m, h0_t, B)
         
         !  --- define snow consentration
         call get_sigma(sigma_r, sigma_w, rho_air, rho_s)
@@ -291,14 +296,9 @@ contains
                       
         deltaS=S_salt-S_mean
     
-         
-        ! --- define Ri-bulk
-        Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2
         Ri_sn = (g * sigma_r * deltaS * h) / U**2
-        !write(*,*) 'ri', S_salt, h_salt, h, w_snow, u_dyn0
-        ! --- define free convection transition zeta = z/L value
-        call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
-                h0_m, h0_t, B)
+        
+        
 
         ! --- get the fluxes
         ! ----------------------------------------------------------------------------
@@ -309,9 +309,9 @@ contains
             !   --- restrict bulk Ri value
             !   *: note that this value is written to output
 !            Rib = min(Rib, Rib_max)
-           
+        do i = 1, maxiters   
         if (surface_type == surface_snow) then    
-            write(*,*) 'sfx_snow1', Ri_sn, Rib, U, deltaS
+            write(*,*) 'sfx_snow1', Ri_sn, Rib, U, deltaS, i
             !write(*,*) 'sfx_snow2', deltaS, u_dyn0, S_mean
             Rib=Rib+Ri_sn 
         endif
@@ -326,8 +326,17 @@ contains
 
             Udyn = kappa * U / (log(h / z0_m) - (psi_m - psi0_m))
             Tdyn = kappa * dT * Pr_t_0_inv / (log(h / z0_t) - (psi_h - psi0_h))
-
-        else if (Rib < Rib_conv_lim) then
+            
+            call get_h_salt(h_salt, Udyn)
+            call get_S_salt(S_salt, Udyn, u_thsnow, g, h_salt)
+            call get_S_mean(S_mean,  S_salt, h_salt, h, w_snow, Udyn)
+                      
+            deltaS=S_salt-S_mean
+    
+            Ri_sn = (g * sigma_r * deltaS * h) / U**2
+            
+            end do
+         else if (Rib < Rib_conv_lim) then
             ! --- strong instability block
 
             call get_psi_convection(psi_m, psi_h, zeta, Rib, &