sfx_log.f90

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

module sfx_log
    !< @brief simple log-roughness surface flux module

    ! modules used
    ! --------------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
    use sfx_common
#endif
    use sfx_data
    use sfx_surface
    use sfx_log_param
    ! --------------------------------------------------------------------------------

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

    ! public interface
    ! --------------------------------------------------------------------------------
    public :: get_surface_fluxes
    public :: get_surface_fluxes_vec
    ! --------------------------------------------------------------------------------

    ! --------------------------------------------------------------------------------
    type, public :: numericsType
        integer :: maxiters_charnock = 10      !< maximum (actual) number of iterations in charnock roughness
    end type
    ! --------------------------------------------------------------------------------

contains

    ! --------------------------------------------------------------------------------
    subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
        !< @brief surface flux calculation for array data
        !< @details contains C/C++ & CUDA interface
        ! ----------------------------------------------------------------------------
        type (sfxDataVecType), intent(inout) :: sfx

        type (meteoDataVecType), intent(in) :: meteo
        type (numericsType), intent(in) :: numerics
        integer, intent(in) :: n
        ! ----------------------------------------------------------------------------

        ! --- local variables
        type (meteoDataType)  meteo_cell
        type (sfxDataType) sfx_cell
        integer i
        ! ----------------------------------------------------------------------------

        do i = 1, n
            meteo_cell = meteoDataType(&
                    h = meteo%h(i), &
                    U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
                    z0_m = meteo%z0_m(i))

            call get_surface_fluxes(sfx_cell, meteo_cell, numerics)

            call push_sfx_data(sfx, sfx_cell, i)
        end do

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

    ! --------------------------------------------------------------------------------
    subroutine get_surface_fluxes(sfx, meteo, numerics)
        !< @brief surface flux calculation for single cell
        !< @details contains C/C++ interface
        ! ----------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
        use ieee_arithmetic
#endif

        type (sfxDataType), intent(out) :: sfx

        type (meteoDataType), intent(in) :: meteo
        type (numericsType), intent(in) :: numerics
        ! ----------------------------------------------------------------------------

        ! --- meteo derived datatype name shadowing
        ! ----------------------------------------------------------------------------
        real :: h       !< constant flux layer height [m]
        real :: U       !< abs(wind speed) at 'h' [m/s]
        real :: dT      !< difference between potential temperature at 'h' and at surface [K]
        real :: Tsemi   !< semi-sum of potential temperature at 'h' and at surface [K]
        real :: dQ      !< difference between humidity at 'h' and at surface [g/g]
        real :: z0_m    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
        ! ----------------------------------------------------------------------------

        ! --- local variables
        ! ----------------------------------------------------------------------------
        real z0_t               !< thermal roughness [m]
        real B                  !< = ln(z0_m / z0_t) [n/d]
        real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]

        real u_dyn0             !< dynamic velocity in neutral conditions [m/s]
        real Re                 !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]

        real Rib                !< bulk Richardson number

        real psi_m, psi_h       !< universal functions (momentum) & (heat) [n/d]
        real phi_m, phi_h       !< stability functions (momentum) & (heat) [n/d]

        real Km                 !< eddy viscosity coeff. at h [m^2/s]
        real Pr_t_inv           !< invese Prandt number [n/d]

        real Cm, Ct             !< transfer coeff. for (momentum) & (heat) [n/d]

        integer surface_type    !< surface type = (ocean || land)


#ifdef SFX_CHECK_NAN
        real NaN
#endif
        ! ----------------------------------------------------------------------------

#ifdef SFX_CHECK_NAN
        ! --- checking if arguments are finite
        if (.not.(is_finite(meteo%U).and.is_finite(meteo%Tsemi).and.is_finite(meteo%dT).and.is_finite(meteo%dQ) &
                .and.is_finite(meteo%z0_m).and.is_finite(meteo%h))) then

            NaN = ieee_value(0.0, ieee_quiet_nan)   ! setting NaN
            sfx = sfxDataType(zeta = NaN, Rib = NaN, &
                    Re = NaN, B = NaN, z0_m = NaN, z0_t = NaN, &
                    Rib_conv_lim = NaN, &
                    Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN)
            return
        end if
#endif

        ! --- shadowing names for clarity
        U = meteo%U
        Tsemi = meteo%Tsemi
        dT = meteo%dT
        dQ = meteo%dQ
        h = meteo%h
        z0_m = meteo%z0_m

        ! --- define surface type
        if (z0_m < 0.0) then
            surface_type = surface_ocean
        else
            surface_type = surface_land
        end if

        if (surface_type == surface_ocean) then
            ! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
            call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
            ! --- define relative height
            h0_m = h / z0_m
        endif
        if (surface_type == surface_land) then
            ! --- define relative height
            h0_m = h / z0_m
            ! --- define dynamic velocity in neutral conditions
            u_dyn0 = U * kappa / log(h0_m)
        end if

        ! --- define thermal roughness & B = log(z0_m / z0_h)
        Re = u_dyn0 * z0_m / nu_air
        call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)

        ! --- define relative height [thermal]
        h0_t = h / z0_t

        ! --- define Ri-bulk
        Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2

        ! --- get the fluxes
        ! ----------------------------------------------------------------------------
        call get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
        ! ----------------------------------------------------------------------------

        phi_m = 1.0
        phi_h = 1.0 / Pr_t_0_inv

        ! --- define transfer coeff. (momentum) & (heat)
        Cm = kappa / psi_m
        Ct = kappa / psi_h

        ! --- define eddy viscosity & inverse Prandtl number
        Km = kappa * Cm * U * h / phi_m
        Pr_t_inv = phi_m / phi_h

        ! --- setting output
        sfx = sfxDataType(zeta = 0.0, Rib = Rib, &
                Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
                Rib_conv_lim = 0.0, &
                Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)

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

    ! universal functions
    ! --------------------------------------------------------------------------------
    subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
        !< @brief universal functions (momentum) & (heat): neutral case
        ! ----------------------------------------------------------------------------
        real, intent(out) :: psi_m, psi_h   !< universal functions

        real, intent(in) :: h0_m, h0_t      !< = z/z0_m, z/z0_h
        real, intent(in) :: B               !< = log(z0_m / z0_h)
        ! ----------------------------------------------------------------------------

        psi_m = log(h0_m)
        psi_h = log(h0_t) / Pr_t_0_inv
        !*: this looks redundant z0_t = z0_m in case |B| ~ 0
        if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv

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

end module sfx_log