Skip to content

Commit

Permalink
Merge pull request #3025 from CliMA/he/isdac-les
Browse files Browse the repository at this point in the history
Forcing for ISDAC
  • Loading branch information
haakon-e authored Aug 27, 2024
2 parents 11e980d + 9bd2a69 commit bd85668
Show file tree
Hide file tree
Showing 11 changed files with 233 additions and 3 deletions.
13 changes: 13 additions & 0 deletions .buildkite/pipeline.yml
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,19 @@ steps:
agents:
slurm_mem: 20GB

- label: ":genie: LES ISDAC in a box"
command: >
julia --color=yes --project=examples examples/hybrid/driver.jl
--config_file $CONFIG_PATH/les_isdac_box.yml
--job_id les_isdac_box
artifact_paths: "les_isdac_box/*"
env:
CLIMACOMMS_DEVICE: "CUDA"
agents:
slurm_mem: 16G
slurm_gpus: 1
soft_fail: true

- group: "Plane Examples"
steps:
- label: ":computer: Agnesi linear hydrostatic mountain experiment (uniform)"
Expand Down
4 changes: 2 additions & 2 deletions config/default_configs/default_config.yml
Original file line number Diff line number Diff line change
Expand Up @@ -179,7 +179,7 @@ surface_temperature:
help: "Prescribed surface temperature functional form ['ZonallySymmetric' (default), 'ZonallyAsymmetric', 'RCEMIPII']"
value: "ZonallySymmetric"
initial_condition:
help: "Initial condition [`DryBaroclinicWave`, `MoistBaroclinicWave`, `DecayingProfile`, `IsothermalProfile`, `Bomex`, `DryDensityCurrentProfile`, `AgnesiHProfile`, `ScharProfile`, `RisingThermalBubbleProfile`]"
help: "Initial condition [`DryBaroclinicWave`, `MoistBaroclinicWave`, `DecayingProfile`, `IsothermalProfile`, `Bomex`, `DryDensityCurrentProfile`, `AgnesiHProfile`, `ScharProfile`, `RisingThermalBubbleProfile`, `ISDAC`]"
value: "DecayingProfile"
perturb_initstate:
help: "Add a perturbation to the initial condition [`false`, `true` (default)]"
Expand Down Expand Up @@ -236,7 +236,7 @@ cfsite_number:
help: "cfsite identifier for single column forcing from `external_forcing_file`, specified as siteN. For site details see Shen et al. 2022 `https://doi.org/10.1029/2021MS002631`. [`site23` (default), `siteN`]"
value: "site23"
subsidence:
help: "Subsidence [`nothing` (default), `Bomex`, `LifeCycleTan2018`, `Rico`, `DYCOMS`]"
help: "Subsidence [`nothing` (default), `Bomex`, `LifeCycleTan2018`, `Rico`, `DYCOMS`, `ISDAC`]"
value: ~
toml:
help: "TOML file(s) used to override model parameters"
Expand Down
30 changes: 30 additions & 0 deletions config/model_configs/les_isdac_box.yml
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
job_id: "les_isdac_box"
initial_condition: "ISDAC"
subsidence: "ISDAC"
surface_setup: "ISDAC"
external_forcing: "ISDAC"
rad: "ISDAC"
moist: "equil"
config: "box"
hyperdiff: "false"
implicit_diffusion: false
ode_algo: "SSP33ShuOsher"
precip_model: "1M"
smagorinsky_lilly: true
c_smag: 0.20
x_max: 3.2e3
y_max: 3.2e3
z_max: 2e3
x_elem: 4
y_elem: 4
z_elem: 10
z_stretch: false
dt: "0.05secs"
t_end: "10mins"
dt_save_state_to_disk: "1mins"
# restart_file: "restart/isdac/day0.0.hdf5"
netcdf_interpolation_num_points: [10, 10, 20]
diagnostics:
- short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl]
reduction: average
period: 1mins
1 change: 1 addition & 0 deletions post_processing/ci_plots.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1094,6 +1094,7 @@ EDMFBoxPlotsWithPrecip = Union{
Val{:diagnostic_edmfx_rico_box},
Val{:diagnostic_edmfx_trmm_box},
Val{:diagnostic_edmfx_trmm_stretched_box},
Val{:les_isdac_box},
}


