Skip to content

Commit cd06566

Browse files
haakon-ehaakon-e
committed
ft: Add Kinematic Driver (KiD) model configuration
temp: commented out various tendencies todo: debug why we have `clw` at surface
1 parent 79028dd commit cd06566

File tree

12 files changed

+251
-37
lines changed

12 files changed

+251
-37
lines changed

‎config/default_configs/default_config.yml‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -419,3 +419,6 @@ use_itime:
419419
1. Surface conditions that explicitly depend on time (e.g. LifeCycleTan2018, TRMM_LBA, etc.),
420420
2. Time dependent forcing/tendencies use time rounded to the nearest unit of time for dt"
421421
value: false
422+
prescribed_flow:
423+
help: "Prescribe a flow field [`nothing` (default), `true`]"
424+
value: ~

‎config/model_configs/kinematic_driver.yml‎

Lines changed: 25 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
## Model
2+
prescribed_flow: true
3+
initial_condition: "ShipwayHill2012"
4+
surface_setup: "ShipwayHill2012"
5+
tracer_upwinding: first_order
6+
moist: "nonequil"
7+
cloud_model: "grid_scale"
8+
call_cloud_diagnostics_per_stage: true
9+
precip_model: "1M"
10+
config: "column"
11+
## Simulation
12+
z_max: 4e3
13+
z_elem: 100
14+
z_stretch: false
15+
dt: "1secs"
16+
t_end: "1hours"
17+
# t_end: "10secs"
18+
dt_save_state_to_disk: "1hours"
19+
toml: [toml/kinematic_driver.toml]
20+
check_nan_every: 1
21+
## Diagnostics
22+
output_default_diagnostics: false
23+
diagnostics:
24+
- short_name: [ta, thetaa, ha, rhoa, wa, hur, hus, cl, clw, cli, husra, hussn]
25+
period: 10secs

