From 27f4413dc88bd4ca6009d5fcccae8217359245a6 Mon Sep 17 00:00:00 2001 From: Kevin Phan <98072684+ph-kev@users.noreply.github.com> Date: Thu, 7 Nov 2024 14:10:30 -0800 Subject: [PATCH] Add reasonable conversion from z to pressure To convert from z to pressure, we use p(z) = P0 * exp(-z / H_EARTH), where P0 = 10000 and H_EARTH = 7000.0. --- NEWS.md | 4 +++ src/Atmos.jl | 22 +++++++-------- test/test_Atmos.jl | 68 +++++++++++++++++++++++++++++++++++++++++----- 3 files changed, 75 insertions(+), 19 deletions(-) diff --git a/NEWS.md b/NEWS.md index 6a788b3c..8583782c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -28,6 +28,10 @@ julia> ClimaAnalysis.global_rmse_pfull(sim_var, obs_var, sim_pressure = pressure ## Bug fixes - `Atmos.to_pressure_coordinates` now works with Unitful units. +- `Atmos.to_pressure_coordinates` now uses reasonable pressure values when `target_pressure` + is not specified. In particular, the vertical dimension is mapped to pressure levels by z + -> P0 * exp(-z / H_EARTH), where P0 = 10000 and H_EARTH = 7000.0, following a simple + hydrostatic model for the atmosphere. v0.5.11 ------- diff --git a/src/Atmos.jl b/src/Atmos.jl index 5d25e63c..65c4e70d 100644 --- a/src/Atmos.jl +++ b/src/Atmos.jl @@ -43,8 +43,8 @@ end Change the vertical dimension of `var` to be in pressure coordinates. If `target_pressure` is nothing, the target pressure levels are computed by -linearly sampling the interval `minimum(pressure), maximum(pressure)`. Then, for -each column in `var`, the values are linearly interpolate onto this new grid. +Plvl(z) = P0 * exp(-z / H_EARTH), where H_EARTH = 7000.0 and P0 = 1e5, following +a simple hydrostatic model for the atmosphere. `target_pressure` can be set to a `Vector` to specify custom pressure levels. @@ -66,18 +66,16 @@ function to_pressure_coordinates( z_index = var.dim2index[z_name] pressure_name = short_name(pressure) - # First, we construct the target pressure grid. For this, we take the - # extrema of pressure and divide the interval linearly with the same number - # of points we originally had in z + # First, we construct the target pressure grid. For this, we use the mapping + # Plvl(z) = P0 * exp(-z / H_EARTH), where H_EARTH = 7000.0 and P0 = 1e5. if isnothing(target_pressure) - # TODO: Pick this more sensibly - # TODO: Make it go from max to min? (This is not supported by Interpolations.jl...) - target_pressure = range( - minimum(pressure.data), - maximum(pressure.data), - length = length(var.dims[z_name]), - ) + H_EARTH = 7000.0 + P0 = 1e5 + Plvl(z) = P0 * exp(-z / H_EARTH) + + # Reverse vector because Interpolations.jl require increasing knots + target_pressure = reverse(Plvl.(var.dims[z_name])) end # Then, we prepare the output variable diff --git a/test/test_Atmos.jl b/test/test_Atmos.jl index 480e92b7..ed0fb91a 100644 --- a/test/test_Atmos.jl +++ b/test/test_Atmos.jl @@ -3,6 +3,12 @@ import ClimaAnalysis import OrderedCollections: OrderedDict @testset "To pressure coordinates" begin + # Define function to make it easier to test by getting equispaced pressure levels + find_pressure_range(var, pressure) = range( + minimum(pressure.data), + maximum(pressure.data), + length = length(var.dims["z"]), + ) # Let start with testing a single column z_alt = 0:100.0 |> collect @@ -20,7 +26,11 @@ import OrderedCollections: OrderedDict ClimaAnalysis.OutputVar(attribs, Dict("z" => z_alt), dim_attribs, pdata) pressure_in_pressure_coordinates = - ClimaAnalysis.Atmos.to_pressure_coordinates(pressure_var, pressure_var) + ClimaAnalysis.Atmos.to_pressure_coordinates( + pressure_var, + pressure_var, + target_pressure = find_pressure_range(pressure_var, pressure_var), + ) @test collect(keys(pressure_in_pressure_coordinates.dims)) == ["pfull"] # reverse because we go from min to max for pressure (to have it increasing @@ -44,8 +54,11 @@ import OrderedCollections: OrderedDict mydata, ) - myvar_in_pressure_coordinates = - ClimaAnalysis.Atmos.to_pressure_coordinates(myvar, pressure_var) + myvar_in_pressure_coordinates = ClimaAnalysis.Atmos.to_pressure_coordinates( + myvar, + pressure_var, + target_pressure = find_pressure_range(myvar, pressure_var), + ) @test collect(keys(myvar_in_pressure_coordinates.dims)) == ["pfull"] @test myvar_in_pressure_coordinates.dims["pfull"] == reverse(pdata) @@ -64,7 +77,11 @@ import OrderedCollections: OrderedDict ) myvar_in_exp_pressure_coordinates = - ClimaAnalysis.Atmos.to_pressure_coordinates(myvar, exp_pressure_var) + ClimaAnalysis.Atmos.to_pressure_coordinates( + myvar, + exp_pressure_var, + target_pressure = find_pressure_range(myvar, exp_pressure_var), + ) # Linear range from min to max expected_range = collect( @@ -113,7 +130,11 @@ import OrderedCollections: OrderedDict pressure3D = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, pdata3D) my3Dvar_in_exp_pressure_coordinates = - ClimaAnalysis.Atmos.to_pressure_coordinates(var3D, pressure3D) + ClimaAnalysis.Atmos.to_pressure_coordinates( + var3D, + pressure3D, + target_pressure = find_pressure_range(var3D, pressure3D), + ) # Linear range from min to max overall_range = @@ -127,6 +148,28 @@ import OrderedCollections: OrderedDict ] @test my3Dvar_in_exp_pressure_coordinates.data ≈ expected_output rtol = 1e-5 + # Test z to pressure level conversion + z_alt = 0:100.0 |> collect + data = copy(z_alt) + + zvar = ClimaAnalysis.OutputVar(Dict("z" => z_alt), data) + + # Fake pressure, linearly decreasing, so that we can check precisely + pressure = 300.0:-2.0:100.0 |> collect + pdata = copy(pressure) + + attribs = Dict("short_name" => "pfull") + dim_attribs = Dict{String, Any}() + pressure_var = + ClimaAnalysis.OutputVar(attribs, Dict("z" => z_alt), dim_attribs, pdata) + + pressure_in_pressure_coordinates = + ClimaAnalysis.Atmos.to_pressure_coordinates(pressure_var, pressure_var) + H_EARTH = 7000.0 + P0 = 1e5 + Plvl(z) = P0 * exp(-z / H_EARTH) + @test reverse(Plvl.(z_alt)) ≈ pressure_in_pressure_coordinates.dims["pfull"] + # Error checking @test_throws ErrorException ClimaAnalysis.Atmos.to_pressure_coordinates( var3D, @@ -141,6 +184,12 @@ import OrderedCollections: OrderedDict end @testset "RMSE of pressure coordinates" begin + find_pressure_range(var, pressure) = range( + minimum(pressure.data), + maximum(pressure.data), + length = length(var.dims["z"]), + ) + lon = 0.5:1.0:359.5 |> collect lat = -89.5:1.0:89.5 |> collect zzz = 1.0:10.0 |> collect @@ -180,13 +229,17 @@ end zero_var3D = ClimaAnalysis.OutputVar(attribs, dims, dim_attribs, zero_data) pressure3D = ClimaAnalysis.OutputVar(pfull_attribs, dims, pfull_dim_attribs, pdata3D) - zero_var = - ClimaAnalysis.Atmos.to_pressure_coordinates(zero_var3D, pressure3D) + zero_var = ClimaAnalysis.Atmos.to_pressure_coordinates( + zero_var3D, + pressure3D, + target_pressure = find_pressure_range(zero_var3D, pressure3D), + ) global_rmse_pfull = ClimaAnalysis.Atmos.global_rmse_pfull( var3D, zero_var, sim_pressure = pressure3D, + target_pressure = find_pressure_range(var3D, pressure3D), ) min_pfull, max_pfull = extrema(zero_var.dims["pfull"]) @@ -217,6 +270,7 @@ end var3D, pfull_obs_var, sim_pressure = pressure3D, + target_pressure = find_pressure_range(var3D, pressure3D), ) @test global_rmse_pfull >= 0.0