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

@@ -5,20 +5,24 @@ Main
 -------
 
 
-v0.28.2
--------
 ### Features
 
-### Add van Leer class operator
+### Read CO2 from file
 
-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
+`ClimaAtmos` now uses data from the Mauna Loa CO2 measurements to set CO2
+concentration. This is currently only relevant for radiation transport with
+RRTGMP.
 
-v0.28.1
+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
+
 ### Read initial conditions from NetCDF files
 
 Added functionality to allow initial conditions to be overwritten by

calibration/experiments/gcm_driven_scm/model_config_diagnostic.yml

@@ -24,7 +24,9 @@ z_max: 40e3
 z_elem: 60
 z_stretch: true
 dz_bottom: 30
-rad: allskywithclear
+rad: clearsky
+co2: fixed
+co2: fixed
 insolation: "gcmdriven"
 dt: "100secs"
 t_end: "72hours"

calibration/experiments/gcm_driven_scm/model_config_prognostic.yml

@@ -24,7 +24,9 @@ z_max: 40e3
 z_elem: 60
 z_stretch: true
 dz_bottom: 30
-rad: allskywithclear
+rad: clearsky
+co2: fixed
+co2: fixed
 insolation: "gcmdriven"
 perturb_initstate: false
 dt: "10secs"
@@ -39,4 +41,4 @@ netcdf_output_at_levels: true
 netcdf_interpolation_num_points: [2, 2, 60]
 diagnostics:
   - short_name: [hus, thetaa, clw, arup, tke, entr, detr, waup, turbentr]
-    period: 10mins
+    period: 10mins

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: "Which source to use for the concentration of prescribed CO2 to use with RRTGMP. Irrelevant when not using RRTGMP. When fixed, the volume mixing ratio is set to 397.547 parts per million. [`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

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

config/model_configs/aquaplanet_diagedmf.yml

@@ -9,7 +9,9 @@ rayleigh_sponge: true
 viscous_sponge: true
 moist: equil
 surface_setup: DefaultMoninObukhov
-rad: allskywithclear
+rad: clearsky
+co2: fixed
+co2: fixed
 insolation: "timevarying"
 dt_rad: 1hours
 dt_cloud_fraction: 1hours

config/model_configs/aquaplanet_progedmf.yml

@@ -9,7 +9,9 @@ rayleigh_sponge: true
 viscous_sponge: true
 moist: equil
 surface_setup: DefaultMoninObukhov
-rad: allskywithclear
+rad: clearsky
+co2: fixed
+co2: fixed
 insolation: "timevarying"
 dt_rad: 1hours
 dt_cloud_fraction: 1hours

config/model_configs/prognostic_edmfx_gcmdriven_column.yml

@@ -33,7 +33,9 @@ toml: [toml/prognostic_edmfx_gcmdriven.toml]
 netcdf_output_at_levels: true
 netcdf_interpolation_num_points: [2, 2, 60]
 output_default_diagnostics: false
-rad: allskywithclear
+rad: clearsky
+co2: fixed
+co2: 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

@@ -2,7 +2,9 @@ surface_setup: DefaultMoninObukhov
 surface_temperature: RCEMIPII
 insolation: rcemipii
 config: box
-rad: allskywithclear
+rad: clearsky
+co2: fixed
+co2: fixed
 turbconv: diagnostic_edmfx
 implicit_diffusion: true
 approximate_linear_solve_iters: 2

config/model_configs/rcemipii_sphere_diagnostic_edmfx.yml

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

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 `MuanaLoaCO2`. 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.MuanaLoaCO2
 ```

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,38 @@ function ozone_cache(::PrescribedOzone, Y, start_date)
     return (; o3, prescribed_o3_timevaryinginput)
 end
 
+co2_cache(_, _, _) = (;)
+function co2_cache(::FixedCO2, Y, start_date)
+    FT = Spaces.undertype(axes(Y.c))
+    # co2 is well mixed, so it is just a number, but we create an array to
+    # update it with evaluate!
+    co2 = FT[zero(FT)]
+
+    years = []
+    months = []
+    CO2_vals = []
+    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 +101,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, ::FixedCO2)
+    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(::MuanaLoaCO2) = 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_o3(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_mode"]) == "fixed"
+        return FixedCO2()
+    elseif lowercase(parsed_args["co2_mode"]) == "maunaloa"
+        return MauanaLoa()
+    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,13 @@ 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 +64,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),