‎docs/bibliography.bib‎

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -283,3 +283,16 @@ @article{Wing2018
283283
url = {https://gmd.copernicus.org/articles/11/793/2018/},
284284
doi = {10.5194/gmd-11-793-2018}
285285
}
286+
287+
288+
@article{ShipwayHill2012,
289+
title = {Diagnosis of systematic differences between multiple parametrizations of warm rain microphysics using a kinematic framework},
290+
volume = {138},
291+
url = {https://onlinelibrary.wiley.com/doi/abs/10.1002/qj.1913},
292+
doi = {10.1002/qj.1913},
293+
pages = {2196--2211},
294+
number = {669},
295+
journaltitle = {Quarterly Journal of the Royal Meteorological Society},
296+
author = {Shipway, B. J. and Hill, A. A.},
297+
date = {2012},
298+
}

‎examples/hybrid/KiD_driver.jl‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
import ClimaComms
2+
import ClimaAtmos as CA
3+
import ClimaCore: Fields
4+
import YAML
5+
import ClimaComms
6+
7+
import Random
8+
# Random.seed!(Random.MersenneTwister())
9+
Random.seed!(1234)
10+
11+
# --> get config
12+
configs_path = joinpath(pkgdir(CA), "config/model_configs/")
13+
pth = joinpath(configs_path, "kinematic_driver.yml");
14+
job_id = "kinematic_driver";
15+
config_dict = YAML.load_file(pth)
16+
# <--
17+
18+
19+
config = CA.AtmosConfig(config_dict; job_id)
20+
simulation = CA.get_simulation(config);
21+
22+
sol_res = CA.solve_atmos!(simulation); # solve!
23+
24+
(; integrator) = simulation;
25+
(; p) = integrator;
26+
(; atmos, params) = p;
27+
28+
# --> Make ci plots
29+
# ]add ClimaAnalysis, ClimaCoreSpectra
30+
include(joinpath(pkgdir(CA), "post_processing", "ci_plots.jl"))
31+
ref_job_id = config.parsed_args["reference_job_id"]
32+
reference_job_id = isnothing(ref_job_id) ? simulation.job_id : ref_job_id
33+
make_plots(Val(Symbol(reference_job_id)), simulation.output_dir)
34+
# <--
35+
36+
# --> ClimaAnalysis
37+
import ClimaAnalysis
38+
# using ClimaAnalysis.Visualize
39+
import ClimaAnalysis.Visualize as viz
40+
using ClimaAnalysis.Utils: kwargs
41+
using CairoMakie;
42+
CairoMakie.activate!();
43+
# using GLMakie; GLMakie.activate!()
44+
simdir = ClimaAnalysis.SimDir(simulation.output_dir);
45+
46+
# entr = get(simdir; short_name = "entr")
47+
# entr.dims # (time, x, y, z)
48+
49+
# fig = Figure();
50+
# viz.plot!(fig, entr, time=0, x=0, y=0, more_kwargs = Dict(:axis => kwargs(dim_on_y = true)))
51+
# viz.plot!(fig, entr, x=0, y=0);
52+
# fig
53+
# <--
54+
55+
56+
57+
#=
58+
%use upwinding for the rain -- yes!
59+
60+
qr, qs, all specific humidities,
61+
62+
take 30% acoustic Courant number c=300m/s, divided by vertical res
63+
- for ~200
64+
65+
ill make some plots
66+
67+
is smag applied to qr??? check this!
68+
69+
most likely precip would cause blow-ups.
70+
=#

‎src/initial_conditions/initial_conditions.jl‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1618,3 +1618,47 @@ function (initial_condition::ISDAC)(params)
16181618
)
16191619
end
16201620
end
1621+
1622+
"""
1623+
ShipwayHill2012
1624+
1625+
The `InitialCondition` described in [ShipwayHill2012](@cite), but with a hydrostatically
1626+
balanced pressure profile.
1627+
1628+
B. J. Shipway and A. A. Hill.
1629+
Diagnosis of systematic differences between multiple parametrizations of warm rain microphysics using a kinematic framework.
1630+
Quarterly Journal of the Royal Meteorological Society 138, 2196-2211 (2012).
1631+
"""
1632+
Base.@kwdef struct ShipwayHill2012 <: InitialCondition end
1633+
1634+
function (initial_condition::ShipwayHill2012)(params)
1635+
FT = eltype(params)
1636+
1637+
## Initialize the profile
1638+
z_0, z_1, z_2 = FT(0), FT(740), FT(3260) # m
1639+
rv_0, rv_1, rv_2 = FT(0.015), FT(0.0138), FT(0.0024) # kg/kg
1640+
θ_0, θ_1, θ_2 = FT(297.9), FT(297.9), FT(312.66) # K
1641+
1642+
# profile of water vapour mixing ratio
1643+
linear_profile(z₀, z₁, x₀, x₁, z) = x₀ + (x₁ - x₀) / (z₁ - z₀) * (z - z₀)
1644+
rv(z) = if z < z_1
1645+
linear_profile(z_0, z_1, rv_0, rv_1, z)
1646+
else
1647+
linear_profile(z_1, z_2, rv_1, rv_2, z)
1648+
end
1649+
q_tot(z) = rv(z) / (1 + rv(z))
1650+
# profile of potential temperature
1651+
θ(z) = z < z_1 ? θ_0 : linear_profile(z_1, z_2, θ_1, θ_2, z)
1652+
##
1653+
thermo_params = CAP.thermodynamics_params(params)
1654+
p_0 = FT(100700) # Pa
1655+
p = hydrostatic_pressure_profile(; thermo_params, p_0, θ, q_tot)
1656+
function local_state(local_geometry)
1657+
(; z) = local_geometry.coordinates
1658+
return LocalState(; params, geometry = local_geometry,
1659+
thermo_state = TD.PhaseEquil_pθq(thermo_params, p(z), θ(z), q_tot(z)),
1660+
precip_state = PrecipStateMassNum(; q_rai = FT(0), q_sno = FT(0)),
1661+
)
1662+
end
1663+
return local_state
1664+
end

‎src/prognostic_equations/advection.jl‎

Lines changed: 30 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -279,7 +279,7 @@ NVTX.@annotate function explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
279279
)
280280
vtt = vertical_transport(ᶜρ, ᶠu³, ᶜh_tot, FT(dt), energy_q_tot_upwinding)
281281
vtt_central = vertical_transport(ᶜρ, ᶠu³, ᶜh_tot, FT(dt), Val(:none))
282-
@. Yₜ.c.ρe_tot += vtt - vtt_central
282+
# @. Yₜ.c.ρe_tot += vtt - vtt_central
283283
end
284284

