Integrate PHS into STEMMUS-SCOPE

src/+equations/Soil_Inertia1.m

@@ -7,18 +7,20 @@
     Sr = SMC / theta_s;
 
     % fss = 0.58; %(sand fraction)
-    gamma_s = 0.27; % (soil texture dependent parameter)
+    gamma_s = 0.96; % (soil texture dependent parameter)
     dels = 1.33; % (shape parameter)
 
     ke = exp(gamma_s * (1 - power(Sr, gamma_s - dels)));
 
     phis  = theta_s0; % (phis == theta_s)
     lambda_d = -0.56 * phis + 0.51;
 
-    QC = 0.20; % (quartz content)
+    QC = 0.60; % (quartz content)
     lambda_qc = 7.7;  % (thermal conductivity of quartz, constant)
 
-    lambda_s = (lambda_qc^(QC)) * lambda_d^(1 - QC);
+    lambda_o = 2.0;
+
+    lambda_s = (lambda_qc^(QC)) * lambda_o^(1 - QC);
     lambda_wtr = 0.57;   % (thermal conductivity of water, W/m.K, constant)
 
     lambda_w = (lambda_s^(1 - phis)) * lambda_wtr^(phis);

src/+init/setBoundaryCondition.m

@@ -1,4 +1,5 @@
-function BoundaryCondition = setBoundaryCondition(SoilVariables, ForcingData, initialLandcoverClass, ModelSettings)
+function BoundaryCondition = setBoundaryCondition(SoilVariables, ForcingData, initialLandcoverClass, ModelSettings, siteName)
+
     % Variables information for boundary condition settings
     % NBCh   Indicator for type of surface boundary condition on mass euqation to be applied;
     %        "1"--Specified matric head;
@@ -66,6 +67,9 @@
     else
         NBChB = 3;
     end
+    if strcmp(siteName,  "CH-HTC")
+        NBChB = 1;
+    end
     BChB = -9e-10;
     if ModelSettings.Thmrlefc == 1
         % NBCT: Energy Surface B.C.:

src/+io/getModelSettings.m

@@ -2,7 +2,7 @@
     %{
     %}
 
-    ModelSettings.R_depth = 350;
+    ModelSettings.R_depth = 30;
 
     % Indicator denotes the index of soil type for choosing soil physical parameters
     ModelSettings.J = 1;
@@ -51,7 +51,7 @@
     ModelSettings.rroot = 1.5 * 1e-3;
     ModelSettings.SFCC = 1;
 
-    ModelSettings.Tot_Depth = 500; % Unit is cm. it should be usually bigger than 0.5m. Otherwise,
+    ModelSettings.Tot_Depth = 100; % Unit is cm. it should be usually bigger than 0.5m. Otherwise,
     ModelSettings.Eqlspace = 0; % Indicator for deciding is the space step equal or not; % the DeltZ would be reset in 50cm by hand;
 
     ModelSettings.NS = 1; % Number of soil types;
@@ -61,7 +61,7 @@
     ModelSettings.KT = 0; % Number of time steps;
 
     % Determination of NL, the number of elments
-    ModelSettings.NL = 100;
+    ModelSettings.NL = 29;
     if ~ModelSettings.Eqlspace
         [DeltZ, DeltZ_R, NL, ML] = Dtrmn_Z(ModelSettings.NL, ModelSettings.Tot_Depth);
     else

src/+io/loadSoilInitialValues.m

@@ -34,12 +34,18 @@
     BtmX = BtmX;
     Tss = Tss;
 
-    InitND1 = 5;    % Unit of it is cm. These variables are used to indicated the depth corresponding to the measurement.
-    InitND2 = 15;
-    InitND3 = 60;
-    InitND4 = 100;
-    InitND5 = 200;
-    InitND6 = 300;
+    % This input of soil layer thickness and initial depth will be revised
+    % later via reading csv file. I read the initial soil depth variable
+    % from .mat file.
+    % ZSong (August 2024)
+    if ~exist("InitND1", "var")
+        InitND1 = 5;    % Unit of it is cm. These variables are used to indicated the depth corresponding to the measurement.
+        InitND2 = 15;
+        InitND3 = 60;
+        InitND4 = 100;
+        InitND5 = 200;
+        InitND6 = 300;
+    end
     if SWCC == 0
         InitT0 = -1.762;  % -1.75estimated soil surface temperature-1.762
         InitT1 = -0.662;
@@ -71,7 +77,8 @@
         if nanmean(Ta_msr) < 0
             BtmT = 0;  % 9 8.1
         else
