Merge pull request #229 from EcoExtreML/save_wsf_psiLeaf

save soil water stress factor and leaf water potential in outputs

CHANGELOG.md

@@ -4,6 +4,8 @@
 
 - Added changes to support groundwater coupling via BMI in
   [#221](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/221)
+- Save water stress factor and water potential into csv files.
+  [#229](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/229)
 
 
 <a name="1.5.0"></a>

src/+io/bin_to_csv.m

@@ -49,6 +49,16 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer)
     write_output(Sim_Temp_names, Sim_Temp_units, ...
                  fnames.Sim_Temp_file, n_col.Sim_Temp, ns, true);
 
+    %% water stress factor
+    waterStressFactor_names = {'simulation_number', 'year', 'DoY', 'soilWaterStressFactor'};
+    waterStressFactor_units = {'', '', '', '-'};
+    write_output(waterStressFactor_names, waterStressFactor_units, fnames.waterStressFactor_file, n_col.waterStressFactor, ns);
+
+    %% water potential
+    waterPotential_names = {'simulation_number', 'year', 'DoY', 'leafWaterPotential'};
+    waterPotential_units = {'', '', '', 'm'};
+    write_output(waterPotential_names, waterPotential_units, fnames.waterPotential_file, n_col.waterPotential, ns);
+
     %% spectrum (added on 19 September 2008)
 
     spectrum_hemis_optical_names = {'hemispherically integrated radiation spectrum'};

src/+io/create_output_files_binary.m

@@ -25,6 +25,9 @@
     fnames.BOC_irradiance_file = fullfile(Output_dir, 'BOC_irradiance.bin');  % reflectance spectrum
     fnames.Sim_Theta_file = fullfile(Output_dir, 'Sim_Theta.bin');  % soil moisture
     fnames.Sim_Temp_file = fullfile(Output_dir, 'Sim_Temp.bin');  % soil temperature
+    fnames.waterStressFactor_file = fullfile(Output_dir, 'waterStressFactor.bin');
+    fnames.waterPotential_file = fullfile(Output_dir, 'waterPotential.bin');
+
     if options.calc_ebal
         fnames.spectrum_obsdir_BlackBody_file = fullfile(Output_dir, 'spectrum_obsdir_BlackBody.bin');  % spectrum observation direction
     end

src/+io/output_data_binary.m

@@ -1,4 +1,4 @@
-function n_col = output_data_binary(f, k, xyt, rad,  canopy, ScopeParameters, vi, vmax, options, fluxes, meteo, iter, thermal, spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap)
+function n_col = output_data_binary(f, k, xyt, rad,  canopy, ScopeParameters, vi, vmax, options, fluxes, meteo, iter, thermal, spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap, WaterStress, WaterPotential)
     %% OUTPUT DATA
     % author C. Van der Tol
     % date:      30 Nov 2019
@@ -32,6 +32,17 @@
     Sim_Temp_out =  [Sim_Temp(k, :)];
     n_col.Sim_Temp = length(Sim_Temp_out);
     fwrite(f.Sim_Temp_file, Sim_Temp_out, 'double');
+
+    %% water stress factor
+    waterStressFactor_out = [k xyt.year(k) xyt.t(k) WaterStress.soil];
+    n_col.waterStressFactor = length(waterStressFactor_out);
+    fwrite(f.waterStressFactor_file, waterStressFactor_out, 'double');
+
+    %% water potential
+    waterPotential_out = [k xyt.year(k) xyt.t(k) WaterPotential.leaf];
+    n_col.waterPotential = length(waterPotential_out);
+    fwrite(f.waterPotential_file, waterPotential_out, 'double');
+
     %% spectrum (added on 19 September 2008)
     spectrum_hemis_optical_out =  rad.Eout_;
     n_col.spectrum_hemis_optical = length(spectrum_hemis_optical_out);

src/STEMMUS_SCOPE.m

@@ -443,7 +443,7 @@
 
                 switch options.calc_ebal
                     case 1
-                        [iter, fluxes, rad, thermal, profiles, soil, RWU, frac] ...
+                        [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, WaterStressFactor, WaterPotential] ...
                             = ebal(iter, options, spectral, rad, gap,  ...
                                    leafopt, angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t, ...
                                    Rl, SoilVariables, VanGenuchten, InitialValues);
@@ -481,6 +481,7 @@
                             end
                         end
                 end
+
                 if options.calc_fluor % total emitted fluorescence irradiance (excluding leaf and canopy re-absorption and scattering)
                     if options.calc_PSI
                         rad.Femtot = 1E3 * (leafbio.fqe(2) * optipar.phiII(spectral.IwlF) * fluxes.aPAR_Cab_eta + leafbio.fqe(1) * optipar.phiI(spectral.IwlF)  * fluxes.aPAR_Cab);
@@ -714,7 +715,7 @@
 
         % Open files for writing
         file_ids = structfun(@(x) fopen(x, 'a'), fnames, 'UniformOutput', false);
-        n_col = io.output_data_binary(file_ids, k, xyt, rad, canopy, ScopeParameters, vi, vmax, options, fluxes, meteo, iter, thermal, spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap);
+        n_col = io.output_data_binary(file_ids, k, xyt, rad, canopy, ScopeParameters, vi, vmax, options, fluxes, meteo, iter, thermal, spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap, WaterStressFactor, WaterPotential);
         fclose("all");
     end
 end

src/ebal.m

@@ -1,4 +1,4 @@
-function [iter, fluxes, rad, thermal, profiles, soil, RWU, frac]             ...
+function [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, WaterStressFactor, WaterPotential]             ...
          = ebal(iter, options, spectral, rad, gap, leafopt,  ...
                 angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t, ...
                 Rl, SoilVariables, VanGenuchten, InitialValues)
@@ -88,6 +88,8 @@ resistance of leaves (or biochemical_MD12: alternative)
        rad         radiation spectra
        profiles    vertical profiles of fluxes
        thermal     temperatures, aerodynamic resistances and friction velocity
+       sfactor     soil water stress factor
+       PSI         leaf water potential
     %}
 
     %% 1. initialisations and other preparations for the iteration loop
@@ -566,4 +568,7 @@ resistance of leaves (or biochemical_MD12: alternative)
     % function Tnew = update(Told, Wc, innovation)
     %     Tnew        = Wc.*innovation + (1-Wc).*Told;
     % return
+
+    WaterStressFactor.soil = sfactor;
+    WaterPotential.leaf = PSI;
 end