Merge pull request #1 from EcoExtreML/add_code_v1_gmd

Add codes of gmd paper, version 1

filenames.m

@@ -0,0 +1,36 @@
+%The following three are always required,
+
+X = {
+'Simulation_Name'	, 'verification_run';
+'soil_file'		, 'soilnew.txt';
+'leaf_file'		, 'Optipar2017_ProspectD.mat'; 
+'atmos_file' 		, 'FLEX-S3_std.atm';
+
+%The following are only for the time series option!
+'Dataset_dir'		, 'for_verification'; 
+'t_file'		, 't_.dat';
+'year_file'		, 'year_.dat';
+'Rin_file'		, 'Rin_.dat';
+'Rli_file'		, 'Rli_.dat';
+'p_file' 		, 'p_.dat';
+'Ta_file'		, 'Ta_.dat';
+'ea_file'		, 'ea_.dat';
+'u_file'		, 'u_.dat';
+
+%optional (leave empty for constant values From inputdata.TXT)
+'CO2_file'		, '';
+'SMC_file'		, '';
+
+% optional (leave empty for calculations based on t_file year timezn)
+'tts_file' 		, '';
+
+%optional two column tables (first column DOY second column value)
+'z_file' 		, '';
+'LAI_file'		, '';
+'hc_file'		, '';
+'Vcmax_file'	, '';
+'Cab_file'		, '';
+
+%optional leaf inclination distribution file with 3 headerlines (see
+%example). It MUST be located in ../data/leafangles/
+'LIDF_file'     , ''};

set_parameter_filenames.m

@@ -0,0 +1,2 @@
+% parameter_file             = { 'setoptions.m', 'filenames.m', 'inputdata.txt'};
+parameter_file            = {'input_data.xlsx'};

setoptions.m

