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,16 +4,27 @@ ClimaAtmos.jl Release Notes
 Main
 -------
 
+v0.28.2
+-------
+
+### Features
+
+### Read CO2 from file
+
+`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
 -------
 ### 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 
+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.
 
 ### Read initial conditions from NetCDF files
 

docs/src/tracers.md

@@ -39,10 +39,18 @@ 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).
+
 ### More docstrings
 
 ```@docs
 ClimaAtmos.AbstractOzone
 ClimaAtmos.IdealizedOzone
 ClimaAtmos.PrescribedOzone
+
+ClimaAtmos.AbstractCO2
+ClimaAtmos.PrescribedCO2
 ```

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,39 @@ function ozone_cache(::PrescribedOzone, Y, start_date)
     return (; o3, prescribed_o3_timevaryinginput)
 end
 
+co2_cache(_, _, _) = (;)
+function co2_cache(::PrescribedCO2, 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 evalaute!
+    co2 = [zero(FT)]
+    extrapolation_bc = (Intp.Periodic(), Intp.Flat(), Intp.Flat())
+
+    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 +102,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, ::PrescribedCO2)
+    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,14 @@ function idealized_ozone(z::FT) where {FT}
     return g1 * p^g2 * exp(-p / g3) * PPMV_TO_VMR
 end
 
+
 function radiation_model_cache(
     Y,
     radiation_mode::RRTMGPI.AbstractRRTMGPMode,
     start_date,
     params,
     ozone,
+    co2,
     aerosol_names,
     insolation_mode;
     interpolation = RRTMGPI.BestFit(),
@@ -160,7 +162,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 = NaN, # initialized in callback
                 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/type_getters.jl

@@ -57,6 +57,7 @@ function get_atmos(config::AtmosConfig, params)
     atmos = AtmosModel(;
         moisture_model,
         ozone,
+        co2 = PrescribedCO2(),
         radiation_mode,
         subsidence = get_subsidence_model(parsed_args, radiation_mode, FT),
         ls_adv = get_large_scale_advection_model(parsed_args, FT),

src/solver/types.jl

@@ -127,6 +127,23 @@ Refer to ClimaArtifacts for more information on how to obtain the artifact.
 """
 struct PrescribedOzone <: AbstractOzone end
 
+"""
+    AbstractCO2
+
+Describe how CO2 concentration should be set.
+"""
+abstract type AbstractCO2 end
+
+"""
+    PrescribedOzone
+
+Implement a time-varying ozone profile as read from disk.
+
+The data from the Mauna Loa CO2 measurements is used. It is a assumed that the
+concentration is constant.
+"""
+struct PrescribedCO2 <: AbstractCO2 end
+
 """
     AbstractCloudInRadiation
 
@@ -456,6 +473,7 @@ Base.@kwdef struct AtmosModel{
     F,
     S,
     OZ,
+    CO2,
     RM,
     LA,
     EXTFORCING,
@@ -486,8 +504,12 @@ Base.@kwdef struct AtmosModel{
     forcing_type::F = nothing
     subsidence::S = nothing
 
+    # Currently only relevant for RRTGMP, but will hopefully become standalone
+    # in the future
     """What to do with ozone for radiation (when using RRTGMP)"""
     ozone::OZ = nothing
+    """What to do with co2 for radiation (when using RRTGMP)"""
+    co2::CO2 = nothing
 
     radiation_mode::RM = nothing
     ls_adv::LA = nothing

src/utils/AtmosArtifacts.jl

@@ -72,7 +72,7 @@ end
 
 Construct the file path for the 60arcsecond orography data NetCDF file.
 
-Downloads the 60arc-second dataset by default. 
+Downloads the 60arc-second dataset by default.
 """
 function earth_orography_file_path(; context = nothing)
     filename = "ETOPO_2022_v1_60s_N90W180_surface.nc"
@@ -82,4 +82,13 @@ function earth_orography_file_path(; context = nothing)
     )
 end
 
+"""
+    co2_concentration_file_path(; context = nothing)
+
+Construct the file path for the co2 concentration CSV file.
+"""
+function co2_concentration_file_path(; context = nothing)
+    return joinpath(@clima_artifact("co2_dataset", context), "co2_mm_mlo.txt")
+end
+
 end