clear all
close all

%opengl software

dirbase = '/discover/nobackup/projects/gmao/smap/SMAP_Nature/';
addpath([dirbase, '/validation/tools'])
addpath  /home/qliu/projects/matlab_code/matlab_base/export_fig

calvalglobalvar(dirbase)

global COREMATCHUPDATAPATH
global PASTSATURDAYstr

use_SMAP_nature = 0;
compare_multiple_sat = 1;

myfontsize = 14;

Nmin = 480;  % mininum number for computing stat

if ~compare_multiple_sat
    if use_SMAP_nature
        CalVal_path = {[COREMATCHUPDATAPATH, '/SMAP_Nature_v05/']};
    else
        CalVal_path = {[COREMATCHUPDATAPATH, '/L4_Vb1010_aup/']};
    end
else
    CalVal_path1= [COREMATCHUPDATAPATH, '/L4_OL4001_aup/'];
    CalVal_path2= [COREMATCHUPDATAPATH, '/L4_Vv4030_aup/'];
    CalVal_path = {CalVal_path1, CalVal_path2};
end

startdate = '20150329';

enddate = '20190330';

%enddate = PASTSATURDAYstr;

start_year = str2num(startdate(1:4));

add_anomR = 1;

do_ST = 0;

if do_ST
    var_tag = {'SMSF'};
else
    var_tag = {'SMSF', 'SMRZ'};
end

QCflagyes = 1;

[col, row, scale, active, site, refpixii, description, refpixids, rootzone, sitealias, lat, lon] = ...
    get_refpix_master_info([], 'SM', 'active');

is_core = zeros(length(refpixids),1);
refpixlist = refpixids;

clear is_36 is_33 is_09 site_vec
for i=1:length(refpixids)
    is_36(i) = strcmp('36', refpixlist{i}(5:6));
    is_33(i) = strcmp('33', refpixlist{i}(5:6));
    is_09(i) = strcmp('09', refpixlist{i}(5:6));
    site_vec{i} = refpixlist{i}(1:4);
% beta release report list    
%     if strcmp(refpixlist{i}, '03013602') || ...
%             strcmp(refpixlist{i}, '03010908') || ...
%             strcmp(refpixlist{i}, '07013601') || ...
%             strcmp(refpixlist{i}, '07010902') || ...
%             strcmp(refpixlist{i}, '09013601') || ...
%             strcmp(refpixlist{i}, '09010906') || ...
%             strcmp(refpixlist{i}, '16013603') || ...
%             strcmp(refpixlist{i}, '16010906') || ...
%             strcmp(refpixlist{i}, '16010907') || ...
%             strcmp(refpixlist{i}, '16023602') || ...
%             strcmp(refpixlist{i}, '16020907') || ...
%             strcmp(refpixlist{i}, '16033602') || ...
%             strcmp(refpixlist{i}, '16030911') || ...
%             strcmp(refpixlist{i}, '16030916') || ...
%             strcmp(refpixlist{i}, '16043602') || ...
%             strcmp(refpixlist{i}, '16040901') || ...
%             strcmp(refpixlist{i}, '16060907') || ...
%             strcmp(refpixlist{i}, '16073602') || ...
%             strcmp(refpixlist{i}, '16070909') || ...
%             strcmp(refpixlist{i}, '16070910') || ...
%             strcmp(refpixlist{i}, '16070911') || ...
%             strcmp(refpixlist{i}, '27013601') || ...
%             strcmp(refpixlist{i}, '27010911') || ...
%             strcmp(refpixlist{i}, '41010906') || ...
%             strcmp(refpixlist{i}, '48013601') || ...
%             strcmp(refpixlist{i}, '48010902') || ...
%             strcmp(refpixlist{i}, '48010911')
% Validated report list
%      if strcmp(refpixlist{i}, '03013602') || ...
%             strcmp(refpixlist{i}, '03010903') || ...
%             strcmp(refpixlist{i}, '03010908') || ...
%             ...strcmp(refpixlist{i}, '04013602') || ...
%             ...strcmp(refpixlist{i}, '04010907') || ...
%             ...strcmp(refpixlist{i}, '04010910') || ...
%             strcmp(refpixlist{i}, '07013601') || ...
%             strcmp(refpixlist{i}, '07010902') || ...
%             strcmp(refpixlist{i}, '07010916') || ...
%             ...strcmp(refpixlist{i}, '07023602') || ...
%             strcmp(refpixlist{i}, '09013601') || ...
%             strcmp(refpixlist{i}, '09010906') || ...
%             ...strcmp(refpixlist{i}, '12043606') || ...
%             strcmp(refpixlist{i}, '16013603') || ...
%             strcmp(refpixlist{i}, '16010906') || ...
%             strcmp(refpixlist{i}, '16010907') || ...
%             strcmp(refpixlist{i}, '16010913') || ...
%             strcmp(refpixlist{i}, '16023602') || ...
%             ...strcmp(refpixlist{i}, '16020905') || ...
%             ...strcmp(refpixlist{i}, '16020906') || ...
%             strcmp(refpixlist{i}, '16020907') || ...
%             strcmp(refpixlist{i}, '16033602') || ...
%             strcmp(refpixlist{i}, '16030911') || ...
%             strcmp(refpixlist{i}, '16030916') || ...
%             strcmp(refpixlist{i}, '16043602') || ...
%             strcmp(refpixlist{i}, '16040901') || ...
%             ...strcmp(refpixlist{i}, '16040906') || ...
%             ...strcmp(refpixlist{i}, '16040911') || ...
%             ...strcmp(refpixlist{i}, '16060906') || ...
%             strcmp(refpixlist{i}, '16060907') || ...
%             strcmp(refpixlist{i}, '16073602') || ...
%             strcmp(refpixlist{i}, '16070909') || ...
%             strcmp(refpixlist{i}, '16070910') || ...
%             strcmp(refpixlist{i}, '16070911') || ...
%             ...strcmp(refpixlist{i}, '19023601') || ...
%             ...strcmp(refpixlist{i}, '19020911') || ...
%             strcmp(refpixlist{i}, '25013601') || ...
%             strcmp(refpixlist{i}, '25010911') || ...
%             strcmp(refpixlist{i}, '27013601') || ...
%             strcmp(refpixlist{i}, '27010910') || ...
%             strcmp(refpixlist{i}, '27010911') || ...
%             strcmp(refpixlist{i}, '41010906') || ...
%             strcmp(refpixlist{i}, '48013601') || ...
%             strcmp(refpixlist{i}, '48010902') || ...
%             strcmp(refpixlist{i}, '48010911')