Expand Down
46 changes: 46 additions & 0 deletions src/initial_conditions/initial_conditions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1196,3 +1196,49 @@ function gcm_initial_conditions(external_forcing_file, cfsite_number)
)
end
end

Base.@kwdef struct ISDAC <: InitialCondition
prognostic_tke::Bool = false
perturb::Bool = false
end

function (initial_condition::ISDAC)(params)
(; prognostic_tke, perturb) = initial_condition
FT = eltype(params)
thermo_params = CAP.thermodynamics_params(params)
p_0 = FT(102000) # 1020 hPa
θ = APL.ISDAC_θ_liq_ice(FT) # K
q_tot = APL.ISDAC_q_tot(FT) # kg/kg
# Note: ISDAC top-of-domain is ~1.5km, but we don't have access to that information here, so we use 5km to be safe
p = hydrostatic_pressure_profile(;
thermo_params,
p_0,
θ,
q_tot,
z_max = 5000,
) # Pa

u = APL.ISDAC_u(FT) # m/s
v = APL.ISDAC_v(FT) # m/s
tke = APL.ISDAC_tke(FT) # m²/s²

# pseudorandom fluctuations with amplitude 0.1 K
θ_pert(z::FT) where {FT} =
perturb && (z < 825) ? FT(0.1) * randn(FT) : FT(0)

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) + θ_pert(z),
q_tot(z),
),
velocity = Geometry.UVVector(u(z), v(z)),
turbconv_state = EDMFState(; tke = prognostic_tke ? tke(z) : FT(0)),
)
end
end
29 changes: 29 additions & 0 deletions src/parameterized_tendencies/radiation/radiation.jl
Original file line number Diff line number Diff line change
Expand Up @@ -363,3 +363,32 @@ function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationTRMM_LBA)
@. Yₜ.c.ρe_tot += ᶜρ * TD.cv_m(thermo_params, ᶜts_gm) * ᶜdTdt_rad
return nothing
end

#####
##### ISDAC radiation
#####

radiation_model_cache(Y, radiation_mode::RadiationISDAC; args...) = (;) # Don't need a cache for ISDAC
function radiation_tendency!(Yₜ, Y, p, t, radiation_mode::RadiationISDAC)
(; F₀, F₁, κ) = radiation_mode
(; params, precomputed) = p
(; ᶜts) = precomputed
thermo_params = CAP.thermodynamics_params(params)

ᶜρq = p.scratch.ᶜtemp_scalar
@. ᶜρq = Y.c.ρ * TD.liquid_specific_humidity(thermo_params, ᶜts)

LWP_zₜ = p.scratch.temp_field_level # column integral of LWP (zₜ = top-of-domain)
Operators.column_integral_definite!(LWP_zₜ, ᶜρq)

LWP_z = p.scratch.ᶠtemp_scalar # column integral of LWP from 0 to z (z = current level)
Operators.column_integral_indefinite!(LWP_z, ᶜρq)

@. Yₜ.c.ρe_tot -= ᶜdivᵥ(
Geometry.WVector(
F₀ * exp(-κ * (LWP_zₜ - LWP_z)) + F₁ * exp(-κ * LWP_z),
),
) # = -∂F/∂z = ρ cₚ ∂T/∂t (longwave radiation)

return nothing
end
55 changes: 55 additions & 0 deletions src/prognostic_equations/forcing/external_forcing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -227,3 +227,58 @@ function external_forcing_tendency!(Yₜ, Y, p, t, ::GCMForcing)

return nothing
end

