% AWS GILL - T1 GoAmazon 2015-2016
% Calcula média 5 min e salva Level2
% Functions used: fmeanwind.m, export_LFA_GoAmazon_jan2015.m
% Output: alldata (avg5min)
% 1-time(utc),2-ws_vectorialmean(m/s),3-wd_vectorialmean(deg),4-press(hPa)
% 5-RH(%),6-temp(C),7-dewpoint(C),8-prt_accumulated(C),9-status_summed
% 10-u(m/s), 11-v(m/s), 12-stpfactor
%
% Parametros de entrada --------------------------------------------
caminho1='C:\A\GOAMAZON\T1\AWS_GILL\raw\'; % arquivos flex (Labview)
caminho2='C:\A\GOAMAZON\T1\AWS_GILL\';
tit='T1_2015-2016';
pref=1013.25; % Reference pressure for STP correction (mbar)
tref=293.15; % Reference temp for STP correction (K)
basicplots=false;
lev2=true;
%
% Lendo rawdata ---------------------------------------------------
files=dir(caminho1);
alldata=zeros(1,11);
disp('-----Reading rawdata files')
for ifile=3:size(files,1)
    nome=getfield(files, {ifile,1}, 'name');        
    disp(nome)    
    fid=fopen(strcat(caminho1,nome));
    textscan(fid,'%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s	%s',1,'Delimiter',',');
    C=textscan(fid,'%*f	%s	%s	%*f	%*s	%f	%*s	%f	%*s	%*f	%*s	%*f	%*s	%*f	%*s	%f	%*s	%f	%*s	%f	%*s	%f	%*s	%f	%*s	%*f	%*s	%*f	%*s	%*f	%*s	%*f	%*s	%*f	%*s	%*f	%*s	%f	%*s	%*s	%*s','Delimiter',',');
    fclose(fid);
    data_(:,2:9)=cell2mat(C(:,3:10));
    if(isnan(data_(1,2)))
        istart=2;
        data_(1,:)=[];
    else
        istart=1;
    end
    for i=istart:length(C{1,1})
        data_(i-istart+1,1)=datenum(strcat(C{1,1}(i),C{1,2}(i)),'yyyy:mm:ddHH:MM:SS');
    end
    %
    % 1-time(utc),2-ws(m/s),3-wd(deg),4-press(hPa),5-RH(%),6-temp(C),7-dewpoint(C),8-prt(C),9-status    
    tline=floor(data_(1,1)):datenum(0,0,0,0,5,0):floor(data_(1,1))+1;
    datam=nan(length(tline),size(data_,2));
    datam(:,1)=tline;
    for i=1:length(tline)-1
        aux=find(data_(:,1)>=tline(i) & data_(:,1)<tline(i+1));
        if(~isempty(aux))
            % mean press,RH,temp,dewpoint
            datam(i,4:size(data_,2)-2)=nanmean(data_(aux,4:size(data_,2)-2));
            % vectorial mean ws,wdir
            [datam(i,2),datam(i,3),datam(i,10),datam(i,11)]=fmeanwind(data_(aux,2),data_(aux,3));
            %[datam(i,3)]=circ_mean(data_(aux,3),[],1)*180/pi; % average wind direction in deg
            % sum prt,status
            datam(i,8)=nansum(data_(aux,8));
            datam(i,9)=nansum(data_(aux,9));
        end        
    end
    alldata=cat(1,alldata,datam);
    clear data_ aux datam tline C    
