changing aerosol sources; using an approximate rate instead of updati…

…ng the state Y
CliMA · Apr 21, 2023 · 5bdee22 · 5bdee22
1 parent bdcabe4
commit 5bdee22
Show file tree

Hide file tree

Showing 3 changed files with 10 additions and 8 deletions.
diff --git a/src/driveKiD/helper_functions.jl b/src/driveKiD/helper_functions.jl
@@ -8,7 +8,7 @@ end
 """
     Returns the number of new activated aerosol particles and updates aerosol number density
 """
-@inline function aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw)
+@inline function aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw, dt)
 
     microphys_params = KP.microphysics_params(params)
     thermo_params = CM.Parameters.thermodynamics_params(microphys_params)
@@ -36,7 +36,7 @@ end
         return (; S_Nl, S_Na)
     end
 
-    S_Nl = max(0, N_act - (N_aer_0 - N_aer))
+    S_Nl = max(0, N_act - (N_aer_0 - N_aer)) / dt
     S_Na = -S_Nl
 
     return (; S_Nl, S_Na)

diff --git a/src/driveKiD/tendency.jl b/src/driveKiD/tendency.jl
@@ -557,9 +557,10 @@ end
         aux.moisture_variables.p,
         aux.moisture_variables.ρ,
         ClimaCore.Operators.InterpolateF2C().(aux.prescribed_velocity.ρw.components.data.:1),
+        aux.TS.dt,
     )
-    @. Y.N_liq += tmp.S_Nl
-    @. Y.N_aer += tmp.S_Na
+    aux.aerosol_variables.S_N_aer = tmp.S_Na
+    aux.precip_sources.S_N_liq = tmp.S_Nl
 
     aux.precip_variables.q_rai = Y.ρq_rai ./ aux.moisture_variables.ρ
     aux.precip_variables.N_liq = Y.N_liq
@@ -579,7 +580,7 @@ end
     aux.precip_sources.S_q_rai = tmp.S_q_rai
     aux.precip_sources.S_q_tot = tmp.S_q_tot
     aux.precip_sources.S_q_liq = tmp.S_q_liq
-    aux.precip_sources.S_N_liq = tmp.S_N_liq
+    @. aux.precip_sources.S_N_liq += tmp.S_N_liq
     aux.precip_sources.S_N_rai = tmp.S_N_rai
 
     If = CC.Operators.InterpolateC2F(bottom = CC.Operators.Extrapolate(), top = CC.Operators.Extrapolate())

diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl
@@ -582,24 +582,25 @@ end
     p = 1e5
     ρ = 1.0
     ρw = 2.0
+    dt = 1.0
 
     #action
-    tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw)
+    tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, ρw, dt)
 
     #test
     @test tmp isa NamedTuple
     @test tmp.S_Nl ≈ -tmp.S_Na
     @test 0 < tmp.S_Nl < N_aer_0
 
     #action
-    tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, 290.0, p, ρ, ρw)
+    tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, 290.0, p, ρ, ρw, dt)
 
     #test
     @test tmp.S_Nl ≈ 0.0 atol = eps(Float64)
     @test tmp.S_Na ≈ 0.0 atol = eps(Float64)
 
     #action
-    tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, 0.0)
+    tmp = KID.aerosol_activation_helper(params, q_tot, q_liq, N_aer, N_aer_0, T, p, ρ, 0.0, dt)
 
     #test
     @test tmp.S_Nl ≈ 0.0 atol = eps(Float64)