Merge pull request #4012 from CliMA/j/era5_warming_scm

Adds ability to run ERA5 simulations with temperature change

Author: Julians42

config/default_configs/default_config.yml

@@ -271,6 +271,9 @@ site_latitude:
 site_longitude:
   help: "Site longitude for single column model. Used for externally driven time varying forcing model to generate the forcing file. Artifact support is currently for eastern Pacific region in July 2007 only. [`-149.0` (default)]"
   value: -149.0
+era5_diurnal_warming:
+  help: Applies a warming scenario by adding a constant temperature offset (in Kelvin). Supported only by the `ReanalysisMonthlyAveragedDiurnal` forcing model. This affects the target surface temperature for the entire simulation, the initial atmospheric temperature profile, and the relaxation temperature profile above `gcmdriven_relaxation_minimum_height`. Specific humidity is increased to hold relative humidity constant in the initial and relaxation profiles.
+  value: ~
 subsidence:
   help: "Subsidence [`nothing` (default), `Bomex`, `LifeCycleTan2018`, `Rico`, `DYCOMS`, `ISDAC`]"
   value: ~

src/solver/model_getters.jl

@@ -466,6 +466,10 @@ function get_external_forcing_model(parsed_args, ::Type{FT}) where {FT}
         @assert parsed_args["config"] == "column" "ReanalysisTimeVarying and ReanalysisMonthlyAveragedDiurnal are only supported in column mode."
         @assert all(reanalysis_required_fields .== "ReanalysisTimeVarying") "All of external_forcing, surface_setup, surface_temperature and initial_condition must be set to ReanalysisTimeVarying."
     end
+    if !isnothing(parsed_args["era5_diurnal_warming"])
+        @assert external_forcing == "ReanalysisMonthlyAveragedDiurnal" "era5_diurnal_warming is only supported for ReanalysisMonthlyAveragedDiurnal."
+        @assert parsed_args["era5_diurnal_warming"] isa Number "era5_diurnal_warming is expected to be a number, but was supplied as a $(typeof(parsed_args["era5_diurnal_warming"]))"
+    end
     return if isnothing(external_forcing)
         nothing
     elseif external_forcing == "GCM"

src/utils/era5_observations_to_forcing_file.jl

@@ -69,6 +69,10 @@ function get_external_monthly_forcing_file_path(
     ),
 )
     start_date = parsed_args["start_date"]
+    warming_amount = parsed_args["era5_diurnal_warming"]
+    warming_amount_str =
+        (warming_amount isa Number) ? "_plus_$(float(warming_amount))K" : ""
+
     # round to era5 quarter degree resolution for site selection
     site_latitude = round(parsed_args["site_latitude"] * 4) / 4
     site_longitude = round(parsed_args["site_longitude"] * 4) / 4
