P3 processes

Project.toml

@@ -22,9 +22,11 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 ProgressLogging = "33c8b6b6-d38a-422a-b730-caa89a2f386c"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TerminalLoggers = "5d786b92-1e48-4d6f-9151-6b4477ca9bed"
+TestEnv = "1e6cf692-eddd-4d53-88a5-2d735e33781b"
 Thermodynamics = "b60c26fb-14c3-4610-9d3e-2d17fe7ff00c"
 
 [compat]
@@ -47,7 +49,6 @@ Plots = "1.29"
 ProgressLogging = "0.1"
 Random = "1"
 SpecialFunctions = "2.3"
-Statistics = "1.10"
 TerminalLoggers = "0.1"
 Thermodynamics = "0.12.6"
 julia = "1.8"

src/Common/initial_condition.jl

@@ -329,6 +329,10 @@ function p3_initial_condition(
     # use approximate temperature values from
     # sounding: kin1d/soundings/Tsfc2.txt
     # (path in p3 fortran github repo)
+    # although I think it could also be grabbed from
+    # this website: https://weather.uwyo.edu/upperair/sounding.html
+    # temperature is a constant gradient with a melting layer at 500 m,
+    # but of course that could easily be changed
     T::FT = -0.004 * (z - 500) + 273.15 # temperature
     p::FT = 990 - (0.1 * z) # pressure
     _q = TD.PhasePartition(q_tot, q_liq, q_ice)

src/Common/ode_utils.jl

@@ -141,6 +141,8 @@ function initialise_aux(
         tmp = similar(ip.q_tot),
         tmp2 = similar(ip.q_tot),
         tmp3 = similar(ip.q_tot),
+        tmp4 = similar(ip.q_tot),
+        tmp5 = similar(ip.q_tot),
         tmp_surface = similar(CC.Fields.level(ip.q_tot, 1), Tuple{FT, FT}),
     )
 
@@ -200,19 +202,31 @@ function initialise_aux(
             q_vap = ip.q_vap,
             ρq_vap = ip.ρq_vap,
         )
+        p3_predicted_eltype = @NamedTuple{
+            F_rim::FT,
+            ρ_rim::FT,
+            F_liq::FT,
+        }
+        p3_predicted = @. p3_predicted_eltype(
+            tuple(
+                copy(ip.ρq_rim / (ip.ρq_ice - ip.ρq_liqonice)),
+                copy(ip.ρq_rim / (ip.ρq_ice - ip.ρq_liqonice) / ip.B_rim),
+                copy(ip.ρq_liqonice / ip.ρq_ice),
+            ),
+        )
         velocities = (;
             term_vel_rai = copy(ip.zero),
             term_vel_ice = copy(ip.zero),
             term_vel_N_rai = copy(ip.zero),
             term_vel_N_ice = copy(ip.zero),
         )
         precip_sources_eltype = @NamedTuple{
-            q_tot::FT,
-            q_liq::FT,
-            q_rai::FT,
-            q_ice::FT,
-            q_rim::FT,
-            q_liqonice::FT,
+            ρq_tot::FT,
+            ρq_liq::FT,
+            ρq_rai::FT,
+            ρq_ice::FT,
+            ρq_rim::FT,
+            ρq_liqonice::FT,
             N_aer::FT,
             N_liq::FT,
             N_rai::FT,
@@ -297,7 +311,7 @@ function initialise_aux(
     if precip isa CloudyPrecip
         aux = merge(aux, (; cloudy_params, cloudy_variables))
     elseif precip isa PrecipitationP3
-        aux = merge(aux, (; p3_boundary_condition))
+        aux = merge(aux, (; p3_boundary_condition, p3_predicted))
     end
 
     return aux