From 9bd2a6926ae8b2aaf5e29a2a14c6e4a9efaebb0c Mon Sep 17 00:00:00 2001 From: Haakon Ludvig Langeland Ervik <45243236+haakon-e@users.noreply.github.com> Date: Mon, 26 Aug 2024 16:42:31 -0700 Subject: [PATCH] Forcing for ISDAC --- .buildkite/pipeline.yml | 13 +++++ config/default_configs/default_config.yml | 4 +- config/model_configs/les_isdac_box.yml | 30 ++++++++++ post_processing/ci_plots.jl | 1 + src/initial_conditions/initial_conditions.jl | 46 ++++++++++++++++ .../radiation/radiation.jl | 29 ++++++++++ .../forcing/external_forcing.jl | 55 +++++++++++++++++++ src/solver/model_getters.jl | 32 ++++++++++- src/solver/type_getters.jl | 6 ++ src/solver/types.jl | 9 +++ src/surface_conditions/surface_setups.jl | 11 ++++ 11 files changed, 233 insertions(+), 3 deletions(-) create mode 100644 config/model_configs/les_isdac_box.yml diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 9625a31c7c..f0540cbcdf 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -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)" diff --git a/config/default_configs/default_config.yml b/config/default_configs/default_config.yml index 57c84b544b..1949635ec4 100644 --- a/config/default_configs/default_config.yml +++ b/config/default_configs/default_config.yml @@ -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)]" @@ -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" diff --git a/config/model_configs/les_isdac_box.yml b/config/model_configs/les_isdac_box.yml new file mode 100644 index 0000000000..255d01dba2 --- /dev/null +++ b/config/model_configs/les_isdac_box.yml @@ -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 diff --git a/post_processing/ci_plots.jl b/post_processing/ci_plots.jl index 607118eaab..4dbc4efa1a 100644 --- a/post_processing/ci_plots.jl +++ b/post_processing/ci_plots.jl @@ -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}, } diff --git a/src/initial_conditions/initial_conditions.jl b/src/initial_conditions/initial_conditions.jl index e5a5c0f064..c3f9edef99 100644 --- a/src/initial_conditions/initial_conditions.jl +++ b/src/initial_conditions/initial_conditions.jl @@ -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 diff --git a/src/parameterized_tendencies/radiation/radiation.jl b/src/parameterized_tendencies/radiation/radiation.jl index bbb90d8f53..f7adb51b35 100644 --- a/src/parameterized_tendencies/radiation/radiation.jl +++ b/src/parameterized_tendencies/radiation/radiation.jl @@ -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 diff --git a/src/prognostic_equations/forcing/external_forcing.jl b/src/prognostic_equations/forcing/external_forcing.jl index 6006c6b960..798fbbcfc8 100644 --- a/src/prognostic_equations/forcing/external_forcing.jl +++ b/src/prognostic_equations/forcing/external_forcing.jl @@ -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 diff --git a/src/solver/model_getters.jl b/src/solver/model_getters.jl index 93093f14d4..3f12da256b 100644 --- a/src/solver/model_getters.jl +++ b/src/solver/model_getters.jl @@ -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) @@ -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 @@ -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 @@ -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" @@ -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 @@ -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 diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 5940a788df..ac4ccdc40f 100644 --- a/src/solver/type_getters.jl +++ b/src/solver/type_getters.jl @@ -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) @@ -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", diff --git a/src/solver/types.jl b/src/solver/types.jl index 552c260fb9..f887e30802 100644 --- a/src/solver/types.jl +++ b/src/solver/types.jl @@ -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 @@ -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} diff --git a/src/surface_conditions/surface_setups.jl b/src/surface_conditions/surface_setups.jl index 47c588d458..a07755d319 100644 --- a/src/surface_conditions/surface_setups.jl +++ b/src/surface_conditions/surface_setups.jl @@ -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