clc
clear all
close all

%-----------------------------------------
% EXPERIMENT DEFINITION
%-----------------------------------------

storm_name = 'Idai';
save_path = ['/discover/nobackup/jkolassa/SMAP_TC_project/Data_analysis/' storm_name '/'];

exp_dir = '/discover/nobackup/projects/smaptc/jkolassa/';
exp_name = 'Idai_LDAS_DA_4D_C360';
exp_resolution = 'CF0360x6C_GLOBAL';
da_interval = 3;   % DA interval in hours

lat_range = [-90 90];
lon_range = [-180 180];

start_year = 2019;
start_month = 2;
start_day = 4;

end_year = 2019;
end_month = 3;
end_day = 22;

%---------------------------------------------
% DATA READ PREPARATION
%---------------------------------------------

leapYears = 1980:4:2022;
monthLength_nly = [31 28 31 30 31 30 31 31 30 31 30 31];
monthLength_ly = [31 29 31 30 31 30 31 31 30 31 30 31];
mat_index = 0;
tot_index = 0;

%---------------------------------------------
% DATA READ
%---------------------------------------------

% loop over years
for y = start_year:end_year

    y_str = num2str(y,'%04i');

    % define month length based on leap year status
    if sum(y==leapYears)>0
        monthLength = monthLength_ly;
    else
        monthLength = monthLength_nly;
    end

    % loop over months
    for m = start_month:end_month

        m_str = num2str(m,'%02i');

        % define start and end day for current month
        day_first = 1;
        day_last = monthLength(m);

        if (m == start_month)
           day_first = start_day;
        end

        if (m == end_month)
           day_last = end_day;
        end

        % loop over days
        for d = day_first:day_last

            mat_index = mat_index + 1;

            d_str = num2str(d,'%02i');
            t_index = 0;

            % loop over time steps
            for t = 1:3:22
                t_index = t_index + 1;
                tot_index = tot_index + 1;

                t_str = num2str(t,'%02i');

                % define increment file name

                exp_incr_file = [exp_dir exp_name '/output/' exp_resolution '/cat/ens_avg/Y' y_str '/M' m_str '/' exp_name ... 
                                '.inst3_2d_lndfcstana_Nx.' y_str m_str d_str '_' t_str '30z.nc4'];

                % on first time step read coordinate information and initialize matrices

                if (mat_index == 1) && (t == 1)

                   lat_vec = double(ncread(exp_incr_file,'lats'));
                   lon_vec = double(ncread(exp_incr_file,'lons'));

                   datenum_vec = NaN*ones(sum(monthLength(start_month:end_month)),1);

                end

                ana_data = double(ncread(exp_incr_file,'SFMC_ANA'));
                fcst_data = double(ncread(exp_incr_file,'SFMC_FCST'));

                AmF = ana_data - fcst_data;

                AmF_mean_vec(tot_index) = nanmean(AmF(:));
                fcst_mean_vec(tot_index) = nanmean(fcst_data(:));
            end % t

            % log date

            datenum_vec(mat_index,1) = datenum(y,m,d);

            % compute daily total increments

         end % d 
     end % m
 end % y

