add new file

CliMA · Sep 18, 2023 · 22bc342 · 22bc342
1 parent b7b4601
commit 22bc342
Show file tree

Hide file tree

Showing 3 changed files with 283 additions and 279 deletions.
diff --git a/src/ClimaAtmos.jl b/src/ClimaAtmos.jl
@@ -80,6 +80,7 @@ include(joinpath("prognostic_equations", "hyperdiffusion.jl"))
 include(joinpath("prognostic_equations", "edmf_coriolis.jl"))
 include(joinpath("prognostic_equations", "buoyancy_gradients.jl"))
 include(joinpath("prognostic_equations", "edmfx_closures.jl"))
+include(joinpath("prognostic_equations", "edmfx_entr_detr.jl"))
 include(joinpath("prognostic_equations", "edmfx_tke.jl"))
 include(joinpath("prognostic_equations", "edmfx_sgs_flux.jl"))
 include(joinpath("prognostic_equations", "edmfx_boundary_condition.jl"))

diff --git a/src/prognostic_equations/edmfx_closures.jl b/src/prognostic_equations/edmfx_closures.jl
@@ -1,5 +1,5 @@
 #####
-##### EDMF entrainment detrainment
+##### EDMF closures (nonhydrostatic pressure drag and mixing length)
 #####
 
 import StaticArrays as SA
@@ -109,284 +109,6 @@ function edmfx_nh_pressure_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMF
     end
 end
 