# ISDAC external forcing (i.e. nudging)
external_forcing_cache(Y, external_forcing::ISDACForcing, params) = (;) # Don't need to cache anything
function external_forcing_tendency!(Yₜ, Y, p, t, ::ISDACForcing)
FT = Spaces.undertype(axes(Y.c))
(; params) = p
thermo_params = CAP.thermodynamics_params(params)
(; ᶜspecific, ᶜts, ᶜh_tot, ᶜp) = p.precomputed

ᶜinv_τ_scalar = APL.ISDAC_inv_τ_scalar(FT) # s⁻¹
ᶜinv_τ_wind = APL.ISDAC_inv_τ_wind(FT) # s⁻¹
θ = APL.ISDAC_θ_liq_ice(FT)
u = APL.ISDAC_u(FT)
v = APL.ISDAC_v(FT)
q_tot = APL.ISDAC_q_tot(FT)

# Convert ISDAC potential temperature to air temperature
ta_ISDAC =
(pres, z) -> TD.air_temperature(
thermo_params,
TD.PhaseEquil_pθq(thermo_params, pres, θ(z), q_tot(z)),
)

ᶜz = Fields.coordinate_field(Y.c).z
ᶜlg = Fields.local_geometry_field(Y.c)
ᶜuₕ_nudge = p.scratch.ᶜtemp_C12
@. ᶜuₕ_nudge = C12(Geometry.UVVector(u(ᶜz), v(ᶜz)), ᶜlg)
@. Yₜ.c.uₕ -= (Y.c.uₕ - ᶜuₕ_nudge) * ᶜinv_τ_wind(ᶜz)

# TODO: May make more sense to use initial ISDAC (hydrostatic) pressure, but would need to add it to cache,
# so for now just use current pressure.
ᶜdTdt_nudging = p.scratch.ᶜtemp_scalar
ᶜdqtdt_nudging = p.scratch.ᶜtemp_scalar_2
@. ᶜdTdt_nudging =
-(TD.air_temperature(thermo_params, ᶜts) - ta_ISDAC(ᶜp, ᶜz)) *
ᶜinv_τ_scalar(ᶜz)
@. ᶜdqtdt_nudging = -(ᶜspecific.q_tot - q_tot(ᶜz)) * ᶜinv_τ_scalar(ᶜz)

T_0 = TD.Parameters.T_0(thermo_params)
Lv_0 = TD.Parameters.LH_v0(thermo_params)
cv_v = TD.Parameters.cv_v(thermo_params)
R_v = TD.Parameters.R_v(thermo_params)
# total energy
@. Yₜ.c.ρe_tot +=
Y.c.ρ * (
TD.cv_m(thermo_params, ᶜts) * ᶜdTdt_nudging +
(
cv_v * (TD.air_temperature(thermo_params, ᶜts) - T_0) + Lv_0 -
R_v * T_0
) * ᶜdqtdt_nudging
)

# total specific humidity
@. Yₜ.c.ρq_tot += Y.c.ρ * ᶜdqtdt_nudging
end
32 changes: 31 additions & 1 deletion src/solver/model_getters.jl
Original file line number Diff line number Diff line change
Expand Up @@ -235,6 +235,7 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
"allskywithclear",
"DYCOMS",
"TRMM_LBA",
"ISDAC",
)
return if radiation_name == "gray"
RRTMGPI.GrayRadiation(add_isothermal_boundary_layer)
Expand Down Expand Up @@ -262,6 +263,8 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
RadiationDYCOMS{FT}()
elseif radiation_name == "TRMM_LBA"
RadiationTRMM_LBA(FT)
elseif radiation_name == "ISDAC"
RadiationISDAC{FT}()
else
nothing
end
Expand Down Expand Up @@ -327,6 +330,8 @@ function get_subsidence_model(parsed_args, radiation_mode, FT)
elseif subsidence == "DYCOMS"
@assert radiation_mode isa RadiationDYCOMS
z -> -z * radiation_mode.divergence
elseif subsidence == "ISDAC"
APL.ISDAC_subsidence(FT)
else
error("Uncaught case")
end
Expand Down Expand Up @@ -372,7 +377,7 @@ end

