fix maxiners

srcF/sfx_sheba_noniterative.f90

@@ -230,6 +230,7 @@ contains
         real sigma_w, sigma_r
         real S_mean, S_salt
         real z0_m1
+        integer i
 #ifdef SFX_CHECK_NAN
         real NaN
 #endif
@@ -288,36 +289,33 @@ contains
                 h0_m, h0_t, B)
         
         !  --- define snow consentration
-        call get_sigma(sigma_r, sigma_w, rho_air, rho_s)
-        call get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
-        call get_h_salt(h_salt, u_dyn0)
-        call get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
-        call get_S_mean(S_mean,  S_salt, h_salt, h, w_snow, u_dyn0)
+        !call get_sigma(sigma_r, sigma_w, rho_air, rho_s)
+        !call get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
+        !call get_h_salt(h_salt, u_dyn0)
+        !call get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
+        !call get_S_mean(S_mean,  S_salt, h_salt, h, w_snow, u_dyn0)
                       
-        deltaS=S_salt-S_mean
+        !deltaS=S_salt-S_mean
     
-        Ri_sn = (g * sigma_r * deltaS * h) / U**2
-        
+        !Ri_sn = (g * sigma_r * deltaS * h) / U**2
+                Ri_sn=0.0
         
-
         ! --- get the fluxes
         ! ----------------------------------------------------------------------------
-        !write(*,*) 'surface_type', surface_type
          if (Rib > 0.0) then
             ! --- stable stratification block
 
             !   --- restrict bulk Ri value
             !   *: note that this value is written to output
 !            Rib = min(Rib, Rib_max)
-        do i = 1, maxiters   
+        do i = 1, numerics%maxiters_convection  
         if (surface_type == surface_snow) then    
-            write(*,*) 'sfx_snow1', Ri_sn, Rib, U, deltaS, i
-            !write(*,*) 'sfx_snow2', deltaS, u_dyn0, S_mean
+            !write(*,*) 'RIsnow', Rib, Ri_sn   
             Rib=Rib+Ri_sn 
         endif
 
             call get_zeta(zeta, Rib, h, z0_m, z0_t)
-
+            
             call get_psi_stable(psi_m, psi_h, zeta, zeta)
             call get_psi_stable(psi0_m, psi0_h, zeta * z0_m / h, zeta * z0_t / h)
 
@@ -327,16 +325,20 @@ 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))
             
+            
+            call get_sigma(sigma_r, sigma_w, rho_air, rho_s)
+            call get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
             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
+            !write(*,*) 'sfx_snow1', Ri_sn, Rib, Udyn, S_salt, S_mean
+        else if (Rib < Rib_conv_lim) then
             ! --- strong instability block
 
             call get_psi_convection(psi_m, psi_h, zeta, Rib, &
@@ -384,10 +386,15 @@ contains
         Pr_t_inv = phi_m / phi_h
 
         ! --- setting output
+        !sfx = sfxDataType(zeta = zeta, Rib = Rib, &
+        !    Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
+        !    Rib_conv_lim = Rib_conv_lim, &
+        !    Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
+
         sfx = sfxDataType(zeta = zeta, Rib = Rib, &
-            Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
-            Rib_conv_lim = Rib_conv_lim, &
-            Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)
+            Re = Ri_sn, B = B, z0_m = z0_m, z0_t = z0_t, &
+            Rib_conv_lim = S_mean, &
+            Cm = Cm, Ct = Ct, Km = S_salt, Pr_t_inv = Udyn)
 
     end subroutine get_surface_fluxes
     ! --------------------------------------------------------------------------------
@@ -418,6 +425,7 @@ contains
         real, intent(in) ::  S_salt, h_salt, z, w_snow, u_dyn0         
         ! ----------------------------------------------------------------------------
         S_mean  = S_salt *  (z/h_salt)**(-w_snow/(kappa*u_dyn0))
+        !write(*,*) 'Smean_Udyn', S_mean, u_dyn0
     end subroutine
 
     subroutine get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
@@ -430,6 +438,7 @@ contains
         if (u_dyn0>u_thsnow) then
             qs  = (u_dyn0*u_dyn0-u_thsnow*u_thsnow)/(Csn*u_dyn0*g*h_salt)
             S_salt=qs/(qs+rho_s/rho_air)
+            !write(*,*) 'S_salt', S_salt
         else
             S_salt=0.0
         endif