Redefine Evap and Trap arrays

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

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);
 

src/+io/bin_to_csv.m

@@ -91,6 +91,10 @@ 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);
 
+    % Comment unnecessary large size files
+    %{
+    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);
     %% Spectrum (added on 19 September 2008)
     spectrum_hemis_optical_names = {'hemispherically integrated radiation spectrum'};
     spectrum_hemis_optical_units = {'W m-2 um-1'};
@@ -112,15 +116,17 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings)
     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({'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); % comment unnecessary large size files
         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'}, ...
@@ -145,6 +151,7 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings)
             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,6 +161,9 @@ 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
+
+        % Comment unnecessary large size files
+        %{
         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'}, ...
@@ -166,9 +176,8 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings)
                      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);
 
     fclose('all');
 

src/+io/create_output_files_binary.m

@@ -20,11 +20,7 @@
     fnames.radiation_file = fullfile(Output_dir, 'radiation.bin');
     fnames.fluorescence_file = fullfile(Output_dir, 'fluorescence.bin');
     fnames.wl_file = fullfile(Output_dir, 'wl.bin');  % wavelength
-    fnames.irradiance_spectra_file = fullfile(Output_dir, 'irradiance_spectra.bin');  % Fluorescence
-    fnames.spectrum_hemis_optical_file = fullfile(Output_dir, 'spectrum_hemis.bin');
-    fnames.spectrum_obsdir_optical_file = fullfile(Output_dir, 'spectrum_obsdir.bin');
     fnames.reflectance_file = fullfile(Output_dir, 'reflectance.bin');  % reflectance spectrum
-    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');
@@ -38,10 +34,23 @@
     fnames.Sim_qva_file = fullfile(Output_dir, 'qva.bin'); % vapour flux due to dry air pressure gradient
     fnames.Sim_qtot_file = fullfile(Output_dir, 'qtot.bin'); % total flux (liquid + vapour)
 
+    % Comment unnecessary large size files
+    %{
+    fnames.BOC_irradiance_file = fullfile(Output_dir, 'BOC_irradiance.bin');  % reflectance spectrum
+    fnames.irradiance_spectra_file = fullfile(Output_dir, 'irradiance_spectra.bin');  % Fluorescence
+    fnames.spectrum_hemis_optical_file = fullfile(Output_dir, 'spectrum_hemis.bin');
+    fnames.spectrum_obsdir_optical_file = fullfile(Output_dir, 'spectrum_obsdir.bin');
+
     if options.calc_ebal
         fnames.spectrum_obsdir_BlackBody_file = fullfile(Output_dir, 'spectrum_obsdir_BlackBody.bin');  % spectrum observation direction
     end
 
+    if options.calc_planck && options.calc_ebal
+        fnames.spectrum_obsdir_thermal_file    = fullfile(Output_dir, 'spectrum_obsdir_thermal.bin');  % spectrum observation direction
+        fnames.spectrum_hemis_thermal_file     = fullfile(Output_dir, 'spectrum_hemis_thermal.bin');   % spectrum hemispherically integrated
+    end
+    %}
+
     % if ~(options.simulation==1)
     fnames.pars_and_input_file           = fullfile(Output_dir, 'pars_and_input.bin');      % wavelength
 
@@ -60,7 +69,7 @@
     %
     %% Optional Output
     if options.calc_vert_profiles
-        fnames.gap_file            = fullfile(Output_dir, 'gap.bin');              % gap
+        % fnames.gap_file            = fullfile(Output_dir, 'gap.bin');            % gap, comment unnecessary large size files
         fnames.leaftemp_file        = fullfile(Output_dir, 'leaftemp.bin');        % leaftemp
         fnames.layer_H_file         = fullfile(Output_dir, 'layer_H.bin');      % vertical profile
         fnames.layer_LE_file        = fullfile(Output_dir, 'layer_lE.bin');         % latent heat
@@ -89,12 +98,15 @@
             fnames.fluorescencePSI_file       = fullfile(Output_dir, 'fluorescencePSI.bin');     % Fluorescence
             fnames.fluorescencePSII_file      = fullfile(Output_dir, 'fluorescencePSII.bin');     % Fluorescence
         end
+        % Comment unnecessary large size files
+        %{
         fnames.fluorescence_hemis_file        = fullfile(Output_dir, 'fluorescence_hemis.bin');     % Fluorescence
         fnames.fluorescence_emitted_by_all_leaves_file         = fullfile(Output_dir, 'fluorescence_emitted_by_all_leaves.bin');
         fnames.fluorescence_emitted_by_all_photosystems_file   = fullfile(Output_dir, 'fluorescence_emitted_by_all_photosystems.bin');
         fnames.fluorescence_sunlit_file       = fullfile(Output_dir, 'fluorescence_sunlit.bin');     % Fluorescence
         fnames.fluorescence_shaded_file       = fullfile(Output_dir, 'fluorescence_shaded.bin');     % Fluorescence
         fnames.fluorescence_scattered_file    = fullfile(Output_dir, 'fluorescence_scattered.bin');      % Fluorescence
+        %}
     else
         delete([Output_dir, 'fluorescence.bin']);
     end
@@ -103,11 +115,6 @@
         delete([Output_dir, 'BRDF/*.bin']);
     end
 
