clear all
close all

%% Pull footprint grids into matlab and use emissions inventory to calculate concentrations.

days = [19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30];
year = [2015,2015,2015,2015,2015,2015,2015,2015,2015,2015,2015,2015];
mont = [4,4,4,4,4,4,4,4,4,4,4,4];

footprint = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/stilt_output/roosevelt_april2015/footprint_05/2015041900_footprint.nc','footprint');

time_model = datenum(year(1),mont(1),days(1),0,0,0);

for d = 1:length(days)
 
% Deal with time to read in the correct file.

   if days(d) < 10
      day = strcat('0',int2str(days(d)));
   elseif days(d) >= 10
      day = int2str(days(d));
   end
   
   if mont(d) < 10
      mon = strcat('0',int2str(mont(d)));
   elseif mont(d) >= 10
      mon = int2str(mont(d));
   end
   
   for h = 0:23
  
      if h <= 9
         hour = strcat('0',int2str(h));
      elseif h >= 10
         hour = int2str(h);
      end
      
      disp(strcat(day,',',hour))
      
      % Footprint.
                                                                                                                    
      footprint_temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/stilt_output/roosevelt_april2015/footprint_05/',int2str(year(d)),mon,day,hour,'_footprint.nc'),'footprint');

      footprint = cat(4,footprint,footprint_temp);      
               
      time_model = cat(1,time_model,datenum(year(d),mont(d),days(d),(h),0,0));
         
   end
end  

% Remove first footprint (it is repeated as artifact of method).

footprint(:,:,:,1) = [];
time_model(1) = [];

% Read in emissions data.

ch4 = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/stilt_output/emissions/regrid_ch4_05.nc','ch4');
lat = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/stilt_output/emissions/regrid_ch4_05.nc','xlat');
lon = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/stilt_output/emissions/regrid_ch4_05.nc','xlon');

%% Calculate hourly concentrations.

dims = size(footprint);

sum_foot = squeeze(sum(footprint,3));

ch4_conc = zeros(1,dims(4));

for t = 1:length(ch4_conc) % Loop over every hour with backward trajs.
   
    ch4_conc(t) = sum(sum(sum_foot(:,:,t)' .* ch4));
    
end

%% Plot some CH4 time series.

figure(1)
plot(time_model,ch4_conc)

NumTicks = ((length(days)*4)+1) ;
L = get(gca,'XLim');
set(gca,'XTick',linspace(L(1),L(2),NumTicks))

datetick('x','DD, HH','keepticks')

data_rsv = readtable('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/stilt_output/emissions/data/roosevelt/2015_04_calibrated.dat');
data_fru = readtable('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/stilt_output/emissions/data/fruitland/2015_04_calibrated.dat');

% Roosevelt observations.

ts_time_rsv = table2array(data_rsv(:,'Time_UTC'));
ts_datenum_rsv = datenum(ts_time_rsv, 'yyyy-mm-dd hh:MM:SS');
ts_ch4_mat_rsv  = table2array(data_rsv(:,'CH4d_ppm_cal'));

% Fruitland observations.

ts_time_fru = table2array(data_fru(:,'Time_UTC'));
ts_datenum_fru = datenum(ts_time_fru, 'yyyy-mm-dd hh:MM:SS');
ts_ch4_fru = table2array(data_fru(:,'CH4d_ppm_cal'));

ts_time_fru_avg = mean(ts_ch4_fru(29992:end));

figure(2)
plot(ts_datenum_rsv,ts_ch4_mat_rsv)
hold on
plot(ts_datenum_fru,ts_ch4_fru,'k-')
hold on
plot(time_model,(ch4_conc+ts_time_fru_avg),'r-','LineWidth',2)
axis([min(time_model),max(time_model),0,10])

NumTicks = length(days)+1 ;
L = get(gca,'XLim');
set(gca,'XTick',linspace(L(1),L(2),NumTicks))

datetick('x','mmm DD','keepticks')

%% Average daily concentration.

ch4_conc_fin = ch4_conc + ts_time_fru_avg;

ch4_conc_24 = reshape(ch4_conc_fin,[24,12]);

ch4_conc_24_avg = mean(ch4_conc_24,2);

ts_datevec_rsv = datevec(ts_datenum_rsv);

sum = 0;
t1 = 0;
t2 = 1;

for i = 1:length(ts_datenum_rsv)-1
   if ts_datevec_rsv(i,4) == ts_datevec_rsv(i+1,4)
       sum = sum + ts_ch4_mat_rsv(i);
       t1 = t1 + 1;
   elseif ts_datevec_rsv(i,4) ~= ts_datevec_rsv(i+1,4)        
       ts_ch4_rsv_24(t2) = sum / t1;
       test(t2) = ts_datevec_rsv(i,4);
       sum = 0;
       t1 = 0;
       t2 = t2 + 1
   end
end
    
rsv_24 = ts_ch4_rsv_24(5:292);
        
rsv_obs_24 = reshape(rsv_24,[24,12]);

rsv_obs_24_avg = mean(rsv_obs_24,2);


figure(3)
plot([rsv_obs_24_avg(7:end); rsv_obs_24_avg(1:6)],'b-','LineWidth',2)
hold on
plot([ch4_conc_24_avg(7:end); ch4_conc_24_avg(1:6)],'b--','LineWidth',2)
title('CH4 April Average Diurnal Cycle: Obs vs Sim')
legend('Horsepool','Horsepool.sim')
xlabel('Hour (MST)')
ylabel('CH4 (ppm)')