diff --git a/srcF/sfx_esm.f90 b/srcF/sfx_esm.f90
index c331669f5697a773d2b4df3cc6b7514ca7b5fc54..cb64fa886303aef682bb8a349f960683a20ae8cc 100644
--- a/srcF/sfx_esm.f90
+++ b/srcF/sfx_esm.f90
@@ -27,6 +27,8 @@ module sfx_esm
     ! --------------------------------------------------------------------------------
     public :: get_surface_fluxes
     public :: get_surface_fluxes_vec
+    integer z0m_id
+    integer z0t_id
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------
@@ -139,7 +141,7 @@ contains
             meteo_cell = meteoDataType(&
                     h = meteo%h(i), &
                     U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
-                    z0_m = meteo%z0_m(i))
+                    z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -174,7 +176,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 :: z0_map    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+        real :: depth   
+        real :: lai
         ! ----------------------------------------------------------------------------
 
         ! --- local variables
@@ -233,34 +236,25 @@ contains
         dQ = meteo%dQ
         h = meteo%h
         z0_m = meteo%z0_m
-        z0_map = meteo%z0_m
+        depth = meteo%depth
+        lai = meteo%lai
+        surface_type=meteo%surface_type
 
+        
+        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, z0m_id)
+        
+        call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, z0m_id)
 
-        ! --- 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)
-            call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_map, 1)
-            ! --- 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)
+        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, 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
 
diff --git a/srcF/sfx_log.f90 b/srcF/sfx_log.f90
index ae5a356209472d3a0d5ee12ad2e1f84d78cd4571..8131c5ca960aa99cf88fe3925225d3e24cbc8ed8 100644
--- a/srcF/sfx_log.f90
+++ b/srcF/sfx_log.f90
@@ -23,6 +23,8 @@ module sfx_log
     ! --------------------------------------------------------------------------------
     public :: get_surface_fluxes
     public :: get_surface_fluxes_vec
+    integer z0m_id
+    integer z0t_id
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------
@@ -55,7 +57,7 @@ contains
             meteo_cell = meteoDataType(&
                     h = meteo%h(i), &
                     U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
-                    z0_m = meteo%z0_m(i))
+                    z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -88,6 +90,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
@@ -138,31 +142,25 @@ contains
         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)
+        depth = meteo%depth
+        lai = meteo%lai
+        surface_type=meteo%surface_type
+
+        
+        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, z0m_id)
+        
+        call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, 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, 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
 
diff --git a/srcF/sfx_most.f90 b/srcF/sfx_most.f90
index f4ce053676793567cd6772a1363c09853d50a45f..a98b46f7611c5479816d2c04a68861114085581e 100644
--- a/srcF/sfx_most.f90
+++ b/srcF/sfx_most.f90
@@ -24,6 +24,8 @@ module sfx_most
     public :: get_surface_fluxes
     public :: get_surface_fluxes_vec
     public :: get_psi
+    integer z0m_id
+    integer z0t_id
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------
@@ -56,7 +58,7 @@ contains
             meteo_cell = meteoDataType(&
                     h = meteo%h(i), &
                     U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
-                    z0_m = meteo%z0_m(i))
+                    z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -89,6 +91,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
@@ -141,32 +145,30 @@ contains
         dQ = meteo%dQ
         h = meteo%h
         z0_m = meteo%z0_m
+        depth = meteo%depth
+        lai = meteo%lai
+        surface_type=meteo%surface_type
+
+        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, z0m_id)
+        
+        call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, 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, z0t_id)
+        
+        Re = u_dyn0 * z0_m / nu_air
 
-        ! --- 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
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
+            
             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)
+        
         Re = u_dyn0 * z0_m / nu_air
-        call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
-
-        ! --- define relative height [thermal]
+       
         h0_t = h / z0_t
 
         ! --- define Ri-bulk
diff --git a/srcF/sfx_most_snow.f90 b/srcF/sfx_most_snow.f90
index 0599d14596010788780bb76d342b6d12ddc8d2fe..e5c8c5568e0b7ba6e1bfd31926de9e3e91277440 100644
--- a/srcF/sfx_most_snow.f90
+++ b/srcF/sfx_most_snow.f90
@@ -26,6 +26,7 @@ module sfx_most_snow
     public :: get_surface_fluxes_vec
     public :: get_psi
     integer z0m_id
