Merge pull request #4094 from CliMA/sa/fix_sedimentatio_detr

fix sedimentation detr

Author: sajjadazimi

src/prognostic_equations/advection.jl

@@ -396,8 +396,6 @@ function edmfx_sgs_vertical_advection_tendency!(
 
     for j in 1:n
         ᶜa = (@. lazy(draft_area(Y.c.sgsʲs.:($$j).ρa, ᶜρʲs.:($$j))))
-        ᶜright_biased_∂a∂z =
-            @. lazy(ᶜprecipdivᵥ(ᶠinterp(ᶜJ) / ᶠJ * ᶠright_bias(Geometry.WVector(ᶜa))))
 
         # Flux form vertical advection of area farction with the grid mean velocity
         vtt =
@@ -462,7 +460,6 @@ function edmfx_sgs_vertical_advection_tendency!(
                 ᶜqʲ = MatrixFields.get_field(Y, qʲ_name)
                 ᶜqʲₜ = MatrixFields.get_field(Yₜ, qʲ_name)
                 ᶜwʲ = MatrixFields.get_field(p.precomputed, wʲ_name)
-                ᶜaqʲ = (@. lazy(ᶜa * ᶜqʲ))
 
                 # Advective form advection of tracers with updraft velocity
                 va = vertical_advection(
@@ -473,21 +470,16 @@ function edmfx_sgs_vertical_advection_tendency!(
                 @. ᶜqʲₜ += va
 
                 # Flux form sedimentation of tracers
-                vtt = vertical_transport_sedimentation(
-                    ᶜρʲs.:($j),
-                    ᶜwʲ,
-                    ᶜaqʲ,
-                    ᶠJ,
-                )
-                sed_detr = sedimentation_detrainment(
+                vtt = updraft_sedimentation(
                     ᶜρʲs.:($j),
                     ᶜwʲ,
+                    ᶜa,
                     ᶜqʲ,
-                    ᶜright_biased_∂a∂z,
+                    ᶠJ,
                 )
-                @. ᶜqʲₜ += ᶜinv_ρ̂ * (vtt + sed_detr)
-                @. Yₜ.c.sgsʲs.:($$j).q_tot += ᶜinv_ρ̂ * (vtt + sed_detr)
-                @. ᶜ∂ρ∂t_sed += (vtt + sed_detr)
+                @. ᶜqʲₜ += ᶜinv_ρ̂ * vtt
+                @. Yₜ.c.sgsʲs.:($$j).q_tot += ᶜinv_ρ̂ * vtt
+                @. ᶜ∂ρ∂t_sed += vtt
 
                 # Flux form sedimentation of energy
                 if name in (@name(q_liq), @name(q_rai))
@@ -501,19 +493,15 @@ function edmfx_sgs_vertical_advection_tendency!(
                 else
                     error("Unsupported moisture tracer variable")
                 end
-                vtt = vertical_transport_sedimentation(
-                    ᶜρʲs.:($j),
-                    ᶜwʲ,
-                    ᶜaqʲ .* ᶜmse_li,
-                    ᶠJ,
-                )
-                sed_detr = sedimentation_detrainment(
+                vtt = updraft_sedimentation(
                     ᶜρʲs.:($j),
                     ᶜwʲ,
+                    ᶜa,
                     ᶜqʲ .* ᶜmse_li,
-                    ᶜright_biased_∂a∂z,
+                    ᶠJ,
                 )
-                @. Yₜ.c.sgsʲs.:($$j).mse += ᶜinv_ρ̂ * (vtt + sed_detr)
+                @. Yₜ.c.sgsʲs.:($$j).mse +=
+                    ᶜinv_ρ̂ * vtt
             end
 
             # Contribution of density variation due to sedimentation
@@ -552,7 +540,6 @@ function edmfx_sgs_vertical_advection_tendency!(
                 ᶜχʲ = MatrixFields.get_field(Y, χʲ_name)
                 ᶜχʲₜ = MatrixFields.get_field(Yₜ, χʲ_name)
                 ᶜwʲ = MatrixFields.get_field(p.precomputed, wʲ_name)
-                ᶜaχʲ = (@. lazy(ᶜa * ᶜχʲ))
 
                 # Advective form advection of tracers with updraft velocity
                 va = vertical_transport(
@@ -563,23 +550,76 @@ function edmfx_sgs_vertical_advection_tendency!(
                 @. ᶜχʲₜ += va
 
                 # Flux form sedimentation of tracers
-                vtt = vertical_transport_sedimentation(
-                    ᶜρʲs.:($j),
-                    ᶜwʲ,
-                    ᶜaχʲ,
-                    ᶠJ,
-                )
-                sed_detr = sedimentation_detrainment(
+                vtt = updraft_sedimentation(
                     ᶜρʲs.:($j),
                     ᶜwʲ,
+                    ᶜa,
                     ᶜχʲ,
-                    ᶜright_biased_∂a∂z,
+                    ᶠJ,
                 )
-                @. ᶜχʲₜ += ᶜinv_ρ̂ * (vtt + sed_detr)
+                @. ᶜχʲₜ += ᶜinv_ρ̂ * vtt
 
                 # Contribution of density variation due to sedimentation
                 @. ᶜχʲₜ -= ᶜinv_ρ̂ * ᶜχʲ * ᶜ∂ρ∂t_sed
             end
         end
     end
 end
+
+"""
+    updraft_sedimentation(ᶜρ, ᶜw, ᶜa, ᶜχ, ᶠJ)
+
+Compute the sedimentation tendency of tracer `χ` within an updraft, including lateral 
+detrainment when the updraft area increases with height.
+
+# Description
+Sedimenting particles fall with velocity `w` through an updraft of fractional area `a(z)`.  
+The vertical flux divergence gives a tendency of ``∂(ρ w a χ)/∂z``.  
+When `∂a/∂z > 0`, some sedimenting mass exits laterally through the expanding sides, 
+producing a detrainment tendency of ``-ρ w χ ∂a/∂z``.  
+The resulting net tendency in this case is ``a * ∂(ρ w χ)/∂z``.
+
+# Equation
+The lateral flux through the updraft side surface `S` within one grid column is  
+``F_side = ∫_S (ρ χ (w · n)) dS ≈ ρ χ (w · n) A_side,``  
+where `n` is the outward unit normal and `A_side` the side area.  
+For predominantly vertical sedimentation,  
+``w·n A_side ≈ w A_grid [a(z+Δz) - a(z)] = w A_grid Δa.``  
+Dividing by the grid column volume `A_grid·Δz` gives the flux divergence (tendency):  
+``tendency ≈ ρ χ w ∂a/∂z.``  
+A negative sign is applied to represent the loss (detrainment) from the updraft:  
+``Dₛ = -ρ w χ ∂a/∂z.``
+
+# Arguments
+- `ᶜρ`: air density  
+- `ᶜw`: sedimentation velocity (positive downward)  
+- `ᶜa`: updraft area fraction  
+- `ᶜχ`: tracer mixing ratio  
+- `ᶠJ`: face Jacobian (grid geometry)
+
+# Returns
+Tracer tendency due to sedimentation and lateral detrainment.
+"""
+function updraft_sedimentation(
+    ᶜρ,
+    ᶜw,
+    ᶜa,
+    ᶜχ,
+    ᶠJ,
+)
+    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
+    ∂a∂z = @. lazy(ᶜprecipdivᵥ(ᶠinterp(ᶜJ) / ᶠJ * ᶠright_bias(Geometry.WVector(ᶜa))))
+    return @. lazy(
+        ifelse(
+            ∂a∂z < 0,
+            -(ᶜprecipdivᵥ(
+                ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠright_bias(Geometry.WVector(-(ᶜw)) * ᶜa * ᶜχ),
+            )),
+            -(
+                ᶜa * ᶜprecipdivᵥ(
+                    ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠright_bias(Geometry.WVector(-(ᶜw)) * ᶜχ),
+                )
+            ),
+        ),
+    )
+end

src/prognostic_equations/edmfx_entr_detr.jl

@@ -526,40 +526,6 @@ function detrainment(
     return max(detr, 0)
 end
 
-"""
-    sedimentation_detrainment(ᶜρ, ᶜw, ᶜχ, ᶜ∂a∂z)
-
-Compute the lateral detrainment tendency of a sedimenting tracer `χ` within a tilted updraft.
-
-# Description
-Sedimenting particles moving downward through an updraft with a height-varying area fraction `a(z)` 
-can exit laterally through its sloping sides. When the updraft area increases with height (`∂a/∂z > 0`), 
-this produces a lateral outflow (detrainment) of tracer mass proportional to the local air density `ρ`, 
-sedimentation velocity `w`, tracer mixing ratio `χ`, and the vertical gradient of area fraction `∂a/∂z`.
-
-# Equation
-The lateral flux through the updraft side surface `S` within one grid column is  
-F_side = ∫_S (ρ χ (w · n)) dS ≈ ρ χ (w · n) A_side,  
-where `n` is the outward unit normal and `A_side` the side area.  
-For predominantly vertical sedimentation,  
-w·n A_side ≈ w A_grid [a(z+Δz) - a(z)] = w A_grid Δa.  
-Dividing by the grid column volume `A_grid·Δz` gives the flux divergence (tendency):  
-tendency ≈ ρ χ w ∂a/∂z. A negative sign is applied to represent the loss (detrainment) from the updraft:  
-Dₛ = -ρ w χ ∂a/∂z.
-
-# Arguments
-- `ᶜρ`: air density at cell center  
-- `ᶜw`: sedimentation velocity (positive downward)  
-- `ᶜχ`: tracer mixing ratio  
-- `ᶜ∂a∂z`: vertical derivative of updraft area fraction  
-
-# Returns
-Lateral detrainment tendency of tracer `χ` due to sedimentation.
-"""
-function sedimentation_detrainment(ᶜρ, ᶜw, ᶜχ, ᶜ∂a∂z)
-    return @. lazy(-1 * ᶜρ * ᶜw * ᶜχ * max(0, ᶜ∂a∂z))
-end
-
 function turbulent_entrainment(turbconv_params, ᶜaʲ)
     turb_entr_param_vec = CAP.turb_entr_param_vec(turbconv_params)
     return max(turb_entr_param_vec[1] * exp(-turb_entr_param_vec[2] * ᶜaʲ), 0)

src/prognostic_equations/implicit/implicit_tendency.jl

@@ -97,18 +97,6 @@ function vertical_transport(ᶜρ, ᶠu³, ᶜχ, dt, ::Val{:third_order})
     return @. lazy(-(ᶜadvdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠupwind3(ᶠu³, ᶜχ))))
 end
 
-function vertical_transport_sedimentation(
-    ᶜρ,
-    ᶜw,
-    ᶜχ,
-    ᶠJ,
-)
-    ᶜJ = Fields.local_geometry_field(axes(ᶜρ)).J
-    return @. lazy(
-        -(ᶜprecipdivᵥ(ᶠinterp(ᶜρ * ᶜJ) / ᶠJ * ᶠright_bias(Geometry.WVector(-(ᶜw)) * ᶜχ))),
-    )
-end
-
 vertical_advection(ᶠu³, ᶜχ, ::Val{:none}) =
     @. lazy(-(ᶜadvdivᵥ(ᶠu³ * ᶠinterp(ᶜχ)) - ᶜχ * ᶜadvdivᵥ(ᶠu³)))
 vertical_advection(ᶠu³, ᶜχ, ::Val{:first_order}) =

src/prognostic_equations/implicit/manual_sparse_jacobian.jl

@@ -871,15 +871,18 @@ function update_jacobian!(alg::ManualSparseJacobian, cache, Y, p, dtγ, t)
                             DiagonalMatrixRow(ᶜχʲ) ⋅ ᶜadvdivᵥ_matrix()
                         ) ⋅ DiagonalMatrixRow(g³³(ᶠgⁱʲ))
 
-                    # precipitation
+                    # sedimentation
                     @. ᶜtridiagonal_matrix_scalar =
-                        dtγ * -(ᶜprecipdivᵥ_matrix()) ⋅
-                        DiagonalMatrixRow(ᶠinterp(ᶜρʲs.:(1) * ᶜJ) / ᶠJ) ⋅
-                        ᶠright_bias_matrix() ⋅
-                        DiagonalMatrixRow(-Geometry.WVector(ᶜwʲ) * ᶜa)
-                    # precipitation detrainment
-                    @. ᶜtridiagonal_matrix_scalar +=
-                        dtγ * -DiagonalMatrixRow(ᶜρʲs.:(1) * ᶜwʲ * max(0, ᶜ∂a∂z))
+                        dtγ * ifelse(ᶜ∂a∂z < 0,
+                            -(ᶜprecipdivᵥ_matrix()) ⋅
+                            DiagonalMatrixRow(ᶠinterp(ᶜρʲs.:(1) * ᶜJ) / ᶠJ) ⋅
+                            ᶠright_bias_matrix() ⋅
+                            DiagonalMatrixRow(-Geometry.WVector(ᶜwʲ) * ᶜa),
+                            -DiagonalMatrixRow(ᶜa) ⋅ ᶜprecipdivᵥ_matrix() ⋅
+                            DiagonalMatrixRow(ᶠinterp(ᶜρʲs.:(1) * ᶜJ) / ᶠJ) ⋅
+                            ᶠright_bias_matrix() ⋅
+                            DiagonalMatrixRow(-Geometry.WVector(ᶜwʲ)),
+                        )
 
                     @. ∂ᶜχʲ_err_∂ᶜχʲ +=
                         DiagonalMatrixRow(ᶜinv_ρ̂) ⋅ ᶜtridiagonal_matrix_scalar