Advection of water with working fluid changes

src/ClimaAtmos.jl

@@ -57,6 +57,7 @@ include(joinpath("prognostic_equations", "zero_tendency.jl"))
 include(joinpath("prognostic_equations", "implicit", "implicit_tendency.jl"))
 include(joinpath("prognostic_equations", "implicit", "implicit_solver.jl"))
 
+include(joinpath("prognostic_equations", "water_advection.jl"))
 include(joinpath("prognostic_equations", "remaining_tendency.jl"))
 include(joinpath("prognostic_equations", "forcing", "large_scale_advection.jl")) # TODO: should this be in tendencies/?
 include(joinpath("prognostic_equations", "forcing", "subsidence.jl"))

src/cache/diagnostic_edmf_precomputed_quantities.jl

@@ -1049,32 +1049,33 @@ NVTX.@annotate function set_diagnostic_edmf_precomputed_quantities_env_precipita
     t,
     precip_model::Microphysics1Moment,
 )
-    thermo_params = CAP.thermodynamics_params(p.params)
-    microphys_1m_params = CAP.microphysics_1m_params(p.params)
-
-    (; ᶜts, ᶜSqₜᵖ⁰, ᶜSeₜᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰) = p.precomputed
-    (; q_tot) = p.precomputed.ᶜspecific
-    (; ᶜqᵣ, ᶜqₛ) = p.precomputed
-
-    ᶜSᵖ = p.scratch.ᶜtemp_scalar
-    ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2
-
-    # Environment precipitation sources (to be applied to grid mean)
-    compute_precipitation_sources!(
-        ᶜSᵖ,
-        ᶜSᵖ_snow,
-        ᶜSqₜᵖ⁰,
-        ᶜSqᵣᵖ⁰,
-        ᶜSqₛᵖ⁰,
-        ᶜSeₜᵖ⁰,
-        Y.c.ρ,
-        ᶜqᵣ,
-        ᶜqₛ,
-        ᶜts,
-        p.core.ᶜΦ,
-        p.dt,
-        microphys_1m_params,
-        thermo_params,
-    )
+    error("Not implemented yet")
+    #thermo_params = CAP.thermodynamics_params(p.params)
+    #microphys_1m_params = CAP.microphysics_1m_params(p.params)
+
+    #(; ᶜts, ᶜSqₜᵖ⁰, ᶜSeₜᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰) = p.precomputed
+    #(; q_tot) = p.precomputed.ᶜspecific
+    #(; ᶜqᵣ, ᶜqₛ) = p.precomputed
+
+    #ᶜSᵖ = p.scratch.ᶜtemp_scalar
+    #ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2
+
+    ## Environment precipitation sources (to be applied to grid mean)
+    #compute_precipitation_sources!(
+    #    ᶜSᵖ,
+    #    ᶜSᵖ_snow,
+    #    ᶜSqₜᵖ⁰,
+    #    ᶜSqᵣᵖ⁰,
+    #    ᶜSqₛᵖ⁰,
+    #    ᶜSeₜᵖ⁰,
+    #    Y.c.ρ,
+    #    ᶜqᵣ,
+    #    ᶜqₛ,
+    #    ᶜts,
+    #    p.core.ᶜΦ,
+    #    p.dt,
+    #    microphys_1m_params,
+    #    thermo_params,
+    #)
     return nothing
 end

src/cache/precipitation_precomputed_quantities.jl

@@ -2,11 +2,7 @@
 ##### Precomputed quantities
 #####
 import CloudMicrophysics.Microphysics1M as CM1