function get_external_forcing_model(parsed_args)
external_forcing = parsed_args["external_forcing"]
@assert external_forcing in (nothing, "GCM")
@assert external_forcing in (nothing, "GCM", "ISDAC")
return if isnothing(external_forcing)
nothing
elseif external_forcing == "GCM"
Expand All @@ -381,6 +386,8 @@ function get_external_forcing_model(parsed_args)
parsed_args["external_forcing_file"],
parsed_args["cfsite_number"],
)
elseif external_forcing == "ISDAC"
ISDACForcing()
end
end

Expand Down Expand Up @@ -491,3 +498,26 @@ function get_tendency_model(parsed_args)
UseAllTendency()
end
end

function check_case_consistency(parsed_args)
# if any flags is ISDAC, check that all are ISDAC
ic = parsed_args["initial_condition"]
subs = parsed_args["subsidence"]
surf = parsed_args["surface_setup"]
rad = parsed_args["rad"]
cor = parsed_args["edmf_coriolis"]
forc = parsed_args["forcing"]
moist = parsed_args["moist"]
ls_adv = parsed_args["ls_adv"]
extf = parsed_args["external_forcing"]

ISDAC_mandatory = (ic, subs, surf, rad, extf)
if "ISDAC" in ISDAC_mandatory
@assert(
allequal(ISDAC_mandatory) &&
all(isnothing, (cor, forc, ls_adv)) &&
moist != "dry",
"ISDAC setup not consistent"
)
end
end
6 changes: 6 additions & 0 deletions src/solver/type_getters.jl
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@ function get_atmos(config::AtmosConfig, params)
(; turbconv_params) = params
(; parsed_args) = config
FT = eltype(config)
check_case_consistency(parsed_args)
moisture_model = get_moisture_model(parsed_args)
precip_model = get_precipitation_model(parsed_args)
cloud_model = get_cloud_model(parsed_args)
Expand Down Expand Up @@ -344,6 +345,11 @@ function get_initial_condition(parsed_args)
return getproperty(ICs, Symbol(parsed_args["initial_condition"]))(
parsed_args["prognostic_tke"],
)
elseif parsed_args["initial_condition"] == "ISDAC"
ICs.ISDAC(
parsed_args["prognostic_tke"],
parsed_args["perturb_initstate"],
)
elseif parsed_args["initial_condition"] in [
"IsothermalProfile",
"AgnesiHProfile",
Expand Down
9 changes: 9 additions & 0 deletions src/solver/types.jl
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,8 @@ struct GCMForcing{FT}
cfsite_number::String
end

struct ISDACForcing end

struct EDMFCoriolis{U, V, FT}
prof_ug::U
prof_vg::V
Expand Down Expand Up @@ -275,6 +277,13 @@ Base.@kwdef struct RadiationDYCOMS{FT}
F0::FT = 70.0
F1::FT = 22.0
end

Base.@kwdef struct RadiationISDAC{FT}
F₀::FT = 72 # W/m²
F₁::FT = 15 # W/m²
κ::FT = 170 # m²/kg
end

import AtmosphericProfilesLibrary as APL

struct RadiationTRMM_LBA{R}
Expand Down
11 changes: 11 additions & 0 deletions src/surface_conditions/surface_setups.jl
Original file line number Diff line number Diff line change
Expand Up @@ -269,3 +269,14 @@ function gcm_surface_conditions(external_forcing_file, cfsite_number)
mean(gcm_driven_timeseries(ds.group[cfsite_number], "ts"))
end
end

struct ISDAC end
function (::ISDAC)(params)
FT = eltype(params)
T = FT(267) # K
p = FT(102000) # Pa
lhf = shf = FT(0) # neglect heat fluxes
z0 = FT(4e-4) # m surface roughness length
parameterization = MoninObukhov(; z0, fluxes = HeatFluxes(; lhf, shf))
return SurfaceState(; parameterization, T, p)
end

0 comments on commit bd85668

Please sign in to comment.