save additional variable into csv file

EcoExtreML · May 7, 2024 · 1cdca8f · 1cdca8f
1 parent ddcf0b2
commit 1cdca8f
Show file tree

Hide file tree

Showing 5 changed files with 31 additions and 6 deletions.
diff --git a/src/+io/bin_to_csv.m b/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'};

diff --git a/src/+io/create_output_files_binary.m b/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

diff --git a/src/+io/output_data_binary.m b/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);

diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m
@@ -443,7 +443,7 @@
 
                 switch options.calc_ebal
                     case 1
-                        [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, sfactor, psiLeaf] ...
+                        [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,8 +481,6 @@
                             end
                         end
                 end
-                sfactorTot(KT) = sfactor;
-                psiLeafTot(KT) = psiLeaf;
 
                 if options.calc_fluor % total emitted fluorescence irradiance (excluding leaf and canopy re-absorption and scattering)
                     if options.calc_PSI
@@ -717,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

diff --git a/src/ebal.m b/src/ebal.m
@@ -1,4 +1,4 @@
-function [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, sfactor, PSI]             ...
+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)
@@ -568,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