Merge #1849

1849: Fix second order rain flux r=charleskawczynski a=charleskawczynski ### Description This PR: - Fixes the second order rain flux, part of which was missed in #1840 - Moves the first order flux definitions to `tendencies_precipitation.jl`, where they all reside. Co-authored-by: Charles Kawczynski <[email protected]>
CliMA · Dec 16, 2020 · 90bcf13 · 90bcf13
2 parents 11c28e6 + 0c64aa3
commit 90bcf13
Show file tree

Hide file tree

Showing 3 changed files with 57 additions and 56 deletions.
diff --git a/src/Atmos/Model/multiphysics_types.jl b/src/Atmos/Model/multiphysics_types.jl
@@ -50,6 +50,17 @@ RemovePrecipitation(b::Bool) = (
     RemovePrecipitation{TotalMoisture}(b),
 )
 
+"""
+    PrecipitationFlux{PV <: Union{Rain, Snow}} <: TendencyDef{Flux{FirstOrder}, PV}
+
+Computes the precipitation flux as a sum of air velocity and terminal velocity
+multiplied by the advected variable.
+"""
+struct PrecipitationFlux{PV <: Union{Rain, Snow}} <:
+       TendencyDef{Flux{FirstOrder}, PV} end
+
+PrecipitationFlux() = (PrecipitationFlux{Rain}(), PrecipitationFlux{Snow}())
+
 function remove_precipitation_sources(
     s::RemovePrecipitation{PV},
     atmos,

diff --git a/src/Atmos/Model/precipitation.jl b/src/Atmos/Model/precipitation.jl
@@ -56,54 +56,6 @@ function compute_gradient_argument!(
 
 source!(::PrecipitationModel, args...) = nothing
 
-"""
-    PrecipitationFlux{PV <: Union{Rain, Snow}} <: TendencyDef{Flux{FirstOrder}, PV}
-
-Computes the precipitation flux as a sum of air velocity and terminal velocity
-multiplied by the advected variable.
-"""
-struct PrecipitationFlux{PV <: Union{Rain, Snow}} <:
-       TendencyDef{Flux{FirstOrder}, PV} end
-
-PrecipitationFlux() = (PrecipitationFlux{Rain()}, PrecipitationFlux{Snow()})
-
-function flux(::PrecipitationFlux{Rain}, m, state, aux, t, ts, direction)
-    FT = eltype(state)
-    u = state.ρu / state.ρ
-    q_rai = state.precipitation.ρq_rai / state.ρ
-
-    v_term_rai::FT = FT(0)
-    if q_rai > FT(0)
-        v_term_rai = terminal_velocity(
-            m.param_set,
-            m.param_set.microphys.rai,
-            state.ρ,
-            q_rai,
-        )
-    end
-
-    k̂ = vertical_unit_vector(m, aux)
-    return state.precipitation.ρq_rai * (u - k̂ * v_term_rai)
-end
-function flux(::PrecipitationFlux{Snow}, m, state, aux, t, ts, direction)
-    FT = eltype(state)
-    u = state.ρu / state.ρ
-    q_sno = state.precipitation.ρq_sno / state.ρ
-
-    v_term_sno::FT = FT(0)
-    if q_sno > FT(0)
-        v_term_sno = terminal_velocity(
-            m.param_set,
-            m.param_set.microphys.sno,
-            state.ρ,
-            q_sno,
-        )
-    end
-
-    k̂ = vertical_unit_vector(m, aux)
-    return state.precipitation.ρq_sno * (u - k̂ * v_term_sno)
-end
-
 """
     NoPrecipitation <: PrecipitationModel
 
@@ -170,18 +122,17 @@ end
 function flux_second_order!(
     precip::RainModel,
     flux::Grad,
+    atmos::AtmosModel,
     state::Vars,
-    diffusive::Vars,
     aux::Vars,
     t::Real,
-    D_t,
+    ts,
+    diffusive::Vars,
+    hyperdiffusive::Vars,
 )
-    d_q_rai = (-D_t) .* diffusive.precipitation.∇q_rai
-
-    flux_second_order!(precip, flux, state, d_q_rai)
-end
-function flux_second_order!(precip::RainModel, flux::Grad, state::Vars, d_q_rai)
-    flux.precipitation.ρq_rai += d_q_rai * state.ρ
+    tend = Flux{SecondOrder}()
+    args = (atmos, state, aux, t, ts, diffusive, hyperdiffusive)
+    flux.precipitation.ρq_rai = Σfluxes(eq_tends(Rain(), atmos, tend), args...)
 end
 
 function source!(m::RainModel, source::Vars, atmos::AtmosModel, args)

diff --git a/src/Atmos/Model/tendencies_precipitation.jl b/src/Atmos/Model/tendencies_precipitation.jl
@@ -4,6 +4,45 @@
 ##### First order fluxes
 #####
 
+function flux(::PrecipitationFlux{Rain}, m, state, aux, t, ts, direction)
+    FT = eltype(state)
+    u = state.ρu / state.ρ
+    q_rai = state.precipitation.ρq_rai / state.ρ
+
+    v_term_rai::FT = FT(0)
+    if q_rai > FT(0)
+        v_term_rai = terminal_velocity(
+            m.param_set,
+            m.param_set.microphys.rai,
+            state.ρ,
+            q_rai,
+        )
+    end
+
+    k̂ = vertical_unit_vector(m, aux)
+    return state.precipitation.ρq_rai * (u - k̂ * v_term_rai)
+end
+
+function flux(::PrecipitationFlux{Snow}, m, state, aux, t, ts, direction)
+    FT = eltype(state)
+    u = state.ρu / state.ρ
+    q_sno = state.precipitation.ρq_sno / state.ρ
+
+    v_term_sno::FT = FT(0)
+    if q_sno > FT(0)
+        v_term_sno = terminal_velocity(
+            m.param_set,
+            m.param_set.microphys.sno,
+            state.ρ,
+            q_sno,
+        )
+    end
+
+    k̂ = vertical_unit_vector(m, aux)
+    return state.precipitation.ρq_sno * (u - k̂ * v_term_sno)
+end
+
+
 #####
 ##### Second order fluxes
 #####