+    integer z0t_id
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------
@@ -60,7 +61,7 @@ contains
             meteo_cell = meteoDataType(&
                     h = meteo%h(i), &
                     U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
-                    z0_m = meteo%z0_m(i))
+                    z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -93,7 +94,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 :: z0_map    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+        real :: depth  
+        real :: lai
         ! ----------------------------------------------------------------------------
 
         ! --- local variables
@@ -156,22 +158,26 @@ contains
         dT = meteo%dT
         dQ = meteo%dQ
         h = meteo%h
-        !z0_m = meteo%z0_m
-        z0_map = meteo%z0_m
+        z0_m = meteo%z0_m
+        depth = meteo%depth
+        lai = meteo%lai
+        surface_type=meteo%surface_type
 
         
-        call get_dynamic_roughness_definition(1, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
+        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, z0m_id)
         
-        call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_map, z0m_id)
+        call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, 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, z0t_id)
+        
+        Re = u_dyn0 * z0_m / nu_air
         
+        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 thermal roughness & B = log(z0_m / z0_h)
-        Re = u_dyn0 * z0_m / nu_air
-        call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
-
         ! --- define relative height [thermal]
         h0_t = h / z0_t
 
diff --git a/srcF/sfx_sheba.f90 b/srcF/sfx_sheba.f90
index f6160dd1af53e6be98bf42204602cfba142579ff..4e51cb8497d802c6b72368ef47e0b5b228f727ff 100644
--- a/srcF/sfx_sheba.f90
+++ b/srcF/sfx_sheba.f90
@@ -29,6 +29,8 @@ module sfx_sheba
     public :: get_surface_fluxes
     public :: get_surface_fluxes_vec
     public :: get_psi
+    integer z0m_id
+    integer z0t_id
     ! --------------------------------------------------------------------------------
 
     ! --------------------------------------------------------------------------------
@@ -138,7 +140,7 @@ contains
             meteo_cell = meteoDataType(&
                     h = meteo%h(i), &
                     U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
-                    z0_m = meteo%z0_m(i))
+                    z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -171,6 +173,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 :: lai
+        real :: depth
         ! ----------------------------------------------------------------------------
 
         ! --- local variables
@@ -223,32 +227,22 @@ contains
         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)
+        depth = meteo%depth
+        lai = meteo%lai
+        surface_type=meteo%surface_type
+
+        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, z0m_id)
+        
+        call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m, 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, z0t_id)
+        
         Re = u_dyn0 * z0_m / nu_air
-        call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)
-
-        ! --- define relative height [thermal]
+        
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
+        
         h0_t = h / z0_t
 
         ! --- define Ri-bulk
diff --git a/srcF/sfx_surface.f90 b/srcF/sfx_surface.f90
index a7c9b37d21db7d002fdc2222a3e663cfa68817d3..8c6c579b212db7024c522b5ae724221f771ee3b1 100644
--- a/srcF/sfx_surface.f90
+++ b/srcF/sfx_surface.f90
@@ -7,6 +7,7 @@ module sfx_surface
     ! --------------------------------------------------------------------------------
     use sfx_phys_const
     use sfx_z0m_all_surface
+    use sfx_z0t_all_surface
     ! --------------------------------------------------------------------------------
 
     ! directives list
@@ -23,6 +24,8 @@ module sfx_surface
     public :: get_thermal_roughness
     public :: get_dynamic_roughness_all
     public :: get_dynamic_roughness_definition
+    public :: get_thermal_roughness_all
+    public :: get_thermal_roughness_definition
     ! --------------------------------------------------------------------------------
 
 
@@ -49,11 +52,38 @@ module sfx_surface
     integer, public, parameter :: z0m_user = 4 
 
 