-"""
-   Return entrainment rate [1/s].
-
-   Inputs (everything defined on cell centers):
-   - params set with model parameters
-   - ᶜz, z_sfc, ᶜp, ᶜρ, - grid-scale height, surface height, grid-scale pressure and density
-   - buoy_flux_surface - buoyancy flux at the surface
-   - ᶜaʲ, ᶜwʲ, ᶜRHʲ, ᶜbuoyʲ - updraft area, physical vertical velocity,
-                                   relative humidity and buoyancy
-   - ᶜw⁰, ᶜRH⁰, ᶜbuoy⁰ - environment physical vertical velocity,
-                              relative humidity and buoyancy
-   - dt - timestep
-"""
-
-function pi_groups_entr(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    dt::FT,
-    ::NoEntrainment,
-) where {FT}
-    return FT(0)
-end
-
-function pi_groups_entr(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    dt::FT,
-    ::PiGroupsEntrainment,
-) where {FT}
-    if ᶜaʲ <= FT(0)
-        return FT(0)
-    else
-        g = CAP.grav(params)
-        turbconv_params = CAP.turbconv_params(params)
-        ᶜaʲ_max = TCP.max_area(turbconv_params)
-
-        # pressure scale height (height where pressure drops by 1/e)
-        ref_H = ᶜp / (ᶜρ * g)
-        # convective velocity
-        w_star = get_wstar(buoy_flux_surface)
-
-        # non-dimensional pi-groups
-        # TODO - using Π₁ blows things up
-        Π₁ =
-            (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) /
-            ((ᶜwʲ - ᶜw⁰)^2 + w_star^2 + eps(FT))
-        Π₃ = sqrt(ᶜaʲ)
-        Π₄ = ᶜRHʲ - ᶜRH⁰
-        Π₆ = (ᶜz - z_sfc) / ref_H
-        entr = max(
-            0,
-            ᶜwʲ / (ᶜz - z_sfc) * (
-                -4.013288 - 0.000968 * Π₁ + 0.356974 * Π₃ - 0.403124 * Π₄ +
-                1.503261 * Π₆
-            ),
-        )
-
-        entr = max(0, min(entr, 1 / dt))
-        return entr
-    end
-end
-
-function pi_groups_entr(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    dt::FT,
-    ::ConstantCoefficientEntrainment,
-) where {FT}
-    if ᶜaʲ <= FT(0)
-        return FT(0)
-    else
-        entr_coeff = CAP.entr_coeff(params)
-        entr = max(0, min(entr_coeff * abs(ᶜwʲ), 1 / dt))
-        return entr
-    end
-end
-
-function pi_groups_entr(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    dt::FT,
-    ::ConstantTimescaleEntrainment,
-) where {FT}
-    if ᶜaʲ <= FT(0)
-        return FT(0)
-    else
-        entr_tau = CAP.entr_tau(params)
-        entr = max(0, min(1 / entr_tau, 1 / dt))
-        return entr
-    end
-end
-
-"""
-   Return detrainment rate [1/s].
-
-   Inputs (everything defined on cell centers):
-   - params set with model parameters
-   - ᶜz, z_sfc, ᶜp, ᶜρ, - grid-scale height, surface height, grid-scale pressure and density
-   - buoy_flux_surface - buoyancy flux at the surface
-   - ᶜaʲ, ᶜwʲ, ᶜRHʲ, ᶜbuoyʲ - updraft area, physical vertical velocity,
-                                   relative humidity and buoyancy
-   - ᶜw⁰, ᶜRH⁰, ᶜbuoy⁰ - environment physical vertical velocity,
-                              relative humidity and buoyancy
-   - dt - timestep
-"""
-function pi_groups_detr(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    dt::FT,
-    ::NoDetrainment,
-) where {FT}
-    return FT(0)
-end
-
-function pi_groups_detr(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    dt::FT,
-    ::PiGroupsDetrainment,
-) where {FT}
-
-    if ᶜaʲ <= FT(0)
-        return FT(0)
-    else
-        g = CAP.grav(params)
-        turbconv_params = CAP.turbconv_params(params)
-        ᶜaʲ_max = TCP.max_area(turbconv_params)
-        max_area_limiter = 0.1 * exp(-10 * (ᶜaʲ_max - ᶜaʲ))
-
-        # pressure scale height (height where pressure drops by 1/e)
-        ref_H = ᶜp / (ᶜρ * g)
-        # convective velocity
-        w_star = get_wstar(buoy_flux_surface)
-
-        # non-dimensional pi-groups
-        # TODO - using Π₁ blows things up
-        Π₁ =
-            (ᶜz - z_sfc) * (ᶜbuoyʲ - ᶜbuoy⁰) /
-            ((ᶜwʲ - ᶜw⁰)^2 + w_star^2 + eps(FT))
-        Π₃ = sqrt(ᶜaʲ)
-        Π₄ = ᶜRHʲ - ᶜRH⁰
-        Π₆ = (ᶜz - z_sfc) / ref_H
-        detr = max(
-            0,
-            ᶜwʲ * (
-                3.535208 + 0.598496 * Π₁ + 1.583348 * Π₃ + 0.046275 * Π₄ -
-                0.344836 * Π₆ + max_area_limiter
-            ),
-        )
-        return detr
-    end
-end
-
-function pi_groups_detr(
-    params,
-    ᶜz::FT,
-    z_sfc::FT,
-    ᶜp::FT,
-    ᶜρ::FT,
-    buoy_flux_surface::FT,
-    ᶜaʲ::FT,
-    ᶜwʲ::FT,
-    ᶜRHʲ::FT,
-    ᶜbuoyʲ::FT,
-    ᶜw⁰::FT,
-    ᶜRH⁰::FT,
-    ᶜbuoy⁰::FT,
-    dt::FT,
-    ::ConstantCoefficientDetrainment,
-) where {FT}
-
-    if ᶜaʲ <= FT(0)
-        return FT(0)
-    else
-        detr_coeff = CAP.detr_coeff(params)
-        detr = max(0, min(detr_coeff * abs(ᶜwʲ), 1 / dt))
-        return detr
-    end
-end
-
-edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model) = nothing
-function edmfx_entr_detr_tendency!(Yₜ, Y, p, t, colidx, turbconv_model::EDMFX)
-
-    n = n_mass_flux_subdomains(turbconv_model)
-    (; ᶜspecificʲs, ᶜh_totʲs, ᶜentrʲs, ᶜdetrʲs) = p
-    (; ᶜu⁰, ᶜspecific⁰, ᶜh_tot⁰) = p
-
-    for j in 1:n
-
-        @. Yₜ.c.sgsʲs.:($$j).ρa[colidx] +=
-            Y.c.sgsʲs.:($$j).ρa[colidx] *
-            (ᶜentrʲs.:($$j)[colidx] - ᶜdetrʲs.:($$j)[colidx])
-
-        @. Yₜ.c.sgsʲs.:($$j).ρae_tot[colidx] +=
-            Y.c.sgsʲs.:($$j).ρa[colidx] * (
-                ᶜentrʲs.:($$j)[colidx] * ᶜh_tot⁰[colidx] -
-                ᶜdetrʲs.:($$j)[colidx] * ᶜh_totʲs.:($$j)[colidx]
-            )
-
-        @. Yₜ.c.sgsʲs.:($$j).ρaq_tot[colidx] +=
-            Y.c.sgsʲs.:($$j).ρa[colidx] * (
-                ᶜentrʲs.:($$j)[colidx] * ᶜspecific⁰.q_tot[colidx] -
-                ᶜdetrʲs.:($$j)[colidx] * ᶜspecificʲs.:($$j).q_tot[colidx]
-            )
-
-        @. Yₜ.f.sgsʲs.:($$j).u₃[colidx] +=
-            ᶠinterp(ᶜentrʲs.:($$j)[colidx] * C3(ᶜu⁰[colidx])) -
-            ᶠinterp(ᶜentrʲs.:($$j)[colidx]) * Y.f.sgsʲs.:($$j).u₃[colidx]
-    end
-    return nothing
-end
-
 # lambert_2_over_e(::Type{FT}) where {FT} = FT(LambertW.lambertw(FT(2) / FT(MathConstants.e)))
 lambert_2_over_e(::Type{FT}) where {FT} = FT(0.46305551336554884) # since we can evaluate