% list prior Tv4000
%     if     strcmp(refpixlist{i}, '03013602') || ... 
%            strcmp(refpixlist{i}, '03010903') || ...
%            strcmp(refpixlist{i}, '03010908') || ...
%            strcmp(refpixlist{i}, '04013603') || ... % new 
%            strcmp(refpixlist{i}, '04010907') || ... % new ref info 
%            strcmp(refpixlist{i}, '04010910') || ... % new ref info
%            strcmp(refpixlist{i}, '07013601') || ...
%            strcmp(refpixlist{i}, '07010902') || ...
%            strcmp(refpixlist{i}, '07010916') || ...
%            ...strcmp(refpixlist{i}, '09013601') || ...
%            strcmp(refpixlist{i}, '09013610') || ...  % replace "09013601"
%            strcmp(refpixlist{i}, '09010906') || ...
%            strcmp(refpixlist{i}, '12033601') || ...   % new Tibet
%            ...strcmp(refpixlist{i}, '12043601') || ...    % new
%            strcmp(refpixlist{i}, '16013603') || ...
%            strcmp(refpixlist{i}, '16010906') || ...
%            strcmp(refpixlist{i}, '16010907') || ...
%            strcmp(refpixlist{i}, '16010913') || ...
%            strcmp(refpixlist{i}, '16023602') || ...
%            ...strcmp(refpixlist{i}, '16020905') || ...  % no data
%            ....strcmp(refpixlist{i}, '16020906') || ... % no data
%            strcmp(refpixlist{i}, '16020907') || ...
%            strcmp(refpixlist{i}, '16033602') || ...
%            strcmp(refpixlist{i}, '16030911') || ...
%            strcmp(refpixlist{i}, '16030916') || ...
%            strcmp(refpixlist{i}, '16043602') || ...
%            strcmp(refpixlist{i}, '16040901') || ...
%            ...strcmp(refpixlist{i}, '16040906') || ... % no data
%            ...strcmp(refpixlist{i}, '16040911') || ...
%            ...strcmp(refpixlist{i}, '16063602') || ... % in situ data quality
%            ...strcmp(refpixlist{i}, '16060906') || ... % new
%            strcmp(refpixlist{i}, '16060907') || ...
%            strcmp(refpixlist{i}, '16073602') || ...
%            strcmp(refpixlist{i}, '16070909') || ...
%            ...strcmp(refpixlist{i}, '16070910') || ... % one sensor data lost
%            strcmp(refpixlist{i}, '16070911') || ...
%            ...strcmp(refpixlist{i}, '17010907') || ... % new
%            strcmp(refpixlist{i}, '19023601') || ... % new 
%            ...strcmp(refpixlist{i}, '19020911') || ...  % new (no show)
%            strcmp(refpixlist{i}, '19020902') || ...  % subset of 19020911 
%            ...strcmp(refpixlist{i}, '24013602') || ...
%            strcmp(refpixlist{i}, '25013601') || ...
%            strcmp(refpixlist{i}, '25010911') || ...
%            strcmp(refpixlist{i}, '27013601') || ...
%            strcmp(refpixlist{i}, '27010910') || ...
%            strcmp(refpixlist{i}, '27010911') || ...
%            strcmp(refpixlist{i}, '41010906') || ...  %new ref info 36 km in situ the same
%            ...strcmp(refpixlist{i}, '44010911') || ...  % new 
%            strcmp(refpixlist{i}, '45013601') || ...  % new
%            strcmp(refpixlist{i}, '45010902') || ...  % new
%            strcmp(refpixlist{i}, '45023601') || ...  % new           
%            strcmp(refpixlist{i}, '45020902') || ...  % new           
%            strcmp(refpixlist{i}, '48013601') || ...
%            strcmp(refpixlist{i}, '48010902') || ...
%            strcmp(refpixlist{i}, '48010911') || ...
%            ...strcmp(refpixlist{i}, '53013601') || ...  %NEW mongolia only 5 sensors
%            strcmp(refpixlist{i}, '67013601') % new
%            ...strcmp(refpixlist{i}, '67010901')  % new, time series too scatter


     if   (     ...strcmp(refpixlist{i}, '03013602') || ... 
           strcmp(refpixlist{i}, '03013302') || ... 
           strcmp(refpixlist{i}, '03010903') || ...
           strcmp(refpixlist{i}, '03010908') || ...
           ...strcmp(refpixlist{i}, '04013602') || ... % 5 sensors only
           strcmp(refpixlist{i}, '04013302') || ... % 7 sensors only
           strcmp(refpixlist{i}, '04010907') || ... % new ref info 
           strcmp(refpixlist{i}, '04010910') || ... % new ref info
           ...strcmp(refpixlist{i}, '07013601') || ...
           strcmp(refpixlist{i}, '07013301') || ...
           strcmp(refpixlist{i}, '07010902') || ...
           strcmp(refpixlist{i}, '07010916') || ...
           ...strcmp(refpixlist{i}, '09013611') || ...  % replace "09013610"
           strcmp(refpixlist{i}, '09013301') || ...  % 
           strcmp(refpixlist{i}, '09010906') || ... % shifted grid (some sensor overlap with 09010902)
           ...strcmp(refpixlist{i}, '12033601') || ...   % new Tibet
           strcmp(refpixlist{i}, '12033301') || ...   % new Tibet
           ...strcmp(refpixlist{i}, '12043601') || ...    % new Twente, Holland
           ...strcmp(refpixlist{i}, '12043306') || ...    % 5 sensors
           ...strcmp(refpixlist{i}, '16013603') || ...
           strcmp(refpixlist{i}, '16013302') || ...
           strcmp(refpixlist{i}, '16010906') || ...
           strcmp(refpixlist{i}, '16010907') || ...
           strcmp(refpixlist{i}, '16010913') || ...
           ...strcmp(refpixlist{i}, '16023602') || ...
           strcmp(refpixlist{i}, '16023302') || ...
           strcmp(refpixlist{i}, '16020905') || ...  % no data
           strcmp(refpixlist{i}, '16020906') || ... % no data
           strcmp(refpixlist{i}, '16020907') || ...
           ...strcmp(refpixlist{i}, '16033602') || ...
           strcmp(refpixlist{i}, '16033302') || ...
           strcmp(refpixlist{i}, '16030911') || ...
           strcmp(refpixlist{i}, '16030916') || ...
           ...strcmp(refpixlist{i}, '16043602') || ...
           strcmp(refpixlist{i}, '16043302') || ...
           strcmp(refpixlist{i}, '16040901') || ...
           strcmp(refpixlist{i}, '16040906') || ... % new
           ...strcmp(refpixlist{i}, '16063602') || ... %***** in situ data quality
           strcmp(refpixlist{i}, '16063302') || ... %***** in situ data quality
           strcmp(refpixlist{i}, '16060907') || ...
           ...strcmp(refpixlist{i}, '16073602') || ...
           strcmp(refpixlist{i}, '16073302') || ...
           strcmp(refpixlist{i}, '16070909') || ...
           strcmp(refpixlist{i}, '16070910') || ... % one sensor data lost
           strcmp(refpixlist{i}, '16070911') || ...
           ...strcmp(refpixlist{i}, '17010907') || ... % strong Bias, not representative. new Sodankyla, Finland
           ...strcmp(refpixlist{i}, '19023601') || ... % MonteBuey
           strcmp(refpixlist{i}, '19023301') || ... %  
           strcmp(refpixlist{i}, '19020902') || ...  % subset of 19020911 
           ...strcmp(refpixlist{i}, '24013602') || ... (only 4 sensors)
           ...strcmp(refpixlist{i}, '25013601') || ...  % Tonzi Ranch
           strcmp(refpixlist{i}, '25013301') || ...
           strcmp(refpixlist{i}, '25010911') || ...
           ...strcmp(refpixlist{i}, '27013601') || ...
           strcmp(refpixlist{i}, '27013301') || ...
           strcmp(refpixlist{i}, '27010910') || ...
           strcmp(refpixlist{i}, '27010911') || ...
           strcmp(refpixlist{i}, '41010906') || ...  %new ref info 36 km in situ the same Valencia, Spain
           ...strcmp(refpixlist{i}, '45013601') || ...  % Niger
           strcmp(refpixlist{i}, '45013301') || ...  % 
           strcmp(refpixlist{i}, '45010902') || ...  % 
           ...strcmp(refpixlist{i}, '45023601') || ...  % Bernin         
           strcmp(refpixlist{i}, '45023301') || ...  %            
           strcmp(refpixlist{i}, '45020902') || ...  %            
           ...strcmp(refpixlist{i}, '48013601') || ...
           strcmp(refpixlist{i}, '48013301') || ...
           strcmp(refpixlist{i}, '48010902') || ...
           strcmp(refpixlist{i}, '48010911') || ...
           ...strcmp(refpixlist{i}, '53013601') || ...  %NEW mongolia only 5 sensors MAHASRI
           ...strcmp(refpixlist{i}, '53013301') || ...  %NEW mongolia only 5 sensors
           ...strcmp(refpixlist{i}, '67013601') || ...% new
           strcmp(refpixlist{i}, '67013301') || ...% new
           strcmp(refpixlist{i}, '67010901') ) % new,          

           is_core(i) = 1;
        
    end
