From fa1e60f54e483826b9a932e12cd047b2189544a2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?=D0=92=D0=B8=D0=BA=D1=82=D0=BE=D1=80=D0=B8=D1=8F=20=D0=A1?=
 =?UTF-8?q?=D1=83=D1=8F=D0=B7=D0=BE=D0=B2=D0=B0?=
 <viktoriasuazova@MacBook-Pro-Viktoria.local>
Date: Wed, 19 Mar 2025 16:54:08 +0300
Subject: [PATCH] mosaic surt tupe added and z0t_map added for TERM

---
 srcF/sfx_api_inmcm.f90          |  5 +++--
 srcF/sfx_api_term.f90           | 12 +++++++++++-
 srcF/sfx_data.f90               |  8 +++++++-
 srcF/sfx_esm.f90                | 19 +++++++++++--------
 srcF/sfx_log.f90                | 20 ++++++++++++--------
 srcF/sfx_most.f90               | 27 ++++++++++++---------------
 srcF/sfx_most_snow.f90          | 20 ++++++++++++--------
 srcF/sfx_run.f90                |  2 ++
 srcF/sfx_sheba.f90              | 22 +++++++++++++---------
 srcF/sfx_sheba_coare.f90        | 21 ++++++++++++---------
 srcF/sfx_sheba_noniterative.f90 | 20 ++++++++++++--------
 srcF/sfx_surface.f90            | 33 +++++++++++++++++++++++++--------
 srcF/sfx_z0t_all_surface.f90    | 21 +++++++++++++++++++--
 13 files changed, 151 insertions(+), 79 deletions(-)

diff --git a/srcF/sfx_api_inmcm.f90 b/srcF/sfx_api_inmcm.f90
index 37f45df..ef58eee 100644
--- a/srcF/sfx_api_inmcm.f90
+++ b/srcF/sfx_api_inmcm.f90
@@ -20,7 +20,7 @@ module sfx_api_inmcm
 contains
     
     ! --------------------------------------------------------------------------------
-    subroutine inmcm_to_sfx_in_cell(meteo, arg, IVEG_sfx, depth_inm, lai_inm)
+    subroutine inmcm_to_sfx_in_cell(meteo, arg, IVEG_sfx, depth_inm, lai_inm, z0t_inm)
         !> @brief converts legacy arg [AR1 INMCM format] array to sfx meteo input
         ! ----------------------------------------------------------------------------
         implicit none
@@ -29,7 +29,7 @@ contains
         integer,intent(in)   :: IVEG_sfx
         real,intent(in) :: depth_inm
         real,intent(in) :: lai_inm
-
+        real, intent(in) :: z0t_inm
         ! ----------------------------------------------------------------------------
         
         
@@ -42,6 +42,7 @@ contains
         meteo%depth = depth_inm
         meteo%lai = lai_inm
         meteo%surface_type = IVEG_sfx
+        meteo%z0_t = z0t_inm
         !write(*,*) 'surface_type, IVEG_sfx', meteo%surface_type, IVEG_sfx
     end subroutine inmcm_to_sfx_in_cell
     ! --------------------------------------------------------------------------------
diff --git a/srcF/sfx_api_term.f90 b/srcF/sfx_api_term.f90
index a7a08c5..e0abb1d 100644
--- a/srcF/sfx_api_term.f90
+++ b/srcF/sfx_api_term.f90
@@ -20,12 +20,17 @@ module sfx_api_term
 contains
     
     ! --------------------------------------------------------------------------------
-    subroutine term_to_sfx_in_cell(meteo, arg)
+    subroutine term_to_sfx_in_cell(meteo, arg, IVEG_sfx, depth_term, lai_term, z0t_term)
         !> @brief converts legacy arg [AR1 INMCM format] array to sfx meteo input
         ! ----------------------------------------------------------------------------
         implicit none
         type (meteoDataType), intent(inout) :: meteo
         real, dimension(6), intent(in)      :: arg
