Merge pull request #4 from EcoExtreML/add_code_v1sth

Add modified codes to  source code v1

src/+equations/Root_properties.m

@@ -0,0 +1,16 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%   Subfunction - Root - Properties         %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%REFERENCES   
+function[Rl]=Root_properties(Rl,rroot,frac,BR)
+%%% INPUTS
+% BR = 10:1:650; %% [gC /m^2 PFT]
+% rroot =  0.5*1e-3 ; % 3.3*1e-4 ;%% [0.5-6 *10^-3] [m] root radius
+%%% OUTPUTS 
+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 
+Delta_Rltot = BR/R_C/root_den/(pi*(rroot^2)); %% %% root length index [m root / m^2 PFT]
+Delta_Rl = Delta_Rltot*frac;
+Rl = Rl + Delta_Rl;
+end 

src/+equations/calc_rssrbs.m

@@ -6,8 +6,7 @@
 %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));
+rss = exp(7.6-1.5*(SMC-0.0875)/(0.25-0.0875));
 %if  rss<70,
   %  rss=70;
 %end

src/+equations/fixedp_brent_ari.m

@@ -328,7 +328,7 @@
     %fprintf('Laggards: %d\n', sum(err_outside_tol));
     counter = counter + 1;
     if (counter > iter_limit)
-        error('iteration limit exceeded');
+        return%error('iteration limit exceeded');
     end
     %-----------------
     %  Now reorganize a, b, c so b is in best

src/+equations/soil_respiration.m

@@ -1,2 +1,2 @@
 function [R] = soil_respiration(Ts)
-R = 0.5+0.14375*Ts;    %umol m-2 s-1
+R = 2.5 + 0.054375*Ts;    %umol m-2 s-1

src/+io/bin_to_csv.m