end

idx_core = find(is_core);
refpixlist = refpixids(idx_core);
[cols, rows, scales, active, sites, refpixii, description, refpixids, rz, sitealias, lat,lon] = get_refpix_master_info(refpixlist);

is_core = is_core(idx_core);
is_36 = is_36(idx_core);
is_33 = is_33(idx_core);
is_09 = is_09(idx_core);
site_vec = site_vec(idx_core);

%% get prposity for the reference pixels
fname_lmc = '//discover/nobackup/projects/gmao/smap/SMAP_L4/L4_SM/Vv3030/lmc/SMAP_L4_SM_lmc_00000000T000000_Vv3030_001.h5';

poros_all = hdf5read(fname_lmc, '/Land-Model-Constants_Data/clsm_poros'); poros_all(poros_all < 0) = NaN;
wp_all = hdf5read(fname_lmc, '/Land-Model-Constants_Data/clsm_wp'); wp_all(wp_all < 0) = NaN;
clay_all = hdf5read(fname_lmc, '/Land-Model-Constants_Data/mwrtm_clay');clay_all(clay_all < 0) = NaN;
sand_all = hdf5read(fname_lmc, '/Land-Model-Constants_Data/mwrtm_sand');sand_all(sand_all < 0) = NaN;

lon_M09 = hdf5read(fname_lmc, '/cell_lon');
lon_M09 = lon_M09(:,1);

