@@ -50,6 +50,7 @@ using ClimaLand.Soil
50
50
using ClimaLand. Soil. Biogeochemistry
51
51
using ClimaLand. Canopy
52
52
using ClimaLand. Canopy. PlantHydraulics
53
+ import ClimaLand. Simulations: LandSimulation, solve!
53
54
import ClimaLand
54
55
import ClimaLand. Parameters as LP
55
56
using DelimitedFiles
@@ -65,7 +66,11 @@ earth_param_set = LP.LandParameters(FT);
65
66
site_ID = " US-MOz"
66
67
# Timezone (offset from UTC in hrs)
67
68
time_offset = 7
68
- start_date = DateTime (2010 ) + Hour (time_offset)
69
+ start_date = DateTime (2010 ) + Hour (time_offset) + Day (150 ) # start mid year
70
+ # Specify the time range and dt value over which to perform the simulation.
71
+ N_days = 100
72
+ end_date = start_date + Day (N_days)
73
+ dt = Float64 (30 )
69
74
70
75
# Site latitude and longitude
71
76
lat = FT (38.7441 ) # degree
@@ -74,11 +79,6 @@ long = FT(-92.2000) # degree
74
79
# Height of the sensor at the site
75
80
atmos_h = FT (32 ) # m
76
81
77
- # Specify the time range and dt value over which to perform the simulation.
78
- t0 = Float64 (150 * 3600 * 24 )# start mid year
79
- N_days = 100
80
- tf = t0 + Float64 (3600 * 24 * N_days)
81
- dt = Float64 (30 )
82
82
83
83
# - We will be using prescribed atmospheric and radiative drivers from the
84
84
# US-MOz tower, which we read in here. We are using prescribed
@@ -186,8 +186,8 @@ soilco2 = Soil.Biogeochemistry.SoilCO2Model{FT}(domain, drivers);
186
186
# Read in prescribed LAI at the site from global MODIS data
187
187
surface_space = domain. space. surface;
188
188
modis_lai_ncdata_path = ClimaLand. Artifacts. modis_lai_multiyear_paths (;
189
- start_date = start_date + Second (t0) ,
190
- end_date = start_date + Second (t0) + Second (tf) ,
189
+ start_date = start_date,
190
+ end_date = end_date ,
191
191
);
192
192
LAI = ClimaLand. prescribed_lai_modis (
193
193
modis_lai_ncdata_path,
@@ -273,50 +273,40 @@ canopy = ClimaLand.Canopy.CanopyModel{FT}(
273
273
# Now we can combine the soil and canopy models into a single combined model.
274
274
land = SoilCanopyModel {FT} (soilco2, soil, canopy);
275
275
276
- # Now we can initialize the state vectors and model coordinates, and initialize
277
- # the explicit/implicit tendencies as usual. The Richard's equation time
278
- # stepping is done implicitly, while the canopy model may be explicitly stepped,
279
- # so we use an IMEX (implicit-explicit) scheme for the combined model.
280
-
281
- Y, p, coords = initialize (land);
282
- exp_tendency! = make_exp_tendency (land);
283
- imp_tendency! = make_imp_tendency (land);
284
- jacobian! = make_jacobian (land);
285
- jac_kwargs =
286
- (; jac_prototype = ClimaLand. FieldMatrixWithSolver (Y), Wfact = jacobian!);
287
-
288
276
# We need to provide initial conditions for the model:
289
- Y. soil. ϑ_l = FT (0.4 )
290
- Y. soil. θ_i = FT (0.0 )
291
- T_0 = FT (288.7 )
292
- ρc_s =
293
- volumetric_heat_capacity .(
294
- Y. soil. ϑ_l,
295
- Y. soil. θ_i,
296
- land. soil. parameters. ρc_ds,
297
- earth_param_set,
298
- )
299
- Y. soil. ρe_int =
300
- volumetric_internal_energy .(Y. soil. θ_i, ρc_s, T_0, earth_param_set)
301
-
302
- Y. soilco2. C .= FT (0.000412 ) # set to atmospheric co2, mol co2 per mol air
303
-
304
- canopy_retention_model = canopy. hydraulics. parameters. retention_model
305
- ψ_stem_0 = FT (- 1e5 / 9800 )
306
- ψ_leaf_0 = FT (- 2e5 / 9800 )
307
-
308
- S_l_ini =
309
- inverse_water_retention_curve .(
310
- canopy_retention_model,
311
- [ψ_stem_0, ψ_leaf_0],
312
- plant_ν,
313
- plant_S_s,
314
- )
315
- for i in 1 : 2
316
- Y. canopy. hydraulics. ϑ_l.:($ i) .=
317
- augmented_liquid_fraction .(plant_ν, S_l_ini[i])
277
+ function set_ic! (Y, p, t0, model)
278
+ Y. soil. ϑ_l = FT (0.4 )
279
+ Y. soil. θ_i = FT (0.0 )
280
+ T_0 = FT (288.7 )
281
+ ρc_s =
282
+ volumetric_heat_capacity .(
283
+ Y. soil. ϑ_l,
284
+ Y. soil. θ_i,
285
+ land. soil. parameters. ρc_ds,
286
+ earth_param_set,
287
+ )
288
+ Y. soil. ρe_int =
289
+ volumetric_internal_energy .(Y. soil. θ_i, ρc_s, T_0, earth_param_set)
290
+
291
+ Y. soilco2. C .= FT (0.000412 ) # set to atmospheric co2, mol co2 per mol air
292
+
293
+ canopy_retention_model = canopy. hydraulics. parameters. retention_model
294
+ ψ_stem_0 = FT (- 1e5 / 9800 )
295
+ ψ_leaf_0 = FT (- 2e5 / 9800 )
296
+
297
+ S_l_ini =
298
+ inverse_water_retention_curve .(
299
+ canopy_retention_model,
300
+ [ψ_stem_0, ψ_leaf_0],
301
+ plant_ν,
302
+ plant_S_s,
303
+ )
304
+ for i in 1 : 2
305
+ Y. canopy. hydraulics. ϑ_l.:($ i) .=
306
+ augmented_liquid_fraction .(plant_ν, S_l_ini[i])
307
+ end
308
+ evaluate! (Y. canopy. energy. T, atmos. T, t0)
318
309
end
319
- evaluate! (Y. canopy. energy. T, atmos. T, t0)
320
310
321
311
# Choose the timestepper and solver needed for the problem.
322
312
timestepper = CTS. ARS343 ()
@@ -328,55 +318,45 @@ ode_algo = CTS.IMEXAlgorithm(
328
318
),
329
319
);
330
320
331
- # Now set the initial values for the cache variables for the combined soil and plant model.
332
- set_initial_cache! = make_set_initial_cache (land)
333
- set_initial_cache! (p, Y, t0);
334
321
335
- # Set the callbacks, which govern
322
+ # Set the callbacks update times , which govern
336
323
# how often we save output, and how often we update
337
324
# the forcing data ("drivers")
338
325
339
326
n = 120
340
- saveat = Array (t0 : (n * dt): tf)
327
+ saveat = Array (start_date : Second (n * dt): end_date);
341
328
sv = (;
342
- t = Array {Float64 } (undef, length (saveat)),
329
+ t = Array {DateTime } (undef, length (saveat)),
343
330
saveval = Array {NamedTuple} (undef, length (saveat)),
344
331
)
345
332
saving_cb = ClimaLand. NonInterpSavingCallback (sv, saveat)
346
- model_drivers = ClimaLand. get_drivers (land)
347
- updatefunc = ClimaLand. make_update_drivers (model_drivers)
348
- updateat = Array (t0: 1800 : tf)
349
- driver_cb = ClimaLand. DriverUpdateCallback (updateat, updatefunc)
350
- cb = SciMLBase. CallbackSet (driver_cb, saving_cb);
351
-
352
- # Carry out the simulation
353
- prob = SciMLBase. ODEProblem (
354
- CTS. ClimaODEFunction (
355
- T_exp! = exp_tendency!,
356
- T_imp! = SciMLBase. ODEFunction (imp_tendency!; jac_kwargs... ),
357
- dss! = ClimaLand. dss!,
358
- ),
359
- Y,
360
- (t0, tf),
361
- p,
362
- );
363
- sol = SciMLBase. solve (
364
- prob,
365
- ode_algo;
366
- dt = dt,
367
- callback = cb,
368
- adaptive = false ,
369
- saveat = saveat,
333
+ updateat = Array (start_date: Second (1800 ): end_date);
334
+
335
+ # Create the LandSimulation object, which will also create and initialize the state vectors,
336
+ # the cache, the driver callbacks, and set the initial conditions.
337
+ simulation = LandSimulation (
338
+ start_date,
339
+ end_date,
340
+ dt,
341
+ land;
342
+ set_ic! = set_ic!,
343
+ updateat = updateat,
344
+ solver_kwargs = (; saveat = deepcopy (saveat)),
345
+ timestepper = ode_algo,
346
+ user_callbacks = (saving_cb,),
347
+ diagnostics = (),
370
348
);
371
349
350
+ sol = solve! (simulation);
351
+
372
352
# # Plotting
373
353
374
354
# Now that we have both a soil and canopy model incorporated together, we will
375
355
# show how to plot some model data demonstrating the time series produced from
376
356
# each of these models. As before, we may plot the GPP of the system as well
377
357
# as transpiration showing fluxes in the canopy.
378
358
379
- daily = sol. t ./ 3600 ./ 24
359
+ daily = FT .( sol. t) ./ 3600 ./ 24
380
360
model_GPP = [
381
361
parent (sv. saveval[k]. canopy. photosynthesis. GPP)[1 ] for
382
362
k in 1 : length (sv. saveval)
0 commit comments