ft: Add Kinematic Driver (KiD) model configuration

temp: commented out various tendencies
todo: debug why we have `clw` at surface

Author: haakon-e

config/default_configs/default_config.yml

@@ -419,3 +419,6 @@ use_itime:
   1. Surface conditions that explicitly depend on time (e.g. LifeCycleTan2018, TRMM_LBA, etc.),
   2. Time dependent forcing/tendencies use time rounded to the nearest unit of time for dt"
   value: false
+prescribed_flow:
+  help: "Prescribe a flow field [`nothing` (default), `true`]"
+  value: ~

config/model_configs/kinematic_driver.yml

@@ -0,0 +1,22 @@
+## Model
+prescribed_flow: true
+initial_condition: "ShipwayHill2012"
+tracer_upwinding: first_order
+moist: "nonequil"
+cloud_model: "grid_scale"
+call_cloud_diagnostics_per_stage: true
+precip_model: "1M"
+config: "column"
+## Simulation
+z_max: 4e3
+z_elem: 100
+z_stretch: false
+dt: "1secs"
+t_end: "1hours"
+dt_save_state_to_disk: "1hours"
+toml: [toml/kinematic_driver.toml]
+## Diagnostics
+output_default_diagnostics: false
+diagnostics:
+  - short_name: [ta, thetaa, ha, rhoa, wa, hur, hus, cl, clw, cli, husra, hussn]
+    period: 10secs

examples/hybrid/KiD_driver.jl

@@ -0,0 +1,70 @@
+import ClimaComms
+import ClimaAtmos as CA
+import ClimaCore: Fields
+import YAML
+import ClimaComms
+
+import Random
+# Random.seed!(Random.MersenneTwister())
+Random.seed!(1234)
+
+# --> get config
+configs_path = joinpath(pkgdir(CA), "config/model_configs/")
+pth = joinpath(configs_path, "kinematic_driver.yml");
+job_id = "kinematic_driver";
+config_dict = YAML.load_file(pth)
+# <--
+
+
+config = CA.AtmosConfig(config_dict; job_id)
+simulation = CA.get_simulation(config);
+
+sol_res = CA.solve_atmos!(simulation);  # solve!
+
+(; integrator) = simulation;
+(; p) = integrator;
+(; atmos, params) = p;
+
+# --> Make ci plots
+# ]add ClimaAnalysis, ClimaCoreSpectra
+include(joinpath(pkgdir(CA), "post_processing", "ci_plots.jl"))
+ref_job_id = config.parsed_args["reference_job_id"]
+reference_job_id = isnothing(ref_job_id) ? simulation.job_id : ref_job_id
+make_plots(Val(Symbol(reference_job_id)), simulation.output_dir)
+# <--
+
+# --> ClimaAnalysis
+import ClimaAnalysis
+# using ClimaAnalysis.Visualize
+import ClimaAnalysis.Visualize as viz
+using ClimaAnalysis.Utils: kwargs
+using CairoMakie;
+CairoMakie.activate!();
+# using GLMakie; GLMakie.activate!()
+simdir = ClimaAnalysis.SimDir(simulation.output_dir);
+
+# entr = get(simdir; short_name = "entr")
+# entr.dims  # (time, x, y, z)
+
+# fig = Figure();
+# viz.plot!(fig, entr, time=0, x=0, y=0, more_kwargs = Dict(:axis => kwargs(dim_on_y = true)))
+# viz.plot!(fig, entr, x=0, y=0);
+# fig
+# <--
+
+
+
+#=
+%use upwinding for the rain -- yes!
+
+qr, qs, all specific humidities,
+
+take 30% acoustic Courant number c=300m/s, divided by vertical res
+- for ~200
+
+ill make some plots
+
+is smag applied to qr??? check this!
+
+most likely precip would cause blow-ups.
+=#

src/initial_conditions/initial_conditions.jl

@@ -1618,3 +1618,49 @@ function (initial_condition::ISDAC)(params)
         )
     end
 end
