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/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,

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

@@ -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);
@@ -223,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

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

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)