lat_M09 = hdf5read(fname_lmc, '/cell_lat');
lat_M09 = lat_M09(1,:);

%uscales =[9, 36];%unique(scales);
%uscales =[9, 33, 36];%unique(scales);
uscales =[9, 33];%unique(scales);
num_uscales = length(uscales);

num_refpix = length(refpixlist);
k = zeros(1,num_refpix);

for n = 1:num_uscales
    
    if n==1
        ind_s = find(is_09);
    elseif n ==2
        ind_s = find(is_33);
    else
        ind_s = find(is_36);
    end
    
    for s = 1:length(ind_s)
        
        k(ind_s(s)) = k(ind_s(s)) + 1;
        
        % note that this depends on the scale!!!!
        if uscales(n) == 9 ||  uscales(n) == 33
            
            if (mod(cols(ind_s(s)),1) == 0 && mod(rows(ind_s(s)),1) == 0)  % regular pixels
                
                poros{ind_s(s)}(k(ind_s(s)), 1) = poros_all(cols(ind_s(s)), rows(ind_s(s)));
                wp{ind_s(s)}(k(ind_s(s)), 1)    = wp_all(cols(ind_s(s)), rows(ind_s(s)));
                clay{ind_s(s)}(k(ind_s(s)), 1)  = clay_all(cols(ind_s(s)), rows(ind_s(s)));
                sand{ind_s(s)}(k(ind_s(s)), 1)  = sand_all(cols(ind_s(s)), rows(ind_s(s)));
                
            else % displaced pixels
                
                cols_1 = floor(cols(ind_s(s)));
                cols_2 = cols_1 + 1;
                
                rows_1 = floor(rows(ind_s(s)));
                rows_2 = rows_1 + 1;
                
                poros_tmp(1) = poros_all(cols_1, rows_1);
                poros_tmp(2) = poros_all(cols_2, rows_1);
                poros_tmp(3) = poros_all(cols_1, rows_2);
                poros_tmp(4) = poros_all(cols_2, rows_2);
                
                wp_tmp(1) = wp_all(cols_1, rows_1);
                wp_tmp(2) = wp_all(cols_2, rows_1);
                wp_tmp(3) = wp_all(cols_1, rows_2);
                wp_tmp(4) = wp_all(cols_2, rows_2);
                
                clay_tmp(1) = clay_all(cols_1, rows_1);
                clay_tmp(2) = clay_all(cols_2, rows_1);
                clay_tmp(3) = clay_all(cols_1, rows_2);
                clay_tmp(4) = clay_all(cols_2, rows_2);
                
                sand_tmp(1) = sand_all(cols_1, rows_1);
                sand_tmp(2) = sand_all(cols_2, rows_1);
                sand_tmp(3) = sand_all(cols_1, rows_2);
                sand_tmp(4) = sand_all(cols_2, rows_2);
                
                weight(1) = (1-abs(cols_1 - cols(ind_s(s))))*(1-abs(rows_1 - rows(ind_s(s))));
                weight(2) = (1-abs(cols_2 - cols(ind_s(s))))*(1-abs(rows_1 - rows(ind_s(s))));
                weight(3) = (1-abs(cols_1 - cols(ind_s(s))))*(1-abs(rows_2 - rows(ind_s(s))));
                weight(4) = (1-abs(cols_2 - cols(ind_s(s))))*(1-abs(rows_2 - rows(ind_s(s))));
                
                if all(isnan(poros_tmp))
                    poros{ind_s(s)}(k(ind_s(s)), 1) = NaN;
                else
                    poros{ind_s(s)}(k(ind_s(s)), 1) = wmean(poros_tmp, weight);
                end
                
                if all(isnan(wp_tmp))
                    wp{ind_s(s)}(k(ind_s(s)), 1) = NaN;
                else
                    wp{ind_s(s)}(k(ind_s(s)), 1) = wmean(wp_tmp, weight);
                end
                
                if all(isnan(clay_tmp))
                    clay{ind_s(s)}(k(ind_s(s)), 1) = NaN;
                else
                    clay{ind_s(s)}(k(ind_s(s)), 1) = wmean(clay_tmp, weight);
                end
                
                if all(isnan(sand_tmp))
                    sand{ind_s(s)}(k(ind_s(s)), 1) = NaN;
                else
                    sand{ind_s(s)}(k(ind_s(s)), 1) = wmean(sand_tmp, weight);
                end
                
            end
            
        elseif uscales(n) == 36 
            
            [col_3_tmp, row_3_tmp, col_9, row_9, col_36_tmp, row_36_tmp] = ...
                EASEv2_3k_9k_36k_index_convert(cols(ind_s(s)), rows(ind_s(s)), 36,1);
            
            clear col_3_tmp row_3_tmp col_36_tmp row_36_tmp
            
            if ( mod(single(col_9),1) ~= 0.5 || mod(single(row_9),1) ~= 0.5)
                
               col_9 = single(round(col_9/0.5)*0.5);
               row_9 = single(round(row_9/0.5)*0.5);
                
            end
            
            poros_9 = poros_all(int32(col_9-1.5):int32(col_9+1.5), int32(row_9-1.5):int32(row_9+1.5));
            wp_9    = wp_all(int32(col_9-1.5):int32(col_9+1.5), int32(row_9-1.5):int32(row_9+1.5));
            clay_9  = clay_all(int32(col_9-1.5):int32(col_9+1.5), int32(row_9-1.5):int32(row_9+1.5));
            sand_9  = sand_all(int32(col_9-1.5):int32(col_9+1.5), int32(row_9-1.5):int32(row_9+1.5));
            
            poros{ind_s(s)}(k(ind_s(s)), 1) = nanmean(poros_9(:));
            wp{ind_s(s)}(k(ind_s(s)), 1) = nanmean(wp_9(:));
            clay{ind_s(s)}(k(ind_s(s)), 1) = nanmean(clay_9(:));
            sand{ind_s(s)}(k(ind_s(s)), 1) = nanmean(sand_9(:));
            
            clear poros_9 wp_9 clay_9 sand_9 col_9 row_9
            
      end
   end
