Part_2_Ex_3.m

%% DA in the Murray-Darling basin
% We now perform DA in the Murray-Darling basin. This river basin was affected 
% by a very severe drought (Millennium Drought) between 2001 and 2009, followed 
% by a notably wet period in 2010-2013 (Forootan et al., 2012; Schumacher et al., 
% 2018).
%% Setting up DA configuration

clc, clear all;
addpath(genpath('./Codes/')); dbstop if error;
load("Config_base.mat");
%% 
% We first need to add some data to the base configuration file, so that the 
% scripts can find the data necessary to perform DA. We need to turn on the DA 
% switch:

config.run.DA = 1;
% config.run.todate = '2007-12-31';
%% 
% We need to provide a directory to the observations, a folder to store the 
% xMinus and xPlus files, and the directory to the design matrix.

config.DAOptions.TWSObservations = 'Data/03_Perturbed_Observations/GRACE_TWS_Observations_Error_2007_2013.mat';
config.DAOptions.DAGeneratedFiles.xMinus = 'Model_States_Output/03_DA_States/xMinus/';
config.DAOptions.DAGeneratedFiles.xPlus = 'Model_States_Output/03_DA_States/xPlus/';
config.DAOptions.OtherData.DesignMatrix = 'Data/03_DA_Files_50km/A_design_matrix/A.mat';
%% 
% For the moment, we will use no localization (this will be used in a later 
% exercise).

config.DAOptions.Localization = 0;
%% 
% Additionally, we will save the states of the DA run on a different folder. 
% We therefore also need to initialize this folder. *Please copy-paste the missing 
% code from exercise 0.*

config.states.directory = 'Model_States_Output/03_DA_States/ens%n/';
config.summary.directory = 'Model_States_Output/03_DA_States/Summary/';
config.summary.nameKeyword = 'DA'; % Keyword for naming summary file
% Make new directory and copy initialization files
%%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%

%% 
% Before we continue, let us save the DA configuration so that we can use it 
% in other exercises.

save('Config_DA.mat','config');
%% DA run
% As we will be performing the DA update step at the end of each month, we first 
% need to separate the ensemble run in monthly chunks.

% Define start and end dates
start_date = datetime(config.run.fromdate);
end_date = datetime(config.run.todate);
% Create array of epochs were DA should be performed
t_array = datetime(year(start_date),month(start_date),1):calmonths(1):datetime(year(end_date),month(end_date),1); % List of months to process.
t_array = t_array + calmonths(1) - caldays(1); % Set list to last day of the month.
t_array = [(start_date - caldays(1)), t_array]; % Add first day.
disp(length(t_array)-1);
%% 
% As we see, our processing period has been divided in 84 1-month chunks. Let 
% us now iterate over these 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. *Please copy-paste the missing 
% code from exercise 1.*

    % Set beginning and ending time for monthly run
    config_m = config;
    config_m.run.fromdate = datestr(t_array(i)+caldays(1),'yyyy-mm-dd');
    config_m.run.todate = datestr(t_array(i+1),'yyyy-mm-dd');
    % Perform run for each ensemble member
    nE = config.ensembleOptions.nE; % Number of ensemble members.
    for j=1:nE
        config_ens = config_m;
        % Assign proper directories to each ensemble member:
        %%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%
        % Run ensemble member.
        %%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%
    end
%% 
% After the monthly run, we will perform the Data Assimilation step. For this 
% purpose, we will use the following functions:
%% 
% * |Main03_computeMonthlyAverage_xminus| retrieves all the model states of 
% the previous months and computes a monthly average for each ensemble. The monthly-averaged 
% states are saved in matrix xMinus.
% * |Main05_performDA| performs the DA step.
% * |Main06_assignUpdateEndOfMonth| adds the monthly update to the model estimates 
% at the last day of the month. This updated state will be used to initialize 
% the next monthly run.

    % DA update step
    t = t_array(i+1); % Set t to the last day of the month
%% Monthly EnKF

    Main03_computeMonthlyAverage_xminus(year(t),month(t),config_ens);
    Main05_performDA(year(t),month(t),config_ens);
    Main06_assignUpdateEndOfMonth(year(t),month(t),config_ens);
end
%% 
% This process is repeated for each month, until we reach the end of the processing 
% period.
%% Result extraction
% Now that the processing is finished, we can compute the ensemble-averaged 
% results. *Please copy-paste the missing code from exercise 1.*

%%%%%%%%%%% COPY-PASTED FROM EX. 1 %%%%%%%%%%%%%%%
%%%%%%%%%%% TO BE FILLED %%%%%%%%%%%%%%%

%% Visualization

F_visualizeEx3(config);
%% Reflection questions
%% 
% # How did the model TWS time series change after DA? Did the representation 
% of the climate (inter-annual) variability change? 
%% 
%% 
% # Which kind of science questions do you think these results can help us to 
% answer? 
%% 
%% 
% # Were all the grids within each sub-basin updated homogeneously? Why? 
%% 
%% 
% # What would have happened if we had assimilated the GRACE fields resampled 
% at 0.5 degrees? 
%%