Merge pull request #39 from EcoExtreML/next_step_work

Next changes to the model to use high resolution soil data

config_file_crib.txt

@@ -1,6 +1,6 @@
 SoilPropertyPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/SoilProperty/
-InputPath=/home/jovyan/private/06_STEMMUS_SCOPE/GitEcoExtreML/data/
-OutputPath=/home/jovyan/private/06_STEMMUS_SCOPE/GitEcoExtreML/data/output/
+InputPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/
+OutputPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/output/
 ForcingPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/Plumber2_data/
-ForcingFileName=FI-Hyy_1996-2014_FLUXNET2015_Met.nc
-DurationSize=NA
+ForcingFileName=AU-DaS_2010-2017_OzFlux_Met.nc
+DurationSize=NA

src/+equations/calc_rssrbs.m

@@ -1,13 +1,5 @@
 function [rss,rbs] = calc_rssrbs(SMC,LAI,rbs)
-
-%rss            = 10*exp(35.63*(0.25-SMC));
-%if rss>1000,
-    %rss=1000;
-%elseif rss<30,
-    %rss=30;
-%end
-rss = exp(7.6-1.5*(SMC-0.0875)/(0.25-0.0875));
-%if  rss<70,
-  %  rss=70;
-%end
+global SaturatedMC ResidualMC fieldMC
+aa=3.8;
+rss = exp((aa+4.1)-aa*(SMC-ResidualMC(1))/(fieldMC(1)-ResidualMC(1)));
 rbs            = rbs*LAI/4.3;

src/+io/output_verification_csv.m

@@ -64,7 +64,7 @@ function output_verification_csv(Output_dir, verification_dir)
                 D0 = dlmread([path0_ info0(i).name],',',1,0);
                 D1 = dlmread([path1_ info1(j).name],',',1,0);
             end
-            if length(D0) ~= length(D1), keyboard, end
+            if length(D0) ~= length(D1), fprintf('\n Warning: FIX_ME: replace keyboard \r'), end
             if (sum(sum(D0-D1).^2))>1E-9
                 fprintf(['\nWarning: data in the output file ' info0(i).name ' are different from the verification output \r '])
                 h0 = textread([path0_ info0(i).name],'%s','bufsize', 1E9); %#ok<DTXTRD>

src/CondL_h.m

@@ -236,7 +236,8 @@
             Se(ML,ND)=0;
         end
         if isnan(Se(ML,ND))==1
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         Theta_II(ML,ND)=(Theta_UU(ML,ND)-Theta_LL(ML,ND))*RHOL/RHOI;  % ice water content
         if Theta_UU(ML,ND)~=0
@@ -340,10 +341,12 @@
                 end
             if isnan(KL_h(ML,ND))==1
                 KL_h(ML,ND)=0;
-                keyboard
+                %%keyboard
+                fprintf('\n Warning: FIX_ME: replace keyboard \r')
             end
             if ~isreal(KL_h(ML,ND))
-                keyboard
+                %%keyboard
+                fprintf('\n Warning: FIX_ME: replace keyboard \r')
             end
             KfL_T(ML,ND)=helpers.heaviside1(TT_CRIT(MN)-(TT(MN)+T0))*L_f*1e4/(g*(T0));   % thermal consider for freezing soil
         else

src/Constants.m

@@ -29,7 +29,7 @@
 uERR=0.02;                  % Maximum desirable change of total water content;                                                                    
 Tot_Depth=500;           % Unit is cm. it should be usually bigger than 0.5m. Otherwise, 
                                  % the DeltZ would be reset in 50cm by hand;                                  
-R_depth=300; %¸ùÉî                                                                               
+R_depth=300; %                                                                               
 Eqlspace=0;                 % Indicator for deciding is the space step equal or not;       
 NL=100;
 if ~Eqlspace                                                                     
@@ -486,12 +486,12 @@
 InitT6=	nanmean(Ta_msr);
 Tss=InitT0;
 BtmT=nanmean(Ta_msr);  %9 8.1
-InitX0=	fieldMC(1)*0.6;  %0.0793
-InitX1=	fieldMC(1)*0.6; % Measured soil liquid moisture content
-InitX2=	fieldMC(2)*0.6; %0.182
-InitX3=	fieldMC(3)*0.6;
-InitX4= fieldMC(4)*0.6; %0.14335
-InitX5=	fieldMC(5)*0.6;
-InitX6=	fieldMC(6)*0.6;
-BtmX=fieldMC(6)*0.6;%0.05;    % The initial moisture content at the bottom of the column.
+InitX0=	fieldMC(1);  %0.0793
+InitX1=	fieldMC(1); % Measured soil liquid moisture content
+InitX2=	fieldMC(2); %0.182
+InitX3=	fieldMC(3);
+InitX4= fieldMC(4); %0.14335
+InitX5=	fieldMC(5);
+InitX6=	fieldMC(6);
+BtmX=fieldMC(6);%0.05;    % The initial moisture content at the bottom of the column.
 end

