diff --git a/docs/bibliography.bib b/docs/bibliography.bib
index 0caca54d..966fdc16 100644
--- a/docs/bibliography.bib
+++ b/docs/bibliography.bib
@@ -2,6 +2,19 @@
 #    Last author name (titlecase), followed by
 #    (no characters in-between) the year.
 
+
+@article{Beljaars1991,
+  title = {Flux Parameterization over Land Surfaces for Atmospheric Models},
+  author = {Beljaars, A. C. M., and A. A. M. Holtslag},
+  journal = {Journal of Applied Meteorology and Climatology},
+  volume = {30},
+  number = {3},
+  doi = {10.1175/1520-0450(1991)030<0327:FPOLSF>2.0.CO;2},
+  pages = {327–341},
+  year = {1991},
+  publisher = {American Meteorological Society}
+}
+
 @article{Nishizawa2018,
   title = {A Surface Flux Scheme Based on the Monin-Obukhov Similarity for Finite Volume Models},
   author = {Nishizawa, S and Kitamura, Y},
diff --git a/docs/src/index.md b/docs/src/index.md
index cf398506..ffda2733 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -37,6 +37,7 @@ SurfaceFluxes.UniversalFunctions
 SurfaceFluxes.UniversalFunctions.Gryanik
 SurfaceFluxes.UniversalFunctions.Grachev
 SurfaceFluxes.UniversalFunctions.Businger
+SurfaceFluxes.UniversalFunctions.Beljaars
 ```
 
 ```@docs
diff --git a/docs/src/plot_universal_functions.jl b/docs/src/plot_universal_functions.jl
index 5173898a..b3cc2e65 100644
--- a/docs/src/plot_universal_functions.jl
+++ b/docs/src/plot_universal_functions.jl
@@ -17,7 +17,7 @@ import Plots
 
 L = FT(10);
 
-ufts = (UF.GryanikType(), UF.BusingerType(), UF.GrachevType())
+ufts = (UF.GryanikType(), UF.BusingerType(), UF.GrachevType(), UF.BeljaarsType())
 
 universal_functions(uft) = UF.universal_func(uft, L, create_uf_parameters(toml_dict, uft))
 
diff --git a/parameters/create_parameters.jl b/parameters/create_parameters.jl
index eaa6c316..9a6b80e8 100644
--- a/parameters/create_parameters.jl
+++ b/parameters/create_parameters.jl
@@ -3,6 +3,42 @@ import SurfaceFluxes as SF
 import SurfaceFluxes.UniversalFunctions as UF
 import Thermodynamics as TD
 
+function create_uf_parameters(toml_dict, ::UF.BeljaarsType)
+    FT = CP.float_type(toml_dict)
+
+    aliases = [
+        "Pr_0_Beljaars",
+        "a_m_Beljaars",
+        "a_h_Beljaars",
+        "b_m_Beljaars",
+        "b_h_Beljaars",
+        "c_h_Beljaars",
+        "c_m_Beljaars",
+        "d_m_Beljaars",
+        "d_h_Beljaars",
+        "ζ_a_Beljaars",
+        "γ_Beljaars",
+    ]
+
+    pairs = CP.get_parameter_values!(toml_dict, aliases, "UniversalFunctions")
+    pairs = (; pairs...) # convert to NamedTuple
+
+    pairs = (;
+        Pr_0 = pairs.Pr_0_Beljaars,
+        a_m = pairs.a_m_Beljaars,
+        a_h = pairs.a_h_Beljaars,
+        b_m = pairs.b_m_Beljaars,
+        b_h = pairs.b_h_Beljaars,
+        c_m = pairs.b_m_Beljaars,
+        c_h = pairs.b_h_Beljaars,
+        d_m = pairs.b_m_Beljaars,
+        d_h = pairs.b_h_Beljaars,
+        ζ_a = pairs.ζ_a_Beljaars,
+        γ = pairs.γ_Beljaars,
+    )
+    return UF.BeljaarsParams{FT}(; pairs...)
+end
+
 function create_uf_parameters(toml_dict, ::UF.GryanikType)
     FT = CP.float_type(toml_dict)
 
diff --git a/src/Parameters.jl b/src/Parameters.jl
index bba79d6b..a2b954fa 100644
--- a/src/Parameters.jl
+++ b/src/Parameters.jl
@@ -38,6 +38,10 @@ a_m(ps::SurfaceFluxesParameters) = a_m(uf_params(ps))
 a_h(ps::SurfaceFluxesParameters) = a_h(uf_params(ps))
 b_m(ps::SurfaceFluxesParameters) = b_m(uf_params(ps))
 b_h(ps::SurfaceFluxesParameters) = b_h(uf_params(ps))
+c_m(ps::SurfaceFluxesParameters) = c_m(uf_params(ps))
+c_h(ps::SurfaceFluxesParameters) = c_h(uf_params(ps))
+d_m(ps::SurfaceFluxesParameters) = d_m(uf_params(ps))
+d_h(ps::SurfaceFluxesParameters) = d_h(uf_params(ps))
 ζ_a(ps::SurfaceFluxesParameters) = ζ_a(uf_params(ps))
 γ(ps::SurfaceFluxesParameters) = γ(uf_params(ps))
 
diff --git a/src/UniversalFunctions.jl b/src/UniversalFunctions.jl
index 3cd3479a..737fcd97 100644
--- a/src/UniversalFunctions.jl
+++ b/src/UniversalFunctions.jl
@@ -7,6 +7,7 @@ universal functions:
  - `Businger`
  - `Gryanik`
  - `Grachev`
+ - `Beljaars`
 """
 module UniversalFunctions
 
@@ -80,6 +81,9 @@ a_h(uf::AUF) = uf.params.a_h
 b_m(uf::AUF) = uf.params.b_m
 b_h(uf::AUF) = uf.params.b_h
 c_h(uf::AUF) = uf.params.c_h
+c_m(uf::AUF) = uf.params.c_m
+d_h(uf::AUF) = uf.params.d_h
+d_m(uf::AUF) = uf.params.d_m
 ζ_a(uf::AUF) = uf.params.ζ_a
 γ(uf::AUF) = uf.params.γ
 
@@ -496,12 +500,197 @@ function psi(uf::Grachev, ζ, tt::HeatTransport)
     end
 end
 
+#####
+##### Beljaars
+#####
+
+Base.@kwdef struct BeljaarsParams{FT} <: AbstractUniversalFunctionParameters{FT}
+    Pr_0::FT
+    a_m::FT
+    a_h::FT
+    b_m::FT
+    b_h::FT
+    c_m::FT
+    c_h::FT
+    d_m::FT
+    d_h::FT
+    ζ_a::FT
+    γ::FT
+end
+
+"""
+    Beljaars <: AbstractUniversalFunction{FT}
+
+# References
+ - [Beljaars1991](@cite)
+
+# Equations in reference:
+
+    `ϕ_m`: Derived from Eq. 28
+    `ϕ_h`: Derived from Eq. 32
+    `ψ_m`: Eq. 28
+    `ψ_h`: Eq. 32
+
+# Beljaars and Holtslag (1991) functions are used in stable conditions
+# In unstable conditions the functions of Businger (1971) are 
+# assigned by default. 
+
+# Fields
+
+$(DSE.FIELDS)
+"""
+struct Beljaars{FT, PS <: BeljaarsParams} <: AbstractUniversalFunction{FT}
+    "Monin-Obhukov Length"
+    L::FT
+    params::PS
+end
+
+struct BeljaarsType <: AbstractUniversalFunctionType end
+Beljaars() = BeljaarsType()
+
+# Nishizawa2018 Eq. A7
+f_momentum(uf::Beljaars, ζ) = sqrt(sqrt(1 - 15 * ζ))
+
+# Nishizawa2018 Eq. A8
+f_heat(uf::Beljaars, ζ) = sqrt(1 - 9 * ζ)
+
+function phi(uf::Beljaars, ζ, tt::MomentumTransport)
+    FT = eltype(uf)
+    if 0 < ζ
+        # Derived from Beljaars1991 Eq. 28
+        _a_m = FT(a_m(uf))
+        _b_m = FT(b_m(uf))
+        _c_m = FT(c_m(uf))
+        _d_m = FT(d_m(uf))
+        _π_group = FT(π_group(uf, tt))
+
+        first_exp = _b_m * exp(-_d_m * ζ)
+        second_exp = -_b_m * _d_m * exp(-_d_m * ζ) * (ζ - _c_m / _d_m)
+
+        return _π_group + ζ * (_a_m + first_exp + second_exp)
+    else
+        # Nishizawa2018 Eq. A1 (L < 0)
+        f_m = f_momentum(uf, ζ)
+        return 1 / f_m
+    end
+end
+
+function phi(uf::Beljaars, ζ, tt::HeatTransport)
+    FT = eltype(uf)
+    if 0 < ζ
+        # Derived from Beljaars1991 Eq. 32
+        _a_h = FT(a_h(uf))
+        _b_h = FT(b_h(uf))
+        _c_h = FT(c_h(uf))
+        _d_h = FT(d_h(uf))
+        _π_group = FT(π_group(uf, tt))
+
+        first_exp = _b_h * exp(-_d_h * ζ)
+        second_exp = -_b_h * _d_h * exp(-_d_h * ζ) * (ζ - _c_h / _d_h)
+        sqrt_term = _a_h * sqrt(2 * _a_h * ζ / 3 + 1)
+
+        return _π_group + ζ * (first_exp + second_exp + sqrt_term)
+    else
+        # Nishizawa2018 Eq. A2 (L < 0)
+        f_h = f_heat(uf, ζ)
+        return 1 / f_h
+    end
+end
+
+function psi(uf::Beljaars, ζ, tt::MomentumTransport)
+    FT = eltype(uf)
+    if 0 < ζ
+        # Beljaars1991 Eq. 28
+        _a_m = FT(a_m(uf))
+        _b_m = FT(b_m(uf))
+        _c_m = FT(c_m(uf))
+        _d_m = FT(d_m(uf))
+
+        return -_a_m * ζ - _b_m * (ζ - _c_m / _d_m) * exp(-_d_m * ζ) - _b_m * _c_m / _d_m
+    else
+        # Nishizawa2018 Eq. A3 (ζ < 0)
+        f_m = f_momentum(uf, ζ)
+        log_term = log((1 + f_m)^2 * (1 + f_m^2) / 8)
+        return log_term - 2 * atan(f_m) + FT(π) / 2
+    end
+end
+
+function psi(uf::Beljaars, ζ, tt::HeatTransport)
+    FT = eltype(uf)
+    if 0 < ζ
+        # Beljaars1991 Eq. 32
+        _a_h = FT(a_h(uf))
+        _b_h = FT(b_h(uf))
+        _c_h = FT(c_h(uf))
+        _d_h = FT(d_h(uf))
+
+        power_term = -(1 + 2 * _a_h * ζ / 3)^(FT(3 / 2))
+        exp_term = -_b_h * (ζ - _c_h / _d_h) * exp(-_d_h * ζ)
+        const_term = -_b_h * _c_h / _d_h
+
+        return power_term + 1 + exp_term + const_term
+
+    else
+        # Nishizawa2018 Eq. A4 (ζ < 0)
+        f_h = f_heat(uf, ζ)
+        return 2 * log((1 + f_h) / 2)
+    end
+end
+
+function Psi(uf::Beljaars, ζ, tt::MomentumTransport)
+    FT = eltype(uf)
+    _a_m = FT(a_m(uf))
+    _b_m = FT(b_m(uf))
+    _c_m = FT(c_m(uf))
+    _d_m = FT(d_m(uf))
+    if (abs(ζ) < eps(FT) && ζ >= 0) || (abs(ζ) >= eps(FT) && ζ < 0)
+        # Volume-averaged form of Beljaars1991 Eq. 28
+        exp_term1 = _c_m * ((1 - exp(-_d_m * ζ)) / (ζ * _d_m^2) - 1 / _d_m)
+        exp_term2 = (exp(-_d_m * ζ) - 1) / (ζ * _d_m^2)
+        exp_term3 = exp(-_d_m * ζ) / _d_m
+        return _b_m * (exp_term1 + exp_term2 + exp_term3) - _a_m * ζ / 2
+    elseif (abs(ζ) < eps(FT) && ζ < 0)
+        # Nishizawa2018 Eq. A13 (ζ < 0)
+        return -FT(15) * ζ / FT(8)
+    else
+        # Nishizawa2018 Eq. A5 (ζ < 0)
+        return -FT(9) * _a_m * ((_b_m * ζ + FT(1))^(FT(4 / 3)) - 1) / ζ / FT(4) / _b_m^FT(2) - FT(1)
+    end
+end
+
+function Psi(uf::Beljaars, ζ, tt::HeatTransport)
+    FT = eltype(uf)
+    _a_h = FT(a_h(uf))
+    _b_h = FT(b_h(uf))
+    _c_h = FT(c_h(uf))
+    _d_h = FT(d_h(uf))
+    if ζ >= 0
+        # Volume-averaged form of Beljaars1991 Eq. 32
+        sqrt_term = ((3^(FT(1 / 2))) * (2 * _a_h * ζ + 3)^(FT(5 / 2)) - 27) / (45 * _a_h)
+        exp_term1 = _b_h * _c_h * (exp(-_d_h * ζ) - 1) / (_d_h^2)
+        linear_term = ζ * (_b_h * _c_h / _d_h - 1)
+        exp_term2 = _b_h * (1 - exp(-_d_h * ζ) * (_d_h * ζ + 1)) / (_d_h^2)
+        return (-1 / ζ) * (sqrt_term + exp_term1 + linear_term + exp_term2)
+    elseif (abs(ζ) < eps(FT) && ζ < 0)
+        # Nishizawa2018 Eq. A14 (ζ < 0)
+        return -9 * ζ / 4
+    else
+        # Nishizawa2018 Eq. A6 (ζ < 0)
+        f_h = f_heat(uf, ζ)
+        log_term = 2 * log((1 + f_h) / 2)
+        return log_term + 2 * (1 - f_h) / (9 * ζ) - 1
+    end
+end
+
+
 universal_func(::BusingerType, L_MO::Real, params::BusingerParams) = Businger(L_MO, params)
 universal_func(::GryanikType, L_MO::Real, params::GryanikParams) = Gryanik(L_MO, params)
 universal_func(::GrachevType, L_MO::Real, params::GrachevParams) = Grachev(L_MO, params)
+universal_func(::BeljaarsType, L_MO::Real, params::BeljaarsParams) = Beljaars(L_MO, params)
 
 universal_func_type(::Type{T}) where {T <: BusingerParams} = BusingerType()
 universal_func_type(::Type{T}) where {T <: GryanikParams} = GryanikType()
 universal_func_type(::Type{T}) where {T <: GrachevParams} = GrachevType()
+universal_func_type(::Type{T}) where {T <: BeljaarsParams} = Beljaars()
 
 end # module
diff --git a/test/test_convergence.jl b/test/test_convergence.jl
index b19e5c6a..43871847 100644
--- a/test/test_convergence.jl
+++ b/test/test_convergence.jl
@@ -212,10 +212,10 @@ function check_over_moist_states(
 end
 
 @testset "Check convergence (dry thermodynamic states): Stable/Unstable" begin
-    for uf_type in [UF.BusingerType(), UF.GryanikType(), UF.GrachevType()]
+    for uf_type in [UF.BusingerType(), UF.GryanikType(), UF.GrachevType(), UF.BeljaarsType()]
         for FT in [Float32, Float64]
             for profile_type in [DryProfiles(), MoistEquilProfiles()]
-                profiles_sfc, profiles_int = generate_profiles(FT, profile_type)
+                profiles_sfc, profiles_int = generate_profiles(FT, profile_type; uf_type = uf_type)
                 scheme = [SF.FVScheme(), SF.FDScheme()]
                 z0_momentum = Array{FT}(range(1e-6, stop = 1e-1, length = 2))
                 z0_thermal = Array{FT}(range(1e-6, stop = 1e-1, length = 2))
@@ -250,6 +250,6 @@ end
             end
         end
     end
-    @info "Tested Businger/Gryanik/Grachev u.f., Float32/Float64, Dry/EquilMoist profiles, non-iterative solver on/off"
+    @info "Tested Businger/Gryanik/Grachev/Beljaars u.f., Float32/Float64, Dry/EquilMoist profiles, non-iterative solver on/off"
 
 end
diff --git a/test/test_universal_functions.jl b/test/test_universal_functions.jl
index 2f08a45e..830efeec 100644
--- a/test/test_universal_functions.jl
+++ b/test/test_universal_functions.jl
@@ -23,7 +23,7 @@ universal_functions(uft, L) = UF.universal_func(uft, L, create_uf_parameters(tom
         FT = Float32
         ζ = FT(-2):FT(0.01):FT(200)
         for L in (-FT(10), FT(10))
-            for uft in (UF.GryanikType(), UF.GrachevType(), UF.BusingerType())
+            for uft in (UF.GryanikType(), UF.GrachevType(), UF.BusingerType(), UF.BeljaarsType())
                 uf = universal_functions(uft, L)
                 for transport in (UF.MomentumTransport(), UF.HeatTransport())
                     ϕ = UF.phi.(uf, ζ, transport)
@@ -38,7 +38,7 @@ universal_functions(uft, L) = UF.universal_func(uft, L, create_uf_parameters(tom
         FT = Float32
         ζ = (-FT(1), FT(0.5) * eps(FT), 2 * eps(FT))
         for L in (-FT(10), FT(10))
-            for uft in (UF.GryanikType(), UF.GrachevType(), UF.BusingerType())
+            for uft in (UF.GryanikType(), UF.GrachevType(), UF.BusingerType(), UF.BeljaarsType())
                 uf = universal_functions(uft, L)
                 for transport in (UF.MomentumTransport(), UF.HeatTransport())
                     ϕ = UF.phi.(uf, ζ, transport)
@@ -53,11 +53,12 @@ universal_functions(uft, L) = UF.universal_func(uft, L, create_uf_parameters(tom
         FT = Float32
         ζ = (-FT(1), -FT(0.5) * eps(FT), FT(0.5) * eps(FT), 2 * eps(FT))
         for L in (-FT(10), FT(10))
-            uft = UF.BusingerType()
-            uf = universal_functions(uft, L)
-            for transport in (UF.MomentumTransport(), UF.HeatTransport())
-                Ψ = UF.Psi.(uf, ζ, transport)
-                @test eltype(Ψ) == FT
+            for uft in (UF.GryanikType(), UF.BusingerType(), UF.BeljaarsType())
+                uf = universal_functions(uft, L)
+                for transport in (UF.MomentumTransport(), UF.HeatTransport())
+                    Ψ = UF.Psi.(uf, ζ, transport)
+                    @test eltype(Ψ) == FT
+                end
             end
         end
     end