end

nn = length(refpixlist);
nv = length(var_tag);
nc = length(CalVal_path);

R   = NaN*ones(nn,nv,nc);
RLO = NaN*ones(nn,nv,nc) ;
RUP = NaN*ones(nn,nv,nc);

if add_anomR
    anomR   = NaN*ones(nn,nv,nc);
    anomRLO = NaN*ones(nn,nv,nc);
    anomRUP = NaN*ones(nn,nv,nc);
end

Bias   = NaN*ones(nn,nv,nc);
BiasLO = NaN*ones(nn,nv,nc);
BiasUP = NaN*ones(nn,nv,nc);

absBias = NaN*ones(nn,nv,nc);
absBiasLO = NaN*ones(nn,nv,nc);
absBiasUP = NaN*ones(nn,nv,nc);

RMSE =   NaN*ones(nn,nv,nc);
ubRMSE =  NaN*ones(nn,nv,nc);
ubRMSELO = NaN*ones(nn,nv,nc);
ubRMSEUP = NaN*ones(nn,nv,nc);

Nsens_max = NaN*ones(nn,nv,nc);
Nsens_min = NaN*ones(nn,nv,nc) ;
Nsens_mean = NaN*ones(nn,nv,nc);
Nsens_mode = NaN*ones(nn,nv,nc);

Nobs_valid = NaN*ones(nn,nv,nc);

