Forcing for ISDAC

.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)"

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"

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

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},
 }
 
 

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

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

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

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

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",

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}

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