Merge pull request #58 from EcoExtreML/nb_era5cli

Use initial soil moisture and temperature from ERA5-land data

config_file_crib.txt

@@ -3,4 +3,5 @@ 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=AU-DaS_2010-2017_OzFlux_Met.nc
+InitialConditionPath=/data/shared/EcoExtreML/STEMMUS_SCOPEv1.0.0/input/Initial_condition/
 DurationSize=NA

src/+io/read_config.m

@@ -1,4 +1,4 @@
-function [SoilPropertyPath, InputPath, OutputPath, ForcingPath, ForcingFileName, DurationSize] = read_config(config_file)
+function [SoilPropertyPath, InputPath, OutputPath, ForcingPath, ForcingFileName, DurationSize, InitialConditionPath] = read_config(config_file)
 
 file_id = fopen(config_file);
 config = textscan(file_id,'%s %s', 'HeaderLines',0, 'Delimiter', '=');
@@ -24,6 +24,9 @@
 indx = find(strcmp(config_vars, 'ForcingFileName'));
 ForcingFileName = config_paths{indx};
 
+indx = find(strcmp(config_vars, 'InitialConditionPath'));
+InitialConditionPath = config_paths{indx};
+
 % value of DurationSize is optional and can be NA
 indx = find(strcmp(config_vars, 'DurationSize'));
 DurationSize = str2double(config_paths{indx});

src/Constants.m

@@ -3,7 +3,7 @@
 global xERR hERR TERR PERR tS uERR
 global KT TIME Delt_t0 DURTN TEND NIT KIT Nmsrmn Eqlspace h_SUR Msrmn_Fitting  
 global Msr_Mois Msr_Temp Msr_Time Tp_t Evap EPCT SAVEDTheta_LLT SAVEDTheta_LLh SAVEDTheta_UUh
-global Ksoil Rl SMC Ztot DeltZ_R Theta_o rroot frac bbx wfrac RWUtot RWUtott RWUtottt Rls Tatot LR PSItot sfactortot LAI_msr R_depth RTB VPD_msr Tss Tsss SUMTIME Tsur FOC FOS FOSL MSOC Coef_Lamda fieldMC DELT Dur_tot Precipp fmax
+global Ksoil Rl SMC Ztot DeltZ_R Theta_o rroot frac bbx wfrac RWUtot RWUtott RWUtottt Rls Tatot LR PSItot sfactortot LAI_msr R_depth RTB VPD_msr Tss Tsss SUMTIME Tsur FOC FOS FOSL MSOC Coef_Lamda fieldMC DELT Dur_tot Precipp fmax sitename
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%% The time and domain information setting. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -453,8 +453,10 @@
 Precip_msr=Mdata(:,6)'*10*DELT;
 Tmin=min(Ts_msr);
 Precip_msr=Precip_msr.*(1-fmax.*exp(-0.5*0.5*Tot_Depth/100));
+% load initial soil moisture and soil temperature from ERA5
+Initial_path=dir(fullfile(InitialConditionPath,[sitename,'*.nc']));
 InitND1=5;    % Unit of it is cm. These variables are used to indicated the depth corresponding to the measurement.
-InitND2=30;
+InitND2=15;
 InitND3=60;
 InitND4=100;
 InitND5=200;
@@ -477,21 +479,46 @@
 InitX6=	0.078;
 BtmX=0.078;%0.078 0.05;    % The initial moisture content at the bottom of the column.
 else
