Part_4_Ex_2.m

%% Ensemble Kalman Smoother (Zaitchik et al., 2008)
% So far, we have been adding the monthly update at the end of each month. Different 
% techniques exist to distribute the monthly update within the month (see, e.g., 
% Zaitchik et al., 2008; Girotto et al., 2016; Tian et al. 2017). In this exercise, 
% we implement the Zaitchik smoother.
%% Setting up
% We will use the same setting we have been using so far for DA. We will also 
% reduce the processing to 2007-2008.

clc, clear all;
addpath(genpath('./Codes/')); dbstop if error;
load("Config_DA.mat");
config.run.todate = '2008-12-31';
%% 
% Let us create a new folder to save the outputs of this run. *Copy-paste the 
% missing code from previous exercises.*

config.states.directory = 'Model_States_Output/03_DA_States_Smoother/ens%n/';
config.summary.directory = 'Model_States_Output/03_DA_States_Smoother/Summary/';
% Make new directory and copy initialization files
%%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%

%% DA run
% First we need to separate the ensemble run in monthly chunks. *Copy-paste 
% the missing code from previous exercises.*

% Define start and end dates
%%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%
% Create array of epochs were DA should be performed
%%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%
%% 
% Let us now iterate over the months:

for i = 1:(length(t_array)-1)
    fprintf('Processing month number %d out of %d\n',i,length(t_array)-1)
%% 
% First we will run an ensemble run for this month. *Copy-paste the missing 
% code from previous exercises.*

    % Set beginning and ending time for monthly run
    %%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%
    % Perform run for each ensemble member
    %%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%
%% 
% We now run the DA step. This is also similar to the regular EnKF update.  
% *Copy-paste the missing code from previous exercises, and omit the last step 
% (|Main06|).*

    % DA update step (OMIT THE LAST STEP, Main06).
    %%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%
    
%     Main06_assignUpdateEndOfMonth(year(t),month(t),config_ens);
%% 
% Now comes the novelty: instead of adding the update to the last day of the 
% month (this was what |Main06| was doing), we will rewind to the beginning of 
% the month. We will then follow this algorithm:
%% 
% # Apply a fraction of the update to first day of the month.
% # Run the model one day ahead.
% # Apply a fraction of the update to second day of the month.
% # Run the model one day ahead.
% # ...
% # Apply a fraction of the update to the last day of the month.
%% 
% More specifically, each day of the month will receive an update of |U/n_days|, 
% where |U| is the total update computed during the DA step, and |n_days| is the 
% amount of days in this month. At the end of the month, the updates applied in 
% all of the days should accumulate, resulting in an update approximately equivalent 
% to |U| in the last day.

    % Smooth update application
    t1 = t_array(i) + caldays(1); % First day of the month;
    % Iterate over day of the month
    for dd = t1:t
        % Setup config to perform a 1-day run
        config_day = config_m;
        config_day.run.fromdate = datestr(dd,'yyyy-mm-dd');
        config_day.run.todate = datestr(dd,'yyyy-mm-dd');
        % Perform run for each ensemble member
        %%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%
%% 
% The daily update process is carried out by the function |Main08_apply_daily_increment|.

        % Apply increment fraction over that day
        Main08_apply_daily_increment(year(dd),month(dd),day(dd),config_ens,'fraction');
    end   
end
%% Result extraction and visualization

F02_ensemble_extract_state_results_W3RA(config);
F_visualizeEx7(config)
%% Reflection questions
%% 
% # How do the TWS dynamics after DA obtained by the EnKS compared to the TWS 
% dynamics of the OL simulation and the DA results obtained by the regular EnKF? 
%% 
%% 
% # What can be the benefits and drawbacks of running an EnKS with respect to 
% an EnKF? 
%%