-            BtmT = nanmean(Ta_msr);
+           BtmT = 14.5; 
+            % BtmT = nanmean(Ta_msr);
         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/loadTimeSeries.m

@@ -122,7 +122,7 @@
         ScopeParameters.SMC = load(fullfile(path_input, SMC_file));
     end
 
-    %% 7. Set Leaf Vcmo, Tparam, m, Type, Rdparam, and leafwidth as time-dependent
+    %% 7. Set Leaf Vcmo, Tparam, m, Type, Rdparam,leafwidth, g1Med and g0Med as time-dependent
     %    parameters. The timeseries of landcover, along with the lookup table values
     %    (lc...) are used to generate the timeseries.
     landcovers = unique(landcoverClass, 'stable');
@@ -134,9 +134,11 @@
         ScopeParameters.Type(timeIndex, 1) = ScopeParameters.lcType(landcoverIndex, 1);
         ScopeParameters.Rdparam(timeIndex, 1) = ScopeParameters.lcRdparam(landcoverIndex, 1);
         ScopeParameters.leafwidth(timeIndex, 1) = ScopeParameters.lcleafwidth(landcoverIndex, 1);
+        % ScopeParameters.g1Med(timeIndex, 1) = ScopeParameters.lcg1Med(landcoverIndex, 1);
+        % ScopeParameters.g0Med(timeIndex, 1) = ScopeParameters.lcg0Med(landcoverIndex, 1);
     end
     % Remove the intermediate variables that were required to generate the timeseries:
-    ScopeParameters = rmfield(ScopeParameters, {'lcVcmo', 'lcTparam', 'lcm', 'lcType', 'lcRdparam', 'lcleafwidth'});
+    ScopeParameters = rmfield(ScopeParameters, {'lcVcmo', 'lcTparam', 'lcm', 'lcType', 'lcRdparam', 'lcleafwidth', 'lcg1Med', 'lcg0Med'});
 
     if ~isempty(Cab_file)
         Cabtable = load(fullfile(path_input, Cab_file));

src/+io/select_input.m

@@ -30,6 +30,8 @@
     leafbio.Tparam  = ScopeParameters.Tparam(digitsVector(14), :); % this is correct (: instead of 14)
     fqe             = ScopeParameters.fqe(digitsVector(15));
     leafbio.Rdparam = ScopeParameters.Rdparam(digitsVector(13));
+    leafbio.g1Med   = ScopeParameters.g1Med(digitsVector(65));
+    leafbio.g0Med   = ScopeParameters.g0Med(digitsVector(66));
 
     leafbio.rho_thermal = ScopeParameters.rho_thermal(digitsVector(7));
     leafbio.tau_thermal = ScopeParameters.tau_thermal(digitsVector(8));
@@ -39,6 +41,8 @@
     leafbio.kNPQs         = ScopeParameters.kNPQs(digitsVector(57));
     leafbio.qLs           = ScopeParameters.qLs(digitsVector(58));
     leafbio.stressfactor  = ScopeParameters.stressfactor(digitsVector(59));
+    leafbio.g1Med         = ScopeParameters.g1Med;
+    leafbio.g0Med         = ScopeParameters.g0Med;
 
     canopy.LAI  = ScopeParameters.LAI(digitsVector(22));
     canopy.hc  = ScopeParameters.hc(digitsVector(23));

src/+parameters/loadParameters.m

@@ -12,7 +12,7 @@
                             'Ca', 'Oa', 'rb', 'Cd', 'CR', 'CD1', 'Psicor', 'CSSOIL', 'rbs', 'rwc', ...
                             'startDOY', 'endDOY', 'LAT', 'LON', 'timezn', 'tts', 'tto', 'psi', 'SMC', ...
                             'Tyear', 'beta', 'kNPQs', 'qLs', 'stressfactor', 'Cant', 'BSMBrightness', ...
-                            'BSMlat', 'BSMlon', 'BallBerry0'
+                            'BSMlat', 'BSMlon', 'BallBerry0', 'g1Med', 'g0Med'
                            };
 
     % create an empty structure with field names of ScopeParametersNames

src/+parameters/setTempParameters.m

@@ -6,7 +6,7 @@
     % where landcoverClass is an array like {"forest", "forest", "forest", "shrubland", ...}
     landcovers = unique(landcoverClass, 'stable');
     for landcoverIndex = 1:length(landcovers)