+        integer,intent(in)   :: IVEG_sfx
+        real,intent(in) :: depth_term
+        real,intent(in) :: lai_term
+        real, intent(in) :: z0t_term
+
         ! ----------------------------------------------------------------------------
         
         
@@ -35,6 +40,11 @@ contains
         meteo%dQ = arg(4)
         meteo%h = arg(5)
         meteo%z0_m = arg(6)
+        meteo%surface_type = IVEG_sfx
+        meteo%depth = depth_term
+        meteo%lai = lai_term
+        meteo%z0_t = z0t_term
+
 
     end subroutine term_to_sfx_in_cell
     ! --------------------------------------------------------------------------------
diff --git a/srcF/sfx_data.f90 b/srcF/sfx_data.f90
index b900a14..c2a46b0 100644
--- a/srcF/sfx_data.f90
+++ b/srcF/sfx_data.f90
@@ -35,6 +35,7 @@ module sfx_data
         real(C_FLOAT) :: Tsemi   !< semi-sum of potential temperature at 'h' and at surface [K]
         real(C_FLOAT) :: dQ      !< difference between humidity at 'h' and at surface [g/g]
         real(C_FLOAT) :: z0_m    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+        real(C_FLOAT) :: z0_t
         real(C_FLOAT) :: depth
         real(C_FLOAT) :: lai
         integer(C_INT) :: surface_type
@@ -49,6 +50,7 @@ module sfx_data
         real, allocatable :: Tsemi(:)   !< semi-sum of potential temperature at 'h' and at surface [K]
         real, allocatable :: dQ(:)      !< difference between humidity at 'h' and at surface [g/g]
         real, allocatable :: z0_m(:)    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+        real, allocatable :: z0_t(:)    !< surface thermal roughness (should be < 0 for water bodies surface)
         real, allocatable :: depth(:)  
         real, allocatable :: lai(:)  
         integer, allocatable :: surface_type(:)  
@@ -63,6 +65,7 @@ module sfx_data
         type(C_PTR) :: Tsemi   !< semi-sum of potential temperature at 'h' and at surface [K]
         type(C_PTR) :: dQ      !< difference between humidity at 'h' and at surface [g/g]
         type(C_PTR) :: z0_m    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
+        type(C_PTR) :: z0_t
         type(C_PTR) :: depth
         type(C_PTR) :: lai
         type(C_PTR) :: surface_type
@@ -167,6 +170,7 @@ contains
         allocate(meteo%Tsemi(n))
         allocate(meteo%dQ(n))
         allocate(meteo%z0_m(n))
+        allocate(meteo%z0_t(n))
         allocate(meteo%depth(n))
         allocate(meteo%lai(n))
         allocate(meteo%surface_type(n))
@@ -186,6 +190,7 @@ contains
         meteo_C%Tsemi = c_loc(meteo%Tsemi)
         meteo_C%dQ = c_loc(meteo%dQ)
         meteo_C%z0_m = c_loc(meteo%z0_m)
+        meteo_C%z0_t = c_loc(meteo%z0_t)
         meteo_C%depth = c_loc(meteo%depth)
         meteo_C%lai = c_loc(meteo%lai)
         meteo_C%surface_type = c_loc(meteo%surface_type)
@@ -206,6 +211,7 @@ contains
         deallocate(meteo%Tsemi)
         deallocate(meteo%dQ)
         deallocate(meteo%z0_m)
+        deallocate(meteo%z0_t)
         deallocate(meteo%depth)
         deallocate(meteo%lai)
         deallocate(meteo%surface_type)
@@ -305,4 +311,4 @@ contains
     end subroutine push_sfx_data
     ! --------------------------------------------------------------------------------
 