clear tmpdata1
% loop through refpix
for kk = find(is_33) % 1:length(refpixlist)   %127 
    
    clear  date_vec1 DATA_L4_std1 DATA_RP_std1 RP_prcp1
    is_root(kk) = 0;
    
    for iv = 1:length(var_tag)
        
        clear tmpdata matchupdir date_vec DATA_RP DATA_L4 DATA_L4_std L4_std RP_std RP_prcp refdata refdata_std DATA data
        
        for cc = 1:length(CalVal_path)
            tmp = dir([CalVal_path{cc}, refpixlist{kk}(1:4), '_*']);
            if isempty(tmp)
                disp(['no data for ', refpixlist{kk}(1:4)])
                matchupdir = [];
                break
            end
            site_name = tmp.name(6:end); clear tmp
            
            matchupdir = [CalVal_path{cc}, refpixlist{kk}(1:4), '_',site_name,'/SMAPL4',var_tag{iv}, '/', ...
                refpixlist{kk}, '/'];
            
            if exist(matchupdir, 'dir')
                
                fullfilenames = findfilenames(matchupdir, startdate, enddate);
                
                if ~isempty(fullfilenames)
                    
                    DATA = [];
                    Yr = [];
                    Mo =[];
                    Day = [];
                    TOY = [];
                    
                    for i = 1:length(fullfilenames)
                        
                        [hdr1_1, cols_1, dtime_1, ALL_1, ID_1, Yr_1, Mo_1, Day_1, Hr_1, Min_1, TOY_1, DATA_1]...
                            = read_data_file(fullfilenames{i});
                        
                        DATA = [DATA; DATA_1];
                        Yr = [Yr; Yr_1];
                        Mo =[Mo; Mo_1];
                        Day = [Day; Day_1];
                        TOY = [TOY; TOY_1];
                        
                    end
                else
                    disp(['no ', var_tag{iv}, ' matchup file found for the period for ', refpixlist{kk}])
                    break
                end
                
            else
                
                disp(['no ', var_tag{iv}, ' matchup file found for the period for ', refpixlist{kk}])
                break
            end
            
            %if ~isempty(strfind(CalVal_path{cc}, 'SMAP_Nature'))
            %   TOY = TOY + 0.0625;
            %end
            
            t2k = resolve_days_since_2000(Yr, TOY);
            
            date_vec{cc} = t2k - resolve_days_since_2000(str2num(startdate(1:4)), 1);
            
            prec = DATA(:,9);
            Nsens = DATA(:,10);
            st = DATA(:,3);
            st_std = DATA(:,7);
            
            Nsens(Nsens < 1) = NaN;
            Nsens_max(kk,iv,cc) = max(Nsens);
            Nsens_min(kk,iv,cc) = min(Nsens);
            Nsens_mean(kk,iv,cc) = nanmean(Nsens);
            Nsens_mode(kk,iv,cc) = mode(Nsens);
            refdata = DATA(:,1); 
            
            refdata_std = DATA(:,5); 
            refdata_std(isnan(refdata)) = NaN;
            
            figure, plot([refdata refdata_std]),title(refpixlist{kk})
            [nanmean(refdata) nanmean(refdata_std)]
            
            if QCflagyes
                QCflag = DATA(:,11);
                cl = 12;
            else
                cl = 11;
            end
            
            if ~do_ST
                data = DATA(:,cl);
                if strcmp(refpixlist{kk}(1:5),'16043') && iv ==1 
                    refdata = 0.55*(0.65*refdata+0.0462)+0.25*(1.5263*refdata+0.1837)+0.20*(1.04*refdata-0.0249);
                end
                DATA_RP = refdata;
                RP_std = refdata_std;
                DATA_L4_std = DATA(:,cl+1);
                DATA_L4 = data;
            else
                data = DATA(:,cl+2);
                DATA_RP = st;
                RP_std = st_std;
                if cc <= 3
                    DATA_L4_std = st_std;
                end
                DATA_L4 = data;
            end
            
            %if QCflagyes
            %   DATA_RP(find(QCflag == 1)) = NaN;
            %end
            
            %DATA_L4(isnan(DATA_RP)) = NaN;
            %DATA_RP(isnan(DATA_L4)) = NaN;
            %DATA_L4_std(isnan(DATA_RP)) = NaN;
            
            tmpdata{cc} = [ DATA_RP DATA_L4];

            %L4_std{cc} = DATA_L4_std;
            
           RP_prcp = prec;
            
        end
 
        % ..........................................................................................

        % cross masking
        if exist('tmpdata', 'var') 
            [datevec, tind1, tind2] = intersect(round(date_vec{1}*8), round(date_vec{2}*8));
            tmpdata{1} = tmpdata{1}(tind1,:);
            tmpdata{2} = tmpdata{2}(tind2,:);

%              [datevec_tmp, tind1, tind2] = intersect(round(date_vec{1}*8), round(date_vec{2}*8));
%              tmpdata{1} = tmpdata{1}(tind1,:);
%              tmpdata{2} = tmpdata{2}(tind2,:);
%              
%             clear tind1 tind2
%               

% 
%          [datevec_tmp, tind1, tind2] = intersect(round(datevec_tmp), round(date_vec{3}*8));
%          tmpdata{1} = tmpdata{1}(tind1,:);
%          tmpdata{2} = tmpdata{2}(tind1,:);
%          tmpdata{3} = tmpdata{3}(tind2,:);
%             
%              clear tind1 tind2
%   
%              [datevec, tind1, tind2] = intersect(round(datevec_tmp), round(date_vec{4}*8));
%              tmpdata{1} = tmpdata{1}(tind1,:);
%              tmpdata{2} = tmpdata{2}(tind1,:);
%              tmpdata{3} = tmpdata{3}(tind1,:);
%              tmpdata{4} = tmpdata{4}(tind2,:);
% 
%            clear tind1 tind2
            
%            [datevec, tind1, tind2] = intersect(round(datevec_tmp), round(date_vec{5}*8));
%            tmpdata{1} = tmpdata{1}(tind1,:);
%            tmpdata{2} = tmpdata{2}(tind1,:);
%            tmpdata{3} = tmpdata{3}(tind1,:);
%            tmpdata{4} = tmpdata{4}(tind1,:);
%            tmpdata{5} = tmpdata{5}(tind2,:);
                     
             DATA_L4_std = DATA_L4_std(tind2);
             RP_prcp = RP_prcp(tind2);
             RP_std = RP_std(tind2);
            
            date_vec = datevec/8.;
            
            clear datevec tind1 tind2
            
           % cross mask in situ data as well
           tmpdata{1}(isnan(tmpdata{2}(:,1)),1) = NaN;
           tmpdata{2}(isnan(tmpdata{1}(:,1)),1) = NaN;
            
           % cross mask L4 data 
           tmpdata{1}(isnan(tmpdata{2}(:,2)),2) = NaN;
           tmpdata{2}(isnan(tmpdata{1}(:,2)),2) = NaN;
            
             % cross mask in situ data as well