-
-# helper function to safely get precipitation from state
-function qₚ(ρqₚ, ρ)
-    return max(zero(ρ), ρqₚ / ρ)
-end
+import CloudMicrophysics.MicrophysicsNonEq as CMNe
 
 """
     set_precipitation_precomputed_quantities!(Y, p, t)
@@ -17,46 +13,49 @@ for the 1-moment microphysics scheme
 function set_precipitation_precomputed_quantities!(Y, p, t)
     @assert (p.atmos.precip_model isa Microphysics1Moment)
 
-    (; ᶜwᵣ, ᶜwₛ, ᶜqᵣ, ᶜqₛ) = p.precomputed
-
+    (; ᶜwᵣ, ᶜwₛ) = p.precomputed
     cmp = CAP.microphysics_1m_params(p.params)
 
-    # compute the precipitation specific humidities
-    @. ᶜqᵣ = qₚ(Y.c.ρq_rai, Y.c.ρ)
-    @. ᶜqₛ = qₚ(Y.c.ρq_sno, Y.c.ρ)
-
     # compute the precipitation terminal velocity [m/s]
     @. ᶜwᵣ = CM1.terminal_velocity(
         cmp.pr,
         cmp.tv.rain,
         Y.c.ρ,
-        abs(Y.c.ρq_rai / Y.c.ρ),
+        max(zero(Y.c.ρ), Y.c.ρq_rai / Y.c.ρ),
     )
     @. ᶜwₛ = CM1.terminal_velocity(
         cmp.ps,
         cmp.tv.snow,
         Y.c.ρ,
-        abs(Y.c.ρq_sno / Y.c.ρ),
+        max(zero(Y.c.ρ), Y.c.ρq_sno / Y.c.ρ),
     )
     return nothing
 end
 
-
 """
     set_sedimentation_precomputed_quantities!(Y, p, t)