-        [Vcmo, Tparam, m, Type, Rdparam, leafwidth] = landcover_variables(landcovers(landcoverIndex), siteName);
+        [Vcmo, Tparam, m, g1Med, g0Med, Type, Rdparam, leafwidth] = landcover_variables(landcovers(landcoverIndex), siteName);
         % The lc... parameters are intermediate values. The timeseries are generated in
         %   io.loadTimeSeries.m.
         ScopeParameters.lcVcmo(landcoverIndex, 1) = Vcmo;
@@ -15,41 +15,53 @@
         ScopeParameters.lcType(landcoverIndex, 1) = Type;
         ScopeParameters.lcRdparam(landcoverIndex, 1) = Rdparam;
         ScopeParameters.lcleafwidth(landcoverIndex, 1) = leafwidth;
+        ScopeParameters.lcg1Med(landcoverIndex, 1) = g1Med;
+        ScopeParameters.lcg0Med(landcoverIndex, 1) = g0Med;
     end
 end
 
 %% lookup table for Vcmo, Tparam, m, Type, Rdparam, and leafwidth
 
-function [Vcmo, Tparam, m, Type, Rdparam, leafwidth] = landcover_variables(landCoverType, siteName)
+function [Vcmo, Tparam, m, g1Med, g0Med, Type, Rdparam, leafwidth] = landcover_variables(landCoverType, siteName)
     Rdparam = [0.015]; % NOTE: Rdparam IS MISSING FOR MOST DEFINTIONS:
     if startsWith(landCoverType, 'Permanent Wetlands')
         Tparam = [0.2 0.3 288 313 328]; % These are five parameters specifying the temperature response.
         Vcmo = [120]; % Vcmax, maximum carboxylation capacity (at optimum temperature)
         m = [9]; % Ball-Berry stomatal conductance parameter
+        g1Med = [14.43]; % Medlyn stomatal conductance parameter: g1
+        g0Med = [0.001]; % Medlyn stomatal conductance parameter: g0
         Type = [0]; % Photochemical pathway: 0=C3, 1=C4
         leafwidth = [0.05]; % leaf width
     elseif startsWith(landCoverType, 'Evergreen Broadleaf')
         Tparam = [0.2 0.3 283 311 328];
         Vcmo = [80];
         m = [9];
+        g1Med = [9];
+        g0Med = [0.001];
         Type = [0];
         leafwidth = [0.05];
     elseif startsWith(landCoverType, 'Deciduous Broadleaf')
         Tparam = [0.2 0.3 283 311 328];
         Vcmo = [80];
         m = [9];
+        g1Med = [6.99];
+        g0Med = [-0.044];
         Type = [0];
         leafwidth = [0.05];
     elseif startsWith(landCoverType, 'Mixed Forests')
         Tparam = [0.2 0.3 281 307 328];
         Vcmo = [80];
         m = [9];
+        g1Med = [5.36];
+        g0Med = [0.024];
         Type = [0];
         leafwidth = [0.04];
     elseif startsWith(landCoverType, 'Evergreen Needleleaf')
         Tparam = [0.2 0.3 278 303 328];
         Vcmo = [80];
         m = [9];
+        g1Med = [1.66];
+        g0Med = [0.033];
         Type = [0];
         leafwidth = [0.01];
     elseif startsWith(landCoverType, 'Croplands')
@@ -67,28 +79,38 @@
             Type = [0];
             leafwidth = [0.03];
         end
+        g1Med = [9];
+        g0Med = [0.001];
     elseif startsWith(landCoverType, 'Open Shrublands')
         Tparam = [0.2 0.3 288 313 328];
         Vcmo = [120];
         m = [9];
+        g1Med = [5.95];
+        g0Med = [0.028];
         Type = [0];
         leafwidth = [0.05];
     elseif startsWith(landCoverType, 'Closed Shrublands')
         Tparam = [0.2 0.3 288 313 328];
         Vcmo = [80];
         m = [9];
+        g1Med = [5.95];
+        g0Med = [0.028];
         Type = [0];
         leafwidth = [0.05];
     elseif startsWith(landCoverType, 'Savannas')
         Tparam = [0.2 0.3 278 313 328];
         Vcmo = [120];
         m = [9];
+        g1Med = [11.23];
+        g0Med = [-0.004];
         Type = [0];
         leafwidth = [0.05];
     elseif startsWith(landCoverType, 'Woody Savannas')
         Tparam = [0.2 0.3 278 313 328];
         Vcmo = [120];
         m = [9];
+        g1Med = [11.23];
+        g0Med = [-0.004];
         Type = [0];
         leafwidth = [0.03];
     elseif startsWith(landCoverType, 'Grassland')
