diff --git a/CHANGELOG.md b/CHANGELOG.md index 20606166..77b23d32 100644 --- 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) diff --git a/docs/getting_started.md b/docs/getting_started.md index ac14bd09..2ce20494 100644 --- 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 index 5bf5133e..fff99aa3 100755 Binary files a/run_model_on_snellius/exe/STEMMUS_SCOPE and b/run_model_on_snellius/exe/STEMMUS_SCOPE differ diff --git a/src/+energy/calculateBoundaryConditions.m b/src/+energy/calculateBoundaryConditions.m index 4e59856c..ee68460b 100644 --- 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 index 77ca0bf6..99fea87b 100644 --- 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 index 8a4baf89..5bc13136 100644 --- 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 index e560d195..d5a918d4 100644 --- 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 index 587f4e99..60f33185 100644 --- 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 index 960c98d4..cd1bc817 100644 --- 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 index 450a5d5d..c865ef66 100644 --- 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 index 8d9a6bde..52ee8544 100644 --- a/src/+soilmoisture/solveSoilMoistureBalance.m +++ b/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 diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 8c6e155a..1db95466 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -39,7 +39,7 @@ if strcmp(bmiMode, "initialize") || strcmp(runMode, "full") % Read the configPath file. Due to using MATLAB compiler, we cannot use run(CFG) disp (['Reading config from ', CFG]); - [InputPath, OutputPath, InitialConditionPath] = io.read_config(CFG); + [InputPath, OutputPath, InitialConditionPath, FullCSVfiles] = io.read_config(CFG); % Prepare forcing and soil data [SiteProperties, SoilProperties, TimeProperties] = io.prepareInputData(InputPath); @@ -80,8 +80,8 @@ SAVE = InitialValues.SAVE; P_gOLD = InitialValues.P_gOLD; gwfluxes = struct(); - Evap = InitialValues.Evap; - EVAP = InitialValues.EVAP; + Evap = []; + Trap = []; RWUs = []; RWUg = []; @@ -653,7 +653,7 @@ GasDispersivity = conductivity.calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g, ModelSettings); % Srt is both input and output - [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, EVAP, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip, RWUs, RWUg, ForcingData] = ... + [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip, RWUs, RWUg, ForcingData] = ... soilmoisture.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); @@ -670,7 +670,7 @@ DSTOR = min(EXCESS, DSTMAX); RS = (EXCESS - DSTOR) / Delt_t; else - AVAIL = AVAIL0 - Evap(KT); + AVAIL = AVAIL0 - Evap; EXCESS = (AVAIL + HBoundaryFlux.QMT) * Delt_t; % (unit converstion from cm/sec to cm/30mins) if abs(EXCESS / Delt_t) <= 1e-10 EXCESS = 0; @@ -701,7 +701,7 @@ AirVariabes, VaporVariables, GasDispersivity, ThermalConductivityCapacity, ... HBoundaryFlux, BoundaryCondition, ForcingData, DRHOVh, DRHOVT, KL_T, ... Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, Delt_t, P_g, P_gg, ... - TOLD, Precip, EVAP, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings); + TOLD, Precip, Evap, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings); end if max(CHK) < 0.1 @@ -822,7 +822,7 @@ io.output_verification(Output_dir); end - io.bin_to_csv(fnames, n_col, k, options, SoilLayer, GroundwaterSettings); + io.bin_to_csv(fnames, n_col, k, options, SoilLayer, GroundwaterSettings, FullCSVfiles); save([Output_dir, 'output.mat']); end diff --git a/src/STEMMUS_SCOPE_exe.m b/src/STEMMUS_SCOPE_exe.m index 1f966975..c090d08c 100644 --- a/src/STEMMUS_SCOPE_exe.m +++ b/src/STEMMUS_SCOPE_exe.m @@ -26,12 +26,13 @@ function STEMMUS_SCOPE_exe(config_file, runMode) 'fluxes', ... % Atmospheric fluxes 'TT', ... % Soil temperature over depth 'SoilVariables', ... % Structure that includes different variables of soil moisture - 'GroundwaterSettings' ... % groundwater settings including input data from MODFLOW - 'gwfluxes' ... % structure that includes groundwater recharge and its individual components - 'EVAP' ... % evaporation - 'RWUs' ... % soil water root uptake - 'RWUg' ... % groundwater root uptake - 'ForcingData' ... % forcing data that includes Dunnian runoff and Hortonian runoff + 'GroundwaterSettings', ... % groundwater settings including input data from MODFLOW + 'gwfluxes', ... % structure that includes groundwater recharge and its individual components + 'Evap', ... % evaporation + 'Trap', ... % transpiration + 'RWUs', ... % soil water root uptake + 'RWUg', ... % groundwater root uptake + 'ForcingData', ... % forcing data that includes Dunnian runoff and Hortonian runoff 'RS' ... % total surface runoff }; %#ok