+
 Updates the sedimentation terminal velocity stored in `p`
 for the non-equilibrium microphysics scheme
 """
 function set_sedimentation_precomputed_quantities!(Y, p, t)
     @assert (p.atmos.moisture_model isa NonEquilMoistModel)
 
     (; ᶜwₗ, ᶜwᵢ) = p.precomputed
+    cmc = CAP.microphysics_cloud_params(p.params)
 
-    FT = eltype(Y)
-
-    # compute the precipitation terminal velocity [m/s]
-    # TODO - the actual parameterization will be added in the next PR
-    @. ᶜwₗ = FT(0)
-    @. ᶜwᵢ = FT(0)
+    # compute sedimentation velocity for cloud condensate [m/s]
+    @. ᶜwₗ = CMNe.terminal_velocity(
+        cmc.liquid,
+        cmc.Ch2022.rain,
+        Y.c.ρ,
+        max(zero(Y.c.ρ), Y.c.ρq_liq / Y.c.ρ),
+    )
+    @. ᶜwᵢ = CMNe.terminal_velocity(
+        cmc.ice,
+        cmc.Ch2022.small_ice,
+        Y.c.ρ,
+        max(zero(Y.c.ρ), Y.c.ρq_ice / Y.c.ρ),
+    )
     return nothing
 end

src/cache/precomputed_quantities.jl

@@ -154,12 +154,7 @@ function precomputed_quantities(Y, atmos)
         end
     precipitation_quantities =
         atmos.precip_model isa Microphysics1Moment ?
-        (;
-            ᶜwᵣ = similar(Y.c, FT),
-            ᶜwₛ = similar(Y.c, FT),
-            ᶜqᵣ = similar(Y.c, FT),
-            ᶜqₛ = similar(Y.c, FT),
-        ) : (;)
+        (; ᶜwᵣ = similar(Y.c, FT), ᶜwₛ = similar(Y.c, FT)) : (;)
     smagorinsky_lilly_quantities =
         if atmos.smagorinsky_lilly isa SmagorinskyLilly
             uvw_vec = UVW(FT(0), FT(0), FT(0))
@@ -441,27 +436,63 @@ NVTX.@annotate function set_precomputed_quantities!(Y, p, t)
     (; ᶜwₜqₜ, ᶜwₕhₜ) = p.precomputed
     @. ᶜwₜqₜ = Geometry.WVector(0)
     @. ᶜwₕhₜ = Geometry.WVector(0)
-    #
-    # TODO - uncomment in the next PR. Right now for the purpose of testing
-    # we want to merge with 0 sedimentation and precipitation
-    #
     if moisture_model isa NonEquilMoistModel
         set_sedimentation_precomputed_quantities!(Y, p, t)
-        #    (; ᶜwₗ, ᶜwᵢ) = p.precomputed
-        #    @. ᶜwₜqₜ += Geometry.WVector(ᶜwₗ * Y.c.ρq_liq + ᶜwᵢ * Y.c.ρq_ice) / Y.c.ρ
-        #    @. ᶜwₕhₜ += Geometry.WVector(
-        #        ᶜwₗ * Y.c.ρq_liq * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ + norm_sqr(Geometry.UVWVector(0, 0, -ᶜwₗ) + Geometry.UVWVector(ᶜu))/2) +
-        #        ᶜwᵢ * Y.c.ρq_ice * (TD.internal_energy_ice(thermo_params, ᶜts)    + ᶜΦ + norm_sqr(Geometry.UVWVector(0, 0, -ᶜwᵢ) + Geometry.UVWVector(ᶜu))/2)
-        #    ) / Y.c.ρ
+        (; ᶜwₗ, ᶜwᵢ) = p.precomputed
+        @. ᶜwₜqₜ +=
+            Geometry.WVector(ᶜwₗ * Y.c.ρq_liq + ᶜwᵢ * Y.c.ρq_ice) / Y.c.ρ
+        @. ᶜwₕhₜ +=
+            Geometry.WVector(
+                ᶜwₗ *
+                Y.c.ρq_liq *
+                (
+                    TD.internal_energy_liquid(thermo_params, ᶜts) +
+                    ᶜΦ +
+                    norm_sqr(
+                        Geometry.UVWVector(0, 0, -(ᶜwₗ)) +
+                        Geometry.UVWVector(ᶜu),
+                    ) / 2
+                ) +
+                ᶜwᵢ *
+                Y.c.ρq_ice *
+                (
+                    TD.internal_energy_ice(thermo_params, ᶜts) +
+                    ᶜΦ +
+                    norm_sqr(
+                        Geometry.UVWVector(0, 0, -(ᶜwᵢ)) +
+                        Geometry.UVWVector(ᶜu),
+                    ) / 2
+                ),
+            ) / Y.c.ρ
     end
     if precip_model isa Microphysics1Moment
         set_precipitation_precomputed_quantities!(Y, p, t)
-        #    (; ᶜwᵣ, ᶜwₛ) = p.precomputed
-        #    @. ᶜwₜqₜ += Geometry.WVector(ᶜwᵣ * Y.c.ρq_rai + ᶜwₛ * Y.c.ρq_sno) / Y.c.ρ
-        #    @. ᶜwₕhₜ += Geometry.WVector(
-        #        ᶜwᵣ * Y.c.ρq_rai * (TD.internal_energy_liquid(thermo_params, ᶜts) + ᶜΦ + norm_sqr(Geometry.UVWVector(0, 0, -ᶜwᵣ) + Geometry.UVWVector(ᶜu))/2) +
-        #        ᶜwₛ * Y.c.ρq_sno * (TD.internal_energy_ice(thermo_params, ᶜts)    + ᶜΦ + norm_sqr(Geometry.UVWVector(0, 0, -ᶜwₛ) + Geometry.UVWVector(ᶜu))/2)
-        #    ) / Y.c.ρ
+        (; ᶜwᵣ, ᶜwₛ) = p.precomputed
+        @. ᶜwₜqₜ +=
+            Geometry.WVector(ᶜwᵣ * Y.c.ρq_rai + ᶜwₛ * Y.c.ρq_sno) / Y.c.ρ
+        @. ᶜwₕhₜ +=
+            Geometry.WVector(
+                ᶜwᵣ *
+                Y.c.ρq_rai *
+                (
+                    TD.internal_energy_liquid(thermo_params, ᶜts) +
+                    ᶜΦ +
+                    norm_sqr(
+                        Geometry.UVWVector(0, 0, -(ᶜwᵣ)) +
+                        Geometry.UVWVector(ᶜu),
+                    ) / 2
+                ) +
+                ᶜwₛ *
+                Y.c.ρq_sno *
+                (
+                    TD.internal_energy_ice(thermo_params, ᶜts) +
+                    ᶜΦ +
+                    norm_sqr(
+                        Geometry.UVWVector(0, 0, -(ᶜwₛ)) +
+                        Geometry.UVWVector(ᶜu),
+                    ) / 2
+                ),
+            ) / Y.c.ρ
     end
 
     if turbconv_model isa PrognosticEDMFX

src/cache/prognostic_edmf_precomputed_quantities.jl

@@ -414,59 +414,60 @@ NVTX.@annotate function set_prognostic_edmf_precomputed_quantities_precipitation
     p,
     ::Microphysics1Moment,
 )
-    @assert !(p.atmos.moisture_model isa DryModel)
-
-    (; params, dt) = p
-    (; ᶜΦ,) = p.core
-    thp = CAP.thermodynamics_params(params)
-    cmp = CAP.microphysics_1m_params(params)
-
-    (; ᶜSeₜᵖʲs, ᶜSqₜᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs, ᶜρʲs, ᶜtsʲs) = p.precomputed
-    (; ᶜSeₜᵖ⁰, ᶜSqₜᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜρ⁰, ᶜts⁰) = p.precomputed
-    (; ᶜqᵣ, ᶜqₛ) = p.precomputed
-
-    # TODO - can I re-use them between js and env?
-    ᶜSᵖ = p.scratch.ᶜtemp_scalar
-    ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2
-
-    n = n_mass_flux_subdomains(p.atmos.turbconv_model)
-
-    # Sources from the updrafts
-    for j in 1:n
-        compute_precipitation_sources!(
-            ᶜSᵖ,
-            ᶜSᵖ_snow,
-            ᶜSqₜᵖʲs.:($j),
-            ᶜSqᵣᵖʲs.:($j),
-            ᶜSqₛᵖʲs.:($j),
-            ᶜSeₜᵖʲs.:($j),
-            ᶜρʲs.:($j),
-            ᶜqᵣ,
-            ᶜqₛ,
-            ᶜtsʲs.:($j),
-            ᶜΦ,
-            dt,
-            cmp,
-            thp,
-        )
-    end
-
-    # Sources from the environment
-    compute_precipitation_sources!(
-        ᶜSᵖ,
-        ᶜSᵖ_snow,
-        ᶜSqₜᵖ⁰,
-        ᶜSqᵣᵖ⁰,
-        ᶜSqₛᵖ⁰,
-        ᶜSeₜᵖ⁰,
-        ᶜρ⁰,
-        ᶜqᵣ,
-        ᶜqₛ,
-        ᶜts⁰,
-        ᶜΦ,
-        dt,
-        cmp,
-        thp,
-    )
+    error("Not implemented yet")
+    #@assert !(p.atmos.moisture_model isa DryModel)
+
+    #(; params, dt) = p
+    #(; ᶜΦ,) = p.core
+    #thp = CAP.thermodynamics_params(params)
+    #cmp = CAP.microphysics_1m_params(params)
+
+    #(; ᶜSeₜᵖʲs, ᶜSqₜᵖʲs, ᶜSqᵣᵖʲs, ᶜSqₛᵖʲs, ᶜρʲs, ᶜtsʲs) = p.precomputed
+    #(; ᶜSeₜᵖ⁰, ᶜSqₜᵖ⁰, ᶜSqᵣᵖ⁰, ᶜSqₛᵖ⁰, ᶜρ⁰, ᶜts⁰) = p.precomputed
+    #(; ᶜqᵣ, ᶜqₛ) = p.precomputed
+
+    ## TODO - can I re-use them between js and env?
+    #ᶜSᵖ = p.scratch.ᶜtemp_scalar
+    #ᶜSᵖ_snow = p.scratch.ᶜtemp_scalar_2
+
+    #n = n_mass_flux_subdomains(p.atmos.turbconv_model)
+
+    ## Sources from the updrafts
+    #for j in 1:n
+    #    compute_precipitation_sources!(
+    #        ᶜSᵖ,
+    #        ᶜSᵖ_snow,
+    #        ᶜSqₜᵖʲs.:($j),
+    #        ᶜSqᵣᵖʲs.:($j),
+    #        ᶜSqₛᵖʲs.:($j),
+    #        ᶜSeₜᵖʲs.:($j),
+    #        ᶜρʲs.:($j),
+    #        ᶜqᵣ,
+    #        ᶜqₛ,
+    #        ᶜtsʲs.:($j),
+    #        ᶜΦ,
+    #        dt,
+    #        cmp,
+    #        thp,
+    #    )
+    #end
+
+    ## Sources from the environment
+    #compute_precipitation_sources!(
+    #    ᶜSᵖ,
+    #    ᶜSᵖ_snow,
+    #    ᶜSqₜᵖ⁰,
+    #    ᶜSqᵣᵖ⁰,
+    #    ᶜSqₛᵖ⁰,
+    #    ᶜSeₜᵖ⁰,
+    #    ᶜρ⁰,
+    #    ᶜqᵣ,
+    #    ᶜqₛ,
+    #    ᶜts⁰,
+    #    ᶜΦ,
+    #    dt,
+    #    cmp,
+    #    thp,
+    #)
     return nothing
 end

src/parameters/create_parameters.jl

@@ -110,6 +110,7 @@ atmos_name_map = (;
 cloud_parameters(FT_or_toml) = (;
     liquid = CM.Parameters.CloudLiquid(FT_or_toml),
     ice = CM.Parameters.CloudIce(FT_or_toml),
+    Ch2022 = CM.Parameters.Chen2022VelType(FT_or_toml),
 )
 
 microphys_1m_parameters(::Type{FT}) where {FT <: AbstractFloat} =