From d4f90fd1a3da5b505cd3dcc150bdfad5bf6c227c Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Wed, 18 Sep 2024 10:53:17 +0200 Subject: [PATCH 01/14] Modifiy nonmean() function to mean() function --- config_file_crib.txt | 12 ++++++------ src/+init/setBoundaryCondition.m | 4 ++-- src/+io/loadSoilInitialValues.m | 4 ++-- src/STEMMUS_SCOPE.m | 2 ++ 4 files changed, 12 insertions(+), 10 deletions(-) diff --git a/config_file_crib.txt b/config_file_crib.txt index 45f5c81d..1859ad40 100644 --- a/config_file_crib.txt +++ b/config_file_crib.txt @@ -1,7 +1,7 @@ -WorkDir=/home/jovyan/private/ +WorkDir=/home/jovyan/private/Test/ SoilPropertyPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/SoilProperty/ ForcingPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/Plumber2_data/ -Location=AU-DaS +Location=US-Bo1 directional=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/directional/ fluspect_parameters=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/fluspect_parameters/ leafangles=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/leafangles/ @@ -9,7 +9,7 @@ radiationdata=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/radiationdata/ soil_spectrum=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/soil_spectrum/ input_data=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/input_data.xlsx InitialConditionPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/Initial_condition/ -StartTime=2001-01-01T00:00 -EndTime=2001-01-02T00:00 -InputPath= -OutputPath= \ No newline at end of file +StartTime=2001-04-22T00:00 +EndTime=2001-04-23T00:00 +InputPath=../../input/US-Bo1_2024-09-18-1011/ +OutputPath=../../output/ diff --git a/src/+init/setBoundaryCondition.m b/src/+init/setBoundaryCondition.m index f65b0c84..fd5a8b2d 100644 --- a/src/+init/setBoundaryCondition.m +++ b/src/+init/setBoundaryCondition.m @@ -82,10 +82,10 @@ % 3 --Zero temperature gradient; NBCTB = 1; - if nanmean(Ta_msr) < 0 + if mean(Ta_msr,'omitnan') < 0 BCTB = 0; % 9 8.1 else - BCTB = nanmean(Ta_msr); + BCTB = mean(Ta_msr,'omitnan'); end end if ModelSettings.Soilairefc == 1 diff --git a/src/+io/loadSoilInitialValues.m b/src/+io/loadSoilInitialValues.m index 94aaa06f..ae5427a4 100644 --- a/src/+io/loadSoilInitialValues.m +++ b/src/+io/loadSoilInitialValues.m @@ -68,10 +68,10 @@ InitT6 = 0; Tss = InitT0; end - if nanmean(Ta_msr) < 0 + if mean(Ta_msr,'omitnan') < 0 BtmT = 0; % 9 8.1 else - BtmT = nanmean(Ta_msr); + BtmT = mean(Ta_msr,'omitnan'); end if InitX0 > SaturatedMC(1) || InitX1 > SaturatedMC(1) || InitX2 > SaturatedMC(2) || ... InitX3 > SaturatedMC(3) || InitX4 > SaturatedMC(4) || InitX5 > SaturatedMC(5) || InitX6 > SaturatedMC(6) diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 351cbd24..1e8c2db2 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -374,6 +374,8 @@ % Will do one timestep in "update mode", and run until the end if in "full run" mode. while KT < endTime KT = KT + 1; % Counting Number of timesteps + fprintf(strcat('\n KT = ',num2str(KT),' \r')); + 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. end From 125e9f1763f3695fa44f0022c16079ddb55c875d Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Wed, 18 Sep 2024 11:24:34 +0200 Subject: [PATCH 02/14] Add a function to read WOFOST plant growth parameters from CropD.crp file --- src/+wofost/WofostRead.m | 52 ++++++++++++++++++++++++++++++++++++++++ src/STEMMUS_SCOPE.m | 12 +++++++--- 2 files changed, 61 insertions(+), 3 deletions(-) create mode 100644 src/+wofost/WofostRead.m diff --git a/src/+wofost/WofostRead.m b/src/+wofost/WofostRead.m new file mode 100644 index 00000000..caa7b694 --- /dev/null +++ b/src/+wofost/WofostRead.m @@ -0,0 +1,52 @@ +function [wofost] = WofostRead(path_input) +%INPUT WOFOST PARAMETERS +%This module is added by Danyang Yu, which loads plant growth parameters +%from file "CropD.crp" in the input folder. + +% Read lines from the file +PlantGrowthFile = [path_input,'plant_growth/CropD.crp']; +fid = fopen(PlantGrowthFile); + +% process it line by line +while true + tline = fgetl(fid); + + % remove empty and note line + if ~isempty(tline) & tline(1) ~= '*' + % judge whether the string contain the table value + if contains(tline,'=') + s = regexp(tline,'=','split'); + vname = strtrim(char(s(1))); + % assign value + if ~isempty(char(s(2))) + values = regexp(char(s(2)),'!','split'); + value = str2num(strtrim(char(values(1)))); + wofost.(vname) = value; % save parameters in wofost struct + end + else + % save tabel value + i = 1; + table_value = []; + while ~isempty(tline) & tline(1) ~= '*' + s = strsplit(strtrim(tline)); + table_value(i,1) = str2num(strtrim(char(s(1)))); + table_value(i,2) = str2num(strtrim(char(s(2)))); + i = i + 1; + tline = fgetl(fid); +% disp(tline); + end + + wofost.(vname) = table_value; + end + end + + %end of file + if contains(tline,'* End of .crp file !') + break; + end +% disp(tline); +end + + +end + diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 1e8c2db2..651e19e2 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -198,7 +198,13 @@ [rho, tau, rs] = deal(zeros(nwlP + nwlT, 1)); - %% 11. load time series data + %% 11. Define plant growth parameters + if options.calc_vegetation_dynamic == 1 + wofostpar = wofost.WofostRead(path_input); + crop_output = zeros(TimeProperties.Dur_tot,12); + end + + %% 12. load time series data ScopeParametersNames = fieldnames(ScopeParameters); if options.simulation == 1 vi = ones(length(ScopeParametersNames), 1); @@ -208,7 +214,7 @@ soil = struct; end - %% 12. preparations + %% 13. preparations if options.simulation == 1 diff_tmin = abs(xyt.t - xyt.startDOY); diff_tmax = abs(xyt.t - xyt.endDOY); @@ -242,7 +248,7 @@ atmfile = [path_input 'radiationdata/' char(F(4).FileName(1))]; atmo.M = helpers.aggreg(atmfile, spectral.SCOPEspec); - %% 13. create output files and + %% 14. 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. From b289e2d7ce2fa08563177ca737871ec4e25d6a16 Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Wed, 18 Sep 2024 13:59:15 +0200 Subject: [PATCH 03/14] add the module of plant growth process --- src/+wofost/afgen.m | 13 ++ src/+wofost/cropgrowth.m | 286 +++++++++++++++++++++++++++++++++++++++ src/STEMMUS_SCOPE.m | 12 ++ 3 files changed, 311 insertions(+) create mode 100644 src/+wofost/afgen.m create mode 100644 src/+wofost/cropgrowth.m diff --git a/src/+wofost/afgen.m b/src/+wofost/afgen.m new file mode 100644 index 00000000..70557e33 --- /dev/null +++ b/src/+wofost/afgen.m @@ -0,0 +1,13 @@ +function [output] = afgen(table,x) +%AFGEN +% intercept the value in the table +if x < table(1,1) | x > table(end,1) + print('the value of developemnt stage is mistaken'); +else + dvslist = table(:,1).'; + loc = discretize([x],[-Inf dvslist Inf]); + slope = (table(loc,2) - table((loc-1),2))/(table(loc,1)-table((loc-1),1)); + output = table(loc-1,2) + slope*(x - table(loc-1,1)); +end +end + diff --git a/src/+wofost/cropgrowth.m b/src/+wofost/cropgrowth.m new file mode 100644 index 00000000..7fa6894e --- /dev/null +++ b/src/+wofost/cropgrowth.m @@ -0,0 +1,286 @@ +function [crop_output] = cropgrowth(meteo,wofostpar,Anet,WaterStressFactor,xyt,KT) +%CROPGROWTH MODULE +%This module is added by Danyang Yu, which aims to simulate the growth of vegetation +%The plant growth process could refer to "SWAP version 3.2. Theory description and user manual" + +%% 0. global variables +global Dur_tot dvs wrt wlv wst wso dwrt dwlv dwst +global sla lvage lv lasum laiexp lai ph crop_output Anet_sum lai_delta + +%% 1. initilization +if KT == wofostpar.CSTART + % initial value of crop parameters + Tsum = 0; + dvs = 0; + tdwi = wofostpar.TDWI; + laiem = wofostpar.LAIEM; + ph = 0; + + frtb = wofostpar.FRTB; + fltb = wofostpar.FLTB; + fstb = wofostpar.FSTB; + fotb = wofostpar.FOTB; + + fr = wofost.afgen(frtb, dvs); + fl = wofost.afgen(fltb, dvs); + fs = wofost.afgen(fstb, dvs); + fo = wofost.afgen(fotb, dvs); + + ssa = wofostpar.SSA; + spa = wofostpar.SPA; + + % initial state variables of the crop + tadw = (1-fr)*tdwi; + wrt = fr*tdwi; + wlv = fl*tadw; + wst = fs*tadw; + wso = fo*tadw; + + dwrt = 0.0; + dwlv = 0.0; + dwst = 0.0; + + sla = zeros(Dur_tot+1,1); + lvage = zeros(Dur_tot+1,1); + lv = zeros(Dur_tot+1,1); + sla(1) = wofost.afgen(wofostpar.SLATB, dvs); + lv(1) = wlv; + lvage(1) = 0.0; + ilvold = 1; + + lasum = laiem; + laiexp = laiem; + glaiex = 0; + laimax = laiem; + lai = lasum+ssa*wst+spa*wso; + Anet_sum = 0; + lai_delta = 0; + +end + +%% 1. phenological development +delt = wofostpar.TSTEP/24; +tav = max(0, meteo.Ta); +dtsum = wofost.afgen(wofostpar.DTSMTB,tav); + +if dvs <= 1 % vegetative stage + dvs = dvs + dtsum*delt/wofostpar.TSUMEA; +else + dvs = dvs + dtsum*delt/wofostpar.TSUMAM; +end + +%% 2. dry matter increase +frtb = wofostpar.FRTB; +fltb = wofostpar.FLTB; +fstb = wofostpar.FSTB; +fotb = wofostpar.FOTB; + +fr = wofost.afgen(frtb, dvs); +fl = wofost.afgen(fltb, dvs); +fs = wofost.afgen(fstb, dvs); +fo = wofost.afgen(fotb, dvs); + +fcheck = fr+(fl+fs+fo)*(1.0-fr) - 1.0; %check on partitions +if abs(fcheck) > 0.0001 + print('The sum of partitioning factors for leaves, stems and storage organs is not equal to one at development stage'); + return; +end + +cvf = 1.0/((fl/wofostpar.CVL+fs/wofostpar.CVS+fo/wofostpar.CVO)*(1.0-fr)+fr/wofostpar.CVR); +asrc = Anet*30/1000000000*10000*3600*wofostpar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; +dmi = cvf * asrc; + +%% 3. Growth rate by root +% root extension +% rr = min(wofostpar.RDM - wofostpar.RD,wofostpar.RRI); +admi = (1.0-fr)*dmi; +grrt = fr*dmi; +drrt = wrt*wofost.afgen (wofostpar.RDRRTB,dvs)*delt; +gwrt = grrt-drrt; + +%% 4. Growth rate by stem +% growth rate stems +grst = fs*admi; +% death rate stems +drst = wofost.afgen (wofostpar.RDRSTB,dvs)*wst*delt; +% net growth rate stems +gwst = grst-drst; + +% growth of plant height +str = wofost.afgen(wofostpar.STRTB, dvs); +phnew = wst/(wofostpar.STN*wofostpar.STD*pi*str*str)*wofostpar.PHCoeff+wofostpar.PHEM; +if phnew >= ph + ph = phnew; +end + +%% 5. Growth rate by organs +gwso = fo*admi; + +%% 6. Growth rate by leave +% weight of new leaves +grlv = fl*admi; + +% death of leaves due to water stress or high lai +if abs(lai) < 0.1 + dslv1 = 0; +else + sfactor = WaterStressFactor.soil; + dslv1 = wlv*(1.0-sfactor)*wofostpar.PERDL*delt; +end + +laicr = wofostpar.LAICR; +dslv2 = wlv*max(0,0.03*(lai-laicr)/laicr); +dslv = max(dslv1,dslv2); + +% death of leaves due to exceeding life span +% first: leaf death due to water stress or high lai is imposed on array +% until no more leaves have to die or all leaves are gone +rest = dslv; +i1 = KT; +iday = ceil(KT*delt); + +while (rest>lv(i1) && i1 >= 1) + rest = rest - lv(i1); + i1 = i1 -1; +end + +% then: check if some of the remaining leaves are older than span, +% sum their weights +dalv = 0.0; +if (lvage(i1)>wofostpar.SPAN && rest>0 && i1>=1) + dalv = lv(i1) - rest; + rest = 0; + i1 = i1 -1; +end + +while (i1>=1 && lvage(i1)>wofostpar.SPAN) + dalv = dalv + lv(i1); + i1 = i1 - 1; +end + +% final death rate +drlv = dslv+dalv; + +% calculation of specific leaf area in case of exponential growth: +slat = wofost.afgen(wofostpar.SLATB, dvs); +if laiexp < 6 + dteff = max(0,tav-wofostpar.TBASE); % effective temperature + glaiex = laiexp*wofostpar.RGRLAI*dteff*delt; % exponential growth + laiexp = laiexp + glaiex; + + glasol = grlv*slat; % source-limited growth + gla = min(glaiex,glasol); +% gla = max(0,gla); + + slat = gla/grlv; + if isnan(slat) + slat = 0; + end +% if grlv>=0 +% slat = gla/grlv; % the actual sla value +% else +% slat = 0; +% end +end + +% update the information of leave +dslvt = dslv; +i1 = KT; +while (dslvt>0 && i1>=1) % water stress and high lai + if dslvt >= lv(i1) + dslvt = dslvt-lv(i1); + lv(i1) = 0.0; + i1 = i1-1; + else + lv(i1) = lv(i1)-dslvt; + dslvt = 0.0; + end +end + +while (lvage(i1)>=wofostpar.SPAN && i1>=1) % leaves older than span die + lv(i1) = 0; + i1 = i1-1; +end + +ilvold = KT; % oldest class with leaves +fysdel = max (0.0d0,(tav-wofostpar.TBASE)/(35.0-wofostpar.TBASE)); % physiologic ageing of leaves per time step + +for i1 = ilvold: -1: 1 % shifting of contents, updating of physiological age + lv(i1+1) = lv(i1); + sla(i1+1) = sla(i1); + lvage(i1+1) = lvage(i1)+fysdel*delt; +end +ilvold = ilvold+1; + +lv(1) = grlv; +sla(1) = slat; +lvage(1) = 0.0; + +% calculation of new leaf area and weight +lasum = 0.0; +wlv = 0.0; +for i1 = 1:ilvold + lasum = lasum + lv(i1)*sla(i1); + wlv = wlv + lv(i1); +end +lasum = max(0,lasum); +wlv = max(0,wlv); + +%% 7. integrals of the crop +% dry weight of living plant organs +wrt = wrt+gwrt; +wst = wst+gwst; +wso = wso+gwso; + +% total above ground biomass +tadw = wlv+wst+wso; + +% dry weight of dead plant organs (roots,leaves & stems) +dwrt = dwrt+drrt; +dwlv = dwlv+drlv; +dwst = dwst+drst; + +% dry weight of dead and living plant organs +twrt = wrt+dwrt; +twlv = wlv+dwlv; +twst = wst+dwst; +cwdm = twlv+twst+wso; + +% leaf area index +lai = wofostpar.LAIEM+lasum+wofostpar.SSA*wst+wofostpar.SPA*wso; + +Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration +if Anet_sum < 0 && KT>1 % the consumed biomass are the assimilated one at daytime + lai_delta = lai_delta + grlv*slat; + lai = lai - lai_delta; +elseif Anet_sum>=0 && Anet>=0 + Anet_sum = 0; + lai_delta = 0; +end + + +% Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration +% if Anet_sum < 0 && KT>1 +% lai_old = crop_output(KT-1,3); +% lai = max(lai,lai_old); +% elseif Anet_sum>=0 && Anet>=0 +% Anet_sum = 0; +% end + + +%% 8. integrals of the crop +crop_output(KT,1) = xyt.t(KT,1); % Day of the year +crop_output(KT,2) = dvs; % Development of stage +crop_output(KT,3) = lai; % LAI +crop_output(KT,4) = ph; % Plant height +crop_output(KT,5) = sfactor; % Water stress +crop_output(KT,6) = wrt; % Dry matter of root +crop_output(KT,7) = wlv; % Dry matter of leaves +crop_output(KT,8) = wst; % Dry matter of stem +crop_output(KT,9) = wso; % Dry matter of organ +crop_output(KT,10) = dwrt; % Dry matter of leaves +crop_output(KT,11) = dwlv; % Dry matter of stem +crop_output(KT,12) = dwst; % Dry matter of organ +end + + diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 651e19e2..02d497dd 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -200,6 +200,7 @@ %% 11. Define plant growth parameters if options.calc_vegetation_dynamic == 1 + global crop_output wofostpar = wofost.WofostRead(path_input); crop_output = zeros(TimeProperties.Dur_tot,12); end @@ -550,6 +551,17 @@ end end end + + % calculate the plant growth process + if options.calc_vegetation_dynamic == 1 && KT >= wofostpar.CSTART && KT <= wofostpar.CEND % vegetation growth process + Anet = fluxes.Actot; + if isnan(Anet) || Anet < -2 % limit value of Anet + Anet = 0; + fluxes.Actot = Anet; + end + [crop_output] = wofost.cropgrowth(meteo,wofostpar,Anet,WaterStressFactor,xyt,KT); + end + if options.simulation == 2 && telmax > 1 vi = helpers.count(nvars, vi, vmax, 1); end From 063c1d79e5d3e938dab4e77da0bba66baf424f37 Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Wed, 18 Sep 2024 14:22:30 +0200 Subject: [PATCH 04/14] output the results of plant growth simulations --- src/+io/bin_to_csv.m | 9 +++++++++ src/+io/create_output_files_binary.m | 5 +++++ src/+io/output_data_binary.m | 9 ++++++++- src/STEMMUS_SCOPE.m | 2 +- 4 files changed, 23 insertions(+), 2 deletions(-) diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m index 8a4baf89..9c46846f 100644 --- a/src/+io/bin_to_csv.m +++ b/src/+io/bin_to_csv.m @@ -169,6 +169,15 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings) 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); + + % output the vegetation dynamic results ydy + if options.calc_vegetation_dynamic + cropgrowth_names = {'DOY','DVS','LAI',... + 'PH','Sfactor','RootDM','LeafDM','StemDM','OrganDM','RootDeath','LeafDeath','StemDeath',}; + cropgrowth_units = {'day','-','m2/m2', ... + 'cm','-','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha'}; + write_output(cropgrowth_names, cropgrowth_units, fnames.cropgrowth_file, n_col.cropgrowth, ns); + end fclose('all'); diff --git a/src/+io/create_output_files_binary.m b/src/+io/create_output_files_binary.m index 76695055..303e5082 100644 --- a/src/+io/create_output_files_binary.m +++ b/src/+io/create_output_files_binary.m @@ -107,6 +107,11 @@ fnames.spectrum_obsdir_thermal_file = fullfile(Output_dir, 'spectrum_obsdir_thermal.bin'); % spectrum observation direction fnames.spectrum_hemis_thermal_file = fullfile(Output_dir, 'spectrum_hemis_thermal.bin'); % spectrum hemispherically integrated end + + % Create cropgrowth file + if options.calc_vegetation_dynamic + fnames.cropgrowth_file = fullfile(Output_dir,'cropgrowth.bin'); % crop growth simulation ydy + end % Create empty files for appending structfun(@(x) fopen(x, 'w'), fnames, 'UniformOutput', false); diff --git a/src/+io/output_data_binary.m b/src/+io/output_data_binary.m index e560d195..d84b206d 100644 --- a/src/+io/output_data_binary.m +++ b/src/+io/output_data_binary.m @@ -1,7 +1,7 @@ function n_col = output_data_binary(f, k, xyt, rad, canopy, ScopeParameters, vi, vmax, options, fluxes, meteo, iter, thermal, ... spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap, WaterStress, WaterPotential, ... Sim_hh, Sim_qlh, Sim_qlt, Sim_qvh, Sim_qvt, Sim_qla, Sim_qva, Sim_qtot, ... - ForcingData, RS, RWUs, RWUg) + ForcingData, RS, RWUs, RWUg,crop_output) %% OUTPUT DATA % author C. Van der Tol @@ -39,6 +39,13 @@ n_col.Sim_Temp = length(Sim_Temp_out); fwrite(f.Sim_Temp_file, Sim_Temp_out, 'double'); + %% Crop growth + if options.calc_vegetation_dynamic == 1 + cropgrowth_out = crop_output(k,:); + n_col.cropgrowth = length(cropgrowth_out); + fwrite(f.cropgrowth_file,cropgrowth_out,'double'); + end + %% Water stress factor waterStressFactor_out = [k xyt.year(k) xyt.t(k) WaterStress.soil]; n_col.waterStressFactor = length(waterStressFactor_out); diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 02d497dd..4ee801be 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -818,7 +818,7 @@ n_col = io.output_data_binary(file_ids, k, xyt, rad, canopy, ScopeParameters, vi, vmax, options, fluxes, ... meteo, iter, thermal, spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, ... Evap, WaterStressFactor, WaterPotential, Sim_hh, Sim_qlh, Sim_qlt, Sim_qvh, ... - Sim_qvt, Sim_qla, Sim_qva, Sim_qtot, ForcingData, RS, RWUs, RWUg); + Sim_qvt, Sim_qla, Sim_qva, Sim_qtot, ForcingData, RS, RWUs, RWUg, crop_output); fclose("all"); end end From 7ebdf1536d9a9660453c3e3b306a91ec2348ff76 Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Wed, 18 Sep 2024 14:48:19 +0200 Subject: [PATCH 05/14] use simulated LAI and PH to substitute the prescribed one --- src/STEMMUS_SCOPE.m | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 4ee801be..feef9a20 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -429,6 +429,18 @@ if options.simulation == 0 vi(vmax == telmax) = k; end + + % using simulated LAI and PH to substitute prescribed LAI + if options.calc_vegetation_dynamic == 1 && KT >= wofostpar.CSTART && KT <= wofostpar.CEND + if KT == wofostpar.CSTART + ScopeParameters.LAI(KT) = wofostpar.LAIEM; % initial LAI + ScopeParameters.hc(KT) = wofostpar.PHEM; % initial PH + elseif KT > wofostpar.CSTART + ScopeParameters.LAI(KT) = crop_output(KT-1,3); % substitute LAI + ScopeParameters.hc(KT) = crop_output(KT-1,4); % substitute PH + end + end + [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); From 88ea0b2faf7e9f61d77326880e9ec721d88f0cdf Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Wed, 18 Sep 2024 15:59:33 +0200 Subject: [PATCH 06/14] Upload the addtional input files in the example_data folder --- example_data/input_data.xlsx | Bin 0 -> 25226 bytes example_data/input_soilLayThick.csv | 61 ---------- example_data/plant_growth/CropD.crp | 183 ++++++++++++++++++++++++++++ 3 files changed, 183 insertions(+), 61 deletions(-) create mode 100644 example_data/input_data.xlsx delete mode 100644 example_data/input_soilLayThick.csv create mode 100644 example_data/plant_growth/CropD.crp diff --git a/example_data/input_data.xlsx b/example_data/input_data.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..99b631a2e4e96bc9296837f4252bab61e19259d6 GIT binary patch literal 25226 zcmeFYV~}l2v?f~iF4r#Gwr$(CZQHhO+ur3}wr$&fb=k-RO_Plr5wE$tOr+!}&1gJy`Pm8U}y{sPd^5nH$OM2PbXQq5=1c@VQ$*?5ZNd-t)O)*Jjh&!J;1uRJi zQ5e($#?S(0RJi-NsMq1*2dRuXp@tCPd4OdMnEYpfgN`XtQ#YKi)@L|7&XFE?>)PeI zy@|vRGF#5FZg&ULp?(J6U1#?AS2%by6}s*3B^`96sZ)grz9KL41=`f9CH&&CHn7OK z-Cz22TPJeUAJ%o%o`p*<3j7QRJY5kHbT1P3wCBgUOH=a$1izA4rX?ZA#)O5j#b({% zp`iWbPMKoNi6t8Md3JD731bCojpICganKB6?J?*i|{+;x8QhW)w<@mr$3rB6uK3{d%7T|2tx*Bnap>p^heF-aJOrjv#BLJplLiv zz*P`yB_CpkY-7kxG7U|ML2x<(7x27vy>_v%O4pi?d)ST^6H>;N2u;PQzr?KyoZT}Q zsN@$D_0{=|oraCQ6Bq*HkZOG9o?KhbDkz8$2WJ3hAH80%fi`s&2|FVBZlsHN9n*C& zz?G`m(2f>7)7NtwPN`0_B2H~Jy#1GuNp5@*3IGNGnESDRa6c4oRn2(s5(Gs2=|I99>PpP`K;8&4)Vc- zHGaFdPP3yLB?hc}UYL0T>Z1Cr1oc_i22?jzGW>?rxvTC1GApEioA9;5K`{Dvi z#a|7{tCc#eYm{iqZ5q3*+hC=$`eMKUvsl!+#dFIQa^zU2a&qi8Cik`~tBfyrkzhzJ zM>Z!XKV&jV7$#f@CElK8^$%o{NAVd_SZRI>%n6UmjG1x-g#P=0D=1rj3aEF4b&Bb>;yShAs9RgP@z1i{19Y6cCnUsQE~(=Zyz@>~U4Ww|kB`7mF@ zP9(pyl2Z9YnceLS)*jRF@T#1J_H5Xr05-?vY)RBvIqnSOq#xVo`hf3x-#??3KK|eU zKq$YnuLg4Y4I!Q>16o>|*)e3L2}~mJ>QI~f9Q`bWzst)MU;|yF1A9-d^>${#d%M}0 z6v9&m7d){XiGrrwHMLbWY+bgc$X`CUL`X$$ydvtnxper3>@sLw#BO4!bMF2M_;Bp@ zVHz3r#zy8ecx-n+I3lsO@$f=mv-8y96IrK%E9bTE#o_bz^YbP5@ZtIiYCjZ^8Si+; zcWnamayvfr`@e*%)w{YB#*c+N0s#QP0Du7ef$M)5x&I!t|IOe5{wEuTsT}jgx zL3HrL*TFr36Wt=4moDU(&oc@umYv>y@S=p3;YM@3J>G$G)pbRvw}@%zciy}v4D49F z>m+YNlFhSNb8P+)%}m18l(5&wsh!+H2SGGKhD;78@jGeTTUa{v2vfQQ%C5-nDAHI- zj+=yrL?|r*2J`t($#W|zY`AL?&H_3@!Bmn7*p0r-c)8{xr2&E5C|lnfu3N9Z)`}_3 zgQkSM8Vgqupos=#nBY{ z4Q`LBaM`5!8l50WG}{@~AxN~*Nq+TR{Kif^$Lw|AuPnxgUEQWNzy9FaKi~zC5Bd=Ny9GHMQVz82OsYSI5z4CjT%LKk!J$++#ti_t|Ha(yMxdIA`e{7#m?lL$o z?RAVsJhvLDwBig=1i8ORngUB^z}LquTe|!nm#svN6fmeasOYr{Yi2JJ`dcGJ`vF+gTjV_R zeI@TcyG6xy!|7w9VkV``x+T1(j>AJ@9TVRaQVI0Xp6Dz$Q4+i@MF%cM#ilC{fpuO% z!P~~Doa(j`{hq^6=* zQ$b+Jb&MD7Q8T4y7%53eJ;lgAXKPOXVF;@yU}FI%w2=(jO&__Tj#wkaG4yaaR-5Y~ zYkJlxV|bp(Zg9%I5Ij`JSHYwkaX(Ke(jN4hEPmbH04wQQ^ROFWA1_VVee2lcbg+Zy z$DRK1ss9y@B+YdjELNCL6vORKoJnt}b-tx7=d50c004K3NH#Iw-VQm!Y(j@_<~Y=Z zf}#P!h&?NuEnlziuk*w6+FDsS-tVvHtK`^Vbnk-+zP>>QjP(A`ozZ!Vb7 zPTkm0SK80d1@4V+ACI>0qwm-FcRSkx7&wJ%JC2RS*Rl+7@~NYTRTAg6)T{A}NB>3+ zI*FkJFw#aoiIQXCwG)K~uxB`B0Br9w+*cGbcNpE;#Mz zd=D#a%#BIAJkjOGvWZI=jE>^9ySjO{h;VhVW+3P%cvYM;w9ft-H>n(dj;WoPN_IK- z8u1+{??FbhUoj8Nq7lb=h$jC~!aX>HRd9LQt2;tAx`9>OuX$#_@oRCV z(nXriV32k!#|putP7_b;p*NJH@PCRuFrS8!1mrOju$_n!ggrCEVev0uAYeTZA`A{O zio!RGfP`|K{?U z+RB(2uJVvrBh1$eQLMA`1>uj!M4mL}9hev%%jWLb!^U-fc|=j%<^;AyG2gwEDN=@w zt4-~fF}u;9xpRfNP9xm|HG|QLWH=O8@H2p~V}5|>Suk7Ay;D!LVy4d?ILXV)z5nwls#IkF;p6^}JUR^;ZzR7!jB z6TaOCROOFiD;LRnRVLgIlJ$rxT%k3Am@Z4RlBqsH7rJ4|d|*Nw!Xe?ys7Y0%hl=Ix ziIBIf7nIFUERsvhf$*Xb(SPI3cy+~9%_(V|x+XiiR77#45Mkq&crB||0l%4kI}J0) z@XZ%;rw}>&=yP=P@9R5#1ZF8M)*x(FL5>8jhbLPe-=W1r(Q@s77+hb zic&&87_HMjmz{W=y^B5G*a>qOLy3&1G;#Y>M0j_ICv?LFPwRsv{L;xhP=O2lwt(P@`!8`v&{<2PS6aBjb_eqOJ$bC2eyPs546JKD|y5x*>pD zTML47`Bx;Gws9v@t8`1`RuJ|BH)(ie75=~$0<%uz&*!D(Kd_D4Ae~V!04^s7s&%EB zzx9tD%n=GJ>qs#(FaLelv}&L5Qr&<|=`QoR^oqXsuM4_&>ALYSIK9vscTi7@OyCPl z)$pI!1Z%%tldKEd_QQIBy3D|719zVM2mQ_`$%;M?_CuPle+J+FCtf5{ay8>$=a z8YsaxhR_RAlkB~n*qW7tzp#NP(ocap`n)H}bYWXCz*V*Syl%bWpLDsSKKt@bAnc8~ zow%zf7IWe;RxAH5FZeT0H6iAISRInfK6d+G&W4)(){+?1fA^2+GMiZ!s8epm$Lplz z-90AqQ5vph*W%1oRqK!4^)J_UBAB-tBe|3TT|6+d8)lVP5()ovmOu6b)Ia}K=p%Do z;B%&4=xf#Rv_+RjVOAmxZ*IdZ^UBb05cfItiu`te%-ip$s#d*`NBWDC&)+NdJf_Tt zs{5-Lhf$opKcvbVylyQKaE(p;-2t3jR45y}=<!4c|&@Ol{xSrj}2n_AK_STh*0Z zy=N_}!*FF*O0AO_84Q?{BAF!TiiAa+Dx{20>QiPjI+bO(SiM?}Jde`I2ga4={#r_N z7p%z1Y`WJmBC+nw6|ZEaB(>*J(Qn-Ng!sGAWxT<0oE zELu*K`;2D58RRdMCjiu3bi$S8mG^?IC-DXhyXYM`R2PS(!9O4{h2O8pyJ+&%ZfF$ZLIFp>z2OC;EX*OV z*};vFo6{S9rs7rLR2hQ9Q0RrXV{34HL41sB0s&L&NarQGTGp5aJFIcd)huCtJlO|N z$!i&IvW*3^i0Z2qWqBuxkk)s-=vc4N!>cQ|uFC7uf2kO%ItpEGEWUn>Xg^=|qT#7h z{jxknyDCwm8g-7|?A54YXg{wMPo63BOYlWFZGYPx4iz-l*3c^A z!TYtB==&#B`A_%5DpndcfDcBbmf;j|Q0O%t0CJ)~nNo(%0-XdXI|&gSW^`sTD=Qwf# zGm>;*=~v-dU3#SzD?OP@i?%8Lw}haaDPq1)@$X3#_k8TK5HXcBN?K3Lq)|#+E{98; zv;jcm`hg!4h2d~B*grgaCvHdD6HQ%!Ec;)&(n6xrM?yaVT0ZRm(v@cT4_A6UVx8_M z7X|GNp8u|vvTuZmWJ=yiwW%O+`eaxAcVA7YK}($7N^-lKkaH9ah`6fG`N8Wx)|9Dt z6H~fnr>ntQ+&w`7sJbcHLY4UQ>r($A1`$_C6Vh6f?3B^Xqtor(q^Hs564Q?_Zssk4 z1>XtOP=`e?NY)CRI%i=C9l`EYWx^F!Xb}^AvVkd^flx^jTwu3Gf2BhO@v3?Q&1?X` zK2MnMvqwIgcP^?6%rr>Twj;i&L~u~S93{t-j9BD%Y~>P1{Q}x^&aps3_W0Y|Gn`GBa5Lxj#%ogZtyO%ynxo`?n31xuz>cD8cl%svGR=v1ptpRz|;aVN7tU?Bt1Q#COF|U)_epFG4vLTUB&Pv}33Z_}dgkvt*Qna1hhssTF zfMyFIoMUmE+c^s(ec7mpR1{04i+6igI}=D*(|ydKztM0@ShwK&eK1JtnY(UuAWkXTHi3LFLA8jNpnqxV_&~*^+r^ND&lEL-$+~?miw%}dKH$mk*m6E! z(HLrSpWhwZd1~)pUcR`wNsc$oKihq(A1NCJmhn z;XYk(`K0v#RLzY@Nrbra{gK)-NR-Cvf-f&y#y0UU_XWC#U^&WZ6fSSi2$~{iVmzQb z`Bx#E*$p{4=lOyZ+L_~S0k0YKB+sT&s!KvbQ5yWC zM`XE^6;;Yg<2obLGq8v;Pu#GOluRTiDO^!gGH{4BUv|oF*(}SnG;nU>8-n3tFnbRt z^-~EH)_{2gVt(_mIgoYkYDcJrgL@TO?FJf-*tQp%=#*(;S#j+U8bjk^ar=zt1QNT- z)=-mIgnF5h2lAYZd}}|%tQmw$a=Jdo*$Q+qHexZ4OR%WHvZa4hhYP1C;D3Fj99JY7(dGlxrn=M6$)^f6ku)o&df z?4#R#m~PWY3oHWjG$b#6LzM4q>$z!eEEhza{)FoPa{&aUac%JSCyXQd^}p0SjQ^>5 zEcXA=Jg9H50?+kc5-TBW#_I|#uHy}%mclqA1PO=9NrLK6o32D7;&^fb(bedlFTLsY zB9WniDn{T-;`7Q(k_r-eNF^zo9R<9;T~a)XdtBERHB#Ux-oT=`DXzJ{Jq=lFK!0NR z>`Jq9Z(h@h@ZU_X;ivfnVun>mIJ%b^?OQ;#g<_N`rKeHmG@KmpM^C>)%rVG!R>P`8 zrr_YABRs3Nn;TpSr6V=%GS(+sav=&mg{q1YzNiJrXAQnOuA&YpgS~MZ5}+OG+>T$6 z4KrRoDH2oqF-$RX@zQ@W>$dt%U!kLt&N+Xx>Ng~)uZ-6Dn<{X-EeRSS&SjM=n4;4gi44$K)cvFQ!Y zyJigo3ccVBHs!;fsJKvXVh~yJC}wezS04%7!v4>yJJ-U5!u8+;T zK0hvPN<%U>KRmy34%abR=iTIQwHFl=Gd4dLCJJ~+&N}VC(}q8Vc@i`M?g#$^+Z8~>5^m6H+A=(QtmTfw6mO)(wkchkwE}Z5$lh*jUn;EMEoX6VFyNV;x?@kX@6SLv|MXOJ1H`iEYhmI zFB$PEVk`wyO0+*)%c#68vkqS^sr-s8S+#4d#W>Vj^7%=17x(`xH37HXX8;Zl0B|Pr zzf==U|Eh_k^_X?MA2mVwh)s262=jDQ#8zohsH}P>tZd!746=^~8U=k&RB^r47n3Mf zKIc>o61MyOg$~BxxnYT79qc@}z)oNvKp^j^_}3%3@Va(7kUk>I)VVLITsQ`{KY!$P z1gqrzS!CnIBvrc7Hlv!IDCa za?nkRzCXLwAKK4)X}$KWDx=|GnMLjmW3(~kmHh=GOSk`{D85|v#X_oA#^nAmM)c4U z-Psw;GE#YyTDb)TUpCA~Q#m|Pvhx6+h?^zwCxjgBy+$833J3%Xs>Guh$MD3T?nTca z@4i-X3)m`LYz(p8mAn#2d=n)F9G=1p??%x?tA`fT?qmZjK@z*MQ6mJ}P<9+f)@^)R zW%r{=-HaYOY;f|l9=}yaL&-Np-^q%6#s$~zz{azUi>1Jv0-u@orf-_!t0sq@QcVHn zV1iPa5mFCK1_F;fS@&`}m~r92i5hVM!EUqJI$hf>xJr;IRbuh?WIG@x8YO%Qf<1R6 z^ZT-Jf%^>DIVcJEDw&99l`x5mi@k;<{OOgnln}K4B2eApI3O}BAib)(q4JI_*PAH zTUJw<8$4}0Y1i|uW|o`29&&59dA{$C?o)Mn*5D7l1$IwEfnl_gp8fZWO0!LeVRP9_ zO)rEq1BOh}*rtq+)9Xjz4GQ0f3c=s415s%ZjUHa}aUsz6<3P(*+tT zuen_*x9M)wS69+ZU>;`qy?>+x*4|ku+%CDpY40zyADQk&aa)l!#wHpimkZrLAX~J; zPz(ssV4(>`=W@uuL5YZFqk#1#f3-^|qKC!yC7#xQC13@$_78{Z#`alHAc^_J@wtTC zGjouMQvWzP8Zq_M^dqb)@Vwv^y>^r4TXQ=EP60Hd2)ZV7Wm#$RWIM?$D+3`*2U-y+ zF2u{p!!6{MiLUEdW%_ZPjT_{VD|8{gsv0UWe>F(bk*$X)iYtL z@x=uF-k-Pw5l|&B)C)7#V6+5leFN#zhhZg^`INtk@fyi-^!?%C)aT&QhOc%gcx15& z7}b%O8U2)bt3ZALPtlP;-%36tO^o#{>R0UVKJ<~@F1!3?yH~qkULm9;tB4Ln7ekKb z%}?gdNyJ7MY}{uXd=(aGmf_TH@Gng0_>9hx19 z=BHQWg!E;N8!;CZgjli|Uej>LSwJGoF4Ixaa!<;AkuR{M=3D%gqqfz-;51uzSoS4h z2Wo+?FBXFI%84-{`bL&29W^fcv+`Sh9M9&$8Noi7F8~H?`9&nk>R+CqOc`f`-HDD@ z^3(yx5Ot3$WGE(-7W?40Fwva9;6Mf$&5s@ClBD_q4ZQlJ@+YmAFpIh6dyOyy)fbDX zjf9$;7u?C_w>pGZ`1925=k%@ETjOd*7xeQ;YF8!Bl<05y9@F)k@ zeDJ0q-ub!l*U+3cd+Q6StHlX3X)93@>9e2V)*=EPozae-XASWiQErJl*bsYBd5jN? z0aJE67QNlsTzBHZ^f=L2SS*e)+X9r2PLtTVr5T=a8J^x2rm`);IU|~Y9y$K4i3Jy6 zi{$xuXUdsITlway$cFODzEZ@_;FdK2jRWYew>7HHMu<;IY4#Du%V-US8=LAdb67l(M|) zy{9^z@1>_C#<(pSmzKXFSyw8T?FQr>cSYjd2kN#^vVP4Vs)V+X@~5} z%jyh()aNzIBc5q(y(1nakJyCM&EXQ{n78W^3eOH7#}5?szj8V(A&)qOiDWC< z2f${?LphOQoqCtQ@NrL8WrTF(7r)TbLg?qJcM5TG3bh!220T~vlq3bl0j_e~0?3HU z_8<^|Lo|#db}8w`We`@+vz>~Uv>e_7W=cRzwwsYs=P;tE<8!6B1#GW3OeBd9FepH3 zXe9%oLSrmLJ#AZ;XYJ~3IY#aIzT$X>` zXvQ7H=}7Qobfr&MTXS61gJ=|+s^QgViPONeB9cDv%D+b#Jry#s-td{;ZduwvOaDhO z-weTqfj7yA@@+5F?wAaRme>-@PJ#)K`*{@sW8?r^5MRoHtj<(in&lI0kF?EpV0>cJ zq3CC_+uLa9Tpse*-^0f5q#O(=(U!D>`VV*HdGr*{9z+3NV_5B@v9(v$N(-X&$krk> z*c`x4L}M#tsQr3lhI6DS$`&r9~8`9AQ3ZT z{v8X5y53N7<;;iwFy4*iuEkK2b}%IWd-wsFfcSYEEmJS`@-vOD7L5;W8#W_ggYZk} zV_t2buQ5^&mthP@#bvq<_l3XJ$I6ArX`~<_l{7Y68?1-t!6Z;YCWdL8042;>b)0$i z&wX%KMTnuPY9CKQGAN$x^a<`uea=-Vb|p8g1yf5VMh(!5Jb)?li$rfqVwXdA{YdAr z{}mOB9q#alRHE9j?6Ykz=-y&LhG@)EfL7I0tQRhIOdBIWEZfY z zV~M@KeQ>vYl7aC^U8R2Z)C;O0#rbi<2M_(y%qb$YK0Mv@91^N1wmQ*VNDvD?^<@5n zx;{i7+ zlqFhCFhGw7mQ1-IaWPntSi>quX&w*j)QU^O-h2ZAg(?!J1@$=3PX)MSEH;}DJ&Ds2 zFA{zwb(!sbi8E`*J@yS@+z5@;F@?!Sz$c$xm%sF379fg!4LcI&^x!$!UXO@VV#C!R zcASeb5fX`Y;rTqk-%=qTg$Vd;-x$6mQ?6L4Q8I!VSKB%U9o^kbi-e(cLl!qCx?!K0 zR5LEruvTG^x`#Mq18$feHpD<;njtF8X1~P=9MKZ*hG(>|5J$2N?TaB{j*eqoH3rJo zUyAXY4$Fo6Yl(s%QXGJKS!VH7Fg7&+wMD~9-)vIbSZn(==6x74Ycz0!JZjh1v{c)Q zb3uKA37}Znrc=9ox%m4524IBV%);2-eCb`Yc`B;}xPe#I+_f^=VsnJ^fmPMW#+f5q zq%vs;ks0s3u_s`=g5|su6KBhhs;54n*43qRvm}%iFVePV=N(RCvLlu?FKO1RAz%si z@>%i>`fpleAf6ZWj^-j_j}Si-Re@w9+PzkFWHU_fUI84yFC(>Cta5e!VYU#mt34|+ zvL%=smcA*>PUP`q;z{UHwHblt!V%DrD786HcHwueVkVJLCh=tAy-T(~?t52m-}}oy zjXuP0pU1dcy{{xusWd*Z2cxLS$Eoyw$Qlt3SL?04Uw1#5W8KfE`AmRMZU0ajFFrXo?rM!J|bpbZK28|PD{Ji9(TQ!8jlipF^*p2Zx6ug1=q z@?;yXx~<>+5!fW&nHt^DqN-$uurx|1&TYQ>+FJwjA8C9PRXLKl<98zsSf9t^cO6c0;*e$26jm?)&X$&qECQ(=t^Y)T#&LZAwOoe{w;Hdkm@%_^c?ay}Me@x&I zwY+7NJu>!yIvz97agfv#WQZu*g2-p7GU-#a@~puWp?-Cy(_?zP)g=oWI1h@txX))Z zpc+j~q^F6yOqnSnwb3+Ag0~u~G%8dxXVz*=hM7bZl&%RF-|SXc)V0--G%-SC;MF-G zg?LkQO0lQK$G3E5|2ejqt-4+iW*8z=coRGCJ}bs(5Ia3oJ#y1mLPjtK4h>N7&A^Wl_Z5M_IVMkQTf<9QkiIn8*r%&a?!X4~j~FwD`YgykUcWc`DhNa?sdJDhSQcn28cQK4;=L5ssxWATQBxNtI2Ep@f}(1WNt z{~CX3?-o(AU^&*sA}U}oSNgAyV}`(xa2toO^6XgYqs$iY4$eC2luh4BouVsb?_)AE z&}p7|Me4CWV(|B0u^O87Dok5uTJL4lJ5La&zB?a<1SU_6p*(WnA0QT#0kywdx&`1D zpRM-&o&1{)?IEHQwVV2i#8*R#FKn*6q60dA(x%!s!Bn4dGS4ByI!vZKoV2L4*F3f+ z7jykN)tfS`dKO}9k>nNYg|B#qScW~Z?YVU?Hz;f#a6G@I(i5Ht`z+figUAPbZ2S&d zry#ITFCMBb`Rp)C8Y!AmcLe(Mw96;f%ItLd1Zh_fY#MOJzuH(CXCc z2b}->N!xES<_Pblp#2+2BDh$B#Q7jlK}UyVCf)?{cSJna^B^9SxL_S>7HuNa%bjz_ zy1&g_HP^4gy>K{G8qAh0Dbj1tIx+rjaajG&AIqteATuJn8!*lkjy<*HV9vO}eb?Yk z6r0qy75k-U8lw4o##{R;uVLLJT>)HB2Gfe7>zIusrwb14ecgF1(Z&;2DQx7tHx-ne z#@y&=dk6k1X}Bsz%1k#m?Nl=cFe8VjKEiRl@Zwt!g!LpLslf$|_=enJF!Oq_v=VP0UD0^sU~*i< zxP$qkSW^b>v<{rwG*r&8z8oL+X~J&LGD0Yu{KPJ8;PfW8tblwY~Q2&p7SFk?-pgZO6XV1M_8#h!vcB*~2OO`4Tj0x}pW z17|XV%~6T~J-~F5*d@E@PXD`rAus-6i`^^!;*Da-BpXC+lE$N-h z70R=e@>>!>|5_$gU^TCk5@*hUC8&KHLNYk1>UA0uZH!md}C+fEbj(dBIWJ6)%N@I-n9lyyDEG{%x*6Wh%S5b>p*a$Rqq3;f(>ak z8i>QWE^-OSbL-Qr{QWq05S>LS%)B;tRPrxKQ1MkR?`i*>S0B^E)`tqKM<71NqsV|}BP(n4*MJ?JC)2S> zOy}xNM^z2Fh{J74BI`ApmhyvG+$u$*Tn|bVNUdqm41=WEyX5Trf3%MMzfXhaJ}&Kf zp#lKRD*^x@|EEzf(|0g7Qgm`Kw=w{IuRT{Scb2WUr z-tKkG>GpZt$%K%roWIY-xwUo;*Ag_=&lS1o`+CKyDUiw#X_R1k=4M&`u#_wN z^|`HMsSWio`Ea|`MTIKr>B*5*sdVPkBZuyjek4pG`C&@GJ|8JCDmI-=F>iFON9Wd_ zI>xq+sIOyuc)N6O3)do)Jk+aXVU`O6rb$}s>~+_{)y}`I)V)P1t8+bzbnnzlOViNs zw~&Ggzw7;8m}JboYKTb&`Kv@Ii@HJ~sU7WCZOwKBQ>+q!hs&8k@3(=kCWMn3^6N31 zLOkbA8~P=S3?jX}!1y!;ViO&sJ!Yz-Pjt*YMLbdb-%7;VhY@*^GP-FAs(>S#c`org z&8m(a(u6t2?tr%>k~RhQtHN}B?q_9;z+GfvTJq)S+q&T~tqmtNR^>a!nDC^~%t|L6v+);fPK2J8<;Hb3gHR5v&85L@U_d9^h?4^4!p+rE`-n_DrJ@${q%_ z@kt|9jzE&pRMKVmc>m1`B&r|M=||s~Wa39NxyrKH6c>yboH1vgs2`85jb@$_L|K#n` z&Wb6#v%f*g@w4`?)?czC%_&ZPOWj&W!i)FH#@oqLaIx(KL@1ghaWM|D@Y_Xl79i1V z8fQj>P!XsMZm2p|{v2XSE(#zkUAf_Ms62gxVqiyZ9o)foj^&M7XjLe8YN=E z_E3w@g2{u3sDd=YD#b{19Nm3WyW^?roN)je5XbyV&<$jv20CV*+!>8B&ocXC0zC91 z-E>R-G&wQXD4Y6(>7+LbqkGWCB`5_p%|%;jtoBbQaaP2uA-CkEnQY_^=+a`>ETm%7 z_XS0qJsgs1iC_1dg7U>(RC}dP&wn=kGl{-XbYx7Hyq|ZkA_2QInkHon_7@|6=~C?j zh3J}V^E2w#y*4A$ft?C@fO}B&>P;vVLo5GH2)WeX@H$D?Z9a9&JtU~rFqm_To5!qM zljCv2jQ7(C;y*3+^mzaNI>s#^9bHN?sM%ZXLY}_Ly`1Tgq7M~sK>lW??B#3=3>du| zYOo<9ZRg%Rdc^=2B;kjP_-KTxBPo`rZF4!$O7Wkq?rZuxpLYuaLm(EQ8-yy%8ch`_On_ zzJYl{Z?DUzj<0zt2824GjSKraY`hXKfM#cPak~psyVH!{JrjKR=R}o$o?`&c`pFwU zhoGn7oBU@G>Ut&{=5xxELU@bE`EQcoOKvYA5WwN(p}K<5U!-Cgqjpm< zqp;B1Z}B7|6=NC-cU^9z`MsmERiQeszE5a(HrAe!e2A?f2FD;%JnX0|h-HFVYD){y zRc_&iu*%v~j+3k#U8{{oD;+sgtDI{j>a#0k4j9dbjM=9P_F$N#4&>#SuU0uXsimlx zIR?KfgUrpLBnFAK@a-e?07MuRT-(~gJx~2x)9>vRqbp4+m|P2VE*3DYAADBWwJuv) zu7$Gx;72f9ay3*qQCyHb^30Pf7$bzT9V85}sBBv8+=FuZ(vfx23W*&Bb;0l7%+(W2 z!&88dL0)ZOlarsN7}PS*4!_Y$9d3pY6cS@SBD2<|9k@zjb_P4(5434T=p-Md&0YD= zbiVh{BqaVq<;Ws+kEcD7BwGRRfUptf7zlFx)b850v(kX*^rYJ;PcRyma+;U?^t|5y zp^46MzGMy-yQB-FAq{XGlwiSoql!#ZmgaP`2thI{PbhcY6LqRyNdeF|v!SuR9HTCU zNHDA)cD>6|_Mq=~zI~D@f&M^Y(ALSyUKc!a_Z<&L{giV&Whn|pdd3f4TWO4efI5T4 zS3Qy`J|jc~`QB-#U5fj%F`_AE?fa0#7p z{9bi~y*bhjMaXpJp($uOf0eq|3TlO7N?dt?s$t^Ne%{SYgMI#rcN5fM9I7->CG(SP z6d+Gy|2UQQbz9*c>l=_K!wIkG66l&?j&Q%%?}Xc2hOxn>0hhcJ>pd3T18zMCk%h3* zdEDbqjV|(ZClW#C05G$vn#@7tL2VSUh~_-t-r!k4nNb-ciqB|1fgHJfJ$KIJa>GqXkO* z4);c|&oEB*9UZB=!Z0(z7~g>rv7zw^qO9Q(1ONNGB{{G3Fu>LwzcC$f9A?kNs@m18 zdiY7{lOg*cXZRJucspP0)tjs5o$vAQk2#S>552&mqMi9xj06+4PXQKkxfYRwD%h%+rwf3G?9@_Z8*li?G*2eojz7B5`?$`D^Cc$ReI{ zWqXa&J^Bv-2|#>NhvMw*66-@1Tj;haMl|5MSy`HZk%?m zS9FGt&8`$3(`^W5jhww0Npl)ZaSij(p?=2S%lBOqj|_pBD*esX z+9Owok-mocz8?crx@kvWKt9omMqE#&%Yc6;mf6LaiIzC#dw1I}m*o53aAP+sYg_T& zwN=kE{ZD?uwGK_;e;&5uhV;kJYT)&%3`Y`yF<5QNdQPdd_aU@q5e!ane_Vz5*=s2q zks{%tYK%a~_-Ze{pAJ0ga_3%dgj5F~;$ER)iHS4-& zpMCb;=bp3g`*+%?4Ln0hdxFAmYsT7|jFm^9~arC&gUCfrJT@pcCtUmEU|qp$>TtIG)UJ5KjUGR0<9kr}dBzR=rlI;^B+% zrQu(VLa$Uf(8VMSjitZ#g{b+@NrKdo&-9smMN96AgQfe*=MTHYxQrys?8LEgD2@`2VSJpk2=+cksc!I`jynE3zjc=9tKe z1JMZ&7Mas0)R)v(cUq2(^~{*+Yv}`uKJCmX4KQu1ohTzErtInxY;Sn$gwyRXkU{if zq(ry;cCd*0!F0%n>NEs{(e)>?%Su0^S=1I+CTUME`gu)N;xB`P{L_3D61&psmCfsf z?(b2%aEfL%N$=(#eNPqdPDe_OyAlcv^NF8K*Dm&r&Y~Z^SJ_WpAW=d2^;o2cytkbO z1fHm65U14H0DN+K#Y=Vi(C{L#bp zGK|8n({pO;wek_5_%K?$^kY{=A~Nfh`imKfdLOnTv;=81wZR7>b)U(ACLZcB1-kuo zd`|;LkbQS_(-!FlHb&#J>rF-Uq82w#cbdpu4(=(%3b&l)=Dw+POFTGadoxcx&OJbY zUX2FjMKNHFWJT>ae}_AMTyZY6ysMI7-IsKpbCpZfY87L!zH5`XKebld<(poMZl=~X zVNTkcyhC7!*08W*ld_sRlZ!GNC8hTmg-MT*kvk(_3+5Dwj-EN5eh9iI2tu2SPTg0~ z){}9w*Y_`sul25Z>>3Uht6%b~0h&WYKs|4lb&WrIJ|H};L3ZdmFv#tWHQYKcwXU-v z6ci__BkeO8XHv!z;(JG5z6xDzI6rg?RCRnzA}sAU%-civgdLIT)8L!)1CG|5Uc<4| zpXxDQgue@K5^NR|9tmdWL=7n_zhior)=cGL7epekQv|p(6J#7F#4sCXw6!2R5vk`g z=M)R3OENZdVoP&q%3qUeIn7?0j%!;&6wNSKjWLGfe(`#!j-=Fk(lZOvMmHnQN+x$iAs zd>+YmR?o`DeDShlu}7ZndEPEM5t|?q!wk}T5B$JS(~nP~*BtfUsFk`ME<}bwM9ko+ z_zhhIjejyno-cklw=MP^_UEEcU^Q{fq$1&OKiAJyg$7q^2?H0V-1ORkBg*XEvkHmr zOY=*~R~N76Fa@g1f_7H#2fBh8jrh8gLUZ|gYCC-KIN@>wo~L&q&sB7a@H63U8Q#Gx z`2-pV`huhx?lXBU7|6OlVCD$z01TxSV0;t z)i52#V63siNYidq9GS%Z*UlK6*GWyCGLeuCs|35ZWARgA@<9UnY|~YuX%pmGLyXNW zE})vgN7f4{ zm%#ZHbbky4Fs&@q^*H-ElW+-KG8B5bvNo3DtaC}vEks3y;w5s5sg!41;@fz6#BP+5 zZU|!?DoaW}u9nBiYPgTY)+gS`_25qEqA&LQaNVPI!sf8+GAtcCGa9PdYK$oc&`rl@3YMa6Q5 zBAneSFv>>^?W@^+hTlBj?_us*pH!v^9d$Y-*GK{PW`i|O}2`?{IK|9nL8RK)Y zAIV*fj@;QltWJD)K&+FrGu2z~AQ_=rzwG<9+N6nTgp};mzi+hJfCKZhvR#$uzQ=wk zDaoA4+QK{QJ1Zz|SC3bc-@6bZNA9H43)I_0f}utn>9o2v{@a+t+)K;0TH0bQ=8aJ8 zLvTW0R40btXWW3Lv8*7=%8eG6E6d%@XLc@{iQJx*)?Wb9o;Kgieu&L@u zb*kjoNN?DrD08^QBi=cPqLesLyJ~y*DP??l2nT3IwXVl@q{S^WJ*B>g7!SmH)~pZ? z)+esSADXRo>u{p{m{BpJ(s3txs}N^xc#(kB)QP&u`r)&Wc`$R_j!JI3u>vx> z3qHVbPQDY+tk`$y&SL|Dr5}qJ-&BxMcGElo6XPL`AX+Y^p=O!Mwsw#2#$0n-Gv%UX zY^q><3OOwNNlpKK&zAmI2eqYG)Rf$&_-D*{iUD|*qgxa>3x{-Eq2W=ptBE>{_Fv4G zd&EKrxb7+1(mr#WAK5r1^~%wSkqwpg>AFx3Sgbm7z&a4FnLs9;0|Ww-{R{mAV@j(C)=_Hw5!p7=JMR`4}~?E;)|D?)}eO zK!rBl7Yv_BTw!gfKpy4T4`KV@0Nl8(BKrREu+Dipc}=SM*K4!&&(IR3cNeV>O40MP ztR=gRB=}7#QS~h!%xsX=5ALuqqC9vCY{QgGKYgqXN(%Q63y<_ zZ^7PsE1~314!qCDv{H95J_u6H@ax(PUPUQuzGUa=)O*uocCp(To&ay(4a`V31^cc< z(U?wYumR1!YRW7fCs{8zZ?zs#A4W+Zf|W<~2>KIJ?TB|uI9w`Sx{F^WFY9V%Iy+^l zI=`qZXCD5RI~L$SSpmiWOyx6Jr6nlIcKlq(1P^6G1lSoK=ahRT`KyfN_S1Sa8Iq%H z#5e7TH=;D@uPs~?h^>i?CB(|^uU%ZR*dDB2E__)8xDWr-PE-7DN(lhaB_&ByWJ9DV zz(_6&hkcacL`BqCmQb&-l;0!4&A{2qBDtG5qM1M6z3#V2`OvYEZ$2Fp=o1g2Xv%`k z?^MIkHKwF&yP>Vwjpp0jx^~Aw3`4wWrWJEMHU?M=FU5^j)TtMmrj@Xac5+zJA1q@& zjVH=qQ~f@?*_VesxqEPg*mC~oX0*>YP<#`jm|X%@zf)lYogDjoNSEN*AyqINS7))D7fBEQr z^PuBGL(T-u$K4<9b9Tm24a>}I4?0=MWK^-pCghr8eNoP+pW;3QV5ACFf7OJZ+JK>2 z6PcunJ-;<2^5e5DHZH#hq2j>y8~rKfcLElpBg~c|d4N#N{zliBZv58BM(NH7ES-UK z!C53MMzvbx^ai2?yM#wi-n*@Txucf&5D2{gjXkH>wnok&27V-6NdY2{@g+B3X!3rq z?Zp1%i(m~i1nOFEN40q8XXzgFXI!|D^{Wr7DQJ4^$)oBf!2Y29xKA5`{&^4A(sP5< z;&IXRo>fN2!UM(YY&6^;vJufaMzz;~i!qT+)GZJe1;9XL+4#rA*Mve@HE8d;?~hN7 z=&sB6x&#g1S8;W45ew_544oZsITu|QlIOY}m4~a4<}E8clD8iul}cFEmhft;5F4|% z8xhqr70Y41vgsWk<{PQ!+rAp-TAuvb)4@GydwIpljomWB&;dF(%RRo@y9#JQz3T`8 zNSoLy!Qp>n z9nod~d1WW{#7_TJajn=MoZ`K^r%|Q0=>^cApsRbOG2JmPk~B+K$FT_H81D^5CA6<) zhM_Z1nJ~Z1EU{HxGT);9d0{R_0&Ocd<-f?U%3!>=&Ihqf{?Mp9o~U)|_e}&s8#GkV zfNn8FCF8_l*olvAPkHg0o9hi1B-Pf~MRP28!3Y3RenVSa+Zm)-X6XFqh^IvP8I%wi6yOBy-KQvPyFaqJ7Ws(l^^9)4|4 z1|z=*kcm)xdY178xD~(S>g31zS&hAa7yNGax9ZujzF9~mv(uD-n#01}eHRLbY(jI9 z)cX>|*hgP+j!f??@wl{5#+E|ZGU4x`3bXu=aV;99(=79QPs#=C$ekmOY!BR$_yI-X+B~StRz0j;9ZJ+ORK<=t4S@=vtiWLBbwXM;T&LGBc1$*t zw6hA_I{Cbd7Mq9?ZMcCN!m`}vipsN}kz-ecj)tC7EgL4pkXVSeE#}JCWR)5oX-?o) z+wdF>6egbBQHZxoBqUVEN*_%GyyG5d$Tv%=6g;H4KmOq^S~=C#mors7szv2qVoj)^ zsFm2oPU_De#yLYQu^C>)3shJT0D-}LFUPV$o4yi^tSjH!gAr{N;a6<#(4~$e!CSw}BfsE?pG=>hcS}aogo~J^qc$7R|r9{Ho5s z?Q%OY|Hg%$` Date: Tue, 31 Dec 2024 15:31:01 +0100 Subject: [PATCH 07/14] Remove all global variables --- src/+wofost/cropgrowth.m | 57 +++++++++++++++++++++++++++++++--------- src/STEMMUS_SCOPE.m | 5 ++-- 2 files changed, 48 insertions(+), 14 deletions(-) diff --git a/src/+wofost/cropgrowth.m b/src/+wofost/cropgrowth.m index 7fa6894e..805cde69 100644 --- a/src/+wofost/cropgrowth.m +++ b/src/+wofost/cropgrowth.m @@ -1,11 +1,11 @@ -function [crop_output] = cropgrowth(meteo,wofostpar,Anet,WaterStressFactor,xyt,KT) +function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,wofostpar,Anet,WaterStressFactor,xyt,KT,Dur_tot) %CROPGROWTH MODULE %This module is added by Danyang Yu, which aims to simulate the growth of vegetation %The plant growth process could refer to "SWAP version 3.2. Theory description and user manual" %% 0. global variables -global Dur_tot dvs wrt wlv wst wso dwrt dwlv dwst -global sla lvage lv lasum laiexp lai ph crop_output Anet_sum lai_delta +% global sla lvage lv lasum laiexp lai ph crop_output Anet_sum lai_delta +% global crop_output %% 1. initilization if KT == wofostpar.CSTART @@ -55,7 +55,29 @@ lai = lasum+ssa*wst+spa*wso; Anet_sum = 0; lai_delta = 0; - + +else + % Unpack state variables + dvs = state_vars.dvs; + wrt = state_vars.wrt; + wlv = state_vars.wlv; + wst = state_vars.wst; + wso = state_vars.wso; + + dwrt = state_vars.dwrt; + dwlv = state_vars.dwlv; + dwst = state_vars.dwst; + + sla = state_vars.sla; + lvage = state_vars.lvage; + lv = state_vars.lv; + lasum = state_vars.lasum; + laiexp = state_vars.laiexp; + lai = state_vars.lai; + + ph = state_vars.ph; + Anet_sum = state_vars.Anet_sum; + lai_delta = state_vars.lai_delta; end %% 1. phenological development @@ -259,14 +281,6 @@ end -% Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration -% if Anet_sum < 0 && KT>1 -% lai_old = crop_output(KT-1,3); -% lai = max(lai,lai_old); -% elseif Anet_sum>=0 && Anet>=0 -% Anet_sum = 0; -% end - %% 8. integrals of the crop crop_output(KT,1) = xyt.t(KT,1); % Day of the year @@ -281,6 +295,25 @@ crop_output(KT,10) = dwrt; % Dry matter of leaves crop_output(KT,11) = dwlv; % Dry matter of stem crop_output(KT,12) = dwst; % Dry matter of organ + +%% 9. Pack updated state variables +state_vars.dvs = dvs; +state_vars.wrt = wrt; +state_vars.wlv = wlv; +state_vars.wst = wst; +state_vars.wso = wso; +state_vars.dwrt = dwrt; +state_vars.dwlv = dwlv; +state_vars.dwst = dwst; +state_vars.sla = sla; +state_vars.lvage = lvage; +state_vars.lv = lv; +state_vars.lasum = lasum; +state_vars.laiexp = laiexp; +state_vars.lai = lai; +state_vars.ph = ph; +state_vars.Anet_sum = Anet_sum; +state_vars.lai_delta = lai_delta; end diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index feef9a20..d002c44e 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -200,9 +200,10 @@ %% 11. Define plant growth parameters if options.calc_vegetation_dynamic == 1 - global crop_output +% global crop_output wofostpar = wofost.WofostRead(path_input); crop_output = zeros(TimeProperties.Dur_tot,12); + state_vars = struct(); end %% 12. load time series data @@ -571,7 +572,7 @@ Anet = 0; fluxes.Actot = Anet; end - [crop_output] = wofost.cropgrowth(meteo,wofostpar,Anet,WaterStressFactor,xyt,KT); + [crop_output,state_vars] = wofost.cropgrowth(crop_output,state_vars,meteo,wofostpar,Anet,WaterStressFactor,xyt,KT,TimeProperties.Dur_tot); end if options.simulation == 2 && telmax > 1 From 9968666e9fa4a96646b0f62504c1eb4c395c9dd0 Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Fri, 3 Jan 2025 16:02:24 +0100 Subject: [PATCH 08/14] Add the explainations in the code, and resolve the confict --- example_data/input_soilLayThick.csv | 61 +++++ src/+io/bin_to_csv.m | 152 ++++++------ src/+wofost/WofostRead.m | 10 +- src/+wofost/afgen.m | 8 +- src/+wofost/cropgrowth.asv | 352 ++++++++++++++++++++++++++++ src/+wofost/cropgrowth.m | 130 ++++++---- src/STEMMUS_SCOPE.m | 16 +- 7 files changed, 600 insertions(+), 129 deletions(-) create mode 100644 example_data/input_soilLayThick.csv create mode 100644 src/+wofost/cropgrowth.asv diff --git a/example_data/input_soilLayThick.csv b/example_data/input_soilLayThick.csv new file mode 100644 index 00000000..0460d3d9 --- /dev/null +++ b/example_data/input_soilLayThick.csv @@ -0,0 +1,61 @@ +NL,LayThick (cm),RootDepth (cm) +1,1,450 +2,1, +3,1, +4,2, +5,2, +6,2, +7,2, +8,2, +9,2, +10,2, +11,2, +12,2, +13,2, +14,2, +15,2.5, +16,2.5, +17,2.5, +18,2.5, +19,5, +20,5, +21,5, +22,5, +23,5, +24,10, +25,10, +26,10, +27,10, +28,10, +29,10, +30,10, +31,10, +32,10, +33,10, +34,10, +35,10, +36,10, +37,10, +38,10, +39,10, +40,10, +41,15, +42,15, +43,20, +44,20, +45,20, +46,20, +47,20, +48,20, +49,20, +50,20, +51,20, +52,20, +53,20, +54,20, +55,25, +56,25, +57,25, +58,25, +59,25, +60,25, diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m index 9c46846f..f5359b5d 100644 --- a/src/+io/bin_to_csv.m +++ b/src/+io/bin_to_csv.m @@ -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,91 @@ 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({'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); + if options.calc_PSI + write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution, for PSI only'}, {'W m-2 um-1 sr-1'}, ... + fnames.fluorescencePSI_file, n_col.fluorescencePSI, ns, true); + 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); + end % output the vegetation dynamic results ydy if options.calc_vegetation_dynamic diff --git a/src/+wofost/WofostRead.m b/src/+wofost/WofostRead.m index caa7b694..c6d19034 100644 --- a/src/+wofost/WofostRead.m +++ b/src/+wofost/WofostRead.m @@ -1,7 +1,11 @@ function [wofost] = WofostRead(path_input) -%INPUT WOFOST PARAMETERS -%This module is added by Danyang Yu, which loads plant growth parameters -%from file "CropD.crp" in the input folder. +%{ + function WofostRead.m loads plant growth parameters + from file "CropD.crp" in the input folder. + + authors: Danyang Yu (yudanyang123@gmail.com) + date: 2025/01/01 +%} % Read lines from the file PlantGrowthFile = [path_input,'plant_growth/CropD.crp']; diff --git a/src/+wofost/afgen.m b/src/+wofost/afgen.m index 70557e33..7aac03a9 100644 --- a/src/+wofost/afgen.m +++ b/src/+wofost/afgen.m @@ -1,6 +1,10 @@ function [output] = afgen(table,x) -%AFGEN -% intercept the value in the table +%{ + function afgen.m intercept the value in the table + authors: Danyang Yu (yudanyang123@gmail.com) + date: 2025/01/01 +%} + if x < table(1,1) | x > table(end,1) print('the value of developemnt stage is mistaken'); else diff --git a/src/+wofost/cropgrowth.asv b/src/+wofost/cropgrowth.asv new file mode 100644 index 00000000..6e828603 --- /dev/null +++ b/src/+wofost/cropgrowth.asv @@ -0,0 +1,352 @@ +function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,wofostpar,Anet,WaterStressFactor,xyt,KT,Dur_tot) +%{ + function cropgrowth.m simulate the growth of vegetation + + authors: Danyang Yu (yudanyang123@gmail.com) + date: 2025/01/01 + reference: + + Table of contents of the function + 1. Initilize crop growth paramters + 2. Calculate the phenological development + 3. Calculate the growth rate of different plant organs + 4. Output the crop growth variables + + Input: + Tsum Temperature sum from emergence to the simulated day + tdwi Initial total crop dry weight + laiem Leaf area index at emergence + f(r,l,s,o) The fraction of dry matter to root, leaf, stem, storage + tadw Initial total crop dry weight + lvage leaf age + lv save leaf dry matter at different steps + spa Specific pod area + ssa Specific stem area + sla specific leaf area + dmi converted dry matter at step KT + admi above-ground dry matter + grrt growth rate of root, similar for the other organs + drrt death rate of root, similar for the other organs + WofostPar Wofost parameters, the definitions are in CropD.crp + + Output: + dvs Developement of stage + lai leaf area index + ph Plant height + wrt dry matter weight of root + wlv dry matter weight of leaves + wst dry matter weight of stem + wso dry matter weight of storage organ + dwrt dry matter weight of dead root + dwlv dry matter weight of dead leaves + dwst dry matter weight of dead stem +%} + +%% 1. initilization +if KT == wofostpar.CSTART + % initial value of crop parameters + Tsum = 0; + dvs = 0; + tdwi = wofostpar.TDWI; + laiem = wofostpar.LAIEM; + ph = 0; + + frtb = wofostpar.FRTB; + fltb = wofostpar.FLTB; + fstb = wofostpar.FSTB; + fotb = wofostpar.FOTB; + + fr = wofost.afgen(frtb, dvs); + fl = wofost.afgen(fltb, dvs); + fs = wofost.afgen(fstb, dvs); + fo = wofost.afgen(fotb, dvs); + + ssa = wofostpar.SSA; + spa = wofostpar.SPA; + + % initial state variables of the crop + tadw = (1-fr)*tdwi; + wrt = fr*tdwi; + wlv = fl*tadw; + wst = fs*tadw; + wso = fo*tadw; + + dwrt = 0.0; + dwlv = 0.0; + dwst = 0.0; + + sla = zeros(Dur_tot+1,1); + lvage = zeros(Dur_tot+1,1); + lv = zeros(Dur_tot+1,1); + sla(1) = wofost.afgen(wofostpar.SLATB, dvs); + lv(1) = wlv; + lvage(1) = 0.0; + ilvold = 1; + + lasum = laiem; + laiexp = laiem; + glaiex = 0; + laimax = laiem; + lai = lasum+ssa*wst+spa*wso; + Anet_sum = 0; + lai_delta = 0; + +else + % Unpack state variables + dvs = state_vars.dvs; + wrt = state_vars.wrt; + wlv = state_vars.wlv; + wst = state_vars.wst; + wso = state_vars.wso; + + dwrt = state_vars.dwrt; + dwlv = state_vars.dwlv; + dwst = state_vars.dwst; + + sla = state_vars.sla; + lvage = state_vars.lvage; + lv = state_vars.lv; + lasum = state_vars.lasum; + laiexp = state_vars.laiexp; + lai = state_vars.lai; + + ph = state_vars.ph; + Anet_sum = state_vars.Anet_sum; + lai_delta = state_vars.lai_delta; +end + +%% 1. phenological development +delt = wofostpar.TSTEP/24; +tav = max(0, meteo.Ta); +dtsum = wofost.afgen(wofostpar.DTSMTB,tav); + +if dvs <= 1 % vegetative stage + dvs = dvs + dtsum*delt/wofostpar.TSUMEA; +else + dvs = dvs + dtsum*delt/wofostpar.TSUMAM; +end + +%% 2. dry matter increase +frtb = wofostpar.FRTB; +fltb = wofostpar.FLTB; +fstb = wofostpar.FSTB; +fotb = wofostpar.FOTB; + +fr = wofost.afgen(frtb, dvs); +fl = wofost.afgen(fltb, dvs); +fs = wofost.afgen(fstb, dvs); +fo = wofost.afgen(fotb, dvs); + +fcheck = fr+(fl+fs+fo)*(1.0-fr) - 1.0; %check on partitions +if abs(fcheck) > 0.0001 + print('The sum of partitioning factors for leaves, stems and storage organs is not equal to one at development stage'); + return; +end + +cvf = 1.0/((fl/wofostpar.CVL+fs/wofostpar.CVS+fo/wofostpar.CVO)*(1.0-fr)+fr/wofostpar.CVR); +asrc = Anet*30/1000000000*10000*3600*wofostpar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; +dmi = cvf * asrc; + +%% 3. Growth rate by root +% root extension +% rr = min(wofostpar.RDM - wofostpar.RD,wofostpar.RRI); +admi = (1.0-fr)*dmi; +grrt = fr*dmi; +drrt = wrt*wofost.afgen (wofostpar.RDRRTB,dvs)*delt; +gwrt = grrt-drrt; + +%% 4. Growth rate by stem +% growth rate stems +grst = fs*admi; +% death rate stems +drst = wofost.afgen (wofostpar.RDRSTB,dvs)*wst*delt; +% net growth rate stems +gwst = grst-drst; + +% growth of plant height +str = wofost.afgen(wofostpar.STRTB, dvs); +phnew = wst/(wofostpar.STN*wofostpar.STD*pi*str*str)*wofostpar.PHCoeff+wofostpar.PHEM; +if phnew >= ph + ph = phnew; +end + +%% 5. Growth rate by organs +gwso = fo*admi; + +%% 6. Growth rate by leave +% weight of new leaves +grlv = fl*admi; + +% death of leaves due to water stress or high lai +if abs(lai) < 0.1 + dslv1 = 0; +else + sfactor = WaterStressFactor.soil; + dslv1 = wlv*(1.0-sfactor)*wofostpar.PERDL*delt; +end + +laicr = wofostpar.LAICR; +dslv2 = wlv*max(0,0.03*(lai-laicr)/laicr); +dslv = max(dslv1,dslv2); + +% death of leaves due to exceeding life span +% first: leaf death due to water stress or high lai is imposed on array +% until no more leaves have to die or all leaves are gone +rest = dslv; +i1 = KT; +iday = ceil(KT*delt); + +while (rest>lv(i1) && i1 >= 1) + rest = rest - lv(i1); + i1 = i1 -1; +end + +% then: check if some of the remaining leaves are older than span, +% sum their weights +dalv = 0.0; +if (lvage(i1)>wofostpar.SPAN && rest>0 && i1>=1) + dalv = lv(i1) - rest; + rest = 0; + i1 = i1 -1; +end + +while (i1>=1 && lvage(i1)>wofostpar.SPAN) + dalv = dalv + lv(i1); + i1 = i1 - 1; +end + +% final death rate +drlv = dslv+dalv; + +% calculation of specific leaf area in case of exponential growth: +slat = wofost.afgen(wofostpar.SLATB, dvs); +if laiexp < 6 + dteff = max(0,tav-wofostpar.TBASE); % effective temperature + glaiex = laiexp*wofostpar.RGRLAI*dteff*delt; % exponential growth + laiexp = laiexp + glaiex; + + glasol = grlv*slat; % source-limited growth + gla = min(glaiex,glasol); +% gla = max(0,gla); + + slat = gla/grlv; + if isnan(slat) + slat = 0; + end +% if grlv>=0 +% slat = gla/grlv; % the actual sla value +% else +% slat = 0; +% end +end + +% update the information of leave +dslvt = dslv; +i1 = KT; +while (dslvt>0 && i1>=1) % water stress and high lai + if dslvt >= lv(i1) + dslvt = dslvt-lv(i1); + lv(i1) = 0.0; + i1 = i1-1; + else + lv(i1) = lv(i1)-dslvt; + dslvt = 0.0; + end +end + +while (lvage(i1)>=wofostpar.SPAN && i1>=1) % leaves older than span die + lv(i1) = 0; + i1 = i1-1; +end + +ilvold = KT; % oldest class with leaves +fysdel = max (0.0d0,(tav-wofostpar.TBASE)/(35.0-wofostpar.TBASE)); % physiologic ageing of leaves per time step + +for i1 = ilvold: -1: 1 % shifting of contents, updating of physiological age + lv(i1+1) = lv(i1); + sla(i1+1) = sla(i1); + lvage(i1+1) = lvage(i1)+fysdel*delt; +end +ilvold = ilvold+1; + +lv(1) = grlv; +sla(1) = slat; +lvage(1) = 0.0; + +% calculation of new leaf area and weight +lasum = 0.0; +wlv = 0.0; +for i1 = 1:ilvold + lasum = lasum + lv(i1)*sla(i1); + wlv = wlv + lv(i1); +end +lasum = max(0,lasum); +wlv = max(0,wlv); + +%% 7. integrals of the crop +% dry weight of living plant organs +wrt = wrt+gwrt; +wst = wst+gwst; +wso = wso+gwso; + +% total above ground biomass +tadw = wlv+wst+wso; + +% dry weight of dead plant organs (roots,leaves & stems) +dwrt = dwrt+drrt; +dwlv = dwlv+drlv; +dwst = dwst+drst; + +% dry weight of dead and living plant organs +twrt = wrt+dwrt; +twlv = wlv+dwlv; +twst = wst+dwst; +cwdm = twlv+twst+wso; + +% leaf area index +lai = wofostpar.LAIEM+lasum+wofostpar.SSA*wst+wofostpar.SPA*wso; + +Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration +if Anet_sum < 0 && KT>1 % the consumed biomass are the assimilated one at daytime + lai_delta = lai_delta + grlv*slat; + lai = lai - lai_delta; +elseif Anet_sum>=0 && Anet>=0 + Anet_sum = 0; + lai_delta = 0; +end + +%% 8. integrals of the crop +crop_output(KT,1) = xyt.t(KT,1); % Day of the year +crop_output(KT,2) = dvs; % Development of stage +crop_output(KT,3) = lai; % LAI +crop_output(KT,4) = ph; % Plant height +crop_output(KT,5) = sfactor; % Water stress +crop_output(KT,6) = wrt; % Dry matter of root +crop_output(KT,7) = wlv; % Dry matter of leaves +crop_output(KT,8) = wst; % Dry matter of stem +crop_output(KT,9) = wso; % Dry matter of organ +crop_output(KT,10) = dwrt; % Dry matter of deadleaves +crop_output(KT,11) = dwlv; % Dry matter of dead stem +crop_output(KT,12) = dwst; % Dry matter of dead organ + +%% 9. Pack updated state variables +state_vars.dvs = dvs; +state_vars.wrt = wrt; +state_vars.wlv = wlv; +state_vars.wst = wst; +state_vars.wso = wso; +state_vars.dwrt = dwrt; +state_vars.dwlv = dwlv; +state_vars.dwst = dwst; +state_vars.sla = sla; +state_vars.lvage = lvage; +state_vars.lv = lv; +state_vars.lasum = lasum; +state_vars.laiexp = laiexp; +state_vars.lai = lai; +state_vars.ph = ph; +state_vars.Anet_sum = Anet_sum; +state_vars.lai_delta = lai_delta; +end + + diff --git a/src/+wofost/cropgrowth.m b/src/+wofost/cropgrowth.m index 805cde69..4d1dcf67 100644 --- a/src/+wofost/cropgrowth.m +++ b/src/+wofost/cropgrowth.m @@ -1,33 +1,69 @@ -function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,wofostpar,Anet,WaterStressFactor,xyt,KT,Dur_tot) -%CROPGROWTH MODULE -%This module is added by Danyang Yu, which aims to simulate the growth of vegetation -%The plant growth process could refer to "SWAP version 3.2. Theory description and user manual" - -%% 0. global variables -% global sla lvage lv lasum laiexp lai ph crop_output Anet_sum lai_delta -% global crop_output +function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,WofostPar,Anet,WaterStressFactor,xyt,KT,Dur_tot) +%{ + function cropgrowth.m simulate the growth of vegetation + + authors: Danyang Yu (yudanyang123@gmail.com) + date: 2025/01/01 + reference: SWAP version 3.2. Theory description and user manual (Page 147-163) + + Table of contents of the function + 1. Initilize crop growth paramters + 2. Calculate the phenological development + 3. Calculate the growth rate of different plant organs + 4. Output the crop growth variables + + Input: + Tsum Temperature sum from emergence to the simulated day + tdwi Initial total crop dry weight + laiem Leaf area index at emergence + f(r,l,s,o) The fraction of dry matter to root, leaf, stem, storage + tadw Initial total crop dry weight + lvage leaf age + lv save leaf dry matter at different steps + spa Specific pod area + ssa Specific stem area + sla specific leaf area + Anet_sum the sum of net assimilation + dmi converted dry matter at step KT + admi above-ground dry matter + grrt growth rate of root, similar for the other organs + drrt death rate of root, similar for the other organs + WofostPar Wofost parameters, the definitions are in CropD.crp + + Output: + dvs Developement of stage + lai leaf area index + ph Plant height + wrt dry matter weight of root + wlv dry matter weight of leaves + wst dry matter weight of stem + wso dry matter weight of storage organ + dwrt dry matter weight of dead root + dwlv dry matter weight of dead leaves + dwst dry matter weight of dead stem +%} %% 1. initilization -if KT == wofostpar.CSTART +if KT == WofostPar.CSTART % initial value of crop parameters Tsum = 0; dvs = 0; - tdwi = wofostpar.TDWI; - laiem = wofostpar.LAIEM; + tdwi = WofostPar.TDWI; + laiem = WofostPar.LAIEM; ph = 0; - frtb = wofostpar.FRTB; - fltb = wofostpar.FLTB; - fstb = wofostpar.FSTB; - fotb = wofostpar.FOTB; + frtb = WofostPar.FRTB; + fltb = WofostPar.FLTB; + fstb = WofostPar.FSTB; + fotb = WofostPar.FOTB; fr = wofost.afgen(frtb, dvs); fl = wofost.afgen(fltb, dvs); fs = wofost.afgen(fstb, dvs); fo = wofost.afgen(fotb, dvs); - ssa = wofostpar.SSA; - spa = wofostpar.SPA; + ssa = WofostPar.SSA; + spa = WofostPar.SPA; % initial state variables of the crop tadw = (1-fr)*tdwi; @@ -43,7 +79,7 @@ sla = zeros(Dur_tot+1,1); lvage = zeros(Dur_tot+1,1); lv = zeros(Dur_tot+1,1); - sla(1) = wofost.afgen(wofostpar.SLATB, dvs); + sla(1) = wofost.afgen(WofostPar.SLATB, dvs); lv(1) = wlv; lvage(1) = 0.0; ilvold = 1; @@ -81,21 +117,21 @@ end %% 1. phenological development -delt = wofostpar.TSTEP/24; +delt = WofostPar.TSTEP/24; tav = max(0, meteo.Ta); -dtsum = wofost.afgen(wofostpar.DTSMTB,tav); +dtsum = wofost.afgen(WofostPar.DTSMTB,tav); if dvs <= 1 % vegetative stage - dvs = dvs + dtsum*delt/wofostpar.TSUMEA; + dvs = dvs + dtsum*delt/WofostPar.TSUMEA; else - dvs = dvs + dtsum*delt/wofostpar.TSUMAM; + dvs = dvs + dtsum*delt/WofostPar.TSUMAM; end %% 2. dry matter increase -frtb = wofostpar.FRTB; -fltb = wofostpar.FLTB; -fstb = wofostpar.FSTB; -fotb = wofostpar.FOTB; +frtb = WofostPar.FRTB; +fltb = WofostPar.FLTB; +fstb = WofostPar.FSTB; +fotb = WofostPar.FOTB; fr = wofost.afgen(frtb, dvs); fl = wofost.afgen(fltb, dvs); @@ -108,29 +144,29 @@ return; end -cvf = 1.0/((fl/wofostpar.CVL+fs/wofostpar.CVS+fo/wofostpar.CVO)*(1.0-fr)+fr/wofostpar.CVR); -asrc = Anet*30/1000000000*10000*3600*wofostpar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; +cvf = 1.0/((fl/WofostPar.CVL+fs/WofostPar.CVS+fo/WofostPar.CVO)*(1.0-fr)+fr/WofostPar.CVR); +asrc = Anet*30/1000000000*10000*3600*WofostPar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; dmi = cvf * asrc; %% 3. Growth rate by root % root extension -% rr = min(wofostpar.RDM - wofostpar.RD,wofostpar.RRI); +% rr = min(WofostPar.RDM - WofostPar.RD,WofostPar.RRI); admi = (1.0-fr)*dmi; grrt = fr*dmi; -drrt = wrt*wofost.afgen (wofostpar.RDRRTB,dvs)*delt; +drrt = wrt*wofost.afgen (WofostPar.RDRRTB,dvs)*delt; gwrt = grrt-drrt; %% 4. Growth rate by stem % growth rate stems grst = fs*admi; % death rate stems -drst = wofost.afgen (wofostpar.RDRSTB,dvs)*wst*delt; +drst = wofost.afgen (WofostPar.RDRSTB,dvs)*wst*delt; % net growth rate stems gwst = grst-drst; % growth of plant height -str = wofost.afgen(wofostpar.STRTB, dvs); -phnew = wst/(wofostpar.STN*wofostpar.STD*pi*str*str)*wofostpar.PHCoeff+wofostpar.PHEM; +str = wofost.afgen(WofostPar.STRTB, dvs); +phnew = wst/(WofostPar.STN*WofostPar.STD*pi*str*str)*WofostPar.PHCoeff+WofostPar.PHEM; if phnew >= ph ph = phnew; end @@ -147,10 +183,10 @@ dslv1 = 0; else sfactor = WaterStressFactor.soil; - dslv1 = wlv*(1.0-sfactor)*wofostpar.PERDL*delt; + dslv1 = wlv*(1.0-sfactor)*WofostPar.PERDL*delt; end -laicr = wofostpar.LAICR; +laicr = WofostPar.LAICR; dslv2 = wlv*max(0,0.03*(lai-laicr)/laicr); dslv = max(dslv1,dslv2); @@ -169,13 +205,13 @@ % then: check if some of the remaining leaves are older than span, % sum their weights dalv = 0.0; -if (lvage(i1)>wofostpar.SPAN && rest>0 && i1>=1) +if (lvage(i1)>WofostPar.SPAN && rest>0 && i1>=1) dalv = lv(i1) - rest; rest = 0; i1 = i1 -1; end -while (i1>=1 && lvage(i1)>wofostpar.SPAN) +while (i1>=1 && lvage(i1)>WofostPar.SPAN) dalv = dalv + lv(i1); i1 = i1 - 1; end @@ -184,10 +220,10 @@ drlv = dslv+dalv; % calculation of specific leaf area in case of exponential growth: -slat = wofost.afgen(wofostpar.SLATB, dvs); +slat = wofost.afgen(WofostPar.SLATB, dvs); if laiexp < 6 - dteff = max(0,tav-wofostpar.TBASE); % effective temperature - glaiex = laiexp*wofostpar.RGRLAI*dteff*delt; % exponential growth + dteff = max(0,tav-WofostPar.TBASE); % effective temperature + glaiex = laiexp*WofostPar.RGRLAI*dteff*delt; % exponential growth laiexp = laiexp + glaiex; glasol = grlv*slat; % source-limited growth @@ -219,13 +255,13 @@ end end -while (lvage(i1)>=wofostpar.SPAN && i1>=1) % leaves older than span die +while (lvage(i1)>=WofostPar.SPAN && i1>=1) % leaves older than span die lv(i1) = 0; i1 = i1-1; end ilvold = KT; % oldest class with leaves -fysdel = max (0.0d0,(tav-wofostpar.TBASE)/(35.0-wofostpar.TBASE)); % physiologic ageing of leaves per time step +fysdel = max (0.0d0,(tav-WofostPar.TBASE)/(35.0-WofostPar.TBASE)); % physiologic ageing of leaves per time step for i1 = ilvold: -1: 1 % shifting of contents, updating of physiological age lv(i1+1) = lv(i1); @@ -269,7 +305,7 @@ cwdm = twlv+twst+wso; % leaf area index -lai = wofostpar.LAIEM+lasum+wofostpar.SSA*wst+wofostpar.SPA*wso; +lai = WofostPar.LAIEM+lasum+WofostPar.SSA*wst+WofostPar.SPA*wso; Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration if Anet_sum < 0 && KT>1 % the consumed biomass are the assimilated one at daytime @@ -280,8 +316,6 @@ lai_delta = 0; end - - %% 8. integrals of the crop crop_output(KT,1) = xyt.t(KT,1); % Day of the year crop_output(KT,2) = dvs; % Development of stage @@ -292,9 +326,9 @@ crop_output(KT,7) = wlv; % Dry matter of leaves crop_output(KT,8) = wst; % Dry matter of stem crop_output(KT,9) = wso; % Dry matter of organ -crop_output(KT,10) = dwrt; % Dry matter of leaves -crop_output(KT,11) = dwlv; % Dry matter of stem -crop_output(KT,12) = dwst; % Dry matter of organ +crop_output(KT,10) = dwrt; % Dry matter of deadleaves +crop_output(KT,11) = dwlv; % Dry matter of dead stem +crop_output(KT,12) = dwst; % Dry matter of dead organ %% 9. Pack updated state variables state_vars.dvs = dvs; diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index d002c44e..2b44770c 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -201,7 +201,7 @@ %% 11. Define plant growth parameters if options.calc_vegetation_dynamic == 1 % global crop_output - wofostpar = wofost.WofostRead(path_input); + WofostPar = wofost.WofostRead(path_input); crop_output = zeros(TimeProperties.Dur_tot,12); state_vars = struct(); end @@ -432,11 +432,11 @@ end % using simulated LAI and PH to substitute prescribed LAI - if options.calc_vegetation_dynamic == 1 && KT >= wofostpar.CSTART && KT <= wofostpar.CEND - if KT == wofostpar.CSTART - ScopeParameters.LAI(KT) = wofostpar.LAIEM; % initial LAI - ScopeParameters.hc(KT) = wofostpar.PHEM; % initial PH - elseif KT > wofostpar.CSTART + if options.calc_vegetation_dynamic == 1 && KT >= WofostPar.CSTART && KT <= WofostPar.CEND + if KT == WofostPar.CSTART + ScopeParameters.LAI(KT) = WofostPar.LAIEM; % initial LAI + ScopeParameters.hc(KT) = WofostPar.PHEM; % initial PH + elseif KT > WofostPar.CSTART ScopeParameters.LAI(KT) = crop_output(KT-1,3); % substitute LAI ScopeParameters.hc(KT) = crop_output(KT-1,4); % substitute PH end @@ -566,13 +566,13 @@ end % calculate the plant growth process - if options.calc_vegetation_dynamic == 1 && KT >= wofostpar.CSTART && KT <= wofostpar.CEND % vegetation growth process + if options.calc_vegetation_dynamic == 1 && KT >= WofostPar.CSTART && KT <= WofostPar.CEND % vegetation growth process Anet = fluxes.Actot; if isnan(Anet) || Anet < -2 % limit value of Anet Anet = 0; fluxes.Actot = Anet; end - [crop_output,state_vars] = wofost.cropgrowth(crop_output,state_vars,meteo,wofostpar,Anet,WaterStressFactor,xyt,KT,TimeProperties.Dur_tot); + [crop_output,state_vars] = wofost.cropgrowth(crop_output,state_vars,meteo,WofostPar,Anet,WaterStressFactor,xyt,KT,TimeProperties.Dur_tot); end if options.simulation == 2 && telmax > 1 From ee71fef1e3c63ffa7f4a927174a8d76533bcc95d Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Fri, 3 Jan 2025 16:16:11 +0100 Subject: [PATCH 09/14] Solve the conficts --- src/+io/bin_to_csv.m | 26 +++++++++----------------- 1 file changed, 9 insertions(+), 17 deletions(-) diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m index f5359b5d..44a21feb 100644 --- a/src/+io/bin_to_csv.m +++ b/src/+io/bin_to_csv.m @@ -102,6 +102,15 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings) end end + % output the vegetation dynamic results ydy + if options.calc_vegetation_dynamic + cropgrowth_names = {'DOY','DVS','LAI',... + 'PH','Sfactor','RootDM','LeafDM','StemDM','OrganDM','RootDeath','LeafDeath','StemDeath',}; + cropgrowth_units = {'day','-','m2/m2', ... + 'cm','-','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha'}; + write_output(cropgrowth_names, cropgrowth_units, fnames.cropgrowth_file, n_col.cropgrowth, ns); + end + % Optional for large output files if FullCSVfiles %% Spectrum (added on 19 September 2008) @@ -161,14 +170,6 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings) 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); - if options.calc_PSI - write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution, for PSI only'}, {'W m-2 um-1 sr-1'}, ... - fnames.fluorescencePSI_file, n_col.fluorescencePSI, ns, true); - 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'}, ... @@ -185,15 +186,6 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings) 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 - - % output the vegetation dynamic results ydy - if options.calc_vegetation_dynamic - cropgrowth_names = {'DOY','DVS','LAI',... - 'PH','Sfactor','RootDM','LeafDM','StemDM','OrganDM','RootDeath','LeafDeath','StemDeath',}; - cropgrowth_units = {'day','-','m2/m2', ... - 'cm','-','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha'}; - write_output(cropgrowth_names, cropgrowth_units, fnames.cropgrowth_file, n_col.cropgrowth, ns); - end fclose('all'); From 9d7f7951565d574739af990f3fe177d19dae883b Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Sun, 5 Jan 2025 17:06:15 +0100 Subject: [PATCH 10/14] Solve the conflict for io.bin_to_csv function.m --- src/+io/bin_to_csv.m | 27 ++- src/+io/read_config.m | 9 +- src/+wofost/cropgrowth.asv | 352 ------------------------------------- src/STEMMUS_SCOPE.m | 4 +- 4 files changed, 23 insertions(+), 369 deletions(-) delete mode 100644 src/+wofost/cropgrowth.asv diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m index 44a21feb..51613bea 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,6 +91,15 @@ 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); + %% output the vegetation dynamic results Danyang Yu + if options.calc_vegetation_dynamic + cropgrowth_names = {'DOY','DVS','LAI',... + 'PH','Sfactor','RootDM','LeafDM','StemDM','OrganDM','RootDeath','LeafDeath','StemDeath',}; + cropgrowth_units = {'day','-','m2/m2', ... + 'cm','-','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha'}; + write_output(cropgrowth_names, cropgrowth_units, fnames.cropgrowth_file, n_col.cropgrowth, ns); + 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); @@ -102,26 +111,17 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings) end end - % output the vegetation dynamic results ydy - if options.calc_vegetation_dynamic - cropgrowth_names = {'DOY','DVS','LAI',... - 'PH','Sfactor','RootDM','LeafDM','StemDM','OrganDM','RootDeath','LeafDeath','StemDeath',}; - cropgrowth_units = {'day','-','m2/m2', ... - 'cm','-','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha'}; - write_output(cropgrowth_names, cropgrowth_units, fnames.cropgrowth_file, n_col.cropgrowth, ns); - end - % 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); @@ -194,7 +194,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) @@ -233,4 +232,4 @@ function write_output(header, units, bin_path, f_n_col, ns, not_header) fprintf(f_csv, '%d,', out_2d(k, 1:end - 1)); fprintf(f_csv, '%d\n', out_2d(k, end)); % saves from extra comma end -end +end \ No newline at end of file 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/+wofost/cropgrowth.asv b/src/+wofost/cropgrowth.asv deleted file mode 100644 index 6e828603..00000000 --- a/src/+wofost/cropgrowth.asv +++ /dev/null @@ -1,352 +0,0 @@ -function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,wofostpar,Anet,WaterStressFactor,xyt,KT,Dur_tot) -%{ - function cropgrowth.m simulate the growth of vegetation - - authors: Danyang Yu (yudanyang123@gmail.com) - date: 2025/01/01 - reference: - - Table of contents of the function - 1. Initilize crop growth paramters - 2. Calculate the phenological development - 3. Calculate the growth rate of different plant organs - 4. Output the crop growth variables - - Input: - Tsum Temperature sum from emergence to the simulated day - tdwi Initial total crop dry weight - laiem Leaf area index at emergence - f(r,l,s,o) The fraction of dry matter to root, leaf, stem, storage - tadw Initial total crop dry weight - lvage leaf age - lv save leaf dry matter at different steps - spa Specific pod area - ssa Specific stem area - sla specific leaf area - dmi converted dry matter at step KT - admi above-ground dry matter - grrt growth rate of root, similar for the other organs - drrt death rate of root, similar for the other organs - WofostPar Wofost parameters, the definitions are in CropD.crp - - Output: - dvs Developement of stage - lai leaf area index - ph Plant height - wrt dry matter weight of root - wlv dry matter weight of leaves - wst dry matter weight of stem - wso dry matter weight of storage organ - dwrt dry matter weight of dead root - dwlv dry matter weight of dead leaves - dwst dry matter weight of dead stem -%} - -%% 1. initilization -if KT == wofostpar.CSTART - % initial value of crop parameters - Tsum = 0; - dvs = 0; - tdwi = wofostpar.TDWI; - laiem = wofostpar.LAIEM; - ph = 0; - - frtb = wofostpar.FRTB; - fltb = wofostpar.FLTB; - fstb = wofostpar.FSTB; - fotb = wofostpar.FOTB; - - fr = wofost.afgen(frtb, dvs); - fl = wofost.afgen(fltb, dvs); - fs = wofost.afgen(fstb, dvs); - fo = wofost.afgen(fotb, dvs); - - ssa = wofostpar.SSA; - spa = wofostpar.SPA; - - % initial state variables of the crop - tadw = (1-fr)*tdwi; - wrt = fr*tdwi; - wlv = fl*tadw; - wst = fs*tadw; - wso = fo*tadw; - - dwrt = 0.0; - dwlv = 0.0; - dwst = 0.0; - - sla = zeros(Dur_tot+1,1); - lvage = zeros(Dur_tot+1,1); - lv = zeros(Dur_tot+1,1); - sla(1) = wofost.afgen(wofostpar.SLATB, dvs); - lv(1) = wlv; - lvage(1) = 0.0; - ilvold = 1; - - lasum = laiem; - laiexp = laiem; - glaiex = 0; - laimax = laiem; - lai = lasum+ssa*wst+spa*wso; - Anet_sum = 0; - lai_delta = 0; - -else - % Unpack state variables - dvs = state_vars.dvs; - wrt = state_vars.wrt; - wlv = state_vars.wlv; - wst = state_vars.wst; - wso = state_vars.wso; - - dwrt = state_vars.dwrt; - dwlv = state_vars.dwlv; - dwst = state_vars.dwst; - - sla = state_vars.sla; - lvage = state_vars.lvage; - lv = state_vars.lv; - lasum = state_vars.lasum; - laiexp = state_vars.laiexp; - lai = state_vars.lai; - - ph = state_vars.ph; - Anet_sum = state_vars.Anet_sum; - lai_delta = state_vars.lai_delta; -end - -%% 1. phenological development -delt = wofostpar.TSTEP/24; -tav = max(0, meteo.Ta); -dtsum = wofost.afgen(wofostpar.DTSMTB,tav); - -if dvs <= 1 % vegetative stage - dvs = dvs + dtsum*delt/wofostpar.TSUMEA; -else - dvs = dvs + dtsum*delt/wofostpar.TSUMAM; -end - -%% 2. dry matter increase -frtb = wofostpar.FRTB; -fltb = wofostpar.FLTB; -fstb = wofostpar.FSTB; -fotb = wofostpar.FOTB; - -fr = wofost.afgen(frtb, dvs); -fl = wofost.afgen(fltb, dvs); -fs = wofost.afgen(fstb, dvs); -fo = wofost.afgen(fotb, dvs); - -fcheck = fr+(fl+fs+fo)*(1.0-fr) - 1.0; %check on partitions -if abs(fcheck) > 0.0001 - print('The sum of partitioning factors for leaves, stems and storage organs is not equal to one at development stage'); - return; -end - -cvf = 1.0/((fl/wofostpar.CVL+fs/wofostpar.CVS+fo/wofostpar.CVO)*(1.0-fr)+fr/wofostpar.CVR); -asrc = Anet*30/1000000000*10000*3600*wofostpar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; -dmi = cvf * asrc; - -%% 3. Growth rate by root -% root extension -% rr = min(wofostpar.RDM - wofostpar.RD,wofostpar.RRI); -admi = (1.0-fr)*dmi; -grrt = fr*dmi; -drrt = wrt*wofost.afgen (wofostpar.RDRRTB,dvs)*delt; -gwrt = grrt-drrt; - -%% 4. Growth rate by stem -% growth rate stems -grst = fs*admi; -% death rate stems -drst = wofost.afgen (wofostpar.RDRSTB,dvs)*wst*delt; -% net growth rate stems -gwst = grst-drst; - -% growth of plant height -str = wofost.afgen(wofostpar.STRTB, dvs); -phnew = wst/(wofostpar.STN*wofostpar.STD*pi*str*str)*wofostpar.PHCoeff+wofostpar.PHEM; -if phnew >= ph - ph = phnew; -end - -%% 5. Growth rate by organs -gwso = fo*admi; - -%% 6. Growth rate by leave -% weight of new leaves -grlv = fl*admi; - -% death of leaves due to water stress or high lai -if abs(lai) < 0.1 - dslv1 = 0; -else - sfactor = WaterStressFactor.soil; - dslv1 = wlv*(1.0-sfactor)*wofostpar.PERDL*delt; -end - -laicr = wofostpar.LAICR; -dslv2 = wlv*max(0,0.03*(lai-laicr)/laicr); -dslv = max(dslv1,dslv2); - -% death of leaves due to exceeding life span -% first: leaf death due to water stress or high lai is imposed on array -% until no more leaves have to die or all leaves are gone -rest = dslv; -i1 = KT; -iday = ceil(KT*delt); - -while (rest>lv(i1) && i1 >= 1) - rest = rest - lv(i1); - i1 = i1 -1; -end - -% then: check if some of the remaining leaves are older than span, -% sum their weights -dalv = 0.0; -if (lvage(i1)>wofostpar.SPAN && rest>0 && i1>=1) - dalv = lv(i1) - rest; - rest = 0; - i1 = i1 -1; -end - -while (i1>=1 && lvage(i1)>wofostpar.SPAN) - dalv = dalv + lv(i1); - i1 = i1 - 1; -end - -% final death rate -drlv = dslv+dalv; - -% calculation of specific leaf area in case of exponential growth: -slat = wofost.afgen(wofostpar.SLATB, dvs); -if laiexp < 6 - dteff = max(0,tav-wofostpar.TBASE); % effective temperature - glaiex = laiexp*wofostpar.RGRLAI*dteff*delt; % exponential growth - laiexp = laiexp + glaiex; - - glasol = grlv*slat; % source-limited growth - gla = min(glaiex,glasol); -% gla = max(0,gla); - - slat = gla/grlv; - if isnan(slat) - slat = 0; - end -% if grlv>=0 -% slat = gla/grlv; % the actual sla value -% else -% slat = 0; -% end -end - -% update the information of leave -dslvt = dslv; -i1 = KT; -while (dslvt>0 && i1>=1) % water stress and high lai - if dslvt >= lv(i1) - dslvt = dslvt-lv(i1); - lv(i1) = 0.0; - i1 = i1-1; - else - lv(i1) = lv(i1)-dslvt; - dslvt = 0.0; - end -end - -while (lvage(i1)>=wofostpar.SPAN && i1>=1) % leaves older than span die - lv(i1) = 0; - i1 = i1-1; -end - -ilvold = KT; % oldest class with leaves -fysdel = max (0.0d0,(tav-wofostpar.TBASE)/(35.0-wofostpar.TBASE)); % physiologic ageing of leaves per time step - -for i1 = ilvold: -1: 1 % shifting of contents, updating of physiological age - lv(i1+1) = lv(i1); - sla(i1+1) = sla(i1); - lvage(i1+1) = lvage(i1)+fysdel*delt; -end -ilvold = ilvold+1; - -lv(1) = grlv; -sla(1) = slat; -lvage(1) = 0.0; - -% calculation of new leaf area and weight -lasum = 0.0; -wlv = 0.0; -for i1 = 1:ilvold - lasum = lasum + lv(i1)*sla(i1); - wlv = wlv + lv(i1); -end -lasum = max(0,lasum); -wlv = max(0,wlv); - -%% 7. integrals of the crop -% dry weight of living plant organs -wrt = wrt+gwrt; -wst = wst+gwst; -wso = wso+gwso; - -% total above ground biomass -tadw = wlv+wst+wso; - -% dry weight of dead plant organs (roots,leaves & stems) -dwrt = dwrt+drrt; -dwlv = dwlv+drlv; -dwst = dwst+drst; - -% dry weight of dead and living plant organs -twrt = wrt+dwrt; -twlv = wlv+dwlv; -twst = wst+dwst; -cwdm = twlv+twst+wso; - -% leaf area index -lai = wofostpar.LAIEM+lasum+wofostpar.SSA*wst+wofostpar.SPA*wso; - -Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration -if Anet_sum < 0 && KT>1 % the consumed biomass are the assimilated one at daytime - lai_delta = lai_delta + grlv*slat; - lai = lai - lai_delta; -elseif Anet_sum>=0 && Anet>=0 - Anet_sum = 0; - lai_delta = 0; -end - -%% 8. integrals of the crop -crop_output(KT,1) = xyt.t(KT,1); % Day of the year -crop_output(KT,2) = dvs; % Development of stage -crop_output(KT,3) = lai; % LAI -crop_output(KT,4) = ph; % Plant height -crop_output(KT,5) = sfactor; % Water stress -crop_output(KT,6) = wrt; % Dry matter of root -crop_output(KT,7) = wlv; % Dry matter of leaves -crop_output(KT,8) = wst; % Dry matter of stem -crop_output(KT,9) = wso; % Dry matter of organ -crop_output(KT,10) = dwrt; % Dry matter of deadleaves -crop_output(KT,11) = dwlv; % Dry matter of dead stem -crop_output(KT,12) = dwst; % Dry matter of dead organ - -%% 9. Pack updated state variables -state_vars.dvs = dvs; -state_vars.wrt = wrt; -state_vars.wlv = wlv; -state_vars.wst = wst; -state_vars.wso = wso; -state_vars.dwrt = dwrt; -state_vars.dwlv = dwlv; -state_vars.dwst = dwst; -state_vars.sla = sla; -state_vars.lvage = lvage; -state_vars.lv = lv; -state_vars.lasum = lasum; -state_vars.laiexp = laiexp; -state_vars.lai = lai; -state_vars.ph = ph; -state_vars.Anet_sum = Anet_sum; -state_vars.lai_delta = lai_delta; -end - - diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 2b44770c..1546c900 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); @@ -854,7 +854,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 From 3d4803786e030a5718dc6c761e1986f7b25d95c1 Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Tue, 7 Jan 2025 16:47:10 +0100 Subject: [PATCH 11/14] Revise according to MISS_HIT style --- src/+init/setBoundaryCondition.m | 4 +- src/+io/bin_to_csv.m | 10 +- src/+io/create_output_files_binary.m | 4 +- src/+io/loadSoilInitialValues.m | 4 +- src/+io/output_data_binary.m | 6 +- src/+wofost/WofostRead.m | 90 ++-- src/+wofost/afgen.m | 28 +- src/+wofost/cropgrowth.asv | 353 ++++++++++++++ src/+wofost/cropgrowth.m | 664 +++++++++++++-------------- src/STEMMUS_SCOPE.m | 26 +- 10 files changed, 771 insertions(+), 418 deletions(-) create mode 100644 src/+wofost/cropgrowth.asv diff --git a/src/+init/setBoundaryCondition.m b/src/+init/setBoundaryCondition.m index fd5a8b2d..b2520967 100644 --- a/src/+init/setBoundaryCondition.m +++ b/src/+init/setBoundaryCondition.m @@ -82,10 +82,10 @@ % 3 --Zero temperature gradient; NBCTB = 1; - if mean(Ta_msr,'omitnan') < 0 + if mean(Ta_msr, 'omitnan') < 0 BCTB = 0; % 9 8.1 else - BCTB = mean(Ta_msr,'omitnan'); + BCTB = mean(Ta_msr, 'omitnan'); end end if ModelSettings.Soilairefc == 1 diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m index 51613bea..0404a4b1 100644 --- a/src/+io/bin_to_csv.m +++ b/src/+io/bin_to_csv.m @@ -92,11 +92,11 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings, write_output(Sim_qtot_names, Sim_qtot_units, fnames.Sim_qtot_file, n_col.Sim_qtot, ns, true); %% output the vegetation dynamic results Danyang Yu - if options.calc_vegetation_dynamic - cropgrowth_names = {'DOY','DVS','LAI',... - 'PH','Sfactor','RootDM','LeafDM','StemDM','OrganDM','RootDeath','LeafDeath','StemDeath',}; - cropgrowth_units = {'day','-','m2/m2', ... - 'cm','-','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha','kg/ha'}; + if options.calc_vegetation_dynamic + cropgrowth_names = {'DOY', 'DVS', 'LAI', ... + 'PH', 'Sfactor', 'RootDM', 'LeafDM', 'StemDM', 'OrganDM', 'RootDeath', 'LeafDeath', 'StemDeath',}; + cropgrowth_units = {'day', '-', 'm2/m2', ... + 'cm', '-', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha'}; write_output(cropgrowth_names, cropgrowth_units, fnames.cropgrowth_file, n_col.cropgrowth, ns); end diff --git a/src/+io/create_output_files_binary.m b/src/+io/create_output_files_binary.m index 303e5082..0e5e2dbc 100644 --- a/src/+io/create_output_files_binary.m +++ b/src/+io/create_output_files_binary.m @@ -107,10 +107,10 @@ fnames.spectrum_obsdir_thermal_file = fullfile(Output_dir, 'spectrum_obsdir_thermal.bin'); % spectrum observation direction fnames.spectrum_hemis_thermal_file = fullfile(Output_dir, 'spectrum_hemis_thermal.bin'); % spectrum hemispherically integrated end - + % Create cropgrowth file if options.calc_vegetation_dynamic - fnames.cropgrowth_file = fullfile(Output_dir,'cropgrowth.bin'); % crop growth simulation ydy + fnames.cropgrowth_file = fullfile(Output_dir, 'cropgrowth.bin'); % crop growth simulation ydy end % Create empty files for appending diff --git a/src/+io/loadSoilInitialValues.m b/src/+io/loadSoilInitialValues.m index ae5427a4..63285173 100644 --- a/src/+io/loadSoilInitialValues.m +++ b/src/+io/loadSoilInitialValues.m @@ -68,10 +68,10 @@ InitT6 = 0; Tss = InitT0; end - if mean(Ta_msr,'omitnan') < 0 + if mean(Ta_msr, 'omitnan') < 0 BtmT = 0; % 9 8.1 else - BtmT = mean(Ta_msr,'omitnan'); + BtmT = mean(Ta_msr, 'omitnan'); end if InitX0 > SaturatedMC(1) || InitX1 > SaturatedMC(1) || InitX2 > SaturatedMC(2) || ... InitX3 > SaturatedMC(3) || InitX4 > SaturatedMC(4) || InitX5 > SaturatedMC(5) || InitX6 > SaturatedMC(6) diff --git a/src/+io/output_data_binary.m b/src/+io/output_data_binary.m index d84b206d..f5db2040 100644 --- a/src/+io/output_data_binary.m +++ b/src/+io/output_data_binary.m @@ -1,7 +1,7 @@ function n_col = output_data_binary(f, k, xyt, rad, canopy, ScopeParameters, vi, vmax, options, fluxes, meteo, iter, thermal, ... spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap, WaterStress, WaterPotential, ... Sim_hh, Sim_qlh, Sim_qlt, Sim_qvh, Sim_qvt, Sim_qla, Sim_qva, Sim_qtot, ... - ForcingData, RS, RWUs, RWUg,crop_output) + ForcingData, RS, RWUs, RWUg, crop_output) %% OUTPUT DATA % author C. Van der Tol @@ -41,9 +41,9 @@ %% Crop growth if options.calc_vegetation_dynamic == 1 - cropgrowth_out = crop_output(k,:); + cropgrowth_out = crop_output(k, :); n_col.cropgrowth = length(cropgrowth_out); - fwrite(f.cropgrowth_file,cropgrowth_out,'double'); + fwrite(f.cropgrowth_file, cropgrowth_out, 'double'); end %% Water stress factor diff --git a/src/+wofost/WofostRead.m b/src/+wofost/WofostRead.m index c6d19034..1166279d 100644 --- a/src/+wofost/WofostRead.m +++ b/src/+wofost/WofostRead.m @@ -1,55 +1,55 @@ function [wofost] = WofostRead(path_input) -%{ - function WofostRead.m loads plant growth parameters - from file "CropD.crp" in the input folder. + %{ + function WofostRead.m loads plant growth parameters + from file "CropD.crp" in the input folder. - authors: Danyang Yu (yudanyang123@gmail.com) - date: 2025/01/01 -%} + authors: Danyang Yu (yudanyang123@gmail.com) + date: 2025/01/01 + %} -% Read lines from the file -PlantGrowthFile = [path_input,'plant_growth/CropD.crp']; -fid = fopen(PlantGrowthFile); - -% process it line by line -while true - tline = fgetl(fid); + % Read lines from the file + PlantGrowthFile = [path_input, 'plant_growth/CropD.crp']; + fid = fopen(PlantGrowthFile); - % remove empty and note line - if ~isempty(tline) & tline(1) ~= '*' - % judge whether the string contain the table value - if contains(tline,'=') - s = regexp(tline,'=','split'); - vname = strtrim(char(s(1))); - % assign value - if ~isempty(char(s(2))) - values = regexp(char(s(2)),'!','split'); - value = str2num(strtrim(char(values(1)))); - wofost.(vname) = value; % save parameters in wofost struct - end - else - % save tabel value - i = 1; - table_value = []; - while ~isempty(tline) & tline(1) ~= '*' - s = strsplit(strtrim(tline)); - table_value(i,1) = str2num(strtrim(char(s(1)))); - table_value(i,2) = str2num(strtrim(char(s(2)))); - i = i + 1; - tline = fgetl(fid); -% disp(tline); + % process it line by line + while true + tline = fgetl(fid); + + % remove empty and note line + if ~isempty(tline) & tline(1) ~= '*' + % judge whether the string contain the table value + if contains(tline, '=') + s = regexp(tline, '=', 'split'); + vname = strtrim(char(s(1))); + % assign value + if ~isempty(char(s(2))) + values = regexp(char(s(2)), '!', 'split'); + value = str2num(strtrim(char(values(1)))); + wofost.(vname) = value; % save parameters in wofost struct + end + else + % save tabel value + i = 1; + table_value = []; + while ~isempty(tline) & tline(1) ~= '*' + s = strsplit(strtrim(tline)); + table_value(i,1) = str2num(strtrim(char(s(1)))); + table_value(i,2) = str2num(strtrim(char(s(2)))); + i = i + 1; + tline = fgetl(fid); + % disp(tline); + end + + wofost.(vname) = table_value; end - - wofost.(vname) = table_value; end + + %end of file + if contains(tline,'* End of .crp file !') + break; + end + % disp(tline); end - - %end of file - if contains(tline,'* End of .crp file !') - break; - end -% disp(tline); -end end diff --git a/src/+wofost/afgen.m b/src/+wofost/afgen.m index 7aac03a9..5e96ec83 100644 --- a/src/+wofost/afgen.m +++ b/src/+wofost/afgen.m @@ -1,17 +1,17 @@ function [output] = afgen(table,x) -%{ - function afgen.m intercept the value in the table - authors: Danyang Yu (yudanyang123@gmail.com) - date: 2025/01/01 -%} - -if x < table(1,1) | x > table(end,1) - print('the value of developemnt stage is mistaken'); -else - dvslist = table(:,1).'; - loc = discretize([x],[-Inf dvslist Inf]); - slope = (table(loc,2) - table((loc-1),2))/(table(loc,1)-table((loc-1),1)); - output = table(loc-1,2) + slope*(x - table(loc-1,1)); -end + %{ + function afgen.m intercept the value in the table + authors: Danyang Yu (yudanyang123@gmail.com) + date: 2025/01/01 + %} + + if x < table(1, 1) || x > table(end, 1) + print('the value of developemnt stage is mistaken'); + else + dvslist = table(:, 1).'; + loc = discretize([x], [-Inf dvslist Inf]); + slope = (table(loc, 2) - table((loc-1), 2)) / (table(loc, 1) - table((loc - 1), 1)); + output = table(loc - 1, 2) + slope * (x - table(loc - 1, 1)); + end end diff --git a/src/+wofost/cropgrowth.asv b/src/+wofost/cropgrowth.asv new file mode 100644 index 00000000..271874a1 --- /dev/null +++ b/src/+wofost/cropgrowth.asv @@ -0,0 +1,353 @@ +function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,WofostPar,Anet,WaterStressFactor,xyt,KT,Dur_tot) + %{ + function cropgrowth.m simulate the growth of vegetation + + authors: Danyang Yu (yudanyang123@gmail.com) + date: 2025/01/01 + reference: SWAP version 3.2. Theory description and user manual (Page 147-163) + + Table of contents of the function + 1. Initilize crop growth paramters + 2. Calculate the phenological development + 3. Calculate the growth rate of different plant organs + 4. Output the crop growth variables + + Input: + Tsum Temperature sum from emergence to the simulated day + tdwi Initial total crop dry weight + laiem Leaf area index at emergence + f(r,l,s,o) The fraction of dry matter to root, leaf, stem, storage + tadw Initial total crop dry weight + lvage leaf age + lv save leaf dry matter at different steps + spa Specific pod area + ssa Specific stem area + sla specific leaf area + Anet_sum the sum of net assimilation + dmi converted dry matter at step KT + admi above-ground dry matter + grrt growth rate of root, similar for the other organs + drrt death rate of root, similar for the other organs + WofostPar Wofost parameters, the definitions are in CropD.crp + + Output: + dvs Developement of stage + lai leaf area index + ph Plant height + wrt dry matter weight of root + wlv dry matter weight of leaves + wst dry matter weight of stem + wso dry matter weight of storage organ + dwrt dry matter weight of dead root + dwlv dry matter weight of dead leaves + dwst dry matter weight of dead stem + %} + + %% 1. initilization + if KT == WofostPar.CSTART + % initial value of crop parameters + Tsum = 0; + dvs = 0; + tdwi = WofostPar.TDWI; + laiem = WofostPar.LAIEM; + ph = 0; + + frtb = WofostPar.FRTB; + fltb = WofostPar.FLTB; + fstb = WofostPar.FSTB; + fotb = WofostPar.FOTB; + + fr = wofost.afgen(frtb, dvs); + fl = wofost.afgen(fltb, dvs); + fs = wofost.afgen(fstb, dvs); + fo = wofost.afgen(fotb, dvs); + + ssa = WofostPar.SSA; + spa = WofostPar.SPA; + + % initial state variables of the crop + tadw = (1 - fr) * tdwi; + wrt = fr * tdwi; + wlv = fl * tadw; + wst = fs * tadw; + wso = fo * tadw; + + dwrt = 0.0; + dwlv = 0.0; + dwst = 0.0; + + sla = zeros(Dur_tot+1, 1); + lvage = zeros(Dur_tot+1, 1); + lv = zeros(Dur_tot+1, 1); + sla(1) = wofost.afgen(WofostPar.SLATB, dvs); + lv(1) = wlv; + lvage(1) = 0.0; + ilvold = 1; + + lasum = laiem; + laiexp = laiem; + glaiex = 0; + laimax = laiem; + lai = lasum + ssa * wst + spa * wso; + Anet_sum = 0; + lai_delta = 0; + + else + % Unpack state variables + dvs = state_vars.dvs; + wrt = state_vars.wrt; + wlv = state_vars.wlv; + wst = state_vars.wst; + wso = state_vars.wso; + + dwrt = state_vars.dwrt; + dwlv = state_vars.dwlv; + dwst = state_vars.dwst; + + sla = state_vars.sla; + lvage = state_vars.lvage; + lv = state_vars.lv; + lasum = state_vars.lasum; + laiexp = state_vars.laiexp; + lai = state_vars.lai; + + ph = state_vars.ph; + Anet_sum = state_vars.Anet_sum; + lai_delta = state_vars.lai_delta; + end + + %% 1. phenological development + delt = WofostPar.TSTEP / 24; + tav = max(0, meteo.Ta); + dtsum = wofost.afgen(WofostPar.DTSMTB, tav); + + if dvs <= 1 % vegetative stage + dvs = dvs + dtsum * delt / WofostPar.TSUMEA; + else + dvs = dvs + dtsum * delt / WofostPar.TSUMAM; + end + + %% 2. dry matter increase + frtb = WofostPar.FRTB; + fltb = WofostPar.FLTB; + fstb = WofostPar.FSTB; + fotb = WofostPar.FOTB; + + fr = wofost.afgen(frtb, dvs); + fl = wofost.afgen(fltb, dvs); + fs = wofost.afgen(fstb, dvs); + fo = wofost.afgen(fotb, dvs); + + fcheck = fr + (fl + fs + fo) * (1.0 - fr) - 1.0; %check on partitions + if abs(fcheck) > 0.0001 + print('The sum of partitioning factors for leaves, stems and storage organs is not equal to one at development stage'); + return; + end + + cvf = 1.0 / ((fl / WofostPar.CVL + fs / WofostPar.CVS + fo / WofostPar.CVO) * (1.0 -fr ) + fr / WofostPar.CVR); + asrc = Anet * 30 / 1000000000 * 10000 * 3600 * WofostPar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; + dmi = cvf * asrc; + + %% 3. Growth rate by root + % root extension + % rr = min(WofostPar.RDM - WofostPar.RD,WofostPar.RRI); + admi = (1.0 - fr) * dmi; + grrt = fr * dmi; + drrt = wrt * wofost.afgen(WofostPar.RDRRTB, dvs) * delt; + gwrt = grrt - drrt; + + %% 4. Growth rate by stem + % growth rate stems + grst = fs * admi; + % death rate stems + drst = wofost.afgen(WofostPar.RDRSTB, dvs) * wst * delt; + % net growth rate stems + gwst = grst - drst; + + % growth of plant height + str = wofost.afgen(WofostPar.STRTB, dvs); + phnew = wst / (WofostPar.STN * WofostPar.STD * pi * str * str) * WofostPar.PHCoeff + WofostPar.PHEM; + if phnew >= ph + ph = phnew; + end + + %% 5. Growth rate by organs + gwso = fo * admi; + + %% 6. Growth rate by leave + % weight of new leaves + grlv = fl * admi; + + % death of leaves due to water stress or high lai + if abs(lai) < 0.1 + dslv1 = 0; + else + sfactor = WaterStressFactor.soil; + dslv1 = wlv * (1.0 - sfactor) * WofostPar.PERDL * delt; + end + + laicr = WofostPar.LAICR; + dslv2 = wlv * max(0, 0.03*(lai - laicr) / laicr); + dslv = max(dslv1, dslv2); + + % death of leaves due to exceeding life span + % first: leaf death due to water stress or high lai is imposed on array + % until no more leaves have to die or all leaves are gone + rest = dslv; + i1 = KT; + iday = ceil(KT*delt); + + while (rest > lv(i1) && i1 >= 1) + rest = rest - lv(i1); + i1 = i1 - 1; + end + + % then: check if some of the remaining leaves are older than span, + % sum their weights + dalv = 0.0; + if (lvage(i1) > WofostPar.SPAN && rest > 0 && i1 >= 1) + dalv = lv(i1) - rest; + rest = 0; + i1 = i1 - 1; + end + + while (i1 >= 1 && lvage(i1) > WofostPar.SPAN) + dalv = dalv + lv(i1); + i1 = i1 - 1; + end + + % final death rate + drlv = dslv + dalv; + + % calculation of specific leaf area in case of exponential growth: + slat = wofost.afgen(WofostPar.SLATB, dvs); + if laiexp < 6 + dteff = max(0, tav - WofostPar.TBASE); % effective temperature + glaiex = laiexp * WofostPar.RGRLAI * dteff * delt; % exponential growth + laiexp = laiexp + glaiex; + + glasol = grlv * slat; % source-limited growth + gla = min(glaiex, glasol); + % gla = max(0,gla); + + slat = gla / grlv; + if isnan(slat) + slat = 0; + end + % if grlv>=0 + % slat = gla/grlv; % the actual sla value + % else + % slat = 0; + % end + end + + % update the information of leave + dslvt = dslv; + i1 = KT; + while (dslvt>0 && i1>=1) % water stress and high lai + if dslvt >= lv(i1) + dslvt = dslvt-lv(i1); + lv(i1) = 0.0; + i1 = i1-1; + else + lv(i1) = lv(i1)-dslvt; + dslvt = 0.0; + end + end + + while (lvage(i1) >= WofostPar.SPAN && i1 >= 1) % leaves older than span die + lv(i1) = 0; + i1 = i1 - 1; + end + + ilvold = KT; % oldest class with leaves + fysdel = max(0.0d0, (tav - WofostPar.TBASE) / (35.0 - WofostPar.TBASE)); % physiologic ageing of leaves per time step + + for i1 = ilvold: -1: 1 % shifting of contents, updating of physiological age + lv(i1+1) = lv(i1); + sla(i1+1) = sla(i1); + lvage(i1+1) = lvage(i1) + fysdel * delt; + end + ilvold = ilvold + 1; + + lv(1) = grlv; + sla(1) = slat; + lvage(1) = 0.0; + + % calculation of new leaf area and weight + lasum = 0.0; + wlv = 0.0; + for i1 = 1:ilvold + lasum = lasum + lv(i1)*sla(i1); + wlv = wlv + lv(i1); + end + lasum = max(0,lasum); + wlv = max(0,wlv); + + %% 7. integrals of the crop + % dry weight of living plant organs + wrt = wrt+gwrt; + wst = wst+gwst; + wso = wso+gwso; + + % total above ground biomass + tadw = wlv+wst+wso; + + % dry weight of dead plant organs (roots,leaves & stems) + dwrt = dwrt+drrt; + dwlv = dwlv+drlv; + dwst = dwst+drst; + + % dry weight of dead and living plant organs + twrt = wrt+dwrt; + twlv = wlv+dwlv; + twst = wst+dwst; + cwdm = twlv+twst+wso; + + % leaf area index + lai = WofostPar.LAIEM+lasum+WofostPar.SSA*wst+WofostPar.SPA*wso; + + Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration + if Anet_sum < 0 && KT>1 % the consumed biomass are the assimilated one at daytime + lai_delta = lai_delta + grlv*slat; + lai = lai - lai_delta; + elseif Anet_sum>=0 && Anet>=0 + Anet_sum = 0; + lai_delta = 0; + end + + %% 8. integrals of the crop + crop_output(KT,1) = xyt.t(KT,1); % Day of the year + crop_output(KT,2) = dvs; % Development of stage + crop_output(KT,3) = lai; % LAI + crop_output(KT,4) = ph; % Plant height + crop_output(KT,5) = sfactor; % Water stress + crop_output(KT,6) = wrt; % Dry matter of root + crop_output(KT,7) = wlv; % Dry matter of leaves + crop_output(KT,8) = wst; % Dry matter of stem + crop_output(KT,9) = wso; % Dry matter of organ + crop_output(KT,10) = dwrt; % Dry matter of deadleaves + crop_output(KT,11) = dwlv; % Dry matter of dead stem + crop_output(KT,12) = dwst; % Dry matter of dead organ + + %% 9. Pack updated state variables + state_vars.dvs = dvs; + state_vars.wrt = wrt; + state_vars.wlv = wlv; + state_vars.wst = wst; + state_vars.wso = wso; + state_vars.dwrt = dwrt; + state_vars.dwlv = dwlv; + state_vars.dwst = dwst; + state_vars.sla = sla; + state_vars.lvage = lvage; + state_vars.lv = lv; + state_vars.lasum = lasum; + state_vars.laiexp = laiexp; + state_vars.lai = lai; + state_vars.ph = ph; + state_vars.Anet_sum = Anet_sum; + state_vars.lai_delta = lai_delta; +end + + diff --git a/src/+wofost/cropgrowth.m b/src/+wofost/cropgrowth.m index 4d1dcf67..e1d4b009 100644 --- a/src/+wofost/cropgrowth.m +++ b/src/+wofost/cropgrowth.m @@ -1,57 +1,133 @@ function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,WofostPar,Anet,WaterStressFactor,xyt,KT,Dur_tot) -%{ - function cropgrowth.m simulate the growth of vegetation - - authors: Danyang Yu (yudanyang123@gmail.com) - date: 2025/01/01 - reference: SWAP version 3.2. Theory description and user manual (Page 147-163) - - Table of contents of the function - 1. Initilize crop growth paramters - 2. Calculate the phenological development - 3. Calculate the growth rate of different plant organs - 4. Output the crop growth variables - - Input: - Tsum Temperature sum from emergence to the simulated day - tdwi Initial total crop dry weight - laiem Leaf area index at emergence - f(r,l,s,o) The fraction of dry matter to root, leaf, stem, storage - tadw Initial total crop dry weight - lvage leaf age - lv save leaf dry matter at different steps - spa Specific pod area - ssa Specific stem area - sla specific leaf area - Anet_sum the sum of net assimilation - dmi converted dry matter at step KT - admi above-ground dry matter - grrt growth rate of root, similar for the other organs - drrt death rate of root, similar for the other organs - WofostPar Wofost parameters, the definitions are in CropD.crp - - Output: - dvs Developement of stage - lai leaf area index - ph Plant height - wrt dry matter weight of root - wlv dry matter weight of leaves - wst dry matter weight of stem - wso dry matter weight of storage organ - dwrt dry matter weight of dead root - dwlv dry matter weight of dead leaves - dwst dry matter weight of dead stem -%} - -%% 1. initilization -if KT == WofostPar.CSTART - % initial value of crop parameters - Tsum = 0; - dvs = 0; - tdwi = WofostPar.TDWI; - laiem = WofostPar.LAIEM; - ph = 0; + %{ + function cropgrowth.m simulate the growth of vegetation + + authors: Danyang Yu (yudanyang123@gmail.com) + date: 2025/01/01 + reference: SWAP version 3.2. Theory description and user manual (Page 147-163) + + Table of contents of the function + 1. Initilize crop growth paramters + 2. Calculate the phenological development + 3. Calculate the growth rate of different plant organs + 4. Output the crop growth variables + + Input: + Tsum Temperature sum from emergence to the simulated day + tdwi Initial total crop dry weight + laiem Leaf area index at emergence + f(r,l,s,o) The fraction of dry matter to root, leaf, stem, storage + tadw Initial total crop dry weight + lvage leaf age + lv save leaf dry matter at different steps + spa Specific pod area + ssa Specific stem area + sla specific leaf area + Anet_sum the sum of net assimilation + dmi converted dry matter at step KT + admi above-ground dry matter + grrt growth rate of root, similar for the other organs + drrt death rate of root, similar for the other organs + WofostPar Wofost parameters, the definitions are in CropD.crp + + Output: + dvs Developement of stage + lai leaf area index + ph Plant height + wrt dry matter weight of root + wlv dry matter weight of leaves + wst dry matter weight of stem + wso dry matter weight of storage organ + dwrt dry matter weight of dead root + dwlv dry matter weight of dead leaves + dwst dry matter weight of dead stem + %} + + %% 1. initilization + if KT == WofostPar.CSTART + % initial value of crop parameters + Tsum = 0; + dvs = 0; + tdwi = WofostPar.TDWI; + laiem = WofostPar.LAIEM; + ph = 0; + + frtb = WofostPar.FRTB; + fltb = WofostPar.FLTB; + fstb = WofostPar.FSTB; + fotb = WofostPar.FOTB; + + fr = wofost.afgen(frtb, dvs); + fl = wofost.afgen(fltb, dvs); + fs = wofost.afgen(fstb, dvs); + fo = wofost.afgen(fotb, dvs); + + ssa = WofostPar.SSA; + spa = WofostPar.SPA; + + % initial state variables of the crop + tadw = (1 - fr) * tdwi; + wrt = fr * tdwi; + wlv = fl * tadw; + wst = fs * tadw; + wso = fo * tadw; + + dwrt = 0.0; + dwlv = 0.0; + dwst = 0.0; + + sla = zeros(Dur_tot+1, 1); + lvage = zeros(Dur_tot+1, 1); + lv = zeros(Dur_tot+1, 1); + sla(1) = wofost.afgen(WofostPar.SLATB, dvs); + lv(1) = wlv; + lvage(1) = 0.0; + ilvold = 1; + + lasum = laiem; + laiexp = laiem; + glaiex = 0; + laimax = laiem; + lai = lasum + ssa * wst + spa * wso; + Anet_sum = 0; + lai_delta = 0; + + else + % Unpack state variables + dvs = state_vars.dvs; + wrt = state_vars.wrt; + wlv = state_vars.wlv; + wst = state_vars.wst; + wso = state_vars.wso; + dwrt = state_vars.dwrt; + dwlv = state_vars.dwlv; + dwst = state_vars.dwst; + + sla = state_vars.sla; + lvage = state_vars.lvage; + lv = state_vars.lv; + lasum = state_vars.lasum; + laiexp = state_vars.laiexp; + lai = state_vars.lai; + + ph = state_vars.ph; + Anet_sum = state_vars.Anet_sum; + lai_delta = state_vars.lai_delta; + end + + %% 1. phenological development + delt = WofostPar.TSTEP / 24; + tav = max(0, meteo.Ta); + dtsum = wofost.afgen(WofostPar.DTSMTB, tav); + + if dvs <= 1 % vegetative stage + dvs = dvs + dtsum * delt / WofostPar.TSUMEA; + else + dvs = dvs + dtsum * delt / WofostPar.TSUMAM; + end + + %% 2. dry matter increase frtb = WofostPar.FRTB; fltb = WofostPar.FLTB; fstb = WofostPar.FSTB; @@ -62,292 +138,216 @@ fs = wofost.afgen(fstb, dvs); fo = wofost.afgen(fotb, dvs); - ssa = WofostPar.SSA; - spa = WofostPar.SPA; + fcheck = fr + (fl + fs + fo) * (1.0 - fr) - 1.0; %check on partitions + if abs(fcheck) > 0.0001 + print('The sum of partitioning factors for leaves, stems and storage organs is not equal to one at development stage'); + return; + end + + cvf = 1.0 / ((fl / WofostPar.CVL + fs / WofostPar.CVS + fo / WofostPar.CVO) * (1.0 -fr ) + fr / WofostPar.CVR); + asrc = Anet * 30 / 1000000000 * 10000 * 3600 * WofostPar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; + dmi = cvf * asrc; + + %% 3. Growth rate by root + % root extension + % rr = min(WofostPar.RDM - WofostPar.RD,WofostPar.RRI); + admi = (1.0 - fr) * dmi; + grrt = fr * dmi; + drrt = wrt * wofost.afgen(WofostPar.RDRRTB, dvs) * delt; + gwrt = grrt - drrt; + + %% 4. Growth rate by stem + % growth rate stems + grst = fs * admi; + % death rate stems + drst = wofost.afgen(WofostPar.RDRSTB, dvs) * wst * delt; + % net growth rate stems + gwst = grst - drst; + + % growth of plant height + str = wofost.afgen(WofostPar.STRTB, dvs); + phnew = wst / (WofostPar.STN * WofostPar.STD * pi * str * str) * WofostPar.PHCoeff + WofostPar.PHEM; + if phnew >= ph + ph = phnew; + end + + %% 5. Growth rate by organs + gwso = fo * admi; + + %% 6. Growth rate by leave + % weight of new leaves + grlv = fl * admi; + + % death of leaves due to water stress or high lai + if abs(lai) < 0.1 + dslv1 = 0; + else + sfactor = WaterStressFactor.soil; + dslv1 = wlv * (1.0 - sfactor) * WofostPar.PERDL * delt; + end + + laicr = WofostPar.LAICR; + dslv2 = wlv * max(0, 0.03*(lai - laicr) / laicr); + dslv = max(dslv1, dslv2); + + % death of leaves due to exceeding life span + % first: leaf death due to water stress or high lai is imposed on array + % until no more leaves have to die or all leaves are gone + rest = dslv; + i1 = KT; + iday = ceil(KT*delt); - % initial state variables of the crop - tadw = (1-fr)*tdwi; - wrt = fr*tdwi; - wlv = fl*tadw; - wst = fs*tadw; - wso = fo*tadw; + while (rest > lv(i1) && i1 >= 1) + rest = rest - lv(i1); + i1 = i1 - 1; + end + + % then: check if some of the remaining leaves are older than span, + % sum their weights + dalv = 0.0; + if (lvage(i1) > WofostPar.SPAN && rest > 0 && i1 >= 1) + dalv = lv(i1) - rest; + rest = 0; + i1 = i1 - 1; + end - dwrt = 0.0; - dwlv = 0.0; - dwst = 0.0; + while (i1 >= 1 && lvage(i1) > WofostPar.SPAN) + dalv = dalv + lv(i1); + i1 = i1 - 1; + end + + % final death rate + drlv = dslv + dalv; + + % calculation of specific leaf area in case of exponential growth: + slat = wofost.afgen(WofostPar.SLATB, dvs); + if laiexp < 6 + dteff = max(0, tav - WofostPar.TBASE); % effective temperature + glaiex = laiexp * WofostPar.RGRLAI * dteff * delt; % exponential growth + laiexp = laiexp + glaiex; + + glasol = grlv * slat; % source-limited growth + gla = min(glaiex, glasol); + % gla = max(0,gla); - sla = zeros(Dur_tot+1,1); - lvage = zeros(Dur_tot+1,1); - lv = zeros(Dur_tot+1,1); - sla(1) = wofost.afgen(WofostPar.SLATB, dvs); - lv(1) = wlv; + slat = gla / grlv; + if isnan(slat) + slat = 0; + end + % if grlv>=0 + % slat = gla/grlv; % the actual sla value + % else + % slat = 0; + % end + end + + % update the information of leave + dslvt = dslv; + i1 = KT; + while (dslvt>0 && i1>=1) % water stress and high lai + if dslvt >= lv(i1) + dslvt = dslvt-lv(i1); + lv(i1) = 0.0; + i1 = i1-1; + else + lv(i1) = lv(i1)-dslvt; + dslvt = 0.0; + end + end + + while (lvage(i1) >= WofostPar.SPAN && i1 >= 1) % leaves older than span die + lv(i1) = 0; + i1 = i1 - 1; + end + + ilvold = KT; % oldest class with leaves + fysdel = max(0.0d0, (tav - WofostPar.TBASE) / (35.0 - WofostPar.TBASE)); % physiologic ageing of leaves per time step + + for i1 = ilvold: -1: 1 % shifting of contents, updating of physiological age + lv(i1+1) = lv(i1); + sla(i1+1) = sla(i1); + lvage(i1+1) = lvage(i1) + fysdel * delt; + end + ilvold = ilvold + 1; + + lv(1) = grlv; + sla(1) = slat; lvage(1) = 0.0; - ilvold = 1; - - lasum = laiem; - laiexp = laiem; - glaiex = 0; - laimax = laiem; - lai = lasum+ssa*wst+spa*wso; - Anet_sum = 0; - lai_delta = 0; - -else - % Unpack state variables - dvs = state_vars.dvs; - wrt = state_vars.wrt; - wlv = state_vars.wlv; - wst = state_vars.wst; - wso = state_vars.wso; - - dwrt = state_vars.dwrt; - dwlv = state_vars.dwlv; - dwst = state_vars.dwst; - - sla = state_vars.sla; - lvage = state_vars.lvage; - lv = state_vars.lv; - lasum = state_vars.lasum; - laiexp = state_vars.laiexp; - lai = state_vars.lai; - - ph = state_vars.ph; - Anet_sum = state_vars.Anet_sum; - lai_delta = state_vars.lai_delta; -end - -%% 1. phenological development -delt = WofostPar.TSTEP/24; -tav = max(0, meteo.Ta); -dtsum = wofost.afgen(WofostPar.DTSMTB,tav); - -if dvs <= 1 % vegetative stage - dvs = dvs + dtsum*delt/WofostPar.TSUMEA; -else - dvs = dvs + dtsum*delt/WofostPar.TSUMAM; -end - -%% 2. dry matter increase -frtb = WofostPar.FRTB; -fltb = WofostPar.FLTB; -fstb = WofostPar.FSTB; -fotb = WofostPar.FOTB; - -fr = wofost.afgen(frtb, dvs); -fl = wofost.afgen(fltb, dvs); -fs = wofost.afgen(fstb, dvs); -fo = wofost.afgen(fotb, dvs); - -fcheck = fr+(fl+fs+fo)*(1.0-fr) - 1.0; %check on partitions -if abs(fcheck) > 0.0001 - print('The sum of partitioning factors for leaves, stems and storage organs is not equal to one at development stage'); - return; -end - -cvf = 1.0/((fl/WofostPar.CVL+fs/WofostPar.CVS+fo/WofostPar.CVO)*(1.0-fr)+fr/WofostPar.CVR); -asrc = Anet*30/1000000000*10000*3600*WofostPar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; -dmi = cvf * asrc; - -%% 3. Growth rate by root -% root extension -% rr = min(WofostPar.RDM - WofostPar.RD,WofostPar.RRI); -admi = (1.0-fr)*dmi; -grrt = fr*dmi; -drrt = wrt*wofost.afgen (WofostPar.RDRRTB,dvs)*delt; -gwrt = grrt-drrt; - -%% 4. Growth rate by stem -% growth rate stems -grst = fs*admi; -% death rate stems -drst = wofost.afgen (WofostPar.RDRSTB,dvs)*wst*delt; -% net growth rate stems -gwst = grst-drst; - -% growth of plant height -str = wofost.afgen(WofostPar.STRTB, dvs); -phnew = wst/(WofostPar.STN*WofostPar.STD*pi*str*str)*WofostPar.PHCoeff+WofostPar.PHEM; -if phnew >= ph - ph = phnew; -end - -%% 5. Growth rate by organs -gwso = fo*admi; - -%% 6. Growth rate by leave -% weight of new leaves -grlv = fl*admi; - -% death of leaves due to water stress or high lai -if abs(lai) < 0.1 - dslv1 = 0; -else - sfactor = WaterStressFactor.soil; - dslv1 = wlv*(1.0-sfactor)*WofostPar.PERDL*delt; -end - -laicr = WofostPar.LAICR; -dslv2 = wlv*max(0,0.03*(lai-laicr)/laicr); -dslv = max(dslv1,dslv2); - -% death of leaves due to exceeding life span -% first: leaf death due to water stress or high lai is imposed on array -% until no more leaves have to die or all leaves are gone -rest = dslv; -i1 = KT; -iday = ceil(KT*delt); - -while (rest>lv(i1) && i1 >= 1) - rest = rest - lv(i1); - i1 = i1 -1; -end - -% then: check if some of the remaining leaves are older than span, -% sum their weights -dalv = 0.0; -if (lvage(i1)>WofostPar.SPAN && rest>0 && i1>=1) - dalv = lv(i1) - rest; - rest = 0; - i1 = i1 -1; -end - -while (i1>=1 && lvage(i1)>WofostPar.SPAN) - dalv = dalv + lv(i1); - i1 = i1 - 1; -end - -% final death rate -drlv = dslv+dalv; - -% calculation of specific leaf area in case of exponential growth: -slat = wofost.afgen(WofostPar.SLATB, dvs); -if laiexp < 6 - dteff = max(0,tav-WofostPar.TBASE); % effective temperature - glaiex = laiexp*WofostPar.RGRLAI*dteff*delt; % exponential growth - laiexp = laiexp + glaiex; - - glasol = grlv*slat; % source-limited growth - gla = min(glaiex,glasol); -% gla = max(0,gla); - - slat = gla/grlv; - if isnan(slat) - slat = 0; + + % calculation of new leaf area and weight + lasum = 0.0; + wlv = 0.0; + for i1 = 1:ilvold + lasum = lasum + lv(i1) * sla(i1); + wlv = wlv + lv(i1); end -% if grlv>=0 -% slat = gla/grlv; % the actual sla value -% else -% slat = 0; -% end -end - -% update the information of leave -dslvt = dslv; -i1 = KT; -while (dslvt>0 && i1>=1) % water stress and high lai - if dslvt >= lv(i1) - dslvt = dslvt-lv(i1); - lv(i1) = 0.0; - i1 = i1-1; - else - lv(i1) = lv(i1)-dslvt; - dslvt = 0.0; + lasum = max(0, lasum); + wlv = max(0, wlv); + + %% 7. integrals of the crop + % dry weight of living plant organs + wrt = wrt + gwrt; + wst = wst + gwst; + wso = wso + gwso; + + % total above ground biomass + tadw = wlv + wst + wso; + + % dry weight of dead plant organs (roots,leaves & stems) + dwrt = dwrt + drrt; + dwlv = dwlv + drlv; + dwst = dwst + drst; + + % dry weight of dead and living plant organs + twrt = wrt + dwrt; + twlv = wlv + dwlv; + twst = wst + dwst; + cwdm = twlv + twst + wso; + + % leaf area index + lai = WofostPar.LAIEM + lasum + WofostPar.SSA * wst + WofostPar.SPA * wso; + + Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration + if Anet_sum < 0 && KT > 1 % the consumed biomass are the assimilated one at daytime + lai_delta = lai_delta + grlv * slat; + lai = lai - lai_delta; + elseif Anet_sum >= 0 && Anet >= 0 + Anet_sum = 0; + lai_delta = 0; end -end - -while (lvage(i1)>=WofostPar.SPAN && i1>=1) % leaves older than span die - lv(i1) = 0; - i1 = i1-1; -end - -ilvold = KT; % oldest class with leaves -fysdel = max (0.0d0,(tav-WofostPar.TBASE)/(35.0-WofostPar.TBASE)); % physiologic ageing of leaves per time step - -for i1 = ilvold: -1: 1 % shifting of contents, updating of physiological age - lv(i1+1) = lv(i1); - sla(i1+1) = sla(i1); - lvage(i1+1) = lvage(i1)+fysdel*delt; -end -ilvold = ilvold+1; - -lv(1) = grlv; -sla(1) = slat; -lvage(1) = 0.0; - -% calculation of new leaf area and weight -lasum = 0.0; -wlv = 0.0; -for i1 = 1:ilvold - lasum = lasum + lv(i1)*sla(i1); - wlv = wlv + lv(i1); -end -lasum = max(0,lasum); -wlv = max(0,wlv); - -%% 7. integrals of the crop -% dry weight of living plant organs -wrt = wrt+gwrt; -wst = wst+gwst; -wso = wso+gwso; - -% total above ground biomass -tadw = wlv+wst+wso; - -% dry weight of dead plant organs (roots,leaves & stems) -dwrt = dwrt+drrt; -dwlv = dwlv+drlv; -dwst = dwst+drst; - -% dry weight of dead and living plant organs -twrt = wrt+dwrt; -twlv = wlv+dwlv; -twst = wst+dwst; -cwdm = twlv+twst+wso; - -% leaf area index -lai = WofostPar.LAIEM+lasum+WofostPar.SSA*wst+WofostPar.SPA*wso; - -Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration -if Anet_sum < 0 && KT>1 % the consumed biomass are the assimilated one at daytime - lai_delta = lai_delta + grlv*slat; - lai = lai - lai_delta; -elseif Anet_sum>=0 && Anet>=0 - Anet_sum = 0; - lai_delta = 0; -end - -%% 8. integrals of the crop -crop_output(KT,1) = xyt.t(KT,1); % Day of the year -crop_output(KT,2) = dvs; % Development of stage -crop_output(KT,3) = lai; % LAI -crop_output(KT,4) = ph; % Plant height -crop_output(KT,5) = sfactor; % Water stress -crop_output(KT,6) = wrt; % Dry matter of root -crop_output(KT,7) = wlv; % Dry matter of leaves -crop_output(KT,8) = wst; % Dry matter of stem -crop_output(KT,9) = wso; % Dry matter of organ -crop_output(KT,10) = dwrt; % Dry matter of deadleaves -crop_output(KT,11) = dwlv; % Dry matter of dead stem -crop_output(KT,12) = dwst; % Dry matter of dead organ - -%% 9. Pack updated state variables -state_vars.dvs = dvs; -state_vars.wrt = wrt; -state_vars.wlv = wlv; -state_vars.wst = wst; -state_vars.wso = wso; -state_vars.dwrt = dwrt; -state_vars.dwlv = dwlv; -state_vars.dwst = dwst; -state_vars.sla = sla; -state_vars.lvage = lvage; -state_vars.lv = lv; -state_vars.lasum = lasum; -state_vars.laiexp = laiexp; -state_vars.lai = lai; -state_vars.ph = ph; -state_vars.Anet_sum = Anet_sum; -state_vars.lai_delta = lai_delta; + + %% 8. integrals of the crop + crop_output(KT, 1) = xyt.t(KT, 1); % Day of the year + crop_output(KT, 2) = dvs; % Development of stage + crop_output(KT, 3) = lai; % LAI + crop_output(KT, 4) = ph; % Plant height + crop_output(KT, 5) = sfactor; % Water stress + crop_output(KT, 6) = wrt; % Dry matter of root + crop_output(KT, 7) = wlv; % Dry matter of leaves + crop_output(KT, 8) = wst; % Dry matter of stem + crop_output(KT, 9) = wso; % Dry matter of organ + crop_output(KT, 10) = dwrt; % Dry matter of deadleaves + crop_output(KT, 11) = dwlv; % Dry matter of dead stem + crop_output(KT, 12) = dwst; % Dry matter of dead organ + + %% 9. Pack updated state variables + state_vars.dvs = dvs; + state_vars.wrt = wrt; + state_vars.wlv = wlv; + state_vars.wst = wst; + state_vars.wso = wso; + state_vars.dwrt = dwrt; + state_vars.dwlv = dwlv; + state_vars.dwst = dwst; + state_vars.sla = sla; + state_vars.lvage = lvage; + state_vars.lv = lv; + state_vars.lasum = lasum; + state_vars.laiexp = laiexp; + state_vars.lai = lai; + state_vars.ph = ph; + state_vars.Anet_sum = Anet_sum; + state_vars.lai_delta = lai_delta; end diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 1546c900..465d84c6 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -198,11 +198,11 @@ [rho, tau, rs] = deal(zeros(nwlP + nwlT, 1)); - %% 11. Define plant growth parameters + %% 11. Define plant growth parameters if options.calc_vegetation_dynamic == 1 -% global crop_output + %global crop_output WofostPar = wofost.WofostRead(path_input); - crop_output = zeros(TimeProperties.Dur_tot,12); + crop_output = zeros(TimeProperties.Dur_tot, 12); state_vars = struct(); end @@ -382,7 +382,7 @@ % Will do one timestep in "update mode", and run until the end if in "full run" mode. while KT < endTime KT = KT + 1; % Counting Number of timesteps - fprintf(strcat('\n KT = ',num2str(KT),' \r')); + fprintf(strcat('\n KT = ', num2str(KT), ' \r')); 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. @@ -431,16 +431,16 @@ vi(vmax == telmax) = k; end - % using simulated LAI and PH to substitute prescribed LAI - if options.calc_vegetation_dynamic == 1 && KT >= WofostPar.CSTART && KT <= WofostPar.CEND + % using simulated LAI and PH to substitute prescribed LAI + if options.calc_vegetation_dynamic == 1 && KT >= WofostPar.CSTART && KT <= WofostPar.CEND if KT == WofostPar.CSTART ScopeParameters.LAI(KT) = WofostPar.LAIEM; % initial LAI ScopeParameters.hc(KT) = WofostPar.PHEM; % initial PH - elseif KT > WofostPar.CSTART - ScopeParameters.LAI(KT) = crop_output(KT-1,3); % substitute LAI - ScopeParameters.hc(KT) = crop_output(KT-1,4); % substitute PH + elseif KT > WofostPar.CSTART + ScopeParameters.LAI(KT) = crop_output(KT - 1, 3); % substitute LAI + ScopeParameters.hc(KT) = crop_output(KT - 1, 4); % substitute PH end - end + end [soil, leafbio, canopy, meteo, angles, xyt] = io.select_input(ScopeParameters, SoilVariables.Theta_LL, vi, canopy, options, xyt, soil); if options.simulation ~= 1 @@ -564,15 +564,15 @@ end end end - + % calculate the plant growth process - if options.calc_vegetation_dynamic == 1 && KT >= WofostPar.CSTART && KT <= WofostPar.CEND % vegetation growth process + if options.calc_vegetation_dynamic == 1 && KT >= WofostPar.CSTART && KT <= WofostPar.CEND % vegetation growth process Anet = fluxes.Actot; if isnan(Anet) || Anet < -2 % limit value of Anet Anet = 0; fluxes.Actot = Anet; end - [crop_output,state_vars] = wofost.cropgrowth(crop_output,state_vars,meteo,WofostPar,Anet,WaterStressFactor,xyt,KT,TimeProperties.Dur_tot); + [crop_output, state_vars] = wofost.cropgrowth(crop_output, state_vars, meteo, WofostPar, Anet, WaterStressFactor, xyt, KT, TimeProperties.Dur_tot); end if options.simulation == 2 && telmax > 1 From 22abb53589528bb1665f638758647b0bec1b41bb Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Tue, 7 Jan 2025 17:05:20 +0100 Subject: [PATCH 12/14] Revise the code based on MISS style --- src/+wofost/cropgrowth.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/+wofost/cropgrowth.m b/src/+wofost/cropgrowth.m index e1d4b009..6936d7ab 100644 --- a/src/+wofost/cropgrowth.m +++ b/src/+wofost/cropgrowth.m @@ -222,7 +222,7 @@ % calculation of specific leaf area in case of exponential growth: slat = wofost.afgen(WofostPar.SLATB, dvs); if laiexp < 6 - dteff = max(0, tav - WofostPar.TBASE); % effective temperature + dteff = max(0, tav - WofostPar.TBASE); % effective temperature glaiex = laiexp * WofostPar.RGRLAI * dteff * delt; % exponential growth laiexp = laiexp + glaiex; From 3119425827c339f0967b6e25ed0b0d70bb21a47d Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Wed, 8 Jan 2025 15:57:37 +0100 Subject: [PATCH 13/14] Revise the code to MISS_HIT style 2 --- src/+io/bin_to_csv.m | 4 +- src/+wofost/WofostRead.m | 21 ++- src/+wofost/cropgrowth.asv | 353 ------------------------------------- src/+wofost/cropgrowth.m | 110 ++++++------ src/STEMMUS_SCOPE.m | 2 +- 5 files changed, 68 insertions(+), 422 deletions(-) delete mode 100644 src/+wofost/cropgrowth.asv diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m index 0404a4b1..16c270da 100644 --- a/src/+io/bin_to_csv.m +++ b/src/+io/bin_to_csv.m @@ -94,9 +94,9 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings, %% output the vegetation dynamic results Danyang Yu if options.calc_vegetation_dynamic cropgrowth_names = {'DOY', 'DVS', 'LAI', ... - 'PH', 'Sfactor', 'RootDM', 'LeafDM', 'StemDM', 'OrganDM', 'RootDeath', 'LeafDeath', 'StemDeath',}; + 'PH', 'Sfactor', 'RootDM', 'LeafDM', 'StemDM', 'OrganDM', 'RootDeath', 'LeafDeath', 'StemDeath'}; cropgrowth_units = {'day', '-', 'm2/m2', ... - 'cm', '-', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha'}; + 'cm', '-', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha', 'kg/ha'}; write_output(cropgrowth_names, cropgrowth_units, fnames.cropgrowth_file, n_col.cropgrowth, ns); end diff --git a/src/+wofost/WofostRead.m b/src/+wofost/WofostRead.m index 1166279d..453fb802 100644 --- a/src/+wofost/WofostRead.m +++ b/src/+wofost/WofostRead.m @@ -10,11 +10,11 @@ % Read lines from the file PlantGrowthFile = [path_input, 'plant_growth/CropD.crp']; fid = fopen(PlantGrowthFile); - + % process it line by line while true tline = fgetl(fid); - + % remove empty and note line if ~isempty(tline) & tline(1) ~= '*' % judge whether the string contain the table value @@ -33,24 +33,23 @@ table_value = []; while ~isempty(tline) & tline(1) ~= '*' s = strsplit(strtrim(tline)); - table_value(i,1) = str2num(strtrim(char(s(1)))); - table_value(i,2) = str2num(strtrim(char(s(2)))); + table_value(i, 1) = str2num(strtrim(char(s(1)))); + table_value(i, 2) = str2num(strtrim(char(s(2)))); i = i + 1; tline = fgetl(fid); % disp(tline); end - + wofost.(vname) = table_value; end end - - %end of file - if contains(tline,'* End of .crp file !') - break; - end + + % end of file + if contains(tline, '* End of .crp file !') + break; + end % disp(tline); end - end diff --git a/src/+wofost/cropgrowth.asv b/src/+wofost/cropgrowth.asv deleted file mode 100644 index 271874a1..00000000 --- a/src/+wofost/cropgrowth.asv +++ /dev/null @@ -1,353 +0,0 @@ -function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,WofostPar,Anet,WaterStressFactor,xyt,KT,Dur_tot) - %{ - function cropgrowth.m simulate the growth of vegetation - - authors: Danyang Yu (yudanyang123@gmail.com) - date: 2025/01/01 - reference: SWAP version 3.2. Theory description and user manual (Page 147-163) - - Table of contents of the function - 1. Initilize crop growth paramters - 2. Calculate the phenological development - 3. Calculate the growth rate of different plant organs - 4. Output the crop growth variables - - Input: - Tsum Temperature sum from emergence to the simulated day - tdwi Initial total crop dry weight - laiem Leaf area index at emergence - f(r,l,s,o) The fraction of dry matter to root, leaf, stem, storage - tadw Initial total crop dry weight - lvage leaf age - lv save leaf dry matter at different steps - spa Specific pod area - ssa Specific stem area - sla specific leaf area - Anet_sum the sum of net assimilation - dmi converted dry matter at step KT - admi above-ground dry matter - grrt growth rate of root, similar for the other organs - drrt death rate of root, similar for the other organs - WofostPar Wofost parameters, the definitions are in CropD.crp - - Output: - dvs Developement of stage - lai leaf area index - ph Plant height - wrt dry matter weight of root - wlv dry matter weight of leaves - wst dry matter weight of stem - wso dry matter weight of storage organ - dwrt dry matter weight of dead root - dwlv dry matter weight of dead leaves - dwst dry matter weight of dead stem - %} - - %% 1. initilization - if KT == WofostPar.CSTART - % initial value of crop parameters - Tsum = 0; - dvs = 0; - tdwi = WofostPar.TDWI; - laiem = WofostPar.LAIEM; - ph = 0; - - frtb = WofostPar.FRTB; - fltb = WofostPar.FLTB; - fstb = WofostPar.FSTB; - fotb = WofostPar.FOTB; - - fr = wofost.afgen(frtb, dvs); - fl = wofost.afgen(fltb, dvs); - fs = wofost.afgen(fstb, dvs); - fo = wofost.afgen(fotb, dvs); - - ssa = WofostPar.SSA; - spa = WofostPar.SPA; - - % initial state variables of the crop - tadw = (1 - fr) * tdwi; - wrt = fr * tdwi; - wlv = fl * tadw; - wst = fs * tadw; - wso = fo * tadw; - - dwrt = 0.0; - dwlv = 0.0; - dwst = 0.0; - - sla = zeros(Dur_tot+1, 1); - lvage = zeros(Dur_tot+1, 1); - lv = zeros(Dur_tot+1, 1); - sla(1) = wofost.afgen(WofostPar.SLATB, dvs); - lv(1) = wlv; - lvage(1) = 0.0; - ilvold = 1; - - lasum = laiem; - laiexp = laiem; - glaiex = 0; - laimax = laiem; - lai = lasum + ssa * wst + spa * wso; - Anet_sum = 0; - lai_delta = 0; - - else - % Unpack state variables - dvs = state_vars.dvs; - wrt = state_vars.wrt; - wlv = state_vars.wlv; - wst = state_vars.wst; - wso = state_vars.wso; - - dwrt = state_vars.dwrt; - dwlv = state_vars.dwlv; - dwst = state_vars.dwst; - - sla = state_vars.sla; - lvage = state_vars.lvage; - lv = state_vars.lv; - lasum = state_vars.lasum; - laiexp = state_vars.laiexp; - lai = state_vars.lai; - - ph = state_vars.ph; - Anet_sum = state_vars.Anet_sum; - lai_delta = state_vars.lai_delta; - end - - %% 1. phenological development - delt = WofostPar.TSTEP / 24; - tav = max(0, meteo.Ta); - dtsum = wofost.afgen(WofostPar.DTSMTB, tav); - - if dvs <= 1 % vegetative stage - dvs = dvs + dtsum * delt / WofostPar.TSUMEA; - else - dvs = dvs + dtsum * delt / WofostPar.TSUMAM; - end - - %% 2. dry matter increase - frtb = WofostPar.FRTB; - fltb = WofostPar.FLTB; - fstb = WofostPar.FSTB; - fotb = WofostPar.FOTB; - - fr = wofost.afgen(frtb, dvs); - fl = wofost.afgen(fltb, dvs); - fs = wofost.afgen(fstb, dvs); - fo = wofost.afgen(fotb, dvs); - - fcheck = fr + (fl + fs + fo) * (1.0 - fr) - 1.0; %check on partitions - if abs(fcheck) > 0.0001 - print('The sum of partitioning factors for leaves, stems and storage organs is not equal to one at development stage'); - return; - end - - cvf = 1.0 / ((fl / WofostPar.CVL + fs / WofostPar.CVS + fo / WofostPar.CVO) * (1.0 -fr ) + fr / WofostPar.CVR); - asrc = Anet * 30 / 1000000000 * 10000 * 3600 * WofostPar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; - dmi = cvf * asrc; - - %% 3. Growth rate by root - % root extension - % rr = min(WofostPar.RDM - WofostPar.RD,WofostPar.RRI); - admi = (1.0 - fr) * dmi; - grrt = fr * dmi; - drrt = wrt * wofost.afgen(WofostPar.RDRRTB, dvs) * delt; - gwrt = grrt - drrt; - - %% 4. Growth rate by stem - % growth rate stems - grst = fs * admi; - % death rate stems - drst = wofost.afgen(WofostPar.RDRSTB, dvs) * wst * delt; - % net growth rate stems - gwst = grst - drst; - - % growth of plant height - str = wofost.afgen(WofostPar.STRTB, dvs); - phnew = wst / (WofostPar.STN * WofostPar.STD * pi * str * str) * WofostPar.PHCoeff + WofostPar.PHEM; - if phnew >= ph - ph = phnew; - end - - %% 5. Growth rate by organs - gwso = fo * admi; - - %% 6. Growth rate by leave - % weight of new leaves - grlv = fl * admi; - - % death of leaves due to water stress or high lai - if abs(lai) < 0.1 - dslv1 = 0; - else - sfactor = WaterStressFactor.soil; - dslv1 = wlv * (1.0 - sfactor) * WofostPar.PERDL * delt; - end - - laicr = WofostPar.LAICR; - dslv2 = wlv * max(0, 0.03*(lai - laicr) / laicr); - dslv = max(dslv1, dslv2); - - % death of leaves due to exceeding life span - % first: leaf death due to water stress or high lai is imposed on array - % until no more leaves have to die or all leaves are gone - rest = dslv; - i1 = KT; - iday = ceil(KT*delt); - - while (rest > lv(i1) && i1 >= 1) - rest = rest - lv(i1); - i1 = i1 - 1; - end - - % then: check if some of the remaining leaves are older than span, - % sum their weights - dalv = 0.0; - if (lvage(i1) > WofostPar.SPAN && rest > 0 && i1 >= 1) - dalv = lv(i1) - rest; - rest = 0; - i1 = i1 - 1; - end - - while (i1 >= 1 && lvage(i1) > WofostPar.SPAN) - dalv = dalv + lv(i1); - i1 = i1 - 1; - end - - % final death rate - drlv = dslv + dalv; - - % calculation of specific leaf area in case of exponential growth: - slat = wofost.afgen(WofostPar.SLATB, dvs); - if laiexp < 6 - dteff = max(0, tav - WofostPar.TBASE); % effective temperature - glaiex = laiexp * WofostPar.RGRLAI * dteff * delt; % exponential growth - laiexp = laiexp + glaiex; - - glasol = grlv * slat; % source-limited growth - gla = min(glaiex, glasol); - % gla = max(0,gla); - - slat = gla / grlv; - if isnan(slat) - slat = 0; - end - % if grlv>=0 - % slat = gla/grlv; % the actual sla value - % else - % slat = 0; - % end - end - - % update the information of leave - dslvt = dslv; - i1 = KT; - while (dslvt>0 && i1>=1) % water stress and high lai - if dslvt >= lv(i1) - dslvt = dslvt-lv(i1); - lv(i1) = 0.0; - i1 = i1-1; - else - lv(i1) = lv(i1)-dslvt; - dslvt = 0.0; - end - end - - while (lvage(i1) >= WofostPar.SPAN && i1 >= 1) % leaves older than span die - lv(i1) = 0; - i1 = i1 - 1; - end - - ilvold = KT; % oldest class with leaves - fysdel = max(0.0d0, (tav - WofostPar.TBASE) / (35.0 - WofostPar.TBASE)); % physiologic ageing of leaves per time step - - for i1 = ilvold: -1: 1 % shifting of contents, updating of physiological age - lv(i1+1) = lv(i1); - sla(i1+1) = sla(i1); - lvage(i1+1) = lvage(i1) + fysdel * delt; - end - ilvold = ilvold + 1; - - lv(1) = grlv; - sla(1) = slat; - lvage(1) = 0.0; - - % calculation of new leaf area and weight - lasum = 0.0; - wlv = 0.0; - for i1 = 1:ilvold - lasum = lasum + lv(i1)*sla(i1); - wlv = wlv + lv(i1); - end - lasum = max(0,lasum); - wlv = max(0,wlv); - - %% 7. integrals of the crop - % dry weight of living plant organs - wrt = wrt+gwrt; - wst = wst+gwst; - wso = wso+gwso; - - % total above ground biomass - tadw = wlv+wst+wso; - - % dry weight of dead plant organs (roots,leaves & stems) - dwrt = dwrt+drrt; - dwlv = dwlv+drlv; - dwst = dwst+drst; - - % dry weight of dead and living plant organs - twrt = wrt+dwrt; - twlv = wlv+dwlv; - twst = wst+dwst; - cwdm = twlv+twst+wso; - - % leaf area index - lai = WofostPar.LAIEM+lasum+WofostPar.SSA*wst+WofostPar.SPA*wso; - - Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration - if Anet_sum < 0 && KT>1 % the consumed biomass are the assimilated one at daytime - lai_delta = lai_delta + grlv*slat; - lai = lai - lai_delta; - elseif Anet_sum>=0 && Anet>=0 - Anet_sum = 0; - lai_delta = 0; - end - - %% 8. integrals of the crop - crop_output(KT,1) = xyt.t(KT,1); % Day of the year - crop_output(KT,2) = dvs; % Development of stage - crop_output(KT,3) = lai; % LAI - crop_output(KT,4) = ph; % Plant height - crop_output(KT,5) = sfactor; % Water stress - crop_output(KT,6) = wrt; % Dry matter of root - crop_output(KT,7) = wlv; % Dry matter of leaves - crop_output(KT,8) = wst; % Dry matter of stem - crop_output(KT,9) = wso; % Dry matter of organ - crop_output(KT,10) = dwrt; % Dry matter of deadleaves - crop_output(KT,11) = dwlv; % Dry matter of dead stem - crop_output(KT,12) = dwst; % Dry matter of dead organ - - %% 9. Pack updated state variables - state_vars.dvs = dvs; - state_vars.wrt = wrt; - state_vars.wlv = wlv; - state_vars.wst = wst; - state_vars.wso = wso; - state_vars.dwrt = dwrt; - state_vars.dwlv = dwlv; - state_vars.dwst = dwst; - state_vars.sla = sla; - state_vars.lvage = lvage; - state_vars.lv = lv; - state_vars.lasum = lasum; - state_vars.laiexp = laiexp; - state_vars.lai = lai; - state_vars.ph = ph; - state_vars.Anet_sum = Anet_sum; - state_vars.lai_delta = lai_delta; -end - - diff --git a/src/+wofost/cropgrowth.m b/src/+wofost/cropgrowth.m index 6936d7ab..c959dd01 100644 --- a/src/+wofost/cropgrowth.m +++ b/src/+wofost/cropgrowth.m @@ -1,17 +1,17 @@ function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,WofostPar,Anet,WaterStressFactor,xyt,KT,Dur_tot) %{ function cropgrowth.m simulate the growth of vegetation - + authors: Danyang Yu (yudanyang123@gmail.com) date: 2025/01/01 reference: SWAP version 3.2. Theory description and user manual (Page 147-163) - + Table of contents of the function 1. Initilize crop growth paramters 2. Calculate the phenological development 3. Calculate the growth rate of different plant organs 4. Output the crop growth variables - + Input: Tsum Temperature sum from emergence to the simulated day tdwi Initial total crop dry weight @@ -29,7 +29,7 @@ grrt growth rate of root, similar for the other organs drrt death rate of root, similar for the other organs WofostPar Wofost parameters, the definitions are in CropD.crp - + Output: dvs Developement of stage lai leaf area index @@ -42,7 +42,7 @@ dwlv dry matter weight of dead leaves dwst dry matter weight of dead stem %} - + %% 1. initilization if KT == WofostPar.CSTART % initial value of crop parameters @@ -51,31 +51,31 @@ tdwi = WofostPar.TDWI; laiem = WofostPar.LAIEM; ph = 0; - + frtb = WofostPar.FRTB; fltb = WofostPar.FLTB; fstb = WofostPar.FSTB; fotb = WofostPar.FOTB; - + fr = wofost.afgen(frtb, dvs); fl = wofost.afgen(fltb, dvs); fs = wofost.afgen(fstb, dvs); fo = wofost.afgen(fotb, dvs); - + ssa = WofostPar.SSA; spa = WofostPar.SPA; - + % initial state variables of the crop tadw = (1 - fr) * tdwi; wrt = fr * tdwi; wlv = fl * tadw; wst = fs * tadw; wso = fo * tadw; - + dwrt = 0.0; dwlv = 0.0; dwst = 0.0; - + sla = zeros(Dur_tot+1, 1); lvage = zeros(Dur_tot+1, 1); lv = zeros(Dur_tot+1, 1); @@ -83,7 +83,7 @@ lv(1) = wlv; lvage(1) = 0.0; ilvold = 1; - + lasum = laiem; laiexp = laiem; glaiex = 0; @@ -91,7 +91,7 @@ lai = lasum + ssa * wst + spa * wso; Anet_sum = 0; lai_delta = 0; - + else % Unpack state variables dvs = state_vars.dvs; @@ -99,55 +99,55 @@ wlv = state_vars.wlv; wst = state_vars.wst; wso = state_vars.wso; - + dwrt = state_vars.dwrt; dwlv = state_vars.dwlv; dwst = state_vars.dwst; - + sla = state_vars.sla; lvage = state_vars.lvage; lv = state_vars.lv; lasum = state_vars.lasum; laiexp = state_vars.laiexp; lai = state_vars.lai; - + ph = state_vars.ph; Anet_sum = state_vars.Anet_sum; lai_delta = state_vars.lai_delta; end - + %% 1. phenological development delt = WofostPar.TSTEP / 24; tav = max(0, meteo.Ta); dtsum = wofost.afgen(WofostPar.DTSMTB, tav); - + if dvs <= 1 % vegetative stage dvs = dvs + dtsum * delt / WofostPar.TSUMEA; else dvs = dvs + dtsum * delt / WofostPar.TSUMAM; end - + %% 2. dry matter increase frtb = WofostPar.FRTB; fltb = WofostPar.FLTB; fstb = WofostPar.FSTB; fotb = WofostPar.FOTB; - + fr = wofost.afgen(frtb, dvs); fl = wofost.afgen(fltb, dvs); fs = wofost.afgen(fstb, dvs); fo = wofost.afgen(fotb, dvs); - + fcheck = fr + (fl + fs + fo) * (1.0 - fr) - 1.0; %check on partitions if abs(fcheck) > 0.0001 print('The sum of partitioning factors for leaves, stems and storage organs is not equal to one at development stage'); return; end - + cvf = 1.0 / ((fl / WofostPar.CVL + fs / WofostPar.CVS + fo / WofostPar.CVO) * (1.0 -fr ) + fr / WofostPar.CVR); asrc = Anet * 30 / 1000000000 * 10000 * 3600 * WofostPar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; dmi = cvf * asrc; - + %% 3. Growth rate by root % root extension % rr = min(WofostPar.RDM - WofostPar.RD,WofostPar.RRI); @@ -155,7 +155,7 @@ grrt = fr * dmi; drrt = wrt * wofost.afgen(WofostPar.RDRRTB, dvs) * delt; gwrt = grrt - drrt; - + %% 4. Growth rate by stem % growth rate stems grst = fs * admi; @@ -163,21 +163,21 @@ drst = wofost.afgen(WofostPar.RDRSTB, dvs) * wst * delt; % net growth rate stems gwst = grst - drst; - + % growth of plant height str = wofost.afgen(WofostPar.STRTB, dvs); phnew = wst / (WofostPar.STN * WofostPar.STD * pi * str * str) * WofostPar.PHCoeff + WofostPar.PHEM; if phnew >= ph ph = phnew; end - + %% 5. Growth rate by organs gwso = fo * admi; - + %% 6. Growth rate by leave % weight of new leaves grlv = fl * admi; - + % death of leaves due to water stress or high lai if abs(lai) < 0.1 dslv1 = 0; @@ -185,23 +185,23 @@ sfactor = WaterStressFactor.soil; dslv1 = wlv * (1.0 - sfactor) * WofostPar.PERDL * delt; end - + laicr = WofostPar.LAICR; dslv2 = wlv * max(0, 0.03*(lai - laicr) / laicr); dslv = max(dslv1, dslv2); - + % death of leaves due to exceeding life span % first: leaf death due to water stress or high lai is imposed on array % until no more leaves have to die or all leaves are gone rest = dslv; i1 = KT; iday = ceil(KT*delt); - + while (rest > lv(i1) && i1 >= 1) rest = rest - lv(i1); i1 = i1 - 1; end - + % then: check if some of the remaining leaves are older than span, % sum their weights dalv = 0.0; @@ -210,26 +210,26 @@ rest = 0; i1 = i1 - 1; end - + while (i1 >= 1 && lvage(i1) > WofostPar.SPAN) dalv = dalv + lv(i1); i1 = i1 - 1; end - + % final death rate drlv = dslv + dalv; - + % calculation of specific leaf area in case of exponential growth: slat = wofost.afgen(WofostPar.SLATB, dvs); if laiexp < 6 dteff = max(0, tav - WofostPar.TBASE); % effective temperature glaiex = laiexp * WofostPar.RGRLAI * dteff * delt; % exponential growth laiexp = laiexp + glaiex; - + glasol = grlv * slat; % source-limited growth gla = min(glaiex, glasol); % gla = max(0,gla); - + slat = gla / grlv; if isnan(slat) slat = 0; @@ -240,7 +240,7 @@ % slat = 0; % end end - + % update the information of leave dslvt = dslv; i1 = KT; @@ -254,26 +254,26 @@ dslvt = 0.0; end end - + while (lvage(i1) >= WofostPar.SPAN && i1 >= 1) % leaves older than span die lv(i1) = 0; i1 = i1 - 1; end - + ilvold = KT; % oldest class with leaves fysdel = max(0.0d0, (tav - WofostPar.TBASE) / (35.0 - WofostPar.TBASE)); % physiologic ageing of leaves per time step - - for i1 = ilvold: -1: 1 % shifting of contents, updating of physiological age - lv(i1+1) = lv(i1); - sla(i1+1) = sla(i1); - lvage(i1+1) = lvage(i1) + fysdel * delt; + + for i1 = ilvold:-1:1 % shifting of contents, updating of physiological age + lv(i1 + 1) = lv(i1); + sla(i1 + 1) = sla(i1); + lvage(i1 + 1) = lvage(i1) + fysdel * delt; end ilvold = ilvold + 1; - + lv(1) = grlv; sla(1) = slat; lvage(1) = 0.0; - + % calculation of new leaf area and weight lasum = 0.0; wlv = 0.0; @@ -283,30 +283,30 @@ end lasum = max(0, lasum); wlv = max(0, wlv); - + %% 7. integrals of the crop % dry weight of living plant organs wrt = wrt + gwrt; wst = wst + gwst; wso = wso + gwso; - + % total above ground biomass tadw = wlv + wst + wso; - + % dry weight of dead plant organs (roots,leaves & stems) dwrt = dwrt + drrt; dwlv = dwlv + drlv; dwst = dwst + drst; - + % dry weight of dead and living plant organs twrt = wrt + dwrt; twlv = wlv + dwlv; twst = wst + dwst; cwdm = twlv + twst + wso; - + % leaf area index lai = WofostPar.LAIEM + lasum + WofostPar.SSA * wst + WofostPar.SPA * wso; - + Anet_sum = Anet_sum + Anet; % assume LAI not changed with the respiration if Anet_sum < 0 && KT > 1 % the consumed biomass are the assimilated one at daytime lai_delta = lai_delta + grlv * slat; @@ -315,7 +315,7 @@ Anet_sum = 0; lai_delta = 0; end - + %% 8. integrals of the crop crop_output(KT, 1) = xyt.t(KT, 1); % Day of the year crop_output(KT, 2) = dvs; % Development of stage @@ -329,7 +329,7 @@ crop_output(KT, 10) = dwrt; % Dry matter of deadleaves crop_output(KT, 11) = dwlv; % Dry matter of dead stem crop_output(KT, 12) = dwst; % Dry matter of dead organ - + %% 9. Pack updated state variables state_vars.dvs = dvs; state_vars.wrt = wrt; diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index b90c8d45..9c1cb299 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -200,7 +200,7 @@ %% 11. Define plant growth parameters if options.calc_vegetation_dynamic == 1 - %global crop_output + % global crop_output WofostPar = wofost.WofostRead(path_input); crop_output = zeros(TimeProperties.Dur_tot, 12); state_vars = struct(); From bee55fc9e06fd42890f232dbb992022a5552becf Mon Sep 17 00:00:00 2001 From: Danyang Yu Date: Wed, 8 Jan 2025 16:31:34 +0100 Subject: [PATCH 14/14] Revise code to MISS_HIT style 3 --- src/+io/bin_to_csv.m | 2 +- src/+wofost/WofostRead.m | 3 +-- src/+wofost/afgen.m | 7 +++-- src/+wofost/cropgrowth.m | 58 +++++++++++++++++++--------------------- 4 files changed, 33 insertions(+), 37 deletions(-) diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m index 16c270da..fb892666 100644 --- a/src/+io/bin_to_csv.m +++ b/src/+io/bin_to_csv.m @@ -232,4 +232,4 @@ function write_output(header, units, bin_path, f_n_col, ns, not_header) fprintf(f_csv, '%d,', out_2d(k, 1:end - 1)); fprintf(f_csv, '%d\n', out_2d(k, end)); % saves from extra comma end -end \ No newline at end of file +end diff --git a/src/+wofost/WofostRead.m b/src/+wofost/WofostRead.m index 453fb802..cd871f8f 100644 --- a/src/+wofost/WofostRead.m +++ b/src/+wofost/WofostRead.m @@ -46,10 +46,9 @@ % end of file if contains(tline, '* End of .crp file !') - break; + break end % disp(tline); end end - diff --git a/src/+wofost/afgen.m b/src/+wofost/afgen.m index 5e96ec83..25520487 100644 --- a/src/+wofost/afgen.m +++ b/src/+wofost/afgen.m @@ -1,17 +1,16 @@ -function [output] = afgen(table,x) +function [output] = afgen(table, x) %{ function afgen.m intercept the value in the table authors: Danyang Yu (yudanyang123@gmail.com) date: 2025/01/01 %} - + if x < table(1, 1) || x > table(end, 1) print('the value of developemnt stage is mistaken'); else dvslist = table(:, 1).'; loc = discretize([x], [-Inf dvslist Inf]); - slope = (table(loc, 2) - table((loc-1), 2)) / (table(loc, 1) - table((loc - 1), 1)); + slope = (table(loc, 2) - table(loc - 1, 2)) / (table(loc, 1) - table(loc - 1, 1)); output = table(loc - 1, 2) + slope * (x - table(loc - 1, 1)); end end - diff --git a/src/+wofost/cropgrowth.m b/src/+wofost/cropgrowth.m index c959dd01..9ed2ed2a 100644 --- a/src/+wofost/cropgrowth.m +++ b/src/+wofost/cropgrowth.m @@ -1,4 +1,4 @@ -function [crop_output,state_vars] = cropgrowth(crop_output,state_vars,meteo,WofostPar,Anet,WaterStressFactor,xyt,KT,Dur_tot) +function [crop_output, state_vars] = cropgrowth(crop_output, state_vars, meteo, WofostPar, Anet, WaterStressFactor, xyt, KT, Dur_tot) %{ function cropgrowth.m simulate the growth of vegetation @@ -33,7 +33,7 @@ Output: dvs Developement of stage lai leaf area index - ph Plant height + ph Plant height wrt dry matter weight of root wlv dry matter weight of leaves wst dry matter weight of stem @@ -56,7 +56,7 @@ fltb = WofostPar.FLTB; fstb = WofostPar.FSTB; fotb = WofostPar.FOTB; - + fr = wofost.afgen(frtb, dvs); fl = wofost.afgen(fltb, dvs); fs = wofost.afgen(fstb, dvs); @@ -76,16 +76,16 @@ dwlv = 0.0; dwst = 0.0; - sla = zeros(Dur_tot+1, 1); - lvage = zeros(Dur_tot+1, 1); - lv = zeros(Dur_tot+1, 1); + sla = zeros(Dur_tot + 1, 1); + lvage = zeros(Dur_tot + 1, 1); + lv = zeros(Dur_tot + 1, 1); sla(1) = wofost.afgen(WofostPar.SLATB, dvs); lv(1) = wlv; lvage(1) = 0.0; ilvold = 1; - lasum = laiem; - laiexp = laiem; + lasum = laiem; + laiexp = laiem; glaiex = 0; laimax = laiem; lai = lasum + ssa * wst + spa * wso; @@ -138,13 +138,13 @@ fs = wofost.afgen(fstb, dvs); fo = wofost.afgen(fotb, dvs); - fcheck = fr + (fl + fs + fo) * (1.0 - fr) - 1.0; %check on partitions + fcheck = fr + (fl + fs + fo) * (1.0 - fr) - 1.0; % check on partitions if abs(fcheck) > 0.0001 print('The sum of partitioning factors for leaves, stems and storage organs is not equal to one at development stage'); - return; + return end - cvf = 1.0 / ((fl / WofostPar.CVL + fs / WofostPar.CVS + fo / WofostPar.CVO) * (1.0 -fr ) + fr / WofostPar.CVR); + cvf = 1.0 / ((fl / WofostPar.CVL + fs / WofostPar.CVS + fo / WofostPar.CVO) * (1.0 - fr) + fr / WofostPar.CVR); asrc = Anet * 30 / 1000000000 * 10000 * 3600 * WofostPar.TSTEP; % [umol m-2 s-1] to [kg CH2O ha-1]; dmi = cvf * asrc; @@ -163,7 +163,7 @@ drst = wofost.afgen(WofostPar.RDRSTB, dvs) * wst * delt; % net growth rate stems gwst = grst - drst; - + % growth of plant height str = wofost.afgen(WofostPar.STRTB, dvs); phnew = wst / (WofostPar.STN * WofostPar.STD * pi * str * str) * WofostPar.PHCoeff + WofostPar.PHEM; @@ -187,7 +187,7 @@ end laicr = WofostPar.LAICR; - dslv2 = wlv * max(0, 0.03*(lai - laicr) / laicr); + dslv2 = wlv * max(0, 0.03 * (lai - laicr) / laicr); dslv = max(dslv1, dslv2); % death of leaves due to exceeding life span @@ -195,9 +195,9 @@ % until no more leaves have to die or all leaves are gone rest = dslv; i1 = KT; - iday = ceil(KT*delt); + iday = ceil(KT * delt); - while (rest > lv(i1) && i1 >= 1) + while rest > lv(i1) && i1 >= 1 rest = rest - lv(i1); i1 = i1 - 1; end @@ -205,13 +205,13 @@ % then: check if some of the remaining leaves are older than span, % sum their weights dalv = 0.0; - if (lvage(i1) > WofostPar.SPAN && rest > 0 && i1 >= 1) + if lvage(i1) > WofostPar.SPAN && rest > 0 && i1 >= 1 dalv = lv(i1) - rest; rest = 0; i1 = i1 - 1; end - while (i1 >= 1 && lvage(i1) > WofostPar.SPAN) + while i1 >= 1 && lvage(i1) > WofostPar.SPAN dalv = dalv + lv(i1); i1 = i1 - 1; end @@ -228,34 +228,34 @@ glasol = grlv * slat; % source-limited growth gla = min(glaiex, glasol); - % gla = max(0,gla); + % gla = max(0,gla); slat = gla / grlv; if isnan(slat) slat = 0; end - % if grlv>=0 - % slat = gla/grlv; % the actual sla value - % else - % slat = 0; - % end + % if grlv>=0 + % slat = gla/grlv; % the actual sla value + % else + % slat = 0; + % end end % update the information of leave dslvt = dslv; i1 = KT; - while (dslvt>0 && i1>=1) % water stress and high lai + while dslvt > 0 && i1 >= 1 % water stress and high lai if dslvt >= lv(i1) - dslvt = dslvt-lv(i1); + dslvt = dslvt - lv(i1); lv(i1) = 0.0; - i1 = i1-1; + i1 = i1 - 1; else - lv(i1) = lv(i1)-dslvt; + lv(i1) = lv(i1) - dslvt; dslvt = 0.0; end end - while (lvage(i1) >= WofostPar.SPAN && i1 >= 1) % leaves older than span die + while lvage(i1) >= WofostPar.SPAN && i1 >= 1 % leaves older than span die lv(i1) = 0; i1 = i1 - 1; end @@ -349,5 +349,3 @@ state_vars.Anet_sum = Anet_sum; state_vars.lai_delta = lai_delta; end - -