Read CO2 from file

Artifacts.toml

@@ -36,7 +36,14 @@ lazy = true
     [[earth_orography_60arcseconds.download]]
     sha256 = "eca66c0701d1c2b9e271742314915ffbf4a0fae92709df611c323f38e019966e"
     url = "https://caltech.box.com/shared/static/4asrxcgl6xsgenfcug9p0wkkyhtqilgk.gz"
-
+
+[co2_dataset]
+git-tree-sha1 = "9c3bd05b68e820fceb43d130ce6b4e86ce788e4e"
+
+    [[co2_dataset.download]]
+    sha256 = "46923ec3e1f9028a11899c47e6b13d675d84daa9db2f37351a19ec9187f87865"
+    url = "https://caltech.box.com/shared/static/fuwajscgyblccy8y9aq01d0pgy91gwut.gz"
+
 [era5_cloud]
 git-tree-sha1 = "10742e0a2e343d13bb04df379e300a83402d4106"
 

NEWS.md

@@ -4,20 +4,23 @@ ClimaAtmos.jl Release Notes
 Main
 -------
 
+### Features
+
+### Read CO2 from file
+
+`ClimaAtmos` now support using data from the Mauna Loa CO2 measurements to set
+CO2 concentration. This is currently only relevant for radiation transfer with
+RRTGMP.
 
 v0.28.2
 -------
 ### Features
 
 ### Add van Leer class operator
 
-Added a new vertical transport option `vanleer_limiter` (for tracer and energy variables)
-which uses methods described in Lin et al. (1994) to apply slope-limited upwinding. Adds
-operator
-
-v0.28.1
--------
-### Features
+Added a new vertical transport option `vanleer_limiter` (for tracer and energy
+variables) which uses methods described in Lin et al. (1994) to apply
+slope-limited upwinding. Adds operator
 
 ### Read initial conditions from NetCDF files
 

calibration/experiments/gcm_driven_scm/model_config_diagnostic.yml

@@ -25,6 +25,7 @@ z_elem: 60
 z_stretch: true
 dz_bottom: 30
 rad: allskywithclear
+co2_model: fixed
 insolation: "gcmdriven"
 dt: "100secs"
 t_end: "72hours"

calibration/experiments/gcm_driven_scm/model_config_prognostic.yml

@@ -25,6 +25,7 @@ z_elem: 60
 z_stretch: true
 dz_bottom: 30
 rad: allskywithclear
+co2_model: fixed
 insolation: "gcmdriven"
 perturb_initstate: false
 dt: "10secs"

calibration/test/model_config.yml

@@ -8,6 +8,7 @@ output_dir: calibration_end_to_end_test
 output_default_diagnostics: false
 dt_rad: 6hours
 rad: clearsky
+co2_model: fixed
 diagnostics:
   - reduction_time: average
     short_name: rsut

config/default_configs/default_config.yml

@@ -289,6 +289,9 @@ restart_file:
 prescribe_ozone:
   help: "Prescribe time and spatially varying ozone from a file [`false` (default), `true`]"
   value: false
+co2_model:
+  help: "What CO2 concentration to use for RRTGMP. When fixed, it is set to 397.547 parts per million. Otherwise, it is read from the MaunaLoa measuraments. [`nothing` (default), `Fixed`, `MaunaLoa`]"
+  value: ~
 detect_restart_file:
   help: "When true, try finding a restart file and use it to restart the simulation. Only works with ActiveLink."
   value: false

config/gpu_configs/gpu_aquaplanet_dyamond_diag_1process.yml

@@ -9,6 +9,7 @@ viscous_sponge: true
 moist: equil
 precip_model: 1M
 rad: allskywithclear
+co2_model: fixed
 idealized_insolation: false
 dt_rad: 1hours
 vert_diff: "DecayWithHeightDiffusion"

config/model_configs/aquaplanet_diagedmf.yml

@@ -10,6 +10,7 @@ viscous_sponge: true
 moist: equil
 surface_setup: DefaultMoninObukhov
 rad: allskywithclear
+co2_model: maunaloa
 insolation: "timevarying"
 dt_rad: 1hours
 dt_cloud_fraction: 1hours

config/model_configs/aquaplanet_progedmf.yml

@@ -10,6 +10,7 @@ viscous_sponge: true
 moist: equil
 surface_setup: DefaultMoninObukhov
 rad: allskywithclear