285285
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)
289289
@. Yₜ.c.ρq_tot += vtt - vtt_central
290290
end
291291

292-
if isnothing(ᶠf¹²)
293-
# shallow atmosphere
294-
@. Yₜ.c.uₕ -=
295-
ᶜinterp(ᶠω¹² × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) / (Y.c.ρ * ᶜJ) +
296-
(ᶜf³ + ᶜω³) × CT12(ᶜu)
297-
@. Yₜ.f.u₃ -= ᶠω¹² × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
298-
for j in 1:n
299-
@. Yₜ.f.sgsʲs.:($$j).u₃ -=
300-
ᶠω¹²ʲs.:($$j) × ᶠinterp(CT12(ᶜuʲs.:($$j))) +
301-
ᶠgradᵥ(ᶜKʲs.:($$j) - ᶜinterp(ᶠKᵥʲs.:($$j)))
302-
end
303-
else
304-
# deep atmosphere
305-
@. Yₜ.c.uₕ -=
306-
ᶜinterp((ᶠf¹² + ᶠω¹²) × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) /
307-
(Y.c.ρ * ᶜJ) + (ᶜf³ + ᶜω³) × CT12(ᶜu)
308-
@. Yₜ.f.u₃ -= (ᶠf¹² + ᶠω¹²) × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
309-
for j in 1:n
310-
@. Yₜ.f.sgsʲs.:($$j).u₃ -=
311-
(ᶠf¹² + ᶠω¹²ʲs.:($$j)) × ᶠinterp(CT12(ᶜuʲs.:($$j))) +
312-
ᶠgradᵥ(ᶜKʲs.:($$j) - ᶜinterp(ᶠKᵥʲs.:($$j)))
313-
end
314-
end
315-
316-
if use_prognostic_tke(turbconv_model) # advect_tke triggers allocations
317-
@. ᶜa_scalar = ᶜtke⁰ * draft_area(ᶜρa⁰, ᶜρ⁰)
318-
vtt = vertical_transport(ᶜρ⁰, ᶠu³⁰, ᶜa_scalar, dt, edmfx_mse_q_tot_upwinding)
319-
@. Yₜ.c.sgs⁰.ρatke += vtt
320-
end
292+
# if isnothing(ᶠf¹²)
293+
# # shallow atmosphere
294+
# @. Yₜ.c.uₕ -=
295+
# ᶜinterp(ᶠω¹² × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) / (Y.c.ρ * ᶜJ) +
296+
# (ᶜf³ + ᶜω³) × CT12(ᶜu)
297+
# @. Yₜ.f.u₃ -= ᶠω¹² × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
298+
# for j in 1:n
299+
# @. Yₜ.f.sgsʲs.:($$j).u₃ -=
300+
# ᶠω¹²ʲs.:($$j) × ᶠinterp(CT12(ᶜuʲs.:($$j))) +
301+
# ᶠgradᵥ(ᶜKʲs.:($$j) - ᶜinterp(ᶠKᵥʲs.:($$j)))
302+
# end
303+
# else
304+
# # deep atmosphere
305+
# @. Yₜ.c.uₕ -=
306+
# ᶜinterp((ᶠf¹² + ᶠω¹²) × (ᶠinterp(Y.c.ρ * ᶜJ) * ᶠu³)) /
307+
# (Y.c.ρ * ᶜJ) + (ᶜf³ + ᶜω³) × CT12(ᶜu)
308+
# @. Yₜ.f.u₃ -= (ᶠf¹² + ᶠω¹²) × ᶠinterp(CT12(ᶜu)) + ᶠgradᵥ(ᶜK)
309+
# for j in 1:n
310+
# @. Yₜ.f.sgsʲs.:($$j).u₃ -=
311+
# (ᶠf¹² + ᶠω¹²ʲs.:($$j)) × ᶠinterp(CT12(ᶜuʲs.:($$j))) +
312+
# ᶠgradᵥ(ᶜKʲs.:($$j) - ᶜinterp(ᶠKᵥʲs.:($$j)))
313+
# end
314+
# end
315+
316+
# if use_prognostic_tke(turbconv_model) # advect_tke triggers allocations
317+
# @. ᶜa_scalar = ᶜtke⁰ * draft_area(ᶜρa⁰, ᶜρ⁰)
318+
# vtt = vertical_transport(ᶜρ⁰, ᶠu³⁰, ᶜa_scalar, dt, edmfx_mse_q_tot_upwinding)
319+
# @. Yₜ.c.sgs⁰.ρatke += vtt
320+
# end
321321
end
322322

