diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE index c924d75a..8195625f 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/+init/applySoilHeteroEffect.m b/src/+init/applySoilHeteroEffect.m index 35fe265e..1a85a68b 100644 --- a/src/+init/applySoilHeteroEffect.m +++ b/src/+init/applySoilHeteroEffect.m @@ -1,7 +1,7 @@ -function [SoilVariables, VanGenuchten] = applySoilHeteroEffect(SoilProperties, SoilConstants, SoilData, SoilVariables, VanGenuchten) +function [SoilVariables, VanGenuchten] = applySoilHeteroEffect(SoilProperties, SoilVariables, VanGenuchten) - initX = SoilData.InitialValues.initX; - initND = SoilData.InitialValues.initND; + initX = SoilVariables.InitialValues.initX; + initND = SoilVariables.InitialValues.initND; SoilVariables.Phi_s = []; % see issue 139, item 5 SoilVariables.Lamda = []; % see issue 139, item 5 @@ -11,6 +11,9 @@ ModelSettings = io.getModelSettings(); Eqlspace = ModelSettings.Eqlspace; + % Get soil constants for StartInit + SoilConstants = io.getSoilConstants(); + for i = 1:6 if ModelSettings.SWCC == 0 if SoilConstants.CHST == 0 @@ -23,7 +26,7 @@ if ~Eqlspace j = ModelSettings.J; for i = 1:length(initX) - SoilData.InitialValues.initH(i) = init.calcInitH(VanGenuchten.Theta_s(j), VanGenuchten.Theta_r(j), initX(i), VanGenuchten.n(j), VanGenuchten.m(j), VanGenuchten.Alpha(j)); + SoilVariables.InitialValues.initH(i) = init.calcInitH(VanGenuchten.Theta_s(j), VanGenuchten.Theta_r(j), initX(i), VanGenuchten.n(j), VanGenuchten.m(j), VanGenuchten.Alpha(j)); end Dmark = []; for i = 1:ModelSettings.NL @@ -31,15 +34,15 @@ InitLnth(i) = ModelSettings.Tot_Depth - SoilConstants.Elmn_Lnth; for subRoutine = 5:-1:1 if abs(InitLnth(i) - initND(subRoutine)) < 1e-10 - [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilConstants, SoilData, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i); - SoilData.InitialValues.initH = initH; + [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i); + SoilVariables.InitialValues.initH = initH; Dmark = i + 2; end end if abs(InitLnth(i)) < 1e-10 subRoutine = 0; - [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilConstants, SoilData, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i); - SoilData.InitialValues.initH = initH; + [SoilVariables, VanGenuchten, initH] = init.soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, i); + SoilVariables.InitialValues.initH = initH; end end else diff --git a/src/+init/applySoilHeteroWithInitialFreezing.m b/src/+init/applySoilHeteroWithInitialFreezing.m index 2c1a5d0d..c77965ab 100644 --- a/src/+init/applySoilHeteroWithInitialFreezing.m +++ b/src/+init/applySoilHeteroWithInitialFreezing.m @@ -1,4 +1,4 @@ -function [SoilVariables] = applySoilHeteroWithInitialFreezing(SoilConstants, SoilVariables) +function [SoilVariables] = applySoilHeteroWithInitialFreezing(LatentHeatOfFreezing, SoilVariables) SoilVariables.ISFT = 0; @@ -6,21 +6,21 @@ ModelSettings = io.getModelSettings(); for i = 1:ModelSettings.NN - SoilVariables.h_frez = init.updateHfreez(i, SoilVariables, SoilConstants); + SoilVariables.h_frez = init.updateHfreez(i, LatentHeatOfFreezing, SoilVariables); SoilVariables.hh_frez(i) = SoilVariables.h_frez(i); SoilVariables.h(i) = SoilVariables.h(i) - SoilVariables.h_frez(i); SoilVariables.hh(i) = SoilVariables.h(i); SoilVariables.SAVEh(i) = SoilVariables.h(i); SoilVariables.SAVEhh(i) = SoilVariables.hh(i); - [SoilVariables.Gama_h, SoilVariables.Gama_hh] = init.updateGamaH(i, SoilConstants, SoilVariables); + [SoilVariables.Gama_h, SoilVariables.Gama_hh] = init.updateGamaH(i, SoilVariables); if ModelSettings.Thmrlefc == 1 SoilVariables.TT(i) = SoilVariables.T(i); end if ModelSettings.Soilairefc == 1 - SoilConstants.P_g(i) = 95197.850; - SoilConstants.P_gg(i) = SoilConstants.P_g(i); + SoilVariables.P_g(i) = 95197.850; + SoilVariables.P_gg(i) = SoilVariables.P_g(i); end if i < ModelSettings.NN SoilVariables.XWRE(i, 1) = 0; diff --git a/src/+init/calculateInitialThermal.m b/src/+init/calculateInitialThermal.m index 556b3404..3406d893 100644 --- a/src/+init/calculateInitialThermal.m +++ b/src/+init/calculateInitialThermal.m @@ -1,4 +1,4 @@ -function ThermalConductivity = calculateInitialThermal(SoilConstants, SoilVariables, VanGenuchten) +function ThermalConductivity = calculateInitialThermal(SoilVariables, VanGenuchten) FEHCAP = []; % see issue 139 HCAP(1) = 0.998 * 4.182; diff --git a/src/+init/defineSoilVariables.m b/src/+init/defineSoilVariables.m new file mode 100644 index 00000000..64828842 --- /dev/null +++ b/src/+init/defineSoilVariables.m @@ -0,0 +1,38 @@ +function SoilVariables = defineSoilVariables(InitialValues, SoilProperties, VanGenuchten) + %{ + Create a SoilVariables structure and define some variables. + SoilVariables is input for initializing Temperature, Matric potential + and soil air pressure. + %} + + % Create a SoilVariables structure + SoilVariables = struct(); + + % Add initial values to SoilVariables + soil_fields = { + 'P_g', 'P_gg', 'h', 'T', 'TT', 'h_frez', 'Theta_L', ... + 'Theta_LL', 'Theta_V', 'Theta_g', 'Se', 'KL_h', ... + 'DTheta_LLh', 'KfL_T', 'Theta_II', 'Theta_I', ... + 'Theta_UU', 'Theta_U', 'TT_CRIT', 'KfL_h', 'DTheta_UUh' + }; + for field = soil_fields + SoilVariables.(field{1}) = InitialValues.(field{1}); + end + + % Calculate some of SoilVariables + SoilVariables.XK = VanGenuchten.Theta_r + 0.02; % 0.11 This is for silt loam; For sand XK=0.025 + SoilVariables.XWILT = equations.van_genuchten(VanGenuchten.Theta_s, VanGenuchten.Theta_r, VanGenuchten.Alpha, -1.5e4, VanGenuchten.n, VanGenuchten.m); + SoilVariables.XCAP = equations.van_genuchten(VanGenuchten.Theta_s, VanGenuchten.Theta_r, VanGenuchten.Alpha, -336, VanGenuchten.n, VanGenuchten.m); + + SoilVariables.VPERSOC = equations.calc_msoc_fraction(SoilProperties.MSOC); % fraction of soil organic matter + SoilVariables.FOSL = 1 - SoilProperties.FOC - SoilProperties.FOS - SoilVariables.VPERSOC; + + % Add variables to SoilVariables, see issue 141 + SoilVariables.POR = SoilProperties.porosity; + SoilVariables.VPERS = SoilProperties.FOS' .* (1 - SoilVariables.POR); + SoilVariables.VPERSL = SoilVariables.FOSL' .* (1 - SoilVariables.POR); + SoilVariables.VPERC = SoilProperties.FOC' .* (1 - SoilVariables.POR); + SoilVariables.Ks = SoilProperties.SaturatedK; + SoilVariables.Vol_qtz = SoilProperties.FOS; + +end diff --git a/src/+init/setBoundaryCondition.m b/src/+init/setBoundaryCondition.m index 2d517449..89dabde1 100644 --- a/src/+init/setBoundaryCondition.m +++ b/src/+init/setBoundaryCondition.m @@ -1,4 +1,4 @@ -function BoundaryCondition = setBoundaryCondition(SoilVariables, SoilConstants, initialLandcoverClass) +function BoundaryCondition = setBoundaryCondition(SoilVariables, ForcingData, initialLandcoverClass) % Variables information for boundary condition settings % NBCh Indicator for type of surface boundary condition on mass euqation to be applied; % "1"--Specified matric head; @@ -47,7 +47,7 @@ IRPT1 = 0; IRPT2 = 0; - Ta_msr = SoilConstants.Ta_msr; + Ta_msr = ForcingData.Ta_msr; % get model settings ModelSettings = io.getModelSettings(); diff --git a/src/+init/setSoilConstants.m b/src/+init/setSoilConstants.m deleted file mode 100644 index 60bc8d79..00000000 --- a/src/+init/setSoilConstants.m +++ /dev/null @@ -1,31 +0,0 @@ -function SoilConstants = setSoilConstants(InitialValues, SoilProperties, ForcingData) - - % create SoilConstants structure holding initial values - SoilConstants = struct(); - soil_fields = { - 'P_g', 'P_gg', 'h', 'T', 'TT', 'h_frez', 'Theta_L', ... - 'Theta_LL', 'Theta_V', 'Theta_g', 'Se', 'KL_h', ... - 'DTheta_LLh', 'KfL_T', 'Theta_II', 'Theta_I', ... - 'Theta_UU', 'Theta_U', 'TT_CRIT', 'KfL_h', 'DTheta_UUh' - }; - for field = soil_fields - SoilConstants.(field{1}) = InitialValues.(field{1}); - end - - SoilConstants.Ta_msr = ForcingData.Ta_msr; - - SoilConstants.hd = -1e7; - SoilConstants.hm = -9899; - SoilConstants.CHST = 0; - SoilConstants.Elmn_Lnth = 0; - SoilConstants.Dmark = 0; - SoilConstants.Vol_qtz = SoilProperties.FOS; - SoilConstants.VPERSOC = equations.calc_msoc_fraction(SoilProperties.MSOC); % fraction of soil organic matter - SoilConstants.FOSL = 1 - SoilProperties.FOC - SoilProperties.FOS - SoilConstants.VPERSOC; - SoilConstants.Phi_S = [-17.9 -17 -17 -19 -10 -10]; - SoilConstants.Phi_soc = -0.0103; - SoilConstants.Lamda_soc = 2.7; - SoilConstants.Theta_soc = 0.6; - SoilConstants.XK = 0.11; % 0.11 This is for silt loam; For sand XK=0.025 - -end diff --git a/src/+init/setSoilVariables.m b/src/+init/setSoilVariables.m deleted file mode 100644 index c25e2a3b..00000000 --- a/src/+init/setSoilVariables.m +++ /dev/null @@ -1,21 +0,0 @@ -function SoilVariables = setSoilVariables(SoilProperties, SoilConstants, VanGenuchten) - - % SoilVariables.XK = 0.09; %0.11 This is for silt loam; For sand XK=0.025 - SoilVariables.XK = VanGenuchten.Theta_r + 0.02; % 0.11 This is for silt loam; For sand XK=0.025 - SoilVariables.XWILT = equations.van_genuchten(VanGenuchten.Theta_s, VanGenuchten.Theta_r, VanGenuchten.Alpha, -1.5e4, VanGenuchten.n, VanGenuchten.m); - SoilVariables.XCAP = equations.van_genuchten(VanGenuchten.Theta_s, VanGenuchten.Theta_r, VanGenuchten.Alpha, -336, VanGenuchten.n, VanGenuchten.m); - - SoilVariables.POR = SoilProperties.porosity; - SoilVariables.VPERS = SoilProperties.FOS' .* (1 - SoilVariables.POR); - SoilVariables.VPERSL = SoilConstants.FOSL' .* (1 - SoilVariables.POR); - SoilVariables.VPERC = SoilProperties.FOC' .* (1 - SoilVariables.POR); - SoilVariables.VPERSOC = SoilConstants.VPERSOC; - - SoilVariables.Ks = SoilProperties.SaturatedK; - - SoilVariables.h = SoilConstants.h; - SoilVariables.T = SoilConstants.T; - SoilVariables.TT = SoilConstants.TT; - SoilVariables.h_frez = SoilConstants.h_frez; - -end diff --git a/src/+init/soilHeteroSubroutine.m b/src/+init/soilHeteroSubroutine.m index 63d020aa..f9d45937 100644 --- a/src/+init/soilHeteroSubroutine.m +++ b/src/+init/soilHeteroSubroutine.m @@ -1,11 +1,11 @@ -function [SoilVariables, VanGenuchten, initH] = soilHeteroSubroutine(subRoutine, SoilConstants, SoilData, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, ML) +function [SoilVariables, VanGenuchten, initH] = soilHeteroSubroutine(subRoutine, SoilProperties, SoilVariables, VanGenuchten, ImpedF, Dmark, ML) %{ considering soil heterogeneity effect %} - initX = SoilData.InitialValues.initX; - initT = SoilData.InitialValues.initT; - initH = SoilData.InitialValues.initH; + initX = SoilVariables.InitialValues.initX; + initT = SoilVariables.InitialValues.initT; + initH = SoilVariables.InitialValues.initH; % get model settings ModelSettings = io.getModelSettings(); @@ -56,14 +56,14 @@ end J = SoilVariables.IS(i); - [SoilVariables, VanGenuchten] = init.updateSoilVariables(SoilVariables, VanGenuchten, SoilConstants, SoilProperties, j, J); + [SoilVariables, VanGenuchten] = init.updateSoilVariables(SoilVariables, VanGenuchten, SoilProperties, j, J); SoilVariables.Imped(i) = ImpedF(J); - initH(indexOfInit) = init.updateInitH(initX(indexOfInit), VanGenuchten, SoilConstants, SoilVariables, j); + initH(indexOfInit) = init.updateInitH(initX(indexOfInit), VanGenuchten, SoilVariables, j); if subRoutine == 5 - Btmh = init.updateBtmh(VanGenuchten, SoilData, SoilVariables, i); - SoilVariables.T(i) = SoilData.BtmT + (i - 1) * (initT(indexOfInit) - SoilData.BtmT) / ML; + Btmh = init.updateBtmh(VanGenuchten, SoilVariables, i); + SoilVariables.T(i) = SoilVariables.BtmT + (i - 1) * (initT(indexOfInit) - SoilVariables.BtmT) / ML; SoilVariables.h(i) = Btmh + (i - 1) * (initH(indexOfInit) - Btmh) / ML; SoilVariables.IH(i) = 1; % Index of wetting history of soil which would be assumed as dry at the first with the value of 1 else diff --git a/src/+init/updateBtmh.m b/src/+init/updateBtmh.m index ab82d23a..cb0c1182 100644 --- a/src/+init/updateBtmh.m +++ b/src/+init/updateBtmh.m @@ -1,10 +1,10 @@ -function Btmh = updateBtmh(VanGenuchten, SoilData, SoilVariables, i) +function Btmh = updateBtmh(VanGenuchten, SoilVariables, i) % get model settings ModelSettings = io.getModelSettings(); if ModelSettings.SWCC == 1 % VG soil water retention model - Btmh = init.calcInitH(VanGenuchten.Theta_s(i), VanGenuchten.Theta_r(i), SoilData.BtmX, VanGenuchten.n(i), VanGenuchten.m(i), VanGenuchten.Alpha(i)); + Btmh = init.calcInitH(VanGenuchten.Theta_s(i), VanGenuchten.Theta_r(i), SoilVariables.BtmX, VanGenuchten.n(i), VanGenuchten.m(i), VanGenuchten.Alpha(i)); else - Btmh = SoilVariables.Phi_s(i) * (SoilData.BtmX / VanGenuchten.Theta_s(i))^(-1 / SoilVariables.Lamda(i)); + Btmh = SoilVariables.Phi_s(i) * (SoilVariables.BtmX / VanGenuchten.Theta_s(i))^(-1 / SoilVariables.Lamda(i)); end end diff --git a/src/+init/updateGamaH.m b/src/+init/updateGamaH.m index f5be1bea..26002264 100644 --- a/src/+init/updateGamaH.m +++ b/src/+init/updateGamaH.m @@ -1,4 +1,7 @@ -function [Gama_h, Gama_hh] = updateGamaH(i, SoilConstants, SoilVariables) +function [Gama_h, Gama_hh] = updateGamaH(i, SoilVariables) + + % get soil constants + SoilConstants = io.getSoilConstants(); hd = SoilConstants.hd; hm = SoilConstants.hm; diff --git a/src/+init/updateHfreez.m b/src/+init/updateHfreez.m index 1d721cd9..0d64ef48 100644 --- a/src/+init/updateHfreez.m +++ b/src/+init/updateHfreez.m @@ -1,10 +1,7 @@ -function h_frez = updateHfreez(i, SoilVariables, SoilConstants) +function h_frez = updateHfreez(i, LatentHeatOfFreezing, SoilVariables) % get Constants Constants = io.define_constants(); - L_f = 3.34 * 1e5; % latent heat of freezing fusion J Kg-1 - T0 = 273.15; % unit K - T = SoilVariables.T; h_frez = SoilVariables.h_frez; h = SoilVariables.h; @@ -14,7 +11,7 @@ ModelSettings = io.getModelSettings(); if T(i) <= 0 - h_frez(i) = L_f * 1e4 * (T(i)) / Constants.g / T0; + h_frez(i) = LatentHeatOfFreezing * 1e4 * (T(i)) / Constants.g / ModelSettings.T0; else h_frez(i) = 0; end diff --git a/src/+init/updateInitH.m b/src/+init/updateInitH.m index e3f4ff4c..6c73c33a 100644 --- a/src/+init/updateInitH.m +++ b/src/+init/updateInitH.m @@ -1,4 +1,4 @@ -function initH = updateInitH(initX, VanGenuchten, SoilConstants, SoilVariables, i) +function initH = updateInitH(initX, VanGenuchten, SoilVariables, i) % get model settings ModelSettings = io.getModelSettings(); diff --git a/src/+init/updateSoilVariables.m b/src/+init/updateSoilVariables.m index 298f5b8a..9273812e 100644 --- a/src/+init/updateSoilVariables.m +++ b/src/+init/updateSoilVariables.m @@ -1,12 +1,15 @@ -function [SoilVariables, VanGenuchten] = updateSoilVariables(SoilVariables, VanGenuchten, SoilConstants, SoilProperties, i, j) +function [SoilVariables, VanGenuchten] = updateSoilVariables(SoilVariables, VanGenuchten, SoilProperties, i, j) + + % Get soil constants for StartInit + SoilConstants = io.getSoilConstants(); SoilVariables.POR(i) = SoilProperties.porosity(j); SoilVariables.Ks(i) = SoilProperties.SaturatedK(j); - SoilVariables.Theta_qtz(i) = SoilConstants.Vol_qtz(j); + SoilVariables.Theta_qtz(i) = SoilVariables.Vol_qtz(j); SoilVariables.VPER(i, 1) = SoilVariables.VPERS(j); SoilVariables.VPER(i, 2) = SoilVariables.VPERSL(j); SoilVariables.VPER(i, 3) = SoilVariables.VPERC(j); - SoilVariables.XSOC(i) = SoilConstants.VPERSOC(j); + SoilVariables.XSOC(i) = SoilVariables.VPERSOC(j); SoilVariables.XK(i) = SoilConstants.XK; % get model settings diff --git a/src/+io/getModelSettings.m b/src/+io/getModelSettings.m index ce022fe3..5560f754 100644 --- a/src/+io/getModelSettings.m +++ b/src/+io/getModelSettings.m @@ -59,7 +59,6 @@ ModelSettings.NS = 1; % Number of soil types; % The time and domain information setting - ModelSettings.KIT = 0; % KIT is used to count the number of iteration in a time step; ModelSettings.NIT = 30; % Desirable number of iterations in a time step; ModelSettings.KT = 0; % Number of time steps; diff --git a/src/+io/getSoilConstants.m b/src/+io/getSoilConstants.m new file mode 100644 index 00000000..f23ca38e --- /dev/null +++ b/src/+io/getSoilConstants.m @@ -0,0 +1,16 @@ +function SoilConstants = getSoilConstants() + + SoilConstants.hd = -1e7; + SoilConstants.hm = -9899; + SoilConstants.CHST = 0; + SoilConstants.Elmn_Lnth = 0; + SoilConstants.Dmark = 0; + SoilConstants.Phi_S = [-17.9 -17 -17 -19 -10 -10]; + SoilConstants.Phi_soc = -0.0103; + SoilConstants.Lamda_soc = 2.7; + SoilConstants.Theta_soc = 0.6; + % XK=0.11 for silt loam; For sand XK=0.025 + % SoilConstants.XK is used in updateSoilVariables + SoilConstants.XK = 0.11; + +end diff --git a/src/+io/loadSoilData.m b/src/+io/loadSoilInitialValues.m similarity index 73% rename from src/+io/loadSoilData.m rename to src/+io/loadSoilInitialValues.m index ce65532c..ebf04fdb 100644 --- a/src/+io/loadSoilData.m +++ b/src/+io/loadSoilInitialValues.m @@ -1,4 +1,4 @@ -function [SoilData] = loadSoilData(InputPath, TimeProperties, Tot_Depth, SoilProperties, ForcingData, SWCC) +function [InitialValues, BtmX, BtmT, Tss] = loadSoilInitialValues(InputPath, TimeProperties, SoilProperties, ForcingData) %{ Producing initial soil moisture and soil temperature profile %} @@ -8,6 +8,11 @@ Dur_tot = TimeProperties.Dur_tot; Ta_msr = ForcingData.Ta_msr; + % Get model settings + ModelSettings = io.getModelSettings(); + Tot_Depth = ModelSettings.Tot_Depth; + SWCC = ModelSettings.SWCC; + % load initial soil moisture and soil temperature from ERA5 load(fullfile(InputPath, "soil_init.mat")); % The following parameters are loaded: @@ -27,8 +32,8 @@ % InitX5 % InitX6 % BtmX - SoilData.BtmX = BtmX; - SoilData.Tss = Tss; + BtmX = BtmX; + Tss = Tss; InitND1 = 5; % Unit of it is cm. These variables are used to indicated the depth corresponding to the measurement. InitND2 = 15; @@ -44,7 +49,7 @@ InitT4 = 4.29; InitT5 = 3.657; InitT6 = 6.033; - SoilData.BtmT = 6.62; % 9 8.1 + BtmT = 6.62; % 9 8.1 InitX0 = 0.088; InitX1 = 0.095; % Measured soil liquid moisture content InitX2 = 0.180; % 0.169 @@ -52,7 +57,7 @@ InitX4 = 0.184; % 0.114 InitX5 = 0.0845; InitX6 = 0.078; - SoilData.BtmX = 0.078; % 0.078 0.05; % The initial moisture content at the bottom of the column. + BtmX = 0.078; % 0.078 0.05; % The initial moisture content at the bottom of the column. else if InitT0 < 0 || InitT1 < 0 || InitT2 < 0 || InitT3 < 0 || InitT4 < 0 || InitT5 < 0 || InitT6 < 0 InitT0 = 0; @@ -62,12 +67,12 @@ InitT4 = 0; InitT5 = 0; InitT6 = 0; - SoilData.Tss = InitT0; + Tss = InitT0; end if nanmean(Ta_msr) < 0 - SoilData.BtmT = 0; % 9 8.1 + BtmT = 0; % 9 8.1 else - SoilData.BtmT = nanmean(Ta_msr); + BtmT = nanmean(Ta_msr); end if InitX0 > SaturatedMC(1) || InitX1 > SaturatedMC(1) || InitX2 > SaturatedMC(2) || ... InitX3 > SaturatedMC(3) || InitX4 > SaturatedMC(4) || InitX5 > SaturatedMC(5) || InitX6 > SaturatedMC(6) @@ -78,12 +83,12 @@ InitX4 = fieldMC(4); % 0.14335 InitX5 = fieldMC(5); InitX6 = fieldMC(6); - SoilData.BtmX = fieldMC(6); + BtmX = fieldMC(6); end end - SoilData.InitialValues.initX = [InitX0, InitX1, InitX2, InitX3, InitX4, InitX5, InitX6]; - SoilData.InitialValues.initND = [InitND1, InitND2, InitND3, InitND4, InitND5, InitND6]; - SoilData.InitialValues.initT = [InitT0, InitT1, InitT2, InitT3, InitT4, InitT5, InitT6]; + InitialValues.initX = [InitX0, InitX1, InitX2, InitX3, InitX4, InitX5, InitX6]; + InitialValues.initND = [InitND1, InitND2, InitND3, InitND4, InitND5, InitND6]; + InitialValues.initT = [InitT0, InitT1, InitT2, InitT3, InitT4, InitT5, InitT6]; end diff --git a/src/CondL_h.m b/src/CondL_h.m index ae1f0ef9..0ade467c 100644 --- a/src/CondL_h.m +++ b/src/CondL_h.m @@ -1,10 +1,13 @@ -function [Theta_LL, Se, KfL_h, KfL_T, DTheta_LLh, hh, hh_frez, Theta_UU, DTheta_UUh, Theta_II, KL_h] = CondL_h(SoilConstants, SoilVariables, Theta_r, Theta_s, Alpha, hh, hh_frez, h_frez, n, m, Ks, NL, Theta_L, h, KIT, TT, Thmrlefc, POR, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, Theta_II, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se) +function [Theta_LL, Se, KfL_h, KfL_T, DTheta_LLh, hh, hh_frez, Theta_UU, DTheta_UUh, Theta_II, KL_h] = CondL_h(SoilVariables, Theta_r, Theta_s, Alpha, hh, hh_frez, h_frez, n, m, Ks, NL, Theta_L, h, KIT, TT, Thmrlefc, POR, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, Theta_II, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se) % load Constants Constants = io.define_constants(); + % get soil constants for StartInit + SoilConstants = io.getSoilConstants(); hd = SoilConstants.hd; hm = SoilConstants.hm; + Gama_hh = SoilVariables.Gama_hh; SFCC = 1; @@ -252,7 +255,6 @@ else Ratio_ice(ML, ND) = 0; end - if KIT MU_WN = Constants.MU_W0 * exp(Constants.MU1 / (8.31441 * (20 + 133.3))); if TT(MN) < -20 diff --git a/src/SOIL1.m b/src/SOIL1.m deleted file mode 100644 index 358c8dec..00000000 --- a/src/SOIL1.m +++ /dev/null @@ -1,50 +0,0 @@ -function SOIL1 - % This subroutine is caled after one time step to update the wetting history.If the change in average moisture content - % of the element during the last time step was in the opposite direction of that during the previous time step, the history - % is updated. As an approximation, only primary scanning curves are used, subject to the constraint that matric - % head and moisture content be continuous in time. - global EX ML NL ND Theta_L XOLD IS J Theta_LL XWRE IH Theta_s XK - - for ML = 1:NL - % Soil type index; - J = ML; - % J=IS(ML); - % The average moisture content of an element; - EX = 0.5 * (Theta_L(ML, 1) + Theta_L(ML, 2)); - % Has average trend of wetting in the element changed? If the trend is - % still in drying, keep running like this. Otherwise, change trend from - % drying to wetting. Then, IH value needs to be changed as 2, and XWRE - % needs to be re-evaluated. However, if the trend is still in wetting, - % keep running like that. Otherwise, change trend from wetting to drying. - if IH(ML) == 1 && XOLD(ML) <= (EX + 1e-12) % IH=1 means wetting. - XOLD(ML) = EX; - return - elseif IH(ML) == 2 % IH=2 means drying. - if XOLD(ML) >= (EX - 1e-12) - XOLD(ML) = EX; - return - else - IH(ML) = 1; - for ND = 1:2 - if (Theta_s(J) - Theta_LL(ML, ND)) < 1e-3 - XWRE(ML, ND) = Theta_s(J); - else - XWRE(ML, ND) = Theta_s(J) * (Theta_L(ML, ND) - Theta_LL(ML, ND)) / (Theta_s(J) - Theta_LL(ML, ND)); - end - end - XOLD(ML) = EX; - return - end - else - IH(ML) = 2; - for ND = 1:2 - if (Theta_LL(J) - XK(J)) < 1e-3 - XWRE(ML, ND) = XK(J); - else - XWRE(ML, ND) = Theta_LL(ML, ND) + Theta_s(J) * (Theta_L(ML, ND) / Theta_LL(ML, ND) - 1); - end - end - XOLD(ML) = EX; - return - end - end diff --git a/src/SOIL2.m b/src/SOIL2.m deleted file mode 100644 index c9962109..00000000 --- a/src/SOIL2.m +++ /dev/null @@ -1,80 +0,0 @@ -function [hh, COR, CORh, Theta_V, Theta_g, Se, KL_h, Theta_LL, DTheta_LLh, KfL_h, KfL_T, hh_frez, Theta_UU, DTheta_UUh, Theta_II] = SOIL2(SoilConstants, SoilVariables, hh, COR, hThmrl, NN, NL, TT, Tr, Hystrs, XWRE, Theta_s, IH, KIT, Theta_r, Alpha, n, m, Ks, Theta_L, h, Thmrlefc, POR, Theta_II, CORh, hh_frez, h_frez, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se) - - if hThmrl == 1 - - for MN = 1:NN - CORh(MN) = 0.0068; - COR(MN) = exp(-1 * CORh(MN) * (TT(MN) - Tr)); % *COR21(MN) - - if COR(MN) == 0 - COR(MN) = 1; - end - end - else - for MN = 1:NN - COR(MN) = 1; - end - end - - for MN = 1:NN - hhU(MN) = COR(MN) * hh(MN); - hh(MN) = hhU(MN); - end - [Theta_LL, Se, KfL_h, KfL_T, DTheta_LLh, hh, hh_frez, Theta_UU, DTheta_UUh, Theta_II, KL_h] = CondL_h(SoilConstants, SoilVariables, Theta_r, Theta_s, Alpha, hh, hh_frez, h_frez, n, m, Ks, NL, Theta_L, h, KIT, TT, Thmrlefc, POR, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, Theta_II, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se); - for MN = 1:NN - hhU(MN) = hh(MN); - hh(MN) = hhU(MN) / COR(MN); - end - - if Hystrs == 0 - for ML = 1:NL - J = ML; - % J=IS(ML); - for ND = 1:2 - Theta_V(ML, ND) = POR(J) - Theta_UU(ML, ND); % -Theta_II(ML,ND); % Theta_LL==>Theta_UU - if Theta_V(ML, ND) <= 1e-14 - Theta_V(ML, ND) = 1e-14; - end - Theta_g(ML, ND) = Theta_V(ML, ND); - end - end - else - for ML = 1:NL - J = ML; - % J=IS(ML); - for ND = 1:2 - if IH(ML) == 2 - if XWRE(ML, ND) < Theta_LL(ML, ND) - Theta_V(ML, ND) = POR(J) - Theta_LL(ML, ND) - Theta_II(ML, ND); - else - XSAVE = Theta_LL(ML, ND); - Theta_LL(ML, ND) = XSAVE * (1 + (XWRE(ML, ND) - Theta_LL(ML, ND)) / Theta_s(J)); - if KIT > 0 - DTheta_LLh(ML, ND) = DTheta_LLh(ML, ND) * (Theta_LL(ML, ND) / XSAVE - XSAVE / Theta_s(J)); - end - Theta_V(ML, ND) = POR(J) - Theta_LL(ML, ND) - Theta_II(ML, ND); - end - end - if IH(ML) == 1 - if XWRE(ML, ND) > Theta_LL(ML, ND) - XSAVE = Theta_LL(ML, ND); - Theta_LL(ML, ND) = (2 - XSAVE / Theta_s(J)) * XSAVE; - if KIT > 0 - DTheta_LLh(ML, ND) = 2 * DTheta_LLh(ML, ND) * (1 - XSAVE / Theta_s(J)); - end - Theta_V(ML, ND) = POR(J) - Theta_LL(ML, ND) - Theta_II(ML, ND); - else - Theta_LL(ML, ND) = Theta_LL(ML, ND) + XWRE(ML, ND) * (1 - Theta_LL(ML, ND) / Theta_s(J)); - if KIT > 0 - DTheta_LLh(ML, ND) = DTheta_LLh(ML, ND) * (1 - XWRE(ML, ND) / Theta_s(J)); - end - Theta_V(ML, ND) = POR(J) - Theta_LL(ML, ND) - Theta_II(ML, ND); - end - end - if Theta_V(ML, ND) <= 1e-14 % consider the negative conditions? - Theta_V(ML, ND) = 1e-14; - end - Theta_g(ML, ND) = Theta_V(ML, ND); - end - end - end diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 70571698..398bcf1e 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -32,40 +32,30 @@ CFG = '../config_file_crib.txt'; end disp (['Reading config from ', CFG]); -global InputPath OutputPath InitialConditionPath [InputPath, OutputPath, InitialConditionPath] = io.read_config(CFG); % Prepare forcing and soil data -global latitude longitude reference_height canopy_height sitename DELT Dur_tot SaturatedMC ResidualMC fieldMC fmax theta_s0 Ks0 +global DELT SaturatedMC ResidualMC fieldMC theta_s0 Ks0 [SiteProperties, SoilProperties, TimeProperties] = io.prepareInputData(InputPath); -landcoverClass = SiteProperties.landcoverClass; -latitude = SiteProperties.latitude; -longitude = SiteProperties.longitude; -reference_height = SiteProperties.reference_height; -canopy_height = SiteProperties.canopy_height; -sitename = SiteProperties.sitename; -SaturatedMC = SoilProperties.SaturatedMC; -ResidualMC = SoilProperties.ResidualMC; -fieldMC = SoilProperties.fieldMC; -fmax = SoilProperties.fmax; -theta_s0 = SoilProperties.theta_s0; -Ks0 = SoilProperties.Ks0; -DELT = TimeProperties.DELT; -Dur_tot = TimeProperties.Dur_tot; +landcoverClass = SiteProperties.landcoverClass; +SaturatedMC = SoilProperties.SaturatedMC; % used in calc_rssrbs +ResidualMC = SoilProperties.ResidualMC; % used in calc_rssrbs +fieldMC = SoilProperties.fieldMC; % used in calc_rssrbs + +theta_s0 = SoilProperties.theta_s0; % used in h_BC +Ks0 = SoilProperties.Ks0; % used in h_BC +DELT = TimeProperties.DELT; % used in h_BC % Load model settings: replacing "run Constants" ModelSettings = io.getModelSettings(); -global J rwuef SWCC Hystrs Thmrlefc Soilairefc hThmrl KIT DURTN KT TIME Delt_t NN ML nD -global W_Chg ThmrlCondCap ThermCond SSUR fc Tr T0 rroot SAVE NL DeltZ Tot_Depth Eqlspace -Tot_Depth = ModelSettings.Tot_Depth; -Eqlspace = ModelSettings.Eqlspace; +global J rwuef SWCC Thmrlefc Soilairefc hThmrl DURTN KT TIME Delt_t NN ML nD +global W_Chg ThmrlCondCap ThermCond SSUR fc T0 rroot SAVE NL DeltZ NL = ModelSettings.NL; DeltZ = ModelSettings.DeltZ; DeltZ_R = ModelSettings.DeltZ_R; J = ModelSettings.J; SWCC = ModelSettings.SWCC; -Hystrs = ModelSettings.Hystrs; Thmrlefc = ModelSettings.Thmrlefc; Soilairefc = ModelSettings.Soilairefc; hThmrl = ModelSettings.hThmrl; @@ -74,24 +64,22 @@ ThermCond = ModelSettings.ThermCond; SSUR = ModelSettings.SSUR; fc = ModelSettings.fc; -Tr = ModelSettings.Tr; T0 = ModelSettings.T0; rwuef = ModelSettings.rwuef; rroot = ModelSettings.rroot; SAVE = ModelSettings.SAVE; -KIT = ModelSettings.KIT; NIT = ModelSettings.NIT; KT = ModelSettings.KT; NN = ModelSettings.NN; ML = ModelSettings.ML; nD = ModelSettings.nD; % defined as global, and used in other scripts -DURTN = TimeProperties.DELT * TimeProperties.Dur_tot; % Duration of simulation period; +DURTN = TimeProperties.DELT * TimeProperties.Dur_tot; % used in Forcing_PARM, Duration of simulation period; TIME = 0 * TimeProperties.DELT; % Time of simulation released; Delt_t = TimeProperties.DELT; % Duration of time step [Unit of second] % load forcing data -ForcingData = io.loadForcingData(InputPath, TimeProperties, fmax, Tot_Depth); +ForcingData = io.loadForcingData(InputPath, TimeProperties, SoilProperties.fmax, ModelSettings.Tot_Depth); % global vars used in Forcing_PARM global Ta_msr RH_msr WS_msr Pg_msr Rn_msr Tmin Rns_msr VPD_msr LAI_msr G_msr Precip_msr @@ -107,45 +95,27 @@ Precip_msr = ForcingData.Precip_msr; Tmin = ForcingData.Tmin; -% load initial soil moisture and soil temperature -SoilData = io.loadSoilData(InputPath, TimeProperties, Tot_Depth, SoilProperties, ForcingData, SWCC); -global Tss -Tss = SoilData.Tss; % global vars used in Forcing_PARM - -global i tS MN ND hOLD TOLD h hh T TT P_g P_gg AVAIL Evap EXCESS QMT hN Trap -global SUMTIME TTT Theta_LLL CHK Theta_LL Theta_L Theta_UUU Theta_UU Theta_U Theta_III Theta_II Theta_I +global MN ND hOLD TOLD h hh T TT P_g P_gg Evap QMT hN Trap +global SUMTIME TTT Theta_LLL CHK Theta_LL Theta_L Theta_UUU Theta_UU Theta_U Theta_III Theta_II global AVAIL0 TIMEOLD SRT ALPHA alpha_h bx Srt CTT_PH CTT_LT CTT_g CTT_Lg c_unsat -global QL QL_h QL_T QV Qa KL_h Chh ChT Khh KhT Resis_a KfL_h KfL_T TT_CRIT T_CRIT TOLD_CRIT +global QL QL_h QL_T QV Qa KL_h Chh ChT Khh KhT Resis_a KfL_h KfL_T TT_CRIT global h_frez L_f CTT EPCT DTheta_LLh DTheta_LLT DTheta_UUh CKT Lambda_eff EfTCON TETCON DDhDZ DhDZ DTDZ DRHOVZ global DEhBAR DRHOVhDz EtaBAR D_Vg DRHOVTDz KLhBAR KLTBAR DTDBAR SAVEDTheta_LLh SAVEDTheta_UUh -global QVT QVH Sa HR QVa QLH QLT DVH DVT Se QL_a DPgDZ k_g V_A Theta_a Theta_V W WW D_Ta Ratio_ice -global KCHK FCHK hCHK SAVEhh_frez SAVETT INDICATOR thermal Xaa XaT Xah KL_T DRHOVT DRHOVh DRHODAt DRHODAz +global QVT QVH Sa HR QVa QLH QLT DVH DVT Se QL_a DPgDZ k_g V_A Theta_V W WW D_Ta Ratio_ice +global thermal Xaa XaT Xah KL_T DRHOVT DRHOVh DRHODAt DRHODAz global Theta_g Alpha_Lg Beta_g D_V D_A Eta ZETA MU_W Ks RHODA RHOV ETCON EHCAP -global uERR L QMTT QMBB Evapo trap RnSOIL PrecipO Beta_gBAR Alpha_LgBAR -global RWU EVAP theta_s0 Ks0 Precip Precipp Tss frac sfactortot sfactor fluxes lEstot lEctot NoTime +global L Evapo Beta_gBAR Alpha_LgBAR +global RWU EVAP theta_s0 Ks0 Precip Tss frac sfactortot sfactor fluxes lEstot lEctot NoTime % Get initial values InitialValues = init.defineInitialValues(TimeProperties.Dur_tot); alpha_h = InitialValues.alpha_h; bx = InitialValues.bx; Srt = InitialValues.Srt; -DTheta_LLh = InitialValues.DTheta_LLh; -DTheta_UUh = InitialValues.DTheta_UUh; SAVEDTheta_UUh = InitialValues.SAVEDTheta_UUh; SAVEDTheta_LLh = InitialValues.SAVEDTheta_LLh; Ratio_ice = InitialValues.Ratio_ice; -KL_h = InitialValues.KL_h; -KfL_h = InitialValues.KfL_h; -KfL_T = InitialValues.KfL_T; KL_T = InitialValues.KL_T; -Theta_L = InitialValues.Theta_L; -Theta_LL = InitialValues.Theta_LL; -Theta_U = InitialValues.Theta_U; -Theta_UU = InitialValues.Theta_UU; -Theta_I = InitialValues.Theta_I; -Theta_II = InitialValues.Theta_II; -Se = InitialValues.Se; -Theta_V = InitialValues.Theta_V; Lambda_eff = InitialValues.Lambda_eff; W = InitialValues.W; WW = InitialValues.WW; @@ -154,7 +124,6 @@ D_V = InitialValues.D_V; Eta = InitialValues.Eta; D_A = InitialValues.D_A; -Theta_g = InitialValues.Theta_g; EHCAP = InitialValues.EHCAP; Chh = InitialValues.Chh; ChT = InitialValues.ChT; @@ -204,10 +173,6 @@ EVAP = InitialValues.EVAP; P_g = InitialValues.P_g; P_gg = InitialValues.P_gg; -h = InitialValues.h; -h_frez = InitialValues.h_frez; -T = InitialValues.T; -TT = InitialValues.TT; T_CRIT = InitialValues.T_CRIT; TT_CRIT = InitialValues.TT_CRIT; EPCT = InitialValues.EPCT; @@ -227,7 +192,7 @@ global f0 L_WT Kha Vvh VvT Chg C1 C2 C3 C4 C5 C6 Cah CaT Caa Kah KaT Kaa Vah VaT Vaa Cag CTh CTa KTh KTT KTa global VTT VTh VTa CTg Kcva Kcah KcaT Kcaa Ccah CcaT Ccaa Ksoil SMC bbx wfrac Ta Ts U HR_a Rns Rnl Rn -global RHOV_s DRHOV_sT LL P_gOLD hh_frez XCAP Tbtm r_a_SOIL Rn_SOIL PSItot Tsur SH MO Zeta_MO TopPg Tp_t DhT RHS C7 C9 +global RHOV_s DRHOV_sT Tbtm r_a_SOIL Rn_SOIL SH MO Zeta_MO TopPg Tp_t RHS C7 C9 f0 = InitialValues.f0; L_WT = InitialValues.L_WT; Kha = InitialValues.Kha; @@ -282,25 +247,19 @@ Zeta_MO = InitialValues.Zeta_MO; TopPg = InitialValues.TopPg; Tp_t = InitialValues.Tp_t; -DhT = InitialValues.DhT; RHS = InitialValues.RHS; C7 = InitialValues.C7; C9 = InitialValues.C9; RHOV_s = InitialValues.RHOV_s; DRHOV_sT = InitialValues.DRHOV_sT; -LL = InitialValues.LL; P_gOLD = InitialValues.P_gOLD; -hh_frez = InitialValues.hh_frez; -XCAP = InitialValues.XCAP; Tbtm = InitialValues.Tbtm; r_a_SOIL = InitialValues.r_a_SOIL; Rn_SOIL = InitialValues.Rn_SOIL; -PSItot = InitialValues.PSItot; -Tsur = InitialValues.Tsur; %% 1. define Constants Constants = io.define_constants(); -global g RHOL RHOI Rv RDA MU_a Lambda1 Lambda2 Lambda3 c_a c_V c_L Hc GWT c_i Gamma0 Gamma_w RHO_bulk Rl +global g RHOL RHOI Rv RDA MU_a Lambda1 Lambda2 Lambda3 c_a c_V c_L Hc c_i Gamma0 Gamma_w RHO_bulk Rl g = Constants.g; RHOL = Constants.RHOL; RHOI = Constants.RHOI; @@ -316,7 +275,6 @@ Hc = Constants.Hc; % used in other scripts not here! -GWT = Constants.GWT; % used in CondL_T Gamma0 = Constants.Gamma0; % used in other scripts! Gamma_w = Constants.Gamma_w; % used in other scripts! c_i = Constants.c_i; % used in EnrgyPARM! @@ -327,19 +285,11 @@ [Rl, Ztot] = Initial_root_biomass(RTB, ModelSettings.DeltZ_R, rroot, ML, SiteProperties.landcoverClass(1)); %% 2. simulation options -path_input = InputPath; % path of all inputs +path_input = InputPath; % path of all inputs path_of_code = cd; -useXLSX = 1; % set it to 1 or 0, the current stemmus-scope does not support useXLSX=0 -if useXLSX == 0 - % parameter_file = { 'setoptions.m', 'filenames.m', 'inputdata.txt'}; - % Read parameter file which is 'input_data.xlsx' and return it as options. - % options = io.setOptions(parameter_file,path_input); - warning("the current stemmus-scope does not support useXLSX=0"); -else - parameter_file = {'input_data.xlsx'}; - options = io.readStructFromExcel([path_input char(parameter_file)], 'options', 2, 1); -end +parameter_file = {'input_data.xlsx'}; +options = io.readStructFromExcel([path_input char(parameter_file)], 'options', 2, 1); if options.simulation > 2 || options.simulation < 0 fprintf('\n simulation option should be between 0 and 2 \r'); @@ -347,14 +297,9 @@ end %% 3. file names -% the current stemmus-scope does not support useXLSX=0 -if useXLSX == 0 - run([path_input parameter_file{2}]); -else - [dummy, X] = xlsread([path_input char(parameter_file)], 'filenames'); - j = find(~strcmp(X(:, 2), {''})); - X = X(j, 1:end); -end +[dummy, X] = xlsread([path_input char(parameter_file)], 'filenames'); +j = find(~strcmp(X(:, 2), {''})); +X = X(j, 1:end); F = struct('FileID', {'Simulation_Name', 'soil_file', 'leaf_file', 'atmos_file'... 'Dataset_dir', 't_file', 'year_file', 'Rin_file', 'Rli_file' ... @@ -368,16 +313,11 @@ end %% 4. input data -% the current stemmus-scope does not support useXLSX=0 -if useXLSX == 0 - X = textread([path_input parameter_file{3}], '%s'); %#ok - N = str2double(X); -else - [N, X] = xlsread([path_input char(parameter_file)], 'inputdata', ''); - X = X(9:end, 1); -end +[N, X] = xlsread([path_input char(parameter_file)], 'inputdata', ''); +X = X(9:end, 1); % Create a structure holding Scope parameters +useXLSX = 1; % set it to 1 or 0, the current stemmus-scope does not support 0 [ScopeParameters, options] = parameters.loadParameters(options, useXLSX, X, F, N); % Define the location information @@ -496,23 +436,36 @@ atmo.M = helpers.aggreg(atmfile, spectral.SCOPEspec); %% 13. create output files and +[Output_dir, fnames] = io.create_output_files_binary(parameter_file, SiteProperties.sitename, path_of_code, path_input, path_output, spectral, options); + %% Initialize Temperature, Matric potential and soil air pressure. -[Output_dir, fnames] = io.create_output_files_binary(parameter_file, sitename, path_of_code, path_input, path_output, spectral, options); +% Define soil variables for StartInit +VanGenuchten = init.setVanGenuchtenParameters(SoilProperties); +SoilVariables = init.defineSoilVariables(InitialValues, SoilProperties, VanGenuchten); -% SoilConstants for init -SoilConstants = init.setSoilConstants(InitialValues, SoilProperties, ForcingData); -[SoilConstants, SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilConstants, SoilProperties, SoilData, SiteProperties); +% Add initial soil moisture and soil temperature +global Tss % global vars used in Forcing_PARM +[SoilInitialValues, BtmX, BtmT, Tss] = io.loadSoilInitialValues(InputPath, TimeProperties, SoilProperties, ForcingData); +SoilVariables.InitialValues = SoilInitialValues; +SoilVariables.BtmX = BtmX; +SoilVariables.BtmT = BtmT; +SoilVariables.Tss = Tss; + +% Run StartInit +[SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten); %% get variables that are defined global and are used by other scripts global hm hd hh_frez XWRE POR IH IS XK XWILT KLT_Switch DVT_Switch KaT_Switch global ISFT Imped XSOC Lamda Phi_s XCAP Gama_hh Gama_h SAVEhh COR CORh global Theta_s Theta_r Theta_f m n Alpha global HCAP SF TCA GA1 GA2 GB1 GB2 HCD ZETA0 CON0 PS1 PS2 FEHCAP -global TCON_dry TPS1 TPS2 TCON0 TCON_s XOLD +global TCON_dry TPS1 TPS2 TCON0 TCON_s +% get soil constants for StartInit +SoilConstants = io.getSoilConstants(); hm = SoilConstants.hm; hd = SoilConstants.hd; -hh_frez = SoilVariables.hh_frez; + XWRE = SoilVariables.XWRE; POR = SoilVariables.POR; IH = SoilVariables.IH; @@ -531,8 +484,6 @@ Gama_hh = SoilVariables.Gama_hh; Gama_h = SoilVariables.Gama_h; SAVEhh = SoilVariables.SAVEhh; -COR = SoilVariables.COR; -CORh = SoilVariables.CORh; Theta_s = VanGenuchten.Theta_s; Theta_r = VanGenuchten.Theta_r; Theta_f = VanGenuchten.Theta_f; @@ -558,33 +509,20 @@ TPS2 = ThermalConductivity.TPS2; FEHCAP = ThermalConductivity.FEHCAP; TCON0 = ThermalConductivity.TCON0; -XOLD = SoilVariables.XOLD; % used in SOIL1 %% these vars are defined as global at the begining of this script %% because they are both input and output of StartInit -KfL_T = SoilConstants.KfL_T; -Theta_II = SoilConstants.Theta_II; -Theta_I = SoilConstants.Theta_I; -Theta_UU = SoilConstants.Theta_UU; -Theta_U = SoilConstants.Theta_U; +Theta_I = SoilVariables.Theta_I; +Theta_U = SoilVariables.Theta_U; +Theta_L = SoilVariables.Theta_L; T = SoilVariables.T; h = SoilVariables.h; TT = SoilVariables.TT; -hh = SoilVariables.hh; Ks = SoilVariables.Ks; h_frez = SoilVariables.h_frez; -Theta_L = SoilVariables.Theta_L; -Theta_LL = SoilVariables.Theta_LL; -Theta_V = SoilVariables.Theta_V; -Theta_g = SoilVariables.Theta_g; -Se = SoilVariables.Se; -KL_h = SoilVariables.KL_h; -DTheta_LLh = SoilVariables.DTheta_LLh; -KfL_h = SoilVariables.KfL_h; -DTheta_UUh = SoilVariables.DTheta_UUh; %% The boundary condition information settings -BoundaryCondition = init.setBoundaryCondition(SoilVariables, SoilConstants, landcoverClass(1)); +BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, landcoverClass(1)); %% get global vars global NBCh NBCT NBChB NBCTB BCh DSTOR DSTOR0 RS NBChh DSTMAX IRPT1 IRPT2 @@ -610,6 +548,16 @@ BCP = BoundaryCondition.BCP; BtmPg = BoundaryCondition.BtmPg; +% Outputs of UpdateSoilWaterContent used in step Run the model +hh = SoilVariables.hh; +hh_frez = SoilVariables.hh_frez; + +KL_h = SoilVariables.KL_h; +KfL_h = SoilVariables.KfL_h; + +% Outputs of UpdateSoilWaterContent used in io.select_input in the loop +Theta_LL = SoilVariables.Theta_LL; + %% 14. Run the model disp('The calculations start now'); calculate = 1; @@ -620,20 +568,20 @@ hCHK = zeros(1, NN); TIMELAST = 0; -% Is the tS(time step) needed to be added with 1? % Cause the start of simulation period is from 0mins, while the input data start from 30mins. tS = DURTN / Delt_t; SAVEtS = tS; kk = 0; % DELT=Delt_t; TimeStep = []; -TEND = TIME + DURTN; % Time to be reached at the end of simulation period; -Delt_t0 = Delt_t; % Duration of last time step; -for i = 1:1:Dur_tot - KT = KT + 1; % Counting Number of timesteps +TEND = TIME + DURTN; % Time to be reached at the end of simulation period +Delt_t0 = Delt_t; % Duration of last time step +TOLD_CRIT = []; +for i = 1:1:TimeProperties.Dur_tot + KT = KT + 1; % Counting Number of timesteps if KT > 1 && Delt_t > (TEND - TIME) - Delt_t = TEND - TIME; % If Delt_t is changed due to excessive change of state variables, the judgement of the last time step is excuted. + Delt_t = TEND - TIME; % If Delt_t is changed due to excessive change of state variables, the judgement of the last time step is excuted. end - TIME = TIME + Delt_t; % The time elapsed since start of simulation + TIME = TIME + Delt_t; % The time elapsed since start of simulation TimeStep(KT, 1) = Delt_t; SUMTIME(KT, 1) = TIME; Processing = TIME / TEND; @@ -644,9 +592,7 @@ k = NoTime(KT); end %%%%% Updating the state variables. %%%%%%%%%%%%%%%%%%%%%%%%%%%% - - L_f = 0; % latent heat of freezing fusion J Kg-1 - T0 = 273.15; + L_f = 0; % ignore Freeze/Thaw, see issue 139 TT_CRIT(NN) = T0; % unit K hOLD_frez = []; if IRPT1 == 0 && IRPT2 == 0 && ISFT == 0 @@ -687,9 +633,7 @@ end if calculate - iter.counter = 0; - LIDF_file = char(F(22).FileName); if ~isempty(LIDF_file) canopy.lidf = dlmread([path_input, 'leafangles/', LIDF_file], '', 3, 0); @@ -708,10 +652,8 @@ leafopt = fversion(spectral, leafbio, optipar); leafbio.V2Z = 1; leafoptZ = fversion(spectral, leafbio, optipar); - IwlP = spectral.IwlP; IwlT = spectral.IwlT; - rho(IwlP) = leafopt.refl; tau(IwlP) = leafopt.tran; rlast = rho(nwlP); @@ -738,16 +680,13 @@ end leafopt.refl = rho; % extended wavelength ranges are stored in structures leafopt.tran = tau; - reflZ = leafopt.refl; tranZ = leafopt.tran; reflZ(1:300) = leafoptZ.refl(1:300); tranZ(1:300) = leafoptZ.tran(1:300); leafopt.reflZ = reflZ; leafopt.tranZ = tranZ; - soil.refl = rs; - soil.Ts = meteo.Ta * ones(2, 1); % initial soil surface temperature if length(F(4).FileName) > 1 && options.simulation == 0 @@ -838,7 +777,7 @@ DSTOR0 = DSTOR; if KT > 1 - run SOIL1; + SoilVariables.XWRE = updateWettingHistory(SoilVariables, VanGenuchten); end for ML = 1:NL @@ -877,11 +816,38 @@ end end run Forcing_PARM; - %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for KIT = 1:NIT % Start the iteration procedure in a time step. [TT_CRIT, hh_frez] = HT_frez(hh, T0, g, L_f, TT, NN, hd, Tmin); - [hh, COR, CORh, Theta_V, Theta_g, Se, KL_h, Theta_LL, DTheta_LLh, KfL_h, KfL_T, hh_frez, Theta_UU, DTheta_UUh, Theta_II] = SOIL2(SoilConstants, SoilVariables, hh, COR, hThmrl, NN, NL, TT, Tr, Hystrs, XWRE, Theta_s, IH, KIT, Theta_r, Alpha, n, m, Ks, Theta_L, h, Thmrlefc, POR, Theta_II, CORh, hh_frez, h_frez, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se); + + % update inputs for UpdateSoilWaterContent + SoilVariables.TT_CRIT = TT_CRIT; + SoilVariables.hh_frez = hh_frez; + SoilVariables.h = h; + SoilVariables.hh = hh; + SoilVariables.KL_h = KL_h; + SoilVariables.KfL_h = KfL_h; + SoilVariables.TT = TT; + SoilVariables.h_frez = h_frez; + SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten); + % these can be removed after refactoring functions below + h = SoilVariables.h; + hh = SoilVariables.hh; + COR = SoilVariables.COR; + CORh = SoilVariables.CORh; + Theta_V = SoilVariables.Theta_V; + Theta_g = SoilVariables.Theta_g; + Theta_LL = SoilVariables.Theta_LL; + Se = SoilVariables.Se; + KL_h = SoilVariables.KL_h; + DTheta_LLh = SoilVariables.DTheta_LLh; + KfL_h = SoilVariables.KfL_h; + KfL_T = SoilVariables.KfL_T; + hh_frez = SoilVariables.hh_frez; + Theta_UU = SoilVariables.Theta_UU; + DTheta_UUh = SoilVariables.DTheta_UUh; + Theta_II = SoilVariables.Theta_II; + [KL_T] = CondL_T(NL); [RHOV, DRHOVh, DRHOVT] = Density_V(TT, hh, g, Rv, NN); [W, WW, MU_W, D_Ta] = CondL_Tdisp(InitialValues, POR, Theta_LL, Theta_L, SSUR, RHOL, TT, Theta_s, h, hh, W_Chg, NL, nD, Delt_t, Theta_g, KLT_Switch); @@ -891,6 +857,7 @@ [k_g] = Condg_k_g(POR, NL, m, Theta_g, g, MU_W, Ks, RHOL, SWCC, Imped, Ratio_ice, Soilairefc, MN); [D_V, Eta, D_A] = CondV_DE(Theta_LL, TT, fc, Theta_s, NL, nD, Theta_g, POR, ThmrlCondCap, ZETA, XK, DVT_Switch, Theta_UU); [D_Vg, V_A, Beta_g, DPgDZ, Beta_gBAR, Alpha_LgBAR] = CondV_DVg(P_gg, Theta_g, Sa, V_A, k_g, MU_a, DeltZ, Alpha_Lg, KaT_Switch, Theta_s, Se, NL, DPgDZ, Beta_gBAR, Alpha_LgBAR, Beta_g); + run h_sub; if NBCh == 1 DSTOR = 0; @@ -931,7 +898,33 @@ KIT; KIT = 0; [TT_CRIT, hh_frez] = HT_frez(hh, T0, g, L_f, TT, NN, hd, Tmin); - [hh, COR, CORh, Theta_V, Theta_g, Se, KL_h, Theta_LL, DTheta_LLh, KfL_h, KfL_T, hh_frez, Theta_UU, DTheta_UUh, Theta_II] = SOIL2(SoilConstants, SoilVariables, hh, COR, hThmrl, NN, NL, TT, Tr, Hystrs, XWRE, Theta_s, IH, KIT, Theta_r, Alpha, n, m, Ks, Theta_L, h, Thmrlefc, POR, Theta_II, CORh, hh_frez, h_frez, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se); + + % updates inputs for UpdateSoilWaterContent + SoilVariables.TT_CRIT = TT_CRIT; + SoilVariables.hh_frez = hh_frez; + SoilVariables.h = h; + SoilVariables.hh = hh; + SoilVariables.TT = TT; + SoilVariables.h_frez = h_frez; + SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten); + + % these can be removed after refactoring codes below + h = SoilVariables.h; + hh = SoilVariables.hh; + COR = SoilVariables.COR; + CORh = SoilVariables.CORh; + Theta_V = SoilVariables.Theta_V; + Theta_g = SoilVariables.Theta_g; + Se = SoilVariables.Se; + KL_h = SoilVariables.KL_h; + Theta_LL = SoilVariables.Theta_LL; + DTheta_LLh = SoilVariables.DTheta_LLh; + KfL_h = SoilVariables.KfL_h; + KfL_T = SoilVariables.KfL_T; + hh_frez = SoilVariables.hh_frez; + Theta_UU = SoilVariables.Theta_UU; + DTheta_UUh = SoilVariables.DTheta_UUh; + Theta_II = SoilVariables.Theta_II; SAVEtS = tS; if IRPT1 == 0 && IRPT2 == 0 @@ -985,10 +978,7 @@ %% soil layer information SoilLayer.thickness = ModelSettings.DeltZ_R; -SoilLayer.depth = Ztot'; +SoilLayer.depth = Ztot'; % used in Initial_root_biomass io.bin_to_csv(fnames, n_col, k, options, SoilLayer); save([Output_dir, 'output.mat']); -% if options.makeplots -% plot.plots(Output_dir) -% end diff --git a/src/StartInit.m b/src/StartInit.m index 297993d7..6e32517f 100644 --- a/src/StartInit.m +++ b/src/StartInit.m @@ -1,94 +1,44 @@ -function [SoilConstants, SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilConstants, SoilProperties, SoilData, SiteProperties) +function [SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten) Ksh = repelem(18 / (3600 * 24), 6); BtmKsh = Ksh(6); Ksh0 = Ksh(1); + LatentHeatOfFreezing = 3.34 * 1e5; % latent heat of freezing fusion i Kg-1 + KIT = 0; % KIT is used to count the number of iteration in a time step; + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Considering soil hetero effect modify date: 20170103 %%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - VanGenuchten = init.setVanGenuchtenParameters(SoilProperties); - SoilVariables = init.setSoilVariables(SoilProperties, SoilConstants, VanGenuchten); - [SoilVariables, VanGenuchten] = init.applySoilHeteroEffect(SoilProperties, SoilConstants, SoilData, SoilVariables, VanGenuchten); + [SoilVariables, VanGenuchten] = init.applySoilHeteroEffect(SoilProperties, SoilVariables, VanGenuchten); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Considering soil hetero effect modify date: 20170103 %%%%%%%%%%%% %%%%%% Perform initial freezing temperature for each soil type.%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - SoilVariables = init.applySoilHeteroWithInitialFreezing(SoilConstants, SoilVariables); + SoilVariables = init.applySoilHeteroWithInitialFreezing(LatentHeatOfFreezing, SoilVariables); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Perform initial thermal calculations for each soil type.%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - ThermalConductivity = init.calculateInitialThermal(SoilConstants, SoilVariables, VanGenuchten); + ThermalConductivity = init.calculateInitialThermal(SoilVariables, VanGenuchten); - % According to hh value get the Theta_LL - % run SOIL2; % For calculating Theta_LL,used in first Balance calculation. + % these will be updated in UpdateSoilWaterContent + SoilVariables.COR = []; + SoilVariables.CORh = []; + SoilVariables = UpdateSoilWaterContent(KIT, LatentHeatOfFreezing, SoilVariables, VanGenuchten); - Theta_L = SoilConstants.Theta_L; - Theta_LL = SoilConstants.Theta_LL; - Theta_V = SoilConstants.Theta_V; - Theta_g = SoilConstants.Theta_g; - Se = SoilConstants.Se; - KL_h = SoilConstants.KL_h; - DTheta_LLh = SoilConstants.DTheta_LLh; - KfL_T = SoilConstants.KfL_T; - Theta_II = SoilConstants.Theta_II; - Theta_I = SoilConstants.Theta_I; - Theta_UU = SoilConstants.Theta_UU; - Theta_U = SoilConstants.Theta_U; - TT_CRIT = SoilConstants.TT_CRIT; - KfL_h = SoilConstants.KfL_h; - DTheta_UUh = SoilConstants.DTheta_UUh; + Theta_L = SoilVariables.Theta_L; + Theta_I = SoilVariables.Theta_I; + Theta_U = SoilVariables.Theta_U; - % get Constants - Constants = io.define_constants(); - g = Constants.g; - RHOI = Constants.RHOI; - RHOL = Constants.RHOL; + Theta_LL = SoilVariables.Theta_LL; + Theta_UU = SoilVariables.Theta_UU; + Theta_II = SoilVariables.Theta_II; - % after refatoring SOIL2, these lines can be removed later, see issue 95! % get model settings ModelSettings = io.getModelSettings(); - NN = ModelSettings.NN; - NL = ModelSettings.NL; - SWCC = ModelSettings.SWCC; - Thmrlefc = ModelSettings.Thmrlefc; - T0 = ModelSettings.T0; - Tr = ModelSettings.Tr; - KIT = ModelSettings.KIT; - hThmrl = ModelSettings.hThmrl; - Hystrs = ModelSettings.Hystrs; - - Theta_s = VanGenuchten.Theta_s; - Theta_r = VanGenuchten.Theta_r; - Alpha = VanGenuchten.Alpha; - n = VanGenuchten.n; - m = VanGenuchten.m; - - POR = SoilVariables.POR; - h = SoilVariables.h; - Lamda = SoilVariables.Lamda; - Phi_s = SoilVariables.Phi_s; - - h_frez = SoilVariables.h_frez; - hh_frez = SoilVariables.hh_frez; - TT = SoilVariables.TT; - hh = SoilVariables.hh; - XWRE = SoilVariables.XWRE; - IH = SoilVariables.IH; - Ks = SoilVariables.Ks; - Imped = SoilVariables.Imped; - XCAP = SoilVariables.XCAP; - - L_f = 3.34 * 1e5; % latent heat of freezing fusion J Kg-1 - T0 = 273.15; % unit K - - % these two are defined inside SOIL2, see issue 95 - COR = []; - CORh = []; - [hh, COR, CORh, Theta_V, Theta_g, Se, KL_h, Theta_LL, DTheta_LLh, KfL_h, KfL_T, hh_frez, Theta_UU, DTheta_UUh, Theta_II] = SOIL2(SoilConstants, SoilVariables, hh, COR, hThmrl, NN, NL, TT, Tr, Hystrs, XWRE, Theta_s, IH, KIT, Theta_r, Alpha, n, m, Ks, Theta_L, h, Thmrlefc, POR, Theta_II, CORh, hh_frez, h_frez, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se); for i = 1:ModelSettings.NL Theta_L(i, 1) = Theta_LL(i, 1); @@ -115,24 +65,7 @@ KLa_Switch = 1; end - % TODO after refatoring SOIL2, these two lines can be removed! see issue 96! - SoilConstants.KfL_T = KfL_T; - SoilConstants.Theta_II = Theta_II; - SoilConstants.Theta_I = Theta_I; - SoilConstants.Theta_UU = Theta_UU; - SoilConstants.Theta_U = Theta_U; - - SoilVariables.hh_frez = hh_frez; - SoilVariables.hh = hh; - SoilVariables.COR = COR; - SoilVariables.CORh = CORh; - SoilVariables.Se = Se; - SoilVariables.KL_h = KL_h; + SoilVariables.Theta_I = Theta_I; + SoilVariables.Theta_U = Theta_U; SoilVariables.Theta_L = Theta_L; - SoilVariables.Theta_LL = Theta_LL; - SoilVariables.DTheta_LLh = DTheta_LLh; - SoilVariables.KfL_h = KfL_h; - SoilVariables.KfL_T = KfL_T; - SoilVariables.Theta_V = Theta_V; - SoilVariables.Theta_g = Theta_g; - SoilVariables.DTheta_UUh = DTheta_UUh; +end diff --git a/src/UpdateSoilWaterContent.m b/src/UpdateSoilWaterContent.m new file mode 100644 index 00000000..61b19ba9 --- /dev/null +++ b/src/UpdateSoilWaterContent.m @@ -0,0 +1,152 @@ +function SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten) + %{ + Update SoilWaterContent i.e. Theta_LL. + %} + + % these lines can be removed after issue 181 + % get model settings + ModelSettings = io.getModelSettings(); + + NN = ModelSettings.NN; + NL = ModelSettings.NL; + SWCC = ModelSettings.SWCC; + Thmrlefc = ModelSettings.Thmrlefc; + T0 = ModelSettings.T0; + Tr = ModelSettings.Tr; + hThmrl = ModelSettings.hThmrl; + Hystrs = ModelSettings.Hystrs; + + Theta_s = VanGenuchten.Theta_s; + Theta_r = VanGenuchten.Theta_r; + Alpha = VanGenuchten.Alpha; + n = VanGenuchten.n; + m = VanGenuchten.m; + + POR = SoilVariables.POR; + h = SoilVariables.h; + Lamda = SoilVariables.Lamda; + Phi_s = SoilVariables.Phi_s; + + h_frez = SoilVariables.h_frez; + hh_frez = SoilVariables.hh_frez; + TT = SoilVariables.TT; + hh = SoilVariables.hh; + COR = SoilVariables.COR; + CORh = SoilVariables.CORh; + XWRE = SoilVariables.XWRE; + IH = SoilVariables.IH; + Ks = SoilVariables.Ks; + Imped = SoilVariables.Imped; + XCAP = SoilVariables.XCAP; + + Theta_L = SoilVariables.Theta_L; + Theta_LL = SoilVariables.Theta_LL; + Theta_V = SoilVariables.Theta_V; + Theta_g = SoilVariables.Theta_g; + Se = SoilVariables.Se; + KL_h = SoilVariables.KL_h; + DTheta_LLh = SoilVariables.DTheta_LLh; + KfL_T = SoilVariables.KfL_T; + Theta_II = SoilVariables.Theta_II; + Theta_I = SoilVariables.Theta_I; + Theta_UU = SoilVariables.Theta_UU; + Theta_U = SoilVariables.Theta_U; + TT_CRIT = SoilVariables.TT_CRIT; + KfL_h = SoilVariables.KfL_h; + DTheta_UUh = SoilVariables.DTheta_UUh; + + % get Constants + Constants = io.define_constants(); + g = Constants.g; + RHOI = Constants.RHOI; + RHOL = Constants.RHOL; + + if hThmrl == 1 + for MN = 1:NN + CORh(MN) = 0.0068; + COR(MN) = exp(-1 * CORh(MN) * (TT(MN) - Tr)); % *COR21(MN) + + if COR(MN) == 0 + COR(MN) = 1; + end + end + else + for MN = 1:NN + COR(MN) = 1; + end + end + + for MN = 1:NN + hhU(MN) = COR(MN) * hh(MN); + hh(MN) = hhU(MN); + end + [Theta_LL, Se, KfL_h, KfL_T, DTheta_LLh, hh, hh_frez, Theta_UU, DTheta_UUh, Theta_II, KL_h] = CondL_h(SoilVariables, Theta_r, Theta_s, Alpha, hh, hh_frez, h_frez, n, m, Ks, NL, Theta_L, h, KIT, TT, Thmrlefc, POR, SWCC, Theta_U, XCAP, Phi_s, RHOI, RHOL, Lamda, Imped, L_f, g, T0, TT_CRIT, Theta_II, KfL_h, KfL_T, KL_h, Theta_UU, Theta_LL, DTheta_LLh, DTheta_UUh, Se); + for MN = 1:NN + hhU(MN) = hh(MN); + hh(MN) = hhU(MN) / COR(MN); + end + if Hystrs == 0 + for i = 1:NL + for ND = 1:2 + Theta_V(i, ND) = POR(i) - Theta_UU(i, ND); % -Theta_II(i,ND); % Theta_LL==>Theta_UU + if Theta_V(i, ND) <= 1e-14 + Theta_V(i, ND) = 1e-14; + end + Theta_g(i, ND) = Theta_V(i, ND); + end + end + else + for i = 1:NL + for ND = 1:2 + if IH(i) == 2 + if XWRE(i, ND) < Theta_LL(i, ND) + Theta_V(i, ND) = POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); + else + XSAVE = Theta_LL(i, ND); + Theta_LL(i, ND) = XSAVE * (1 + (XWRE(i, ND) - Theta_LL(i, ND)) / Theta_s(i)); + if KIT > 0 + DTheta_LLh(i, ND) = DTheta_LLh(i, ND) * (Theta_LL(i, ND) / XSAVE - XSAVE / Theta_s(i)); + end + Theta_V(i, ND) = POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); + end + end + if IH(i) == 1 + if XWRE(i, ND) > Theta_LL(i, ND) + XSAVE = Theta_LL(i, ND); + Theta_LL(i, ND) = (2 - XSAVE / Theta_s(i)) * XSAVE; + if KIT > 0 + DTheta_LLh(i, ND) = 2 * DTheta_LLh(i, ND) * (1 - XSAVE / Theta_s(i)); + end + Theta_V(i, ND) = POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); + else + Theta_LL(i, ND) = Theta_LL(i, ND) + XWRE(i, ND) * (1 - Theta_LL(i, ND) / Theta_s(i)); + if KIT > 0 + DTheta_LLh(i, ND) = DTheta_LLh(i, ND) * (1 - XWRE(i, ND) / Theta_s(i)); + end + Theta_V(i, ND) = POR(i) - Theta_LL(i, ND) - Theta_II(i, ND); + end + end + if Theta_V(i, ND) <= 1e-14 % consider the negative conditions? + Theta_V(i, ND) = 1e-14; + end + Theta_g(i, ND) = Theta_V(i, ND); + end + end + end + SoilVariables.KfL_T = KfL_T; + SoilVariables.Theta_II = Theta_II; + SoilVariables.Theta_UU = Theta_UU; + + SoilVariables.hh_frez = hh_frez; + SoilVariables.hh = hh; + SoilVariables.COR = COR; + SoilVariables.CORh = CORh; + SoilVariables.Se = Se; + SoilVariables.KL_h = KL_h; + SoilVariables.Theta_LL = Theta_LL; + SoilVariables.DTheta_LLh = DTheta_LLh; + SoilVariables.KfL_h = KfL_h; + SoilVariables.Theta_V = Theta_V; + SoilVariables.Theta_g = Theta_g; + SoilVariables.DTheta_UUh = DTheta_UUh; +end diff --git a/src/updateWettingHistory.m b/src/updateWettingHistory.m new file mode 100644 index 00000000..f80db6c8 --- /dev/null +++ b/src/updateWettingHistory.m @@ -0,0 +1,61 @@ +function XWRE = updateWettingHistory(SoilVariables, VanGenuchten) + %{ + This subroutine is caled after one time step to update the wetting + history. If the change in average moisture content of the element during the + last time step was in the opposite direction of that during the previous + time step, the history is updated. As an approximation, only primary + scanning curves are used, subject to the constraint that matric head and + moisture content be continuous in time. + %} + % get model settings + ModelSettings = io.getModelSettings(); + + XOLD = SoilVariables.XOLD; + Theta_L = SoilVariables.Theta_L; + Theta_LL = SoilVariables.Theta_LL; + XWRE = SoilVariables.XWRE; + IH = SoilVariables.IH; + XK = SoilVariables.XK; + Theta_s = VanGenuchten.Theta_s; + + for i = 1:ModelSettings.NL + % The average moisture content of an element; + EX = 0.5 * (Theta_L(i, 1) + Theta_L(i, 2)); + % Has average trend of wetting in the element changed? If the trend is + % still in drying, keep running like this. Otherwise, change trend from + % drying to wetting. Then, IH value needs to be changed as 2, and XWRE + % needs to be re-evaluated. However, if the trend is still in wetting, + % keep running like that. Otherwise, change trend from wetting to drying. + if IH(i) == 1 && XOLD(i) <= (EX + 1e-12) % IH=1 means wetting. + XOLD(i) = EX; + return + elseif IH(i) == 2 % IH=2 means drying. + if XOLD(i) >= (EX - 1e-12) + XOLD(i) = EX; + return + else + IH(i) = 1; + for j = 1:2 + if (Theta_s(i) - Theta_LL(i, j)) < 1e-3 + XWRE(i, j) = Theta_s(i); + else + XWRE(i, j) = Theta_s(i) * (Theta_L(i, j) - Theta_LL(i, j)) / (Theta_s(i) - Theta_LL(i, j)); + end + end + XOLD(i) = EX; + return + end + else + IH(i) = 2; + for j = 1:2 + if (Theta_LL(i) - XK(i)) < 1e-3 + XWRE(i, j) = XK(i); + else + XWRE(i, j) = Theta_LL(i, j) + Theta_s(i) * (Theta_L(i, j) / Theta_LL(i, j) - 1); + end + end + XOLD(i) = EX; + return + end + end +end