function [] = get_cdf_match( varname,                                   ...
    run_months, exp_path, exp_run, domain, start_year, end_year, ...
    species, hscale,  inv_cdf_match, time_of_day_in_hours )

%
% get_cdf_match.m
%
% Match cdfs from stats files (produced with get_model_and_obs_stats.m),
% compute polynomial fit, write files that contain cdf-matching parameters.
%
% reichle, 22 Aug 2005
%
% GDL, 16 Nov 2010
% GDL, 03 May 2010 - added option to find an polynomial fit for
%                    `inverse' CDF matching, i.e. to transform model
%                    output to an observed climatology (inv_cdf_match)
% CSD, 31 Jan 2014. 
%
% -------------------------------------------------------------------
% begin user-defined inputs
% -------------------------------------------------------------------

% minimum number of data points to include in statistics

N_data_min = 100;        %Was 100 for global!!!

% no-data-value for points that don't have good statistics

no_data_stats = -9999.;

% option = 1 - write lookup table
% option = 2 - write cubic interpolation of lookup table

%option = [ 1 2 ];              % poly interpolation and lookup table
option = 2;                      % poly interpolation only

% order of polynomial for interpolation of lookup table

N_poly = 3;   % 

% -------------------------------------------------------------------
% end user-defined inputs
% -------------------------------------------------------------------

if ~exist('hscale','var')
  
  hscale = 0;
  
end

if ~exist('inv_cdf_match','var')
  
  inv_cdf_match = 0;
  
end

% assemble input and output path

filepath  = [ exp_path, '/', exp_run, '/', domain ];

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

% assemble input and output file names

month_string = {'Jan'; 'Feb'; 'Mar'; 'Apr'; 'May'; 'Jun'; ...
		'Jul'; 'Aug'; 'Sep'; 'Oct'; 'Nov'; 'Dec'};

if(strcmp(varname,'sfmc'))
    spec_tag = 'SM';
elseif(strcmp(varname,'''Tb'''))
    spec_tag = 'Tb';
elseif(strcmp(varname,'Tskin'))
    spec_tag = 'Tskin';
else
   disp('unknown varname') 
end

for s=1:length(species)
    spec_tag = [spec_tag,'_',num2str(species(s))];
end

% tmpstr = [ 'species' ];
% 
% for i=1:length(species)
%   
%   tmpstr = [ tmpstr, '_', num2str(species(i)) ];
%   
% end

fname_in  = [ filepath, '/stats/', exp_run, '.stats',... %.', ...%tmpstr,           ...
	      '.hscale_', num2str(hscale,'%2.2f'), '_',                  ...
	      num2str(start_year), '-', num2str(end_year), '.' ];

fname_out_lookup = [ filepath, '/stats/', exp_run, '.cdf_match.lookup.', ...
		    '.hscale_', num2str(hscale,'%2.2f'), '_',            ...
		    num2str(start_year), '-', num2str(end_year), '.' ];

% shortened "cdf_match" to "cdf" for RedArkOSSE - reichle, 26 Sep 2006

fname_out_poly = [ filepath, '/stats/', exp_run, '.cdf.poly',     ...
		   '.hscale_', num2str(hscale,'%2.2f'), '_',             ...
		   num2str(start_year), '-', num2str(end_year), '.' ];

if length(run_months)==12
  fname_in         = [fname_in,         'all_months' ];
  fname_out_lookup = [fname_out_lookup, 'all_months' ];
  fname_out_poly   = [fname_out_poly,   'all_months' ];
else
  for i=1:length(run_months)
    fname_in         = [fname_in,  month_string{run_months(i)} ];
    fname_out_lookup = [fname_out_lookup, month_string{run_months(i)} ];
    fname_out_poly   = [fname_out_poly,   month_string{run_months(i)} ];
  end
end

fname_in         = [fname_in,         spec_tag];
fname_out_lookup = [fname_out_lookup, spec_tag];
fname_out_poly   = [fname_out_poly,   spec_tag];


if exist( 'time_of_day_in_hours', 'var')
  
  fname_out_lookup = [fname_out_lookup, '_', num2str(time_of_day_in_hours,'%2.2d'), 'z'];
  fname_out_poly   = [fname_out_poly, '_', num2str(time_of_day_in_hours,'%2.2d'), 'z'];
  fname_in   = [fname_in, '_', num2str(time_of_day_in_hours,'%2.2d'), 'z'];

end

fname_in         = [fname_in,         '.dat'];

if (inv_cdf_match ~= 1)
fname_out_lookup = [fname_out_lookup, '.dat'];
fname_out_poly   = [fname_out_poly,   '.dat'];
else
fname_out_lookup = [fname_out_lookup, '.inv.dat'];
fname_out_poly   = [fname_out_poly,   '.inv.dat'];    
end

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

% load catchment coordinates

fname = [ exp_path '/' exp_run '/output/' domain '/rc_out', '/', exp_run, '.ldas_tilecoord.bin'];

[ tile_coord ] = read_tilecoord( fname );

N_tile = length(tile_coord.tile_id);
  
% -------------------------------------------------------------

% load stats file

disp(['loading stats file ', fname_in])

stats_data = load( fname_in );

disp('done loading stats file')

if ( any(stats_data(:,8)~=stats_data(1,8)) |      ...
     any(stats_data(:,9)~=stats_data(1,9))    )
  input('something wrong (1), ctrl-c now')
end

edge_min = stats_data(1,8);
edge_max = stats_data(1,9);

edge_dx  = (edge_max - edge_min) / (size(stats_data,2) - 9);

edges = edge_min:edge_dx:edge_max;


% initialize to no-data-value

obs_prime = no_data_stats*ones(N_tile,length(edges));

std_obs   = no_data_stats*ones(N_tile,1);
std_mod   = no_data_stats*ones(N_tile,1);

min_obs   = no_data_stats*ones(N_tile,1);
max_obs   = no_data_stats*ones(N_tile,1);

fit_coeff = no_data_stats*ones(N_tile,N_poly+1);

disp('processing statistics...')

for i=1:N_tile
  
  %if (i==2539)
    
  obs_ind = 2*i-1;
  mod_ind = 2*i;
  
  if (inv_cdf_match == 1) 
     mod_ind = 2*i-1;
     obs_ind = 2*i; 
  end
  
  % check tile id
  
  if (stats_data(obs_ind,1)~=stats_data(mod_ind,1))
    input('something wrong (2), ctrl-c now')
  end
  
  N_obs = sum(stats_data(obs_ind,10:end));
  N_mod = sum(stats_data(mod_ind,10:end));
  
  if (N_obs~=N_mod)
    disp(num2str([tile_coord.tile_id(i), N_obs, N_mod]))
    input('something wrong (3), ctrl-c now')
  end
  
  if (N_obs>N_data_min)
    
    obs_cdf = cumsum(stats_data(obs_ind,10:end))/N_obs;
    mod_cdf = cumsum(stats_data(mod_ind,10:end))/N_mod;  
    
    std_obs(i,1) = sqrt(stats_data(obs_ind,4));  
    std_mod(i,1) = sqrt(stats_data(mod_ind,4));  
    
    min_obs(i,1) = stats_data(obs_ind,6);  
    max_obs(i,1) = stats_data(obs_ind,7);  
    
    % fill in lookup table
    
    obs_prime(i,1) = edge_min;
   
    for j=1:length(obs_cdf) % each j is a bin
      
      tmpdist = abs( obs_cdf(j) - mod_cdf );
      
      ind = find( tmpdist==min(tmpdist) );
      
      obs_prime(i,j+1) = edges(ind(1));
      
    end
    
    % polynomial interpolation of lookup table
    
    % min and max obs determine acceptable range of polynomial fit
    
    ind = find( edges>=min_obs(i,1)  &  edges<=max_obs(i,1) );

    if (length(ind) >= 2*(N_poly+1) )
      
      tmp_fit_coeff = polyfit( edges(ind), obs_prime(i,ind), N_poly );

      fit_coeff( i, 1:N_poly+1 ) = tmp_fit_coeff;
      
    else
      
      disp(['get_cdf_match(): obs range insufficient for polyfit in tile ', ...
	    num2str( tile_coord.tile_id(i) ) ])
      disp('CHOOSE FINER EDGES IN get_model_and_obs_stats.m')
      
    end
    
  end
  
 %   end
  
end

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

% write output - lookup table

if any(option==1)
  
  disp(' ')
  disp(['writing ', fname_out_lookup])
  
  ofp = fopen(fname_out_lookup, 'w');
  
  format_string = ['%10d%13.5e%13.5e'];
  for i=1:length(obs_prime)
    format_string = [ format_string, '%11.3e'];
  end
  format_string = [ format_string, '\n'];
  
  % header line
  
  fprintf(ofp, format_string, no_data_stats, no_data_stats*[1 1], edges );
  
  % scaled obs information for each tile
  
  for k=1:N_tile
    
    fprintf(ofp, format_string,                             ...
	    tile_coord.tile_id(k), std_obs(k), std_mod(k),  ...
	    obs_prime(k,:) );
    
  end
  
  fclose(ofp);
  
end
  
% ----------------------------------------------

% write output - polynomial fit of lookup table

if any(option==2)
  
  disp(' ')
  disp(['writing ', fname_out_poly])
  
  ofp = fopen(fname_out_poly, 'w');
  
  format_string = ['%10d%13.5e%13.5e%13.5e%13.5e'];
  for i=1:size(fit_coeff,2)
    format_string = [ format_string, '%13.5e'];
  end
  format_string = [ format_string, '\n'];
  
  % header line 
  
  fprintf(ofp, '%10d%10d%13.5e%13.5e%13.5e', ...
	  N_tile, N_poly, edge_min, edge_max, edge_dx );
  
  % add a few no_data values so resulting files contains a complete
  % matrix that can be read easily with matlab
  
  for i=1:size(fit_coeff,2)
    fprintf(ofp, '%13.5e', no_data_stats);
  end
  fprintf(ofp, '\n');
  
  
  % scaled obs information for each tile
  
  for k=1:N_tile
    
    fprintf(ofp, format_string,                             ...
	    tile_coord.tile_id(k), std_obs(k), std_mod(k),  ...
	    min_obs(k), max_obs(k), fit_coeff(k,:) );
    
  end
  
  fclose(ofp);
  
end
  
% ============== EOF ===============================================
