rename to sfx_z0m/z0t 2

srcF/sfx_z0m_all_surface.f90

+module sfx_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  sfx_z0m_all_surface
\ No newline at end of file

srcF/sfx_z0t_all_surface.f90

+module sfx_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 sfx_z0t_all_surface