src/Initial_root_biomass.m

@@ -4,7 +4,7 @@
 root_den = 250*1000; %% [gDM / m^3] Root density  Jackson et al., 1997 
 R_C = 0.488; %% [gC/gDM] Ratio Carbon-Dry Matter in root   Jackson et al.,  1997 
 if     IGBP_veg_long(1:2)== ['P'; 'e']     %['Permanent Wetlands'] 
-    beta = 0.993; %% ? is a plant-dependent root distribution parameter adopted from Jackson et al. (1996) 
+    beta = 0.993; %% beta is a plant-dependent root distribution parameter adopted from Jackson et al. (1996), it is refered to CLM5.0 document. 
 elseif IGBP_veg_long(1:2)== ['E'; 'v']     %['Evergreen Broadleaf'] 
     beta = 0.993; 
 elseif IGBP_veg_long(1:2)== ['D'; 'e']     %['Deciduous Broadleaf'] 
@@ -16,15 +16,15 @@
 elseif IGBP_veg_long(1:2)== ['C'; 'r']     %['Croplands'] 
     beta = 0.943; 
 elseif IGBP_veg_long(1:2)== ['O'; 'p']     %['Open Shrublands']
-    beta = 0.943; 
+    beta = 0.966; 
 elseif IGBP_veg_long(1:2)== ['C'; 'l']     %['Closed Shrublands'] 
-    beta = 0.943; 
+    beta = 0.966; 
 elseif IGBP_veg_long(1:2)== ['S'; 'a']     %['Savannas'] 
     beta = 0.943; 
 elseif IGBP_veg_long(1:2)== ['W'; 'o']     %['Woody Savannas'] 
     beta = 0.943; 
 else %IGBP_veg_long(1)==['Grasslands'] 
-    beta = 0.993; 
+    beta = 0.943; 
 end
 
 Rltot = RTB/R_C/root_den/(pi*(rroot^2)); %% %% root length index [m root / m^2 PFT]

src/STEMMUS_SCOPE.m

@@ -54,7 +54,7 @@
 global HCAP SF TCA GA1 GA2 GB1 GB2 HCD ZETA0 CON0 PS1 PS2 XWILT FEHCAP QMTT QMBB Evapo trap RnSOIL PrecipO
 global constants
 global RWU EVAP theta_s0 Ks0
-global HR Precip Precipp Tss frac sfactortot sfactor Tsss fluxes lEstot lEctot NoTime DELT IGBP_veg_long latitude longitude reference_height canopy_height sitename Dur_tot Tmin fmax
+global HR Precip Precipp Tss frac sfactortot sfactor fluxes lEstot lEctot NoTime DELT IGBP_veg_long latitude longitude reference_height canopy_height sitename Dur_tot Tmin fmax
 
 %% 1. define constants
 [constants] = io.define_constants();
@@ -199,52 +199,91 @@
 V(23).Val=canopy_height;
 V(55).Val=mean(Ta_msr);
 %Input T parameters for different vegetation type
+    sitename1=cellstr(sitename);
 if     IGBP_veg_long(1:2)== ['P'; 'e']     %['Permanent Wetlands'] 
-   V(14).Val= [0.2 0.3 288 313 328];
+    V(14).Val= [0.2 0.3 288 313 328]; % These are five parameters specifying the temperature response.
+    V(9).Val= [120]; % Vcmax, maximum carboxylation capacity (at optimum temperature)
+    V(10).Val= [9]; % Ball-Berry stomatal conductance parameter
+    V(11).Val= [0]; % Photochemical pathway: 0=C3, 1=C4
+    V(28).Val= [0.05]; % leaf width
 elseif IGBP_veg_long(1:2)== ['E'; 'v']     %['Evergreen Broadleaf'] 
     V(14).Val= [0.2 0.3 283 311 328];
+    V(9).Val= [80];
+    V(10).Val= [9];
+    V(11).Val= [0];
+    V(28).Val= [0.05];
 elseif IGBP_veg_long(1:2)== ['D'; 'e']     %['Deciduous Broadleaf'] 
     V(14).Val= [0.2 0.3 283 311 328];
+    V(9).Val= [80];
+    V(10).Val= [9];
+    V(11).Val= [0];  
+    V(28).Val= [0.05];
 elseif IGBP_veg_long(1:2)== ['M'; 'i']     %['Mixed Forests'] 
     V(14).Val= [0.2 0.3 281 307 328];
