Address review

	modified:   NEWS.md

.buildkite/pipeline.yml

@@ -281,6 +281,15 @@ steps:
         artifact_paths: "deep_sphere_baroclinic_wave_rhoe_equilmoist/output_active/*"
         agents:
           slurm_constraint: icelake|cascadelake|skylake|epyc
+
+      - label: ":computer: test DYAMOND interpolated initial conditions"
+        command: >
+          julia --color=yes --project=examples examples/hybrid/driver.jl
+          --config_file $GPU_CONFIG_PATH/gpu_aquaplanet_dyamond_summer.yml
+          --job_id gpu_aquaplanet_dyamond_summer
+        artifact_paths: "gpu_aquaplanet_dyamond_summer/output_active/*"
+        agents:
+          slurm_constraint: icelake|cascadelake|skylake|epyc
 
       # Add this back when we figure out what to do with SSP and zalesak
       # - label: ":computer: SSP zalesak tracer & energy upwind baroclinic wave (ρe_tot) equilmoist"

NEWS.md

@@ -5,6 +5,12 @@ Main
 -------
 ### Features
 
+### DYAMOND-summer initial conditions
+
+Added functionality to allow initial conditions to be overwritten by 
+interpolated datasets. Currently supports equilibrium moisture models
+for DYAMOND runs, assuming a simulation start date of 1 August 2016.
+
 ### Write diagnostics to text files
 
 Added functionality to write diagnostics in DictWriter to text files.

config/gpu_configs/gpu_aquaplanet_dyamond_summer.yml

@@ -17,14 +17,14 @@ vert_diff: "FriersonDiffusion"
 implicit_diffusion: true
 approximate_linear_solve_iters: 2
 surface_setup: "DefaultMoninObukhov"
-dt: "30secs"
-t_end: "120secs"
+dt: "5secs"
+t_end: "30secs"
 toml: [toml/longrun_aquaplanet.toml]
 prescribe_ozone: true
 aerosol_radiation: true
 prescribed_aerosols: ["CB1", "CB2", "DST01", "OC1", "OC2", "SO4", "SSLT01"]
 initial_condition: "DYAMONDSummer"
 topography: "Earth"
 diagnostics:
-  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]
+  - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr, rv]
     period: "30secs"

post_processing/ci_plots.jl

@@ -679,6 +679,7 @@ end
 
 SphereOrographyPlots = Union{
     Val{:sphere_baroclinic_wave_rhoe_topography_dcmip_rs},
+    Val{:gpu_aquaplanet_dyamond_summer},
     Val{:sphere_baroclinic_wave_rhoe_hughes2023},
 }
 

src/initial_conditions/initial_conditions.jl

@@ -161,10 +161,11 @@ function (initial_condition::DecayingProfile)(params)
 end
 
 """
-    DYAMONDSummer(; perturb = true)
+    DYAMONDSummer()
 
-An `InitialCondition` with a decaying temperature profile, and with an optional
-perturbation to the temperature.
+This function assigns a default initial condition, populating the `LocalState`
+with `NaN`, prior to applying the `overwrite_initial_conditions` method, which
+interpolates values from NetCDF data onto the `ExtrudedFiniteDifferenceSpace`.
 """
 struct DYAMONDSummer <: InitialCondition end
 
@@ -174,7 +175,7 @@ function (initial_condition::DYAMONDSummer)(params)
         grav = CAP.grav(params)
         thermo_params = CAP.thermodynamics_params(params)
 