323323
"""

‎src/prognostic_equations/dss.jl‎

Lines changed: 17 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -77,8 +77,22 @@ dss!(Y_state, params, t_current)
7777
"""
7878

7979
NVTX.@annotate function dss!(Y, p, t)
80-
if do_dss(axes(Y.c))
81-
Spaces.weighted_dss!(Y.c => p.ghost_buffer.c, Y.f => p.ghost_buffer.f)
82-
end
80+
# if do_dss(axes(Y.c))
81+
# Spaces.weighted_dss!(Y.c => p.ghost_buffer.c, Y.f => p.ghost_buffer.f)
82+
# end
83+
prescribe_flow!(Y, p, t, p.atmos.prescribed_flow)
8384
return nothing
8485
end
86+
87+
prescribe_flow!(Y, p, t, ::Nothing) = nothing
88+
function prescribe_flow!(Y, p, t, flow::PrescribedFlow)
89+
FT = eltype(p.params)
90+
# (; uₕ, u₃, ρ, p) = flow
91+
(; prescribed_u₃) = flow
92+
lg = Fields.local_geometry_field(Y.f)
93+
@. Y.f.u₃ = C3(Geometry.WVector(prescribed_u₃(FT, t)), lg)
94+
# @. Y.c.uₕ = uₕ(t)
95+
# @. Y.c.ρ = ρ₀ # move to IC
96+
# @. Y.c.p = p₀ # move to IC
97+
return nothing
98+
end

‎src/prognostic_equations/remaining_tendency.jl‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -64,10 +64,10 @@ Returns:
6464
NVTX.@annotate function remaining_tendency!(Yₜ, Yₜ_lim, Y, p, t)
6565
Yₜ_lim .= zero(eltype(Yₜ_lim))
6666
Yₜ .= zero(eltype(Yₜ))
67-
horizontal_tracer_advection_tendency!(Yₜ_lim, Y, p, t)
67+
# horizontal_tracer_advection_tendency!(Yₜ_lim, Y, p, t)
6868
fill_with_nans!(p) # TODO: would be better to limit this to debug mode (e.g., if p.debug_mode...)
69-
horizontal_dynamics_tendency!(Yₜ, Y, p, t)
70-
hyperdiffusion_tendency!(Yₜ, Yₜ_lim, Y, p, t)
69+
# horizontal_dynamics_tendency!(Yₜ, Y, p, t)
70+
# hyperdiffusion_tendency!(Yₜ, Yₜ_lim, Y, p, t)
7171
explicit_vertical_advection_tendency!(Yₜ, Y, p, t)
7272
vertical_advection_of_water_tendency!(Yₜ, Y, p, t)
7373
additional_tendency!(Yₜ, Y, p, t)

‎src/solver/type_getters.jl‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -114,6 +114,18 @@ function get_atmos(config::AtmosConfig, params)
114114
FT,
115115
)
116116

