diff --git a/srcF/sfx_sheba_noit_param.f90 b/srcF/sfx_sheba_noit_param.f90
index 6461dde50e1fa347b7f1cfe86934cdc22b3be970..76807354a96aa7cd4383691edb72474b042d927f 100644
--- a/srcF/sfx_sheba_noit_param.f90
+++ b/srcF/sfx_sheba_noit_param.f90
@@ -33,8 +33,17 @@ module sfx_sheba_noit_param
     real, parameter :: a_h = 5.0
     real, parameter :: b_h = 5.0
     real, parameter :: c_h = 3.0
-
+    
     real, parameter  ::  gamma = 2.91, zeta_a = 3.6 ! for stable psi
 
+        	!< snow parameters
+    real, parameter :: rho_s =900
+    real, parameter :: d_s=0.0000886
+    real, parameter :: h_salt_const=0.05
+    real, parameter :: u_thsnow=0.35   
+    real, parameter :: S_salt_const=0.000004   
+    real, parameter :: rho_air=1.2
+    real, parameter :: Csn=3.25
+
 
 end module sfx_sheba_noit_param
diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90
index d40165edff3988ee3b2b5a69a4b8ed6066740326..35e3b7ca50b6bac1b5d51bb2343d01e1f291f9fb 100644
--- a/srcF/sfx_sheba_noniterative.f90
+++ b/srcF/sfx_sheba_noniterative.f90
@@ -28,6 +28,8 @@ module sfx_sheba_noniterative
     public :: get_surface_fluxes
     public :: get_surface_fluxes_vec
     public :: get_psi_stable
+    integer z0m_id
+    integer z0t_id
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------
@@ -184,6 +186,8 @@ contains
         real :: Tsemi   !< semi-sum of potential temperature at 'h' and at surface [K]
         real :: dQ      !< difference between humidity at 'h' and at surface [g/g]
         real :: z0_m    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+        real :: depth   
+        real :: lai
         ! ----------------------------------------------------------------------------
 
         ! --- local variables
@@ -193,10 +197,11 @@ contains
         real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
 
         real u_dyn0             !< dynamic velocity in neutral conditions [m/s]
-        real Re                 !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
+        real Rib, Ri_sn         !< bulk Richardson number, snow Richardson number
+
 
         real zeta               !< = z/L [n/d]
-        real Rib                !< bulk Richardson number
+        real Re                 !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
 
         real zeta_conv_lim      !< z/L critical value for matching free convection limit [n/d]
         real Rib_conv_lim       !< Ri-bulk critical value for matching free convection limit [n/d]
@@ -220,7 +225,10 @@ contains
 
         real fval               !< just a shortcut for partial calculations
 
-
+        real w_snow, deltaS, h_salt
+        real sigma_w, sigma_r
+        real S_mean, S_salt
+        real z0_m1
 #ifdef SFX_CHECK_NAN
         real NaN
 #endif
@@ -246,50 +254,65 @@ contains
         dT = meteo%dT
         dQ = meteo%dQ
         h = meteo%h
-        z0_m = meteo%z0_m
-
-        ! --- define surface type
-        if (z0_m < 0.0) then
-            surface_type = surface_ocean
-        else
-            surface_type = surface_land
-        end if
-
-        if (surface_type == surface_ocean) then
-            ! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
-            call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
-            ! --- define relative height
-            h0_m = h / z0_m
-        endif
-        if (surface_type == surface_land) then
-            ! --- define relative height
-            h0_m = h / z0_m
-            ! --- define dynamic velocity in neutral conditions
-            u_dyn0 = U * kappa / log(h0_m)
-        end if
-
-        ! --- define thermal roughness & B = log(z0_m / z0_h)
+        z0_m1 = meteo%z0_m
+        depth = meteo%depth
+        lai = meteo%lai
+        surface_type=meteo%surface_type
+       
+        
+        !write (*,*) surface_type, 'esm'
+        call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
+        forest_z0m_id, usersf_z0m_id, ice_z0m_id, z0m_id)
+        
+        
+        call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
+        
+        call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
+        forest_z0t_id, usersf_z0t_id, ice_z0t_id, z0t_id)
+        
+        
         Re = u_dyn0 * z0_m / nu_air