+    integer, public, parameter :: z0t_kl_water = 0 
+    integer, public, parameter :: z0t_kl_land = 1   
+    integer, public, parameter :: z0t_re = 2      
+    integer, public, parameter :: z0t_zi = 3         
+    integer, public, parameter :: z0t_ca = 4    
+    integer, public, parameter :: z0t_cz = 5      
+    integer, public, parameter :: z0t_br = 6      
+    integer, public, parameter :: z0t_ot = 7      
+    integer, public, parameter :: z0t_du = 8      
+    integer, public, parameter :: z0t_zm = 9      
+    integer, public, parameter :: z0t_mix = 10      
+    integer, public, parameter :: z0t_user = 11   
+
+
     character(len = 16), parameter :: z0m_tag_ch = 'charnock'
     character(len = 16), parameter :: z0m_tag_fe = 'fetch'
     character(len = 16), parameter :: z0m_tag_ow = 'owen'
     character(len = 16), parameter :: z0m_tag_map = 'map'
     character(len = 16), parameter :: z0m_tag_user = 'user'
+
+    character(len = 16), parameter :: z0t_tag_kl_water = 'kl_water'
+    character(len = 16), parameter :: z0t_tag_kl_land = 'kl_land'
+    character(len = 16), parameter :: z0t_tag_re = 're'
+    character(len = 16), parameter :: z0t_tag_zi = 'zi'
+    character(len = 16), parameter :: z0t_tag_ca = 'ca'
+    character(len = 16), parameter :: z0t_tag_cz = 'cz'
+    character(len = 16), parameter :: z0t_tag_br = 'br'
+    character(len = 16), parameter :: z0t_tag_ot = 'ot'
+    character(len = 16), parameter :: z0t_tag_du = 'du'
+    character(len = 16), parameter :: z0t_tag_zm = 'zm'
+    character(len = 16), parameter :: z0t_tag_mix = 'mix'
+    character(len = 16), parameter :: z0t_tag_user = 'zt_user'
     
 
     integer, public, parameter :: ocean_z0m_id = z0m_ch     !< ocean surface
@@ -63,7 +93,13 @@ module sfx_surface
     integer, public, parameter :: forest_z0m_id = z0m_map     !< forest csurface  
     integer, public, parameter :: usersf_z0m_id = z0m_map       !< user surface  
 
-    real, public, parameter :: depth = 1.0
+    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_re        !< lake surface
+    integer, public, parameter :: snow_z0t_id = z0t_ca        !< snow covered surface    
+    integer, public, parameter :: forest_z0t_id = z0t_du      !< forest csurface  
+    integer, public, parameter :: usersf_z0t_id = z0t_mix       !< user surface  
+    
     ! --------------------------------------------------------------------------------
     real, parameter, private :: kappa = 0.40         !< von Karman constant [n/d]
     ! --------------------------------------------------------------------------------
@@ -81,21 +117,21 @@ module sfx_surface
     ! --------------------------------------------------------------------------------
 
     !< Re fully roughness minimum value [n/d]
-    real, parameter ::  Re_rough_min = 16.3
+    !real, parameter ::  Re_rough_min = 16.3
 
     !< roughness model coeff. [n/d]
     !< --- transitional mode
     !<     B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
-    real, parameter :: B1_rough = 5.0 / 6.0
-    real, parameter :: B2_rough = 0.45
-    real, parameter :: B3_rough = kappa * Pr_m
+    !real, parameter :: B1_rough = 5.0 / 6.0
+    !real, parameter :: B2_rough = 0.45
+    !real, parameter :: B3_rough = kappa * Pr_m
     !< --- fully rough mode (Re > Re_rough_min)
     !<     B = B4 * Re^(B2)
-    real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
+    !real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
 
-    real, parameter :: B_max_land = 2.0
-    real, parameter :: B_max_ocean = 8.0
-    real, parameter :: B_max_lake = 8.0
+    !real, parameter :: B_max_land = 2.0
+    !real, parameter :: B_max_ocean = 8.0
+    !real, parameter :: B_max_lake = 8.0
 
 contains
 
@@ -185,6 +221,74 @@ contains
         end if 
      end function
 
