%function isfs_fluxes
addpath(genpath('/home/bbutterworth/mfiles/'))

clear
fl=listfiles('/air/incoming/CHEESEHEAD/ISFS/high-rate','isfs_geo');
%fl=listfiles('/Users/brian/Downloads','isfs_geo');

sf = 20;%sample frequency (Hz)
Nmin_flux = 30;%length of flux interval
Nsamp_per_int = Nmin_flux * 60 * sf;%number of samples per flux interval

% Specific heat
% Dry air Cp - not a moist air Cp
Cpd = 1004.67 ; %  J/Kg K;
Cpv = 1875;     %  J/Kg K                               
Cp = Cpd ;  %mean(Cpd*(ones(size(Xw))-Xw) + Cpv.*Xw)
                 
% conversion from umol/m2/s to mol/m2/yr
% fc_factor = 24 * 3600 * 365.25 / 1e6;

%initialize variables     
wind_stats=[]; 
hsensible_stats=[];
hlatent_stats=[];
fco2_stats=[];

u=NaN(2*length(fl(:,1)),17);   
v=NaN(2*length(fl(:,1)),17);   
ws=NaN(2*length(fl(:,1)),17);   
wd=NaN(2*length(fl(:,1)),17);   
wind_uw=NaN(2*length(fl(:,1)),17);   
wind_vw=NaN(2*length(fl(:,1)),17);   
ustar=NaN(2*length(fl(:,1)),17);   

% create variable names
uheader={['u_30m_nw1'];
['u_10m_nw2'];
['u_2m_nw3'];
['u_30m_nw4'];
['u_30m_ne1'];
['u_30m_ne2'];
['u_30m_ne3'];
['u_30m_ne4'];
['u_30m_sw1'];
['u_25m_sw2'];
['u_30m_sw3'];
['u_30m_sw4'];
['u_30m_se2'];
['u_30m_se3'];
['u_3m_se4'];
['u_10m_se5'];
['u_30m_se6'];};

for j=1:length(uheader);
    vheader{j}=['v' uheader{j}(2:end)];
    wheader{j}=['w' uheader{j}(2:end)];
    h2oheader{j}=['h2o' uheader{j}(2:end)];
    co2header{j}=['co2' uheader{j}(2:end)];
    Theader{j}=['tc' uheader{j}(2:end)];
    Pheader{j}=['P' uheader{j}(2:end)];

    sites{j}=uheader{j}(end-2:end);
end