-    if options.calc_planck && options.calc_ebal
-        fnames.spectrum_obsdir_thermal_file    = fullfile(Output_dir, 'spectrum_obsdir_thermal.bin');  % spectrum observation direction
-        fnames.spectrum_hemis_thermal_file     = fullfile(Output_dir, 'spectrum_hemis_thermal.bin');   % spectrum hemispherically integrated
-    end
-
     % Create empty files for appending
     structfun(@(x) fopen(x, 'w'), fnames, 'UniformOutput', false);
     fclose("all");

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');
 
@@ -83,6 +83,8 @@
     n_col.Sim_qtot = length(Sim_qtot_out);
     fwrite(f.Sim_qtot_file, Sim_qtot_out, 'double');
 
+    % Comment unnecessary large size files
+    %{
     %% Spectrum (added on 19 September 2008)
     spectrum_hemis_optical_out =  rad.Eout_;
     n_col.spectrum_hemis_optical = length(spectrum_hemis_optical_out);
@@ -115,6 +117,11 @@
     n_col.irradiance_spectra = length(irradiance_spectra_out);
     fwrite(f.irradiance_spectra_file, irradiance_spectra_out, 'double');
 
+    BOC_irradiance_out  =  [rad.Emin_(canopy.nlayers + 1, :), rad.Emin_(canopy.nlayers + 1, :) + (rad.Esun_ * gap.Ps(canopy.nlayers + 1)')'];
+    n_col.BOC_irradiance = length(BOC_irradiance_out);
+    fwrite(f.BOC_irradiance_file, BOC_irradiance_out, 'double');
+    %}
+
     reflectance = pi * rad.Lo_ ./ (rad.Esun_ + rad.Esky_);
     reflectance(spectral.wlS > 3000) = NaN;
     reflectance_out =  [reflectance'];
@@ -140,9 +147,9 @@
     if options.calc_vert_profiles
 
         % gap
-        gap_out   =  [gap.Ps gap.Po gap.Pso];
-        n_col.gap = numel(gap_out(:));
-        fwrite(f.gap_file, gap_out, 'double');
+        % gap_out   =  [gap.Ps gap.Po gap.Pso];
+        % n_col.gap = numel(gap_out(:));
+        % fwrite(f.gap_file, gap_out, 'double'); % comment unnecessary large size files
 
         layer_aPAR_out   =  [1E6 * profiles.Pn1d' 0];
         n_col.layer_aPAR = length(layer_aPAR_out);
@@ -202,7 +209,8 @@
                 n_col.fluorescencePSII = length(fluorescencePSII_out);
                 fwrite(f.fluorescencePSII_file, fluorescencePSII_out, 'double');
             end
-
+            % Comment unnecessary large size files
+            %{
             fluorescence_hemis_out  =  [rad.Fhem_];
             n_col.fluorescence_hemis = length(fluorescence_hemis_out);
             fwrite(f.fluorescence_hemis_file, fluorescence_hemis_out, 'double');
@@ -226,12 +234,9 @@
             fluorescence_scattered_out  =  [sum(rad.LoF_scattered, 2) + sum(rad.LoF_soil, 2)];
             n_col.fluorescence_scattered = length(fluorescence_scattered_out);
             fwrite(f.fluorescence_scattered_file, fluorescence_scattered_out, 'double');
-
+    %}
         end
     end
-    BOC_irradiance_out  =  [rad.Emin_(canopy.nlayers + 1, :), rad.Emin_(canopy.nlayers + 1, :) + (rad.Esun_ * gap.Ps(canopy.nlayers + 1)')'];
-    n_col.BOC_irradiance = length(BOC_irradiance_out);
-    fwrite(f.BOC_irradiance_file, BOC_irradiance_out, 'double');
 
     %%
     if options.calc_directional && options.calc_ebal

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;

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

src/+soilmoisture/solveSoilMoistureBalance.m

@@ -1,5 +1,5 @@
-function [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip, RWUs, RWUg, ForcingData] = solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, ...
-                                                                                                                                                                                        BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, ModelSettings, GroundwaterSettings)
+function [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip, RWUs, RWUg, ForcingData] = solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, ...
+                                                                                                                                                                                  BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, ModelSettings, GroundwaterSettings)
     %{
         Solve the soil moisture balance equation with the Thomas algorithm to
         update the soil matric potential 'hh', the finite difference
@@ -13,14 +13,13 @@
 
     [RHS, HeatMatrices, boundaryFluxArray] = soilmoisture.calculateTimeDerivatives(InitialValues, SoilVariables, HeatMatrices, Delt_t, P_gg, ModelSettings, GroundwaterSettings);
 
-    % When BoundaryCondition.NBCh == 3, otherwise Rn_SOIL, Evap, EVAP, Trap,
+    % When BoundaryCondition.NBCh == 3, otherwise Rn_SOIL, Evap, Trap,
     % r_a_SOIL, Srt will be empty arrays! see issue 98, item 2
     if BoundaryCondition.NBCh == 3
-        [Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, RWUs, RWUg] = soilmoisture.calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt, ModelSettings, GroundwaterSettings);
+        [Rn_SOIL, Evap, Trap, r_a_SOIL, Srt, RWUs, RWUg] = soilmoisture.calculateEvapotranspiration(InitialValues, ForcingData, SoilVariables, KT, RWU, fluxes, Srt, ModelSettings, GroundwaterSettings);
     else
         Rn_SOIL = InitialValues.Rn_SOIL;
-        Evap = InitialValues.Evap;
-        EVAP = InitialValues.EVAP;
+        Evap = [];
         Trap = [];
         r_a_SOIL = InitialValues.r_a_SOIL;
     end