end
alldata(1,:)=[];
%
% Exclusões -----------------------------------------------------
origsize=length(alldata);
alldata(alldata(:,9)>0,2:12)=NaN; % exclui status>0
alldata(:,8:9)=[]; % excluí prt, pois não parece ser chuva
%
% Organiza alldata em base de tempo de 5min sem buracos
first = floor(alldata(1,1));   % first averaging time
last = datenum(2016,4,26);  % last averaging time
delta = datenum(0,0,0,0,5,0);     % 5min
tol=0; %
[~,stat1]=fmedia_qualquer_coisa2(alldata,first,last,delta,tol);
alldata=stat1(:,1:10);
%
% STP factor -----------------------------------
% Calcula ciclo diurno médio de T e P.
% Esses valores médios serão usados para substituir buracos na base de dados de STP.
%
%alldata(:,12)=(pref*(alldata(:,6)+273.15))./(tref*alldata(:,4));
tcol=6; pcol=4;
tvec=datevec(alldata(:,1));
deltat=[0 12; 12 15; 15 19; 19 24]; %night(0-12utc),morning(12-15utc),afternoon1(15-19utc),afternoon2(19-24utc)
for i=1:4
    aux=find(mod(alldata(:,1),1)>=deltat(i,1)/24 & mod(alldata(:,1),1)<deltat(i,2)/24);
    lookup(i,:)=[nanmean(alldata(aux,tcol))+273.15 nanmean(alldata(aux,pcol))]; % substituir pontos isolados sem dados
end
lookup(:,3)=(pref*lookup(:,1))./(lookup(:,2)*tref);
first = 0;   % first averaging time
last = 1;  % last averaging time
delta = datenum(0,0,0,0,5,0);     % 5min
tol=0; %
[~,stat1]=fmedia_qualquer_coisa2([mod(alldata(:,1),1),alldata(:,tcol)+273.15,alldata(:,pcol)],first,last,delta,tol);
mediatp=stat1(:,1:3); % ciclo diurno médio T,P 
mediatp(length(mediatp),:)=[];
for j=2:3
    aux=find(isnan(mediatp(:,j)));
    for i=1:length(aux)
        mediatp(aux(i),j)=mediatp(aux(i)-1,j); % substituir dias inteiros sem dados
    end
end
mediatp(:,4)=(pref*mediatp(:,2))./(mediatp(:,3)*tref); % stp mean
clear stat1 
%
% Calcula fator de correção STP, completando pontos faltantes com o ciclo diurno médio
stp=[alldata(:,1) alldata(:,tcol)+273.15 alldata(:,pcol)]; % matlabtime, temp(K), pressure
stp(:,4)=(pref*stp(:,2))./(stp(:,3)*tref); % stp when there is both T and P
%
% Substitui dias inteiros sem dados
for i=floor(stp(1,1)):1:floor(stp(length(stp),1))
    aux1=find(floor(stp(:,1))==i);
    if(sum(isnan(stp(aux1,4)))>280)
        %disp([sum(isnan(stp(aux1,4))) tvec(aux1(1),:)])
        stp(aux1,4)=mediatp(:,4);       
    end
    clear aux1
end
% Substitui pontos isolados sem dados
for i=1:length(stp)
    if(sum(isnan(stp(i,2:3)))==2) % there is no T nor P
        if(tvec(i,4)>=0 && tvec(i,4)<12) % night (0-12 utc)
            stp(i,4)=lookup(1,3);
        elseif(tvec(i,4)>=12 && tvec(i,4)<15) % morning (12-15 utc)
            stp(i,4)=lookup(2,3);
        elseif(tvec(i,4)>=15 && tvec(i,4)<19) % afternoon1 (15-19 utc)
            stp(i,4)=lookup(3,3);
        elseif(tvec(i,4)>=19 && tvec(i,4)<24) % afternoon2 (19-24 utc)
            stp(i,4)=lookup(4,3);
        end
    elseif(isnan(stp(i,2))) % there is P but no T
        if(tvec(i,4)>=0 && tvec(i,4)<12) % night (0-12 utc)
            stp(i,4)=(pref*lookup(1,1))/(stp(i,3)*tref);
        elseif(tvec(i,4)>=12 && tvec(i,4)<15) % morning (12-15 utc)
            stp(i,4)=(pref*lookup(2,1))/(stp(i,3)*tref);
        elseif(tvec(i,4)>=15 && tvec(i,4)<19) % afternoon1 (15-19 utc)
            stp(i,4)=(pref*lookup(3,1))/(stp(i,3)*tref);
        elseif(tvec(i,4)>=19 && tvec(i,4)<24) % afternoon2 (19-24 utc)
            stp(i,4)=(pref*lookup(4,1))/(stp(i,3)*tref);
        end
    elseif(isnan(stp(i,3))) % there is T but no P
        if(tvec(i,4)>=0 && tvec(i,4)<12) % night (0-12 utc)
            stp(i,4)=(pref*stp(i,2))/(lookup(1,2)*tref);
        elseif(tvec(i,4)>=12 && tvec(i,4)<15) % morning (12-15 utc)
            stp(i,4)=(pref*stp(i,2))/(lookup(2,2)*tref);
        elseif(tvec(i,4)>=15 && tvec(i,4)<19) % afternoon1 (15-19 utc)
            stp(i,4)=(pref*stp(i,2))/(lookup(3,2)*tref);
        elseif(tvec(i,4)>=19 && tvec(i,4)<24) % afternoon2 (19-24 utc)
            stp(i,4)=(pref*stp(i,2))/(lookup(4,2)*tref);
        end
    end