-end module sfx_data
\ No newline at end of file
+end module sfx_data
diff --git a/srcF/sfx_esm.f90 b/srcF/sfx_esm.f90
index 52c59f1..6de4e93 100644
--- a/srcF/sfx_esm.f90
+++ b/srcF/sfx_esm.f90
@@ -141,7 +141,8 @@ 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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+                    z0_m = meteo%z0_m(i), z0_t = meteo%z0_t(i), &
+                    depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -176,6 +177,7 @@ 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_t    !< surface thermal roughness (should be < 0 for water bodies surface)
         real :: depth   
         real :: lai
         integer :: surface_type
@@ -183,7 +185,7 @@ contains
 
         ! --- local variables
         ! ----------------------------------------------------------------------------
-        real z0_t               !< thermal roughness [m]
+        !real z0_t               !< thermal roughness [m]
         real B                  !< = ln(z0_m / z0_t) [n/d]
         real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
 
@@ -209,7 +211,7 @@ contains
 
 
         real fval               !< just a shortcut for partial calculations
-        real z0_m1
+        real z0_m1, z0_t1
 #ifdef SFX_CHECK_NAN
         real NaN
 #endif
@@ -236,27 +238,28 @@ contains
         dQ = meteo%dQ
         h = meteo%h
         z0_m1 = meteo%z0_m
+        z0_t1 = meteo%z0_t
         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, ice_z0m_id, z0m_id)
+        forest_z0m_id, usersf_z0m_id, ice_z0m_id, mosaic_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)
+        forest_z0t_id, usersf_z0t_id, ice_z0t_id, mosaic_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)
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0_t1, z0t_id)
        
        ! --- define relative height
-            h0_m = h / z0_m
+        h0_m = h / z0_m
         ! --- define relative height [thermal]
         h0_t = h / z0_t
 
@@ -551,4 +554,4 @@ contains
     end subroutine
     ! --------------------------------------------------------------------------------
 
-end module sfx_esm
\ No newline at end of file
+end module sfx_esm
diff --git a/srcF/sfx_log.f90 b/srcF/sfx_log.f90
index 6a51cf3..02399f6 100644
--- a/srcF/sfx_log.f90
+++ b/srcF/sfx_log.f90
@@ -57,7 +57,9 @@ 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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+                    z0_m = meteo%z0_m(i), z0_t = meteo%z0_t(i), &
+                    depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -92,11 +94,12 @@ contains
         real :: z0_m    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
         real :: depth
         real :: lai
+        real :: z0_t
         ! ----------------------------------------------------------------------------
 
         ! --- local variables
         ! ----------------------------------------------------------------------------
-        real z0_t               !< thermal roughness [m]
+        !real z0_t               !< thermal roughness [m]
         real B                  !< = ln(z0_m / z0_t) [n/d]
         real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
 
@@ -112,7 +115,7 @@ contains
         real Pr_t_inv           !< invese Prandt number [n/d]
 
         real Cm, Ct             !< transfer coeff. for (momentum) & (heat) [n/d]
-
+        real z0_m1, z0_t1
         integer surface_type    !< surface type = (ocean || land)
 
 
@@ -142,22 +145,23 @@ contains
         dQ = meteo%dQ
         h = meteo%h
         z0_m = meteo%z0_m
+        z0_t1 = meteo%z0_t
         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, ice_z0m_id, z0m_id)
+        forest_z0m_id, usersf_z0m_id, ice_z0m_id, mosaic_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_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)
+        forest_z0t_id, usersf_z0t_id, ice_z0t_id, mosaic_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)
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0_t1, z0t_id)
        
        ! --- define relative height
             h0_m = h / z0_m
@@ -211,4 +215,4 @@ contains
     end subroutine
     ! --------------------------------------------------------------------------------
 
-end module sfx_log
\ No newline at end of file
+end module sfx_log
diff --git a/srcF/sfx_most.f90 b/srcF/sfx_most.f90
index 968693e..bb4fcbe 100644
--- a/srcF/sfx_most.f90
+++ b/srcF/sfx_most.f90
@@ -58,7 +58,9 @@ 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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+                    z0_m = meteo%z0_m(i), z0_t = meteo%z0_t(i), &
+                    depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -91,13 +93,14 @@ 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_t
         real :: depth
         real lai
         ! ----------------------------------------------------------------------------
 
         ! --- local variables
         ! ----------------------------------------------------------------------------
