Merge pull request #3458 from CliMA/sk/dyamond_summer_ics

`DYAMOND Summer` low-resolution simulation

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

Artifacts.toml

@@ -43,3 +43,11 @@ git-tree-sha1 = "10742e0a2e343d13bb04df379e300a83402d4106"
     [[era5_cloud.download]]
     sha256 = "bb51e2f2d315b487e05a8d38944d4ad937ee4a40c43b68541220c5d54425e24a"
     url = "https://caltech.box.com/shared/static/b6ur4ap4vo04j09vdulem96z9fxqlgyn.gz"
+
+[DYAMOND_SUMMER_ICS_p98deg]
+git-tree-sha1 = "46aa28a9177dd7e4cbb813eeba3e02597c3ba071"
+lazy = true
+
+    [[DYAMOND_SUMMER_ICS_p98deg.download]]
+    sha256 = "9a5efd20d68d9b954af2fd7803bccdda3fc7ef4ff60421ed13d18f768312d2e3"
+    url = "https://caltech.box.com/shared/static/whn01xr83v0s4p8tc59ds3nvfsfe6ebs.gz"

NEWS.md

@@ -5,10 +5,29 @@ Main
 -------
 ### Features
 
+### Read initial conditions from NetCDF files
+
+Added functionality to allow initial conditions to be overwritten by
+interpolated NetCDF datasets.
+
+To use this feature from the YAML interface, just pass the path of the file.
+We expect the file to contain the following variables:
+- `p`, for pressure,
+- `t`, for temperature,
+- `q`, for humidity,
+- `u, v, w`, for velocity,
+- `cswc, crwc` for snow and rain water content (for 1 moment microphysics).
+
+For example, to use the DYAMONDSummer initial condition, set
+```
+initial_condition: "artifact\"DYAMONDSummer\"/DYAMOND_SUMMER_ICS_p98deg.nc"
+```
+in your configuration file.
+
 ### Write diagnostics to text files
 
 Added functionality to write diagnostics in DictWriter to text files.