-        call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
 
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
+        ! --- define relative height
+            h0_m = h / z0_m
         ! --- define relative height [thermal]
         h0_t = h / z0_t
-
+      
+       
+        
+        !  --- 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)
+                      
+        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
         ! ----------------------------------------------------------------------------
-        if (Rib > 0.0) then
+        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)
+           
+        if (surface_type == surface_snow) then    
+            write(*,*) 'snow', Ri_sn, Rib
+            Rib=Rib+Ri_sn 
+        endif
 
             call get_zeta(zeta, Rib, h, z0_m, z0_t)
 
@@ -357,6 +380,63 @@ contains
 
     end subroutine get_surface_fluxes
     ! --------------------------------------------------------------------------------
+    !--------------------------------------snow functions-----------------------------
+    subroutine get_sigma(sigma_r,  sigma_w, rho_air, rho_s)
+        !< @brief function for 
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: sigma_r, sigma_w  !< s
+        real, intent(in) ::  rho_air, rho_s            
+        ! ----------------------------------------------------------------------------
+        sigma_r  = ((rho_s/rho_air)-1)
+        sigma_w  = (rho_s-rho_air)/rho_air
+    end subroutine
+
+    subroutine get_h_salt(h_salt, u_dyn0)
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: h_salt   
+        real, intent(in) :: u_dyn0
+        ! ----------------------------------------------------------------------------
+        h_salt  = 0.08436*u_dyn0**1.27
+    end subroutine
+
+    
+    subroutine get_S_mean(S_mean, S_salt, h_salt, z, w_snow, u_dyn0)
+        !< @brief function for snow consentration
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: S_mean   !< snow consentration
+        real, intent(in) ::  S_salt, h_salt, z, w_snow, u_dyn0         
+        ! ----------------------------------------------------------------------------
+        S_mean  = S_salt *  (z/h_salt)**(-w_snow/(kappa*u_dyn0))
+    end subroutine
+
+    subroutine get_S_salt(S_salt, u_dyn0, u_thsnow, g, h_salt)
+        !< @brief function for snow consentration
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: S_salt   !< snow consentration
+        real, intent(in) ::  u_dyn0,u_thsnow, g, h_salt         
+        ! ----------------------------------------------------------------------------
+        real qs
+        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)
+        else
+            S_salt=0.0
+        endif
+    end subroutine
+
+    subroutine get_w_snow(w_snow, sigma_w, g, d_s, nu_air)
+        !< @brief function for smow velosity
+        ! ----------------------------------------------------------------------------
+        real, intent(out) :: w_snow   !< 
+        real, intent(in) ::  sigma_w, g, d_s, nu_air            
+        ! ----------------------------------------------------------------------------
+        w_snow  = (sigma_w * g * d_s * d_s) / (18.0 * nu_air);
+
+    end subroutine
+         !-----------------------------------------------------------------------------
+
+
+
 
     subroutine get_zeta(zeta, Rib, h, z0_m, z0_t)
     real,intent(out) :: zeta
diff --git a/srcF/sfx_surface.f90 b/srcF/sfx_surface.f90
index e356c31688fba2c09ad0cc8909b3c0c22f1c766c..70b10b2379513ee729f4612603236e7feea0ffd1 100644
--- a/srcF/sfx_surface.f90
+++ b/srcF/sfx_surface.f90
@@ -20,8 +20,8 @@ module sfx_surface
 #if defined(INCLUDE_CXX)
     public :: set_c_struct_sfx_surface_param_values
 #endif
-    public :: get_charnock_roughness
-    public :: get_thermal_roughness
+    !public :: get_charnock_roughness
+    !public :: get_thermal_roughness
     public :: get_dynamic_roughness_all
     public :: get_dynamic_roughness_definition
     public :: get_thermal_roughness_all