-        real z0_t               !< thermal roughness [m]
+        !real z0_t               !< thermal roughness [m]
         real B                  !< = ln(z0_m / z0_t) [n/d]
         real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
 
@@ -117,7 +120,7 @@ contains
         real Cm, Ct             !< transfer coeff. for (momentum) & (heat) [n/d]
 
         integer surface_type    !< surface type = (ocean || land)
-        real z0_m1
+        real z0_m1, z0_t1
 
 #ifdef SFX_CHECK_NAN
         real NaN
@@ -145,30 +148,24 @@ contains
         dQ = meteo%dQ
         h = meteo%h
         z0_m1 = meteo%z0_m
+        z0_t1 = meteo%z0_t
         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, ice_z0m_id, z0m_id)
+        forest_z0m_id, usersf_z0m_id, ice_z0m_id, mosaic_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)
+        forest_z0t_id, usersf_z0t_id, ice_z0t_id, mosaic_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)
-            
-            h0_m = h / z0_m
-    
-            h0_m = h / z0_m
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0_t1, z0t_id)
             
-
-        
-        Re = u_dyn0 * z0_m / nu_air
-       
+        h0_m = h / z0_m
         h0_t = h / z0_t
 
         ! --- define Ri-bulk
@@ -347,4 +344,4 @@ contains
     end subroutine
     ! --------------------------------------------------------------------------------
 
-end module sfx_most
\ No newline at end of file
+end module sfx_most
diff --git a/srcF/sfx_most_snow.f90 b/srcF/sfx_most_snow.f90
index cc7d250..5890fdc 100644
--- a/srcF/sfx_most_snow.f90
+++ b/srcF/sfx_most_snow.f90
@@ -61,7 +61,9 @@ 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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+                    z0_m = meteo%z0_m(i), z0_t = meteo%z0_t(i), &
+                    depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -94,13 +96,14 @@ 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_t
         real :: depth  
         real :: lai
         ! ----------------------------------------------------------------------------
 
         ! --- local variables
         ! ----------------------------------------------------------------------------
-        real z0_t               !< thermal roughness [m]
+        !real z0_t               !< thermal roughness [m]
         real B                  !< = ln(z0_m / z0_t) [n/d]
         real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
 
@@ -118,7 +121,7 @@ contains
         real Pr_t_inv           !< invese Prandt number [n/d]
 
         real Cm, Ct             !< transfer coeff. for (momentum) & (heat) [n/d]
-        real z0_m1
+        real z0_m1, z0_t1
         ! --- local additional variables
         ! ----------------------------------------------------------------------------
 		!real phi_m, phi_m
@@ -160,25 +163,26 @@ contains
         dQ = meteo%dQ
         h = meteo%h
         z0_m1 = meteo%z0_m
+        z0_t1 = meteo%z0_t
         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, ice_z0m_id, z0m_id)
+        forest_z0m_id, usersf_z0m_id, ice_z0m_id, mosaic_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)
+        forest_z0t_id, usersf_z0t_id, ice_z0t_id, mosaic_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)
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0_t1, z0t_id)
         
        ! --- define relative height
-            h0_m = h / z0_m
+        h0_m = h / z0_m
         ! --- define relative height [thermal]
         h0_t = h / z0_t
 
@@ -453,4 +457,4 @@ contains
 
         end subroutine
 
-end module sfx_most_snow
\ No newline at end of file
+end module sfx_most_snow
diff --git a/srcF/sfx_run.f90 b/srcF/sfx_run.f90
index a7cdf1a..5f93c14 100644
--- a/srcF/sfx_run.f90
+++ b/srcF/sfx_run.f90
@@ -132,6 +132,7 @@ contains
         ! --- setting height & roughness
         meteo_cell%h = dataset%h
         meteo_cell%z0_m = dataset%z0_m
