Upload New File

init_virtual.m

+%
+% Copyright 2019-2020 LIST (Luxembourg Institute of Science and Technology), all right reserved.
+%
+% Author: Nicolas Rodriguez (nicolas.bjorn.rodriguez@gmail.com)
+%
+function [meas,model_input,param,options]=init_virtual(options)
+
+if options.display==1
+    disp('Initializing the model...');
+    tic
+end
+    
+    SASQ_p1=NaN;
+    SASQ_p2=NaN;
+    SASQ_p3=NaN;
+    SASQ_p4=NaN;
+    SASQ_p5=NaN;
+    SASQ_p6=NaN;
+    SASQ_p7=NaN;
+    SASQ_p8=NaN;
+    SASQ_p9=NaN;
+    SASET_p1=NaN;
+    SASET_p2=NaN;
+    SASET_p3=NaN;
+    SASET_p4=NaN;
+    SASET_p5=NaN;
+    SASET_p6=NaN;
+    SASET_p7=NaN;
+    SASET_p8=NaN;
+    SASET_p9=NaN;
+    
+    S_i=1000; % mm
+    S_ET=850; % mm
+    
+    Psi=1; % no unit
+    m=20; % no unit
+   
+    T0 = 1.7; % years
+    n=1e-2; % no unit
+    
+    switch options.tracer
+        case 1
+            C_i=-50; % per mil
+    end
+    
+% assign parameters to variables
+    % fixed parameters
+        for k=1:options.param.fix.nb
+            temp=options.param.fix.name{k};
+            options.param.fix.val(k)=eval(temp);
+        end
+    % calibration parameters
+        for k=1:options.param.cal.nb
+            temp=options.param.cal.name{k};
+            param(k)=eval(temp);
+        end
+        
+% Load data
+   
+t=datenum('01-Oct-1917'):1:datenum('01-Oct-2017'); % time vector, 100 years of daily points
+load('J_virtual.mat'); % precipitation in mm per time step
+
+% randomly remove the many small values in J with a probabibilty of 80%
+sum_lost=0;
+for k=1:length(J)
+    if J(k)<=1 % values small than 1 mm/d
+        r=rand();
+        if r<=0.8
+            sum_lost=sum_lost+J(k); % count how much precip is removed
+            J(k)=0; % remove the corresponding precip amounts
+        end
+    end
+end
+gain=sum_lost/sum(J~=0); % then distribute the amount lost to all other days with precipitation
+for k=1:length(J)
+    if J(k)~=0
+        J(k)=J(k)+gain;
+    end
+end
+% and add the large precipitation event corresponding to the pesticide injection, this corresponds to the date '01-Oct-2009 00:00'
+J(end-8*365.25)=100; % mm
+
+% Q, PET, isotope in precip and streamflow
+load('Q_virtual.mat');
+load('ET_virtual.mat');
+load('Cp_virtual.mat');
+load('C_Q.mat');
+    
+% assignment of various model run periods
+options.periods.all=1:length(t);
+options.periods.warmup=[];
+options.periods.run=1:length(t);
+
+% save the data in the input structure
+model_input.t=t;
+model_input.tfull=t;
+model_input.J=J;
+model_input.PET=ET/24; % mm/d -> mm/h (ET*dt needs to be unitless and dt is in hours)
+
+model_input.Cp=Cp;
+model_input.Cpfull=Cp;
+model_input.size_true_warmup=0;
+model_input.Q=Q/24; % mm/d -> mm/h (Q*dt needs to be unitless and dt is in hours)
+
+meas.Cq=C_Q;
+
+if options.display==1
+    toc()
+end
+
+end
\ No newline at end of file