-InitT0= nanmean(Ta_msr);  %-1.75estimated soil surface temperature-1.762
-InitT1=	nanmean(Ta_msr);
-InitT2=	nanmean(Ta_msr);
-InitT3=	nanmean(Ta_msr);
-InitT4=	nanmean(Ta_msr);
-InitT5=	nanmean(Ta_msr);
-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.
+      InitT0=   ncread([InitialConditionPath,Initial_path(1).name],'skt')-273.15;  %-1.75estimated soil surface temperature-1.762
+      InitT1=	ncread([InitialConditionPath,Initial_path(2).name],'stl1')-273.15;
+      InitT2=	ncread([InitialConditionPath,Initial_path(3).name],'stl2')-273.15;
+      InitT3=	ncread([InitialConditionPath,Initial_path(4).name],'stl3')-273.15;
+      InitT4=	ncread([InitialConditionPath,Initial_path(5).name],'stl4')-273.15;
+      InitT5=	ncread([InitialConditionPath,Initial_path(5).name],'stl4')-273.15;
+      InitT6=	ncread([InitialConditionPath,Initial_path(5).name],'stl4')-273.15;
+      Tss = InitT0;
+    if InitT0 < 0 || InitT1 < 0 || InitT2 < 0 || InitT3 < 0 || InitT4 < 0 || InitT5 < 0 || InitT6 < 0 
+      InitT0=   0;
+      InitT1=	0;
+      InitT2=	0;
+      InitT3=	0;
+      InitT4=	0;
+      InitT5=	0;
+      InitT6=	0;
+      Tss = InitT0;
+    end
+if nanmean(Ta_msr)<0
+    BtmT  = 0;  %9 8.1
+else
+    BtmT  =  nanmean(Ta_msr);
+end
+InitX0=	ncread([InitialConditionPath,Initial_path(6).name],'swvl1');  %0.0793
+InitX1=	ncread([InitialConditionPath,Initial_path(6).name],'swvl1'); % Measured soil liquid moisture content
+InitX2=	ncread([InitialConditionPath,Initial_path(7).name],'swvl2'); %0.182
+InitX3=	ncread([InitialConditionPath,Initial_path(8).name],'swvl3');
+InitX4= ncread([InitialConditionPath,Initial_path(9).name],'swvl4'); %0.14335
+InitX5=	ncread([InitialConditionPath,Initial_path(9).name],'swvl4');
+InitX6=	ncread([InitialConditionPath,Initial_path(9).name],'swvl4');
+BtmX  = ncread([InitialConditionPath,Initial_path(9).name],'swvl4');%0.05;    % The initial moisture content at the bottom of the column.
+if InitX0 > SaturatedMC(1) || InitX1 > SaturatedMC(1) ||InitX2 > SaturatedMC(2) ||...
+InitX3 > SaturatedMC(3) || InitX4 > SaturatedMC(4) || InitX5 > SaturatedMC(5) || InitX6 > SaturatedMC(6)
+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);    
+end
 end

src/StartInit.m

@@ -12,7 +12,7 @@
 global HCAP TCON SF TCA GA1 GA2 GB1 GB2 S1 S2 HCD TARG1 TARG2 GRAT VPER 
 global TERM ZETA0 CON0 PS1 PS2 i KLT_Switch DVT_Switch KaT_Switch
 global Kaa_Switch DVa_Switch KLa_Switch
-global KfL_T Theta_II Theta_I Theta_UU Theta_U T_CRIT L_f g T0 TT_CRIT XUOLD XIOLD ISFT TCON_s TCON_dry TCON_L RHo_bulk Imped TPS1 TPS2 FEHCAP TS1 TS2 TCON0 TCON_min Theta_qtz XSOC
+global KfL_T Theta_II Theta_I Theta_UU Theta_U  L_f g T0 TT_CRIT XUOLD XIOLD ISFT TCON_s TCON_dry TCON_L RHo_bulk Imped TPS1 TPS2 FEHCAP TS1 TS2 TCON0 TCON_min Theta_qtz XSOC
 global Lamda Phi_s SWCC XCAP ThermCond Gama_hh Theta_a Gama_h hm hd SAVEhh SAVEh
 global COR CORh Theta_V Theta_g Se KL_h DTheta_LLh KfL_h DTheta_UUh hThmrl Tr Hystrs KIT RHOI RHOL FOC FOS FOSL MSOC Coef_Lamda fieldMC Theta_f Ta_msr IGBP_veg_long
 
@@ -543,9 +543,13 @@
 BChB=-9e-10; 
 if Thmrlefc==1
     NBCT=1;  % Energy Surface B.C.: 1 --> Specified temperature (BCT); 2 --> Specified heat flux (BCT); 3 --> Atmospheric forcing;