+    V(9).Val= [80];
+    V(10).Val= [9];
+    V(11).Val= [0]; 
+    V(28).Val= [0.04];
 elseif IGBP_veg_long(11:12)== ['N'; 'e']     %['Evergreen Needleleaf'] 
     V(14).Val= [0.2 0.3 278 303 328];
-elseif IGBP_veg_long(1:2)== ['C'; 'r']     %['Croplands'] 
-    V(14).Val= [0.2 0.3 278 303 328];
+    V(9).Val= [80];
+    V(10).Val= [9];
+    V(11).Val= [0];   
+    V(28).Val= [0.01];
+elseif IGBP_veg_long(1:2)== ['C'; 'r']     %['Croplands']    
+    if isequal(sitename1,{'ES-ES2'})||isequal(sitename1,{'FR-Gri'})||isequal(sitename1,{'US-ARM'})||isequal(sitename1,{'US-Ne1'})
+        V(14).Val= [0.2 0.3 278 303 328];
+        V(9).Val= [50];
+        V(10).Val= [4];
+        V(11).Val= [1]; 
+        V(13).Val= [0.025]; % Respiration = Rdparam*Vcmcax
+        V(28).Val= [0.03];
+    else 
+        V(14).Val= [0.2 0.3 278 303 328];
+        V(9).Val= [120];
+        V(10).Val= [9];
+        V(11).Val= [0]; 
+        V(28).Val= [0.03];    
+    end
 elseif IGBP_veg_long(1:2)== ['O'; 'p']     %['Open Shrublands']
     V(14).Val= [0.2 0.3 288 313 328];
+    V(9).Val= [120];
+    V(10).Val= [9];
+    V(11).Val= [0];  
+    V(28).Val= [0.05];
 elseif IGBP_veg_long(1:2)== ['C'; 'l']     %['Closed Shrublands'] 
     V(14).Val= [0.2 0.3 288 313 328];
+    V(9).Val= [80];
+    V(10).Val= [9];
+    V(11).Val= [0];
+    V(28).Val= [0.05];
 elseif IGBP_veg_long(1:2)== ['S'; 'a']     %['Savannas'] 
     V(14).Val= [0.2 0.3 278 313 328];
+    V(9).Val= [120];
+    V(10).Val= [9];
+    V(11).Val= [0];
+    V(28).Val= [0.05];
 elseif IGBP_veg_long(1:2)== ['W'; 'o']     %['Woody Savannas'] 
     V(14).Val= [0.2 0.3 278 313 328];
-else %IGBP_veg_long(1)==['Grasslands'] 
+    V(9).Val= [120];
+    V(10).Val= [9];
+    V(11).Val= [0];
+    V(28).Val= [0.03];
+else %IGBP_veg_long==['Grasslands'] 
     V(14).Val= [0.2 0.3 288 303 328];
+    if isequal(sitename1,{'AR-SLu'})||isequal(sitename1,{'AU-Ync'})||isequal(sitename1,{'CH-Oe1'})||isequal(sitename1,{'DK-Lva'})||isequal(sitename1,{'US-AR1'})||isequal(sitename1,{'US-AR2'})||isequal(sitename1,{'US-Aud'})||isequal(sitename1,{'US-SRG'})
+        V(9).Val= [120];
+        V(10).Val= [4];
+        V(11).Val= [1];
+        V(13).Val= [0.025]; 
+    else
+        V(9).Val= [120];
+        V(10).Val= [9];
+        V(11).Val= [0];
+    end
+    V(28).Val= [0.02];
 end
-if     IGBP_veg_long(1:2)== ['P'; 'e']     %['Permanent Wetlands'] 
-   V(28).Val= [0.05];
-elseif IGBP_veg_long(1:2)== ['E'; 'v']     %['Evergreen Broadleaf'] 
-    V(28).Val= [0.05];
-elseif IGBP_veg_long(1:2)== ['D'; 'e']     %['Deciduous Broadleaf'] 
-    V(28).Val= [0.05];
-elseif IGBP_veg_long(1:2)== ['M'; 'i']     %['Mixed Forests'] 
-    V(28).Val= [0.05];
-elseif IGBP_veg_long(11:12)== ['N'; 'e']     %['Evergreen Needleleaf'] 
-    V(28).Val= [0.05];
-elseif IGBP_veg_long(1:2)== ['C'; 'r']     %['Croplands'] 
-    V(28).Val= [0.05];
-elseif IGBP_veg_long(1:2)== ['O'; 'p']     %['Open Shrublands']
-    V(28).Val= [0.05];
-elseif IGBP_veg_long(1:2)== ['C'; 'l']     %['Closed Shrublands'] 
-    V(28).Val= [0.05];
-elseif IGBP_veg_long(1:2)== ['S'; 'a']     %['Savannas'] 
-    V(28).Val= [0.05];
-elseif IGBP_veg_long(1:2)== ['W'; 'o']     %['Woody Savannas'] 
-    V(28).Val= [0.05];
-else %IGBP_veg_long(1)==['Grasslands'] 
-    V(28).Val= [0.05];
-end
+
 TZ=fix(longitude/15);
 TZZ=mod(longitude,15);
 if longitude>0