+     function get_surface_z0t_id(tag_z0t) result(z0t_id)
+        implicit none
+        character(len=*), intent(in) :: tag_z0t
+        integer :: z0t_id
+
+        z0t_id = - 1
+        if (trim(tag_z0t) == trim(z0t_tag_kl_water)) then
+            z0t_id = z0t_kl_water
+        else if (trim(tag_z0t) == trim(z0t_tag_kl_land)) then
+            z0t_id = z0t_kl_land
+        else if (trim(tag_z0t) == trim(z0t_tag_re)) then
+            z0t_id = z0t_re
+        else if (trim(tag_z0t) == trim(z0t_tag_zi)) then
+            z0t_id = z0t_zi
+        else if (trim(tag_z0t) == trim(z0t_tag_ca)) then
+            z0t_id = z0t_ca
+        else if (trim(tag_z0t) == trim(z0t_tag_cz)) then
+            z0t_id = z0t_cz
+        else if (trim(tag_z0t) == trim(z0t_tag_br)) then
+            z0t_id = z0t_br
+        else if (trim(tag_z0t) == trim(z0t_tag_ot)) then
+            z0t_id = z0t_ot
+        else if (trim(tag_z0t) == trim(z0t_tag_du)) then
+            z0t_id = z0t_du
+        else if (trim(tag_z0t) == trim(z0t_tag_zm)) then
+            z0t_id = z0t_zm
+        else if (trim(tag_z0t) == trim(z0t_tag_mix)) then
+            z0t_id = z0t_mix
+        else if (trim(tag_z0t) == trim(z0t_tag_user)) then
+            z0t_id = z0t_user
+        end if
+
+    end function
+
+    function get_surface_z0t_tag(z0t_id) result(tag_z0t)
+        implicit none
+        integer :: z0t_id
+        character(len=:), allocatable :: tag_z0t
+
+        tag_z0t = 'undefined'
+        if (z0t_id == z0t_kl_water) then
+            tag_z0t = z0t_tag_kl_water
+        else if (z0t_id == z0t_kl_land) then
+            tag_z0t =  z0t_tag_kl_land
+        else if (z0t_id == z0t_re) then
+            tag_z0t = z0t_tag_re
+        else if (z0t_id == z0t_zi) then
+            tag_z0t = z0t_tag_zi
+        else if (z0t_id == z0t_ca) then
+            tag_z0t = z0t_tag_ca
+        else if (z0t_id == z0t_cz) then
+            tag_z0t = z0t_tag_cz
+        else if (z0t_id == z0t_br) then
+            tag_z0t = z0t_tag_br
+        else if (z0t_id == z0t_ot) then
+            tag_z0t = z0t_tag_ot
+        else if (z0t_id == z0t_du) then
+            tag_z0t = z0t_tag_du
+        else if (z0t_id == z0t_zm) then
+            tag_z0t = z0t_tag_zm
+        else if (z0t_id == z0t_mix) then
+            tag_z0t = z0t_tag_mix
+        else if (z0t_id == z0t_user) then
+            tag_z0t = z0t_tag_user
+        end if 
+     end function
+
+
 
 #if defined(INCLUDE_CXX)
     subroutine set_c_struct_sfx_surface_param_values(surface_param)
@@ -278,7 +382,86 @@ end if
 end subroutine
 ! --------------------------------------------------------------------------------  
 
+! thermal roughness definition
+    ! --------------------------------------------------------------------------------
+subroutine get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
+    forest_z0t_id, usersf_z0t_id, z0t_id)
+! ----------------------------------------------------------------------------
+    integer, intent(out) :: z0t_id              
 