-    BCT=Ta_msr(1);  % surface temperature
+    BCT= Ta_msr(1);  % surface temperature
     NBCTB=1;% Energy Bottom B.C.: 1 --> Specified temperature (BCTB); 2 --> Specified heat flux (BCTB); 3 --> Zero temperature gradient;
-    BCTB=nanmean(Ta_msr);
+    if nanmean(Ta_msr)<0
+        BCTB  = 0;  %9 8.1
+    else
+        BCTB  =  nanmean(Ta_msr);
+    end
 end
 if Soilairefc==1
     NBCP=2; % Soil air pressure B.C.: 1 --> Ponded infiltration caused a specified pressure value; 

src/ebal.m

@@ -364,10 +364,10 @@
     % 2.7. New estimates of soil (s) and leaf (c) temperatures, shaded (h) and sunlit (1) 
     Tch         = Tch + Wc*EBerch./((rhoa*cp)./rac + rhoa*lambdah*e_to_q.*sh./(rac+rcwh)+ 4*leafbio.emis*sigmaSB*(Tch+273.15).^3);
     Tcu         = Tcu + Wc*EBercu./((rhoa*cp)./rac + rhoa*lambdau*e_to_q.*su./(rac+rcwu)+ 4*leafbio.emis*sigmaSB*(Tcu+273.15).^3);
-    Ts          = Ts + Wc*EBers./(rhoa*cp./ras + rhoa*lambdas*e_to_q.*ss/(ras+rss)+ 4*(1-soil.rs_thermal)*sigmaSB*(Ts+273.15).^3);  
+    Ts          = Ts + Wc*EBers./(rhoa*cp./ras + rhoa*lambdas*e_to_q.*ss/(ras+rss)+ 4*(1-soil.rs_thermal)*sigmaSB*(Ts+273.15).^3); % Ts contains shaded soil temperature and sunlit soil temperature 
     Tch(abs(Tch)>100) = Ta;
     Tcu(abs(Tcu)>100) = Ta;
-
+    Ts(abs(Ts)>100) = Ta;
     if (any(isnan(Tch)) || any(isnan(Tcu(:)))), warning('Canopy temperature gives NaNs'), end
     if any(isnan(Ts)), warning('Soil temperature gives NaNs'), end
 

src/filesread.m

@@ -9,7 +9,7 @@
 
 %% Read the CFG file. Due to using MATLAB compiler, we cannot use run(CFG)
 disp (['Reading config from ',CFG])
-[SoilPropertyPath, InputPath, OutputPath, ForcingPath, ForcingFileName, DurationSize] = io.read_config(CFG);
+[SoilPropertyPath, InputPath, OutputPath, ForcingPath, ForcingFileName, DurationSize, InitialConditionPath] = io.read_config(CFG);
 
 %%%%%%% Prepare input files. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 global DELT IGBP_veg_long latitude longitude reference_height canopy_height sitename
@@ -127,7 +127,7 @@
 LAI=ncread(ForcingFilePath,'LAI');
 LAIL=length(LAI);
 LAIa=reshape(LAI,LAIL,1);
-
+LAIa(LAIa<0.01)=0.01;
 
 LAI_alternative=ncread(ForcingFilePath,'LAI_alternative');
 LAI_alternativeL=length(LAI_alternative);
@@ -162,4 +162,4 @@
 save([InputPath, 'Mdata.txt'], '-ascii', 'Meteodata') %save meteorological data for STEMMUS
 %Lacationdata=[latitude;longitude;reference_height;canopy_height;elevation]';
 %save([InputPath, 'Lacationdata.txt'], '-ascii', 'Lacationdata')
-clearvars -except SoilPropertyPath InputPath OutputPath DELT Dur_tot IGBP_veg_long latitude longitude reference_height canopy_height sitename
+clearvars -except SoilPropertyPath InputPath OutputPath InitialConditionPath DELT Dur_tot IGBP_veg_long latitude longitude reference_height canopy_height sitename