% ==== Main loop ====
for n=1:length(fl(:,1))

    
    for j=1:length(uheader)
    
    eval([ uheader{j} '=ncread(' '''' fl(n,:) '''' ',' '''' uheader{j}  '''' ');'])
    eval([ vheader{j} '=ncread(' '''' fl(n,:) '''' ',' '''' vheader{j}  '''' ');'])
    eval([ wheader{j} '=ncread(' '''' fl(n,:) '''' ',' '''' wheader{j}  '''' ');'])
    eval([ h2oheader{j} '=ncread(' '''' fl(n,:) '''' ',' '''' h2oheader{j}  '''' ');'])
    eval([ co2header{j} '=ncread(' '''' fl(n,:) '''' ',' '''' co2header{j}  '''' ');'])
    eval([ Theader{j} '=ncread(' '''' fl(n,:) '''' ',' '''' Theader{j}  '''' ');'])
    eval([ Pheader{j} '=ncread(' '''' fl(n,:) '''' ',' '''' Pheader{j}  '''' ');'])
    Pheader{15}='P_2m_se4';%NCAR had a typo

    time=ncread(fl(n,:),'time');
       
    ts1(n)=datenum(str2num(fl(n,end-13:end-10)),...
       str2num(fl(n,end-9:end-8)),...
       str2num(fl(n,end-7:end-6)),...
       str2num(fl(n,end-4:end-3)),0,0);

    ts=datenum(str2num(fl(n,end-13:end-10)),...
       str2num(fl(n,end-9:end-8)),...
       str2num(fl(n,end-7:end-6)),...
       str2num(fl(n,end-4:end-3)),0,0:1/sf:3600);
    ts=ts(1:end-1);  

    ts1hz=datenum(str2num(fl(n,end-13:end-10)),...
       str2num(fl(n,end-9:end-8)),...
       str2num(fl(n,end-7:end-6)),...
       str2num(fl(n,end-4:end-3)),0,time);
   
    % Reshape matrixes from 3600x20 to 72000x1:
    ux=reshape(eval(uheader{j}),[],1);  
    uy=reshape(eval(vheader{j}),[],1);       
    uz=reshape(eval(wheader{j}),[],1);      % vertical wind
    xh2o=reshape(eval(h2oheader{j}),[],1);  % h2o mixing ratio
    xco2=reshape(eval(co2header{j}),[],1);  % co2 mixing ratio
    T=reshape(eval(Theader{j}),[],1);       % temperature
    P=reshape(eval(Pheader{j}),[],1);       % pressure (most of the site were recorded at 1Hz; not sw2)
    RH=reshape(eval(Theader{j}),[],1);      % relative humidity

    % Gap fill NaNs:
    [ux,nan01_U]=nanfill(ux');
    [uy,nan01_V]=nanfill(uy');
    [uz,nan01_W]=nanfill(uz');
    [xh2o,nan01_H2O]=nanfill(xh2o');
    [xco2,nan01_CO2]=nanfill(xco2');
    [T,nan01_T]=nanfill(T');
    [P,nan01_P]=nanfill(P');
    [RH,nan01_RH]=nanfill(RH');
    
    % Coordinate Rotation:
    um=[ux;uy;uz];
    uc=coordrot(um,0);
    
    Td=(T+273.15)/(1+0.32*nanmean(1e-3*xh2o));
      
    if numel(find(isnan(xco2)==1))<36000
         xco2=((62.4.*(Td).*(xco2))./nanmean(P.*0.750062.*44.01));%Right now is using hourly mean P. Might want to use 30-minute mean. But it would require some rearrangement of the code.
    end
            


     % calculate a reference frequency axis (stop at 1 Hz to save memory)
     ptmp=hspec_ss([1:Nsamp_per_int],0,sf);
     fref=ptmp(1:length(ptmp)); 

     for ii=1:2
        if ii==1;
            i20hz=find( ts>=ts1(n) & ts<ts1(n) + (1/48) );
            i1hz=find( ts1hz>=ts1(n) & ts1hz<ts1(n) + (1/48) );
        else
            i20hz=find( ts>=ts1(n) + (1/48)  & ts<ts1(n) + (1/24) );
            i1hz=find( ts1hz>=ts1(n) + (1/48)  & ts1hz<ts1(n) + (1/24) );
        end

        iflux=n*2-2+ii;

        tmaster(iflux)=ts1(n) + (ii/48 - 1/48);

        % Median values:
        u(iflux,j)=nanmedian(ux(i20hz));
        v(iflux,j)=nanmedian(uy(i20hz));
        w(iflux,j)=nanmedian(uz(i20hz));
        H2O_mean(iflux,j)=nanmedian(xh2o(i20hz));
        CO2_mean(iflux,j)=nanmedian(xco2(i20hz));
        T_mean(iflux,j)=nanmedian(T(i20hz));
        P_mean(iflux,j)=nanmedian(P(i1hz));%1hz
        RH_mean(iflux,j)=nanmedian(RH(i20hz));
        Td_mean(iflux,j)=nanmedian(Td(i20hz));
        
        % Percentage of NaNs in the flux intervals:
        nanU(iflux,j)=sum(nan01_U(i20hz))/numel(i20hz);
        nanV(iflux,j)=sum(nan01_V(i20hz))/numel(i20hz);
        nanW(iflux,j)=sum(nan01_W(i20hz))/numel(i20hz);
        nanH2O(iflux,j)=sum(nan01_H2O(i20hz))/numel(i20hz);
        nanCO2(iflux,j)=sum(nan01_CO2(i20hz))/numel(i20hz);
        nanT(iflux,j)=sum(nan01_T(i20hz))/numel(i20hz);
        nanP(iflux,j)=sum(nan01_P(i1hz))/numel(i1hz);%1hz
        nanRH(iflux,j)=sum(nan01_RH(i20hz))/numel(i20hz);

        % Wind Direction calculation:
        ws(iflux,j)=sqrt(u(iflux,j)^2+v(iflux,j)^2);
        wd(iflux,j)=-atan2(v(iflux,j),u(iflux,j))*180/pi+270;
        if wd(iflux,j)<0; wd(iflux,j)=wd(iflux,j)+360;end
        if wd(iflux,j)>360; wd(iflux,j)=wd(iflux,j)-360;end
        
        % Air Density
         rho(iflux,j)=1e2/8.314*P_mean(iflux,j) * (1-nanmean(1e-3* xh2o(i20hz))) / Td_mean(iflux,j); 

        % Lv - Latent heat of vaporization (J/Kg) from Stull p 641
         Lvap(iflux,j)=(2.501-0.00237*(Td_mean(iflux,j)-273.15))*10^6;

         
        % ---------------------------------------------------------------
        %                       Flux Calculations
        % ---------------------------------------------------------------
        
        % Tau - Momentum Flux:    
         wind_stats=[];
         [wind_stats,wind_header]=...
          velocitystats(uc(:,i20hz),wind_stats,'sonic',[],...
          sf,fref,1,'detrend');

        wind_uw(iflux,j)=wind_stats(8);
        wind_vw(iflux,j)=wind_stats(9);
        ustar(iflux,j)=sqrt(sqrt(wind_stats(8).^2+wind_stats(9).^2));
        sigmav(iflux,j)=nanstd(uc(2,i20hz));% Needed for flux footprint

        % Hs - Sensible Heat Flux:
        [hsensible_stats,hsensible_header]=...
        scalarstats(uc(3,i20hz),28.966/1000*rho(iflux,j)*Cp*Td(i20hz),hsensible_stats,...
        'sonic',[],sf,fref,iflux,'detrend');
        hs(iflux,j)=hsensible_stats(5,iflux);

        % Hl - Latent Heat Flux:
        [hlatent_stats,hlatent_header]=...
        scalarstats(uc(3,i20hz),18/1000/1000*Lvap(iflux,j)*rho(iflux,j)*xh2o(i20hz),...
        hlatent_stats,'sonic',[],sf,fref,iflux,'detrend');
        hl(iflux,j)=hlatent_stats(5,iflux);

        % FCO2 - CO2 Flux:
        [fco2_stats,fco2_header]=...
        scalarstats(uc(3,i20hz),rho(iflux,j)*xco2(i20hz),[],...
        'sonic',[],sf,fref,[],'detrend');
        fco2(iflux,j)=fco2_stats(5,:);
       end

    end
   
   disp([num2str(n) ' - ' datestr(ts(n))])
   
   savepath='/home/bbutterworth/proc/isfs_fluxes_testX.mat';
   if rem(iflux,400)==0
      save(savepath);
   end
   
end

h2o=H2O_mean;clear H2O_mean;
co2=CO2_mean;clear CO2_mean;
T=T_mean;clear T_mean;
P=P_mean;clear P_mean;
RH=RH_mean;clear RH_mean;
Td=Td_mean;clear Td_mean;

save(savepath);


% Copy mfiles used to process the fluxes to archive
% Create a run_id based on date code was processes

mfilepath = [mfilename('fullname') '.m'];
[fList,pList] = matlab.codetools.requiredFilesAndProducts(mfilepath);
[d1 d2 d3]=datevec(date);run_id=[sprintf('%04d',d1) '_' sprintf('%02d',d2) '_' sprintf('%02d',d3) ];
archive_path=['/home/bbutterworth/proc/archive/' run_id '/'];
    if ~exist(archive_path, 'dir')
       mkdir(archive_path)
    end

for k=1:length(fList)
    k
    strbeg=strfind(fList{k},'/');strbeg=strbeg(end);strend=strfind(fList{k},'.m');
filename=fList{k}(strbeg+1:strend-1);
copyfile(fList{k},[archive_path filename '_' run_id '.m']);

end



% clean the output to include only the 30 minute variables (not the 20hz):
varnames=who;
cnt=0;
for j=1:length(varnames)
    len1=length(ustar(:,1));
    len2=length(ustar(1,:));
    if ~isempty(eval( varnames{j}))
    if length(eval([ varnames{j} '(:,1)']))==len1 |...
            length(eval([ varnames{j} '(1,:)']))==len1 |...
            length(eval([ varnames{j} '(:,1)']))==len2 |...
            length(eval([ varnames{j} '(1,:)']))==len2 
        cnt=cnt+1;
        varkeep{cnt}=varnames{j}
    end
    end
end
keep(varkeep{:},'run_id')

clear co2header h2oheader Theader hlatent_stats hsensible_stats wheader vheader uheader Pheader 
save('/home/bbutterworth/proc/isfs_fluxes_trimX.mat');
