diff --git a/CHANGELOG.md b/CHANGELOG.md index a6443d4c..37e6f30e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,8 @@ [197](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/197) - Remove `ObservationPoints` script [200](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/200) +- Remove all global variables from `STEMMUS_SCOPE` script +[201](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/201) **Fixed:** @@ -20,6 +22,7 @@ - issue [#99](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/99) - issue [#100](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/100) - issue [#101](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/101) +- issue [#195](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/195) # [1.3.0](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.3.0) - 22 Jun 2023 diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE index 29138cca..04eb7a8b 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/+equations/Soil_Inertia1.m b/src/+equations/Soil_Inertia1.m index af74ab36..ce3499f7 100644 --- a/src/+equations/Soil_Inertia1.m +++ b/src/+equations/Soil_Inertia1.m @@ -1,5 +1,4 @@ -function [GAM] = Soil_Inertia1(SMC) - global theta_s0 +function [GAM] = Soil_Inertia1(SMC, theta_s0) % soil inertia method by Murray and Verhoef ( % % parameters @@ -32,3 +31,4 @@ Hc = (Hcw * SMC) + (1 - theta_s) * Hcs; GAM = sqrt(lambdas .* Hc); +end diff --git a/src/+equations/calc_rssrbs.m b/src/+equations/calc_rssrbs.m index 6be64a5d..f9738078 100644 --- a/src/+equations/calc_rssrbs.m +++ b/src/+equations/calc_rssrbs.m @@ -1,5 +1,7 @@ -function [rss, rbs] = calc_rssrbs(SMC, LAI, rbs) - global SaturatedMC ResidualMC fieldMC +function [rss, rbs] = calc_rssrbs(SoilProperties, LAI, rbs) + + % an alias + SP = SoilProperties; aa = 3.8; - rss = exp((aa + 4.1) - aa * (SMC - ResidualMC(1)) / (fieldMC(1) - ResidualMC(1))); + rss = exp((aa + 4.1) - aa * (SP.SMC - SP.ResidualMC(1)) / (SP.fieldMC(1) - SP.ResidualMC(1))); rbs = rbs * LAI / 4.3; diff --git a/src/+io/select_input.m b/src/+io/select_input.m index 11100207..b91810e4 100644 --- a/src/+io/select_input.m +++ b/src/+io/select_input.m @@ -77,12 +77,12 @@ %% derived input if options.soil_heat_method == 1 - SoilProperties.GAM = equations.Soil_Inertia1(SoilProperties.SMC); + SoilProperties.GAM = equations.Soil_Inertia1(SoilProperties.SMC, SoilProperties.theta_s0); else SoilProperties.GAM = equations.Soil_Inertia0(SoilProperties.cs, SoilProperties.rhos, SoilProperties.lambdas); end if options.calc_rss_rbs - [SoilProperties.rss, SoilProperties.rbs] = equations.calc_rssrbs(SoilProperties.SMC, canopy.LAI, SoilProperties.rbs); + [SoilProperties.rss, SoilProperties.rbs] = equations.calc_rssrbs(SoilProperties, canopy.LAI, SoilProperties.rbs); end if leafbio.Type diff --git a/src/Max_Rootdepth.m b/src/Max_Rootdepth.m index 05dccc40..2373de6f 100644 --- a/src/Max_Rootdepth.m +++ b/src/Max_Rootdepth.m @@ -3,7 +3,7 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% REFERENCES -function [bbx] = Max_Rootdepth(bbx, NL, KT, TT) +function [bbx] = Max_Rootdepth(bbx) % get model settings ModelSettings = io.getModelSettings(); @@ -18,11 +18,11 @@ Tot_Depth = ModelSettings.Tot_Depth; DeltZ = ModelSettings.DeltZ; Elmn_Lnth = 0; - for ML = 1:NL + for ML = 1:ModelSettings.NL Elmn_Lnth = Elmn_Lnth + DeltZ(ML); if Elmn_Lnth < Tot_Depth - R_depth bbx(ML) = 0; - elseif Elmn_Lnth >= Tot_Depth - R_depth && Elmn_Lnth <= Tot_Depth - 5 % && TT(ML)>0 + elseif Elmn_Lnth >= Tot_Depth - R_depth && Elmn_Lnth <= Tot_Depth - 5 bbx(ML) = 1; else bbx(ML) = 0; diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 60b95488..29238a64 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -23,95 +23,40 @@ pkg load statistics io; end -% Read the configPath file. Due to using MATLAB compiler, we cannot use run(CFG) -global CFG -if isempty(CFG) +% set CFG to a path if it is not defined +if exist('CFG', 'var') == 0 CFG = '../config_file_crib.txt'; end + +% 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); % Prepare forcing and soil data -global SaturatedMC ResidualMC fieldMC theta_s0 [SiteProperties, SoilProperties, TimeProperties] = io.prepareInputData(InputPath); -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 select_input % Load model settings: replacing "run Constants" ModelSettings = io.getModelSettings(); - -global J KT Delt_t NN ML rwuef -global T0 rroot NL DeltZ -NL = ModelSettings.NL; -DeltZ = ModelSettings.DeltZ; -DeltZ_R = ModelSettings.DeltZ_R; -J = ModelSettings.J; -T0 = ModelSettings.T0; -rroot = ModelSettings.rroot; -NIT = ModelSettings.NIT; -KT = ModelSettings.KT; NN = ModelSettings.NN; -ML = ModelSettings.ML; -rwuef = ModelSettings.rwuef; - -% defined as global, and used in other scripts -TIME = 0; % Time of simulation released; -Delt_t = TimeProperties.DELT; % Duration of time step [Unit of second] % load forcing data ForcingData = io.loadForcingData(InputPath, TimeProperties, SoilProperties.fmax, ModelSettings.Tot_Depth); -global LAI_msr -LAI_msr = ForcingData.LAI_msr; % used in Root_properties - -global ND h hh T TT P_g P_gg RWU frac TTT Theta_LLL Theta_LL Theta_UUU -global Theta_III SRT TT_CRIT L_f HR Se W thermal Xaa -global XaT Xah DRHOVT DRHOVh DRHODAt DRHODAz Ks RHODA RHOV L sfactor fluxes - % Get initial values InitialValues = init.defineInitialValues(TimeProperties.Dur_tot); -DhDZ = InitialValues.DhDZ; -frac = InitialValues.frac; T_CRIT = InitialValues.T_CRIT; TT_CRIT = InitialValues.TT_CRIT; -HR = InitialValues.HR; % used in Density_V -RHOV = InitialValues.RHOV; -DRHOVh = InitialValues.DRHOVh; -DRHOVT = InitialValues.DRHOVT; -RHODA = InitialValues.RHODA; -DRHODAt = InitialValues.DRHODAt; -DRHODAz = InitialValues.DRHODAz; -Xaa = InitialValues.Xaa; -XaT = InitialValues.XaT; -Xah = InitialValues.Xah; -L = InitialValues.L; hOLD = InitialValues.hOLD; TOLD = InitialValues.TOLD; SAVE = InitialValues.SAVE; - -global C4 SMC bbx Ta Ts RHOV_s DRHOV_sT -SMC = InitialValues.SMC; -bbx = InitialValues.bbx; -Ta = InitialValues.Ta; -Ts = InitialValues.Ts; -RHOV_s = InitialValues.RHOV_s; -DRHOV_sT = InitialValues.DRHOV_sT; P_gOLD = InitialValues.P_gOLD; %% 1. define Constants Constants = io.define_constants(); -global g Rv RDA Hc Rl -g = Constants.g; -Rv = Constants.Rv; -RDA = Constants.RDA; -Hc = Constants.Hc; RTB = 1000; % initial root total biomass (g m-2) % Rl used in ebal -[Rl, Ztot] = Initial_root_biomass(RTB, ModelSettings.DeltZ_R, rroot, ML, SiteProperties.landcoverClass(1)); +[Rl, Ztot] = Initial_root_biomass(RTB, ModelSettings.DeltZ_R, ModelSettings.rroot, ModelSettings.ML, SiteProperties.landcoverClass(1)); %% 2. simulation options path_input = InputPath; % path of all inputs @@ -220,11 +165,10 @@ [rho, tau, rs] = deal(zeros(nwlP + nwlT, 1)); %% 11. load time series data -Theta_LL = InitialValues.Theta_LL; ScopeParametersNames = fieldnames(ScopeParameters); if options.simulation == 1 vi = ones(length(ScopeParametersNames), 1); - [soil, leafbio, canopy, meteo, angles, xyt] = io.select_input(ScopeParameters, Theta_LL, vi, canopy, options, SiteProperties, SoilProperties); + [soil, leafbio, canopy, meteo, angles, xyt] = io.select_input(ScopeParameters, InitialValues.Theta_LL, vi, canopy, options, SiteProperties, SoilProperties); [ScopeParameters, xyt, canopy] = io.loadTimeSeries(ScopeParameters, leafbio, soil, canopy, meteo, F, xyt, path_input, options, SiteProperties.landcoverClass); else soil = struct; @@ -282,58 +226,25 @@ % Run StartInit [SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten); -%% get variables that are defined global and are used by other scripts -global hd hh_frez POR -global XCAP m n Alpha -global Theta_s Theta_r Theta_f - -% get soil constants -SoilConstants = io.getSoilConstants(); -hd = SoilConstants.hd; - -POR = SoilVariables.POR; -XK = SoilVariables.XK; -XCAP = SoilVariables.XCAP; -Theta_s = VanGenuchten.Theta_s; -Theta_r = VanGenuchten.Theta_r; -Theta_f = VanGenuchten.Theta_f; -Alpha = VanGenuchten.Alpha; -n = VanGenuchten.n; -m = VanGenuchten.m; - -%% these vars are defined as global at the begining of this script -%% because they are both input and output of StartInit -Theta_I = SoilVariables.Theta_I; -Theta_U = SoilVariables.Theta_U; -Theta_L = SoilVariables.Theta_L; +% Set SoilVariables that are used in the loop T = SoilVariables.T; h = SoilVariables.h; TT = SoilVariables.TT; -Ks = SoilVariables.Ks; h_frez = SoilVariables.h_frez; +hh = SoilVariables.hh; +hh_frez = SoilVariables.hh_frez; -%% The boundary condition information settings -BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, landcoverClass(1)); +% get soil constants +SoilConstants = io.getSoilConstants(); -%% get global vars -global NBCT NBCTB DSTOR0 BCTB BCT -NBCT = BoundaryCondition.NBCT; -NBCTB = BoundaryCondition.NBCTB; +%% The boundary condition information settings +BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, SiteProperties.landcoverClass(1)); DSTOR = BoundaryCondition.DSTOR; DSTOR0 = BoundaryCondition.DSTOR0; RS = BoundaryCondition.RS; DSTMAX = BoundaryCondition.DSTMAX; IRPT1 = BoundaryCondition.IRPT1; IRPT2 = BoundaryCondition.IRPT2; -BCTB = BoundaryCondition.BCTB; -BCT = BoundaryCondition.BCT; - -% Outputs of UpdateSoilWaterContent used in step Run the model -hh = SoilVariables.hh; -hh_frez = SoilVariables.hh_frez; - -% 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'); @@ -346,7 +257,10 @@ TIMELAST = 0; % the start of simulation period is from 0mins, while the input data start from 30mins. -kk = 0; % DELT=Delt_t; +TIME = 0; % Time of simulation released; +Delt_t = TimeProperties.DELT; % Duration of time step [Unit of second] +KT = ModelSettings.KT; +kk = 0; TimeStep = []; TEND = TIME + TimeProperties.DELT * TimeProperties.Dur_tot; % Time to be reached at the end of simulation period Delt_t0 = Delt_t; % Duration of last time step @@ -380,7 +294,7 @@ end %%%%% Updating the state variables. %%%%%%%%%%%%%%%%%%%%%%%%%%%% L_f = 0; % ignore Freeze/Thaw, see issue 139 - TT_CRIT(NN) = T0; % unit K + TT_CRIT(NN) = ModelSettings.T0; % unit K hOLD_frez = []; if IRPT1 == 0 && IRPT2 == 0 && SoilVariables.ISFT == 0 for MN = 1:NN @@ -391,7 +305,7 @@ hOLD(MN) = h(MN); h(MN) = hh(MN); - DDhDZ(MN, KT) = DhDZ(MN); + DDhDZ(MN, KT) = InitialValues.DhDZ(MN); if ModelSettings.Thmrlefc == 1 TOLD(MN) = T(MN); T(MN) = TT(MN); @@ -411,7 +325,7 @@ if options.simulation == 0 vi(vmax == telmax) = k; end - [soil, leafbio, canopy, meteo, angles, xyt] = io.select_input(ScopeParameters, Theta_LL, vi, canopy, options, xyt, soil); + [soil, leafbio, canopy, meteo, angles, xyt] = io.select_input(ScopeParameters, SoilVariables.Theta_LL, vi, canopy, options, xyt, soil); if options.simulation ~= 1 fprintf('simulation %i ', k); fprintf('of %i \n', telmax); @@ -488,7 +402,8 @@ case 1 [iter, fluxes, rad, thermal, profiles, soil, RWU, frac] ... = ebal(iter, options, spectral, rad, gap, ... - leafopt, angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t); + leafopt, angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t, ... + Rl, SoilVariables, VanGenuchten, InitialValues); if options.calc_fluor if options.calc_vert_profiles [rad, profiles] = RTMf(spectral, rad, soil, leafopt, canopy, gap, angles, profiles); @@ -563,8 +478,8 @@ SoilVariables.XWRE = updateWettingHistory(SoilVariables, VanGenuchten); end - for ML = 1:NL - DDhDZ(ML, KT) = DhDZ(ML); + for i = 1:ModelSettings.NL + DDhDZ(i, KT) = InitialValues.DhDZ(i); end end if Delt_t ~= Delt_t0 @@ -603,8 +518,6 @@ end end - Ts(KT) = Tss; % Tss is calculated above - % Start the iteration procedure in a time step. SoilVariables.h = h; SoilVariables.hh = hh; @@ -613,8 +526,8 @@ SoilVariables.h_frez = h_frez; SoilVariables.Tss(KT) = Tss; - for KIT = 1:NIT - [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, T0, g, L_f, SoilVariables.TT, NN, hd, ForcingData.Tmin); + for KIT = 1:ModelSettings.NIT + [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, NN, SoilConstants.hd, ForcingData.Tmin); % update inputs for UpdateSoilWaterContent SoilVariables.TT_CRIT = TT_CRIT; SoilVariables.hh_frez = hh_frez; @@ -625,7 +538,7 @@ % see issue 181, item 4 KL_T = InitialValues.KL_T; - [RHOV, DRHOVh, DRHOVT] = Density_V(SoilVariables.TT, SoilVariables.hh, g, Rv, NN); + [RHOV, DRHOVh, DRHOVT] = Density_V(SoilVariables.TT, SoilVariables.hh, Constants.g, Constants.Rv, NN); TransportCoefficient = conductivity.calculateTransportCoefficient(InitialValues, SoilVariables, VanGenuchten, Delt_t); W = TransportCoefficient.W; @@ -633,7 +546,7 @@ [L] = Latent(SoilVariables.TT, NN); % DRHODAt unused! - [Xaa, XaT, Xah, DRHODAt, DRHODAz, RHODA] = Density_DA(SoilVariables.T, RDA, P_g, Rv, DeltZ, SoilVariables.h, SoilVariables.hh, SoilVariables.TT, P_gg, Delt_t, NL, NN, DRHOVT, DRHOVh, RHOV); + [Xaa, XaT, Xah, DRHODAt, DRHODAz, RHODA] = Density_DA(SoilVariables.T, Constants.RDA, P_g, Constants.Rv, ModelSettings.DeltZ, SoilVariables.h, SoilVariables.hh, SoilVariables.TT, P_gg, Delt_t, ModelSettings.NL, NN, DRHOVT, DRHOVh, RHOV); ThermalConductivityCapacity = conductivity.calculateThermalConductivityCapacity(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, DRHOVT, L, RHOV); @@ -702,75 +615,58 @@ KIT; KIT = 0; - [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, T0, g, L_f, SoilVariables.TT, NN, hd, ForcingData.Tmin); + [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, NN, SoilConstants.hd, ForcingData.Tmin); % updates inputs for UpdateSoilWaterContent SoilVariables.TT_CRIT = TT_CRIT; SoilVariables.hh_frez = hh_frez; SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten); - % these can be removed after refactoring codes below - h = SoilVariables.h; - hh = SoilVariables.hh; - T = SoilVariables.T; - TT = SoilVariables.TT; - Se = SoilVariables.Se; - Theta_LL = SoilVariables.Theta_LL; - DTheta_LLh = SoilVariables.DTheta_LLh; - hh_frez = SoilVariables.hh_frez; - h_frez = SoilVariables.h_frez; - Theta_UU = SoilVariables.Theta_UU; - DTheta_UUh = SoilVariables.DTheta_UUh; - Theta_II = SoilVariables.Theta_II; - if IRPT1 == 0 && IRPT2 == 0 if KT % In case last time step is not convergent and needs to be repeated. - MN = 0; - for ML = 1:NL - for ND = 1:2 - MN = ML + ND - 1; - Theta_LLL(ML, ND, KT) = Theta_LL(ML, ND); - Theta_L(ML, ND) = Theta_LL(ML, ND); - Theta_UUU(ML, ND, KT) = Theta_UU(ML, ND); - Theta_U(ML, ND) = Theta_UU(ML, ND); - Theta_III(ML, ND, KT) = Theta_II(ML, ND); - Theta_I(ML, ND) = Theta_II(ML, ND); + for i = 1:ModelSettings.NL + for j = 1:ModelSettings.nD + Theta_LLL(i, j, KT) = SoilVariables.Theta_LL(i, j); + SoilVariables.Theta_L(i, j) = SoilVariables.Theta_LL(i, j); + Theta_UUU(i, j, KT) = SoilVariables.Theta_UU(i, j); + SoilVariables.Theta_U(i, j) = SoilVariables.Theta_UU(i, j); + Theta_III(i, j, KT) = SoilVariables.Theta_II(i, j); + SoilVariables.Theta_I(i, j) = SoilVariables.Theta_II(i, j); end end - % update SoilVariables - SoilVariables.Theta_L = Theta_L; - SoilVariables.Theta_U = Theta_U; - SoilVariables.Theta_I = Theta_I; % replace run ObservationPoints, see issue 101 Sim_Theta_U(KT, 1:length(monitorDepthSoilMoisture)) = Theta_UUU(monitorDepthSoilMoisture, 1, KT); Sim_Temp(KT, 1:length(monitorDepthTemperature)) = TTT(monitorDepthTemperature, KT); end if (TEND - TIME) < 1E-3 - for MN = 1:NN - hOLD(MN) = h(MN); - h(MN) = hh(MN); + for i = 1:NN + hOLD(i) = SoilVariables.h(i); + SoilVariables.h(i) = SoilVariables.hh(i); if ModelSettings.Thmrlefc == 1 - TOLD(MN) = T(MN); - T(MN) = TT(MN); - TTT(MN, KT) = TT(MN); - TOLD_CRIT(MN) = T_CRIT(MN); - T_CRIT(MN) = TT_CRIT(MN); - hOLD_frez(MN) = h_frez(MN); - h_frez(MN) = hh_frez(MN); + TOLD(i) = SoilVariables.T(i); + SoilVariables.T(i) = SoilVariables.TT(i); + TTT(i, KT) = SoilVariables.TT(i); + TOLD_CRIT(i) = T_CRIT(i); + T_CRIT(i) = TT_CRIT(i); + hOLD_frez(i) = SoilVariables.h_frez(i); + SoilVariables.h_frez(i) = SoilVariables.hh_frez(i); end if ModelSettings.Soilairefc == 1 - P_gOLD(MN) = P_g(MN); - P_g(MN) = P_gg(MN); + P_gOLD(i) = P_g(i); + P_g(i) = P_gg(i); end end end end - % Update SoilVariables - SoilVariables.h = h; - SoilVariables.T = T; - SoilVariables.h_frez = h_frez; + % set SoilVariables for the rest of the loop + h = SoilVariables.h; + hh = SoilVariables.hh; + T = SoilVariables.T; + TT = SoilVariables.TT; + hh_frez = SoilVariables.hh_frez; + h_frez = SoilVariables.h_frez; kk = k; diff --git a/src/biochemical.m b/src/biochemical.m index 4f5dcaf0..33e3aaa5 100644 --- a/src/biochemical.m +++ b/src/biochemical.m @@ -1,6 +1,5 @@ -function biochem_out = biochemical(biochem_in, Ci_input) - % - global sfactor KT TT +function biochem_out = biochemical(biochem_in, sfactor, Ci_input) + if isnan(sfactor) sfactor = 1; end diff --git a/src/ebal.m b/src/ebal.m index 4899d72e..a257dfbb 100644 --- a/src/ebal.m +++ b/src/ebal.m @@ -1,94 +1,97 @@ function [iter, fluxes, rad, thermal, profiles, soil, RWU, frac] ... = ebal(iter, options, spectral, rad, gap, leafopt, ... - angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t) - global Rl DeltZ Ks Theta_s Theta_r Theta_LL bbx NL KT sfactor PSItot sfactortot Theta_f - global m n Alpha TT - % function ebal.m calculates the energy balance of a vegetated surface - % - % authors: Christiaan van der Tol (tol@itc.nl) - % Joris Timmermans (j_timmermans@itc.nl) - % date 26 Nov 2007 (CvdT) - % updates 29 Jan 2008 (JT & CvdT) converted into a function - % 11 Feb 2008 (JT & CvdT) improved soil heat flux and temperature calculation - % 14 Feb 2008 (JT) changed h in to hc (as h=Avogadro`s constant) - % 31 Jul 2008 (CvdT) Included Pntot in output - % 19 Sep 2008 (CvdT) Converted F0 and F1 from units per aPAR into units per iPAR - % 07 Nov 2008 (CvdT) Changed layout - % 18 Sep 2012 (CvdT) Changed Oc, Cc, ec - % Feb 2012 (WV) introduced structures for variables - % Sep 2013 (JV, CvT) introduced additional biochemical model - % - % parent: master.m (script) - % uses: - % RTMt_sb.m, RTMt_planck.m (optional), RTMf.m (optional) - % resistances.m - % heatfluxes.m - % biochemical.m - % soil_respiration.m - % - % Table of contents of the function - % - % 1. Initialisations for the iteration loop - % intial values are attributed to variables - % 2. Energy balance iteration loop - % iteration between thermal RTM and surface fluxes - % 3. Write warnings whenever the energy balance did not close - % 4. Calculate vertical profiles (optional) - % 5. Calculate spectrally integrated energy, water and CO2 fluxes - % - % The energy balance iteration loop works as follows: - % - % RTMo More or less the classic SAIL model for Radiative - % Transfer of sun and sky light (no emission by the vegetation) - % While continue Here an iteration loop starts to close the energy - % balance, i.e. to match the micro-meteorological model - % and the radiative transfer model - % RTMt_sb A numerical Radiative Transfer Model for thermal - % radiation emitted by the vegetation - % resistances Calculates aerodynamic and boundary layer resistances - % of vegetation and soil (the micro-meteorological model) - % biochemical Calculates photosynthesis, fluorescence and stomatal - % resistance of leaves (or biochemical_MD12: alternative) - % heatfluxes Calculates sensible and latent heat flux of soil and - % vegetation - % Next soil heat flux is calculated, the energy balance - % is evaluated, and soil and leaf temperatures adjusted - % to force energy balance closure - % end {while continue} - % - % meanleaf Integrates the fluxes over all leaf inclinations - % azimuth angles and layers, integrates over the spectrum - % - % usage: - % [iter,fluxes,rad,profiles,thermal] ... - % = ebal(iter,options,spectral,rad,gap,leafopt, ... - % angles,meteo,soil,canopy,leafbio) - % - % The input and output are structures. These structures are further - % specified in a readme file. - % - % Input: - % - % iter numerical parameters used in the iteration for energy balance closure - % options calculation options - % spectral spectral resolutions and wavelengths - % rad incident radiation - % gap probabilities of direct light penetration and viewing - % leafopt leaf optical properties - % angles viewing and observation angles - % soil soil properties - % canopy canopy properties - % leafbio leaf biochemical parameters - % - % Output: - % - % iter numerical parameters used in the iteration for energy balance closure - % fluxes energy balance, turbulent, and CO2 fluxes - % rad radiation spectra - % profiles vertical profiles of fluxes - % thermal temperatures, aerodynamic resistances and friction velocity + angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t, ... + Rl, SoilVariables, VanGenuchten, InitialValues) + + %{ + function ebal.m calculates the energy balance of a vegetated surface + + authors: Christiaan van der Tol (tol@itc.nl) + Joris Timmermans (j_timmermans@itc.nl) + date 26 Nov 2007 (CvdT) + updates 29 Jan 2008 (JT & CvdT) converted into a function + 11 Feb 2008 (JT & CvdT) improved soil heat flux and temperature calculation + 14 Feb 2008 (JT) changed h in to hc (as h=Avogadro`s constant) + 31 Jul 2008 (CvdT) Included Pntot in output + 19 Sep 2008 (CvdT) Converted F0 and F1 from units per aPAR into units per iPAR + 07 Nov 2008 (CvdT) Changed layout + 18 Sep 2012 (CvdT) Changed Oc, Cc, ec + Feb 2012 (WV) introduced structures for variables + Sep 2013 (JV, CvT) introduced additional biochemical model + + parent: master.m (script) + uses: + RTMt_sb.m, RTMt_planck.m (optional), RTMf.m (optional) + resistances.m + heatfluxes.m + biochemical.m + soil_respiration.m + + Table of contents of the function + + 1. Initialisations for the iteration loop + intial values are attributed to variables + 2. Energy balance iteration loop + iteration between thermal RTM and surface fluxes + 3. Write warnings whenever the energy balance did not close + 4. Calculate vertical profiles (optional) + 5. Calculate spectrally integrated energy, water and CO2 fluxes + + The energy balance iteration loop works as follows: + + RTMo More or less the classic SAIL model for Radiative + Transfer of sun and sky light (no emission by the vegetation) + While continue Here an iteration loop starts to close the energy + balance, i.e. to match the micro-meteorological model + and the radiative transfer model + RTMt_sb A numerical Radiative Transfer Model for thermal + radiation emitted by the vegetation + resistances Calculates aerodynamic and boundary layer resistances + of vegetation and soil (the micro-meteorological model) + biochemical Calculates photosynthesis, fluorescence and stomatal + resistance of leaves (or biochemical_MD12: alternative) + heatfluxes Calculates sensible and latent heat flux of soil and + vegetation + Next soil heat flux is calculated, the energy balance + is evaluated, and soil and leaf temperatures adjusted + to force energy balance closure + end {while continue} + + meanleaf Integrates the fluxes over all leaf inclinations + azimuth angles and layers, integrates over the spectrum + + usage: + [iter,fluxes,rad,profiles,thermal] ... + = ebal(iter,options,spectral,rad,gap,leafopt, ... + angles,meteo,soil,canopy,leafbio) + + The input and output are structures. These structures are further + specified in a readme file. + + Input: + + iter numerical parameters used in the iteration for energy balance closure + options calculation options + spectral spectral resolutions and wavelengths + rad incident radiation + gap probabilities of direct light penetration and viewing + leafopt leaf optical properties + angles viewing and observation angles + soil soil properties + canopy canopy properties + leafbio leaf biochemical parameters + + Output: + + iter numerical parameters used in the iteration for energy balance closure + fluxes energy balance, turbulent, and CO2 fluxes + rad radiation spectra + profiles vertical profiles of fluxes + thermal temperatures, aerodynamic resistances and friction velocity + %} %% 1. initialisations and other preparations for the iteration loop + ModelSettings = io.getModelSettings(); counter = 0; % Iteration counter of ebal maxit = iter.maxit; @@ -167,11 +170,11 @@ LAI = canopy.LAI; PSI = 0; - % [bbx]=Max_Rootdepth(bbx,TIME,NL,KT); - [bbx] = Max_Rootdepth(bbx, NL, KT, TT); - [PSIs, rsss, rrr, rxx] = calc_rsoil(Rl, DeltZ, Ks, Theta_s, Theta_r, Theta_LL, bbx, m, n, Alpha); - [sfactor] = calc_sfactor(Rl, Theta_s, Theta_r, Theta_LL, bbx, Ta, Theta_f); - PSIss = PSIs(NL, 1); + + [bbx] = Max_Rootdepth(InitialValues.bbx); + [PSIs, rsss, rrr, rxx] = calc_rsoil(Rl, ModelSettings.DeltZ, SoilVariables.Ks, VanGenuchten.Theta_s, VanGenuchten.Theta_r, SoilVariables.Theta_LL, bbx, VanGenuchten.m, VanGenuchten.n, VanGenuchten.Alpha); + [sfactor] = calc_sfactor(Rl, VanGenuchten.Theta_s, VanGenuchten.Theta_r, SoilVariables.Theta_LL, bbx, Ta, VanGenuchten.Theta_f); + PSIss = PSIs(ModelSettings.NL, 1); %% 2. Energy balance iteration loop % 'Energy balance loop (Energy balance and radiative transfer) @@ -234,14 +237,13 @@ biochem_in.Rdparam = leafbio.Rdparam; if options.Fluorescence_model == 2 % specific for the v.Caemmerer-Magnani model - b = @biochemical_MD12; biochem_in.Tyear = leafbio.Tyear; biochem_in.beta = leafbio.beta; biochem_in.qLs = leafbio.qLs; biochem_in.NPQs = leafbio.kNPQs; biochem_in.stressfactor = leafbio.stressfactor; else - b = @biochemical; % specific for Berry-v.d.Tol model + % specific for Berry-v.d.Tol model biochem_in.tempcor = options.apply_T_corr; biochem_in.Tparams = leafbio.Tparam; biochem_in.stressfactor = SMCsf; @@ -254,9 +256,15 @@ biochem_in.Cs = Cch; biochem_in.Q = rad.Pnh_Cab * 1E6; - biochem_out = b(biochem_in); + if options.Fluorescence_model == 2 + biochem_out = biochemical_MD12(biochem_in); + else + Ci_input = []; + biochem_out = biochemical(biochem_in, sfactor, Ci_input); + end + Ah = biochem_out.A; - Ahh = biochem_out.Ag; + Ahh = biochem_out.Ag; Cih = biochem_out.Ci; Fh = biochem_out.eta; rcwh = biochem_out.rcw; @@ -270,7 +278,12 @@ biochem_in.Cs = Ccu; biochem_in.Q = rad.Pnu_Cab * 1E6; - biochem_out = b(biochem_in); + if options.Fluorescence_model == 2 + biochem_out = biochemical_MD12(biochem_in); + else + Ci_input = []; + biochem_out = biochemical(biochem_in, sfactor, Ci_input); + end Au = biochem_out.A; % Ag? or A? Auu = biochem_out.Ag; % GPP calculation. @@ -335,7 +348,7 @@ end PSI = (PSI + PSI1) / 2; end - PSItot(KT) = PSI; + %%%%%%% if SoilHeatMethod == 2 G = 0.30 * Rns; @@ -553,3 +566,4 @@ % function Tnew = update(Told, Wc, innovation) % Tnew = Wc.*innovation + (1-Wc).*Told; % return +end