+        meteo_cell%z0_t = dataset%z0_h
         meteo_cell%depth = dataset%depth
         meteo_cell%lai = dataset%lai
         meteo_cell%surface_type = dataset%surface_type
@@ -156,6 +157,7 @@ contains
             meteo%Tsemi(i) = meteo_cell%Tsemi
             meteo%dQ(i) = meteo_cell%dQ
             meteo%z0_m(i) = meteo_cell%z0_m
+            meteo%z0_t(i) = meteo_cell%z0_t
 
         enddo
         close(io)
diff --git a/srcF/sfx_sheba.f90 b/srcF/sfx_sheba.f90
index ad61fd9..89a16f7 100644
--- a/srcF/sfx_sheba.f90
+++ b/srcF/sfx_sheba.f90
@@ -140,7 +140,9 @@ 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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+                    z0_m = meteo%z0_m(i), z0_t = meteo%z0_t(i), &
+                    depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -173,13 +175,14 @@ 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_t 
         real :: lai
         real :: depth
         ! ----------------------------------------------------------------------------
 
         ! --- local variables
         ! ----------------------------------------------------------------------------
-        real z0_t               !< thermal roughness [m]
+        !real z0_t               !< thermal roughness [m]
         real B                  !< = ln(z0_m / z0_t) [n/d]
         real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
 
@@ -197,7 +200,7 @@ contains
         real Pr_t_inv           !< invese Prandt number [n/d]
 
         real Cm, Ct             !< transfer coeff. for (momentum) & (heat) [n/d]
-
+        real z0_m1, z0_t1
         integer surface_type    !< surface type = (ocean || land)
 
 
@@ -226,22 +229,23 @@ contains
         dT = meteo%dT
         dQ = meteo%dQ
         h = meteo%h
-        z0_m = meteo%z0_m
+        z0_m1 = meteo%z0_m
+        z0_t1 = meteo%z0_t
         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, ice_z0m_id, z0m_id)
+        forest_z0m_id, usersf_z0m_id, ice_z0m_id, mosaic_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_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)
+        forest_z0t_id, usersf_z0t_id, ice_z0t_id, mosaic_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)
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0_t1, z0t_id)
         
         h0_t = h / z0_t
 
@@ -456,4 +460,4 @@ contains
     end subroutine
     ! --------------------------------------------------------------------------------
 
-end module sfx_sheba
\ No newline at end of file
+end module sfx_sheba
diff --git a/srcF/sfx_sheba_coare.f90 b/srcF/sfx_sheba_coare.f90
index 6c7eb23..f120190 100644
--- a/srcF/sfx_sheba_coare.f90
+++ b/srcF/sfx_sheba_coare.f90
@@ -156,7 +156,9 @@ 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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+                    z0_m = meteo%z0_m(i), z0_t = meteo%z0_t(i), &
+                    depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+
 
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
@@ -189,7 +191,7 @@ 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 :: hpbl
+        real :: z0_t
         real :: depth   
         real :: lai
         integer :: surface_type
@@ -197,7 +199,7 @@ contains
 
         ! --- local variables
         ! ----------------------------------------------------------------------------
-        real z0_t               !< thermal roughness [m]
+        !real z0_t               !< thermal roughness [m]
         real B                  !< = ln(z0_m / z0_t) [n/d]
         real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
 
@@ -215,7 +217,7 @@ contains
 
         real psi_m, psi_h       !< universal functions (momentum) & (heat) [n/d]
         real psi0_m, psi0_h       !< universal functions (momentum) & (heat) [n/d]
-        real z0_m1
+        real z0_m1, z0_t1
 !        real psi_m_BD, psi_h_BD       !< universal functions (momentum) & (heat) [n/d]
 !        real psi0_m_BD, psi0_h_BD       !< universal functions (momentum) & (heat) [n/d]
 !        real psi_m_conv, psi_h_conv       !< universal functions (momentum) & (heat) [n/d]
