modified reckenholz study

.gitignore

@@ -146,15 +146,11 @@ examples/*/*/*/RK4/
 examples/*/*/hydrus_benchmark/*.txt
 examples/*/*/hydrus_benchmark/*.py
 input_examples/
-examples/plot_scale/reckenholz/svat_crop_sobol
 examples/plot_scale/reckenholz/svat_bromide
 examples/plot_scale/reckenholz/svat_nitrate
 examples/plot_scale/reckenholz/svat_crop_bromide
-examples/plot_scale/reckenholz/svat_crop_nitrate
 examples/plot_scale/reckenholz/svat_crop_bromide_sobol
-examples/plot_scale/reckenholz/svat_crop_nitrate_sobol
 examples/plot_scale/reckenholz/svat_crop_bromide_monte_carlo
-examples/plot_scale/reckenholz/svat_crop_nitrate_monte_carlo
 examples/plot_scale/reckenholz/svat_nitrate_monte_carlo
 examples/plot_scale/reckenholz/input
 examples/plot_scale/reckenholz/output

examples/plot_scale/boadkh/svat_crop_nitrate/svat_crop_nitrate.py

@@ -51,12 +51,13 @@ class SVATCROPNITRATESetup(RogerSetup):
         """A SVAT-CROP transport model for nitrate."""
 
         _base_path = Path(__file__).parent
-        if tmp_dir:
-            # read fluxes and states from local SSD on cluster node
-            _input_dir = Path(tmp_dir)
-        else:
-            _input_dir = _base_path.parent / "output" / "svat_crop"
+        # if tmp_dir:
+        #     # read fluxes and states from local SSD on cluster node
+        #     _input_dir = Path(tmp_dir)
+        # else:
+        #     _input_dir = _base_path.parent / "output" / "svat_crop"
         # _input_dir = Path("/Volumes/LaCie/roger/examples/plot_scale/boadkh") / "output" / "svat_crop"
+        _input_dir = _base_path.parent / "output" / "svat_crop"
 
         def _read_var_from_nc(self, var, path_dir, file):
             nc_file = path_dir / file
@@ -289,7 +290,7 @@ def set_parameters_setup(self, state):
             vs.sas_params_re_rg = update(vs.sas_params_re_rg, at[2:-2, 2:-2, 0], 6)
             vs.sas_params_re_rg = update(vs.sas_params_re_rg, at[2:-2, 2:-2, 1], 0.5)
             vs.sas_params_re_rl = update(vs.sas_params_re_rl, at[2:-2, 2:-2, 0], 6)
-            vs.sas_params_re_rl = update(vs.sas_params_re_rl, at[2:-2, 2:-2, 1], 10)
+            vs.sas_params_re_rl = update(vs.sas_params_re_rl, at[2:-2, 2:-2, 1], 3)
 
             # denitrification parameters
             vs.km_denit_rz = update(vs.km_denit_rz, at[2:-2, 2:-2], self._read_var_from_nc("km_denit", self._base_path, "parameters.nc"))
@@ -445,8 +446,8 @@ def set_initial_conditions(self, state):
             )
 
             # initial nitrate concentration (in mg/l)
-            vs.C_rz = update(vs.C_rz, at[2:-2, 2:-2, :vs.taup1], 10)
-            vs.C_ss = update(vs.C_ss, at[2:-2, 2:-2, :vs.taup1], 10)
+            vs.C_rz = update(vs.C_rz, at[2:-2, 2:-2, :vs.taup1], 5.)
+            vs.C_ss = update(vs.C_ss, at[2:-2, 2:-2, :vs.taup1], 5.)
             # initial mineral soil nitrogen
             vs.Nmin_rz = update(vs.Nmin_rz, at[2:-2, 2:-2, :vs.taup1, :], (100 / settings.ages) * settings.dx * settings.dy * 100)
             vs.Nmin_ss = update(vs.Nmin_ss, at[2:-2, 2:-2, :vs.taup1, :], (0 / settings.ages) * settings.dx * settings.dy * 100)
@@ -666,7 +667,6 @@ def set_forcing(self, state):
                     )
                     # set fertilization rate
                     row_no = _get_row_no(vs.lut_fert1[:, 0], i)
-                    print(lut_fert[row_no, 1].dtype)
                     vs.doy_fert1 = update(
                         vs.doy_fert1,
                         at[2:-2, 2:-2],

examples/plot_scale/boadkh/svat_crop_nitrate/write_parameters.py

@@ -27,7 +27,7 @@
 
 _VALS = {
     "alpha_transp": 0.55,
-    "alpha_q": 0.7,
+    "alpha_q": 0.55,
     "c2_transp": 0.6,
     "c2_q_rz": 1.5,
     "c2_q_ss": 1.5,

examples/plot_scale/reckenholz/calculate_clay_content.py

@@ -1,7 +1,4 @@
-
-from pathlib import Path
 import numpy as np
-import pandas as pd
 
 
 theta_pwp = 0.189

examples/plot_scale/reckenholz/svat_crop/compare_simulations_and_observations.py

@@ -47,6 +47,7 @@
         path_fig = base_path_figs / file_str
         fig.savefig(path_fig, dpi=300)
         # plot cumulated observed and simulated time series
+        df_eval = df_eval.dropna()
         fig = eval_utils.plot_obs_sim_cum(df_eval, y_lab=labs._Y_LABS_CUM[var_sim], x_lab='Time [year]')
         file_str = '%s_cum_%s.pdf' % (var_sim, lys_experiment)
         path_fig = base_path_figs / file_str
@@ -56,13 +57,13 @@
         path_fig = base_path_figs / file_str
         fig.savefig(path_fig, dpi=300)
 
-    vars_obs = ['WEIGHT']
-    vars_sim = ['S']
+    vars_obs = ['THETA']
+    vars_sim = ['theta']
     for var_obs, var_sim in zip(vars_obs, vars_sim):
         obs_vals = ds_obs[var_obs].isel(x=0, y=0).values
         df_obs = pd.DataFrame(index=date_obs, columns=['obs'])
         df_obs.loc[:, 'obs'] = obs_vals
-        sim_vals = ds_sim[var_sim].isel(x=0, y=0).values
+        sim_vals = onp.nanmean(ds_sim[var_sim].isel(x=0, y=0).values, axis=-1)
         # join observations on simulations
         df_eval = eval_utils.join_obs_on_sim(date_sim, sim_vals, df_obs)
         # plot observed and simulated time series

examples/plot_scale/reckenholz/svat_crop/parallel_jobs.sh

@@ -0,0 +1,7 @@
+#!/bin/bash
+
+python svat_crop.py -b jax -lys lys2 &
+python svat_crop.py -b jax -lys lys3 &
+python svat_crop.py -b jax -lys lys8 &
+# python svat_crop.py -b jax -lys lys2_bromide &
+# python svat_crop.py -b jax -lys lys8_bromide &

examples/plot_scale/reckenholz/svat_crop/svat_crop.py

@@ -7,7 +7,7 @@
 
 
 @click.option("-lys", "--lys-experiment", type=click.Choice(["lys1", "lys2", "lys3", "lys4", "lys8", "lys9", "lys2_bromide", "lys8_bromide", "lys9_bromide"]), default="lys8")
-@click.option("-td", "--tmp-dir", type=str, default=Path(__file__).parent)
+@click.option("-td", "--tmp-dir", type=str, default=Path(__file__).parent.parent / "output" / "svat_crop")
 @roger_base_cli
 def main(lys_experiment, tmp_dir):
     from roger import RogerSetup, roger_routine, roger_kernel, KernelOutput
@@ -145,6 +145,8 @@ def set_topography(self, state):
                 "z_root_crop",
                 "root_growth_scale",
                 "canopy_growth_scale",
+                "basal_crop_coeff_scale",
+                "pet_weight",
             ],
         )
         def set_parameters_setup(self, state):
@@ -162,15 +164,17 @@ def set_parameters_setup(self, state):
 
             vs.lu_id = update(vs.lu_id, at[2:-2, 2:-2], vs.crop_type[2:-2, 2:-2, 0])
             vs.z_soil = update(vs.z_soil, at[2:-2, 2:-2], 1350)
-            vs.dmpv = update(vs.dmpv, at[2:-2, 2:-2], 25)
+            vs.dmpv = update(vs.dmpv, at[2:-2, 2:-2], 50)
             vs.lmpv = update(vs.lmpv, at[2:-2, 2:-2], 300)
-            vs.theta_ac = update(vs.theta_ac, at[2:-2, 2:-2], 0.1062)
-            vs.theta_ufc = update(vs.theta_ufc, at[2:-2, 2:-2], 0.1247)
+            vs.theta_ac = update(vs.theta_ac, at[2:-2, 2:-2], 0.0962)
+            vs.theta_ufc = update(vs.theta_ufc, at[2:-2, 2:-2], 0.1047)
             vs.theta_pwp = update(vs.theta_pwp, at[2:-2, 2:-2], 0.1869)
-            vs.ks = update(vs.ks, at[2:-2, 2:-2], 1.02)
+            vs.ks = update(vs.ks, at[2:-2, 2:-2], 5.02)
             vs.kf = update(vs.kf, at[2:-2, 2:-2], 2500)
             vs.root_growth_scale = update(vs.root_growth_scale, at[2:-2, 2:-2], 1)
             vs.canopy_growth_scale = update(vs.canopy_growth_scale, at[2:-2, 2:-2], 1)
+            vs.basal_crop_coeff_scale = update(vs.basal_crop_coeff_scale, at[2:-2, 2:-2], 1.5)
+            vs.pet_weight = update(vs.pet_weight, at[2:-2, 2:-2], 1.4)
 
         @roger_routine
         def set_parameters(self, state):
@@ -233,7 +237,7 @@ def set_forcing(self, state):
                 vs.ta_day = update(vs.ta_day, at[:, :, :], vs.TA[npx.newaxis, npx.newaxis, vs.itt_forc:vs.itt_forc+6*24])
                 vs.ta_min = update(vs.ta_min, at[:, :], npx.min(vs.TA_MIN[npx.newaxis, npx.newaxis, vs.itt_forc:vs.itt_forc+6*24], axis=-1))
                 vs.ta_max = update(vs.ta_max, at[:, :], npx.max(vs.TA_MAX[npx.newaxis, npx.newaxis, vs.itt_forc:vs.itt_forc+6*24], axis=-1))
-                vs.pet_day = update(vs.pet_day, at[:, :, :], vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc:vs.itt_forc+6*24])
+                vs.pet_day = update(vs.pet_day, at[:, :, :], vs.PET[npx.newaxis, npx.newaxis, vs.itt_forc:vs.itt_forc+6*24] * vs.pet_weight[:, :, npx.newaxis])
                 vs.itt_forc = vs.itt_forc + 6 * 24
 
             if (vs.year[1] != vs.year[0]) & (vs.itt > 1):