diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4c4486280e7b79d953059ce30274d1cbf7efe1e9..ea81f2abc374b714f09ec1ad3ad0840f61b7db77 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -83,10 +83,15 @@ set(SOURCES_F
     srcF/sfx_log_param.f90
     srcF/sfx_run.f90
     srcF/sfx_phys_const.f90
+    srcF/sfx_z0m_all_surface
+    srcF/sfx_z0t_all_surface
+    srcF/sfx_z0m_all
+    srcF/sfx_z0t_all
     srcF/sfx_surface.f90
-    srcF/sfx_thermal_roughness.f90
     srcF/sfx_most.f90
     srcF/sfx_most_param.f90
+    srcF/sfx_most_snow.f90
+    srcF/sfx_most_snow_param.f90
     srcF/sfx_sheba.f90
     srcF/sfx_sheba_noniterative.f90
     srcF/sfx_sheba_param.f90
diff --git a/srcF/sfx_config.f90 b/srcF/sfx_config.f90
index 3fde2a0995da10b521591349c02afe7ba8e98d57..0528a049482fdfca385dd99767857f61bbfe02a2 100644
--- a/srcF/sfx_config.f90
+++ b/srcF/sfx_config.f90
@@ -21,12 +21,16 @@ module sfx_config
     integer, parameter :: model_sheba = 3           !< SHEBA model
     integer, parameter :: model_sheba_noit = 4           !< SHEBA model noniterative
 
+    integer, parameter :: model_most_snow = 5       !< MOST_SNOW model
+    
     character(len = 16), parameter :: model_esm_tag = 'esm'
     character(len = 16), parameter :: model_log_tag = 'log'
     character(len = 16), parameter :: model_most_tag = 'most'
     character(len = 16), parameter :: model_sheba_tag = 'sheba'
     character(len = 16), parameter :: model_sheba_noit_tag = 'sheba_noit'
 
+    character(len = 16), parameter :: model_most_snow_tag = 'most_snow'
+    
     !> @brief dataset enum def.
     integer, parameter :: dataset_mosaic = 1        !< MOSAiC campaign
     integer, parameter :: dataset_irgason = 2       !< IRGASON data
@@ -72,7 +76,9 @@ contains
             id = model_sheba
         else if (trim(tag) == trim(model_sheba_noit_tag)) then
             id = model_sheba_noit
-        end if
+       else if (trim(tag) == trim(model_most_snow_tag)) then
+            id = model_most_snow 
+       end if
 
     end function
 
@@ -92,6 +98,8 @@ contains
             tag = model_sheba_tag
         else if (id == model_sheba_noit) then
             tag = model_sheba_noit_tag
+        else if (id == model_most_snow) then
+            tag = model_most_snow_tag
         end if 
 
     end function
diff --git a/srcF/sfx_run.f90 b/srcF/sfx_run.f90
index eec3f6ea7eb0196907298ea8892b6a21113a3482..7fd4880e2e686c9d61763d8ff4ea6b06ad442b8b 100644
--- a/srcF/sfx_run.f90
+++ b/srcF/sfx_run.f90
@@ -43,6 +43,9 @@ contains
         use sfx_sheba_noniterative, only: &
                 get_surface_fluxes_vec_sheba_noit => get_surface_fluxes_vec, &
                 numericsType_sheba_noit => numericsType
+        use sfx_most_snow, only: &
+                get_surface_fluxes_vec_most_snow => get_surface_fluxes_vec, &
+                numericsType_most_snow => numericsType
         ! --------------------------------------------------------------------------------
     
         character(len=*), intent(in) :: filename_out
@@ -63,7 +66,8 @@ contains
         type(numericsType_most) :: numerics_most    !< surface flux module (MOST) numerics parameters
         type(numericsType_sheba) :: numerics_sheba  !< surface flux module (SHEBA) numerics parameters
         type(numericsType_sheba_noit) :: numerics_sheba_noit  !< surface flux module (SHEBA) numerics parameters
-    
+        type(numericsType_most_snow) :: numerics_most_snow    !< surface flux module (MOST_SNOW) numerics parameters
+
         integer :: num                          !< number of 'cells' in input
         ! --------------------------------------------------------------------------------
     
@@ -153,6 +157,8 @@ contains
             call get_surface_fluxes_vec_sheba(sfx, meteo, numerics_sheba, num)
         else if (model == model_sheba_noit) then
             call get_surface_fluxes_vec_sheba_noit(sfx, meteo, numerics_sheba_noit, num)
+        else if (model == model_most_snow) then
+            call get_surface_fluxes_vec_most_snow(sfx, meteo, numerics_most_snow, num)
         end if
     
     
@@ -239,7 +245,7 @@ contains
                 write(*, *) ' --help'
                 write(*, *) '    print usage options'
                 write(*, *) ' --model [key]'
-                write(*, *) '    key = esm (default) || log || most || sheba || sheba_noit'
+                write(*, *) '    key = esm (default) || log || most || sheba || sheba_noit || most_snow'
                 write(*, *) ' --dataset [key]'
                 write(*, *) '    key = mosaic (default) || irgason || sheba'
                 write(*, *) '        = lake || papa || toga || user [filename] [param]'
diff --git a/srcF/sfx_z0m_all.f90 b/srcF/sfx_z0m_all.f90
index 89fcaf3df9d8a4e24fc115342e1821309fd9c2b7..c325e4d427b0e5d9b059731750bdd33c083e82bd 100644
--- a/srcF/sfx_z0m_all.f90
+++ b/srcF/sfx_z0m_all.f90
@@ -2,7 +2,7 @@ module sfx_z0m_all
    
   !< @brief surface thermal roughness parameterizations for all type surface
 
-    use z0m_all_surface
+    use sfx_z0m_all_surface
     
 
     implicit none
diff --git a/srcF/sfx_z0t_all.f90 b/srcF/sfx_z0t_all.f90
index 8f615b145bb998a1455afbc62914041dea867202..67958e5209a2facc883c8343865885f8069aa075 100644
--- a/srcF/sfx_z0t_all.f90
+++ b/srcF/sfx_z0t_all.f90
@@ -2,7 +2,7 @@ module sfx_z0t_all
    
   !< @brief surface thermal roughness parameterizations for all type surface
 
-    use z0t_all_surface
+    use sfx_z0t_all_surface
     
 
     implicit none
diff --git a/srcF/z0m_all_surface.f90 b/srcF/z0m_all_surface.f90
deleted file mode 100644
index 46262d6e5798afda7f596b5607a1ab852cbc2f70..0000000000000000000000000000000000000000
--- a/srcF/z0m_all_surface.f90
+++ /dev/null
@@ -1,207 +0,0 @@
-module z0m_all_surface
-   !< @brief surface roughness parameterizations 
-
-    use sfx_phys_const
-    implicit none
-    public 
-    ! --------------------------------------------------------------------------------
-
-
-    
-
-    ! --------------------------------------------------------------------------------
-    real, parameter, private :: kappa = 0.40         !< von Karman constant [n/d]
-    ! --------------------------------------------------------------------------------
-
-    !< Charnock parameters
-    !<     z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
-    ! --------------------------------------------------------------------------------
-    
-    real, parameter :: gamma_c = 0.0144
-    real, parameter :: Re_visc_min = 0.111
-
-    real, parameter :: h_charnock = 10.0
-    real, parameter :: c1_charnock = log(h_charnock * (g / gamma_c))
-    real, parameter :: c2_charnock = Re_visc_min * nu_air * c1_charnock
-    real, parameter :: gamma_min = 0.01
-    real, parameter :: gamma_max = 0.11
-    real, parameter :: f_c = 100
-    real, parameter :: eps = 1
-    ! --------------------------------------------------------------------------------
-
-    
-
-   
- contains
-
-    ! charnock roughness definition
-    ! --------------------------------------------------------------------------------
-    subroutine get_dynamic_roughness_ch(z0_m, u_dyn0, U, h, maxiters)
-        ! ----------------------------------------------------------------------------
-        real, intent(out) :: z0_m           !< aerodynamic roughness [m]
-        real, intent(out) :: u_dyn0         !< dynamic velocity in neutral conditions [m/s]
-
-        real, intent(in) :: h               !< constant flux layer height [m]
-        real, intent(in) :: U               !< abs(wind speed) [m/s]
-        integer, intent(in) :: maxiters     !< maximum number of iterations
-        ! ----------------------------------------------------------------------------
-
-        ! --- local variables
-        real :: Uc                  ! wind speed at h_charnock [m/s]
-
-        real :: a, b, c, c_min
-        real :: f
-
-        integer :: i, j
-        ! ----------------------------------------------------------------------------
-
-        Uc = U
-        a = 0.0
-        b = 25.0
-        c_min = log(h_charnock) / kappa
-
-        do i = 1, maxiters
-            f = c1_charnock - 2.0 * log(Uc)
-            do j = 1, maxiters
-                c = (f + 2.0 * log(b)) / kappa
-                ! looks like the check should use U10 instead of U
-                !   but note that a1 should be set = 0 in cycle beforehand
-                if (U <= 8.0e0) a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
-                c = max(c - a, c_min)
-                b = c
-            end do
-            z0_m = h_charnock * exp(-c * kappa)
-            z0_m = max(z0_m, 0.000015e0)
-            Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
-        end do
-
-        ! --- define dynamic velocity in neutral conditions
-        u_dyn0 = Uc / c
-
-    end subroutine
-    ! --------------------------------------------------------------------------------
-    
-      subroutine get_dynamic_roughness_ow(z0_m, u_dyn0, U, h, maxiters)
-      !Owen 1964
-        ! ----------------------------------------------------------------------------
-        real, intent(out) :: z0_m           !< aerodynamic roughness [m]
-        real, intent(out) :: u_dyn0         !< dynamic velocity in neutral conditions [m/s]
-
-        real, intent(in) :: h               !< constant flux layer height [m]
-        real, intent(in) :: U               !< abs(wind speed) [m/s]
-        integer, intent(in) :: maxiters     !< maximum number of iterations
-        ! ----------------------------------------------------------------------------
-
-        ! --- local variables
-        real :: Uc                  ! wind speed at h_charnock [m/s]
-
-        real :: b1, b2, Cuz, betta_u, nu_m, C_z0,c
-        real :: f
-
-        integer :: i, j
-        ! ----------------------------------------------------------------------------
-        Uc=U
-        C_z0=0.007
-        betta_u=0.111
-        nu_m=0.0000133
-        b1=log(h*g/C_z0)
-        b2=betta_u*nu_m*g/C_z0
-        Cuz=25.0
-        do i = 1, maxiters
-        f = c1_charnock - 2.0 * log(Uc)
-        c = (f + 2.0 * log(Cuz)) / kappa
-        Cuz=(1.0/kappa)*(b1+log(U/Cuz)-log(b2+(U/Cuz)*(U/Cuz)))
-        if(Cuz==0.0) exit
-        z0_m=h*exp(-kappa*Cuz)
-        end do
-        
-
-        u_dyn0 = Uc / c
-
-        end subroutine
-     ! --------------------------------------------------------------------------------
-   
-        subroutine get_dynamic_roughness_fetch(z0_m, u_dyn0, U, depth, h, maxiters)
-        ! ----------------------------------------------------------------------------
-        real, intent(out) :: z0_m           !< aerodynamic roughness [m]
-        real, intent(out) :: u_dyn0         !< dynamic velocity in neutral conditions [m/s]
-
-        real, intent(in) :: U               !< abs(wind speed) [m/s]
-        real, intent(in) :: depth           !< depth [m]
-        real, intent(in) :: h               !< constant flux layer height [m]
-        integer, intent(in) :: maxiters     !< maximum number of iterations
-        ! ----------------------------------------------------------------------------
-
-        ! --- local variables
-        real :: Uc                  ! wind speed at h_charnock [m/s]
-
-        real :: a, b, c, c_min
-        real :: f
-
-        real :: A_lake, B_lake, gamma_c, fetch, c1_charnock_lake, c2_charnock_lake
-
-        integer :: i, j
-        ! ----------------------------------------------------------------------------
-
-        Uc = U
-        a = 0.0
-        b = 25.0
-        c_min = log(h_charnock) / kappa
-        
-        fetch = 25.0 * depth !25.0 * depth
-
-        !<     z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
-        !<     gamma_c = gamma_min + (gamma_max - gamma_min) * exp(-min(A_lake, B_lake))
-        !<     А_lake = (fetch * g / U^2)^(1/3) / f_c
-        !<     B_lake = eps (sqrt(depth * g)/U)
-
-        do i = 1, maxiters
-            A_lake = ((fetch * g / (U)**2)**(1/3)) / f_c
-            B_lake = eps * (sqrt(depth * g)/U)
-            gamma_c =  gamma_min + (gamma_max - gamma_min) * exp(-min(A_lake, B_lake))
-            !write(*,*) A_lake
-            !write(*,*) B_lake
-            c1_charnock_lake = log(h_charnock * (g / gamma_c))
-            c2_charnock_lake = Re_visc_min * nu_air * c1_charnock_lake
-            f = c1_charnock_lake - 2.0 * log(Uc)
-            do j = 1, maxiters
-                c = (f + 2.0 * log(b)) / kappa
-                if (U <= 8.0e0) a = log(1.0 + c2_charnock_lake * ((b / Uc)**3)) / kappa
-                c = max(c - a, c_min)
-                b = c
-            end do
-            z0_m = h_charnock * exp(-c * kappa)
-            z0_m = max(z0_m, 0.000015e0)
-            Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
-        end do
-
-        ! --- define dynamic velocity in neutral conditions
-        u_dyn0 = Uc / c
-
-     end subroutine
-   
-   
-subroutine get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map)
-      
-        ! ----------------------------------------------------------------------------
-        real, intent(out) :: z0_m           !< aerodynamic roughness [m]
-        real, intent(out) :: u_dyn0         !< dynamic velocity in neutral conditions [m/s]
-
-        real, intent(in) :: h               !< constant flux layer height [m]
-        real, intent(in) :: z0m_map           !< aerodynamic roughness from map[m]
-        real, intent(in) :: U               !< abs(wind speed) [m/s]
-        ! ----------------------------------------------------------------------------
-        real :: h0_m
-        
-        z0_m=z0m_map  
-        h0_m = h / z0_m
-        u_dyn0 = U * kappa / log(h0_m)
-
-        end subroutine
-     ! --------------------------------------------------------------------------------
-
-
-
-
-   
-    end module  z0m_all_surface
\ No newline at end of file
diff --git a/srcF/z0t_all_surface.f90 b/srcF/z0t_all_surface.f90
deleted file mode 100644
index ca3a8b4f8b1d9baa4a84c1c7698b2313e82a9e78..0000000000000000000000000000000000000000
--- a/srcF/z0t_all_surface.f90
+++ /dev/null
@@ -1,285 +0,0 @@
-module z0t_all_surface
-    !< @brief surface thermal roughness parameterizations 
-
-    implicit none
-    public
-   
-    
-
-    ! --------------------------------------------------------------------------------
-    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
-    !< 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
-    !< --- 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
-    
-
-    contains
-
-    
-     ! thermal roughness definition by Kazakov, Lykosov
-     ! --------------------------------------------------------------------------------
-       subroutine get_thermal_roughness_kl_land(z0_t, B, &
-            z0_m, Re)
-        ! ----------------------------------------------------------------------------
-        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]
-        ! ----------------------------------------------------------------------------
-
-         !--- define B = log(z0_m / z0_t)
-        if (Re <= Re_rough_min) then
-            B = B1_rough * alog(B3_rough * Re) + B2_rough
-        else
-            ! *: B4 takes into account Re value at z' ~ O(10) z0
-            B = B4_rough * (Re**B2_rough)
-        end if
-         
-            B = min(B, B_max_land)
-       
-       z0_t = z0_m / exp(B)
-
-     end subroutine
-    ! --------------------------------------------------------------------------------  
-        
-! --------------------------------------------------------------------------------
-       subroutine get_thermal_roughness_kl_water(z0_t, B, &
-            z0_m, Re)
-        ! ----------------------------------------------------------------------------
-        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]
-        ! ----------------------------------------------------------------------------
-
-         !--- define B = log(z0_m / z0_t)
-        if (Re <= Re_rough_min) then
-            B = B1_rough * alog(B3_rough * Re) + B2_rough
-        else
-            ! *: B4 takes into account Re value at z' ~ O(10) z0
-            B = B4_rough * (Re**B2_rough)
-        end if
-         
-            B = min(B, B_max_ocean)
-       
-       z0_t = z0_m / exp(B)
-
-     end subroutine
-        
-
-    ! thermal roughness definition by Chen, F., Zhang, Y., 2009. 
-    ! --------------------------------------------------------------------------------  
-    subroutine get_thermal_roughness_cz(z0_t, B, &
-            z0_m, Re)
-        ! ----------------------------------------------------------------------------
-        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]
-       
-        
-       
-        B=(kappa*10.0**(-0.4*z0_m/0.07))*(Re**0.45)        !Chen and Zhang
-        
-         ! --- define roughness [thermal]
-        z0_t = z0_m / exp(B)
-
-    end subroutine
-      ! --------------------------------------------------------------------------------
-        ! thermal roughness definition by Zilitinkevich, S., 1995.
-        ! --------------------------------------------------------------------------------
-        subroutine get_thermal_roughness_zi(z0_t, B, &
-            z0_m, Re)
-        ! ----------------------------------------------------------------------------
-        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]
-       
-       
-        
-       
-        B=0.1*kappa*(Re**0.5)         !6-Zilitinkevich
-       
-       
-
-        ! --- define roughness [thermal]
-        z0_t = z0_m / exp(B)
-
-    end subroutine
-    ! --------------------------------------------------------------------------------
-     ! thermal roughness definition by Cahill, A.T., Parlange, M.B., Albertson, J.D., 1997.
-     ! It is better to use for dynamic surfaces such as sand
-    ! --------------------------------------------------------------------------------
-        subroutine get_thermal_roughness_ca(z0_t, B, &
-            z0_m, Re)
-    ! ----------------------------------------------------------------------------
-        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]
-       
-       
-        
-        
-        B=2.46*(Re**0.25)-3.8         !4-Cahill et al. 
-       
-        ! --- define roughness [thermal]
-        z0_t = z0_m / exp(B)
-
-    end subroutine
-
-
-    ! --------------------------------------------------------------------------------
-     ! thermal roughness definition by Brutsaert W., 2003.
-    ! --------------------------------------------------------------------------------
-        subroutine get_thermal_roughness_br(z0_t, B, &
-            z0_m, Re)
-    ! ----------------------------------------------------------------------------
-        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]
-       
-       
-        
-        
-        B=2.46*(Re**0.25)-2.0       !Brutsaert
-       
-        ! --- define roughness [thermal]
-        z0_t = z0_m / exp(B)
-
-    end subroutine
-
-
-
-    ! --------------------------------------------------------------------------------
-     ! thermal roughness definition by Owen P. R., Thomson W. R., 1963. 
-    ! --------------------------------------------------------------------------------
-        subroutine get_thermal_roughness_ot(z0_t, B, &
-            z0_m, Re)
-    ! ----------------------------------------------------------------------------
-        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]
-    
-        
-        
-        B=kappa*(Re**0.45)         !Owen P. R., Thomson W. R. 
-       
-        ! --- define roughness [thermal]
-        z0_t = z0_m / exp(B)
-
-    end subroutine
-    ! --------------------------------------------------------------------------------
-    ! thermal roughness definition by Duynkerke P. G., 1992. 
-    !It is better to use for surfaces wiht forest
-    ! --------------------------------------------------------------------------------
-        subroutine get_thermal_roughness_du(z0_t, B, &
-            z0_m, u_dyn, LAI)
-    ! ----------------------------------------------------------------------------
-        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) :: u_dyn              !< dynamic velocity [m/s]
-        real, intent(in) :: LAI                 !< leaf-area index
-       
-       
-        
-        
-        B=(13*u_dyn**0.4)/LAI+0.85         !Duynkerke P. G., 1992. 
-       
-        ! --- define roughness [thermal]
-        z0_t = z0_m / exp(B)
-
-    end subroutine
-    ! --------------------------------------------------------------------------------
-    ! thermal roughness definition z0_t = C*z0_m  
-    ! --------------------------------------------------------------------------------
-        subroutine get_thermal_roughness_zm(z0_t, B, &
-            z0_m, Czm)
-    ! ----------------------------------------------------------------------------
-        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) :: Czm                 !< proportionality coefficient
-        
-       
-
-        z0_t =Czm*z0_m
-        B=log(z0_m / z0_t)
-    end subroutine
-    ! --------------------------------------------------------------------------------
-    ! thermal roughness definition by  Chen and Zhang  and Zilitinkevich
-    ! --------------------------------------------------------------------------------
-        subroutine get_thermal_roughness_mix(z0_t, B, &
-            z0_m, u_dyn, Re)
-    ! ----------------------------------------------------------------------------
-        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) :: u_dyn              !< dynamic velocity [m/s]
-        real, intent(in) :: Re                  !< roughness Reynolds number [n/d]
-    
-       
-        real, parameter :: u_dyn_th=0.17         !< dynamic velocity treshhold [m/s]
-        
-        if (u_dyn <= u_dyn_th) then
-          B=0.1*kappa*(Re**0.5)  !Zilitinkevich
-        else
-            B=(kappa*10.0**(-0.4*z0_m/0.07))*(Re**0.45) !Chen and Zhang
-        end if
-       
-       
-        ! --- define roughness [thermal]
-        z0_t = z0_m / exp(B)
-
-    end subroutine
-    
-! --------------------------------------------------------------------------------
-        subroutine get_thermal_roughness_re(z0_t, B, &
-            z0_m, Re)
-    ! ----------------------------------------------------------------------------
-        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]
-       
-       
-      
-       
-        
-        B=alog(-0.56*(4.0*(Re)**(0.5)-3.4))     !Repina, 2023
-
-
-       ! --- define roughness [thermal]
-        z0_t = z0_m / exp(B)   
-       
-    end subroutine
- 
-
-
-end module z0t_all_surface