+co2_model: fixed
 insolation: "timevarying"
 dt_rad: 1hours
 dt_cloud_fraction: 1hours

config/model_configs/diagnostic_edmfx_aquaplanet.yml

@@ -4,6 +4,7 @@ dz_bottom: 50.0
 rayleigh_sponge: true
 surface_setup: DefaultMoninObukhov
 rad: clearsky
+co2_model: fixed
 turbconv: diagnostic_edmfx
 implicit_diffusion: true
 approximate_linear_solve_iters: 2

config/model_configs/diagnostic_edmfx_aquaplanet_gpu.yml

@@ -3,6 +3,7 @@ z_elem: 31
 dz_bottom: 50.0
 surface_setup: DefaultMoninObukhov
 rad: clearsky
+co2_model: fixed
 turbconv: diagnostic_edmfx
 implicit_diffusion: true
 approximate_linear_solve_iters: 2

config/model_configs/prognostic_edmfx_aquaplanet.yml

@@ -5,6 +5,7 @@
 # rayleigh_sponge: true
 surface_setup: DefaultMoninObukhov
 rad: clearsky
+co2_model: fixed
 turbconv: prognostic_edmfx
 prognostic_tke: true
 edmfx_upwinding: first_order

config/model_configs/prognostic_edmfx_gcmdriven_column.yml

@@ -34,6 +34,7 @@ netcdf_output_at_levels: true
 netcdf_interpolation_num_points: [2, 2, 60]
 output_default_diagnostics: false
 rad: allskywithclear
+co2_model: fixed
 insolation: "gcmdriven"
 diagnostics:
   - short_name: [ts, ta, thetaa, ha, pfull, rhoa, ua, va, wa, hur, hus, cl, clw, cli, hussfc, evspsbl, pr]

config/model_configs/rcemipii_box_diagnostic_edmfx.yml

@@ -3,6 +3,7 @@ surface_temperature: RCEMIPII
 insolation: rcemipii
 config: box
 rad: allskywithclear
+co2_model: fixed
 turbconv: diagnostic_edmfx
 implicit_diffusion: true
 approximate_linear_solve_iters: 2

config/model_configs/rcemipii_sphere_diagnostic_edmfx.yml

@@ -1,6 +1,7 @@
 surface_setup: DefaultMoninObukhov
 surface_temperature: RCEMIPII
 rad: allskywithclear
+co2_model: fixed
 turbconv: diagnostic_edmfx
 implicit_diffusion: true
 approximate_linear_solve_iters: 2

config/perf_configs/flame_perf_target_prognostic_edmfx_aquaplanet.yml

@@ -7,6 +7,7 @@ dt_save_to_sol: "Inf"
 log_progress: false
 surface_setup: DefaultExchangeCoefficients
 rad: gray
+co2_model: fixed
 vert_diff: false
 turbconv: prognostic_edmfx 
 implicit_diffusion: true

docs/src/tracers.md

@@ -39,10 +39,21 @@ We interpolate the data from file in time every time radiation is called. The
 interpolation used is the `LinerPeriodFilling` from `ClimaUtilities`. This is a
 linear period-aware interpolation that preserves the annual cycle.
 
