Feature/correct dummy climate

docs/source/data_recipes/CDS_toolbox_template.md

@@ -48,7 +48,7 @@ and optionally:
 We recommend the following data sets to force the virtual rainforest microclimate
 simulations:
 
-* ERA5
+* ERA5 / ERA5-Land
 
   ERA5 is the fifth generation ECMWF reanalysis for the global climate and weather for
   the past 4 to 7 decades. This reanalysis dataset combines model data with
@@ -61,6 +61,14 @@ simulations:
   [here for hourly data](https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels?tab=overview)
   and [here for monthly data](https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single-levels-monthly-means?tab=overview)
 
+  ERA5-Land is a reanalysis dataset providing a consistent view of the evolution of land
+  variables over several decades at an enhanced resolution compared to ERA5 (0.1 x 0.1
+  deg resolution).
+
+  The full documentation and download link can be accessed
+  [here for hourly data](https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land?tab=overview)
+  and [here for monthly data](https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land-monthly-means?tab=overview)
+
 * WFDE5
 
   This global dataset provides bias-corrected reconstruction of near-surface

docs/source/data_recipes/ERA5_preprocessing_example.md

@@ -64,11 +64,11 @@ dataset = xr.open_dataset("./ERA5_land.nc")
 dataset
 ```
 
-### 2. Convert temperatures to Celsius
+### 2. Convert temperatures units
 
 The standard output unit of ERA5-Land temperatures is Kelvin which we need to convert
-to degree Celsius for the Virtual Rainforest. This includes `2m air temperature` and
-`2m dewpoint temperature` which are used to calculate relative humidity in next step.
+to degree Celsius for the Virtual Rainforest. This includes 2m air temperature and
+2m dewpoint temperature which are used to calculate relative humidity in next step.
 
 ```{code-cell} ipython3
 dataset["t2m_C"] = dataset["t2m"]-273.15 # 2m air temperature
@@ -92,24 +92,43 @@ dataset["rh2m"] = (
     )
 ```
 
-### 4. Clean dataset and rename variables
+### 4. Convert precipitation units
 
-In this step, we delete the initial temperature variables (K) and rename the remaining
-variables according to the Virtual Rainforest naming convention (see
-[here](../../../virtual_rainforest/data_variables.toml) ).
+The standard output unit for total precipitation in ERA5-Land is meters which we need to
+convert to millimeters.
 
 ```{code-cell} ipython3
-dataset_cleaned = dataset.drop_vars(["d2m",'d2m_C',"t2m"])
+dataset["tp_mm"] = dataset["tp"]*1000
+```
+
+### 5. Convert surface pressure units
+
+The standard output unit for surface pressure in ERA5-Land is Pascal (Pa) which we need
+to convert to Kilopascal (kPa).
+
+```{code-cell} ipython3
+dataset["sp_kPa"] = dataset["sp"]/1000
+```
+
+### 6. Clean dataset and rename variables
+
+In this step, we delete the initial temperature variables (K), precipitation (m), and
+surface pressure(Pa) and rename the remaining variables according to the Virtual
+Rainforest naming convention
+(see [here](../../../virtual_rainforest/data_variables.toml) ).
+
+```{code-cell} ipython3
+dataset_cleaned = dataset.drop_vars(["d2m",'d2m_C',"t2m","tp","sp"])
 dataset_renamed = dataset_cleaned.rename({
-    'sp':'atmospheric_pressure_ref',
-    'tp':'precipitation',
+    'sp_kPa':'atmospheric_pressure_ref',
+    'tp_mm':'precipitation',
     't2m_C':'air_temperature_ref',
     'rh2m': 'relative_humidity_ref',
     })
 dataset_renamed.data_vars
 ```
 
-### 5. Add further required variables
+### 7. Add further required variables
 
 In addition to the variables from the ERA5-Land datasset, a time series of atmospheric
 $\ce{CO_{2}}$ is needed. We add this here as a constant field. Mean annual temperature
@@ -127,7 +146,7 @@ dataset_renamed['mean_annual_temperature'] = (
 dataset_renamed.data_vars
 ```
 
-### 7. Change coordinates to x-y in meters
+### 8. Change coordinates to x-y in meters
 
 The following code segment changes the coordinate names from `longitude/latitude` to
 `x/y` and the units from `minutes` to `meters`. The ERA5-Land coordinates are treated as
@@ -139,7 +158,7 @@ dataset_xy = dataset_renamed.rename_dims({'longitude':'x','latitude':'y'}).assig
 dataset_xy.coords
 ```
 
-### 8. Scale to 90 m resolution
+### 9. Scale to 90 m resolution
 
 The Virtual Rainforest is run on a 90 x 90 m grid. This means that some form of spatial
 downscaling has to be applied to the dataset, for example by spatially interpolating
@@ -157,7 +176,7 @@ dataset_xy_dummy = dataset_xy_100.isel(x=np.arange(9),y=np.arange(9))
 dataset_xy_dummy.coords
 ```
 
