clc
clear all
close all

% make selections
start_date = '20150401';
stop_date = '20201231';
lat_range = [0 33];      % [southern_border northern_border];
lon_range = [-110 -70];  % [western_border eastern_border]

% convert selections
smap_l4_path = '/css/smapl4/public/L4_Products/L4_SM/';

start_year = str2double(start_date(1:4));
start_month = str2double(start_date(5:6));
start_day = str2double(start_date(7:8));
stop_year = str2double(stop_date(1:4));
stop_month = str2double(stop_date(5:6));
stop_day = str2double(stop_date(7:8));

total_months = (stop_year-start_year+1)*12;
total_time_steps = total_months*31*8;

t_vec = {'0000' '0300' '0600' '0900' '1200' '1500' '1800' '2100'};
month_length_nly = [31 28 31 30 31 30 31 31 30 31 30 31];
month_length_ly = [31 29 31 30 31 30 31 31 30 31 30 31];
leap_years = 2012:4:2024;

stream_datenum_1 = datenum(2020,8,27,0,0,0);
mat_index = 0;

for y = start_year:stop_year

    y_str = num2str(y);
    if (sum(y==leap_years)>0)                                                                                            
        month_length = month_length_ly;                                                                                 
    else                                                                                                              
        month_length = month_length_nly;                                                                                
    end 

    m_first = 1;
    m_last = 12;
    if (y==start_year)
        m_first = start_month;
    end

    if (y==stop_year)
        m_last = stop_month;
    end

    for m = m_first:m_last

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

        day_first = 1;
        day_last = month_length(m);
        if (y==start_year)&(m==start_month)
           day_first = start_day;
        end
        if (y==stop_year)&(m==stop_month)
            day_last = stop_day;
        end

        for d = day_first:day_last

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

            for t = 1:length(t_vec)

                curr_datenum = datenum(y,m,d,str2double(t_vec{t}(1:2)),0,0);
                if (curr_datenum<=stream_datenum_1)
                   stream_str = 'Vv4030';
                else
                   stream_str = 'Vv5014';
                end

                mat_index = mat_index + 1;

                smap_l4_file = [smap_l4_path stream_str '/aup/Y' y_str '/M' m_str '/D' d_str ...
                           '/SMAP_L4_SM_aup_' y_str m_str d_str 'T' t_vec{t} '00_' stream_str '_001.h5'];

                if (mat_index==1)
                   lat_data = double(h5read(smap_l4_file,'/cell_lat'));
                   lon_data = double(h5read(smap_l4_file,'/cell_lon'));
                   lat_index = (lat_data(1,:)>=lat_range(1)) + (lat_data(1,:)<=lat_range(2));   
                   lat_index = lat_index == 2;                                    
                   lon_index = (lon_data(:,1)>=lon_range(1)) + (lon_data(:,1)<=lon_range(2));   
                   lon_index = lon_index == 2;                                    
                   domain_lats = lat_data(lon_index,lat_index);                   
                   domain_lons = lon_data(lon_index,lat_index);  
                   sm_incr_mat = NaN*ones(size(domain_lats,1),size(domain_lats,2),total_time_steps);
                end
                sm_ana_data = double(h5read(smap_l4_file,'/Analysis_Data/sm_surface_analysis'));
                sm_fcst_data = double(h5read(smap_l4_file,'/Forecast_Data/sm_surface_forecast'));

                sm_incr_mat(:,:,mat_index) = sm_ana_data(lon_index,lat_index) - sm_fcst_data(lon_index,lat_index);
                doy = sum(month_length_ly(1:(m-1))+d);
                time_vec(mat_index,1) = doy;
                time_vec(mat_index,2) = str2double(t_vec{t});

            end

        end

    end
end


