module module_z0t
    !< @brief surface thermal roughness parameterizations for land

    implicit none
    
   
    public :: get_thermal_roughness_kl
    public :: get_thermal_roughness_cz
    public :: get_thermal_roughness_zi
    public :: get_thermal_roughness_ca
    public :: get_thermal_roughness_zm
    public :: get_thermal_roughness_ot
    public :: get_thermal_roughness_du
    public :: get_thermal_roughness_mix



    ! --------------------------------------------------------------------------------
    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_land = 2.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_land)
       
       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 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


end module module_z0t