@@ -0,0 +1,160 @@
+function bin_to_csv(fnames, V, vmax, n_col, ns, options, Ztot)
+%% flu
+flu_names = {'simulation_number','nu_iterations','year','DoY',...
+    'Rntot','lEtot','Htot', ...
+    'Rnctot','lEctot','Hctot','Actot', ...
+    'Rnstot','lEstot','Hstot','Gtot','Resp', ...
+   'aPAR','aPAR_Cab','aPAR/rad.PAR','aPAR_Wm2','PAR','rad.Eoutf','rad.Eoutf./fluxes.aPAR_Wm2','Trap','Evap','ET'};
+flu_units = {'','','','',  ...
+    'W m-2','W m-2','W m-2',...
+    'W m-2','W m-2','W m-2','umol m-2 s-1',...
+    'W m-2',' W m-2','W m-2','W m-2','umol m-2 s-1',...
+    'umol m-2 s-1',' umol m-2 s-1','umol m-2 s-1','W m-2','umol m-2 s-1','W m-2','W m-2','cm s-1','cm s-1','cm s-1'};
+write_output(flu_names, flu_units, fnames.flu_file, n_col.flu, ns)
+
+%% surftemp
+surftemp_names = {'simulation_number','year','DoY',...
+    'Ta','Ts(1)','Ts(2)','Tcave','Tsave'};
+surftemp_units = {'','','',  ...
+    'C','C','C','C','C'};
+write_output(surftemp_names, surftemp_units, fnames.surftemp_file, n_col.surftemp, ns)
+
+%% aerodyn
+aerodyn_names = {'simulation_number','year','DoY',...
+    'raa','rawc','raws','ustar','rac','ras'};
+aerodyn_units = {'','','',  ...
+    's m-1','s m-1','s m-1', 'm s-1','s m-1','s m-1'};
+write_output(aerodyn_names, aerodyn_units, fnames.aerodyn_file, n_col.aerodyn, ns)
+
+%% radiation
+radiation_names = {'simulation_number','year','DoY',...
+    'Rin','Rli','HemisOutShort','HemisOutLong','HemisOutTot','Netshort','Netlong','Rntot'};
+radiation_units = {'','','',  ...
+    'W m-2','W m-2','W m-2', 'W m-2','W m-2','W m-2','W m-2','W m-2'};
+write_output(radiation_names, radiation_units, fnames.radiation_file, n_col.radiation, ns)
+%% Soil temperature and moisture
+Sim_Theta_names = num2str(Ztot');
+Sim_Theta_names=cellstr(Sim_Theta_names);
+Sim_Theta_units ={'m-3 m-3'};
+write_output(Sim_Theta_names, Sim_Theta_units, ...
+fnames.Sim_Theta_file, n_col.Sim_Theta, ns, true)
+Sim_Temp_names = num2str(Ztot');
+Sim_Temp_names=cellstr(Sim_Temp_names);
+Sim_Temp_units ={'oC'};
+write_output(Sim_Temp_names, Sim_Temp_units, ...
+fnames.Sim_Temp_file, n_col.Sim_Temp, ns, true)
+%% spectrum (added on 19 September 2008)
+
+spectrum_hemis_optical_names = {'hemispherically integrated radiation spectrum'};
+spectrum_hemis_optical_units = {'W m-2 um-1'};
+write_output(spectrum_hemis_optical_names, spectrum_hemis_optical_units, fnames.spectrum_hemis_optical_file, n_col.spectrum_hemis_optical, ns, true)
+
+spectrum_obsdir_optical_names = {'radiance spectrum in observation direction'};
+spectrum_obsdir_optical_units = {'W m-2 sr-1 um-1'};
+write_output(spectrum_obsdir_optical_names, spectrum_obsdir_optical_units, fnames.spectrum_obsdir_optical_file, n_col.spectrum_obsdir_optical, ns, true)
+
+if options.calc_ebal
+    write_output({'thermal BlackBody emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ...
+        fnames.spectrum_obsdir_BlackBody_file, n_col.spectrum_obsdir_BlackBody, ns, true)
+     if options.calc_planck
+         write_output({'thermal emission spectrum in hemispherical direction'}, {'W m-2 sr-1 um-1'}, ...
+         fnames.spectrum_hemis_thermal_file, n_col.spectrum_hemis_thermal, ns, true)
+         write_output({'thermal emission spectrum in observation direction'}, {'W m-2 sr-1 um-1'}, ...
+         fnames.spectrum_obsdir_thermal_file, n_col.spectrum_obsdir_thermal, ns, true)
+     end
+end
+         write_output({'irradiance'}, {'W m-2 um-1'}, ...
+         fnames.irradiance_spectra_file, n_col.irradiance_spectra, ns, true)
+         write_output({'reflectance'}, {'fraction of radiation in observation direction *pi / irradiance'}, ...
+         fnames.reflectance_file, n_col.reflectance, ns, true)
+%% input and parameter values (added June 2012)
+%write_output(fnames.pars_and_input_file, true)
+%write_output(fnames.pars_and_input_short_file, true)
+%% Optional Output
+if options.calc_vert_profiles
+    write_output({'Fraction leaves in the sun, fraction of observed, fraction of observed&visible per layer'}, {'rows: simulations or time steps, columns: layer numbers'}, ...
+    fnames.gap_file, n_col.gap, ns, true)
+    write_output({'aPAR per leaf layer'}, {'umol m-2 s-1'}, ...
+    fnames.layer_aPAR_file, n_col.layer_aPAR, ns, true)
+    write_output({'aPAR by Cab per leaf layer'}, {'umol m-2 s-1'}, ...
+    fnames.layer_aPAR_Cab_file, n_col.layer_aPAR_Cab, ns, true)
+    if options.calc_ebal
+    write_output({'leaf temperature of sunlit leaves, shaded leaves, and weighted average leaf temperature per layer'}, {'^oC ^oC ^oC'}, ...
+    fnames.leaftemp_file, n_col.leaftemp, ns, true)
+    write_output({'sensible heat flux per layer'}, {'W m-2'}, ...
+    fnames.layer_H_file, n_col.layer_H, ns, true)
+    write_output({'latent heat flux per layer'}, {'W m-2'}, ...
+    fnames.layer_LE_file, n_col.layer_LE, ns, true)
+    write_output({'photosynthesis per layer'}, {'umol m-2 s-1'}, ...
+    fnames.layer_A_file, n_col.layer_A, ns, true)
+    write_output({'average NPQ = 1-(fm-fo)/(fm0-fo0), per layer'}, {''}, ...
+    fnames.layer_NPQ_file, n_col.layer_NPQ, ns, true)
+    write_output({'net radiation per leaf layer'}, {'W m-2'}, ...
+    fnames.layer_Rn_file, n_col.layer_Rn, ns, true)
+    end
+if options.calc_fluor
+fluorescence_names = {'supward fluorescence per layer'};
+fluorescence_units = {'W m-2'};
+write_output(fluorescence_names, fluorescence_units, fnames.layer_fluorescence_file, n_col.layer_fluorescence, ns, true)
+end
+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)
+    if options.calc_PSI
+        write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution, for PSI only'}, {'W m-2 um-1 sr-1'}, ...
+        fnames.fluorescencePSI_file, n_col.fluorescencePSI, ns, true)
+        write_output({'fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution, for PSII only'}, {'W m-2 um-1 sr-1'}, ...
+        fnames.fluorescencePSII_file, n_col.fluorescencePSII, ns, true)
+    end
+    write_output({'hemispherically integrated fluorescence per simulation for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ...
+    fnames.fluorescence_hemis_file, n_col.fluorescence_hemis, ns, true)
+    write_output({'total emitted fluorescence by all leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ...
+    fnames.fluorescence_emitted_by_all_leaves_file, n_col.fluorescence_emitted_by_all_leaves, ns, true)
+    write_output({'total emitted fluorescence by all photosystems for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1'}, ...
+    fnames.fluorescence_emitted_by_all_photosystems_file, n_col.fluorescence_emitted_by_all_photosystems, ns, true)
+    write_output({'TOC fluorescence contribution from sunlit leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
+    fnames.fluorescence_sunlit_file, n_col.fluorescence_sunlit, ns, true)
+    write_output({'TOC fluorescence contribution from shaded leaves for wavelengths of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
+    fnames.fluorescence_shaded_file, n_col.fluorescence_shaded, ns, true)
+    write_output({'TOC fluorescence contribution from from leaves and soil after scattering for wavelenghts of 640 to 850 nm, with 1 nm resolution'}, {'W m-2 um-1 sr-1'}, ...
+    fnames.fluorescence_scattered_file, n_col.fluorescence_scattered, ns, true)
+end
+write_output({'Bottom of canopy irradiance in the shaded fraction, and average BOC irradiance'}, {'First 2162 columns: shaded fraction. Last 2162 columns: average BOC irradiance. Unit: Wm-2 um-1'}, ...
+fnames.BOC_irradiance_file, n_col.BOC_irradiance, ns, true)
+
+fclose('all');
+
+%% deleting .bin
+structfun(@delete, fnames)
+end
+
+function write_output(header, units, bin_path, f_n_col, ns, not_header)
+    if nargin == 5
+        not_header = false;
+    end
+    n_csv = strrep(bin_path, '.bin', '.csv');
+
+    f_csv = fopen(n_csv, 'w');
+    header_str = [strjoin(header, ','), '\n'];
+    if not_header
+        header_str = ['#' header_str];
+    else
+        % it is a header => each column must have one
+        assert(length(header) == f_n_col, 'Less headers than lines `%s` or n_col is wrong', bin_path)
+    end
+    fprintf(f_csv, header_str);
+    fprintf(f_csv, ['#' strjoin(units, ','), '\n']);
+
+    f_bin = fopen(bin_path, 'r');
+    out = fread(f_bin, 'double');
+%     fclose(f_bin);  % + some useconds to execution
+
+    out_2d = reshape(out, f_n_col, ns)';
+%     dlmwrite(n_csv, out_2d, '-append', 'precision', '%d'); % SLOW!
+    for k=1:ns
+        fprintf(f_csv, '%d,', out_2d(k, 1:end-1));
+        fprintf(f_csv, '%d\n', out_2d(k, end));  % saves from extra comma
+    end
+%     fclose(f_csv);
+end

src/+io/create_output_files.m

@@ -33,7 +33,7 @@
 
 %% Log File
 for i = 1:length(parameter_file)  
-    copyfile(['../' parameter_file{i}],[Output_dir,'Parameters/', parameter_file{i}],'f')
+    copyfile(['../input/' parameter_file{i}],[Output_dir,'Parameters/', parameter_file{i}],'f')
 end
 fidpath          = fopen([Output_dir,'Parameters/SCOPEversion.txt'],'w');      % complete path of the SCOPE code
 fprintf(fidpath,'%s', path_of_code);

src/+io/create_output_files_binary.m

@@ -0,0 +1,114 @@
+function [Output_dir, f, fnames] = create_output_files_binary(parameter_file, sitename, path_of_code, input_path, spectral, options)
+%% Create Output dir
+string          = clock;
+simulation_name = char(sitename);
+outdir_name     = sprintf('%s_%4.0f-%02.0f-%02.0f-%02.0f%02.0f', simulation_name, string(1:5));
+Output_dir = [fullfile('../output/', outdir_name) filesep];
+
+warning('off','MATLAB:DELETE:FileNotFound')
+if any(~exist(Output_dir,'dir'))
+    mkdir(Output_dir)
+    mkdir([Output_dir,'Parameters' filesep])
+end
+
+%% Log File
+for i = 1:length(parameter_file)
+    copy_name = [strrep(parameter_file{i}, '.csv', '') '_' outdir_name '.csv'];
+    copyfile([input_path parameter_file{i}],[Output_dir,'Parameters/', copy_name],'f')
+end
+fidpath          = fopen([Output_dir,'Parameters/SCOPEversion.txt'],'w');      % complete path of the SCOPE code
+fprintf(fidpath,'%s', path_of_code);
+
+%% Filenames, will become .csv if options is on
+fnames.flu_file                = fullfile(Output_dir,'fluxes.bin');
+fnames.surftemp_file           = fullfile(Output_dir,'surftemp.bin');
+fnames.aerodyn_file       = fullfile(Output_dir,'aerodyn.bin');
+fnames.radiation_file          = fullfile(Output_dir,'radiation.bin');
+fnames.fluorescence_file       = fullfile(Output_dir,'fluorescence.bin');
+fnames.wl_file       = fullfile(Output_dir,'wl.bin');              % wavelength
+fnames.irradiance_spectra_file        = fullfile(Output_dir,'irradiance_spectra.bin');     % Fluorescence
+fnames.spectrum_hemis_optical_file            = fullfile(Output_dir,'spectrum_hemis.bin');
+fnames.spectrum_obsdir_optical_file           = fullfile(Output_dir,'spectrum_obsdir.bin');
+fnames.reflectance_file           = fullfile(Output_dir,'reflectance.bin');   % reflectance spectrum
+fnames.BOC_irradiance_file           = fullfile(Output_dir,'BOC_irradiance.bin');   % reflectance spectrum
+fnames.Sim_Theta_file           = fullfile(Output_dir,'Sim_Theta.bin');   % soil moisture
+fnames.Sim_Temp_file           = fullfile(Output_dir,'Sim_Temp.bin');   % soil temperature
+if options.calc_ebal
+    fnames.spectrum_obsdir_BlackBody_file           = fullfile(Output_dir,'spectrum_obsdir_BlackBody.bin');  % spectrum observation direction
+end
+
+%if ~(options.simulation==1)
+fnames.pars_and_input_file           = fullfile(Output_dir,'pars_and_input.bin');      % wavelength
+
+%for j = 1:length(V) 
+ %   fprintf(fidv,'%s\t',V(j).Name);
+%end
+%fprintf(fidv,'\r');
+%end
+
+%if ~(options.simulation==1)
+fnames.pars_and_input_short_file           = fullfile(Output_dir,'pars_and_input_short.bin');      % wavelength
+%for j = find(vmax>1)
+  %  fprintf(fidvs,'%s\t',V(vmax>1).Name);
+%end
+%fprintf(fidvs,' \r');
+%
+%% Optional Output
+if options.calc_vert_profiles
+    fnames.gap_file            = fullfile(Output_dir,'gap.bin');              % gap
+    fnames.leaftemp_file        = fullfile(Output_dir,'leaftemp.bin');        % leaftemp
+    fnames.layer_H_file         = fullfile(Output_dir,'layer_H.bin');      % vertical profile
+    fnames.layer_LE_file        = fullfile(Output_dir,'layer_lE.bin');         % latent heat
+    fnames.layer_A_file         = fullfile(Output_dir,'layer_A.bin');        %
+    fnames.layer_aPAR_file      = fullfile(Output_dir,'layer_aPAR.bin');    %
+    fnames.layer_aPAR_Cab_file  = fullfile(Output_dir,'layer_aPAR_Cab.bin');     %
+    fnames.layer_Rn_file        = fullfile(Output_dir,'layer_Rn.bin');        %
+    if options.calc_fluor
+        fnames.layer_fluorescence_file = fullfile(Output_dir,'layer_fluorescence.bin');        %
+        fnames.layer_fluorescenceEm_file = fullfile(Output_dir,'layer_fluorescenceEm.bin');        %
+        fnames.layer_NPQ_file = fullfile(Output_dir,'layer_NPQ.bin');        %
+    end
+
+else
+    delete([Output_dir,'../output/leaftemp.bin'])
+    delete([Output_dir,'../output/layer_H.bin'])
+    delete([Output_dir,'../output/layer_lE.bin'])
+    delete([Output_dir,'../output/layer_A.bin'])
+    delete([Output_dir,'../output/layer_aPAR.bin'])
+    delete([Output_dir,'../output/layer_Rn.bin'])
+end
+
+if options.calc_fluor
+    fnames.fluorescence_file              = fullfile(Output_dir,'fluorescence.bin');      % Fluorescence
+    if options.calc_PSI
+        fnames.fluorescencePSI_file       = fullfile(Output_dir,'fluorescencePSI.bin');     % Fluorescence
+        fnames.fluorescencePSII_file      = fullfile(Output_dir,'fluorescencePSII.bin');     % Fluorescence
+    end
+    fnames.fluorescence_hemis_file        = fullfile(Output_dir,'fluorescence_hemis.bin');     % Fluorescence
+    fnames.fluorescence_emitted_by_all_leaves_file         = fullfile(Output_dir,'fluorescence_emitted_by_all_leaves.bin');
+    fnames.fluorescence_emitted_by_all_photosystems_file   = fullfile(Output_dir,'fluorescence_emitted_by_all_photosystems.bin');
+    fnames.fluorescence_sunlit_file       = fullfile(Output_dir,'fluorescence_sunlit.bin');     % Fluorescence
+    fnames.fluorescence_shaded_file       = fullfile(Output_dir,'fluorescence_shaded.bin');     % Fluorescence
+    fnames.fluorescence_scattered_file    = fullfile(Output_dir,'fluorescence_scattered.bin');      % Fluorescence
+else
+    delete([Output_dir,'fluorescence.bin'])
+end
+
+if options.calc_directional
+    delete([Output_dir,'BRDF/*.bin'])
+end
+
+if options.calc_planck && options.calc_ebal
+    fnames.spectrum_obsdir_thermal_file    = fullfile(Output_dir,'spectrum_obsdir_thermal.bin');  % spectrum observation direction
+    fnames.spectrum_hemis_thermal_file     = fullfile(Output_dir,'spectrum_hemis_thermal.bin');   % spectrum hemispherically integrated
+end
+%% Open files for writing
+f = structfun(@(x) fopen(x, 'w'), fnames, 'UniformOutput',false);
+
+%% write wl
+wlS = spectral.wlS; %#ok<*NASGU>
+wlF = spectral.wlF;
+
+save([Output_dir 'wlS.txt'], 'wlS', '-ascii');
+save([Output_dir 'wlF.txt'], 'wlF', '-ascii');
+end

src/+io/initialize_output_structures.m

@@ -1,7 +1,6 @@
 function [rad,thermal,fluxes] = initialize_output_structures(spectral)
 
-[iter.counter                                        ,...
-    fluxes.Rntot,     fluxes.lEtot,      fluxes.Htot,   fluxes.Atot            ,...
+[    fluxes.Rntot,     fluxes.lEtot,      fluxes.Htot,   fluxes.Atot            ,...
     fluxes.Rnctot,    fluxes.lEctot,     fluxes.Hctot,  fluxes.Actot           ,...
     fluxes.Rnstot,    fluxes.lEstot,     fluxes.Hstot,  fluxes.Gtot,   fluxes.Resp    ,...
     thermal.Tcave,    thermal.Tsave   ,...

src/+io/load_timeseries.m

@@ -1,6 +1,6 @@
 function [V,xyt,canopy]  = load_timeseries(V,leafbio,soil,canopy,meteo,constants,F,xyt,path_input,options)
 
-Dataset_dir         = ['dataset ' char(F(5).FileName)];
+Dataset_dir         = '';
 t_file              = char(F(6).FileName);
 year_file           = char(F(7).FileName);
 Rin_file            = char(F(8).FileName);
@@ -78,7 +78,7 @@
 end
 if  ~isempty(LAI_file)
     LAItable        = load([path_input,Dataset_dir,'/',LAI_file]);
-    V(22).Val         = interp1(LAItable(:,1),LAItable(:,2),t_);
+    V(22).Val       = LAItable(:,2);
 else
     V(22).Val          = canopy.LAI*ones(size(time_));
 end