upload the original code of PG model

config_file_crib.txt

@@ -1,15 +1,15 @@
-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/
 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=
+StartTime=2001-04-22T00:00
+EndTime=2001-04-23T00:00
+InputPath=../../input/US-Bo1_2024-09-18-1011/
+OutputPath=../../output/

example_data/plant_growth/CropD.crp

@@ -0,0 +1,183 @@
+***********************************************************************************************
+* Filename: Maize.CRP
+* Contents: SWAP 3.2 - Data for detailed crop model
+***********************************************************************************************
+*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     =   950.00    ! Temperature sum from emergence to anthesis, [0..10000 C, R]
+  TSUMAM     =   880.00     ! Temperature sum from anthesis to maturity  [0..10000 C, R]
+
+* List increase in temperature sum [0..60 C, R] as function of daily average temp. [0..100 C, R]
+*         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.3 ! 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   =   35.00 ! Life span under leaves under optimum conditions, [0..366 d, R]
+  LAICR   =  6.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 [0..1 ha/kg, R] as function of devel. stage [0..2, R]
+
+*         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
+
+  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 [0..1 m, R] as function of devel. stage [0..2, R]
+
+*         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 [kg/kg, R]
+* as function of development stage [0..2 -, R]
+*          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 [kg/kg, R]
+* as function of development stage [0..2 -, R]
+
+*          DVS     FL   (maximum 15 records)
+  FLTB = 
+           0.0000   0.6000
+           0.7000   0.5000
+           0.7500   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 [kg/kg, R]
+* as function of development stage [0..2 -, R]
+
+*          DVS    FS   (maximum 15 records)
+  FSTB = 
+           0.0000   0.4000
+           0.7000   0.5000
+           0.7500   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 [kg/kg, R]
+* as function of development stage [0..2 -, R]
+
+*          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 [kg/kg/d] as function of dev. stage [0..2 -, R]
+*          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 [kg/kg/d] as function of dev. stage [0..2 -, R]
+*          DVS     RDRS    (maximum 15 records)
+  RDRSTB = 
+          0.0000 0.0000
+          0.9000 0.0020
+          1.2500 0.0020
+          2.0000 0.0030
+* End of table
+***********************************************************************************************
+
+***********************************************************************************************
+* Part 8: Root density distribution and root growth               
+*
+* List relative root density [0..1 -, R], as function of rel. rooting depth [0..1 -, R]:
+*          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 !

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

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');
 

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);

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)

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);

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
+