diff --git a/srcF/module_z0t_ocean.f90 b/srcF/module_z0t_ocean.f90
new file mode 100644
index 0000000000000000000000000000000000000000..bd0d6e5243ccd012250dbe6a003b0ac63cc09ceb
--- /dev/null
+++ b/srcF/module_z0t_ocean.f90
@@ -0,0 +1,125 @@
+module module_z0t_ocean
+    !< @brief surface thermal roughness parameterizations for ocean
+
+    implicit none
+    
+   
+    public :: get_thermal_roughness_kl
+    public :: get_thermal_roughness_ca
+    public :: get_thermal_roughness_zm
+    public :: get_thermal_roughness_br
+    
+
+
+
+    ! --------------------------------------------------------------------------------
+    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
+    
+
+    contains
+
+    
+     ! thermal roughness definition by Kazakov, Lykosov
+     ! --------------------------------------------------------------------------------
+       subroutine get_thermal_roughness_kl(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 Cahill, A.T., Parlange, M.B., Albertson, J.D., 1997.
+    ! --------------------------------------------------------------------------------
+        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 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 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
+   
+end module module_z0t_ocean
diff --git a/srcF/module_z0t_snow.f90 b/srcF/module_z0t_snow.f90
new file mode 100644
index 0000000000000000000000000000000000000000..a8545a9409a2217920b5f7bd0d7e9163e0566256
--- /dev/null
+++ b/srcF/module_z0t_snow.f90
@@ -0,0 +1,125 @@
+module module_z0t_snow
+    !< @brief surface thermal roughness parameterizations for snow
+
+    implicit none
+    
+   
+    public :: get_thermal_roughness_kl
+    public :: get_thermal_roughness_ca
+    public :: get_thermal_roughness_zm
+    public :: get_thermal_roughness_br
+    
+
+
+
+    ! --------------------------------------------------------------------------------
+    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_snow = 8.0
+    
+
+    contains
+
+    
+     ! thermal roughness definition by Kazakov, Lykosov
+     ! --------------------------------------------------------------------------------
+       subroutine get_thermal_roughness_kl(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_snow)
+       
+       z0_t = z0_m / exp(B)
+
+     end subroutine
+    
+        
+
+    ! --------------------------------------------------------------------------------
+     ! thermal roughness definition by Cahill, A.T., Parlange, M.B., Albertson, J.D., 1997.
+    ! --------------------------------------------------------------------------------
+        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 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 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
+   
+end module module_z0t_snow