@@ -261,26 +263,27 @@ contains
         dQ = meteo%dQ
         h = meteo%h
         z0_m1 = meteo%z0_m
+        z0_t1 = meteo%z0_t
         depth = meteo%depth
         lai = meteo%lai
         surface_type=meteo%surface_type
         !hpbl = meteo%hpbl
 
         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)
+        forest_z0m_id, usersf_z0m_id, ice_z0m_id, mosaic_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)
+        forest_z0t_id, usersf_z0t_id, ice_z0t_id, mosaic_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)
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0_t1, z0t_id)
         ! --- define relative height
-            h0_m = h / z0_m
+        h0_m = h / z0_m
         ! --- define relative height [thermal]
         h0_t = h / z0_t
       
@@ -802,4 +805,4 @@ contains
     end subroutine
 
 
-end module sfx_sheba_coare
\ No newline at end of file
+end module sfx_sheba_coare
diff --git a/srcF/sfx_sheba_noniterative.f90 b/srcF/sfx_sheba_noniterative.f90
index aee16a7..1c2c039 100644
--- a/srcF/sfx_sheba_noniterative.f90
+++ b/srcF/sfx_sheba_noniterative.f90
@@ -153,7 +153,9 @@ 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), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+            z0_m = meteo%z0_m(i), z0_t = meteo%z0_t(i), &
+            depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))
+
             call get_surface_fluxes(sfx_cell, meteo_cell, numerics)
 
             call push_sfx_data(sfx, sfx_cell, i)
@@ -184,7 +186,8 @@ contains
         real :: dT      !< difference between potential temperature at 'h' and at surface [K]
         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_m
+        real :: z0_t    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
         real :: depth   
         real :: lai
         integer :: surface_type
@@ -192,7 +195,7 @@ contains
 
         ! --- local variables
         ! ----------------------------------------------------------------------------
-        real z0_t               !< thermal roughness [m]
+        !real z0_t               !< thermal roughness [m]
         real B                  !< = ln(z0_m / z0_t) [n/d]
         real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
 
@@ -228,7 +231,7 @@ contains
         real w_snow, deltaS, h_salt
         real sigma_w, sigma_r
         real S_mean, S_salt, Linv
-        real z0_m1
+        real z0_m1, z0_t1
         integer i
 #ifdef SFX_CHECK_NAN
         real NaN
@@ -256,24 +259,25 @@ contains
         dQ = meteo%dQ
         h = meteo%h
         z0_m1 = meteo%z0_m
+        z0_t1 = meteo%z0_t
         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, ice_z0m_id, z0m_id)
+        forest_z0m_id, usersf_z0m_id, ice_z0m_id, mosaic_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)
+        forest_z0t_id, usersf_z0t_id, ice_z0t_id, mosaic_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)
+        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0_t1, z0t_id)
         ! --- define relative height
             h0_m = h / z0_m
         ! --- define relative height [thermal]
@@ -698,4 +702,4 @@ contains
     end subroutine
     ! --------------------------------------------------------------------------------
 
-end module sfx_sheba_noniterative
\ No newline at end of file
+end module sfx_sheba_noniterative
diff --git a/srcF/sfx_surface.f90 b/srcF/sfx_surface.f90
index 6cd033c..8471130 100644
--- a/srcF/sfx_surface.f90
+++ b/srcF/sfx_surface.f90
@@ -37,6 +37,7 @@ module sfx_surface
     integer, public, parameter :: surface_forest = 4      !< forest covered surface
     integer, public, parameter :: surface_user = 5      !< coast covered surface
     integer, public, parameter :: surface_ice = 6      !< ice covered surface
+    integer, public, parameter :: surface_mosaic = 7      !< mosaic covered surface
 
     character(len = 16), parameter :: surface_ocean_tag = 'ocean'
     character(len = 16), parameter :: surface_land_tag = 'land'
