% Calculate scaling files for SMAP L4 Tb-DA based on SMOS Tb
%
% GDL,  11 Sep 2012
% QLiu, 12 2016
% sqz,  2 2020
% =====================================================================

clear

addpath /home/qliu/projects/matlab_code/matlab_base

%run_months = [1:12 1 2]; %loop around months get complete pentads
run_months = [1:12 1 2];

exp_path = '/discover/nobackup/qzhang/';
exp_run = {'C90omfe24'};
domain   = 'CF0090x6C_GLOBAL';

prefix_out = 'SMAPc90e24_zscore_stats_';

%%%  C90omfe24 
start_year = [repmat(2016,1,3) repmat(2015,1,9) repmat(2016,1,2)]; %corresp to [1:12 1 2]
end_year   = [repmat(2020,1,3) repmat(2019,1,9) repmat(2020,1,2)];


% Run Asc and Desc separately
%run one and two in two runs 
orbit    = [ 1 ]; %1=A, 2=D   !DO *NOT* USE ASC AND DESC TOGETHER!

% Don't change
pol      = [ 1 2 ]; %1=H, 2=V
inc_ang  = [ 40.0 ];


dt_assim   = 3*60*60;

%======

obs_param_fname = [exp_path, '/', exp_run{1}, '/output/', domain, '/rc_out/', ...
     '/Y2015/M05/',exp_run{1}, '.ldas_obsparam.20150501_2100z.txt'];

var_name = {'Tb'};

descr = 'SMAP_L1C' ; % for identify obs species

%======
if (length(orbit) > 1)
    error('ONLY pick one orbit!')
end

if (orbit(1) == 1) int_Asc = 1; end %Asc
if (orbit(1) == 2) int_Asc = 0; end %Desc


%Spatial sampling
hscale = 0.0;          % degrees lat/lon

%Temporal sampling (days), current hard coded and need to be divisive by 5 and be an odd number
w_days    = 45; %95;

Ndata_min = 20 ;

if (mod(w_days,10) == 0)
    disp('w_days should be 5, 15, 25, 35, ...')
    error('Need an odd number of pentads |xxxxx|xxXxx|xxxxx|')
end
if (mod(w_days, 5) > 0)
    error('Aiming at pentad files')
end

% ------------------------------------------------------------------------

[N_obs_param, obs_param ] = read_obsparam_2018(obs_param_fname);

species =[];

for oo=1:length(orbit)
    for pp=1:length(pol)
        for aa=1:length(inc_ang)
            
            add_species = obs_param(strcmp(var_name,{obs_param.varname}) & ...
                orbit(oo) == [obs_param.orbit] & ...
                inc_ang(aa) == [obs_param.ang] & ...
                pol(pp) == [obs_param.pol] & ...
                ~cellfun(@isempty, strfind({obs_param.descr},descr))).species;
            
            species = union(species,add_species);
            
        end
    end
end

% ------------------

for n=1:length(exp_run)
    
    if (exist('convert_grid','var'))
        
        if exist('time_of_day_in_hours','var')
            
            
            for j=1:length(time_of_day_in_hours)
                
                for k=1:length(run_months)
                    
                    get_model_and_obs_clim_stats( var_name,               ...
                        run_months{k}, exp_path, exp_run{n}, domain,     ...
                        start_year, end_year, ...
                        dt_assim, species, obs_param, ...
                        hscale, inc_ang, int_Asc, w_days, Ndata_min, prefix_out,...
                        convert_grid, time_of_day_in_hours(j)  );
                    
                end
                
            end
            
        else
            
            get_model_and_obs_clim_stats_GEOSldas( var_name,                              ...
                run_months, exp_path, exp_run{n}, domain, start_year, end_year, ...
                dt_assim, species, obs_param, ...
                hscale, inc_ang, int_Asc, w_days, Ndata_min, prefix_out,...
                convert_grid );
            
        end
    else
        
        if exist('time_of_day_in_hours','var')
            
            
            for j=1:length(time_of_day_in_hours)
                
                for k=1:length(run_months)
                    
                    get_model_and_obs_clim_stats( var_name,               ...
                        run_months{k}, exp_path, exp_run{n}, domain,     ...
                        start_year, end_year, ...
                        dt_assim, species, obs_param, ...
                        hscale, inc_ang, int_Asc, w_days, Ndata_min, prefix_out,...
                        time_of_day_in_hours(j)  );
                    
                end
                
            end
            
        else
 
%SMAP L1 here            
            get_model_and_obs_clim_stats( var_name,                              ...
                run_months, exp_path, exp_run{n}, domain, start_year, end_year, ...
                dt_assim, species, obs_param, ...
                hscale, inc_ang, int_Asc, w_days, Ndata_min, prefix_out);
            
        end
        
    end
end


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


