sfx_surface.f90

#include "../includeF/sfx_def.fi"

module sfx_surface
    !< @brief surface roughness parameterizations

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

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

    ! public interface
    ! --------------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
    public :: set_c_struct_sfx_surface_param_values
#endif
    public :: get_charnock_roughness
    public :: get_thermal_roughness
    ! --------------------------------------------------------------------------------


    ! *: add surface type as input value
    integer, public, parameter :: surface_ocean = 0     !< ocean surface
    integer, public, parameter :: surface_land = 1      !< land surface
    integer, public, parameter :: surface_lake = 2      !< lake surface
    integer, public, parameter :: surface_snow = 3      !< snow covered surface

    character(len = 16), parameter :: surface_ocean_tag = 'ocean'
    character(len = 16), parameter :: surface_land_tag = 'land'
    character(len = 16), parameter :: surface_lake_tag = 'lake'
    character(len = 16), parameter :: surface_snow_tag = 'snow'

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

    !< 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
    real, parameter :: B_max_ocean = 8.0
    real, parameter :: B_max_lake = 8.0

contains

    ! surface type definition
    ! --------------------------------------------------------------------------------
    function get_surface_id(tag) result(id)
        implicit none
        character(len=*), intent(in) :: tag
        integer :: id

        id = - 1
        if (trim(tag) == trim(surface_ocean_tag)) then
            id = surface_ocean
        else if (trim(tag) == trim(surface_land_tag)) then
            id = surface_land
        else if (trim(tag) == trim(surface_lake_tag)) then
            id = surface_lake
        else if (trim(tag) == trim(surface_snow_tag)) then
            id = surface_snow
        end if

    end function

    function get_surface_tag(id) result(tag)
        implicit none
        integer :: id
        character(len=:), allocatable :: tag

        tag = 'undefined'
        if (id == surface_ocean) then
            tag = surface_ocean_tag
        else if (id == surface_land) then
            tag = surface_land_tag
        else if (id == surface_lake) then
            tag = surface_lake_tag
        else if (id == surface_snow) then
            tag = surface_snow_tag
        end if 

    end function

#if defined(INCLUDE_CXX)
    subroutine set_c_struct_sfx_surface_param_values(surface_param)
        use sfx_data
        implicit none
        type (sfx_surface_param), intent(inout) :: surface_param
        surface_param%surface_ocean = surface_ocean
        surface_param%surface_land = surface_land
        surface_param%surface_lake = surface_lake

        surface_param%kappa = kappa
        surface_param%gamma_c = gamma_c
        surface_param%Re_visc_min = Re_visc_min
        surface_param%h_charnock = h_charnock
        surface_param%c1_charnock = c1_charnock
        surface_param%c2_charnock = c2_charnock
        surface_param%Re_rough_min = Re_rough_min
        surface_param%B1_rough = B1_rough
        surface_param%B2_rough = B2_rough
        surface_param%B3_rough = B3_rough
        surface_param%B4_rough = B4_rough
        surface_param%B_max_lake = B_max_lake
        surface_param%B_max_ocean = B_max_ocean
        surface_param%B_max_land = B_max_land
    end subroutine set_c_struct_sfx_surface_param_values
#endif
    ! 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
    ! --------------------------------------------------------------------------------

    ! thermal roughness definition
    ! --------------------------------------------------------------------------------
    subroutine get_thermal_roughness(z0_t, B, &
            z0_m, Re, surface_type)
        ! ----------------------------------------------------------------------------
        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]
        integer, intent(in) :: surface_type     !< = [ocean] || [land] || [lake]
        ! ----------------------------------------------------------------------------

        ! --- local variables
        ! ----------------------------------------------------------------------------

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

        !   --- apply max restriction based on surface type
        if (surface_type == surface_ocean) then
            B = min(B, B_max_ocean)
        else if (surface_type == surface_lake) then
            B = min(B, B_max_lake)
        else if (surface_type == surface_land) then
            B = min(B, B_max_land)
        end if

        ! --- define roughness [thermal]
        z0_t = z0_m / exp(B)

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

end module sfx_surface