@@ -566,27 +605,23 @@
            Acc=fluxes.Actot;
            lEstot =fluxes.lEstot;
            lEctot =fluxes.lEctot;
-           Tsss=Tss;
            Tss=thermal.Tsave; 
         else
            Acc=0;
            lEstot =0;
            lEctot =0;
-           Tsss=Tss;
            Tss=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;
            lEstot =fluxes.lEstot;
            lEctot =fluxes.lEctot;
-           Tsss=Tss;
            Tss=thermal.Tsave; 
         else
            Acc=0;
            lEstot =0;
            lEctot =0;
-           Tsss=Tss;
            Tss=Ta_msr(KT); 
         end
     end

src/TimestepCHK.m

@@ -1,7 +1,8 @@
 function [KT,TIME,Delt_t,tS,NBChh,IRPT1,IRPT2]=TimestepCHK(NBCh,BCh,NBChh,DSTOR,DSTOR0,EXCESS,QMT,Precip,Evap,hh,IRPT1,NN,SAVEhh,TT_CRIT,T0,TT,T,EPCT,h,T_CRIT,xERR,hERR,TERR,Theta_LL,Theta_L,Theta_UU,Theta_U,KT,TIME,Delt_t,NL,Thmrlefc,NBChB,NBCT,NBCTB,tS,uERR)
 
     if Delt_t<1.0e-10
-        keyboard  
+        %%keyboard
+        fprintf('\n Warning: FIX_ME: replace keyboard \r')
     end
 % testing this function %%%
 if NBCh==1

src/calc_rsoil.m

@@ -1,11 +1,7 @@
 function [PSIs,rsss,rrr,rxx] = calc_rsoil(Rl,DeltZ,Ks,Theta_s,Theta_r,Theta_LL,bbx,m,n,Alpha)
-%global SSMM RRll ZZZ KT
 DeltZ0=DeltZ';
 l=0.5;
 SMC=Theta_LL(1:54,2);
-%SSMM(KT)=sum(SMC.*(DeltZ/100).*bbx)/sum((DeltZ/100).*bbx);
-%RRll(KT)=sum(Rl.*bbx.*(DeltZ/100))/sum((DeltZ/100).*bbx);
-%ZZZ(KT)=sum((DeltZ/100).*bbx);
 Se  = (SMC-Theta_r')./(Theta_s'-Theta_r');
 Ksoil=Ks'.*Se.^l.*(1-(1-Se.^(1./m')).^(m')).^2;
 PSIs=-((Se.^(-1./m')-1).^(1./n'))./(Alpha*100)'.*bbx;

src/ebal.m

@@ -321,7 +321,6 @@
     BB1=AA1(~isnan(AA1));
     BB2=AA2(~isinf(AA2));
     PSI1 = (sum(BB1)-Trans)/sum(BB2);
-   % ��λҪͳһ�����뻹�ǰ�Сʱ���������Ƿ�Ҫ����
     if isnan(PSI1)
     PSI1 = -1; 
     end

src/hPARM.m

@@ -59,34 +59,44 @@
             VvT(ML,ND)=0;
         end
         if isnan(Chh(ML,ND))
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         if isnan(Khh(ML,ND))==1
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         if isnan(Chg(ML,ND))
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         if isnan(ChT(ML,ND))
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         if isnan(KhT(ML,ND))
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         if ~isreal(Chh(ML,ND))
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         if ~isreal(Khh(ML,ND))==1
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         if ~isreal(Chg(ML,ND))
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         if ~isreal(ChT(ML,ND))
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
         if ~isreal(KhT(ML,ND))
-            keyboard
+            %%keyboard
+            fprintf('\n Warning: FIX_ME: replace keyboard \r')
         end
     end
 end 

src/hh_Solve.m

@@ -17,8 +17,10 @@
 end
 
 if isnan(SAVEhh)==1
-    keyboard
+    %%keyboard
+    fprintf('\n Warning: FIX_ME: replace keyboard \r')
 end
 if ~isreal(SAVEhh)
-    keyboard
+    %%keyboard
+    fprintf('\n Warning: FIX_ME: replace keyboard \r')
 end