diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 00000000..a277f006 Binary files /dev/null and b/.DS_Store differ 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/example_data/.DS_Store b/example_data/.DS_Store new file mode 100644 index 00000000..5392c9fa Binary files /dev/null and b/example_data/.DS_Store differ diff --git a/example_data/input_data.xlsx b/example_data/input_data.xlsx new file mode 100644 index 00000000..99b631a2 Binary files /dev/null and b/example_data/input_data.xlsx differ diff --git a/example_data/input_soilLayThick.csv b/example_data/input_soilLayThick.csv index bca17029..0460d3d9 100644 --- a/example_data/input_soilLayThick.csv +++ b/example_data/input_soilLayThick.csv @@ -1,61 +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, +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/example_data/plant_growth/CropD.crp b/example_data/plant_growth/CropD.crp new file mode 100644 index 00000000..0e4e9fce --- /dev/null +++ b/example_data/plant_growth/CropD.crp @@ -0,0 +1,200 @@ +*********************************************************************************************** +* Filename: Maize.CRP +* Contents: SWAP 3.2 - Data for detailed crop model +* Authors: Danyang Yu (yudanyang123@gmail.com) +*********************************************************************************************** +*c Grain maize (Zea mays L.) +*********************************************************************************************** + +*** PLANT GROWTH SECTION *** + +*********************************************************************************************** +* Part 1: Crop development + + CSTART = 1 ! Start time step of crop growth + CEND = 9073 ! End time step of crop growth + TSTEP = 0.5 ! Time step of the simulation, [0..1 hour] + TSUMEA = 900.00 ! Temperature sum from emergence to anthesis, [0..10000 C, R] + TSUMAM = 980.00 ! Temperature sum from anthesis to maturity [0..10000 C, R] + +* List increase in temperature sum [DTSM, 0..60 C, R] as function of average temp. [TAV, 0..100 C, R] +* The intermediate values obtained through linear interpolation. For example, DSTM=DSTM1+(TAV-TAV1)/(TAV2-TAV1)*(DSTM2-DSTM1) + +* TAV DTSM (maximum 15 records) + DTSMTB = + 0.00 0.00 + 10.00 0.00 + 30.00 20.00 + 45.00 20.00 +* End of Table + + DVSEND = 2.00 ! development stage at harvest [-] +*********************************************************************************************** + +*********************************************************************************************** +* Part 2: Initial values + + TDWI = 20.000 ! Initial total crop dry weight [0..10000 kg/ha, R] + LAIEM = 0.1 ! Leaf area index at emergence [0..10 m2/m2, R] + PHEM = 0.1 ! plant height at emergence [0..10 m, R] + RGRLAI = 0.05! Maximum relative increase in LAI [0..1 m2/m2/d, R] +*********************************************************************************************** + +*********************************************************************************************** +* Part 3: Green surface area + + SPA = 0.0000 ! Specific pod area [0..1 ha/kg, R] + SSA = 0.0000 ! Specific stem area [0..1 ha/kg, R] + SPAN = 30.00 ! Life span under leaves under optimum conditions, [0..366 d, R] + LAICR = 10.00 ! The criteriation of LAI which affect the potential death rate due to self-shading, [0..10 m2/m2, R] + TBASE = 10.00 ! Lower threshold temperature for ageing of leaves ,[-10..30 C, R] + +* List specific leaf area [SLA, 0..1 ha/kg, R] as function of development of stage [DVS, 0..2, R] +* The intermediate values obtained through linear interpolation. Similar with the first table + +* DVS SLA (maximum 15 records) + SLATB = + 0.0000 0.0015 + 0.3000 0.0015 + 1.0000 0.0015 + 1.2500 0.0015 + 2.0000 0.0015 +* End of Table +*********************************************************************************************** + +*********************************************************************************************** +* Part 4: Plant height growth + + PHOpt = 0 ! Options of PH simulations; 0 using the prescribed values, and 1 for simulation + STN = 40000 ! Number of stems [1..10000 /ha, I] + STD = 25.0 ! Density of stems [0.1..10 kg/m3, R] + PHCoeff = 1.22 ! Coeffcient of plant height, [1..2 , R] + +* List specific stems radius [STR, 0..1 m, R] as function of devel. stage [DVS, 0..2, R] +* The intermediate values obtained through linear interpolation. Similar with the first table + +* DVS STR (maximum 15 records) + STRTB = + 0.0000 0.020 + 0.6000 0.021 + 1.3000 0.025 + 2.0000 0.025 +* End of Table +*********************************************************************************************** + +*********************************************************************************************** +* Part 5: Conversion of assimilates into biomass +* + CVL = 0.6900 ! Efficiency of conversion into leaves, [0..1 kg/kg, R] + CVO = 0.7500 ! Efficiency of conversion into storage organs, [0..1 kg/kg, R] + CVR = 0.6900 ! Efficiency of conversion into roots, [0..1 kg/kg, R] + CVS = 0.7200 ! Efficiency of conversion into stems, [0..1 kg/kg, R] +*********************************************************************************************** + +*********************************************************************************************** +* Part 6: Partitioning + +* List fraction of total dry matter increase partitioned to the roots [FR, kg/kg, R] +* as function of development stage [DVS, 0..2 -, R] +* The intermediate values obtained through linear interpolation. Similar with the first table + +* DVS FR (maximum 15 records) + FRTB = + 0.00 0.90 + 0.30 0.50 + 1.00 0.00 + 2.00 0.00 +* End of table + +* List fraction of total above ground dry matter incr. part. to the leaves [FL, kg/kg, R] +* as function of development stage [DVS, 0..2 -, R] +* The intermediate values obtained through linear interpolation. Similar with the first table + +* DVS FL (maximum 15 records) + FLTB = + 0.0000 0.6000 + 0.7500 0.5000 + 0.8000 0.1000 + 0.9000 0.1000 + 1.0000 0.1000 + 1.0500 0.0000 + 2.0000 0.0000 +* End of table + +* List fraction of total above ground dry matter incr. part. to the stems [FS, kg/kg, R] +* as function of development stage [DVS, 0..2 -, R] +* The intermediate values obtained through linear interpolation. Similar with the first table + +* DVS FS (maximum 15 records) + FSTB = + 0.0000 0.4000 + 0.7500 0.5000 + 0.8000 0.9000 + 0.9000 0.9000 + 1.0000 0.9000 + 1.0500 0.0000 + 2.0000 0.0000 +* End of table + +* List fraction of total above ground dry matter incr. part. to the st. organs [FO, kg/kg, R] +* as function of development stage [DVS, 0..2 -, R] +* The intermediate values obtained through linear interpolation. Similar with the first table + +* DVS FO (maximum 15 records) + FOTB = + 0.0000 0.0000 + 0.9000 0.0000 + 1.0000 0.0000 + 1.0500 1.0000 + 2.0000 1.0000 +* End of table +*********************************************************************************************** + +*********************************************************************************************** +* Part 7: Death rates + + PERDL = 0.010 ! Maximum rel. death rate of leaves due to water stress [0..3 /d, R] + +* List relative death rates of roots [RDRR, kg/kg/d] as function of dev. stage [DVS, 0..2 -, R] +* The intermediate values obtained through linear interpolation. Similar with the first table + +* DVS RDRR (maximum 15 records) + RDRRTB = + 0.0000 0.0000 + 1.5000 0.0000 + 1.5001 0.0200 + 2.0000 0.0200 +* End of table + +* List relative death rates of stems [RDRS, kg/kg/d] as function of dev. stage [DVS, 0..2 -, R] +* The intermediate values obtained through linear interpolation. Similar with the first table + +* DVS RDRS (maximum 15 records) + RDRSTB = + 0.0000 0.0000 + 1.0000 0.0000 + 1.2500 0.0020 + 2.0000 0.0020 +* End of table +*********************************************************************************************** + +*********************************************************************************************** +* Part 8: Root density distribution and root growth +* +* List relative root density [RDC, 0..1 -, R], as function of rel. rooting depth [RD, 0..1 -, R]: +* The intermediate values obtained through linear interpolation. Similar with the first table + +* RD RDC (maximum 11 records) + RDCTB = + 0.00 1.00 + 1.00 1.00 +* End of table +* + RDI = 10.00 ! Initial rooting depth, [0..1000 cm, R] + RRI = 1.20 ! Maximum daily increase in rooting depth, [0..100 cm/d, R] + RDC = 80.00 ! Maximum rooting depth crop/cultivar, [0..1000 cm, R] +* +************************************************************************************ + + +* End of .crp file ! \ No newline at end of file diff --git a/src/+init/setBoundaryCondition.m b/src/+init/setBoundaryCondition.m index f65b0c84..b2520967 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/bin_to_csv.m b/src/+io/bin_to_csv.m index 5bc13136..fb892666 100644 --- a/src/+io/bin_to_csv.m +++ b/src/+io/bin_to_csv.m @@ -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); diff --git a/src/+io/create_output_files_binary.m b/src/+io/create_output_files_binary.m index 76695055..0e5e2dbc 100644 --- a/src/+io/create_output_files_binary.m +++ b/src/+io/create_output_files_binary.m @@ -108,6 +108,11 @@ 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); fclose("all"); diff --git a/src/+io/loadSoilInitialValues.m b/src/+io/loadSoilInitialValues.m index 94aaa06f..63285173 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/+io/output_data_binary.m b/src/+io/output_data_binary.m index d5a918d4..ab07a6a4 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/+wofost/WofostRead.m b/src/+wofost/WofostRead.m new file mode 100644 index 00000000..cd871f8f --- /dev/null +++ b/src/+wofost/WofostRead.m @@ -0,0 +1,54 @@ +function [wofost] = WofostRead(path_input) + %{ + 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']; + 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/+wofost/afgen.m b/src/+wofost/afgen.m new file mode 100644 index 00000000..25520487 --- /dev/null +++ b/src/+wofost/afgen.m @@ -0,0 +1,16 @@ +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 +end diff --git a/src/+wofost/cropgrowth.m b/src/+wofost/cropgrowth.m new file mode 100644 index 00000000..9ed2ed2a --- /dev/null +++ b/src/+wofost/cropgrowth.m @@ -0,0 +1,351 @@ +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/STEMMUS_SCOPE.asv b/src/STEMMUS_SCOPE.asv new file mode 100644 index 00000000..e43667f5 --- /dev/null +++ b/src/STEMMUS_SCOPE.asv @@ -0,0 +1,870 @@ +%% STEMMUS-SCOPE.m (script) +% STEMMUS-SCOPE is a model for Integrated modeling of canopy photosynthesis, fluorescence, +% and the transfer of energy, mass, and momentum in the soil-plant-atmosphere continuum +% Copyright (C) 2021 Yunfei Wang, Lianyu Yu, Yijian Zeng, Christiaan Van der Tol, Bob Su +% Contact: y.wang-3@utwente.nl; l.yu@utwente.nl; y.zeng@utwente.nl; c.vandertol@utwente.nl; z.su@utwente.nl +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . + +% Load in required Octave packages if STEMMUS-SCOPE is being run in Octave: +if exist('OCTAVE_VERSION', 'builtin') ~= 0 + disp('Loading Octave packages...'); + pkg load statistics io; +end + +% set CFG to a path if it is not defined +if exist('CFG', 'var') == 0 + CFG = '../config_file_crib.txt'; +end + +% set runMode to "full" if it is not defined +if exist('runMode', 'var') == 0 + runMode = "full"; + bmiMode = "none"; +end + +% Initialization routine +start_time = clock; +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, FullCSVfiles] = io.read_config(CFG); + + % Prepare forcing and soil data + [SiteProperties, SoilProperties, TimeProperties] = io.prepareInputData(InputPath); + + % Load model settings: replacing "run Constants" + ModelSettings = io.getModelSettings(); + + % Load soil layers settings + soilLayersFile = [InputPath, 'input_soilLayThick.csv']; + if isfile(soilLayersFile) + SoilLayerSettings = io.readSoilLayerSettings(soilLayersFile); + % Override the default settings + fields = fieldnames(SoilLayerSettings); + for i = 1:numel(fields) + ModelSettings.(fields{i}) = SoilLayerSettings.(fields{i}); + end + end + + % Calculate other ModelSettings + ModelSettings.NN = ModelSettings.NL + 1; % Number of nodes; + ModelSettings.mN = ModelSettings.NN + 1; + ModelSettings.mL = ModelSettings.NL + 1; % Number of elements. Preventing the exceeds of arrays size; + ModelSettings.nD = 2; + + % Load groundwater settings + GroundwaterSettings = groundwater.initializeGroundwaterSettings(); + GroundwaterSettings.soilThick = io.calculateSoilLayerThickness(ModelSettings); + + % load forcing data + ForcingData = io.loadForcingData(InputPath, TimeProperties, SoilProperties, ModelSettings.Tot_Depth, GroundwaterSettings); + + % Get initial values + InitialValues = init.defineInitialValues(TimeProperties.Dur_tot, ModelSettings); + T_CRIT = InitialValues.T_CRIT; + TT_CRIT = InitialValues.TT_CRIT; + hOLD = InitialValues.hOLD; + TOLD = InitialValues.TOLD; + SAVE = InitialValues.SAVE; + P_gOLD = InitialValues.P_gOLD; + gwfluxes = struct(); + Evap = []; + Trap = []; + RWUs = []; + RWUg = []; + + %% 1. define Constants + Constants = io.define_constants(); + + RTB = 1000; % initial root total biomass (g m-2) + % Rl used in ebal + [Rl, Ztot] = Initial_root_biomass(RTB, ModelSettings.DeltZ_R, ModelSettings.rroot, ModelSettings.ML, SiteProperties.landcoverClass(1)); + + %% 2. simulation options + path_input = InputPath; % path of all inputs + path_of_code = cd; + + parameter_file = {'input_data.xlsx'}; + options = io.readStructFromExcel([path_input char(parameter_file)], 'options', 2, 1); + + if options.simulation > 2 || options.simulation < 0 + fprintf('\n simulation option should be between 0 and 2 \r'); + return + end + + %% 3. file names + [dummy, X] = xlsread([path_input char(parameter_file)], 'filenames'); + j = find(~strcmp(X(:, 2), {''})); + X = X(j, 1:end); + + F = struct('FileID', {'Simulation_Name', 'soil_file', 'leaf_file', 'atmos_file'... + 'Dataset_dir', 't_file', 'year_file', 'Rin_file', 'Rli_file' ... + 'p_file', 'Ta_file', 'ea_file', 'u_file', 'CO2_file', 'z_file', 'tts_file' ... + 'LAI_file', 'hc_file', 'SMC_file', 'Vcmax_file', 'Cab_file', 'LIDF_file'}); + for i = 1:length(F) + k = find(strcmp(F(i).FileID, strtok(X(:, 1)))); + if ~isempty(k) + F(i).FileName = strtok(X(k, 2)); + end + end + + %% 4. input data + [N, X] = xlsread([path_input char(parameter_file)], 'inputdata', ''); + X = X(9:end, 1); + + % Create a structure holding Scope parameters + useXLSX = 1; % set it to 1 or 0, the current stemmus-scope does not support 0 + [ScopeParameters, options] = parameters.loadParameters(options, useXLSX, X, F, N); + + % Define the location information + ScopeParameters.LAT = SiteProperties.latitude; % latitude + ScopeParameters.LON = SiteProperties.longitude; % longitude + ScopeParameters.BSMlat = SiteProperties.latitude; % latitude of BSM model + ScopeParameters.BSMlon = SiteProperties.longitude; % longitude of BSM model + ScopeParameters.z = SiteProperties.reference_height; % reference height + ScopeParameters.hc = SiteProperties.canopy_height; % canopy height + ScopeParameters.Tyear = mean(ForcingData.Ta_msr); % calculate mean air temperature + + % calculate the time zone based on longitude + ScopeParameters.timezn = helpers.calculateTimeZone(SiteProperties.longitude); + + % Input T parameters for different vegetation type + [ScopeParameters] = parameters.setTempParameters(ScopeParameters, SiteProperties.sitename, SiteProperties.landcoverClass); + + %% 5. Declare paths + path_input = InputPath; % path of all inputs + path_output = OutputPath; + %% 6. Numerical parameters (iteration stops etc) + iter.maxit = 100; % maximum number of iterations + iter.maxEBer = 5; % [W m-2] maximum accepted error in energy bal. + iter.Wc = 1; % Weight coefficient for iterative calculation of Tc + + %% 7. Load spectral data for leaf and soil + load([path_input, 'fluspect_parameters/', char(F(3).FileName)]); + rsfile = load([path_input, 'soil_spectrum/', char(F(2).FileName)]); % file with soil reflectance spectra + + %% 8. Load directional data from a file + directional = struct; + if options.calc_directional + anglesfile = load([path_input, 'directional/brdf_angles2.dat']); % Multiple observation angles in case of BRDF calculation + directional.tto = anglesfile(:, 1); % [deg] Observation zenith Angles for calcbrdf + directional.psi = anglesfile(:, 2); % [deg] Observation zenith Angles for calcbrdf + directional.noa = length(directional.tto); % Number of Observation Angles + end + + %% 9. Define canopy structure + canopy.nlayers = 30; + nl = canopy.nlayers; + canopy.x = (-1 / nl:-1 / nl:-1)'; % a column vector + canopy.xl = [0; canopy.x]; % add top level + canopy.nlincl = 13; + canopy.nlazi = 36; + canopy.litab = [5:10:75 81:2:89]'; % a column, never change the angles unless 'ladgen' is also adapted + canopy.lazitab = (5:10:355); % a row + + %% 10. Define spectral regions + [spectral] = io.define_bands(); + + wlS = spectral.wlS; % SCOPE 1.40 definition + wlP = spectral.wlP; % PROSPECT (fluspect) range + wlT = spectral.wlT; % Thermal range + wlF = spectral.wlF; % Fluorescence range + + I01 = find(wlS < min(wlF)); % zero-fill ranges for fluorescence + I02 = find(wlS > max(wlF)); + N01 = length(I01); + N02 = length(I02); + + nwlP = length(wlP); + nwlT = length(wlT); + + nwlS = length(wlS); + + spectral.IwlP = 1:nwlP; + spectral.IwlT = nwlP + 1:nwlP + nwlT; + spectral.IwlF = (640:850) - 399; + + [rho, tau, rs] = deal(zeros(nwlP + nwlT, 1)); + + %% 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); + state_vars = struct(); + end + + %% 12. load time series data + ScopeParametersNames = fieldnames(ScopeParameters); + if options.simulation == 1 + vi = ones(length(ScopeParametersNames), 1); + [soil, leafbio, canopy, meteo, angles, xyt] = io.select_input(ScopeParameters, InitialValues.Theta_LL, vi, canopy, options, SiteProperties, SoilProperties); + [ScopeParameters, xyt, canopy] = io.loadTimeSeries(ScopeParameters, leafbio, soil, canopy, meteo, F, xyt, path_input, options, SiteProperties.landcoverClass); + else + soil = struct; + end + + %% 13. preparations + if options.simulation == 1 + diff_tmin = abs(xyt.t - xyt.startDOY); + diff_tmax = abs(xyt.t - xyt.endDOY); + I_tmin = find(min(diff_tmin) == diff_tmin); + I_tmax = find(min(diff_tmax) == diff_tmax); + if options.soil_heat_method < 2 + meteo.Ta = ForcingData.Ta_msr(1); + soil.Tsold = meteo.Ta * ones(12, 2); + end + end + ScopeParametersNames = fieldnames(ScopeParameters); + nvars = length(ScopeParametersNames); + vmax = ones(nvars, 1); + for i = 1:nvars + name = ScopeParametersNames{i}; + vmax(i) = length(ScopeParameters.(name)); + end + vmax(27, 1) = 1; % these are Tparam and LIDFb + vi = ones(nvars, 1); + switch options.simulation + case 0 + telmax = max(vmax); + [xyt.t, xyt.year] = deal(zeros(telmax, 1)); + case 1 + telmax = size(xyt.t, 1); + case 2 + telmax = prod(double(vmax)); + [xyt.t, xyt.year] = deal(zeros(telmax, 1)); + end + [rad, thermal, fluxes] = io.initialize_output_structures(spectral); + atmfile = [path_input 'radiationdata/' char(F(4).FileName(1))]; + atmo.M = helpers.aggreg(atmfile, spectral.SCOPEspec); + + %% 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. + % Define soil variables for StartInit + VanGenuchten = init.setVanGenuchtenParameters(SoilProperties); + SoilVariables = init.defineSoilVariables(InitialValues, SoilProperties, VanGenuchten); + + % Add initial soil moisture and soil temperature + [SoilInitialValues, BtmX, BtmT, Tss] = io.loadSoilInitialValues(InputPath, TimeProperties, SoilProperties, ForcingData, ModelSettings); + SoilVariables.InitialValues = SoilInitialValues; + SoilVariables.BtmX = BtmX; + SoilVariables.BtmT = BtmT; + SoilVariables.Tss = Tss; + + % Run StartInit + [SoilVariables, VanGenuchten, ThermalConductivity] = StartInit(SoilVariables, SoilProperties, VanGenuchten, ModelSettings); + + % Set SoilVariables that are used in the loop + T = SoilVariables.T; + h = SoilVariables.h; + TT = SoilVariables.TT; + h_frez = SoilVariables.h_frez; + hh = SoilVariables.hh; + hh_frez = SoilVariables.hh_frez; + + % get soil constants + SoilConstants = io.getSoilConstants(); + + %% The boundary condition information settings + BoundaryCondition = init.setBoundaryCondition(SoilVariables, ForcingData, SiteProperties.landcoverClass(1), ModelSettings); + DSTOR = BoundaryCondition.DSTOR; + DSTOR0 = BoundaryCondition.DSTOR0; + RS = BoundaryCondition.RS; + DSTMAX = BoundaryCondition.DSTMAX; + IRPT1 = BoundaryCondition.IRPT1; + IRPT2 = BoundaryCondition.IRPT2; + + %% 14. Initialize matrices + calculate = 1; + TIMEOLD = 0; + SAVEhh_frez = zeros(ModelSettings.NN + 1, 1); + FCHK = zeros(1, ModelSettings.NN); + KCHK = zeros(1, ModelSettings.NN); + hCHK = zeros(1, ModelSettings.NN); + TIMELAST = 0; + + % the start of simulation period is from 0mins, while the input data start from 30mins. + TIME = 0; % Time of simulation released; + Delt_t = TimeProperties.DELT; % Duration of time step [Unit of second] + KT = ModelSettings.KT; + kk = 0; + TimeStep = []; + TEND = TIME + TimeProperties.DELT * TimeProperties.Dur_tot; % Time to be reached at the end of simulation period + Delt_t0 = Delt_t; % Duration of last time step + TOLD_CRIT = []; + + % for soil moisture and temperature outputs + monitorDepthTemperature = ModelSettings.NL:-1:1; + monitorDepthSoilMoisture = ModelSettings.NL:-1:1; + Sim_Theta_U = []; + Sim_Temp = []; + + % for matric potential and fluxes outputs, added by Mostafa + Sim_hh = []; + Sim_qlh = []; + Sim_qlt = []; + Sim_qvh = []; + Sim_qvt = []; + Sim_qla = []; + Sim_qva = []; + + % Srt, root water uptake; + Srt = InitialValues.Srt; % will be updated! + P_gg = InitialValues.P_gg; % will be updated! + P_g = InitialValues.P_g; % will be updated! + + %% soil layer information + SoilLayer.thickness = ModelSettings.DeltZ_R; + SoilLayer.depth = Ztot'; % used in Initial_root_biomass + + % NOTE: workspace will be saved, this code block is after the update step. + % (this is to not repeat the save-workspace code). + disp('Finished model initialization'); +end + +% If the runMode is update, retrieve the (possibly) updated state. +% The state can be modified by the STEMMUS_SCOPE BMI. See PyStemmusScope. +if strcmp(bmiMode, 'update') + load([OutputPath, 'STEMMUS_SCOPE_state.mat']); % Load the workspace to be able to (continue) running the model + + if KT + 1 >= TimeProperties.Dur_tot + bmiMode = 'none'; % Ensure the model does not try to update. + disp("Finished running the model. Updating won't do anything!"); + else + endTime = KT + 1; + end +elseif strcmp(runMode, 'full') + endTime = TimeProperties.Dur_tot; + disp('The calculations start now'); +end + +% Actually run the model +if strcmp(bmiMode, 'update') || strcmp(runMode, 'full') + + if GroundwaterSettings.GroundwaterCoupling == 1 % Groundwater coupling is enabled + BoundaryCondition.NBChB = 1; + + % Update GroundwaterSettings.headBotmLayer and GroundwaterSettings.tempBotm, from MODFLOW through BMI + GroundwaterSettings.gw_Dep = groundwater.calculateGroundWaterDepth(GroundwaterSettings.topLevel, GroundwaterSettings.headBotmLayer, ModelSettings.Tot_Depth); + + % Update Dunnian runoff and ForcingData.Precip_msr + [ForcingData.R_Dunn, ForcingData.Precip_msr] = groundwater.updateDunnianRunoff(ForcingData.Precip_msr, GroundwaterSettings.gw_Dep); + + % Calculate the index of the bottom layer level + [GroundwaterSettings.indxBotmLayer, GroundwaterSettings.indxBotmLayer_R] = groundwater.calculateIndexBottomLayer(GroundwaterSettings.soilThick, GroundwaterSettings.gw_Dep, ModelSettings); + + % Check the position of the groundwater table + if any(GroundwaterSettings.gw_Dep > ModelSettings.Tot_Depth) || any(isempty(GroundwaterSettings.indxBotmLayer)) + warning('Groundwater table is below the end of the soil column. Please enlarge the total soil thickness!'); + return + end + + % Check if no available groundwater temperature data + if isnan(GroundwaterSettings.tempBotm) + GroundwaterSettings.tempBotm = SoilVariables.TT(GroundwaterSettings.indxBotmLayer); + end + end + + % 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 + TIME = TIME + Delt_t; % The time elapsed since start of simulation + TimeStep(KT, 1) = Delt_t; + SUMTIME(KT, 1) = TIME; + Processing = TIME / TEND; + NoTime(KT) = fix(SUMTIME(KT) / TimeProperties.DELT); + if NoTime(KT) == 0 + k = 1; + else + k = NoTime(KT); + end + %%%%% Updating the state variables. %%%%%%%%%%%%%%%%%%%%%%%%%%%% + L_f = 0; % ignore Freeze/Thaw, see issue 139 + TT_CRIT(ModelSettings.NN) = ModelSettings.T0; % unit K + hOLD_frez = []; + if IRPT1 == 0 && IRPT2 == 0 && SoilVariables.ISFT == 0 + for MN = 1:ModelSettings.NN + hOLD_frez(MN) = h_frez(MN); + h_frez(MN) = hh_frez(MN); + TOLD_CRIT(MN) = T_CRIT(MN); + T_CRIT(MN) = TT_CRIT(MN); + + hOLD(MN) = h(MN); + h(MN) = hh(MN); + DDhDZ(MN, KT) = InitialValues.DhDZ(MN); + if ModelSettings.Thmrlefc == 1 + TOLD(MN) = T(MN); + T(MN) = TT(MN); + TTT(MN, KT) = TT(MN); + end + if ModelSettings.Soilairefc == 1 + P_gOLD(MN) = P_g(MN); + P_g(MN) = P_gg(MN); + end + if ModelSettings.rwuef == 1 + SRT(MN, KT) = Srt(MN, 1); + end + end + if options.simulation == 1 + vi(vmax > 1) = k; + end + 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 + if WofostPar.PHOpt == 0 % + ScopeParameters.hc(KT) = WofostPar.PHEM; % initial PH + end + elseif KT > WofostPar.CSTART + ScopeParameters.LAI(KT) = crop_output(KT - 1, 3); % substitute LAI + if WofostPar.PHOpt == 0 + 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 + fprintf('simulation %i ', k); + fprintf('of %i \n', telmax); + else + calculate = 1; + end + + if calculate + iter.counter = 0; + LIDF_file = char(F(22).FileName); + if ~isempty(LIDF_file) + canopy.lidf = dlmread([path_input, 'leafangles/', LIDF_file], '', 3, 0); + else + canopy.lidf = equations.leafangles(canopy.LIDFa, canopy.LIDFb); % This is 'ladgen' in the original SAIL model, + end + + leafbio.emis = 1 - leafbio.rho_thermal - leafbio.tau_thermal; + + if options.calc_PSI + fversion = @fluspect_B_CX; + else + fversion = @fluspect_B_CX_PSI_PSII_combined; + end + leafbio.V2Z = 0; + leafopt = fversion(spectral, leafbio, optipar); + leafbio.V2Z = 1; + leafoptZ = fversion(spectral, leafbio, optipar); + IwlP = spectral.IwlP; + IwlT = spectral.IwlT; + rho(IwlP) = leafopt.refl; + tau(IwlP) = leafopt.tran; + rlast = rho(nwlP); + tlast = tau(nwlP); + + if options.soilspectrum == 0 + rs(IwlP) = rsfile(:, soil.spectrum + 1); + else + soilemp.SMC = 25; % empirical parameter (fixed) + soilemp.film = 0.015; % empirical parameter (fixed) + rs(IwlP) = BSM(soil, optipar, soilemp); + end + rslast = rs(nwlP); + + switch options.rt_thermal + case 0 + rho(IwlT) = ones(nwlT, 1) * leafbio.rho_thermal; + tau(IwlT) = ones(nwlT, 1) * leafbio.tau_thermal; + rs(IwlT) = ones(nwlT, 1) * soil.rs_thermal; + case 1 + rho(IwlT) = ones(nwlT, 1) * rlast; + tau(IwlT) = ones(nwlT, 1) * tlast; + rs(IwlT) = ones(nwlT, 1) * rslast; + end + leafopt.refl = rho; % extended wavelength ranges are stored in structures + leafopt.tran = tau; + reflZ = leafopt.refl; + tranZ = leafopt.tran; + reflZ(1:300) = leafoptZ.refl(1:300); + tranZ(1:300) = leafoptZ.tran(1:300); + leafopt.reflZ = reflZ; + leafopt.tranZ = tranZ; + soil.refl = rs; + soil.Ts = meteo.Ta * ones(2, 1); % initial soil surface temperature + + if length(F(4).FileName) > 1 && options.simulation == 0 + atmfile = [path_input 'radiationdata/' char(F(4).FileName(k))]; + atmo.M = helpers.aggreg(atmfile, spectral.SCOPEspec); + end + atmo.Ta = meteo.Ta; + + [rad, gap, profiles] = RTMo(spectral, atmo, soil, leafopt, canopy, angles, meteo, rad, options); + + switch options.calc_ebal + case 1 + [iter, fluxes, rad, thermal, profiles, soil, RWU, frac, WaterStressFactor, WaterPotential] ... + = ebal(iter, options, spectral, rad, gap, ... + leafopt, angles, meteo, soil, canopy, leafbio, xyt, k, profiles, Delt_t, ... + Rl, SoilVariables, VanGenuchten, InitialValues, ModelSettings, GroundwaterSettings); + if options.calc_fluor + if options.calc_vert_profiles + [rad, profiles] = RTMf(spectral, rad, soil, leafopt, canopy, gap, angles, profiles); + else + [rad] = RTMf(spectral, rad, soil, leafopt, canopy, gap, angles, profiles); + end + end + if options.calc_xanthophyllabs + [rad] = RTMz(spectral, rad, soil, leafopt, canopy, gap, angles, profiles); + end + + if options.calc_planck + rad = RTMt_planck(spectral, rad, soil, leafopt, canopy, gap, angles, thermal.Tcu, thermal.Tch, thermal.Ts(2), thermal.Ts(1), 1); + end + + if options.calc_directional + directional = calc_brdf(options, directional, spectral, angles, rad, atmo, soil, leafopt, canopy, meteo, profiles, thermal); + end + + otherwise + Fc = (1 - gap.Ps(1:end - 1))' / nl; % Matrix containing values for Ps of canopy + fluxes.aPAR = canopy.LAI * (Fc * rad.Pnh + equations.meanleaf(canopy, rad.Pnu, 'angles_and_layers', gap.Ps)); % net PAR leaves + fluxes.aPAR_Cab = canopy.LAI * (Fc * rad.Pnh_Cab + equations.meanleaf(canopy, rad.Pnu_Cab, 'angles_and_layers', gap.Ps)); % net PAR leaves + [fluxes.aPAR_Wm2, fluxes.aPAR_Cab_eta] = deal(canopy.LAI * (Fc * rad.Rnh_PAR + equations.meanleaf(canopy, rad.Rnu_PAR, 'angles_and_layers', gap.Ps))); % net PAR leaves + if options.calc_fluor + profiles.etah = ones(nl, 1); + profiles.etau = ones(13, 36, nl); + if options.calc_vert_profiles + [rad, profiles] = RTMf(spectral, rad, soil, leafopt, canopy, gap, angles, profiles); + else + [rad] = RTMf(spectral, rad, soil, leafopt, canopy, gap, angles, profiles); + end + end + end + + if options.calc_fluor % total emitted fluorescence irradiance (excluding leaf and canopy re-absorption and scattering) + if options.calc_PSI + rad.Femtot = 1E3 * (leafbio.fqe(2) * optipar.phiII(spectral.IwlF) * fluxes.aPAR_Cab_eta + leafbio.fqe(1) * optipar.phiI(spectral.IwlF) * fluxes.aPAR_Cab); + else + rad.Femtot = 1E3 * leafbio.fqe * optipar.phi(spectral.IwlF) * fluxes.aPAR_Cab_eta; + 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, state_vars] = wofost.cropgrowth(crop_output, state_vars, meteo, WofostPar, Anet, WaterStressFactor, xyt, KT, TimeProperties.Dur_tot); + end + + if options.simulation == 2 && telmax > 1 + vi = helpers.count(nvars, vi, vmax, 1); + end + if KT == 1 + if isreal(fluxes.Actot) && isreal(thermal.Tsave) && isreal(fluxes.lEstot) && isreal(fluxes.lEctot) + Acc = fluxes.Actot; + Tss = thermal.Tsave; + else + Acc = 0; + fluxes.lEstot = 0; + fluxes.lEctot = 0; + Tss = ForcingData.Ta_msr(KT); + end + elseif NoTime(KT) > NoTime(KT - 1) + if isreal(fluxes.Actot) && isreal(thermal.Tsave) && isreal(fluxes.lEstot) && isreal(fluxes.lEctot) + Acc = fluxes.Actot; + Tss = thermal.Tsave; + else + Acc = 0; + fluxes.lEstot = 0; + fluxes.lEctot = 0; + Tss = ForcingData.Ta_msr(KT); + end + end + + DSTOR0 = DSTOR; + BoundaryCondition.DSTOR0 = DSTOR0; + + if KT > 1 + SoilVariables.XWRE = updateWettingHistory(SoilVariables, VanGenuchten, ModelSettings); + end + + for i = 1:ModelSettings.NL + DDhDZ(i, KT) = InitialValues.DhDZ(i); + end + end + if Delt_t ~= Delt_t0 + for MN = 1:ModelSettings.NN + hh(MN) = h(MN) + (h(MN) - hOLD(MN)) * Delt_t / Delt_t0; + TT(MN) = T(MN) + (T(MN) - TOLD(MN)) * Delt_t / Delt_t0; + end + end + hSAVE = hh(ModelSettings.NN); + TSAVE = TT(ModelSettings.NN); + + % set "hN" empty when the "if statement" below does not happen, see issue 98, + % item 5 + hN = []; + if BoundaryCondition.NBCh == 1 + hN = BoundaryCondition.BCh; + hh(ModelSettings.NN) = hN; + hSAVE = hN; + elseif BoundaryCondition.NBCh == 2 + if BoundaryCondition.NBChh ~= 2 + if BoundaryCondition.BCh < 0 + hN = DSTOR0; + hh(ModelSettings.NN) = hN; + hSAVE = hN; + else + hN = -1e6; + hh(ModelSettings.NN) = hN; + hSAVE = hN; + end + end + else + if BoundaryCondition.NBChh ~= 2 + hN = DSTOR0; + hh(ModelSettings.NN) = hN; + hSAVE = hN; + end + end + + % Start the iteration procedure in a time step. + SoilVariables.h = h; + SoilVariables.hh = hh; + SoilVariables.TT = TT; + SoilVariables.T = T; + SoilVariables.h_frez = h_frez; + SoilVariables.Tss(KT) = Tss; + + for KIT = 1:ModelSettings.NIT + [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, ModelSettings.NN, SoilConstants.hd, ForcingData.Tmin); + % update inputs for UpdateSoilWaterContent + SoilVariables.TT_CRIT = TT_CRIT; + SoilVariables.hh_frez = hh_frez; + + SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten, ModelSettings); + + % Reset KL_T here. CondL_T script is replaced by this line + % see issue 181, item 4 + KL_T = InitialValues.KL_T; + + [RHOV, DRHOVh, DRHOVT] = Density_V(SoilVariables.TT, SoilVariables.hh, Constants.g, Constants.Rv, ModelSettings.NN); + + TransportCoefficient = conductivity.calculateTransportCoefficient(InitialValues, SoilVariables, VanGenuchten, Delt_t, ModelSettings); + W = TransportCoefficient.W; + D_Ta = TransportCoefficient.D_Ta; + + [L] = Latent(SoilVariables.TT, ModelSettings.NN); + % DRHODAt unused! + [Xaa, XaT, Xah, DRHODAt, DRHODAz, RHODA] = Density_DA(SoilVariables.T, Constants.RDA, P_g, Constants.Rv, ModelSettings.DeltZ, SoilVariables.h, SoilVariables.hh, SoilVariables.TT, P_gg, Delt_t, ModelSettings.NL, ModelSettings.NN, DRHOVT, DRHOVh, RHOV); + + ThermalConductivityCapacity = conductivity.calculateThermalConductivityCapacity(InitialValues, ThermalConductivity, SoilVariables, VanGenuchten, ModelSettings, DRHOVT, L, RHOV); + + k_g = conductivity.calculateGasConductivity(InitialValues, TransportCoefficient, VanGenuchten, SoilVariables, ModelSettings); + + VaporVariables = conductivity.calculateVaporVariables(InitialValues, SoilVariables, VanGenuchten, ModelSettings, ThermalConductivityCapacity, SoilVariables.TT); + + GasDispersivity = conductivity.calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g, ModelSettings); + + % Srt is both input and output + [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip, RWUs, RWUg, ForcingData] = ... + soilmoisture.solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, ... + TimeProperties, SoilProperties, BoundaryCondition, Delt_t, RHOV, DRHOVh, ... + DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, ModelSettings, GroundwaterSettings); + + if BoundaryCondition.NBCh == 1 + DSTOR = 0; + RS = 0; + elseif BoundaryCondition.NBCh == 2 + AVAIL = -BoundaryCondition.BCh; + EXCESS = (AVAIL + HBoundaryFlux.QMT) * Delt_t; + if abs(EXCESS / Delt_t) <= 1e-10 + EXCESS = 0; + end + DSTOR = min(EXCESS, DSTMAX); + RS = (EXCESS - DSTOR) / Delt_t; + else + AVAIL = AVAIL0 - Evap; + EXCESS = (AVAIL + HBoundaryFlux.QMT) * Delt_t; % (unit converstion from cm/sec to cm/30mins) + if abs(EXCESS / Delt_t) <= 1e-10 + EXCESS = 0; + end + DSTOR = min(EXCESS, DSTMAX); % Depth of depression storage at end of current time step + % Next line is commented and Surface runoff is re-calculated using different approach (the following 3 lines) + % RS(KT) = (EXCESS - DSTOR) / Delt_t; % surface runoff, (unit conversion from cm/30mins to cm/sec) + RS(KT) = ForcingData.R_Hort(KT) + ForcingData.R_Dunn(KT); % total surface runoff (cm/sec) + end + + if ModelSettings.Soilairefc == 1 + [AirVariabes, RHS, SAVE, P_gg] = dryair.solveDryAirEquations(SoilVariables, GasDispersivity, TransportCoefficient, InitialValues, VaporVariables, ... + BoundaryCondition, ForcingData, P_gg, P_g, Xah, XaT, Xaa, RHODA, KT, Delt_t, ModelSettings, GroundwaterSettings); + else + AirVariabes.KLhBAR = InitialValues.KLhBAR; + AirVariabes.KLTBAR = InitialValues.KLTBAR; + AirVariabes.DDhDZ = DDhDZ; % DDhDZ is not defined as InitialValues, see issue 100, item 3 + AirVariabes.DhDZ = InitialValues.DhDZ; + AirVariabes.DTDZ = InitialValues.DTDZ; + AirVariabes.Kaa = InitialValues.Kaa; + AirVariabes.Vaa = InitialValues.Vaa; + AirVariabes.QL = InitialValues.QL; + end + + if ModelSettings.Thmrlefc == 1 + % CHK will be updated + [RHS, SAVE, CHK, SoilVariables, EnergyVariables] = energy.solveEnergyBalanceEquations(InitialValues, SoilVariables, HeatVariables, TransportCoefficient, ... + AirVariabes, VaporVariables, GasDispersivity, ThermalConductivityCapacity, ... + HBoundaryFlux, BoundaryCondition, ForcingData, DRHOVh, DRHOVT, KL_T, ... + Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, Delt_t, P_g, P_gg, ... + TOLD, Precip, Evap, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings); + end + + if max(CHK) < 0.1 + break + end + hSAVE = SoilVariables.hh(ModelSettings.NN); + TSAVE = SoilVariables.TT(ModelSettings.NN); + end + + TIMEOLD = KT; + KIT; + KIT = 0; + + [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, ModelSettings.NN, SoilConstants.hd, ForcingData.Tmin); + % updates inputs for UpdateSoilWaterContent + SoilVariables.TT_CRIT = TT_CRIT; + SoilVariables.hh_frez = hh_frez; + SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten, ModelSettings); + + if IRPT1 == 0 && IRPT2 == 0 + if KT % In case last time step is not convergent and needs to be repeated. + for i = 1:ModelSettings.NL + hhh(i, KT) = SoilVariables.hh(i); + qlh(i, KT) = transpose(EnergyVariables.QL_h(i)); + qlt(i, KT) = transpose(EnergyVariables.QL_T(i)); + qla(i, KT) = transpose(EnergyVariables.QL_a(i)); + qvh(i, KT) = transpose(EnergyVariables.QVH(i)); + qvt(i, KT) = transpose(EnergyVariables.QVT(i)); + qva(i, KT) = transpose(EnergyVariables.QVa(i)); + qtot(i, KT) = transpose(EnergyVariables.Qtot(i)); + for j = 1:ModelSettings.nD + Theta_LLL(i, j, KT) = SoilVariables.Theta_LL(i, j); + SoilVariables.Theta_L(i, j) = SoilVariables.Theta_LL(i, j); + Theta_UUU(i, j, KT) = SoilVariables.Theta_UU(i, j); + SoilVariables.Theta_U(i, j) = SoilVariables.Theta_UU(i, j); + Theta_III(i, j, KT) = SoilVariables.Theta_II(i, j); + SoilVariables.Theta_I(i, j) = SoilVariables.Theta_II(i, j); + end + end + + % replace run ObservationPoints, see issue 101 + Sim_Theta_U(KT, 1:length(monitorDepthSoilMoisture)) = Theta_UUU(monitorDepthSoilMoisture, 1, KT); + Sim_Temp(KT, 1:length(monitorDepthTemperature)) = TTT(monitorDepthTemperature, KT); + Sim_hh(KT, 1:length(monitorDepthSoilMoisture)) = hhh(monitorDepthSoilMoisture, KT); + Sim_qlh(KT, 1:length(monitorDepthSoilMoisture)) = qlh(monitorDepthSoilMoisture, KT); + Sim_qlt(KT, 1:length(monitorDepthSoilMoisture)) = qlt(monitorDepthSoilMoisture, KT); + Sim_qla(KT, 1:length(monitorDepthSoilMoisture)) = qla(monitorDepthSoilMoisture, KT); + Sim_qvh(KT, 1:length(monitorDepthSoilMoisture)) = qvh(monitorDepthSoilMoisture, KT); + Sim_qvt(KT, 1:length(monitorDepthSoilMoisture)) = qvt(monitorDepthSoilMoisture, KT); + Sim_qva(KT, 1:length(monitorDepthSoilMoisture)) = qva(monitorDepthSoilMoisture, KT); + Sim_qtot(KT, 1:length(monitorDepthSoilMoisture)) = qtot(monitorDepthSoilMoisture, KT); + end + if (TEND - TIME) < 1E-3 + for i = 1:ModelSettings.NN + hOLD(i) = SoilVariables.h(i); + SoilVariables.h(i) = SoilVariables.hh(i); + if ModelSettings.Thmrlefc == 1 + TOLD(i) = SoilVariables.T(i); + SoilVariables.T(i) = SoilVariables.TT(i); + TTT(i, KT) = SoilVariables.TT(i); + TOLD_CRIT(i) = T_CRIT(i); + T_CRIT(i) = TT_CRIT(i); + hOLD_frez(i) = SoilVariables.h_frez(i); + SoilVariables.h_frez(i) = SoilVariables.hh_frez(i); + end + if ModelSettings.Soilairefc == 1 + P_gOLD(i) = P_g(i); + P_g(i) = P_gg(i); + end + end + end + end + + % Recharge calculations, added by Mostafa + if GroundwaterSettings.GroundwaterCoupling == 1 % Groundwater coupling is enabled + gwfluxes = groundwater.calculateGroundwaterRecharge(EnergyVariables, SoilVariables, KT, ModelSettings, GroundwaterSettings); + if GroundwaterSettings.gw_Dep <= 1 % soil is fully saturated + gwfluxes.recharge = 0; + end + else + gwfluxes.recharge = 0; + gwfluxes.indxRchrg = NaN; + end + + % set SoilVariables for the rest of the loop + h = SoilVariables.h; + hh = SoilVariables.hh; + T = SoilVariables.T; + TT = SoilVariables.TT; + hh_frez = SoilVariables.hh_frez; + h_frez = SoilVariables.h_frez; + kk = k; + % Open files for writing + file_ids = structfun(@(x) fopen(x, 'a'), fnames, 'UniformOutput', false); + 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, crop_output); + fclose("all"); + end +end + +if strcmp(bmiMode, 'initialize') || strcmp(bmiMode, 'update') + % Save the required variables to the model state file. + % NOTE: bmiVarNames are defined in STEMMUS_SCOPE_exe.m + % save & load is incompatible with Octave. Octave only supports full run mode, not BMI mode. + save([Output_dir, 'STEMMUS_SCOPE_state.mat'], bmiVarNames{:}, "-v7.3", "-nocompression"); +end + +if strcmp(bmiMode, 'finalize') + % Load the workspace to be able to finalize the model. + load([OutputPath, 'STEMMUS_SCOPE_state.mat']); +end + +if strcmp(bmiMode, 'finalize') || strcmp(runMode, 'full') + disp('Finalizing STEMMUS_SCOPE'); + if options.verify + io.output_verification(Output_dir); + end + + io.bin_to_csv(fnames, n_col, k, options, SoilLayer, GroundwaterSettings, FullCSVfiles); + save([Output_dir, 'output.mat']); +end + +% Calculate the total simulatiom time, added by mostafa +end_time = clock; +simtime_min = etime(end_time, start_time) / 60; +simtime_hr = simtime_min / 24; +disp(['Simulation time is : ' num2str(simtime_hr) ' hrs (' num2str(simtime_min) ' minutes)']); diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 1db95466..440df004 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -198,7 +198,15 @@ [rho, tau, rs] = deal(zeros(nwlP + nwlT, 1)); - %% 11. load time series data + %% 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); + state_vars = struct(); + end + + %% 12. load time series data ScopeParametersNames = fieldnames(ScopeParameters); if options.simulation == 1 vi = ones(length(ScopeParametersNames), 1); @@ -208,7 +216,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 +250,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. @@ -374,6 +382,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 @@ -420,6 +430,22 @@ 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 + if WofostPar.PHOpt == 1 % Options of PH simulations; 0 using the prescribed values, and 1 for simulation + ScopeParameters.hc(KT) = WofostPar.PHEM; % initial PH + end + elseif KT > WofostPar.CSTART + ScopeParameters.LAI(KT) = crop_output(KT - 1, 3); % substitute LAI + if WofostPar.PHOpt == 1 + 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 fprintf('simulation %i ', k); @@ -542,6 +568,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, state_vars] = wofost.cropgrowth(crop_output, state_vars, meteo, WofostPar, Anet, WaterStressFactor, xyt, KT, TimeProperties.Dur_tot); + end + if options.simulation == 2 && telmax > 1 vi = helpers.count(nvars, vi, vmax, 1); end @@ -799,7 +836,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