%              tmpdata{4}(isnan(tmpdata{1}(:,1)),1) = NaN;
%              tmpdata{4}(isnan(tmpdata{2}(:,1)),1) = NaN;
%              tmpdata{4}(isnan(tmpdata{3}(:,1)),1) = NaN;
% %             tmpdata{4}(isnan(tmpdata{5}(:,1)),1) = NaN;
%              tmpdata{1}(isnan(tmpdata{4}(:,1)),1) = NaN;
%              tmpdata{2}(isnan(tmpdata{4}(:,1)),1) = NaN;
%              tmpdata{3}(isnan(tmpdata{4}(:,1)),1) = NaN;
% %             tmpdata{5}(isnan(tmpdata{4}(:,1)),1) = NaN;
%             
%              % cross mask L4 data 
%              tmpdata{4}(isnan(tmpdata{1}(:,1)),1) = NaN;
%              tmpdata{4}(isnan(tmpdata{2}(:,1)),1) = NaN;
%              tmpdata{4}(isnan(tmpdata{3}(:,1)),1) = NaN;
% %             tmpdata{4}(isnan(tmpdata{5}(:,1)),1) = NaN;
%              tmpdata{1}(isnan(tmpdata{4}(:,2)),2) = NaN;
%              tmpdata{2}(isnan(tmpdata{4}(:,2)),2) = NaN;
%              tmpdata{3}(isnan(tmpdata{4}(:,2)),2) = NaN;
% %             tmpdata{5}(isnan(tmpdata{4}(:,2)),2) = NaN;
        end

      
        for cc = 1:nc
            if exist('tmpdata','var')
                % Compute statistical metrics
                
                clear stats stats_tags
                [stats, stats_tags] = get_validation_stats(...
                    tmpdata{cc}, 1, 'complete', 1, [1:2], Nmin, 1 );
                
                nodata_val = -9999;
                nodata_tol =  1E-4;
                
                stats(abs(stats-nodata_val)<nodata_tol ) = NaN;
                stats(stats == 0 ) = NaN;
                
                
                R(kk,iv,cc)   = stats(strcmp(stats_tags,'R'));
                RLO(kk,iv,cc) = stats(strcmp(stats_tags,'RLO')) - R(kk,iv,cc) ;
                
                RUP(kk,iv,cc) = stats(strcmp(stats_tags,'RUP')) - R(kk,iv,cc);
                
                Bias(kk,iv,cc)   = -stats(strcmp(stats_tags,'bias'));
                BiasLO(kk,iv,cc) = -stats(strcmp(stats_tags,'CI_bias')) ;
                BiasUP(kk,iv,cc) =  stats(strcmp(stats_tags,'CI_bias')) ;
                
                absBias(kk,iv,cc) =  abs(stats(strcmp(stats_tags,'bias')));
                absBiasLO(kk,iv,cc) = -stats(strcmp(stats_tags,'CI_bias')) ;
                absBiasUP(kk,iv,cc) =  stats(strcmp(stats_tags,'CI_bias')) ;
                
                RMSE(kk,iv,cc) =     sqrt( stats(strcmp(stats_tags,'MSE')) );
                ubRMSE(kk,iv,cc) =   sqrt(stats(strcmp(stats_tags,'MSE')) - stats(strcmp(stats_tags,'bias')).^2);
                ubRMSELO(kk,iv,cc) = sqrt(stats(strcmp(stats_tags,'ubMSELO'))) -ubRMSE(kk,iv,cc);
                ubRMSEUP(kk,iv,cc) = sqrt(stats(strcmp(stats_tags,'ubMSEUP'))) -ubRMSE(kk,iv,cc); 
                
                Nobs_valid(kk,iv,cc) = length(find(~isnan(tmpdata{cc}(:,1))));
                %tmpdata{cc}(find(date_vec < 91),:) = NaN; % remove data before mission start time 03/31/2015
                tmpdata1{kk,cc,iv} = tmpdata{cc};
                date_vec1{iv} = date_vec;
                DATA_L4_std1{cc,iv} = DATA_L4_std;
                DATA_RP_std1{iv} = RP_std;
                RP_prcp1{iv} = RP_prcp;
                
                if add_anomR
                    
                    insitu_data = tmpdata{cc}(:,1);
                    model_data  = tmpdata{cc}(:,2);
                    
                    % cross mask data
                    insitu_data(isnan(model_data)) = NaN;
                    model_data(isnan(insitu_data)) = NaN;
                    
                    month_vec = NaN *ones(size(model_data));
                    year_vec = NaN * ones(size(model_data));
                    
                    disp('computing anomalies...')
                    Nmin_mon = 80;
                    Nmin_year = 3;
                    
                    % compute monthly mean
                    day_in_month = [31 28 31 30 31 30 31 31 30 31 30 31];
                    day_in_month_leap = [31 29 31 30 31 30 31 31 30 31 30 31];
                    
                    for year = 1: ceil(date_vec(end)/365)
                        
                        this_year = start_year + year - 1;
                        
                        for month = 1:12
                            
                            if this_year < 2016
                                
                                if month > 1
                                    ind = intersect(find((date_vec-(year-1)*365) > sum(day_in_month(1:month-1))), ...
                                        find((date_vec-(year-1)*365) <= sum(day_in_month(1:month))));
                                else
                                    
                                    ind = intersect(find((date_vec-(year-1)*365) <=sum(day_in_month(1))), ...
                                       find((date_vec-(year-1)*365) > 0)) ;
                                end
                                
                            elseif this_year == 2016
                                
                                if month > 1
                                    ind = intersect(find((date_vec-(year-1)*365) > sum(day_in_month_leap(1:month-1))), ...
                                        find((date_vec-(year-1)*365) <= sum(day_in_month_leap(1:month))));
                                else
                                    
                                    ind = intersect(find((date_vec-(year-1)*365) <=sum(day_in_month_leap(1))), ...
                                        find((date_vec-(year-1)*365) > 0));
                                end
                            else
                                if month > 1
                                    ind = intersect(find((date_vec-(year-1)*365-1) > sum(day_in_month(1:month-1))), ...
                                        find((date_vec-(year-1)*365-1) <= sum(day_in_month(1:month))));
                                else
                                    
                                    ind = intersect(find((date_vec-(year-1)*365-1) <=sum(day_in_month(1))), ...
                                        find((date_vec-(year-1)*365-1) > 0));
                                end
                            end
                            
                            if ~isempty(ind)
                                month_vec(ind) = month;
                                year_vec(ind) = year;
                            end
                            
                            tmp = insitu_data(ind);
                            tmp = tmp(~isnan(tmp));
                            if length(tmp) >= Nmin_mon
                                tmp = mean(tmp);
                            else
                                tmp = NaN;
                            end
                            insitu_monthlymean(year,month) = tmp; clear tmp
                            
                            tmp = model_data(ind);
                            tmp = tmp(~isnan(tmp));
                            if length(tmp) >= Nmin_mon
                                tmp = mean(tmp);
                            else
                                tmp = NaN;
                            end
                            model_monthlymean(year,month) = tmp; clear tmp
                            
                        end
                    end
                    
                    for month = 1:12
                        tmp = insitu_monthlymean(:,month);
                        tmp = tmp(~isnan(tmp));
                        if length(tmp) >= Nmin_year
                            tmp = mean(tmp);
                        else
                            tmp = NaN;
                        end
                        insitu_clim(month) = tmp;
                        
                        tmp = model_monthlymean(:,month);
                        tmp = tmp(~isnan(tmp));
                        if length(tmp) >= Nmin_year
                            tmp = mean(tmp);
                        else
                            tmp = NaN;
                        end
                        model_clim(month) = tmp; clear tmp
                        
                    end
                    
                    insitu_anom = NaN *ones(size(insitu_data));
                    model_anom = NaN * ones(size(model_data));
                    
                    for month = 1:12
                        ind = find(month_vec == month);
                        if ~isempty(ind)
                            insitu_anom(ind) = insitu_data(ind) - insitu_clim(month);
                            model_anom(ind) = model_data(ind) - model_clim(month);
                        end
                    end
                    
                    clear stats stats_tags
                    [stats, stats_tags] = get_validation_stats(...
                        [insitu_anom model_anom], 1, 'complete', 1, [1:2], Nmin, 1 );
                    
                    nodata_val = -9999;
                    nodata_tol =  1E-4;
                    
                    stats(abs(stats-nodata_val)<nodata_tol ) = NaN;
                    stats(stats == 0 ) = NaN;
                    
                    
                    anomR(kk,iv,cc)   = stats(strcmp(stats_tags,'R'));
                    anomRLO(kk,iv,cc) = stats(strcmp(stats_tags,'RLO')) - anomR(kk,iv,cc) ;
                    anomRUP(kk,iv,cc) = stats(strcmp(stats_tags,'RUP')) - anomR(kk,iv,cc);
                    
                    tmp = corrcoef([insitu_monthlymean(:) model_monthlymean(:)],'rows','complete');
                    R_monthlymean(kk,iv,cc) = tmp(2,1); clear tmp
                end
        
                clear stats stats_tags
            else
                tmpdata1{kk,cc,iv} = [NaN NaN];
            end
        end
        
        
    end
        
end

clear sfmc
for kk = 1:size(tmpdata1)
sfmc{kk}(:,1) = tmpdata1{kk,1,1}(:,2);
sfmc{kk}(:,2) = tmpdata1{kk,2,1}(:,2);
% sfmc{kk}(:,3) = tmpdata1{kk,3,1}(:,2);
% sfmc{kk}(:,4) = tmpdata1{kk,4,1}(:,2);
% sfmc{kk}(:,5) = tmpdata1{kk,4,1}(:,1);
% sfmc{kk}(isnan(sfmc{kk}(:,1)),5) = NaN;
% sfmc{kk}(isnan(sfmc{kk}(:,5)),1:4) = NaN;
end


% =================== EOF ======================================
