Merge pull request #273 from EcoExtreML/add_Trap_to_BMI

MostafaGomaa93 · web-flow · commit 8d3ef7532577 · 2024-12-17T14:39:39.000+01:00
Redefine Evap and Trap arrays
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -15,9 +15,13 @@ The 1.6.1 release comes with minor changes as:
 **Fixed:**
 - The bug in the calculation of the soil type index discussed in
   [#252](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/252)
+- The EVAP is removed discussed in [#274](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/274) and fixed in [#273](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/273)
+- The shape of Evap variable is changed to match BMI requirements discussed in [#274](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/274) and fixed in [#273](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/273)
 
 **Added:**
 - The recharge index is exposed to BMI in [#257](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/257)
+- The Trap variable (after changing its shape) is exposed to BMI discussed in [#271](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/271) and added in [#273](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/273)
+- The FullCSVfiles option is added in [#273](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/273)
 
 
 <a name="1.6.0"></a>
diff --git a/docs/getting_started.md b/docs/getting_started.md
@@ -76,6 +76,13 @@ EndTime=2001-01-02T00:00
 InputPath=/path_to_model_input_folder/
 OutputPath=/path_to_model_output_folder/
 ```
+The configuration file could also contain the optional key `FullCSVfiles`. If
+`FullCSVfiles=0`, the model will **not** store some **large** binary files in
+csv format. The default value is `FullCSVfiles=1`. To know which files are
+stored in csv format, see the function `bin_to_csv()` in `src/+io` folder.
+```text
+FullCSVfiles=0
+```
 
 See example configuration files below:
 
diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE
diff --git a/src/+energy/calculateBoundaryConditions.m b/src/+energy/calculateBoundaryConditions.m
@@ -1,5 +1,5 @@
 function [RHS, EnergyMatrices] = calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, SoilVariables, ...
-                                                             Precip, EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings)
+                                                             Precip, Evap, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings)
     %{
         Determine the boundary condition for solving the energy equation, see
         STEMMUS Technical Notes.
@@ -46,6 +46,6 @@
     else
         L_ts = L(n);
         SH = 0.1200 * Constants.c_a * (SoilVariables.T(n) - ForcingData.Ta_msr(KT)) / r_a_SOIL(KT);
-        RHS(n) = RHS(n) + 100 * Rn_SOIL(KT) / 1800 - Constants.RHOL * L_ts * EVAP(KT) - SH + Constants.RHOL * Constants.c_L * (ForcingData.Ta_msr(KT) * Precip(KT) + BoundaryCondition.DSTOR0 * SoilVariables.T(n) / Delt_t);    % J cm-2 s-1
+        RHS(n) = RHS(n) + 100 * Rn_SOIL(KT) / 1800 - Constants.RHOL * L_ts * Evap - SH + Constants.RHOL * Constants.c_L * (ForcingData.Ta_msr(KT) * Precip(KT) + BoundaryCondition.DSTOR0 * SoilVariables.T(n) / Delt_t);    % J cm-2 s-1
     end
 end
diff --git a/src/+init/defineInitialValues.m b/src/+init/defineInitialValues.m
@@ -178,8 +178,7 @@
 
     %% Structure 5: variables with zeros(Nmsrmn / 10, 1)
     fields = {
-              'Tp_t', 'Evap', 'Tbtm', 'r_a_SOIL', 'Rn_SOIL', 'EVAP', ...
-              'PSItot', 'sfactortot', 'Tsur'
+              'Tbtm', 'r_a_SOIL', 'Rn_SOIL', 'PSItot', 'sfactortot', 'Tsur'
              };
     structures{5} = helpers.createStructure(zeros(Dur_tot, 1), fields);
 
diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m
@@ -1,4 +1,4 @@
-function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings)
+function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings, FullCSVfiles)
     %% flu
     flu_names = {'simulation_number', 'nu_iterations', 'year', 'DoY', ...
                  'Rntot', 'lEtot', 'Htot', ...
@@ -91,60 +91,6 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings)
     Sim_qtot_units = repelem({'cm s-1'}, length(depth));
     write_output(Sim_qtot_names, Sim_qtot_units, fnames.Sim_qtot_file, n_col.Sim_qtot, ns, true);
 
-    %% Spectrum (added on 19 September 2008)
-    spectrum_hemis_optical_names = {'hemispherically integrated radiation spectrum'};
-    spectrum_hemis_optical_units = {'W m-2 um-1'};
-    write_output(spectrum_hemis_optical_names, spectrum_hemis_optical_units, fnames.spectrum_hemis_optical_file, n_col.spectrum_hemis_optical, ns, true);
-
-    spectrum_obsdir_optical_names = {'radiance spectrum in observation direction'};
-    spectrum_obsdir_optical_units = {'W m-2 sr-1 um-1'};
-    write_output(spectrum_obsdir_optical_names, spectrum_obsdir_optical_units, fnames.spectrum_obsdir_optical_file, n_col.spectrum_obsdir_optical, ns, true);
-
-    if options.calc_ebal
-        write_output({'thermal BlackBody emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ...
-                     fnames.spectrum_obsdir_BlackBody_file, n_col.spectrum_obsdir_BlackBody, ns, true);
-        if options.calc_planck
-            write_output({'thermal emission spectrum in hemispherical direction'}, {'W m-2 sr-1 um-1'}, ...
-                         fnames.spectrum_hemis_thermal_file, n_col.spectrum_hemis_thermal, ns, true);
-            write_output({'thermal emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ...
-                         fnames.spectrum_obsdir_thermal_file, n_col.spectrum_obsdir_thermal, ns, true);
-        end
-    end
-    write_output({'irradiance'}, {'W m-2 um-1'}, ...
-                 fnames.irradiance_spectra_file, n_col.irradiance_spectra, ns, true);
-    write_output({'reflectance'}, {'fraction of radiation in observation direction *pi / irradiance'}, ...
-                 fnames.reflectance_file, n_col.reflectance, ns, true);
-    %% input and parameter values (added June 2012)
-    % write_output(fnames.pars_and_input_file, true)
-    % write_output(fnames.pars_and_input_short_file, true)
-    %% Optional Output
-    if options.calc_vert_profiles
-        write_output({'Fraction leaves in the sun, fraction of observed, fraction of observed&visible per layer'}, {'rows: simulations or time steps, columns: layer numbers'}, ...
-                     fnames.gap_file, n_col.gap, ns, true);
-        write_output({'aPAR per leaf layer'}, {'umol m-2 s-1'}, ...
-                     fnames.layer_aPAR_file, n_col.layer_aPAR, ns, true);
-        write_output({'aPAR by Cab per leaf layer'}, {'umol m-2 s-1'}, ...
-                     fnames.layer_aPAR_Cab_file, n_col.layer_aPAR_Cab, ns, true);
-        if options.calc_ebal
-            write_output({'leaf temperature of sunlit leaves, shaded leaves, and weighted average leaf temperature per layer'}, {'^oC ^oC ^oC'}, ...
-                         fnames.leaftemp_file, n_col.leaftemp, ns, true);
-            write_output({'sensible heat flux per layer'}, {'W m-2'}, ...
-                         fnames.layer_H_file, n_col.layer_H, ns, true);
-            write_output({'latent heat flux per layer'}, {'W m-2'}, ...
-                         fnames.layer_LE_file, n_col.layer_LE, ns, true);
-            write_output({'photosynthesis per layer'}, {'umol m-2 s-1'}, ...
-                         fnames.layer_A_file, n_col.layer_A, ns, true);
-            write_output({'average NPQ = 1-(fm-fo)/(fm0-fo0), per layer'}, {''}, ...
-                         fnames.layer_NPQ_file, n_col.layer_NPQ, ns, true);
-            write_output({'net radiation per leaf layer'}, {'W m-2'}, ...
-                         fnames.layer_Rn_file, n_col.layer_Rn, ns, true);
-        end
-        if options.calc_fluor
-            fluorescence_names = {'supward fluorescence per layer'};
-            fluorescence_units = {'W m-2'};
-            write_output(fluorescence_names, fluorescence_units, fnames.layer_fluorescence_file, n_col.layer_fluorescence, ns, true);
-        end
-    end
     if options.calc_fluor
         write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
                      fnames.fluorescence_file, n_col.fluorescence, ns, true);
@@ -154,21 +100,83 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings)
             write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution, for PSII only'}, {'W m-2 um-1 sr-1'}, ...
                          fnames.fluorescencePSII_file, n_col.fluorescencePSII, ns, true);
         end
-        write_output({'hemispherically integrated fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ...
-                     fnames.fluorescence_hemis_file, n_col.fluorescence_hemis, ns, true);
-        write_output({'total emitted fluorescence by all leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ...
-                     fnames.fluorescence_emitted_by_all_leaves_file, n_col.fluorescence_emitted_by_all_leaves, ns, true);
-        write_output({'total emitted fluorescence by all photosystems for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ...
-                     fnames.fluorescence_emitted_by_all_photosystems_file, n_col.fluorescence_emitted_by_all_photosystems, ns, true);
-        write_output({'TOC fluorescence contribution from sunlit leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
-                     fnames.fluorescence_sunlit_file, n_col.fluorescence_sunlit, ns, true);
-        write_output({'TOC fluorescence contribution from shaded leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
-                     fnames.fluorescence_shaded_file, n_col.fluorescence_shaded, ns, true);
-        write_output({'TOC fluorescence contribution from from leaves and soil after scattering for wavelenghts of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
-                     fnames.fluorescence_scattered_file, n_col.fluorescence_scattered, ns, true);
     end
-    write_output({'Bottom of canopy irradiance in the shaded fraction, and average BOC irradiance'}, {'First 2162 columns: shaded fraction. Last 2162 columns: average BOC irradiance. Unit: Wm-2 um-1'}, ...
-                 fnames.BOC_irradiance_file, n_col.BOC_irradiance, ns, true);
+
+    % Optional for large output files
+    if FullCSVfiles
+        %% Spectrum (added on 19 September 2008)
+        spectrum_hemis_optical_names = {'hemispherically integrated radiation spectrum'};
+        spectrum_hemis_optical_units = {'W m-2 um-1'};
+        write_output(spectrum_hemis_optical_names, spectrum_hemis_optical_units, fnames.spectrum_hemis_optical_file, n_col.spectrum_hemis_optical, ns, true);
+
+        spectrum_obsdir_optical_names = {'radiance spectrum in observation direction'};
+        spectrum_obsdir_optical_units = {'W m-2 sr-1 um-1'};
+        write_output(spectrum_obsdir_optical_names, spectrum_obsdir_optical_units, fnames.spectrum_obsdir_optical_file, n_col.spectrum_obsdir_optical, ns, true);
+
+        if options.calc_ebal
+            write_output({'thermal BlackBody emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ...
+                         fnames.spectrum_obsdir_BlackBody_file, n_col.spectrum_obsdir_BlackBody, ns, true);
+            if options.calc_planck
+                write_output({'thermal emission spectrum in hemispherical direction'}, {'W m-2 sr-1 um-1'}, ...
+                             fnames.spectrum_hemis_thermal_file, n_col.spectrum_hemis_thermal, ns, true);
+                write_output({'thermal emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ...
+                             fnames.spectrum_obsdir_thermal_file, n_col.spectrum_obsdir_thermal, ns, true);
+            end
+        end
+        write_output({'irradiance'}, {'W m-2 um-1'}, ...
+                     fnames.irradiance_spectra_file, n_col.irradiance_spectra, ns, true);
+        write_output({'reflectance'}, {'fraction of radiation in observation direction *pi / irradiance'}, ...
+                     fnames.reflectance_file, n_col.reflectance, ns, true);
+
+        %% input and parameter values (added June 2012)
+        % write_output(fnames.pars_and_input_file, true)
+        % write_output(fnames.pars_and_input_short_file, true)
+
+        %% Optional Output
+        if options.calc_vert_profiles
+            write_output({'Fraction leaves in the sun, fraction of observed, fraction of observed&visible per layer'}, {'rows: simulations or time steps, columns: layer numbers'}, ...
+                         fnames.gap_file, n_col.gap, ns, true);
+            write_output({'aPAR per leaf layer'}, {'umol m-2 s-1'}, ...
+                         fnames.layer_aPAR_file, n_col.layer_aPAR, ns, true);
+            write_output({'aPAR by Cab per leaf layer'}, {'umol m-2 s-1'}, ...
+                         fnames.layer_aPAR_Cab_file, n_col.layer_aPAR_Cab, ns, true);
+            if options.calc_ebal
+                write_output({'leaf temperature of sunlit leaves, shaded leaves, and weighted average leaf temperature per layer'}, {'^oC ^oC ^oC'}, ...
+                             fnames.leaftemp_file, n_col.leaftemp, ns, true);
+                write_output({'sensible heat flux per layer'}, {'W m-2'}, ...
+                             fnames.layer_H_file, n_col.layer_H, ns, true);
+                write_output({'latent heat flux per layer'}, {'W m-2'}, ...
+                             fnames.layer_LE_file, n_col.layer_LE, ns, true);
+                write_output({'photosynthesis per layer'}, {'umol m-2 s-1'}, ...
+                             fnames.layer_A_file, n_col.layer_A, ns, true);
+                write_output({'average NPQ = 1-(fm-fo)/(fm0-fo0), per layer'}, {''}, ...
+                             fnames.layer_NPQ_file, n_col.layer_NPQ, ns, true);
+                write_output({'net radiation per leaf layer'}, {'W m-2'}, ...
+                             fnames.layer_Rn_file, n_col.layer_Rn, ns, true);
+            end
+            if options.calc_fluor
+                fluorescence_names = {'supward fluorescence per layer'};
+                fluorescence_units = {'W m-2'};
+                write_output(fluorescence_names, fluorescence_units, fnames.layer_fluorescence_file, n_col.layer_fluorescence, ns, true);
+            end
+        end
+        if options.calc_fluor
+            write_output({'hemispherically integrated fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ...
+                         fnames.fluorescence_hemis_file, n_col.fluorescence_hemis, ns, true);
+            write_output({'total emitted fluorescence by all leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ...
+                         fnames.fluorescence_emitted_by_all_leaves_file, n_col.fluorescence_emitted_by_all_leaves, ns, true);
+            write_output({'total emitted fluorescence by all photosystems for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ...
+                         fnames.fluorescence_emitted_by_all_photosystems_file, n_col.fluorescence_emitted_by_all_photosystems, ns, true);
+            write_output({'TOC fluorescence contribution from sunlit leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
+                         fnames.fluorescence_sunlit_file, n_col.fluorescence_sunlit, ns, true);
+            write_output({'TOC fluorescence contribution from shaded leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
+                         fnames.fluorescence_shaded_file, n_col.fluorescence_shaded, ns, true);
+            write_output({'TOC fluorescence contribution from from leaves and soil after scattering for wavelenghts of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
+                         fnames.fluorescence_scattered_file, n_col.fluorescence_scattered, ns, true);
+        end
+        write_output({'Bottom of canopy irradiance in the shaded fraction, and average BOC irradiance'}, {'First 2162 columns: shaded fraction. Last 2162 columns: average BOC irradiance. Unit: Wm-2 um-1'}, ...
+                     fnames.BOC_irradiance_file, n_col.BOC_irradiance, ns, true);
+    end
 
     fclose('all');
 
@@ -177,7 +185,6 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings)
     for k = 1:numel(fn)
         delete(fnames.(fn{k}));
     end
-
 end
 
 function write_output(header, units, bin_path, f_n_col, ns, not_header)
diff --git a/src/+io/output_data_binary.m b/src/+io/output_data_binary.m
@@ -10,7 +10,7 @@
     %% Fluxes product
     flu_out = [k iter.counter xyt.year(k) xyt.t(k) fluxes.Rntot fluxes.lEtot fluxes.Htot fluxes.Rnctot fluxes.lEctot, ...
                fluxes.Hctot fluxes.Actot fluxes.Rnstot fluxes.lEstot fluxes.Hstot fluxes.Gtot fluxes.Resp 1E6 * fluxes.aPAR, ...
-               1E6 * fluxes.aPAR_Cab fluxes.aPAR / rad.PAR fluxes.aPAR_Wm2 1E6 * rad.PAR rad.Eoutf rad.Eoutf ./ fluxes.aPAR_Wm2 Trap(k) * 10 Evap(k) * 10 Trap(k) * 10 + Evap(k) * 10 fluxes.GPP fluxes.NEE];
+               1E6 * fluxes.aPAR_Cab fluxes.aPAR / rad.PAR fluxes.aPAR_Wm2 1E6 * rad.PAR rad.Eoutf rad.Eoutf ./ fluxes.aPAR_Wm2 Trap * 10 Evap * 10 Trap * 10 + Evap * 10 fluxes.GPP fluxes.NEE];
     n_col.flu = length(flu_out);
     fwrite(f.flu_file, flu_out, 'double');
 
diff --git a/src/+io/read_config.m b/src/+io/read_config.m
@@ -1,4 +1,4 @@
-function [InputPath, OutputPath, InitialConditionPath] = read_config(config_file)
+function [InputPath, OutputPath, InitialConditionPath, FullCSVfiles] = read_config(config_file)
 
     file_id = fopen(config_file);
     config = textscan(file_id, '%s %s', 'HeaderLines', 0, 'Delimiter', '=');
@@ -17,3 +17,10 @@
 
     indx = find(strcmp(config_vars, 'InitialConditionPath'));
     InitialConditionPath = config_paths{indx};
+
+    indx = find(strcmp(config_vars, 'FullCSVfiles'));
+    if isempty(indx)
+        FullCSVfiles = 1;
+    else
+        FullCSVfiles = str2num(config_paths{indx});
+    end
diff --git a/src/+soilmoisture/calculateBoundaryConditions.m b/src/+soilmoisture/calculateBoundaryConditions.m
@@ -108,7 +108,7 @@
             C4(n - 1, 2) = 0;
             C4_a(n - 1) = 0;
         else
-            RHS(n) = RHS(n) + AVAIL0 - Evap(KT);
+            RHS(n) = RHS(n) + AVAIL0 - Evap;
         end
     end
     HeatMatrices.C4 = C4;
diff --git a/src/+soilmoisture/calculateEvapotranspiration.m b/src/+soilmoisture/calculateEvapotranspiration.m
@@ -1,4 +1,4 @@
-function [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, RWUs, RWUg] = calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt, ModelSettings, GroundwaterSettings)
+function [Rn_SOIL, Evap, Trap, r_a_SOIL, Srt, RWUs, RWUg] = calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt, ModelSettings, GroundwaterSettings)
 
     if ~GroundwaterSettings.GroundwaterCoupling  % no Groundwater coupling, added by Mostafa
         indxBotm = 1; % index of bottom layer is 1, STEMMUS calculates from bottom to top
@@ -14,22 +14,18 @@
     U = 100 .* (ForcingData.WS_msr(KT));
     AR = soilmoisture.calculateAerodynamicResistance(U);
     r_a_SOIL = AR.soil;
-
     Ta = ForcingData.Ta_msr(KT);
-    Evap = InitialValues.Evap;
-    EVAP = InitialValues.EVAP;
 
     if fluxes.lEctot < 1000 && fluxes.lEstot < 800 && fluxes.lEctot > -300 && fluxes.lEstot > -300 && any(SoilVariables.TT > 5)
         lambda1      = (2.501 - 0.002361 * Ta) * 1E6;
         lambda2      = (2.501 - 0.002361 * SoilVariables.Tss(KT)) * 1E6;
-        Evap(KT) = fluxes.lEstot / lambda2 * 0.1; % transfer to second value unit: cm s-1
-        EVAP(KT, 1) = Evap(KT);
-        Tp_t = fluxes.lEctot / lambda1 * 0.1; % transfer to second value
+        Evap = fluxes.lEstot / lambda2 * 0.1; % transfer to second value unit: cm s-1
+        Trap = fluxes.lEctot / lambda1 * 0.1; % transfer to second value
         Srt1 = RWU * 100 ./ ModelSettings.DeltZ';
     else
-        Evap(KT) = 0; % transfer to second value unit: cm s-1
-        EVAP(KT, 1) = Evap(KT);
-        Tp_t = 0; % transfer to second value
+        Evap = 0; % transfer to second value unit: cm s-1
+        Trap = 0; % transfer to second value
+        RWU = 0 * RWU;
         Srt1 = 0 ./ ModelSettings.DeltZ';
     end
 
@@ -38,7 +34,6 @@
             Srt(i, j) = Srt1(i, 1);
         end
     end
-    Trap(KT) = Tp_t;   % root water uptake integration by ModelSettings.DeltZ;
 
     % Calculate root water uptake from soil water (RWUs) and groundwater (RWUg)
     RWU = RWU * 100; % unit conversion from m/sec to cm/sec
diff --git a/src/+soilmoisture/solveSoilMoistureBalance.m b/src/+soilmoisture/solveSoilMoistureBalance.m
diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m
diff --git a/src/STEMMUS_SCOPE_exe.m b/src/STEMMUS_SCOPE_exe.m
