Linvers for stable case, noiter

srcF/sfx_sheba_noniterative.f90

@@ -228,7 +228,7 @@ contains
 
         real w_snow, deltaS, h_salt
         real sigma_w, sigma_r
-        real S_mean, S_salt
+        real S_mean, S_salt, Linv
         real z0_m1
         integer i
 #ifdef SFX_CHECK_NAN
@@ -315,29 +315,50 @@ contains
         !endif
 
             call get_zeta(zeta, Rib, h, z0_m, z0_t)
+            !write(*,*) 'get_zeta', zeta, h
             
             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)
 
             phi_m = 1.0 + (a_m * zeta * (1.0 + zeta)**(1.0 / 3.0)) / (1.0 + b_m * zeta)
             phi_h = 1.0 + (a_h * zeta + b_h * zeta * zeta) / (1.0 + c_h * zeta + zeta * zeta)
-
+            S_mean=0.0
             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))
+            !write(*,*) 'sfx_before_snow', Udyn, zeta, S_mean
+
+            if (surface_type==3.or.surface_type==6.) then
+                if (Udyn>u_thsnow) then
             
-            
-         !   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)
+            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
-            
+            Linv=Linv*((1-S_mean)/(1+sigma_w*S_mean))+(g*w_snow*sigma_w*S_mean/(Udyn**3.0))/(1+sigma_w*S_mean)
+            zeta = h * Linv 
+            !write(*,*) 'sfx_snow1', Udyn, Linv, zeta, S_mean
+                  
+                 
+             
+            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)
+
+            phi_m = 1.0 + (a_m * zeta * (1.0 + zeta)**(1.0 / 3.0)) / (1.0 + b_m * zeta)
+            phi_h = 1.0 + (a_h * zeta + b_h * zeta * zeta) / (1.0 + c_h * zeta + zeta * zeta)
+
+            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))
+
+            write(*,*) 'sfx_snow2', Udyn, zeta, S_mean, Linv
+            endif 
+            endif
+        
            ! Ri_sn = (g * sigma_r * deltaS * h) / U**2
             
             !end do
-            !write(*,*) 'sfx_snow', Ri_sn, Rib, Udyn, surface_type
+            
         else if (Rib < Rib_conv_lim) then
             ! --- strong instability block
 
@@ -386,15 +407,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 = 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)
+        !    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)
+        !write(*,*) 'Smean_0ut', S_mean
+        sfx = sfxDataType(zeta = zeta, Rib = Rib, &
+            Re = Linv, 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
     ! --------------------------------------------------------------------------------
@@ -425,7 +446,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
+        !write(*,*) 'Smean_func', S_mean, u_dyn0
     end subroutine
 
     subroutine get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)

srcF/sfx_z0m_all_surface.f90

@@ -84,7 +84,7 @@ module sfx_z0m_all_surface
 
         ! --- define dynamic velocity in neutral conditions
         u_dyn0 = Uc / c
-        write(*,*) 'ch', z0_m, u_dyn0
+        !write(*,*) 'ch', z0_m, u_dyn0
     end subroutine
     ! --------------------------------------------------------------------------------
     
@@ -204,7 +204,7 @@ subroutine get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map)
         z0_m=z0m_map  
         h0_m = h / z0_m
         u_dyn0 = U * kappa / log(h0_m)
-        write(*,*) 'map', z0_m, u_dyn0
+        !write(*,*) 'map', z0_m, u_dyn0
         end subroutine
      ! --------------------------------------------------------------------------------
 

srcF/sfx_z0t_all_surface.f90

@@ -49,7 +49,7 @@ module sfx_z0t_all_surface
             B = min(B, B_max_land)
        
        z0_t = z0_m / exp(B)
-       write(*,*) 'kl_land', z0_t, B
+       !write(*,*) 'kl_land', z0_t, B
      end subroutine
     ! --------------------------------------------------------------------------------