Forcing for ISDAC

.buildkite/pipeline.yml

@@ -158,6 +158,16 @@ steps:
           --job_id box_density_current_test
         artifact_paths: "box_density_current_test/output_active/*"
 
+      - 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/*"
+        agents:
+          slurm_mem: 16G
+          slurm_gpus: 1
+
   - group: "Plane Examples"
     steps:
       - label: ":computer: Agnesi linear hydrostatic mountain experiment (uniform)"

config/model_configs/les_isdac_box.yml

@@ -0,0 +1,31 @@
+device: CPUSingleThreaded
+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: 64
+y_elem: 64
+z_elem: 200
+z_stretch: false
+dt: "0.1secs" 
+t_end: "1hour" 
+dt_save_state_to_disk: "10mins"
+# restart_file: "restart/isdac/day0.0.hdf5"
+netcdf_interpolation_num_points: [64, 64, 200]
+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

@@ -1087,6 +1087,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

@@ -72,7 +72,7 @@ function column_indefinite_integral(
         boundary_names = (:bottom, :top),
     )
     z_mesh = Meshes.IntervalMesh(z_domain; nelems)
-    context = ClimaComms.SingletonCommsContext()
+    context = ClimaComms.SingletonCommsContext(ClimaComms.CPUSingleThreaded())  # TODO: ClimaComms v0.6 fixes this
     z_topology = Topologies.IntervalTopology(context, z_mesh)
     cspace = Spaces.CenterFiniteDifferenceSpace(z_topology)
     fspace = Spaces.FaceFiniteDifferenceSpace(z_topology)
@@ -1183,3 +1183,48 @@ function gcm_initial_conditions(external_forcing_file)
         )
     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

@@ -343,3 +343,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

@@ -184,3 +184,57 @@ 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) = (;)  # Don't need to cache anything
+function external_forcing_tendency!(Yₜ, Y, p, t, ::ISDACForcing)
+    (; 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

@@ -215,6 +215,7 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
         "allskywithclear",
         "DYCOMS_RF01",
         "TRMM_LBA",
+        "ISDAC",
     )
     return if radiation_name == "clearsky"
         RRTMGPI.ClearSkyRadiation(
@@ -248,6 +249,8 @@ function get_radiation_mode(parsed_args, ::Type{FT}) where {FT}
         RadiationDYCOMS_RF01{FT}()
     elseif radiation_name == "TRMM_LBA"
         RadiationTRMM_LBA(FT)
+    elseif radiation_name == "ISDAC"
+        RadiationISDAC{FT}()
     else
         nothing
     end
@@ -311,6 +314,8 @@ function get_subsidence_model(parsed_args, radiation_mode, FT)
     elseif subsidence == "DYCOMS"
         @assert radiation_mode isa RadiationDYCOMS_RF01
         z -> -z * radiation_mode.divergence
+    elseif subsidence == "ISDAC"
+        APL.ISDAC_subsidence(FT)
     else
         error("Uncaught case")
     end
@@ -356,12 +361,14 @@ 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"
         DType = Float64  # TODO: Read from `parsed_args`
         GCMForcing{DType}(parsed_args["external_forcing_file"])
+    elseif external_forcing == "ISDAC"
+        ISDACForcing()
     end
 end
 
@@ -470,3 +477,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

@@ -17,6 +17,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)
@@ -343,6 +344,8 @@ 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

@@ -119,6 +119,8 @@ struct GCMForcing{FT}
     external_forcing_file::String
 end
 
+struct ISDACForcing end
+
 struct EDMFCoriolis{U, V, FT}
     prof_ug::U
     prof_vg::V
@@ -262,6 +264,13 @@ Base.@kwdef struct RadiationDYCOMS_RF01{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

@@ -268,3 +268,14 @@ function gcm_surface_conditions(external_forcing_file, imin)
         )
     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