+    integer, intent(in) :: surface_type             
+
+    integer, intent(in) :: ocean_z0t_id            
+    integer, intent(in) :: land_z0t_id             
+    integer, intent(in) :: lake_z0t_id
+    integer, intent(in) :: snow_z0t_id
+    integer, intent(in) :: forest_z0t_id
+    integer, intent(in) :: usersf_z0t_id
+
+! ---------------------------------------------------------------------------
+
+    if (surface_type == surface_ocean) then
+    z0t_id = ocean_z0t_id
+    else if (surface_type == surface_land) then
+    z0t_id = land_z0t_id
+    else if (surface_type == surface_snow) then
+    z0t_id = snow_z0t_id
+    else if (surface_type == surface_lake) then
+    z0t_id = lake_z0t_id
+    else if (surface_type == surface_forest) then
+    z0t_id = forest_z0t_id
+    else if (surface_type == surface_user) then
+    z0t_id  = usersf_z0t_id
+    end if 
+
+    end subroutine
+
+!---------------------------------------------------------------
+    subroutine get_thermal_roughness_all(z0_t, B, &
+    z0_m, Re, u_dyn, LAI, z0t_id)
+! ----------------------------------------------------------------------------
+    real, intent(out) :: z0_t               !< thermal roughness [m]
+    real, intent(out) :: B                  !< = log(z0_m / z0_t) [n/d]
+
+    real, intent(in) :: z0_m                !< aerodynamic roughness [m]
+    real, intent(in) :: Re                  !< roughness Reynolds number [n/d]
+    real, intent(in) :: LAI                 !< leaf-area index
+    real, intent(in) :: u_dyn              !< dynamic velocity [m/s]
+    
+
+    integer, intent(in) :: z0t_id
+! ---------------------------------------------------------------------------
+    real  Czm                 !< proportionality coefficient z0_t =Czm*z0_m
+    Czm=1
+
+    if (z0t_id == z0t_kl_water) then
+    call get_thermal_roughness_kl_water(z0_t, B, z0_m, Re)
+    else if (z0t_id == z0t_kl_land) then
+    call get_thermal_roughness_kl_land(z0_t, B, z0_m, Re)
+    else if (z0t_id == z0t_re) then
+    call get_thermal_roughness_re(z0_t, B, z0_m, Re)
+    else if (z0t_id == z0t_zi) then
+    call get_thermal_roughness_zi(z0_t, B, z0_m, Re)
+    else if (z0t_id == z0t_ca) then
+    call get_thermal_roughness_ca(z0_t, B, z0_m, Re)
+    else if (z0t_id == z0t_cz) then
+    call get_thermal_roughness_cz(z0_t, B, z0_m, Re)
+    else if (z0t_id == z0t_br) then
+    call get_thermal_roughness_br(z0_t, B, z0_m, Re)
+    else if (z0t_id == z0t_ot) then
+    call get_thermal_roughness_ot(z0_t, B, z0_m, Re)
+    else if (z0t_id == z0t_du) then
+    call get_thermal_roughness_du(z0_t, B, z0_m, u_dyn, LAI)
+    else if (z0t_id == z0t_zm) then
+    call get_thermal_roughness_zm(z0_t, B, z0_m, Czm)
+    else if (z0t_id == z0t_mix) then
+    call get_thermal_roughness_mix(z0_t, B, z0_m, u_dyn, Re)
+    else if (z0t_id == z0t_user) then
+    write(*, *) 'z0t_user'
+    end if
+
+
+end subroutine
 
 
 
@@ -355,7 +538,7 @@ end subroutine
         if (surface_type == surface_ocean) then
             B = min(B, B_max_ocean)
         else if (surface_type == surface_lake) then
-            B = min(B, B_max_lake)
+            B = min(B, B_max_ocean)
         else if (surface_type == surface_land) then
             B = min(B, B_max_land)
         end if
diff --git a/srcF/sfx_z0t_all_surface.f90 b/srcF/sfx_z0t_all_surface.f90
index 6670babb1f6bcc4d9fb0a3f9467b6b4b48de62e0..cee7517a5afcb6d663b7a77769798c0cb5803b40 100644
--- a/srcF/sfx_z0t_all_surface.f90
+++ b/srcF/sfx_z0t_all_surface.f90
@@ -10,18 +10,18 @@ module sfx_z0t_all_surface
     real, parameter, private :: kappa = 0.40         !< von Karman constant [n/d]
     real, parameter, private :: Pr_m = 0.71              !< molecular Prandtl number (air) [n/d]
     !< Re fully roughness minimum value [n/d]
-    real, parameter ::  Re_rough_min = 16.3
+    real, public, parameter ::  Re_rough_min = 16.3
     !< roughness model coeff. [n/d]
     !< --- transitional mode
     !<     B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
-    real, parameter :: B1_rough = 5.0 / 6.0
-    real, parameter :: B2_rough = 0.45
+    real, public, parameter :: B1_rough = 5.0 / 6.0
+    real, public, parameter :: B2_rough = 0.45
     real, parameter :: B3_rough = kappa * Pr_m
     !< --- fully rough mode (Re > Re_rough_min)
     !<     B = B4 * Re^(B2)
-    real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
-     real, parameter :: B_max_ocean = 8.0
-    real, parameter :: B_max_land = 2.0
+    real, public, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
+    real, public, parameter :: B_max_ocean = 8.0
+    real, public, parameter :: B_max_land = 2.0
     
 
     contains