-### 9. Add time_index
+### 10. Add time_index
 
 The dummy model iterates over time indices rather than real datetime. Therefore, we add
 a `time_index` coordinate to the dataset:
@@ -167,7 +186,7 @@ dataset_xy_timeindex = dataset_xy_dummy.assign_coords({'time_index':np.arange(0,
 dataset_xy_timeindex.coords
 ```
 
-### 10. Save netcdf
+### 11. Save netcdf
 
 Once we confirmed that our dataset is complete and our calculations are correct, we save
 it as a new netcdf file. This can then be fed into the code data loading system here

tests/conftest.py

@@ -254,7 +254,7 @@ def dummy_climate_data(layer_roles_fixture):
         dims=["cell_id", "time_index"],
     )
     data["atmospheric_pressure_ref"] = DataArray(
-        np.full((3, 3), 96000),
+        np.full((3, 3), 96),
         dims=["cell_id", "time_index"],
     )
     data["atmospheric_co2_ref"] = DataArray(

virtual_rainforest/models/abiotic_simple/microclimate.py

@@ -1,9 +1,9 @@
 r"""The ``models.abiotic_simple.microclimate`` module uses linear regressions from
 :cite:t:`hardwick_relationship_2015` and :cite:t:`jucker_canopy_2018` to predict
-atmospheric temperature and relative humidity at ground level (2m) given the above
-canopy conditions and leaf area index of intervening canopy. A within canopy profile is
-then interpolated using a logarithmic curve between the above canopy observation and
-ground level prediction.
+atmospheric temperature, relative humidity, and vapour pressure deficit at ground level
+(1.5 m) given the above canopy conditions and leaf area index of intervening canopy. A
+within canopy profile is then interpolated using a logarithmic curve between the above
+canopy observation and ground level prediction.
 Soil temperature is interpolated between the surface layer and the soil temperature at
 1 m depth which equals the mean annual temperature.
 The module also provides a constant vertical profile of atmospheric pressure and
@@ -63,17 +63,19 @@ def run_microclimate(
     r"""Calculate simple microclimate.
 
     This function uses empirical relationships between leaf area index (LAI) and
-    atmospheric temperature and relative humidity to derive logarithmic
-    profiles of atmospheric temperature and humidity from external climate data such as
-    regional climate models or satellite observations. For below canopy values (1.5 m),
+    atmospheric temperature, relative humidity, and vapour pressure deficit to derive
+    logarithmic profiles of these variables from external climate data such as
+    regional climate models or satellite observations. Note that these sources provide
+    data at different heights and with different underlying assumptions which lead to
+    different biases in the model output. For below canopy values (1.5 m),
     the implementation is based on :cite:t:`hardwick_relationship_2015` as
 
     :math:`y = m * LAI + c`
 
     where :math:`y` is the variable of interest, math:`m` is the gradient
     (:data:`~virtual_rainforest.models.abiotic_simple.microclimate.MicroclimateGradients`)
-    and :math:`c` is the intersect which we set to the
-    external data values. We assume that the gradient remains constant.
+    and :math:`c` is the intersect which we set to the external data values. We assume
+    that the gradient remains constant.
 
     The other atmospheric layers are calculated by logaritmic regression and
     interpolation between the input at the top of the canopy and the 1.5 m values.
@@ -84,7 +86,7 @@ def run_microclimate(
 
     The `layer_roles` list is composed of the following layers (index 0 above canopy):
 
-    * above canopy (canopy height + reference measurement height, typically 2m)
+    * above canopy (canopy height)
     * canopy layers (maximum of ten layers, minimum one layers)
     * subcanopy (1.5 m)
     * surface layer
@@ -95,7 +97,7 @@ def run_microclimate(
     * air_temperature_ref [C]
     * relative_humidity_ref []
     * vapour_pressure_deficit_ref [kPa]
-    * atmospheric_pressure_ref [Pa]
+    * atmospheric_pressure_ref [kPa]
     * atmospheric_co2_ref [ppm]
     * leaf_area_index [m m-1]
     * layer_heights [m]
@@ -116,7 +118,6 @@ def run_microclimate(
         atmospheric :math:`\ce{CO2}` [ppm]
     """
 
-    # TODO correct gap between 1.5 m and 2m reference height for LAI = 0
     # TODO make sure variables are representing correct time interval, e.g. mm per day
     output = {}
 
@@ -138,7 +139,7 @@ def run_microclimate(
 
     # Mean atmospheric pressure profile, [kPa]
     output["atmospheric_pressure"] = (
-        (data["atmospheric_pressure_ref"] / 1000)
+        (data["atmospheric_pressure_ref"])
         .isel(time_index=time_index)
         .where(output["air_temperature"].coords["layer_roles"] != "soil")
         .rename("atmospheric_pressure")
@@ -267,7 +268,7 @@ def calculate_saturation_vapour_pressure(
         factor3: factor 3 in saturation vapour pressure calculation
 
     Returns:
-        saturation vapour pressure, kPa
+        saturation vapour pressure, [kPa]
     """
 
     return DataArray(