From dfc391a20dc13fbf70e0fcd61480a70320bb1ff7 Mon Sep 17 00:00:00 2001 From: "Katherine M. Deck" Date: Wed, 3 Apr 2024 11:33:42 -0700 Subject: [PATCH] force bucket and richards model with ERA5 --- experiments/Artifacts.toml | 15 ++ experiments/Manifest.toml | 6 +- .../Bucket/global_bucket_staticmap.jl | 173 ++++++++++++------ experiments/standalone/Soil/Artifacts.toml | 10 + .../standalone/Soil/richards_runoff.jl | 140 +++++++++++--- 5 files changed, 267 insertions(+), 77 deletions(-) create mode 100644 experiments/Artifacts.toml create mode 100644 experiments/standalone/Soil/Artifacts.toml diff --git a/experiments/Artifacts.toml b/experiments/Artifacts.toml new file mode 100644 index 0000000000..7fb85504f8 --- /dev/null +++ b/experiments/Artifacts.toml @@ -0,0 +1,15 @@ +[lwp_pallardy_etal2018] +git-tree-sha1 = "29ed37c288601ab353f4f5dc13cf652ec4aec305" + +["processed_topographic_index 2.5 deg"] +git-tree-sha1 = "31d4110d008fcfb1ecfd77eb9675c983b650eedf" + +[era5_land_forcing_data2021] +git-tree-sha1 = "ec424296df6b60cfe273ac8f981701fbbed0bd8a" + +[soil_params_Gupta2020_2022] +git-tree-sha1 = "ffc4c3fa76ea632ba46ae5d58671d4161124af29" + + [[soil_params_Gupta2020_2022.download]] + sha256 = "636db48ecb387a9aa59afbff8af86f0ccc092690c688f0cbb52041c133361e22" + url = "https://caltech.box.com/shared/static/f2y23qx0lggjskftzgh7ht7fsbh36gmm.gz" diff --git a/experiments/Manifest.toml b/experiments/Manifest.toml index 4737455c03..b1b75c571d 100644 --- a/experiments/Manifest.toml +++ b/experiments/Manifest.toml @@ -336,10 +336,10 @@ uuid = "595c0a79-7f3d-439a-bc5a-b232dc3bde79" version = "0.7.19" [[deps.ClimaUtilities]] -deps = ["Artifacts", "CFTime", "Dates"] -git-tree-sha1 = "d1b30a07829248b325ffffbb739d0f6fe79dfe7c" +deps = ["Artifacts", "Dates"] +git-tree-sha1 = "617ef474d78e2fd713626e22db8e9fb4881e90ee" uuid = "b3f4f4ca-9299-4f7f-bd9b-81e1242a7513" -version = "0.1.3" +version = "0.1.5" weakdeps = ["Adapt", "ClimaComms", "ClimaCore", "ClimaCoreTempestRemap", "Interpolations", "NCDatasets"] [deps.ClimaUtilities.extensions] diff --git a/experiments/standalone/Bucket/global_bucket_staticmap.jl b/experiments/standalone/Bucket/global_bucket_staticmap.jl index db563dab87..d564c76c59 100644 --- a/experiments/standalone/Bucket/global_bucket_staticmap.jl +++ b/experiments/standalone/Bucket/global_bucket_staticmap.jl @@ -20,9 +20,10 @@ using CairoMakie using Dates using DelimitedFiles using Statistics - +using ClimaUtilities.ClimaArtifacts +import Interpolations import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput - +import ClimaUtilities.ClimaArtifacts: @clima_artifact import ClimaTimeSteppers as CTS import NCDatasets using ClimaCore @@ -55,7 +56,14 @@ function compute_extrema(v) return (minimum(mins), maximum(maxes)) end -anim_plots = false +function compute_clims(v) + means = [mean(u) for u in v] + sigmas = [std(u) for u in v] + return (minimum(means) - maximum(sigmas), maximum(means) + maximum(sigmas)) +end + +anim_plots = true +regridder_type = :InterpolationsRegridder FT = Float64; context = ClimaComms.context() earth_param_set = LP.LandParameters(FT); @@ -70,11 +78,11 @@ soil_depth = FT(3.5); bucket_domain = ClimaLand.Domains.SphericalShell(; radius = FT(6.3781e6), depth = soil_depth, - nelements = (50, 10), + nelements = (30, 10), npolynomial = 1, dz_tuple = FT.((1.0, 0.05)), ); -ref_time = DateTime(2005); +ref_time = DateTime(2021); # Set up parameters σS_c = FT(0.2); @@ -85,43 +93,92 @@ z_0b = FT(1e-3); ρc_soil = FT(2e6); τc = FT(3600); t0 = 0.0; -tf = 7 * 86400; -Δt = 3600.0; +tf = 14 * 86400; +Δt = 3600.0 / 3; # Construct albedo parameter object using static map # Use separate regridding directory for CPU and GPU runs to avoid race condition device_suffix = typeof(ClimaComms.context().device) <: ClimaComms.CPUSingleThreaded ? "cpu" : "gpu" -regrid_dir = joinpath( - pkgdir(ClimaLand), - "experiments/standalone/Bucket/$device_suffix/regrid-static/", -) -!ispath(regrid_dir) && mkpath(regrid_dir) +t_start = t0 surface_space = bucket_domain.space.surface α_snow = FT(0.8) albedo = PrescribedBaregroundAlbedo{FT}(α_snow, surface_space); bucket_parameters = BucketModelParameters(FT; albedo, z_0m, z_0b, τc); +# Forcing data +# TODO: Change with era5_artifact_path = @clima_artifact("era5_land_forcing_data2021", context) +era5_artifact_path = "/groups/esm/ClimaArtifacts/artifacts/era5_land_forcing_data2021" + +# Below, the preprocess_func argument is used to +# 1. Convert precipitation to be negative (as it is downwards) +# 2. Convert accumulations over an hour to a rate per second # Precipitation: -precip = (t) -> 0; -snow_precip = (t) -> -5e-7 * (t < 1 * 86400); -# Diurnal temperature variations: -T_atmos = (t) -> 275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2); -# Constant otherwise: -u_atmos = (t) -> 3.0; -q_atmos = (t) -> 0.001; -h_atmos = FT(2); -P_atmos = (t) -> 101325; +precip = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "rf", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,), +) + +snow_precip = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "sf", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,), +) + +u_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "ws", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, +) +q_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25_clima.nc"), + "q", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, +) +P_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "sp", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, +) + +T_atmos = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "t2m", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, +) +h_atmos = FT(10); + bucket_atmos = PrescribedAtmosphere( - TimeVaryingInput(precip), - TimeVaryingInput(snow_precip), - TimeVaryingInput(T_atmos), - TimeVaryingInput(u_atmos), - TimeVaryingInput(q_atmos), - TimeVaryingInput(P_atmos), + precip, + snow_precip, + T_atmos, + u_atmos, + q_atmos, + P_atmos, ref_time, h_atmos, ); @@ -130,14 +187,26 @@ bucket_atmos = PrescribedAtmosphere( # peak at local noon, and a prescribed downwelling LW radiative # flux, assuming the air temperature is on average 275 degrees # K with a diurnal amplitude of 5 degrees K: -SW_d = (t) -> max(1361 * sin(2π * t / 86400 - π / 2), 0.0); -LW_d = (t) -> 5.67e-8 * (275.0 + 5.0 * sin(2.0 * π * t / 86400 - π / 2))^4; -bucket_rad = PrescribedRadiativeFluxes( - FT, - TimeVaryingInput(SW_d), - TimeVaryingInput(LW_d), - ref_time, -); +SW_d = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "ssrd", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,), +) +LW_d = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "strd", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> data / 3600,), +) + +bucket_rad = PrescribedRadiativeFluxes(FT, SW_d, LW_d, ref_time); model = BucketModel( parameters = bucket_parameters, @@ -214,18 +283,22 @@ F_sfc = [ ) for k in 1:length(sol.t) ]; +sw_forcing = [ + Remapping.interpolate(remapper, saved_values.saveval[k].drivers.SW_d) + for k in 1:length(sol.t) +]; + # save prognostic state to CSV - for comparison between GPU and CPU output open(joinpath(outdir, "tf_state_$device_suffix.txt"), "w") do io writedlm(io, hcat(T_sfc[end][:], W[end][:], Ws[end][:], σS[end][:]), ',') end; # animation settings -nframes = length(W) -framerate = 2 +nframes = length(T_sfc) # hourly data fig_ts = Figure(size = (1000, 1000)) for (i, (field_ts, field_name)) in enumerate( zip( - [W, σS, T_sfc, evaporation, F_sfc], - ["W", "σS", "T_sfc", "evaporation", "F_sfc"], + [W, σS, T_sfc, evaporation, F_sfc, sw_forcing], + ["W", "σS", "T_sfc", "evaporation", "F_sfc", "SW forcing"], ), ) if anim_plots @@ -236,13 +309,18 @@ for (i, (field_ts, field_name)) in enumerate( ylabel = "Latitude", title = field_name, ) - clims = compute_extrema(field_ts) + clims = compute_clims(field_ts) CairoMakie.Colorbar(fig[:, end + 1], colorrange = clims) outfile = joinpath( outdir, string("anim_$(device_suffix)_", field_name, ".mp4"), ) - record(fig, outfile, 1:nframes; framerate = framerate) do frame + record( + fig, + outfile, + (nframes - 7 * 24):2:nframes; + framerate = 3, + ) do frame CairoMakie.heatmap!( longpts, latpts, @@ -251,18 +329,7 @@ for (i, (field_ts, field_name)) in enumerate( ) end end - # Plot the timeseries of the mean value as well. - xlabel = i == 5 ? "Time (days)" : "" - ax2 = Axis( - fig_ts[i, 1], - xlabel = xlabel, - ylabel = field_name, - title = "Global bucket with static map albedo", - ) - CairoMakie.lines!(ax2, sol.t ./ 3600 ./ 24, [mean(x) for x in field_ts]) + end outfile = joinpath(outdir, string("ts_$device_suffix.png")) CairoMakie.save(outfile, fig_ts) - -# delete regrid directory -rm(regrid_dir, recursive = true) diff --git a/experiments/standalone/Soil/Artifacts.toml b/experiments/standalone/Soil/Artifacts.toml new file mode 100644 index 0000000000..5512e5fb0b --- /dev/null +++ b/experiments/standalone/Soil/Artifacts.toml @@ -0,0 +1,10 @@ +["processed_topographic_index 2.5 deg"] +git-tree-sha1 = "0c252917b24ce3c201fa62d639342ab4d872b388" + +[soil_params_Gupta2020_2022] +git-tree-sha1 = "ffc4c3fa76ea632ba46ae5d58671d4161124af29" + + [[soil_params_Gupta2020_2022.download]] + sha256 = "636db48ecb387a9aa59afbff8af86f0ccc092690c688f0cbb52041c133361e22" + url = "https://caltech.box.com/shared/static/f2y23qx0lggjskftzgh7ht7fsbh36gmm.gz" + diff --git a/experiments/standalone/Soil/richards_runoff.jl b/experiments/standalone/Soil/richards_runoff.jl index 41e4e4837e..3e06c73450 100644 --- a/experiments/standalone/Soil/richards_runoff.jl +++ b/experiments/standalone/Soil/richards_runoff.jl @@ -1,10 +1,17 @@ using CairoMakie using Statistics using ArtifactWrappers +using Dates import SciMLBase import ClimaTimeSteppers as CTS using ClimaCore +using ClimaUtilities.ClimaArtifacts +import Interpolations + +import ClimaUtilities.TimeVaryingInputs: TimeVaryingInput import ClimaUtilities.SpaceVaryingInputs: SpaceVaryingInput +import ClimaUtilities.Regridders: InterpolationsRegridder +import ClimaUtilities.ClimaArtifacts: @clima_artifact import NCDatasets import ClimaParams as CP using ClimaComms @@ -14,7 +21,7 @@ using ClimaLand.Soil import ClimaLand import ClimaLand.Parameters as LP - +regridder_type = :InterpolationsRegridder context = ClimaComms.context() outdir = joinpath(pkgdir(ClimaLand), "experiments/standalone/Soil/artifacts") !ispath(outdir) && mkpath(outdir) @@ -43,8 +50,13 @@ infile_path = joinpath(get_data_folder(topmodel_dataset), "means_2.5_new.nc") outfile_root = joinpath(pkgdir(ClimaLand), "experiments/standalone/Soil/static_data_cgll") -f_max = SpaceVaryingInput(infile_path, "fmax", surface_space) -mask = SpaceVaryingInput(infile_path, "landsea_mask", surface_space) +f_max = SpaceVaryingInput(infile_path, "fmax", surface_space; regridder_type) +mask = SpaceVaryingInput( + infile_path, + "landsea_mask", + surface_space; + regridder_type, +) oceans_to_zero(field, mask) = mask > 0.5 ? field : eltype(field)(0) f_over = FT(3.28) # 1/m @@ -54,14 +66,81 @@ runoff_model = ClimaLand.Soil.Runoff.TOPMODELRunoff{FT}(; f_max = f_max, R_sb = R_sb, ) -vg_α = ClimaCore.Fields.ones(subsurface_space) .* FT(0.2) -hydrology_cm = map(vg_α) do (α) - FT = typeof(α) - ClimaLand.Soil.vanGenuchten{FT}(; α = α, n = α + FT(2)) +# TODO: Change with @clima_artifact("soil_params_Gupta2020_2022", context) +soil_params_artifact_path = "/groups/esm/ClimaArtifacts/artifacts/soil_params_Gupta2020_2022" + +extrapolation_bc = + (Interpolations.Periodic(), Interpolations.Flat(), Interpolations.Flat()) +# We use this mask to set values of these parameters over the ocean, in order +# to keep them in the physical range +function mask_vg(var, value) + if var < 1e-8 + return value + else + return var + end +end +vg_α = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "vGalpha_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "α", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) +vg_α .= mask_vg.(vg_α, 1e-3) +# We use this mask to set values of this parameter over the ocean, in order +# to keep it in the physical range +function mask_vg_n(var, value) + if var < 1 + return value + else + return var + end end -θ_r = ClimaCore.Fields.zeros(subsurface_space) -ν = ClimaCore.Fields.zeros(subsurface_space) .+ FT(0.5) -K_sat = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-6) +vg_n = SpaceVaryingInput( + joinpath(soil_params_artifact_path, "vGn_map_gupta_etal2020_2.5x2.5x4.nc"), + "n", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) +vg_n .= mask_vg_n.(vg_n, 1.001) +vg_fields_to_hcm_field(α::FT, n::FT) where {FT} = + ClimaLand.Soil.vanGenuchten{FT}(; @NamedTuple{α::FT, n::FT}((α, n))...) +hydrology_cm = vg_fields_to_hcm_field.(vg_α, vg_n) + +θ_r = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "residual_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "θ_r", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) + +ν = SpaceVaryingInput( + joinpath( + soil_params_artifact_path, + "porosity_map_gupta_etal2020_2.5x2.5x4.nc", + ), + "ν", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) +ν .= mask_vg.(ν, 1.0) +K_sat = SpaceVaryingInput( + joinpath(soil_params_artifact_path, "ksat_map_gupta_etal2020_2.5x2.5x4.nc"), + "Ksat", + subsurface_space; + regridder_type, + regridder_kwargs = (; extrapolation_bc,), +) S_s = ClimaCore.Fields.zeros(subsurface_space) .+ FT(1e-3) soil_params = ClimaLand.Soil.RichardsParameters(; @@ -72,10 +151,24 @@ soil_params = ClimaLand.Soil.RichardsParameters(; θ_r = θ_r, ) -function precip_function(t) - return -1e-6 -end -precip = TimeVaryingInput(precip_function) +# TODO: Change with era5_artifact_path = @clima_artifact("era5_land_forcing_data2021", context) +era5_artifact_path = "/groups/esm/ClimaArtifacts/artifacts/era5_land_forcing_data2021" + +# Below, the preprocess_func argument is used to +# 1. Convert precipitation to be negative (as it is downwards) +# 2. Convert accumulations over an hour to a rate per second +ref_time = DateTime(2021); +t_start = 0.0 +# Precipitation: +precip = TimeVaryingInput( + joinpath(era5_artifact_path, "era5_2021_0.9x1.25.nc"), + "tp", + surface_space; + reference_date = ref_time, + t_start, + regridder_type, + file_reader_kwargs = (; preprocess_func = (data) -> -data / 3600,), +) atmos = ClimaLand.PrescribedPrecipitation{FT, typeof(precip)}(precip) bottom_bc = ClimaLand.Soil.WaterFluxBC((p, t) -> 0.0) bc = (; @@ -116,9 +209,9 @@ function hydrostatic_profile( return FT(ϑ_l) end t0 = 0.0 -tf = 3600.0 * 24 +tf = 3600.0 * 24 * 2 dt = 1800.0 -Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_α .+ 2, S_s, f_max) +Y.soil.ϑ_l .= hydrostatic_profile.(lat, z, ν, θ_r, vg_α, vg_n, S_s, f_max) @. Y.soil.ϑ_l = oceans_to_zero(Y.soil.ϑ_l, mask) set_initial_cache! = make_set_initial_cache(model) exp_tendency! = make_exp_tendency(model); @@ -158,7 +251,7 @@ sv = (; saveval = Array{NamedTuple}(undef, length(saveat)), ) saving_cb = ClimaLand.NonInterpSavingCallback(sv, saveat) -updateat = deepcopy(saveat) +updateat = Array(t0:dt:tf) updatefunc = ClimaLand.make_update_drivers(atmos, nothing) driver_cb = ClimaLand.DriverUpdateCallback(updateat, updatefunc) cb = SciMLBase.CallbackSet(driver_cb, saving_cb) @@ -223,14 +316,16 @@ CairoMakie.heatmap!(ax, longpts, latpts, R_ss_end, colorrange = clims) Colorbar(fig[:, end + 1], colorrange = clims) outfile = joinpath(outdir, string("heatmap_R_ss.png")) CairoMakie.save(outfile, fig) +field_to_error(field) = field < 1 & ~isnan(field) ? field : eltype(field)(-0.1) θ_sfc_end = ClimaCore.Remapping.interpolate( remapper, ClimaLand.Soil.get_top_surface_field( - oceans_to_zero.(sol.u[end].soil.ϑ_l, mask), + oceans_to_zero.(field_to_error.(sol.u[end].soil.ϑ_l), mask), surface_space, ), ) + fig = Figure(size = (1000, 400)) ax = Axis(fig[1, 1], xlabel = "Longitude", ylabel = "Latitude", title = "θ_sfc") clims1 = extrema(θ_sfc_end) @@ -240,7 +335,10 @@ Colorbar(fig[1, 2], colorrange = clims1) Δθ_sfc = ClimaCore.Remapping.interpolate( remapper, ClimaLand.Soil.get_top_surface_field( - oceans_to_zero.(sol.u[end].soil.ϑ_l .- sol.u[1].soil.ϑ_l, mask), + oceans_to_zero.( + field_to_error.(sol.u[end].soil.ϑ_l .- sol.u[1].soil.ϑ_l), + mask, + ), surface_space, ), ) @@ -255,7 +353,7 @@ CairoMakie.save(outfile, fig) int_ϑ_end = similar(sv.saveval[1].soil.h∇) ClimaCore.Operators.column_integral_definite!( int_ϑ_end, - oceans_to_zero.(sol.u[end].soil.ϑ_l, mask), + oceans_to_zero.(field_to_error.(sol.u[end].soil.ϑ_l), mask), ) normalized_int_ϑ_end = ClimaCore.Remapping.interpolate(remapper, int_ϑ_end ./ 50.0) @@ -279,7 +377,7 @@ Colorbar(fig[1, 2], colorrange = clims1) int_ϑ_1 = similar(sv.saveval[1].soil.h∇) ClimaCore.Operators.column_integral_definite!( int_ϑ_1, - oceans_to_zero.(sol.u[1].soil.ϑ_l, mask), + oceans_to_zero.(field_to_error.(sol.u[1].soil.ϑ_l), mask), ) Δ_normalized_int_ϑ = ClimaCore.Remapping.interpolate(remapper, (int_ϑ_end .- int_ϑ_1) ./ 50.0)