117+
if parsed_args["prescribed_flow"]
118+
function prescribed_u₃(FT::Type{<:Real}, _t)
119+
t = FT(_t)
120+
w1 = FT(2)
121+
t1 = FT(600) # 10 minutes
122+
return t < t1 ? w1 * sin(π * t / t1) : FT(0)
123+
end
124+
prescribed_flow = PrescribedFlow(; prescribed_u₃)
125+
else
126+
prescribed_flow = nothing
127+
end
128+
117129
atmos = AtmosModel(;
118130
# AtmosWater - Moisture, Precipitation & Clouds
119131
moisture_model,
@@ -130,6 +142,9 @@ function get_atmos(config::AtmosConfig, params)
130142
advection_test,
131143
scm_coriolis = get_scm_coriolis(parsed_args, FT),
132144

145+
# PrescribedFlow
146+
prescribed_flow,
147+
133148
# AtmosRadiation
134149
radiation_mode = final_radiation_mode,
135150
ozone,
@@ -433,6 +448,8 @@ function get_initial_condition(parsed_args, atmos)
433448
parsed_args["start_date"],
434449
parsed_args["era5_initial_condition_dir"],
435450
)
451+
elseif parsed_args["initial_condition"] == "ShipwayHill2012"
452+
return ICs.ShipwayHill2012()
436453
else
437454
error(
438455
"Unknown `initial_condition`: $(parsed_args["initial_condition"])",

‎src/solver/types.jl‎

Lines changed: 12 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,11 @@ import Dates
66
import ClimaUtilities.ClimaArtifacts: @clima_artifact
77
import LazyArtifacts
88

9+
abstract type AbstractPrescribedFlow end
10+
@kwdef struct PrescribedFlow
11+
prescribed_u₃::Function
12+
end
13+
914
abstract type AbstractMoistureModel end
1015
abstract type AbstractMoistModel <: AbstractMoistureModel end
1116
struct DryModel <: AbstractMoistureModel end
@@ -681,11 +686,12 @@ Base.broadcastable(x::AtmosGravityWave) = tuple(x)
681686
Base.broadcastable(x::AtmosSponge) = tuple(x)
682687
Base.broadcastable(x::AtmosSurface) = tuple(x)
683688

684-
struct AtmosModel{W, SCM, R, TC, GW, VD, SP, SU, NU}
689+
struct AtmosModel{W, SCM, R, TC, PF, GW, VD, SP, SU, NU}
685690
water::W
686691
scm_setup::SCM
687692
radiation::R
688693
turbconv::TC
694+
prescribed_flow::PF
689695
gravity_wave::GW
690696
vertical_diffusion::VD
691697
sponge::SP
@@ -701,6 +707,7 @@ const ATMOS_MODEL_GROUPS = (
701707
(AtmosWater, :water),
702708
(AtmosRadiation, :radiation),
703709
(AtmosTurbconv, :turbconv),
710+
(PrescribedFlow, :prescribed_flow),
704711
(AtmosGravityWave, :gravity_wave),
705712
(AtmosSponge, :sponge),
706713
(AtmosSurface, :surface),
@@ -918,11 +925,14 @@ function AtmosModel(; kwargs...)
918925
disable_surface_flux_tendency =
919926
get(atmos_model_kwargs, :disable_surface_flux_tendency, false)
920927

928+
prescribed_flow = get(atmos_model_kwargs, :prescribed_flow, nothing)
929+
921930
return AtmosModel{
922931
typeof(water),
923932
typeof(scm_setup),
924933
typeof(radiation),
925934
typeof(turbconv),
935+
typeof(prescribed_flow),
926936
typeof(gravity_wave),
927937
typeof(vertical_diffusion),
928938
typeof(sponge),
@@ -933,6 +943,7 @@ function AtmosModel(; kwargs...)
933943
scm_setup,
934944
radiation,
935945
turbconv,
946+
prescribed_flow,
936947
gravity_wave,
937948
vertical_diffusion,
938949
sponge,

0 commit comments

Comments
 (0)