From d40b6dbaf2d1d691e7c6aee2352d896a111c1da9 Mon Sep 17 00:00:00 2001 From: Gabriele Bozzola Date: Thu, 16 Jan 2025 17:51:48 -0800 Subject: [PATCH] Read CO2 from file --- Artifacts.toml | 9 ++++- NEWS.md | 17 +++++++-- docs/src/tracers.md | 8 ++++ src/cache/cache.jl | 9 ++++- src/cache/tracer_cache.jl | 37 ++++++++++++++++++- src/callbacks/callbacks.jl | 11 ++++++ .../radiation/radiation.jl | 4 +- src/solver/type_getters.jl | 1 + src/solver/types.jl | 22 +++++++++++ src/utils/AtmosArtifacts.jl | 11 +++++- 10 files changed, 121 insertions(+), 8 deletions(-) diff --git a/Artifacts.toml b/Artifacts.toml index 1bbe36cfbc..0945016de4 100644 --- a/Artifacts.toml +++ b/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" diff --git a/NEWS.md b/NEWS.md index 0b54101004..1dfd023d12 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,17 @@ 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 ------- @@ -11,9 +22,9 @@ v0.28.1 ### 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 diff --git a/docs/src/tracers.md b/docs/src/tracers.md index ec73782201..faa8b78f7e 100644 --- a/docs/src/tracers.md +++ b/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 ``` diff --git a/src/cache/cache.jl b/src/cache/cache.jl index abf60dd48d..b81badc008 100644 --- a/src/cache/cache.jl +++ b/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) diff --git a/src/cache/tracer_cache.jl b/src/cache/tracer_cache.jl index 064f397258..225c57a4df 100644 --- a/src/cache/tracer_cache.jl +++ b/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 diff --git a/src/callbacks/callbacks.jl b/src/callbacks/callbacks.jl index 1334d397ae..36bc693d7b 100644 --- a/src/callbacks/callbacks.jl +++ b/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) diff --git a/src/parameterized_tendencies/radiation/radiation.jl b/src/parameterized_tendencies/radiation/radiation.jl index cda8f81837..1ab6f025a8 100644 --- a/src/parameterized_tendencies/radiation/radiation.jl +++ b/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"), diff --git a/src/solver/type_getters.jl b/src/solver/type_getters.jl index 082508c161..098eb0f513 100644 --- a/src/solver/type_getters.jl +++ b/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), diff --git a/src/solver/types.jl b/src/solver/types.jl index 6059303dc0..e2c96b9bc6 100644 --- a/src/solver/types.jl +++ b/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 diff --git a/src/utils/AtmosArtifacts.jl b/src/utils/AtmosArtifacts.jl index 2a3092406c..9ae2e90ce8 100644 --- a/src/utils/AtmosArtifacts.jl +++ b/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