@@ -0,0 +1,21 @@
+N=[
+1;		%calc_ebal		calculate the complete energy balance
+0;		%calc_vert_profiles		calculate vertical profiles of fluxes and temperatures
+1;		%calc_fluor		calculate chlorophyll fluorescence
+0;		%calc_planck		calculate spectrum of thermal radiation with spectral emissivity instead of broadband
+0;		%calc_directional		calculate BRDF and directional temperature for many angles specified in a file. Be patient, this takes some time
+1;      %calc_xanthophyllabs		calculate dynamic xanthopyll absorption (zeaxanthin)
+0;      %calc_PSI       0 (recommended): treat the whole fluorescence spectrum as one spectrum (new calibrated optipar), 1: differentiate PSI and PSII with Franck et al. spectra (of SCOPE 1.62 and older)
+0;		%rt_thermal		0: provide emissivity values as input. 1: use values from fluspect and soil at 2400 nm for the TIR range
+0;		%calc_zo		0: use the zo and d values provided in the inputdata, 1: calculate zo and d from the LAI, canopy height, CD1, CR, CSSOIL (recommended if LAI changes in time series)
+0;      %0: use soil spectrum from a file, 1: simulate soil spectrum with the BSM model
+0;		%SoilHeatMethod		0: standard calculation of thermal inertia from soil characteristics, 1: empiricaly calibrated formula (make function), 2: as constant fraction of soil net radiation
+1;		%Fluorescence_model		0: empirical, with sustained NPQ (fit to Flexas' data); 1: empirical, with sigmoid for Kn; 2: Magnani 2012 model
+0;		%calc_rss_rbs		0: use resistance rss and rbs as provided in inputdata. 1:  calculate rss from soil moisture content and correct rbs for LAI (calc_rssrbs.m)
+1;		%applTcorr		correct Vcmax and rate constants for temperature in biochemical.m
+1;		%verify		verifiy the results (compare to saved 'standard' output) to test the code for the firstt ime
+1;		%saveheaders		write header lines in output files
+0;		%makeplots 		plot the results
+1];		%simulation		0: individual runs. Specify one value for constant input, and an equal number (>1) of values for all input that varies between the runs.
+		%		1: time series (uses text files with meteo input as time series)
+		%		2: Lookup-Table (specify the values to be included. All possible combinations of inputs will be  used)

src/+equations/Planck.m

@@ -0,0 +1,11 @@
+function Lb = Planck(wl,Tb,em)
+
+    c1 = 1.191066e-22;
+    c2 = 14388.33;
+    if nargin<3
+        em = ones(size(Tb));
+    end
+
+    Lb = em.* c1*(wl*1e-9).^(-5)./(exp(c2./(wl*1e-3*Tb))-1);
+
+end    

src/+equations/Soil_Inertia0.m

@@ -0,0 +1,3 @@
+function [GAM]  = Soil_Inertia0(cs,rhos,lambdas)
+% soil thermal inertia 
+GAM             = sqrt(cs*rhos*lambdas);                            % soil thermal intertia

src/+equations/Soil_Inertia1.m

@@ -0,0 +1,35 @@
+function [GAM] = Soil_Inertia1(SMC)
+
+%soil inertia method by Murray and Verhoef (
+
+%% parameters
+
+theta_s = 0.42; %(saturated water content, m3/m3)
+Sr = SMC/theta_s;
+
+%fss = 0.58; %(sand fraction)
+gamma_s = 0.27; %(soil texture dependent parameter)
+dels = 1.33; %(shape parameter)
+
+
+ke = exp(gamma_s*(1- power(Sr,(gamma_s - dels))));
+
+phis  = 0.5; %(phis == theta_s)
+lambda_d = -0.56*phis + 0.51;
+
+QC = 0.20; %(quartz content)
+lambda_qc = 7.7;  %(thermal conductivity of quartz, constant)
+
+lambda_s = (lambda_qc^(QC))*lambda_d^(1-QC);
+lambda_wtr = 0.57;   %(thermal conductivity of water, W/m.K, constant)
+
+lambda_w = (lambda_s^(1-phis))*lambda_wtr^(phis);
+
+lambdas = ke*(lambda_w - lambda_d) + lambda_d;
+
+Hcs = 2.0*10^6;
+Hcw = 4.2*10^6;
+
+Hc = (Hcw * SMC)+ (1-theta_s)*Hcs;
+
+GAM = sqrt(lambdas.*Hc);

src/+equations/calc_rssrbs.m

@@ -0,0 +1,14 @@
+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 =11.2*exp(0.3563*100.0*(0.38-SMC));
+%rss = exp(7.9-1.6*(SMC-0.0008)/(0.38-0.0008));
+%if  rss<70,
+  %  rss=70;
+%end
+rbs            = rbs*LAI/4.3;

src/+equations/calczenithangle.m

@@ -0,0 +1,64 @@
+function [Fi_s,Fi_gs,Fi_g,Omega_s] = calczenithangle(Doy,t,Omega_g,Fi_gm,Long,Lat)
+%
+% author: Christiaan van der Tol ([email protected])
+% date:     Jan 2003
+% update:   Oct 2008 by Joris Timmermans ([email protected]): 
+%               - corrected equation of time
+%           Oct 2012 (CvdT) comment: input time is GMT, not local time!
+%
+% function [Fi_s,Fi_gs,Fi_g]= calczenithangle(Doy,t,Omega_g,Fi_gm,Long,Lat)
+%
+% calculates pi/2-the angle of the sun with the slope of the surface.
+%
+% input:
+% Doy       day of the year
+% t         time of the day (hours, GMT)
+% Omega_g   slope azimuth angle (deg)
+% Fi_gm     slope of the surface (deg)
+% Long      Longitude (decimal)
+% Lat       Latitude (decimal)
+%
+% output:
+% Fi_s      'classic' zenith angle: perpendicular to horizontal plane
+% Fi_gs     solar angle perpendicular to surface slope
+% Fi_g      projected slope of the surface in the plane through the solar beam and the vertical
+% 
+
+%parameters (if not already supplied)
+if nargin<6
+    Long        =   13.75;                      % longitude
+    Lat         =   45.5;                       % latitude
+    if (nargin<4)
+        Omega_g =   210;                        % aspect
+        Fi_gm   =   30;                         % slope angle
+    end
+end
+
+%convert angles into radials
+G               =   (Doy-1)/365*2*pi;           % converts day of year to radials
+Omega_g         =   Omega_g/180*pi;             % converts direction of slope to radials
+Fi_gm           =   Fi_gm/180*pi;               % converts maximum slope to radials
+Lat             =   Lat/180*pi;                 % converts latitude to radials
+
+%computes the declination of the sun
+d               =   0.006918-0.399912*cos(G  )+ 0.070247*sin(G  )- ...
+                     0.006758*cos(2*G)+ 0.000907*sin(2*G)- ...
+                     0.002697*cos(3*G)+ 0.00148*sin(3*G);
+
+%Equation of time
+Et              =   0.017 + .4281 * cos(G) - 7.351 * sin(G) - 3.349 * cos(2*G) - 9.731 * sin(2*G);
+
+%computes the time of the day when the sun reaches its highest angle                                
+tm              =   12+(4*(-Long)-Et)/60;      % de Pury and Farquhar (1997), Iqbal (1983)
+
+%computes the hour angle of the sun
+Omega_s         =   (t-tm)/12*pi;
+
+%computes the zenithangle (equation 3.28 in De Bruin)
+Fi_s            =   acos(sin(d)*sin(Lat)+cos(d)*cos(Lat).*cos(Omega_s));
+
+%computes the slope of the surface Fi_g in the same plane as the solar beam
+Fi_g            =   atan(tan(Fi_gm).*cos(Omega_s-Omega_g));
+
+%computes the angle of the sun with the vector perpendicular to the surface
+Fi_gs           =   Fi_s + Fi_g;