module sfx_z0m_all_surface
   !< @brief surface roughness parameterizations 

    use sfx_phys_const
    implicit none
    public :: get_dynamic_roughness_ch
    public :: get_dynamic_roughness_map
    public :: get_dynamic_roughness_ow
    public :: get_dynamic_roughness_fetch
    public :: get_dynamic_roughness_and
    public :: get_dynamic_roughness_coast1
    public :: get_dynamic_roughness_coast2
    public :: get_dynamic_roughness_coast3
    ! --------------------------------------------------------------------------------


    

    ! --------------------------------------------------------------------------------
    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
        !!write(*,*) 'ch', z0_m, u_dyn0
    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
        !write(*,*) 'ow', z0_m, u_dyn0
        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)
        !write(*,*) 'map', z0_m, u_dyn0
        end subroutine
     ! --------------------------------------------------------------------------------

        subroutine get_dynamic_roughness_and(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, u_star_prev, nu, g
            
            
            real :: tolerance
            integer :: max_iterations, iter
        
            
            nu = 1.7e-5  
            g = 9.81     
            u_dyn0 = 0.2 
            
            
            tolerance = 1.0e-5  
            max_iterations = 10
        
            
            do iter = 1, max_iterations
                
                u_star_prev = u_dyn0
        
                
                z0_m = (0.135 * nu / u_dyn0) + (0.035 * u_dyn0**2 / g) * &
                      (1.0 + exp(-((u_dyn0 - 0.18) / 0.1)**2))
        
                      h0_m = h / z0_m
                      u_dyn0 =  U * kappa / log(h0_m) 
        
                
                if (abs(u_dyn0 - u_star_prev) < tolerance) exit
            end do

            u_dyn0 = U * kappa / log(h0_m)

           
             
            end subroutine
         ! --------------------------------------------------------------------------------

            subroutine get_dynamic_roughness_coast1(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
                ! ----------------------------------------------------------------------------
        ! Taylor&Yelland formulation
                ! --- local variables
                real :: Uc, c, c_min, a, b                  ! wind speed at h_charnock [m/s] & Cd^(-1)
        
               real :: pi=3.14159265
               real :: hs, Tp, Lp
               integer :: i, j
                ! ----------------------------------------------------------------------------
        
                Uc = U
                a = 0
                b = 0
                c_min = log(h_charnock) / kappa
        
                do i = 1, maxiters
                  hs = 0.0248*(Uc**2.) !hs is the significant wave height
                  Tp = 0.729*max(Uc,0.1) !Tp dominant wave period
                  Lp = g*Tp**2/(2*pi) !Lp is the wavelength of the dominant wave
                    do j = 1, maxiters
                        z0_m = 1200.*hs*(hs/Lp)**4.5
                        c = log(h_charnock / z0_m) / kappa
                        if (Uc <= 8.0e0) then
                        a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
                        c = max(c - a, c_min)
                        b = c
                        z0_m = h_charnock * exp(-kappa*c)
                        end if
                    end do
                z0_m = max( z0_m, 1.27e-7)  !These max/mins were suggested by
                z0_m = min( z0_m, 2.85e-3)  !Davis et al. (2008)
                Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
                end do
                u_dyn0 = Uc * kappa / log (h_charnock / z0_m)
        !       c = Uc/u_dyn0
        !       write (*,*) 'out1', u_dyn0
        
            end subroutine
        ! Hwang
            subroutine get_dynamic_roughness_coast2(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, c, c1, a, b, c_min                  ! wind speed at h_charnock [m/s] & Cd^(-1)
               integer :: i, j
                ! ----------------------------------------------------------------------------
                Uc = U
                a = 0
                b = 0
                c_min = log(h_charnock) / kappa
        
                do i = 1, maxiters
                    c1 = 0.016 * Uc**2
                    do j = 1, maxiters
                        c = 1.0 / (0.01*sqrt(8.058 + 0.967 * Uc - c1))
                        if (Uc <= 8.0e0) then
                        a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
                        c = max(c - a, c_min)
                        b = c
                        end if
                    end do
                    z0_m = h_charnock * exp(-kappa*c)
                    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
        !       write (*,*) 'out1', u_dyn0
            ! --------------------------------------------------------------------------------
            end subroutine
        ! Zhao et al 2015 10^3*z0=15.6*u_dyn0**2/g + 10**(-2)
            subroutine get_dynamic_roughness_coast3(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, u_dyn, z0_m0                  ! wind speed at h_charnock [m/s]
        
                real :: c, c_min, a, b
                real :: f1
        
                integer :: i, j
                ! ----------------------------------------------------------------------------
                Uc = U
                u_dyn = U / 28.0
                a = 0
                b = 0
                z0_m0=0.082
                c_min = log(h_charnock) / kappa
                    if (Uc >8.00 .AND. Uc <= 35.0e0) then
                    do i = 1, maxiters
                        f1 = 15.6*u_dyn**2/g
                        z0_m0 = 0.001 * (f1 + 0.01)
                        c = log(h_charnock / z0_m0) / kappa
                    end do
                    end if
                    if (Uc <= 8.0e0) then
                    do j = 1, maxiters
                        a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
                        c = max(c - a, c_min)
                        b = c
                        z0_m0 = h_charnock * exp(-c * kappa)
                    end do
                    else
                        z0_m0 = 2.82
                    end if
                    z0_m = max(z0_m0, 0.000015e0)
                    Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
        
                ! --- define dynamic velocity in neutral conditions
        !        u_dyn0 = Uc / c
                u_dyn0 = Uc * kappa / log (h_charnock / z0_m)
        
            end subroutine
        
   
    end module  sfx_z0m_all_surface