%start dealing with the daily CERES surface data

clear all; close all;

dir ='D:\ceres_daily\extracted\'; %change to the directory where you keep the extracted data

fils={'CERES_SYN1deg-Day_Terra-Aqua-MODIS_Ed4A_Subset_20000301-20020902.nc';  
      'CERES_SYN1deg-Day_Terra-Aqua-MODIS_Ed4A_Subset_20020903-20050306.nc '; 
      'CERES_SYN1deg-Day_Terra-Aqua-MODIS_Ed4A_Subset_20050307-20070908.nc '; 
      'CERES_SYN1deg-Day_Terra-Aqua-MODIS_Ed4A_Subset_20070909-20100312.nc';  
      'CERES_SYN1deg-Day_Terra-Aqua-MODIS_Ed4A_Subset_20100313-20120913.nc';  
      'CERES_SYN1deg-Day_Terra-Aqua-MODIS_Ed4A_Subset_20120914-20150318.nc';  
      'CERES_SYN1deg-Day_Terra-Aqua-MODIS_Ed4A_Subset_20150319-20170919.nc';  
      'CERES_SYN1deg-Day_Terra-Aqua-MODIS_Ed4A_Subset_20170920-20171031.nc';}  
  
  daysperyear=365.25;
  dayoff = monthday(3)+1;
  off=0;
  for n=1:1:length(fils); %loop over the files
      n 
      tic 
      FIL = char([char(dir) char(fils(n,:))]);
      %ncdisp(FIL)
      t = ncread(FIL, 'time', [1], [Inf]); NT =length(t);
      lon =ncread(FIL, 'lon', [1], [Inf]); N =length(lon);
      lat = ncread(FIL, 'lat', [1], [Inf]); NN =length(lat);
      year = 2000+ rounddown((t+dayoff)/daysperyear); %rounddown should be the same as matlab's built in floor
      day = (t+dayoff)-(year-2000)*daysperyear; % day of year
      
      DAY(off+(1:NT))= day;  % make it one continous vector
      YEAR(off+(1:NT))= year; %same
      
      
      DIR(off+(1:NT),1:NN,1:N) = permute(ncread(FIL, 'adj_sfc_sw_direct_all_daily', [1 1 1] , [Inf Inf Inf]), [3 2 1]); %direct solar radiation
      DIFF(off+(1:NT),1:NN,1:N) = permute(ncread(FIL, 'adj_sfc_sw_diff_all_daily', [1 1 1] , [Inf Inf Inf]), [3 2 1]); % diffuse solar 
      
      off=off+NT; % update the index offset for the next iteration of the file loop
      
      toc
  end
  
  
lat=double(lat);
lon=double(lon);
DAY=double(DAY);

%%%% start trying to sort the surface solar phasing
 day_out=0:0.25:365; %this is the date for the smoothed seasonal cyle
 
 NT = length(DAY);
 for nt=1:NT; %loop over times
     do_smooth(nt) = find(day_out == DAY(nt)); %index fot the smoothed seasonal cycle
 end
     

for nn=1:NN; % loop over latitude
    nn
    tic
    for n=1:N; % loop over longitude
       ts = squeeze(DIR(:,nn,n)+DIFF(:,nn,n)); % the time series at a single grid point
       pp_periodic = splinefit(DAY, ts, 6, 'p');
       y3 = ppval(pp_periodic, day_out);
       %%%%%%%%%%%%%%%%%%%%%%%%%%%% deseasonalize the data by removing the
       %%%%%%%%%%%%%%%%%%%%%%%%%%%% smooth seasonal cycle (y3) at the DAY
       %%%%%%%%%%%%%%%%%%%%%%%%%%%% whic is given by index do_smooth
       
       aDS(1:NT,nn,n) = ts'-y3(do_smooth); %%%!!! remove the smoothes seasonal cycle -- this if for the total downwelling shortwave at surface, you could do individually for the DIR and DIFF 
       
       
       %% dates of max and min surface insolation (this was leftover from what I was doing prstd(aDSeviously
       MAX_DAY(nn,n) = day_out(max(find(y3==max(y3))));  %% 
       MIN_DAY(nn,n) = day_out(min(find(y3==min(y3))));
    end
    toc
end

clear DIR DIFF %clear up some spcae
save('CERES_deseas.mat', 'lat', 'lon','aDS', 'DAY', 'YEAR') 