@@ -104,12 +126,16 @@
             Type = [0];
         end
         leafwidth = [0.02];
+        g1Med = [11.23];
+        g0Med = [-0.004];
     else
         Tparam = [0.2 0.3 288 313 328];
         Vcmo = [80];
         m = [9];
         Type = [0];
         leafwidth = [0.05];
+        g1Med = [9];
+        g0Med = [0.001];
         warning('landCover name unknown, "%s" is not recognized. ', landCoverType);
     end
 end

src/Dtrmn_Z.m

@@ -3,51 +3,75 @@
         The determination of the element length
     %}
     Elmn_Lnth = 0;
+    sitename = 'CH-HTC';
+    if strcmp(sitename,'CH-HTC')
+        for ML=1:6
+            DeltZ_R(ML)=1;%4
+        end
+    %         DeltZ_R(6)=1;%4
+        for ML=7:13
+            DeltZ_R(ML)=2;%5
+        end
+        for ML=14:21
+            DeltZ_R(ML)=2.5;%5
+        end
+        for ML=22:25
+            DeltZ_R(ML)=5;
+        end
+        for ML=26:29
+            DeltZ_R(ML)=10;
+        end
+        NL= ML;
+        for ML=1:NL
+            MML=NL-ML+1;
+            DeltZ(ML)=DeltZ_R(MML);
+        end
+    else
+        for ML = 1:3
+            DeltZ_R(ML) = 1; % 4
+        end
+        DeltZ_R(4) = 2; % 4
+        for ML = 5:14
+            DeltZ_R(ML) = 2; % 5
+        end
+        for ML = 15:18
+            DeltZ_R(ML) = 2.5; % 5
+        end
+        for ML = 19:23
+            DeltZ_R(ML) = 5;
+        end
+        for ML = 24:31
+            DeltZ_R(ML) = 10;
+        end
+        for ML = 32:40
+            DeltZ_R(ML) = 10;
+        end
+        for ML = 41:42
+            DeltZ_R(ML) = 15; % 5
+        end
+        % Sum of element lengths and compared to the total lenght, so that judge
+        % can be made to determine the length of rest elements.
 
-    for ML = 1:3
-        DeltZ_R(ML) = 1; % 4
-    end
-    DeltZ_R(4) = 2; % 4
-    for ML = 5:14
-        DeltZ_R(ML) = 2; % 5
-    end
-    for ML = 15:18
-        DeltZ_R(ML) = 2.5; % 5
-    end
-    for ML = 19:23
-        DeltZ_R(ML) = 5;
-    end
-    for ML = 24:31
-        DeltZ_R(ML) = 10;
-    end
-    for ML = 32:40
-        DeltZ_R(ML) = 10;
-    end
-    for ML = 41:42
-        DeltZ_R(ML) = 15; % 5
-    end
-    % Sum of element lengths and compared to the total lenght, so that judge
-    % can be made to determine the length of rest elements.
-
-    for ML = 1:42
-        Elmn_Lnth = Elmn_Lnth + DeltZ_R(ML);
-    end
+        for ML = 1:42
+            Elmn_Lnth = Elmn_Lnth + DeltZ_R(ML);
+        end
 
-    % If the total sum of element lenth is over the predefined depth, stop the
-    % for loop, make the ML, at which the element lenth sumtion is over defined
-    % depth, to be new NL.
-    DeltZ = [];
-    for ML = 43:NL
-        DeltZ_R(ML) = 20;
-        Elmn_Lnth = Elmn_Lnth + DeltZ_R(ML);
-        if Elmn_Lnth >= Tot_Depth
-            DeltZ_R(ML) = Tot_Depth - Elmn_Lnth + DeltZ_R(ML);
-            NL = ML;
-            for ML = 1:NL
-                MML = NL - ML + 1;
-                DeltZ(ML) = DeltZ_R(MML);
+        % If the total sum of element lenth is over the predefined depth, stop the
+        % for loop, make the ML, at which the element lenth sumtion is over defined
+        % depth, to be new NL.
+        DeltZ = [];
+        for ML = 43:NL
+            DeltZ_R(ML) = 20;
+            Elmn_Lnth = Elmn_Lnth + DeltZ_R(ML);
+            if Elmn_Lnth >= Tot_Depth
+                DeltZ_R(ML) = Tot_Depth - Elmn_Lnth + DeltZ_R(ML);
+                NL = ML;
+                for ML = 1:NL
+                    MML = NL - ML + 1;
+                    DeltZ(ML) = DeltZ_R(MML);
+                end
+                return
             end
-            return
         end
     end
 end