Merge pull request #3928 from CliMA/zs/updraft_rayleigh_sponge

add rayleigh sponge for edmf variables

Author: szy21

Project.toml

@@ -47,13 +47,13 @@ ClimaComms = "0.6.9"
 ClimaCore = "0.14.37"
 ClimaDiagnostics = "0.2.12"
 ClimaInterpolations = "0.1.0"
-ClimaParams = "0.10.35, 0.11"
+ClimaParams = "0.11.1"
 ClimaTimeSteppers = "0.8.2"
 ClimaUtilities = "0.1.23"
-CloudMicrophysics = "0.26.1"
+CloudMicrophysics = "0.26.2"
 Dates = "1"
-ForwardDiff = "0.10.38, 1"
-Insolation = "0.9.2"
+ForwardDiff = "1"
+Insolation = "0.9.5"
 Interpolations = "0.15.1, 0.16"
 LazyArtifacts = "1"
 LazyBroadcast = "1"
@@ -62,14 +62,14 @@ Logging = "1"
 NCDatasets = "0.14.2"
 NVTX = "0.3, 1"
 NullBroadcasts = "0.1"
-RRTMGP = "0.21.3"
+RRTMGP = "0.21.4"
 Random = "1"
-SciMLBase = "2.34"
+SciMLBase = "2.108"
 SparseMatrixColorings = "0.4.20"
 StaticArrays = "1.9"
 Statistics = "1"
-SurfaceFluxes = "0.11, 0.12"
-Thermodynamics = "0.12.14, 0.13"
+SurfaceFluxes = "0.12.1"
+Thermodynamics = "0.13.1"
 UnrolledUtilities = "0.1.9"
 YAML = "0.4"
 julia = "1.9"

config/model_configs/prognostic_edmfx_aquaplanet.yml

@@ -2,8 +2,6 @@
 z_max: 30000.0
 z_elem: 10
 dz_bottom: 500.0
-rayleigh_sponge: false
-viscous_sponge: false
 dt: 20secs
 surface_setup: DefaultMoninObukhov
 rad: clearsky

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-253
+254
 
 # **README**
 #
@@ -20,6 +20,12 @@
 
 
 #=
+254
+- Use rayleigh and viscous sponges in the prognostic_edmfx_aquaplanet
+
+253
+- Update deps, specifically ClimaParams.jl and Thermodynamics.jl
+
 252
 - Update deps, specifically RootSolvers.jl
 

src/parameterized_tendencies/sponge/rayleigh_sponge.jl

@@ -11,10 +11,26 @@ import ClimaCore.Fields as Fields
     αₘ(s, z, s.α_uₕ) * ζ_rayleigh(s, z, zmax)
 β_rayleigh_w(s::RayleighSponge{FT}, z, zmax) where {FT} =
     αₘ(s, z, s.α_w) * ζ_rayleigh(s, z, zmax)
+β_rayleigh_sgs_tracer(s::RayleighSponge{FT}, z, zmax) where {FT} =
+    αₘ(s, z, s.α_sgs_tracer) * ζ_rayleigh(s, z, zmax)
 
 function rayleigh_sponge_tendency_uₕ(ᶜuₕ, s)
     s isa Nothing && return NullBroadcasted()
     (; ᶜz, ᶠz) = z_coordinate_fields(axes(ᶜuₕ))
     zmax = z_max(axes(ᶠz))
     return @. lazy(-β_rayleigh_uₕ(s, ᶜz, zmax) * ᶜuₕ)
 end
+
+function rayleigh_sponge_tendency_sgs_tracer(ᶜχ, s)
+    s isa Nothing && return NullBroadcasted()
+    (; ᶜz, ᶠz) = z_coordinate_fields(axes(ᶜχ))
+    zmax = z_max(axes(ᶠz))
+    return @. lazy(-β_rayleigh_sgs_tracer(s, ᶜz, zmax) * ᶜχ)
+end
+
+function rayleigh_sponge_tendency_sgs_tracer(ᶜχʲ, ᶜχ, s)
+    s isa Nothing && return NullBroadcasted()
+    (; ᶜz, ᶠz) = z_coordinate_fields(axes(ᶜχ))
+    zmax = z_max(axes(ᶠz))
+    return @. lazy(-β_rayleigh_sgs_tracer(s, ᶜz, zmax) * (ᶜχʲ - ᶜχ))
+end

src/parameters/Parameters.jl