-This is useful for outputing scalar diagnostics, such as total mass of
+This is useful for outputting scalar diagnostics, such as total mass of
 the atmosphere. PR [3476](https://github.com/CliMA/ClimaAtmos.jl/pull/3476)
 
 v0.28.0

config/default_configs/default_config.yml

@@ -182,7 +182,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`, `ISDAC`]"
+  help: "Initial condition [`DryBaroclinicWave`, `MoistBaroclinicWave`, `DecayingProfile`, `IsothermalProfile`, `Bomex`, `DryDensityCurrentProfile`, `AgnesiHProfile`, `ScharProfile`, `RisingThermalBubbleProfile`, `ISDAC`], or a file path for a NetCDF file (read documentation about requirements)."
   value: "DecayingProfile"
 perturb_initstate:
   help: "Add a perturbation to the initial condition [`false`, `true` (default)]"

config/gpu_configs/gpu_aquaplanet_dyamond_summer.yml

@@ -0,0 +1,31 @@
+dt_save_state_to_disk: "Inf"
+dt_save_to_sol: "Inf"
+output_default_diagnostics: false
+h_elem: 16
+z_max: 60000.0
+z_elem: 63
+dz_bottom: 30.0
+rayleigh_sponge: true
+viscous_sponge: true
+moist: "equil"
+precip_model: "1M"
+rad: "allskywithclear"
+insolation: "timevarying"
+dt_rad: "1hours"
+dt_cloud_fraction: "1hours"
+vert_diff: "FriersonDiffusion"
+implicit_diffusion: true
+approximate_linear_solve_iters: 2
+surface_setup: "DefaultMoninObukhov"
+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"]
+start_date: "20160801"
+initial_condition: "artifact\"DYAMOND_SUMMER_ICS_p98deg\"/DYAMOND_SUMMER_ICS_p98deg.nc"
+topography: "Earth"
+diagnostics:
+  - 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/cache/tracer_cache.jl

@@ -1,4 +1,5 @@
 import Dates: Year
+import ClimaUtilities.TimeVaryingInputs
 import ClimaUtilities.TimeVaryingInputs:
     TimeVaryingInput, LinearPeriodFillingInterpolation
 import Interpolations as Intp
@@ -14,7 +15,10 @@ function ozone_cache(::PrescribedOzone, Y, start_date)
         reference_date = start_date,
         regridder_type = :InterpolationsRegridder,
         regridder_kwargs = (; extrapolation_bc),
-        method = LinearPeriodFillingInterpolation(Year(1)),
+        method = LinearPeriodFillingInterpolation(
+            Year(1),
+            TimeVaryingInputs.Flat(),
+        ),
     )
     return (; o3, prescribed_o3_timevaryinginput)
 end

src/initial_conditions/InitialConditions.jl

@@ -9,8 +9,11 @@ import ..Microphysics0Moment
 import ..Microphysics1Moment
 import ..PrescribedSurfaceTemperature
 import ..PrognosticSurfaceTemperature
+import ..ᶜinterp
+import ..ᶠinterp
 import ..C3
 import ..C12
+import ..compute_kinetic!
 import ..PrognosticEDMFX
 import ..DiagnosticEDMFX
 import ..n_mass_flux_subdomains
@@ -30,6 +33,7 @@ import SciMLBase
 import Interpolations as Intp
 import NCDatasets as NC
 import Statistics: mean
+import ClimaUtilities.SpaceVaryingInputs
 
 include("local_state.jl")
 include("atmos_state.jl")

src/initial_conditions/initial_conditions.jl

@@ -160,6 +160,33 @@ function (initial_condition::DecayingProfile)(params)
     return local_state
 end
 
+"""
+    MoistFromFile(file_path)
+
+This function assigns an empty initial condition for , populating the `LocalState` with
+`NaN`, and later overwriting it with the content of the given file
+"""
+struct MoistFromFile <: InitialCondition
+    file_path::String
+end
+
+function (initial_condition::MoistFromFile)(params)
+    function local_state(local_geometry)
+        FT = eltype(params)
+        grav = CAP.grav(params)
+        thermo_params = CAP.thermodynamics_params(params)
+
+        T, p = FT(NaN), FT(NaN) # placeholder values
+
+        return LocalState(;
+            params,
+            geometry = local_geometry,
+            thermo_state = TD.PhaseDry_pT(thermo_params, p, T),
+        )
+    end
+    return local_state
+end
+
 """
     AgnesiHProfile(; perturb = false)
 
@@ -351,10 +378,117 @@ 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. 
+"""
+function overwrite_initial_conditions!(
+    initial_condition::InitialCondition,
+    args...,
+)
+    return nothing
+end
+
+"""
+    overwrite_initial_conditions!(initial_condition::MoistFromFile, 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`, this function
+overwrites the default initial condition and populates prognostic variables with
+interpolated values using the `SpaceVaryingInputs` tool. To mitigate issues related to
+unbalanced states following the interpolation operation, we recompute vertical pressure
+levels assuming hydrostatic balance, given the surface pressure.
+
+We expect the file to contain the following variables:
+- `p`, for pressure,
+- `t`, for temperature,
+- `q`, for humidity,
+- `u, v, w`, for velocity,
+- `cswc, crwc` for snow and rain water content (for 1 moment microphysics).
+"""
+function overwrite_initial_conditions!(
+    initial_conditions::MoistFromFile,
+    Y,
+    thermo_params,
+)
+    file_path = initial_conditions.file_path
+    isfile(file_path) || error("$(file_path) is not a file")
+    @info "Overwriting initial conditions with data from file $(file_path)"
+    center_space = Fields.axes(Y.c)
+    face_space = Fields.axes(Y.f)
+    # Using surface pressure, air temperature and specific humidity 
+    # from the dataset, compute air pressure. 
+    p_sfc = Fields.level(
+        SpaceVaryingInputs.SpaceVaryingInput(file_path, "p", face_space),
+        Fields.half,
+    )
+    ᶜ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),
+            SpaceVaryingInputs.SpaceVaryingInput(file_path, "v", center_space),
+            SpaceVaryingInputs.SpaceVaryingInput(file_path, "w", center_space),
+        )
+    Y.c.uₕ .= C12.(Geometry.UVVector.(vel))
+    Y.f.u₃ .= ᶠinterp.(C3.(Geometry.WVector.(vel)))
+    e_kin = similar(ᶜT)
+    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 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 hasproperty(Y.c, :ρq_sno) && hasproperty(Y.c, :ρq_rai)
+        Y.c.ρq_sno .=
+            SpaceVaryingInputs.SpaceVaryingInput(
+                file_path,
+                "cswc",
+                center_space,
+            ) .* Y.c.ρ
+        Y.c.ρq_rai .=
+            SpaceVaryingInputs.SpaceVaryingInput(
+                file_path,
+                "crwc",
+                center_space,
+            ) .* Y.c.ρ
+    end
+    return nothing
+end
+
 ##
 ## Baroclinic Wave
 ##
-
 function shallow_atmos_baroclinic_wave_values(z, ϕ, λ, params, perturb)
     FT = eltype(params)
     R_d = CAP.R_d(params)

src/solver/type_getters.jl

@@ -296,6 +296,8 @@ function get_initial_condition(parsed_args)
         "PrecipitatingColumn",
     ]
         return getproperty(ICs, Symbol(parsed_args["initial_condition"]))()
+    elseif isfile(parsed_args["initial_condition"])
+        return ICs.MoistFromFile(parsed_args["initial_condition"])
     elseif parsed_args["initial_condition"] == "GCM"
         @assert parsed_args["prognostic_tke"] == true
         return ICs.GCMDriven(
@@ -623,7 +625,6 @@ end
 function get_simulation(config::AtmosConfig)
     params = create_parameter_set(config)
     atmos = get_atmos(config, params)
-
     sim_info = get_sim_info(config)
     job_id = sim_info.job_id
     output_dir = sim_info.output_dir
@@ -652,7 +653,6 @@ function get_simulation(config::AtmosConfig)
 
     initial_condition = get_initial_condition(config.parsed_args)
     surface_setup = get_surface_setup(config.parsed_args)
-
     if !sim_info.restart
         s = @timed_str begin
             Y = ICs.atmos_state(
@@ -664,6 +664,17 @@ function get_simulation(config::AtmosConfig)
             t_start = Spaces.undertype(axes(Y.c))(0)
         end
         @info "Allocating Y: $s"
+
+        # In instances where we wish to interpolate existing datasets, e.g.
+        # NetCDF files containing spatially varying thermodynamic properties,
+        # this call to `overwrite_initial_conditions` overwrites the variables
+        # in `Y` (specific to `initial_condition`) with those computed using the
+        # `SpaceVaryingInputs` tool.
+        CA.InitialConditions.overwrite_initial_conditions!(
+            initial_condition,
+            Y,
+            params.thermodynamics_params,
+        )
     end
 
     tracers = get_tracers(config.parsed_args)

test/runtests.jl

@@ -24,6 +24,7 @@ using Test
 @safetestset "Topography tests" begin @time include("topography.jl") end
 @safetestset "Restarts" begin @time include("restart.jl") end
 @safetestset "Reproducibility infra" begin @time include("unit_reproducibility_infra.jl") end
+@safetestset "Init with file" begin @time include("test_init_with_file.jl") end
 
 #! format: on
 

test/test_init_with_file.jl

@@ -0,0 +1,14 @@
+import ClimaAtmos as CA
+
+simulation = CA.get_simulation(
+    CA.AtmosConfig(
+        Dict(
+            "initial_condition" => "artifact\"DYAMOND_SUMMER_ICS_p98deg\"/DYAMOND_SUMMER_ICS_p98deg.nc",
+            "moist" => "equil",
+        ),
+        job_id = "test_init_with_file_dyamond",
+    ),
+)
+
+# Just a small test to see that we got here
+@test maximum(simulation.integrator.u.c.ρ) > 0