@@ -45,6 +46,7 @@ module sfx_surface
     character(len = 16), parameter :: surface_forest_tag = 'forest'
     character(len = 16), parameter :: surface_user_tag = 'user'
     character(len = 16), parameter :: surface_ice_tag = 'ice'
+    character(len = 16), parameter :: surface_mosaic_tag = 'mosaic'
 
 
     integer, public, parameter :: z0m_ch = 0 
@@ -69,7 +71,7 @@ module sfx_surface
     integer, public, parameter :: z0t_zm = 9      
     integer, public, parameter :: z0t_mix = 10      
     integer, public, parameter :: z0t_user = 11   
-
+    integer, public, parameter :: z0t_map_id = 12
 
     character(len = 16), parameter :: z0m_tag_ch = 'charnock'
     character(len = 16), parameter :: z0m_tag_fe = 'fetch'
@@ -94,6 +96,7 @@ module sfx_surface
     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'
+    character(len = 16), parameter :: z0t_tag_map = 'zt_map'
     
 
     integer, public, parameter :: ocean_z0m_id = z0m_ch     !< ocean surface
@@ -103,7 +106,7 @@ module sfx_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 :: ice_z0m_id = z0m_map_id       !< ice surface  
-
+    integer, public, parameter :: mosaic_z0m_id = z0m_map_id
 
     integer, public, parameter :: ocean_z0t_id = z0t_kl_water     !< ocean surface
     integer, public, parameter :: land_z0t_id = z0t_kl_land       !< land surface
@@ -112,7 +115,7 @@ module sfx_surface
     integer, public, parameter :: forest_z0t_id = z0t_mix      !< forest csurface  
     integer, public, parameter :: usersf_z0t_id = z0t_kl_water      !< user surface  
     integer, public, parameter :: ice_z0t_id = z0t_kl_land     !< user surface  
-    
+    integer, public, parameter :: mosaic_z0t_id = z0t_map_id
     ! --------------------------------------------------------------------------------
     real, parameter, private :: kappa = 0.40         !< von Karman constant [n/d]
     ! --------------------------------------------------------------------------------
@@ -170,6 +173,8 @@ contains
             id = surface_user
         else if (trim(tag) == trim(surface_ice_tag)) then
             id = surface_ice
+        else if (trim(tag) == trim(surface_mosaic_tag)) then
+            id = surface_mosaic
         end if
 
     end function
@@ -194,6 +199,8 @@ contains
             tag = surface_user_tag
         else if (id == surface_ice) then
             tag = surface_ice_tag
+        else if (id == surface_mosaic) then
+            tag = surface_mosaic_tag
         end if 
 
     end function
@@ -280,6 +287,8 @@ contains
             z0t_id = z0t_mix
         else if (trim(tag_z0t) == trim(z0t_tag_user)) then
             z0t_id = z0t_user
+        else if (trim(tag_z0t) == trim(z0t_tag_map)) then
+            z0t_id = z0t_map_id
         end if
         
     end function
@@ -350,7 +359,7 @@ contains
     ! dynamic roughness definition
     ! --------------------------------------------------------------------------------
       subroutine 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)
+        forest_z0m_id, usersf_z0m_id, ice_z0m_id, mosaic_z0m_id, z0m_id)
     ! ----------------------------------------------------------------------------
         integer, intent(out) :: z0m_id              
 
@@ -363,6 +372,7 @@ contains
         integer, intent(in) :: forest_z0m_id
         integer, intent(in) :: usersf_z0m_id
         integer, intent(in) :: ice_z0m_id
+        integer, intent(in) :: mosaic_z0m_id
  ! ---------------------------------------------------------------------------
 
     if (surface_type == surface_ocean) then
@@ -379,6 +389,8 @@ contains
      z0m_id  = usersf_z0m_id
     else if (surface_type == surface_ice) then
     z0m_id  = ice_z0m_id
