boundary?

CliMA · Aug 12, 2024 · ff48555 · ff48555
1 parent 02bb16b
commit ff48555
Show file tree

Hide file tree

Showing 3 changed files with 70 additions and 9 deletions.
diff --git a/src/K1DModel/tendency.jl b/src/K1DModel/tendency.jl
@@ -240,7 +240,38 @@ end
 
     return dY
 end
-@inline function advection_tendency!(::Union{CO.NonEquilibriumMoisture, CO.MoistureP3}, dY, Y, aux, t)
+@inline function advection_tendency!(::CO.NonEquilibriumMoisture, dY, Y, aux, t)
+    FT = eltype(Y.ρq_tot)
+
+    If = CC.Operators.InterpolateC2F()
+
+    ∂_qt = CC.Operators.DivergenceF2C(
+        bottom = CC.Operators.SetValue(CC.Geometry.WVector(aux.prescribed_velocity.ρw0 * aux.q_surf)),
+        top = CC.Operators.Extrapolate(),
+    )
+    ∂_ql = CC.Operators.DivergenceF2C(
+        bottom = CC.Operators.SetValue(CC.Geometry.WVector(aux.prescribed_velocity.ρw0 * FT(0.0))),
+        top = CC.Operators.Extrapolate(),
+    )
+    ∂_qi = CC.Operators.DivergenceF2C(
+        bottom = CC.Operators.SetValue(CC.Geometry.WVector(aux.prescribed_velocity.ρw0 * FT(0.0))),
+        top = CC.Operators.Extrapolate(),
+    )
+
+    @. dY.ρq_tot += -∂_qt(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) * If(Y.ρq_tot))
+    @. dY.ρq_liq += -∂_ql(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) * If(Y.ρq_liq))
+    @. dY.ρq_ice += -∂_qi(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ) * If(Y.ρq_ice))
+
+    fcc = CC.Operators.FluxCorrectionC2C(bottom = CC.Operators.Extrapolate(), top = CC.Operators.Extrapolate())
+    if Bool(aux.kid_params.qtot_flux_correction)
+        @. dY.ρq_tot += fcc(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ), Y.ρq_tot)
+    end
+    @. dY.ρq_liq += fcc(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ), Y.ρq_liq)
+    @. dY.ρq_ice += fcc(aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ), Y.ρq_ice)
+
+    return dY
+end
+@inline function advection_tendency!(::CO.MoistureP3, dY, Y, aux, t)
     FT = eltype(Y.ρq_tot)
 
     If = CC.Operators.InterpolateC2F()
@@ -401,21 +432,47 @@ end
 
     # TODO add compat with 2M rain scheme
 
+    _q_flux = FT(0.65e-4)
+    _N_flux = FT(25000)
+    _F_rim = FT(0)
+    _F_liq = FT(0)
+    _ρ_r_init = FT(900)
+
     If = CC.Operators.InterpolateC2F()
