|
| 1 | +# # Global run of land model at low resolution |
| 2 | + |
| 3 | +# The code sets up and runs ClimaLand v1, which |
| 4 | +# includes soil, canopy, and snow, on a spherical domain, |
| 5 | +# using low resolution ERA5 data as forcing. |
| 6 | + |
| 7 | +# Simulation Setup |
| 8 | +# Number of spatial elements: 30 in horizontal, 15 in vertical |
| 9 | +# Soil depth: 50 m |
| 10 | +# Simulation duration: 730 d |
| 11 | +# Timestep: 450 s |
| 12 | +# Timestepper: ARS111 |
| 13 | +# Fixed number of iterations: 3 |
| 14 | +# Jacobian update: every new Newton iteration |
| 15 | +# Atmos forcing update: every 3 hours |
| 16 | + |
| 17 | +import ClimaComms |
| 18 | +ClimaComms.@import_required_backends |
| 19 | +using ClimaUtilities.ClimaArtifacts |
| 20 | +import ClimaUtilities.TimeManager: ITime, date |
| 21 | + |
| 22 | +import ClimaDiagnostics |
| 23 | +import ClimaUtilities |
| 24 | + |
| 25 | +import ClimaUtilities.TimeVaryingInputs: |
| 26 | + TimeVaryingInput, LinearInterpolation, PeriodicCalendar |
| 27 | +import ClimaUtilities.ClimaArtifacts: @clima_artifact |
| 28 | +import ClimaParams as CP |
| 29 | +using ClimaCore |
| 30 | +using ClimaLand |
| 31 | +using ClimaLand.Snow |
| 32 | +using ClimaLand.Soil |
| 33 | +using ClimaLand.Canopy |
| 34 | +import ClimaLand |
| 35 | +import ClimaLand.Parameters as LP |
| 36 | +import ClimaLand.Simulations: LandSimulation, solve! |
| 37 | + |
| 38 | +using Dates |
| 39 | + |
| 40 | +using CairoMakie, GeoMakie, Poppler_jll, ClimaAnalysis |
| 41 | +import ClimaLand.LandSimVis as LandSimVis |
| 42 | + |
| 43 | +const FT = Float64; |
| 44 | +context = ClimaComms.context() |
| 45 | +ClimaComms.init(context) |
| 46 | +device = ClimaComms.device() |
| 47 | +device_suffix = device isa ClimaComms.CPUSingleThreaded ? "cpu" : "gpu" |
| 48 | +root_path = "lowres_snowy_land_longrun_temp_$(device_suffix)" |
| 49 | +diagnostics_outdir = joinpath(root_path, "global_diagnostics") |
| 50 | +outdir = |
| 51 | + ClimaUtilities.OutputPathGenerator.generate_output_path(diagnostics_outdir) |
| 52 | + |
| 53 | +function setup_model(FT, start_date, stop_date, Δt, domain, earth_param_set) |
| 54 | + surface_domain = ClimaLand.Domains.obtain_surface_domain(domain) |
| 55 | + surface_space = domain.space.surface |
| 56 | + # Forcing data |
| 57 | + era5_ncdata_path = ClimaLand.Artifacts.era5_land_forcing_data2008_path(; |
| 58 | + context, |
| 59 | + lowres = true, |
| 60 | + ) |
| 61 | + atmos, radiation = ClimaLand.prescribed_perturbed_temperature_era5( |
| 62 | + era5_ncdata_path, |
| 63 | + surface_space, |
| 64 | + start_date, |
| 65 | + earth_param_set, |
| 66 | + 10.0, |
| 67 | + FT; |
| 68 | + max_wind_speed = 25.0, |
| 69 | + time_interpolation_method = LinearInterpolation(PeriodicCalendar()), |
| 70 | + ) |
| 71 | + forcing = (; atmos, radiation) |
| 72 | + |
| 73 | + # Read in LAI from MODIS data |
| 74 | + modis_lai_ncdata_path = ClimaLand.Artifacts.modis_lai_multiyear_paths(; |
| 75 | + context = nothing, |
| 76 | + start_date, |
| 77 | + end_date = stop_date, |
| 78 | + ) |
| 79 | + LAI = ClimaLand.prescribed_lai_modis( |
| 80 | + modis_lai_ncdata_path, |
| 81 | + surface_space, |
| 82 | + start_date; |
| 83 | + time_interpolation_method = LinearInterpolation(), |
| 84 | + ) |
| 85 | + |
| 86 | + # Overwrite some defaults for the canopy model |
| 87 | + # Energy model |
| 88 | + ac_canopy = FT(2.5e3) |
| 89 | + energy = Canopy.BigLeafEnergyModel{FT}(; ac_canopy) |
| 90 | + |
| 91 | + # Plant hydraulics |
| 92 | + a = FT(0.2 * 0.0098) # 1/m |
| 93 | + retention_model = Canopy.PlantHydraulics.LinearRetentionCurve{FT}(a) |
| 94 | + hydraulics = |
| 95 | + Canopy.PlantHydraulicsModel{FT}(surface_domain, LAI; retention_model) |
| 96 | + |
| 97 | + # Roughness lengths |
| 98 | + h_canopy = hydraulics.compartment_surfaces[end] |
| 99 | + z0_m = FT(0.13) * h_canopy |
| 100 | + z0_b = FT(0.1) * z0_m |
| 101 | + |
| 102 | + ground = ClimaLand.PrognosticGroundConditions{FT}() |
| 103 | + canopy_forcing = (; atmos, radiation, ground) |
| 104 | + |
| 105 | + canopy = ClimaLand.Canopy.CanopyModel{FT}( |
| 106 | + surface_domain, |
| 107 | + canopy_forcing, |
| 108 | + LAI, |
| 109 | + earth_param_set; |
| 110 | + prognostic_land_components = (:canopy, :snow, :soil, :soilco2), |
| 111 | + energy, |
| 112 | + hydraulics, |
| 113 | + z_0m = z0_m, |
| 114 | + z_0b = z0_b, |
| 115 | + ) |
| 116 | + |
| 117 | + # Snow model setup |
| 118 | + # Set β = 0 in order to regain model without density dependence |
| 119 | + α_snow = Snow.ZenithAngleAlbedoModel( |
| 120 | + FT(0.64), |
| 121 | + FT(0.06), |
| 122 | + FT(2); |
| 123 | + β = FT(0.4), |
| 124 | + x0 = FT(0.2), |
| 125 | + ) |
| 126 | + horz_degree_res = |
| 127 | + sum(ClimaLand.Domains.average_horizontal_resolution_degrees(domain)) / 2 # mean of resolution in latitude and longitude, in degrees |
| 128 | + scf = Snow.WuWuSnowCoverFractionModel( |
| 129 | + FT(0.08), |
| 130 | + FT(1.77), |
| 131 | + FT(1.0), |
| 132 | + horz_degree_res, |
| 133 | + ) |
| 134 | + snow = Snow.SnowModel( |
| 135 | + FT, |
| 136 | + surface_domain, |
| 137 | + forcing, |
| 138 | + earth_param_set, |
| 139 | + Δt; |
| 140 | + prognostic_land_components = (:canopy, :snow, :soil, :soilco2), |
| 141 | + α_snow, |
| 142 | + scf, |
| 143 | + ) |
| 144 | + |
| 145 | + # Construct the land model with all default components except for snow |
| 146 | + land = |
| 147 | + LandModel{FT}(forcing, LAI, earth_param_set, domain, Δt; snow, canopy) |
| 148 | + return land |
| 149 | +end |
| 150 | +# Note that since the Northern hemisphere's winter season is defined as DJF, |
| 151 | +# we simulate from and until the beginning of |
| 152 | +# March so that a full season is included in seasonal metrics. |
| 153 | +start_date = DateTime("2008-03-01") |
| 154 | +stop_date = DateTime("2010-03-02") |
| 155 | +Δt = 450.0 |
| 156 | +nelements = (30, 15) |
| 157 | +domain = ClimaLand.Domains.global_domain( |
| 158 | + FT; |
| 159 | + context, |
| 160 | + nelements, |
| 161 | + mask_threshold = FT(0.99), |
| 162 | +) |
| 163 | +params = LP.LandParameters(FT) |
| 164 | +model = setup_model(FT, start_date, stop_date, Δt, domain, params) |
| 165 | +user_callbacks = ( |
| 166 | + ClimaLand.NaNCheckCallback( |
| 167 | + Dates.Month(6), |
| 168 | + start_date, |
| 169 | + ITime(Δt, epoch = start_date), |
| 170 | + mask = ClimaLand.Domains.landsea_mask(ClimaLand.get_domain(model)), |
| 171 | + ), |
| 172 | + ClimaLand.ReportCallback(10000), |
| 173 | +) |
| 174 | +simulation = |
| 175 | + LandSimulation(FT, start_date, stop_date, Δt, model; user_callbacks, outdir) |
| 176 | +@info "Run: Global Soil-Canopy-Snow Model" |
| 177 | +@info "Resolution: $nelements" |
| 178 | +@info "Timestep: $Δt s" |
| 179 | +@info "Start Date: $start_date" |
| 180 | +@info "Stop Date: $stop_date" |
| 181 | +ClimaLand.Simulations.solve!(simulation) |
| 182 | + |
| 183 | +LandSimVis.make_annual_timeseries(simulation; savedir = root_path) |
| 184 | +LandSimVis.make_heatmaps(simulation; savedir = root_path, date = stop_date) |
| 185 | +LandSimVis.make_leaderboard_plots(simulation; savedir = root_path) |
0 commit comments