@@ -79,7 +83,7 @@ function get_external_monthly_forcing_file_path(
     end
     return joinpath(
         data_dir,
-        "monthly_diurnal_cycle_forcing_$(site_latitude)_$(site_longitude)_$(start_date).nc",
+        "monthly_diurnal_cycle_forcing_$(site_latitude)_$(site_longitude)_$(start_date)$(warming_amount_str).nc",
     )
 end
 
@@ -340,6 +344,40 @@ function generate_external_forcing_file(
         )
     end
 
+    warming_amount = parsed_args["era5_diurnal_warming"]
+    if warming_amount isa Number
+        # Get the relative humidity before warming; this will be used to scale the specific humidity after warming
+        param_set = TD.Parameters.ThermodynamicsParameters(FT)
+        phase_partition =
+            TD.PhasePartition.(sim_forcing["hus"], sim_forcing["clw"], sim_forcing["cli"])
+        # expand dimension and convert pressure levels from hPa (ERA5 raw) to Pa (needed for thermodynamics)
+        p_expanded =
+            FT.(
+                repeat(
+                    sim_forcing["pressure_level"] .* 100,
+                    1,
+                    size(sim_forcing["hus"], 2),
+                )
+            )
+        relative_humidity =
+            TD.relative_humidity.(
+                param_set,
+                sim_forcing["ta"],
+                p_expanded,
+                TD.PhaseEquil{FT},
+                phase_partition,
+            )
+
+        # Warming the temperature
+        sim_forcing["ta"] .+= warming_amount
+        ρ = TD.air_density.(param_set, sim_forcing["ta"], p_expanded, phase_partition)
+
+        # Get the saturation specific humidity at the warmed temperature
+        saturation_humidity_at_warming =
+            TD.q_vap_saturation.(param_set, sim_forcing["ta"], ρ, TD.PhaseEquil{FT})
+        sim_forcing["hus"] = saturation_humidity_at_warming .* relative_humidity
+    end
+
     sim_forcing["z"] =
         sim_forcing["zg"] / external_tv_params.gravitational_acceleration # height in meters
 
@@ -501,6 +539,11 @@ function generate_external_forcing_file(
         smooth_amount = smooth_amount,
     )
 
+    # warm the surface temperature if warming amount is specified
+    if warming_amount isa Number
+        skt .+= warming_amount
+    end
+
     for time_ind in 1:ds.dim["time"]
         ds["coszen"][:, :, :, time_ind] .= coszen_list[time_ind][1]
         ds["rsdt"][:, :, :, time_ind] .= coszen_list[time_ind][2]
@@ -557,6 +600,7 @@ function generate_multiday_era5_external_forcing_file(
             "start_date" => Dates.format(dd, "yyyymmdd"),
             "site_latitude" => parsed_args["site_latitude"],
             "site_longitude" => parsed_args["site_longitude"],
+            "era5_diurnal_warming" => parsed_args["era5_diurnal_warming"],
         )
         single_file_path = get_external_daily_forcing_file_path(
             single_parsed_args;

test/test_helpers.jl

@@ -180,6 +180,7 @@ function create_mock_era5_datasets(
     num_pressure = 37,
     base_date = "20000101",
 )
+    params = CA.ClimaAtmosParameters(FT)
     column_data_path = joinpath(
         temporary_dir,
         "forcing_and_cloud_hourly_profiles_$(start_date).nc",
@@ -214,7 +215,7 @@ function create_mock_era5_datasets(
     tvforcing["latitude"][:] = collect(-2.0:lat_step:2.0)
     tvforcing["longitude"][:] = collect(-2.0:lon_step:2.0)
     tvforcing["pressure_level"][:] =
-        10 .^ (range(1, stop = 4, length = num_pressure))
+        10 .^ (range(1, stop = 3, length = num_pressure))
     # Time values are hours since base_date, not since this file's date
     tvforcing["valid_time"][:] =
         collect(hours_offset:(hours_offset + hours - 1))
@@ -234,11 +235,19 @@ function create_mock_era5_datasets(
     tvforcing["u"][:, :, :, :] .= ones(FT, size(tvforcing["u"]))
     tvforcing["v"][:, :, :, :] .= ones(FT, size(tvforcing["v"]))
     tvforcing["w"][:, :, :, :] .= ones(FT, size(tvforcing["w"]))
-    tvforcing["t"][:, :, :, :] .= ones(FT, size(tvforcing["t"]))
-    tvforcing["q"][:, :, :, :] .= ones(FT, size(tvforcing["q"]))
-    tvforcing["z"][:, :, :, :] .= ones(FT, size(tvforcing["z"]))
-    tvforcing["clwc"][:, :, :, :] .= ones(FT, size(tvforcing["clwc"]))
-    tvforcing["ciwc"][:, :, :, :] .= ones(FT, size(tvforcing["ciwc"]))
+    tvforcing["clwc"][:, :, :, :] .= zeros(FT, size(tvforcing["clwc"]))
+    tvforcing["ciwc"][:, :, :, :] .= zeros(FT, size(tvforcing["ciwc"]))
+
+    tvforcing["z"] .= reshape(
+        tvforcing["pressure_level"] .* 100,
+        (1, 1, length(tvforcing["pressure_level"]), 1),
+    )
+    tvforcing["z"] .=
+        geopotential_from_pressure(tvforcing["z"]; R_d = CA.Parameters.R_d(params))
+    tvforcing["t"][:, :, :, :] .=
+        temperature_from_geopotential(tvforcing["z"]; g = CA.Parameters.grav(params))
+    tvforcing["q"][:, :, :, :] .=
+        shum_from_geopotential(tvforcing["z"]; g = CA.Parameters.grav(params))
     close(tvforcing)
 
     # Create accumulated dataset (3D: lon, lat, time)
@@ -287,3 +296,65 @@ function create_mock_era5_datasets(
 
     return column_data_path, accum_data_path, inst_data_path
 end
+
+"""
+    geopotential_from_pressure(pressure_levels; R_d = 287.05, sp = 101325, T_avg = 250)
+
+Convert pressure levels to geopotential levels using the hypsometric equation.
+R_d is the gas constant for dry air, sp is the surface pressure in Pa, and T_avg is the average temperature in the troposphere in K.
+"""
+function geopotential_from_pressure(pressure_levels; R_d = 287.05, sp = 101325, T_avg = 250)
+    return R_d .* T_avg .* log.(sp ./ pressure_levels)
+end
+
+"""
+    temperature_from_geopotential(geopotential_levels; Γ = 0.0065, T_surf = 300, T_min = 200, g = 9.81)
+
+Produce a typical temperature profile from geopotential levels. 
+Γ is the lapse rate in K/m, T_surf is the surface temperature in K, and T_min is the minimum temperature in K, which simulates the tropopause.
+"""
+function temperature_from_geopotential(
+    geopotential_levels;
+    Γ = 0.0065,
+    T_surf = 300,
+    T_min = 200,
+    g = 9.81,
+)
+    return max.(T_min, T_surf .- Γ .* geopotential_levels ./ g)
+end
+
+"""
+    shum_from_geopotential(geopotential_levels; q0 = 0.02, shum_scale_height = 2e3, g = 9.81)
+
+Produce a typical specific humidity profile from geopotential levels.
+q0 is the specific humidity at the surface and shum_scale_height is the scale height in geopotential coordinates.
+The default values are a 2% mixing ratio at the surface and a 2km scale height.
+"""
+function shum_from_geopotential(
+    geopotential_levels;
+    q0 = 0.02,
+    shum_scale_height = 2e3,
+    g = 9.81,
+)
+    return q0 .* exp.(.-geopotential_levels ./ (shum_scale_height * g))
+end
+
+"""
+    monotonic_decreasing(A, dim)
+
+Check if an array is monotonically decreasing along a given dimension.
+"""
+function monotonic_decreasing(A, dim)
+    all_diffs = mapslices(x -> diff(x), A; dims = dim)
+    return all(all_diffs .<= 0)
+end
+
+"""
+    monotonic_increasing(A, dim)
+
+Check if an array is monotonically increasing along a given dimension.
+"""
+function monotonic_increasing(A, dim)
+    all_diffs = mapslices(x -> diff(x), A, dims = dim)
+    return all(all_diffs .>= 0)
+end

test/utilities.jl

@@ -309,6 +309,7 @@ end
         "site_latitude" => 0.0,
         "site_longitude" => 0.0,
         "t_end" => "5hours",
+        "era5_diurnal_warming" => Nothing,
     )
 
     temporary_dir = mktempdir()
@@ -345,10 +346,21 @@ end
     processed_data = NCDataset(sim_forcing_daily, "r")
 
     # Test fixed variables - this tests that the variables are copied correctly
-    for clima_var in ["hus", "ta", "ua", "va", "wap", "zg", "clw", "cli", "ts"]
+    for clima_var in ["ua", "va", "wap", "ts"]
         @test all(isapprox.(processed_data[clima_var][:], 1, atol = 1e-10))
     end
 
+    for clima_var in ["clw", "cli"]
+        @test all(isapprox.(processed_data[clima_var][:], 0, atol = 1e-10))
+    end
+
+    # data is stored from top of atmosphere to surface 
+    @test monotonic_decreasing(processed_data["zg"], 3)
+    @test monotonic_increasing(processed_data["ta"], 3)
+    @test monotonic_increasing(processed_data["hus"], 3)
+    @test all(processed_data["hus"] .>= 0)
+    @test all(processed_data["ta"] .>= 200) # 200K is the minimum temperature set in the helper function
+
     # Test accumulated variables - note that the sign is flipped because of differences between ecmwf and clima
     for clima_var in ["hfls", "hfss"]
         @test all(
@@ -387,6 +399,7 @@ end
         "site_latitude" => 0.0,
         "site_longitude" => 0.0,
         "t_end" => "2days",
+        "era5_diurnal_warming" => Nothing,
     )
 
     input_dir = mktempdir()
@@ -436,7 +449,7 @@ end
     @test length(processed_data["time"][:]) == expected_time_steps
 
     # Test that data is consistent across time
-    @test all(x -> all(isapprox.(x, 1, atol = 1e-10)), processed_data["ta"][:])
+    @test monotonic_increasing(processed_data["ta"], 3)
     @test all(
         x -> all(isapprox.(x, -1 / time_resolution, atol = 1e-10)),
         processed_data["hfls"][:],
@@ -457,6 +470,7 @@ end
     # compute the vertical temperature gradient
     vertical_temperature_gradient =
         CA.get_vertical_tendencies(vert_partial_ds, "ta")
+    @test all(vertical_temperature_gradient .<= 0)
 
     close(processed_data)
 end
@@ -558,3 +572,67 @@ end
         close(test_ds)
     end
 end
+
+@testset "ERA5 diurnal warming" begin
+    FT = Float32
+    temporary_dir = mktempdir()
+
+    parsed_args_0K = Dict(
+        "start_date" => "20000506",
+        "site_latitude" => 0.0,
+        "site_longitude" => 0.0,
+        "t_end" => "5hours",
+        "era5_diurnal_warming" => Nothing,
+    )
+
+    parsed_args_4K = Dict(
+        "start_date" => "20000506",
+        "site_latitude" => 0.0,
+        "site_longitude" => 0.0,
+        "t_end" => "5hours",
+        "era5_diurnal_warming" => 4,
+    )
+
+    create_mock_era5_datasets(temporary_dir, parsed_args_0K["start_date"], FT)
+
+    sim_forcing_0K = CA.get_external_monthly_forcing_file_path(
+        parsed_args_0K,
+        data_dir = temporary_dir,
+    )
+    sim_forcing_4K = CA.get_external_monthly_forcing_file_path(
+        parsed_args_4K,
+        data_dir = temporary_dir,
+    )
+
+    @test basename(sim_forcing_0K) == "monthly_diurnal_cycle_forcing_0.0_0.0_20000506.nc"
+    @test basename(sim_forcing_4K) ==
+          "monthly_diurnal_cycle_forcing_0.0_0.0_20000506_plus_4.0K.nc"
+
+    CA.generate_external_forcing_file(
+        parsed_args_0K,
+        sim_forcing_0K,
+        FT,
+        input_data_dir = temporary_dir,
+    )
+
+    CA.generate_external_forcing_file(
+        parsed_args_4K,
+        sim_forcing_4K,
+        FT,
+        input_data_dir = temporary_dir,
+    )
+
+    # open the datasets and check the temperature and specific humidity profiles have been adjusted
+    processed_data_0K = NCDataset(sim_forcing_0K, "r")
+    processed_data_4K = NCDataset(sim_forcing_4K, "r")
+
+    # check that air and surface temperatures have increased by the +4K amount
+    @test all(isapprox.(processed_data_0K["ta"] .+ 4, processed_data_4K["ta"]))
+    @test isapprox(processed_data_0K["ts"] .+ 4, processed_data_4K["ts"])
+
+    # check that the specific humidity has increased
+    @test all(processed_data_0K["hus"] .< processed_data_4K["hus"])
+
+    close(processed_data_0K)
+    close(processed_data_4K)
+end