+### Prescribed CO2 Profile
+
+In addition to ozone, `ClimaAtmos` can prescribe CO2 concentration using data
+from [Mauna Loa CO2 measurements](https://gml.noaa.gov/ccgg/trends/data.html).
+This option is enabled with `MaunaLoaCO2`. Alternatively, `FixedCO2`
+utilizes a constant value that can be prescribed (by default, 397.547 ppm).
+
 ### More docstrings
 
 ```@docs
 ClimaAtmos.AbstractOzone
 ClimaAtmos.IdealizedOzone
 ClimaAtmos.PrescribedOzone
+
+ClimaAtmos.AbstractCO2
+ClimaAtmos.FixedCO2
+ClimaAtmos.MaunaLoaCO2
 ```

reproducibility_tests/ref_counter.jl

@@ -1,4 +1,4 @@
-201
+202
 
 # **README**
 #
@@ -20,6 +20,9 @@
 
 
 #=
+202
+- Slightly changed CO2 prescription
+
 201
 - Updated ClimaTimeSteppers to 0.7.39, slightly improving conservation properties
 

src/cache/cache.jl

@@ -147,7 +147,14 @@ function build_cache(Y, atmos, params, surface_setup, sim_info, aerosol_names)
 
     radiation_args =
         atmos.radiation_mode isa RRTMGPI.AbstractRRTMGPMode ?
-        (start_date, params, atmos.ozone, aerosol_names, atmos.insolation) : ()
+        (
+            start_date,
+            params,
+            atmos.ozone,
+            atmos.co2,
+            aerosol_names,
+            atmos.insolation,
+        ) : ()
 
     hyperdiff = hyperdiffusion_cache(Y, atmos)
     precipitation = precipitation_cache(Y, atmos)

src/cache/tracer_cache.jl

@@ -1,4 +1,5 @@
 import Dates: Year
+import ClimaUtilities
 import ClimaUtilities.TimeVaryingInputs
 import ClimaUtilities.TimeVaryingInputs:
     TimeVaryingInput, LinearPeriodFillingInterpolation
@@ -20,6 +21,43 @@ function ozone_cache(::PrescribedOzone, Y, start_date)
     return (; o3, prescribed_o3_timevaryinginput)
 end
 
+co2_cache(_, _, _) = (;)
+function co2_cache(::MaunaLoaCO2, Y, start_date)
+    FT = Spaces.undertype(axes(Y.c))
+    # ClimaUtilities < v0.1.21 can only write to Arrays that are on the same
+    # device as the space
+    ArrayType =
+        pkgversion(ClimaUtilities) < v"0.1.21" ? ClimaComms.array_type(Y.c) :
+        Array
+    # co2 is well mixed, so it is just a number, but we create a mutable object
+    # to update it with `evaluate!`
+    co2 = ArrayType([zero(FT)])
+
+    years = Int[]
+    months = Int[]
+    CO2_vals = FT[]
+    open(
+        AA.co2_concentration_file_path(; context = ClimaComms.context(Y.c)),
+        "r",
+    ) do file
+        for line in eachline(file)
+            # Skip comments
+            startswith(line, '#') && continue
+            parts = split(line)
+            push!(years, parse(Int, parts[1]))
+            push!(months, parse(Int, parts[2]))
+            # convert from ppm to fraction, data is in fourth column of the text file
+            push!(CO2_vals, parse(Float64, parts[4]) / 1_000_000)
+        end
+    end
+    # The text file only has month and year, so we set the day to 15th of the month
+    CO2_dates = Dates.DateTime.(years, months, 15)
+    CO2_times =
+        ClimaUtilities.Utils.period_to_seconds_float.(CO2_dates .- start_date)
+    prescribed_co2_timevaryinginput = TimeVaryingInput(CO2_times, CO2_vals)
+    return (; co2, prescribed_co2_timevaryinginput)
+end
+
 function tracer_cache(Y, atmos, prescribed_aerosol_names, start_date)
     if !isempty(prescribed_aerosol_names)
         target_space = axes(Y.c)
@@ -68,5 +106,6 @@ function tracer_cache(Y, atmos, prescribed_aerosol_names, start_date)
         aerosol_cache = (;)
     end
     o3_cache = ozone_cache(atmos.ozone, Y, start_date)
-    return (; aerosol_cache..., o3_cache...)
+    co2_cache_nt = co2_cache(atmos.co2, Y, start_date)
+    return (; aerosol_cache..., o3_cache..., co2_cache_nt...)
 end

src/callbacks/callbacks.jl

@@ -59,6 +59,13 @@ function update_o3!(p, t, ::PrescribedOzone)
     return nothing
 end
 
+update_co2!(_, _, _) = nothing
+function update_co2!(p, t, ::MaunaLoaCO2)
+    evaluate!(p.tracers.co2, p.tracers.prescribed_co2_timevaryinginput, t)
+    return nothing
+end
+
+
 NVTX.@annotate function rrtmgp_model_callback!(integrator)
     Y = integrator.u
     p = integrator.p
@@ -71,6 +78,7 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
 
     # If we have prescribed ozone or aerosols, we need to update them
     update_o3!(p, t, p.atmos.ozone)
+    update_co2!(p, t, p.atmos.co2)
     if :prescribed_aerosols_field in propertynames(p.tracers)
         for (key, tv) in pairs(p.tracers.prescribed_aerosol_timevaryinginputs)
             field = getproperty(p.tracers.prescribed_aerosols_field, key)
@@ -255,6 +263,9 @@ NVTX.@annotate function rrtmgp_model_callback!(integrator)
             )
             @. ᶜvmr_o3 = p.tracers.o3
         end