-    ∂ = CC.Operators.DivergenceF2C(
-        bottom = CC.Operators.Extrapolate(),
-        top = CC.Operators.SetValue(CC.Geometry.WVector(0.0)),
+    ∂q_ice = CC.Operators.DivergenceF2C(
+        bottom = CC.Operators.SetValue(CC.Geometry.WVector(aux.prescribed_velocity.ρw0 * aux.q_surf)),
+        top = CC.Operators.SetValue(CC.Geometry.WVector(FT(-0.05) * _q_flux)),
+    )
+
+    ∂q_rim = CC.Operators.DivergenceF2C(
+        bottom = CC.Operators.SetValue(CC.Geometry.WVector(aux.prescribed_velocity.ρw0 * aux.q_surf)),
+        top = CC.Operators.SetValue(CC.Geometry.WVector(FT(-0.05) * _F_rim * _q_flux)),
+    )
+
+    ∂q_liqonice = CC.Operators.DivergenceF2C(
+        bottom = CC.Operators.SetValue(CC.Geometry.WVector(aux.prescribed_velocity.ρw0 * aux.q_surf)),
+        top = CC.Operators.SetValue(CC.Geometry.WVector(FT(-0.05) * _F_liq * _q_flux)),
+    )
+
+    ∂B_rim = CC.Operators.DivergenceF2C(
+        bottom = CC.Operators.SetValue(CC.Geometry.WVector(aux.prescribed_velocity.ρw0 * aux.q_surf)),
+        top = CC.Operators.SetValue(CC.Geometry.WVector(FT(-0.05) * _F_rim * _q_flux / _ρ_r_init)),
+    )
+
+    ∂N_ice = CC.Operators.DivergenceF2C(
+        bottom = CC.Operators.SetValue(CC.Geometry.WVector(aux.prescribed_velocity.ρw0 * aux.q_surf)),
+        top = CC.Operators.SetValue(CC.Geometry.WVector(FT(-0.05) * _N_flux)),
     )
 
-    @. dY.ρq_ice += ∂(
+    @. dY.ρq_ice += ∂q_ice(
         FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) +
         (CC.Geometry.WVector(If(aux.velocities.term_vel_ice))) * If(Y.ρq_ice),
     )
-    @. dY.ρq_rim += ∂(
+    @. dY.ρq_rim += ∂q_rim(
         FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) +
         (CC.Geometry.WVector(If(aux.velocities.term_vel_ice))) * If(Y.ρq_rim),
     )
-    @. dY.B_rim += ∂(
+    @. dY.B_rim += ∂B_rim(
         FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) +
         (CC.Geometry.WVector(If(aux.velocities.term_vel_ice))) * If(Y.B_rim),
     )
@@ -427,7 +484,7 @@ end
     #             CC.Geometry.WVector(If(aux.velocities.term_vel_ice))
     #         ) * If(Y.ρq_liqonice),
     #     )
-    @. dY.N_ice += ∂(
+    @. dY.N_ice += ∂N_ice(
         FT(-1) * (aux.prescribed_velocity.ρw / If(aux.thermo_variables.ρ)) +
         (CC.Geometry.WVector(If(aux.velocities.term_vel_N_ice))) * If(Y.N_ice),
     )

diff --git a/test/unit_tests/common_unit_test.jl b/test/unit_tests/common_unit_test.jl
@@ -132,7 +132,7 @@ end
     @test LA.norm(state.N_rai) == 0
     @test LA.norm(state.N_aer) == 0
 
-    state = CO.initialise_state(moisture_p3, precip_p3, initial_profiles)
+    state = CO.initialise_state(p3_moist, precip_p3, initial_profiles)
     @test state isa CC.Fields.FieldVector
     @test LA.norm(state.ρq_tot) == 0
     @test LA.norm(state.ρq_liq) == 0

diff --git a/test/unit_tests/unit_test.jl b/test/unit_tests/unit_test.jl
@@ -30,9 +30,12 @@ thermo_params = create_thermodynamics_parameters(toml_dict)
 kid_params = create_kid_parameters(toml_dict)
 air_params = CMP.AirProperties(toml_dict)
 activation_params = CMP.AerosolActivationParameters(toml_dict)
+p3 = CMP.ParametersP3(FT)
+Chen2022 = CMP.Chen2022VelType(FT)
 # ... for cloud condensate options ...
 equil_moist = CO.EquilibriumMoisture()
 nequil_moist = CO.NonEquilibriumMoisture(CMP.CloudLiquid(toml_dict), CMP.CloudIce(toml_dict))
+p3_moist = CO.MoistureP3()
 # ... and precipitation options
 no_precip = CO.NoPrecipitation()
 precip_0m = CO.Precipitation0M(CMP.Parameters0M(toml_dict))
@@ -46,6 +49,7 @@ precip_1m = CO.Precipitation1M(
     CMP.Blk1MVelType(toml_dict),
 )
 precip_2m = CO.Precipitation2M(CMP.SB2006(toml_dict), CMP.SB2006VelType(toml_dict))
+precip_p3 = CO.PrecipitationP3(p3, Chen2022)
 
 # common unit tests
 include("./common_unit_test.jl")