diff --git a/docs/bibliography.bib b/docs/bibliography.bib
index 0caca54d..78889b95 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{Holtslag1988,
+  title = {Applied Modeling of the Nighttime Surface Energy Balance over Land},
+  author = {Holtslag, A. A. M., and H. A. R. De Bruin},
+  journal = {Journal of Applied Meteorology and Climatology},
+  volume = {27},
+  number = {6},
+  doi = {10.1175/1520-0450(1988)027<0689:AMOTNS>2.0.CO;2},
+  pages = {689–704},
+  year = {1988},
+  publisher = {Wiley Online Library}
+}
+
 @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..dcfb3456 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.Holtslag
 ```
 
 ```@docs
diff --git a/docs/src/plot_universal_functions.jl b/docs/src/plot_universal_functions.jl
index 5173898a..84e4a249 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.HoltslagType())
 
 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..2f10e19d 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.HoltslagType)
+    FT = CP.float_type(toml_dict)
+
+    aliases = [
+        "Pr_0_Holtslag",
+        "a_m_Holtslag",
+        "a_h_Holtslag",
+        "b_m_Holtslag",
+        "b_h_Holtslag",
+        "c_m_Holtslag",
+        "c_h_Holtslag",
+        "d_m_Holtslag",
+        "d_h_Holtslag",
+        "ζ_a_Holtslag",
+        "γ_Holtslag",
+    ]
+
+    pairs = CP.get_parameter_values!(toml_dict, aliases, "UniversalFunctions")
+    pairs = (; pairs...) # convert to NamedTuple
+
+    pairs = (;
+        Pr_0 = pairs.Pr_0_Holtslag,
+        a_m = pairs.a_m_Holtslag,
+        a_h = pairs.a_h_Holtslag,
+        b_m = pairs.b_m_Holtslag,
+        b_h = pairs.b_h_Holtslag,
+        c_m = pairs.c_m_Holtslag,
+        c_h = pairs.c_h_Holtslag,
+        d_m = pairs.d_m_Holtslag,
+        d_h = pairs.d_h_Holtslag,
+        ζ_a = pairs.ζ_a_Holtslag,
+        γ = pairs.γ_Holtslag,
+    )
+    return UF.HoltslagParams{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..8a145413 100644
--- a/src/UniversalFunctions.jl
+++ b/src/UniversalFunctions.jl
@@ -7,6 +7,7 @@ universal functions:
  - `Businger`
  - `Gryanik`
  - `Grachev`
+ - `Holtslag`
 """
 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,194 @@ function psi(uf::Grachev, ζ, tt::HeatTransport)
     end
 end
 
+#####
+##### Holtslag
+#####
+
+Base.@kwdef struct HoltslagParams{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
+
+"""
+    Holtslag <: AbstractUniversalFunction{FT}
+
+# References
+ - [Holtslag1988](@cite)
+
+# Equations in reference:
+
+    `ϕ_m`: Derived from Eq. 12
+    `ϕ_h`: Derived from Eq. 12
+    `ψ_m`: Eq. 12
+    `ψ_h`: Eq. 12
+
+# Holtslag and Bruin 1988 functions are used in stable conditions
+# In unstable conditions the functions of Businger (1971) are 
+# assigned by default. 
+
+# Fields
+
+$(DSE.FIELDS)
+"""
+struct Holtslag{FT, PS <: HoltslagParams} <: AbstractUniversalFunction{FT}
+    "Monin-Obhukov Length"
+    L::FT
+    params::PS
+end
+
+struct HoltslagType <: AbstractUniversalFunctionType end
+Holtslag() = HoltslagType()
+
+# Nishizawa2018 Eq. A7
+f_momentum(uf::Holtslag, ζ) = sqrt(sqrt(1 - 15 * ζ))
+
+# Nishizawa2018 Eq. A8
+f_heat(uf::Holtslag, ζ) = sqrt(1 - 9 * ζ)
+
+function phi(uf::Holtslag, ζ, tt::MomentumTransport)
+    FT = eltype(uf)
+    if 0 < ζ
+        # Derived from Holtslag1988 Eq. 12
+        _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::Holtslag, ζ, tt::HeatTransport)
+    FT = eltype(uf)
+    if 0 < ζ
+        # Derived from Holtslag1988 Eq. 12
+        _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)
+
+        return _π_group + ζ * (_a_h + first_exp + second_exp)
+    else
+        # Nishizawa2018 Eq. A2 (L < 0)
+        f_h = f_heat(uf, ζ)
+        return 1 / f_h
+    end
+end
+
+function psi(uf::Holtslag, ζ, tt::MomentumTransport)
+    FT = eltype(uf)
+    if 0 < ζ
+        # Holtslag1988 Eq. 12
+        _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::Holtslag, ζ, tt::HeatTransport)
+    FT = eltype(uf)
+    if 0 < ζ
+        # Holtslag1988 Eq. 12
+        _a_h = FT(a_h(uf))
+        _b_h = FT(b_h(uf))
+        _c_h = FT(c_h(uf))
+        _d_h = FT(d_h(uf))
+
+        return -_a_h * ζ - _b_h * (ζ - _c_h / _d_h) * exp(-_d_h * ζ) - _b_h * _c_h / _d_h
+
+    else
+        # Nishizawa2018 Eq. A4 (ζ < 0)
+        f_h = f_heat(uf, ζ)
+        return 2 * log((1 + f_h) / 2)
+    end
+end
+
+function Psi(uf::Holtslag, ζ, 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 ζ >= 0
+        # Derived from Holtslag1988 Eq. 12
+        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::Holtslag, ζ, 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
+        # Derived from Holtslag1988 Eq. 12
+        exp_term1 = _c_h * ((1 - exp(-_d_h * ζ)) / (ζ * _d_h^2) - 1 / _d_h)
+        exp_term2 = (exp(-_d_h * ζ) - 1) / (ζ * _d_h^2)
+        exp_term3 = exp(-_d_h * ζ) / _d_h
+        return _b_h * (exp_term1 + exp_term2 + exp_term3) - _a_h * ζ / 2
+
+    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(::HoltslagType, L_MO::Real, params::HoltslagParams) = Holtslag(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 <: HoltslagParams} = HoltslagType()
 
 end # module
diff --git a/test/test_convergence.jl b/test/test_convergence.jl
index b19e5c6a..a7dfa61a 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.HoltslagType()]
         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/Holtslag 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..53ce66c5 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.HoltslagType())
                 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.HoltslagType())
                 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.HoltslagType())
+                uf = universal_functions(uft, L)
+                for transport in (UF.MomentumTransport(), UF.HeatTransport())
+                    Ψ = UF.Psi.(uf, ζ, transport)
+                    @test eltype(Ψ) == FT
+                end
             end
         end
     end