+        if :co2 in propertynames(p.tracers)
+            @. rrtmgp_model.volume_mixing_ratio_co2 = p.tracers.co2
+        end
     end
 
     set_surface_albedo!(Y, p, t, p.atmos.surface_albedo)

src/parameterized_tendencies/radiation/radiation.jl

@@ -86,12 +86,24 @@ function idealized_ozone(z::FT) where {FT}
     return g1 * p^g2 * exp(-p / g3) * PPMV_TO_VMR
 end
 
+#######
+# CO2 #
+#######
+
+function center_vmr_co2(co2::FixedCO2)
+    return co2.value
+end
+
+# Initialized in callback
+center_vmr_co2(::MaunaLoaCO2) = NaN
+
 function radiation_model_cache(
     Y,
     radiation_mode::RRTMGPI.AbstractRRTMGPMode,
     start_date,
     params,
     ozone,
+    co2,
     aerosol_names,
     insolation_mode;
     interpolation = RRTMGPI.BestFit(),
@@ -150,6 +162,10 @@ function radiation_model_cache(
         else
             center_volume_mixing_ratio_o3 = center_vmr_o3(ozone, Y)
 
+            # FT is needed in case FixedCO2 is being used with an inconsistent
+            # floating point type
+            center_volume_mixing_ratio_co2 = FT(center_vmr_co2(co2))
+
             # the first value for each global mean volume mixing ratio is the
             # present-day value
             input_vmr(name) =
@@ -160,7 +176,7 @@ function radiation_model_cache(
                 center_volume_mixing_ratio_h2o = NaN, # initialize in tendency
                 center_relative_humidity = NaN, # initialized in callback
                 center_volume_mixing_ratio_o3,
-                volume_mixing_ratio_co2 = input_vmr("carbon_dioxide_GM"),
+                volume_mixing_ratio_co2 = center_volume_mixing_ratio_co2,
                 volume_mixing_ratio_n2o = input_vmr("nitrous_oxide_GM"),
                 volume_mixing_ratio_co = input_vmr("carbon_monoxide_GM"),
                 volume_mixing_ratio_ch4 = input_vmr("methane_GM"),

src/solver/model_getters.jl

@@ -308,6 +308,18 @@ function get_ozone(parsed_args)
     return parsed_args["prescribe_ozone"] ? PrescribedOzone() : IdealizedOzone()
 end
 
+function get_co2(parsed_args)
+    if isnothing(parsed_args["co2_model"])
+        return nothing
+    elseif lowercase(parsed_args["co2_model"]) == "fixed"
+        return FixedCO2()
+    elseif lowercase(parsed_args["co2_model"]) == "maunaloa"
+        return MaunaLoaCO2()
+    else
+        error("The CO2 models supported are $(subtypes(AbstractCO2))")
+    end
+end
+
 function get_cloud_in_radiation(parsed_args)
     isnothing(parsed_args["prescribe_clouds_in_radiation"]) && return nothing
     return parsed_args["prescribe_clouds_in_radiation"] ?

src/solver/type_getters.jl

@@ -30,6 +30,16 @@ function get_atmos(config::AtmosConfig, params)
         @warn "prescribe_ozone is set to nothing with an RRTMGP model. Resetting to IdealizedOzone. This behavior will stop being supported in some future release"
         ozone = IdealizedOzone()
     end
+    co2 = get_co2(parsed_args)
+    with_rrtgmp = radiation_mode isa RRTMGPI.AbstractRRTMGPMode
+    if with_rrtgmp && isnothing(co2)
+        @warn (
+            "co2_model set to nothing with an RRTGMP model. Resetting to FixedCO2"
+        )
+        co2 = FixedCO2()
+    end
+    (isnothing(co2) && !with_rrtgmp) &&
+        @warn ("$(co2) does nothing if RRTGMP is not used")
 
     diffuse_momentum = !(forcing_type isa HeldSuarezForcing)
 
@@ -57,6 +67,7 @@ function get_atmos(config::AtmosConfig, params)
     atmos = AtmosModel(;
         moisture_model,
         ozone,
+        co2,
         radiation_mode,
         subsidence = get_subsidence_model(parsed_args, radiation_mode, FT),
         ls_adv = get_large_scale_advection_model(parsed_args, FT),