+
+"""
+    ShipwayHill2012
+
+The `InitialCondition` described in [ShipwayHill2012](@cite), but with a hydrostatically
+balanced pressure profile.
+
+B. J. Shipway and A. A. Hill. 
+Diagnosis of systematic differences between multiple parametrizations of warm rain microphysics using a kinematic framework. 
+Quarterly Journal of the Royal Meteorological Society 138, 2196-2211 (2012).
+"""
+Base.@kwdef struct ShipwayHill2012 <: InitialCondition end
+
+function (initial_condition::ShipwayHill2012)(params)
+    FT = eltype(params)
+
+    ## Initialize the profile
+    z_0 = FT(0.0)   # 0.0
+    z_1 = FT(740.0)   # 740.0
+    z_2 = FT(3260.0)   # 3260.0
+    rv_0 = FT(0.015) # 0.015
+    rv_1 = FT(0.0138) # 0.0138
+    rv_2 = FT(0.0024) # 0.0024
+    θ_0 = FT(297.9)   # 297.9
+    θ_1 = FT(297.9)   # 297.9
+    θ_2 = FT(312.66)   # 312.66
+
+    # profile of water vapour mixing ratio
+    rv(z) = z < z_1 ? rv_0 + (rv_1 - rv_0) / (z_1 - z_0) * (z - z_0) : rv_1 + (rv_2 - rv_1) / (z_2 - z_1) * (z - z_1)
+    q_tot(z) = rv(z) / (1 + rv(z))
+
+    # profile of potential temperature
+    θ(z) = z < z_1 ? θ_0 : θ_1 + (θ_2 - θ_1) / (z_2 - z_1) * (z - z_1)
+    ## 
+    thermo_params = CAP.thermodynamics_params(params)
+    p_0 = FT(100700) # Pa
+    p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot)
+    function local_state(local_geometry)
+        (; z) = local_geometry.coordinates
+        return LocalState(; params, geometry = local_geometry,
+            thermo_state = TD.PhaseEquil_pθq(thermo_params, p(z), θ(z), q_tot(z)),
+            precip_state = PrecipStateMassNum(; q_rai = FT(0), q_sno = FT(0)),
+        )
+    end
+    return local_state
+end

src/prognostic_equations/advection.jl

@@ -279,7 +279,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         )
         vtt = vertical_transport(ᶜρ, ᶠu³, ᶜh_tot, FT(dt), energy_q_tot_upwinding)
         vtt_central = vertical_transport(ᶜρ, ᶠu³, ᶜh_tot, FT(dt), Val(:none))
-        @. Yₜ.c.ρe_tot += vtt - vtt_central
+        # @. Yₜ.c.ρe_tot += vtt - vtt_central
     end
 
     if !(p.atmos.moisture_model isa DryModel) && energy_q_tot_upwinding != Val(:none)
@@ -289,35 +289,35 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
         @. Yₜ.c.ρq_tot += vtt - vtt_central
     end
 