@@ -104,6 +104,7 @@ Base.@kwdef struct ClimaAtmosParameters{
     # Sponge
     alpha_rayleigh_w::FT
     alpha_rayleigh_uh::FT
+    alpha_rayleigh_sgs_tracer::FT
     zd_viscous::FT
     zd_rayleigh::FT
     kappa_2_sponge::FT

src/parameters/create_parameters.jl

@@ -109,6 +109,7 @@ atmos_name_map = (;
     :c_smag => :c_smag,
     :alpha_rayleigh_w => :alpha_rayleigh_w,
     :alpha_rayleigh_uh => :alpha_rayleigh_uh,
+    :alpha_rayleigh_sgs_tracer => :alpha_rayleigh_sgs_tracer,
     :astronomical_unit => :astro_unit,
     :held_suarez_T_equator_dry => :T_equator_dry,
     :drag_layer_vertical_extent => :σ_b,

src/prognostic_equations/remaining_tendency.jl

@@ -139,12 +139,13 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     ᶜuₕ = Y.c.uₕ
     ᶠu₃ = Y.f.u₃
     ᶜρ = Y.c.ρ
-    (; radiation_mode, moisture_model, turbconv_model) = p.atmos
+    (; radiation_mode, moisture_model, turbconv_model, microphysics_model) =
+        p.atmos
     (; rayleigh_sponge, viscous_sponge) = p.atmos
     (; ls_adv, scm_coriolis) = p.atmos
     (; params) = p
     thermo_params = CAP.thermodynamics_params(params)
-    (; ᶜp, sfc_conditions, ᶜts) = p.precomputed
+    (; ᶜp, sfc_conditions, ᶜts, ᶜK) = p.precomputed
 
     ᶜh_tot = @. lazy(
         TD.total_specific_enthalpy(
@@ -157,6 +158,55 @@ NVTX.@annotate function additional_tendency!(Yₜ, Y, p, t)
     vst_u₃ = viscous_sponge_tendency_u₃(ᶠu₃, viscous_sponge)
     vst_ρe_tot = viscous_sponge_tendency_ρe_tot(ᶜρ, ᶜh_tot, viscous_sponge)
     rst_uₕ = rayleigh_sponge_tendency_uₕ(ᶜuₕ, rayleigh_sponge)
+
+    if use_prognostic_tke(turbconv_model)
+        rst_ρatke =
+            rayleigh_sponge_tendency_sgs_tracer(Y.c.sgs⁰.ρatke, rayleigh_sponge)
+        @. Yₜ.c.sgs⁰.ρatke += rst_ρatke
+    end
+    if turbconv_model isa PrognosticEDMFX
+        ᶜmse = @. lazy(ᶜh_tot - ᶜK)
+        ᶜq_tot = @. lazy(specific(Y.c.ρq_tot, Y.c.ρ))
+        n = n_mass_flux_subdomains(p.atmos.turbconv_model)
+        for j in 1:n
+            rst_sgs_mse = rayleigh_sponge_tendency_sgs_tracer(
+                Y.c.sgsʲs.:($j).mse,
+                ᶜmse,
+                rayleigh_sponge,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).mse += rst_sgs_mse
+            rst_sgs_q_tot = rayleigh_sponge_tendency_sgs_tracer(
+                Y.c.sgsʲs.:($j).q_tot,
+                ᶜq_tot,
+                rayleigh_sponge,
+            )
+            @. Yₜ.c.sgsʲs.:($$j).q_tot += rst_sgs_q_tot
+        end
+        if moisture_model isa NonEquilMoistModel &&
+           microphysics_model isa Microphysics1Moment
+            # TODO: This doesn't work for multiple updrafts
+            moisture_species = (
+                (@name(c.sgsʲs.:(1).q_liq), @name(c.ρq_liq)),
+                (@name(c.sgsʲs.:(1).q_ice), @name(c.ρq_ice)),
+                (@name(c.sgsʲs.:(1).q_rai), @name(c.ρq_rai)),
+                (@name(c.sgsʲs.:(1).q_sno), @name(c.ρq_sno)),
+            )
+            MatrixFields.unrolled_foreach(
+                moisture_species,
+            ) do (sgs_q_name, ρq_name)
+                ᶜρq = MatrixFields.get_field(Y, ρq_name)
+                ᶜq = @. lazy(specific(ᶜρq, Y.c.ρ))
+                ᶜsgs_q = MatrixFields.get_field(Y, sgs_q_name)
+                ᶜsgs_qₜ = MatrixFields.get_field(Yₜ, sgs_q_name)
+                rst_sgs_q = rayleigh_sponge_tendency_sgs_tracer(
+                    ᶜsgs_q,
+                    ᶜq,
+                    rayleigh_sponge,
+                )
+                @. ᶜsgs_qₜ += rst_sgs_q
+            end
+        end
+    end
     # For HeldSuarezForcing, the radiation_mode is used as the forcing parameter
     forcing = radiation_mode isa HeldSuarezForcing ? radiation_mode : nothing
     hs_args = (ᶜuₕ, ᶜp, params, sfc_conditions.ts, moisture_model, forcing)

src/solver/model_getters.jl

@@ -183,7 +183,8 @@ function get_rayleigh_sponge_model(parsed_args, params, ::Type{FT}) where {FT}
         zd = params.zd_rayleigh
         α_uₕ = params.alpha_rayleigh_uh
         α_w = params.alpha_rayleigh_w
-        RayleighSponge{FT}(; zd, α_uₕ, α_w)
+        α_sgs_tracer = params.alpha_rayleigh_sgs_tracer
+        RayleighSponge{FT}(; zd, α_uₕ, α_w, α_sgs_tracer)
     else
         error("Uncaught rayleigh sponge model `$rs_name`.")
     end

src/solver/types.jl

@@ -238,6 +238,7 @@ Base.@kwdef struct RayleighSponge{FT} <: AbstractSponge
     zd::FT
     α_uₕ::FT
     α_w::FT
+    α_sgs_tracer::FT
 end
 
 abstract type AbstractGravityWave end

test/parameterized_tendencies/sponge/rayleigh_sponge.jl

@@ -35,7 +35,12 @@ using ClimaCore.CommonSpaces
     @. ᶜuₕ.components.data.:1 = 1
     @. ᶜuₕ.components.data.:2 = 1
     ### Component test begins here
-    rs = CA.RayleighSponge(; zd = FT(0), α_uₕ = FT(1), α_w = FT(1))
+    rs = CA.RayleighSponge(;
+        zd = FT(0),
+        α_uₕ = FT(1),
+        α_w = FT(1),
+        α_sgs_tracer = FT(1),
+    )
     expected = @. sin(FT(π) / 2 * ᶜz / zmax)^2
     computed = CA.rayleigh_sponge_tendency_uₕ(ᶜuₕ, rs)
     @test CA.β_rayleigh_uₕ.(rs, ᶜz, zmax) == expected