+    else if (surface_type == surface_mosaic) then
+    z0m_id  = mosaic_z0m_id
     end if 
     !write (*,*) z0m_id, surface_type
  end subroutine
@@ -426,7 +438,7 @@ 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, ice_z0t_id, z0t_id)
+    forest_z0t_id, usersf_z0t_id, ice_z0t_id, mosaic_z0t_id, z0t_id)
 ! ----------------------------------------------------------------------------
     integer, intent(out) :: z0t_id              
 
@@ -439,7 +451,7 @@ subroutine get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t
     integer, intent(in) :: forest_z0t_id
     integer, intent(in) :: usersf_z0t_id
     integer, intent(in) :: ice_z0t_id
-
+    integer, intent(in) :: mosaic_z0t_id
 ! ---------------------------------------------------------------------------
 
     if (surface_type == surface_ocean) then
@@ -456,13 +468,16 @@ subroutine get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t
     z0t_id  = usersf_z0t_id
     else if (surface_type == surface_ice) then
     z0t_id  = ice_z0t_id
+    else if (surface_type == surface_mosaic) then
+    z0t_id  = mosaic_z0t_id
+
     end if 
     
     end subroutine
 
 !---------------------------------------------------------------
     subroutine get_thermal_roughness_all(z0_t, B, &
-    z0_m, Re, u_dyn, LAI, z0t_id)
+    z0_m, Re, u_dyn, LAI, z0t_map, z0t_id)
 ! ----------------------------------------------------------------------------
     real, intent(out) :: z0_t               !< thermal roughness [m]
     real, intent(out) :: B                  !< = log(z0_m / z0_t) [n/d]
@@ -471,7 +486,7 @@ subroutine get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t
     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]
-    
+    real, intent(in) :: z0t_map
 
     integer, intent(in) :: z0t_id
 ! ---------------------------------------------------------------------------
@@ -502,6 +517,8 @@ subroutine get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t
     call get_thermal_roughness_mix(z0_t, B, z0_m, u_dyn, Re)
     else if (z0t_id == z0t_user) then
     write(*, *) 'z0t_user'
+    else if (z0t_id == z0t_map_id) then
+    call get_thermal_roughness_map(z0_t, B, z0_m, z0t_map)
     end if
     
 
diff --git a/srcF/sfx_z0t_all_surface.f90 b/srcF/sfx_z0t_all_surface.f90
index 503eeeb..2bba114 100644
--- a/srcF/sfx_z0t_all_surface.f90
+++ b/srcF/sfx_z0t_all_surface.f90
@@ -277,10 +277,10 @@ module sfx_z0t_all_surface
         z0_t = z0_m*(-0.56*(4.0*(Re1)**(0.5)-3.4))
         if (z0_t<0)  z0_t = z0_m*0.9
         if (z0_m>z0_t) then
-        B=alog(z0_m/z0_t) 
+        B=log(z0_m/z0_t) 
         else
         z0_t = z0_m*0.9 
-        B=alog(z0_m/z0_t)   
+        B=log(z0_m/z0_t)   
         end if
 
        ! --- define roughness [thermal]
@@ -288,6 +288,23 @@ module sfx_z0t_all_surface
         !write(*,*)  z0_m, z0_t, B, Re1, 're'
     end subroutine
  
+    ! --------------------------------------------------------------------------------
+    ! thermal roughness definition z0_t = z0t_map
+    ! --------------------------------------------------------------------------------
+    subroutine get_thermal_roughness_map(z0_t, B, z0_m, z0t_map)
+    
+        ! ----------------------------------------------------------------------------
+        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) :: z0t_map             !< thermal roughness from map [m] (already known value)
+        ! ----------------------------------------------------------------------------
+       
+        z0_t=z0t_map
+        B=log(z0_m / z0_t)
+        
+    end subroutine
 
 
 end module sfx_z0t_all_surface
-- 
GitLab