-    if isnothing(ᶠf¹²)
-        # shallow atmosphere
-        @. Yₜ.c.uₕ -=
-            ᶜinterp(ᶠω¹² × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) / (Y.c.ρ * ᶜJ) +
-            (ᶜf³ + ᶜω³) × CT12(ᶜu)
-        @. Yₜ.f.u₃ -= ᶠω¹² × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
-        for j in 1:n
-            @. Yₜ.f.sgsʲs.:($$j).u₃ -=
-                ᶠω¹²ʲs.:($$j) × ᶠinterp(CT12(ᶜuʲs.:($$j))) +
-                ᶠgradᵥ(ᶜKʲs.:($$j) - ᶜinterp(ᶠKᵥʲs.:($$j)))
-        end
-    else
-        # deep atmosphere
-        @. Yₜ.c.uₕ -=
-            ᶜinterp((ᶠf¹² + ᶠω¹²) × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) /
-            (Y.c.ρ * ᶜJ) + (ᶜf³ + ᶜω³) × CT12(ᶜu)
-        @. Yₜ.f.u₃ -= (ᶠf¹² + ᶠω¹²) × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
-        for j in 1:n
-            @. Yₜ.f.sgsʲs.:($$j).u₃ -=
-                (ᶠf¹² + ᶠω¹²ʲs.:($$j)) × ᶠinterp(CT12(ᶜuʲs.:($$j))) +
-                ᶠgradᵥ(ᶜKʲs.:($$j) - ᶜinterp(ᶠKᵥʲs.:($$j)))
-        end
-    end
-
-    if use_prognostic_tke(turbconv_model) # advect_tke triggers allocations
-        @. ᶜa_scalar = ᶜtke⁰ * draft_area(ᶜρa⁰, ᶜρ⁰)
-        vtt = vertical_transport(ᶜρ⁰, ᶠu³⁰, ᶜa_scalar, dt, edmfx_mse_q_tot_upwinding)
-        @. Yₜ.c.sgs⁰.ρatke += vtt
-    end
+    # if isnothing(ᶠf¹²)
+    #     # shallow atmosphere
+    #     @. Yₜ.c.uₕ -=
+    #         ᶜinterp(ᶠω¹² × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) / (Y.c.ρ * ᶜJ) +
+    #         (ᶜf³ + ᶜω³) × CT12(ᶜu)
+    #     @. Yₜ.f.u₃ -= ᶠω¹² × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
+    #     for j in 1:n
+    #         @. Yₜ.f.sgsʲs.:($$j).u₃ -=
+    #             ᶠω¹²ʲs.:($$j) × ᶠinterp(CT12(ᶜuʲs.:($$j))) +
+    #             ᶠgradᵥ(ᶜKʲs.:($$j) - ᶜinterp(ᶠKᵥʲs.:($$j)))
+    #     end
+    # else
+    #     # deep atmosphere
+    #     @. Yₜ.c.uₕ -=
+    #         ᶜinterp((ᶠf¹² + ᶠω¹²) × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) /
+    #         (Y.c.ρ * ᶜJ) + (ᶜf³ + ᶜω³) × CT12(ᶜu)
+    #     @. Yₜ.f.u₃ -= (ᶠf¹² + ᶠω¹²) × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
+    #     for j in 1:n
+    #         @. Yₜ.f.sgsʲs.:($$j).u₃ -=
+    #             (ᶠf¹² + ᶠω¹²ʲs.:($$j)) × ᶠinterp(CT12(ᶜuʲs.:($$j))) +
+    #             ᶠgradᵥ(ᶜKʲs.:($$j) - ᶜinterp(ᶠKᵥʲs.:($$j)))
+    #     end
+    # end
+
+    # if use_prognostic_tke(turbconv_model) # advect_tke triggers allocations
+    #     @. ᶜa_scalar = ᶜtke⁰ * draft_area(ᶜρa⁰, ᶜρ⁰)
+    #     vtt = vertical_transport(ᶜρ⁰, ᶠu³⁰, ᶜa_scalar, dt, edmfx_mse_q_tot_upwinding)
+    #     @. Yₜ.c.sgs⁰.ρatke += vtt
+    # end
 end
 
 """

src/prognostic_equations/dss.jl

@@ -77,8 +77,23 @@ dss!(Y_state, params, t_current)
 """
 
 NVTX.@annotate function dss!(Y, p, t)
-    if do_dss(axes(Y.c))
-        Spaces.weighted_dss!(Y.c => p.ghost_buffer.c, Y.f => p.ghost_buffer.f)
-    end
+    # if do_dss(axes(Y.c))
+    #     Spaces.weighted_dss!(Y.c => p.ghost_buffer.c, Y.f => p.ghost_buffer.f)
+    # end
+    prescribe_flow!(Y, p, t, p.atmos.prescribed_flow)
     return nothing
 end
+
+prescribe_flow!(Y, p, t, ::Nothing) = nothing
+function prescribe_flow!(Y, p, t, flow::PrescribedFlow)
+    FT = eltype(p.params)
+    # (; uₕ, u₃, ρ, p) = flow
+    (; prescribed_u₃) = flow
+    lg = Fields.local_geometry_field(Y.f)
+    # Main.@infiltrate
+    @. Y.f.u₃ = C3(Geometry.WVector(prescribed_u₃(FT, t)), lg)
+    # @. Y.c.uₕ = uₕ(t)
+    # @. Y.c.ρ = ρ₀  # move to IC
+    # @. Y.c.p = p₀  # move to IC
+    return nothing
+end