-        T, p = FT(300), FT(100000) # placeholder values
+        T, p = FT(NaN), FT(NaN) # placeholder values
 
         return LocalState(;
             params,
@@ -376,22 +377,27 @@ function (initial_condition::RisingThermalBubbleProfile)(params)
     return local_state
 end
 
+"""
+    overwrite_initial_conditions!(initial_condition, args...)
+Do-nothing fallback method for the operation overwriting initial conditions 
+(this functionality required in instances where we interpolate initial conditions from NetCDF files). 
+Future work may revisit this design choice. 
+"""
 overwrite_initial_conditions!(initial_condition::InitialCondition, args...) =
     (return nothing)
 
 """
     overwrite_initial_conditions!(initial_condition, Y, thermo_params, config)
 Given a prognostic state `Y`, an `initial condition` (specifically, where
-initial values are assigned from interpolations of existing datasets), a `thermo_state`,
-and the `config::AtmosConfig` object, this function overwrites the default initial condition
+initial values are assigned from interpolations of existing datasets), a `thermo_state`, this function overwrites the default initial condition
 and populates prognostic variables with interpolated values using the `SpaceVaryingInputs`
-tool. 
+tool. To mitigate issues related to unbalanced states following the interpolation operation, 
+we recompute vertical pressure levels assuming hydrostatic balance, given the surface pressure
 """
 function overwrite_initial_conditions!(
     initial_condition::DYAMONDSummer,
     Y,
     thermo_params,
-    config,
 )
     # Get file from AtmosArtifacts
     file_path = AA.dyamond_summer_artifact_path(; context = config.comms_ctx)
@@ -405,14 +411,29 @@ function overwrite_initial_conditions!(
     )
     ᶜT = SpaceVaryingInputs.SpaceVaryingInput(file_path, "t", center_space)
     ᶜq_tot = SpaceVaryingInputs.SpaceVaryingInput(file_path, "q", center_space)
+
+    # With the known temperature (ᶜT) and moisture (ᶜq_tot) profile, 
+    # recompute the pressure levels assuming hydrostatic balance is maintained.
+    # Uses the ClimaCore `column_integral_indefinite!` function to solve 
+    # ∂(ln𝑝)/∂z = -g/(Rₘ(q)T), where
+    # p is the local pressure
+    # g is the gravitational constant
+    # q is the specific humidity
+    # Rₘ is the gas constant for moist air
+    # T is the air temperature
+    # p is then updated with the integral result, given p_sfc,
+    # following which the thermodynamic state is constructed. 
     ᶜ∂lnp∂z = @. -thermo_params.grav /
        (TD.gas_constant_air(thermo_params, TD.PhasePartition(ᶜq_tot)) * ᶜT)
     ᶠlnp_over_psfc = zeros(face_space)
     Operators.column_integral_indefinite!(ᶠlnp_over_psfc, ᶜ∂lnp∂z)
     ᶠp = p_sfc .* exp.(ᶠlnp_over_psfc)
     ᶜts = TD.PhaseEquil_pTq.(thermo_params, ᶜinterp.(ᶠp), ᶜT, ᶜq_tot)
+
     # Assign prognostic variables from equilibrium moisture models
     Y.c.ρ .= TD.air_density.(thermo_params, ᶜts)
+    # Velocity is first assigned on cell-centers and then interpolated onto
+    # cell faces.
     vel =
         Geometry.UVWVector.(
             SpaceVaryingInputs.SpaceVaryingInput(file_path, "u", center_space),
@@ -425,12 +446,14 @@ function overwrite_initial_conditions!(
     compute_kinetic!(e_kin, Y.c.uₕ, Y.f.u₃)
     e_pot = Fields.coordinate_field(Y.c).z .* thermo_params.grav
     Y.c.ρe_tot .= TD.total_energy.(thermo_params, ᶜts, e_kin, e_pot) .* Y.c.ρ
-    if config.parsed_args["moist"] == "dry"
-        error("`dry` configurations are incompatible with the interpolated DYAMOND-summer initial conditions.")
-    else
+    if hasproperty(Y.c, :ρq_tot)
         Y.c.ρq_tot .= ᶜq_tot .* Y.c.ρ
+    else
+        error(
+            "`dry` configurations are incompatible with the interpolated initial conditions.",
+        )
     end
-    if config.parsed_args["precip_model"] == "1M"
+    if hasproperty(Y.c, :ρq_sno) && hasproperty(Y.c, :ρq_rai)
         Y.c.ρq_sno .=
             SpaceVaryingInputs.SpaceVaryingInput(
                 file_path,

src/solver/type_getters.jl

@@ -265,8 +265,6 @@ function get_initial_condition(parsed_args)
         return getproperty(ICs, Symbol(parsed_args["initial_condition"]))(
             parsed_args["perturb_initstate"],
         )
-    elseif parsed_args["initial_condition"] in ["DYAMONDSummer"]
-        return getproperty(ICs, Symbol(parsed_args["initial_condition"]))()
     elseif parsed_args["initial_condition"] in [
         "Nieuwstadt",
         "GABLS",
@@ -296,6 +294,7 @@ function get_initial_condition(parsed_args)
         "RisingThermalBubbleProfile",
         "ScharProfile",
         "PrecipitatingColumn",
+        "DYAMONDSummer",
     ]
         return getproperty(ICs, Symbol(parsed_args["initial_condition"]))()
     elseif parsed_args["initial_condition"] == "GCM"
@@ -668,11 +667,16 @@ function get_simulation(config::AtmosConfig)
 
     tracers = get_tracers(config.parsed_args)
 
+    # In instances where we wish to interpolate existing datasets,
+    # e.g. NetCDF files containing spatially varying thermodynamic properties, 
+    # this call to `overwrite_initial_conditions`
+    # accesses the appropriate artifact for a given `initial_condition`
+    # and then updates some predetermined state `Y` (specific to `initial_condition`) 
+    # with those computed using the `SpaceVaryingInputs` tool.
     CA.InitialConditions.overwrite_initial_conditions!(
         initial_condition,
         Y,
         params.thermodynamics_params,
-        config,
     )
 
     s = @timed_str begin

src/utils/AtmosArtifacts.jl

@@ -82,6 +82,13 @@ function earth_orography_file_path(; context = nothing)
     )
 end
 
+"""
+    dyamond_summer_artifact_path(; context = nothing)
+
+Construct the file path for the 0.98deg DYAMOND initial conditions: 
+the artifact contains initial conditions for DYAMOND simulations 
+beginning on 1 August 2016.
+"""
 function dyamond_summer_artifact_path(; context = nothing)
     filename = "DYAMOND_SUMMER_ICS_p98deg.nc"
     return joinpath(