@@ -96,22 +96,22 @@ module sfx_surface
     character(len = 16), parameter :: z0t_tag_user = 'zt_user'
     
 
-    integer, public, parameter :: ocean_z0m_id = z0m_fe     !< ocean surface
+    integer, public, parameter :: ocean_z0m_id = z0m_ch     !< ocean surface
     integer, public, parameter :: land_z0m_id = z0m_map_id        !< land surface
-    integer, public, parameter :: lake_z0m_id = z0m_fe        !< lake surface
-    integer, public, parameter :: snow_z0m_id = z0m_and        !< snow covered surface    
+    integer, public, parameter :: lake_z0m_id = z0m_map_id        !< lake surface
+    integer, public, parameter :: snow_z0m_id = z0m_map_id       !< snow covered surface    
     integer, public, parameter :: forest_z0m_id = z0m_map_id     !< forest csurface  
-    integer, public, parameter :: usersf_z0m_id = z0m_ch       !< user surface  
+    integer, public, parameter :: usersf_z0m_id = z0m_map_id      !< user surface  
     integer, public, parameter :: ice_z0m_id = z0m_map_id       !< ice surface  
 
 
-    integer, public, parameter :: ocean_z0t_id = z0t_ca     !< ocean surface
-    integer, public, parameter :: land_z0t_id = z0t_ot        !< land surface
-    integer, public, parameter :: lake_z0t_id = z0t_ot        !< lake surface
-    integer, public, parameter :: snow_z0t_id = z0t_zi         !< snow covered surface    
-    integer, public, parameter :: forest_z0t_id = z0t_ot      !< forest csurface  
-    integer, public, parameter :: usersf_z0t_id = z0t_kl_water      !< user surface  
-    integer, public, parameter :: ice_z0t_id = z0t_ca     !< user surface  
+    integer, public, parameter :: ocean_z0t_id = z0t_kl_water     !< ocean surface
+    integer, public, parameter :: land_z0t_id = z0t_kl_land        !< land surface
+    integer, public, parameter :: lake_z0t_id = z0t_kl_land        !< lake surface
+    integer, public, parameter :: snow_z0t_id = z0t_kl_land         !< snow covered surface    
+    integer, public, parameter :: forest_z0t_id = z0t_kl_land      !< forest csurface  
+    integer, public, parameter :: usersf_z0t_id = z0t_kl_land      !< user surface  
+    integer, public, parameter :: ice_z0t_id = z0t_kl_land     !< user surface  
     
     ! --------------------------------------------------------------------------------
     real, parameter, private :: kappa = 0.40         !< von Karman constant [n/d]
diff --git a/srcF/sfx_z0m_all_surface.f90 b/srcF/sfx_z0m_all_surface.f90
index 29185386a7c3a7f204f4094aab59ee9ec59cd9f7..593f663250ee0905416efcf5129116e14cbad955 100644
--- a/srcF/sfx_z0m_all_surface.f90
+++ b/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
     end subroutine
     ! --------------------------------------------------------------------------------
     
@@ -125,7 +125,7 @@ module sfx_z0m_all_surface
         
 
         u_dyn0 = Uc / c
-        !  write(*,*) z0_m, 'oewn'
+        write(*,*) 'ow', 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
         end subroutine
      ! --------------------------------------------------------------------------------
 
diff --git a/srcF/sfx_z0t_all_surface.f90 b/srcF/sfx_z0t_all_surface.f90
index 8b9eb71d96b67fce85b188d1a2610362adf3ca80..c6581f3b82f83bed460de75bfed918af5a7be352 100644
--- a/srcF/sfx_z0t_all_surface.f90
+++ b/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
      end subroutine
     ! --------------------------------------------------------------------------------  
         
@@ -75,7 +75,7 @@ module sfx_z0t_all_surface
             B = min(B, B_max_ocean)
        
        z0_t = z0_m / exp(B)
-
+       write(*,*) 'kl_water', z0_t, B
      end subroutine