end
alldata(:,10)=stp(:,4);
%
% Basic plots --------------------------------------------------------
% 1-time(utc),2-ws_vectorialmean(m/s),3-wd_vectorialmean(deg),4-press(hPa)
% 5-RH(%),6-temp(C),7-dewpoint(C),8-u(m/s), 9-v(m/s), 10-stpfactor
if(basicplots)
    figure
    subplot(3,1,1),plot(alldata(:,1),alldata(:,6),'-r.'),grid on, dynamicDateTicks    
    ylabel('Temp(oC)'),axis tight
    subplot(3,1,2),plot(alldata(:,1),alldata(:,4),'-g.'),grid on, dynamicDateTicks   
    ylabel('Pressure(hPa)'),axis tight
    subplot(3,1,3),plot(alldata(:,1),alldata(:,10),'-k.'),grid on, dynamicDateTicks
    ylabel('STP factor'),axis tight
    saveas(gcf,[caminho2,tit,'_temp_press_stp.jpg'],'jpg')
    %
    figure
    subplot(3,1,1),plot(alldata(:,1),alldata(:,2),'-b.'),grid on, dynamicDateTicks    
    ylabel('WS(m/s)'),axis tight
    subplot(3,1,2),plot(alldata(:,1),alldata(:,3),'-r.'),grid on, dynamicDateTicks   
    ylabel('WDIR(deg)'),axis tight    
    subplot(3,1,3),plot(alldata(:,1),alldata(:,5),'-c.'),grid on, dynamicDateTicks    
    ylabel('RH(%)'),axis tight   
    saveas(gcf,[caminho2,tit,'_ws_wdir_RH.jpg'],'jpg')
    %
    figure
    plot(mod(alldata(:,1),1),alldata(:,10),'k.'),grid on
    ylabel('STP factor'),axis tight
    saveas(gcf,[caminho2,tit,'_stp_dcycle.jpg'],'jpg')
end
%
% ------------ Level 2 ------------------------------------------------
% 1-time(utc),2-ws_vectorialmean(m/s),3-wd_vectorialmean(deg),4-press(hPa)
% 5-RH(%),6-temp(C),7-dewpoint(C),8-u(m/s), 9-v(m/s), 10-stpfactor
if(lev2)                
    %
    Instrument='AWS';
    Model='Gill';
    serial_num=' ';
    Site=tit;
    Val_level='2';
    Proc_vers='1.0';
    labels={'WS_vector';'WDIR_vector';'Press';'RH';'Temp';'DewPoint';'u';'v';'STP(20oC,1atm)'};
    units={'(m/s)';'(deg)';'(hPa)';'(prct)';'(oC)';'(oC)';'(m/s)';'(m/s)';'(adim)'};    
    export_LFA_GoAmazon_jan2015(labels,units,Instrument,Model,serial_num,Site,Val_level,Proc_vers,origsize,alldata(:,1),alldata(:,2:size(alldata,2)))
end