module sfx_roughness
    !> @brief surface roughness parameterizations

    ! macro defs.
    ! --------------------------------------------------------------------------------
    ! --------------------------------------------------------------------------------

    ! modules used
    ! --------------------------------------------------------------------------------
    use sfx_phys_const
    ! --------------------------------------------------------------------------------

    ! directives list
    ! --------------------------------------------------------------------------------
    implicit none
    private
    ! --------------------------------------------------------------------------------

    ! --------------------------------------------------------------------------------
    real, parameter :: kappa = 0.40         !> von Karman constant [n/d]
    ! --------------------------------------------------------------------------------

    ! public interface
    ! --------------------------------------------------------------------------------
    public :: get_charnock_roughness
    ! --------------------------------------------------------------------------------

    !> 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
    ! --------------------------------------------------------------------------------

contains

    ! charnock roughness definition
    ! --------------------------------------------------------------------------------
    subroutine get_charnock_roughness(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
    ! --------------------------------------------------------------------------------

end module sfx_roughness
