clc
clear all
close all

exp_name = 'IDS_fcst_CTRL_M09_2007_ens1';
exp_path = ['/discover/nobackup/elee15/GEOSldas_4_0/sims/IDS/fcst_ctrl/' exp_name '/output/SMAP_EASEv2_M09_GLOBAL/cat/ens0000/'];
start_year = 2007;
stop_year = 2007;

% read first file for lat/lon data

data_file = [exp_path 'Y' num2str(start_year) '/M01/' exp_name '.tavg24_1d_lnd_Nt.monthly.' num2str(start_year) '01.nc4'];
lat_data = double(ncread(data_file,'lat'));
lon_data = double(ncread(data_file,'lon'));

% compute indices on 0.5 degree lat/lon grid

R_earth = 6371;

% compute threshold latitude outside of which 36km > 0.5deg
delta_lon_thresh = 0.5/(180/pi);
%lat_thresh = acos(36/(delta_lon_thresh*R_earth));
lat_thresh = acos(9/(delta_lon_thresh*R_earth));
lat_thresh = lat_thresh*180/pi;

lat_data_rad = lat_data./(180/pi);
lon_data_rad = lon_data./(180/pi);

lat_ll_deg = -89.75:0.5:89.75;
lat_ll_rad = lat_ll_deg./(180/pi);
lon_ll_deg = -179.75:0.5:179.75;
lon_ll_rad = lon_ll_deg./(180/pi);

im = 720;
jm = 360;

for n = 1:length(lat_data)

    if abs(lat_data(n))>lat_thresh
       lat_dist = (lat_ll_deg - lat_data(n));
       lat_index = abs(lat_dist) <= 0.5;
       jx{n} = find(lat_index==1);

       %delta_lon = 36/(R_earth*cos(lat_data_rad(n)));  % convert 36km into distance in deg lon
       delta_lon = 9/(R_earth*cos(lat_data_rad(n)));  % convert 9km into distance in deg lon
       delta_lon = delta_lon*180/pi;
       lon_dist = (lon_ll_deg-lon_data(n));
       lon_index = abs(lon_dist) <= delta_lon;
       ix{n} = find(lon_index==1);
     else
       ix{n} = floor((lon_data(n)+180)*(im/360)) + 1;
       jx{n} = floor((lat_data(n)+90)*(jm/180)) + 1;
     end

end

disp('entering year loop')
for y = start_year%:stop_year
    year_str = num2str(y);
    y

    %for m = 1:9
    for m = 1%:6
        if m < 10
           month_str = ['0' num2str(m)];
        else
           month_str = num2str(m);
        end

        data_file = [exp_path 'Y' year_str '/M' month_str '/' exp_name '.tavg24_1d_lnd_Nt.monthly.' year_str month_str '.nc4'];
        finfo = ncinfo(data_file);
        var_list = finfo.Variables;
        var_list(1:4) = [];
        
        regrid_mat = zeros(im,jm,length(var_list));
        
        for v = 1:length(var_list)
            
            var_name = var_list(v).Name;

            data_vec = double(ncread(data_file,var_name));
 
            data_mat = zeros(im,jm);
            count_mat = zeros(im,jm);

            for n = 1:length(lat_data)
                for i = 1:length(ix{n})
                    for j = 1:length(jx{n})                
                        data_mat(ix{n}(i),jx{n}(j)) = data_mat(ix{n}(i),jx{n}(j)) + data_vec(n);
                        count_mat(ix{n}(i),jx{n}(j)) = count_mat(ix{n}(i),jx{n}(j)) + 1;
                    end 
                end
            end

            data_mat_new = NaN*ones(size(data_mat));
            data_index = count_mat~=0;
            data_mat_new(data_index) = data_mat(data_index)./count_mat(data_index);

            regrid_mat(:,:,v) = data_mat_new;
        end
 
        save_file = [exp_path 'Y' year_str '/M' month_str '/monthly_0.5deg_lnd_' year_str month_str '.mat']
        save(save_file,'regrid_mat','var_list');
    end
end

