%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% This program calculates the mean hourly values %%%%
%%%% and the diurnal cycle of all variables %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Estimate cloud fraction and its diurnal cycle %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%keep cr
%clear all;
%warning off;

        
ceilo_path = {'C:\THESIS\DATA\EPIC\data\ceilo\ETL_web\epic2_ceilo_15s.txt';'C:\THESIS\DATA\PACS03\DATA\ceilo\STRATUS03_ceilo_30s_a.txt';'C:\THESIS\DATA\Stratus04\ceilo\STRATUS04_ceilo_30s_a.txt'};
jdc = [273 304 335];


load C:\THESIS\Processing_Results\Stratus04\combined\stratus04_clfr_correction

%cr(3).ceilo.hr.mday      = cf_mday;
cr(3).ceilo.hr.clfr       = cf_hr_corr'/100;
cr(3).ceilo.hr.clfr_init  = cf_hr_init'/100;
cr(3).ceilo.diu.clfr      = cf_diu_corr/100;
cr(3).ceilo.diu.clfr_init = cf_diu_init/100;

cr(3).ceilo.hr.clfr(380) = 0;


keep cr ceilo_path jdc

for i=1:2    %%%% 1 corresponds to EPIC, 2 to Stratus03 - Stratus04 (i=3) is being processed separately above the loop due to the need of the cloud fraction correction 
             %%%% (permanent mat file created in the separate code "all_hourly_diurnal_mat")   
    
    ceilo_data = load(char(ceilo_path(i)));

    lam = 0:24;  		  %%%% hours
    jdx = lam/24;		  %%%% fraction of day

    jdd = ceilo_data(:,1);  %%%% date vector
    jcc = ceilo_data(:,2);  %%%% cloud code vector
    
    jd0 = floor(jdd(1));   %%%% first Julian Day of data series 
    jxm = floor(jdd(end)); %%%% last Julian Day of data series   
    jx=1;				   %%%% counts number of hours of data

    while jd0<=jxm        %%%% the current day is l.t.e. last day
       ii = find(jdd>jd0 & jdd<(jd0+1));  %%%% find all data in the current day

       jj=1;		      %%%% start first hour calc.
       while jj<25        %%%% not the end of the day
            fpd = jdd(ii)-floor(jdd(ii));           %%%% fractional part of day
            ij = find(fpd>jdx(jj) & fpd<jdx(jj+1)); %%%% select hour
            
            alldata(jx) = length(fpd(ij));		    %%%% all data in selected hour            
            iscloud     = find(jcc(ii(ij))>0 & jcc(ii(ij))<5);   %%%% Consider the 3 cloud layers + obscure (correct, especially for Stratus04) 
            cloud(jx)   = length(iscloud);
            cf_hr(jx)   = cloud(jx)/alldata(jx);
            
            %iscloud2   = find(jcc(ii(ij))>0 & jcc(ii(ij))<4);  %%%% Consider only the 3 cloud layers                                  
            %cloud2(jx) = length(iscloud2);
            %clfr2(jx) = cloud2(jx)/alldata(jx);
            
            %israin = find(ceilo1(ii(ij),3)<=500);    % Consider 500m as threshold for precipation (only cloud layer 1) 
            %rain(jx) = length(israin);
            %rainfr_ovcl(jx)  = rain(jx)/cloud(jx);
            %rainfr_ovall(jx) = rain(jx)/alldata(jx);

            jdh(jx) = jd0+jdx(jj);
            jj = jj+1;
            jx = jx+1;
        end;        %end while jj
        jd0 = jd0+1;
    end;            

    %rainfr_ovcl(find(isnan(rainfr_ovcl)))=0;
    %rainfr_ovcl(find(isinf(rainfr_ovcl)))=0;
    %rainfr_ovall(find(isnan(rainfr_ovall)))=0;
    %rainfr_ovall(find(isinf(rainfr_ovall)))=0;
    %avg_rf_ovcl=mean(rainfr_ovcl,2);
    %avg_rf_ovall=mean(rainfr_ovall,2);
    
   
    %cr(i).ceilo.hr.mday  = jdh-jdc(i);
    cr(i).ceilo.hr.clfr  = cf_hr';
    cr(i).ceilo.diu.clfr = nanmean(reshape(cf_hr,24,[]),2);
   
    keep cr i ceilo_path jdc
end


%%%% Addition: Mean & Std Values for Cloud Fraction %%%%

ds = [16 16 11.2]; 
de = [22 21 16.2];

for i=1:3

    ind_cf = find(cr(i).ceilo.hr.mday>=ds(i) & cr(i).ceilo.hr.mday<de(i));
    
    if i==3
        cr(i).ceilo.hr.mean.clfr = nanmean(cr(i).ceilo.hr.clfr(ind_cf));
        cr(i).ceilo.hr.std.clfr  = nanstd(cr(i).ceilo.hr.clfr(ind_cf));
        cr(i).ceilo.hr.mean.clfr_init = nanmean(cr(i).ceilo.hr.clfr_init(ind_cf));
        cr(i).ceilo.hr.std.clfr_init  = nanstd(cr(i).ceilo.hr.clfr_init(ind_cf));
    else
        cr(i).ceilo.hr.mean.clfr = nanmean(cr(i).ceilo.hr.clfr(ind_cf));
        cr(i).ceilo.hr.std.clfr  = nanstd(cr(i).ceilo.hr.clfr(ind_cf));
    end
    clear ind_cf
end


keep cr


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Estimate drizzle fraction and its diurnal cycle %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%clear all

radar_path = {'C:\THESIS\Processing_Results\EPIC2001\MMCR\EPIC_mmcr_reader.mat';'C:\THESIS\Processing_Results\Stratus03\MMCR\PACS_mmcr_reader.mat';'C:\THESIS\Processing_Results\Stratus04\radars\stratus04_fmcw_drfr_-10dBZ.mat'};
jdc = [273 304 335];

i=3;
load(char(radar_path(i)));

%cr(i).radar.hr.mday  = df_mday;
cr(i).radar.hr.drfr  = [df_hr (1:24)*NaN];
cr(i).radar.diu.drfr = df_diu;

clear df* i

for i=1:2    %%%% 1 corresponds to EPIC and 2 to Stratus03 - The Stratus04 respective variables (i=3) are already calculated and are just added to the structure above the loop
    
    cr(i).radar.data = load(char(radar_path(i)));
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% Get Max Refl of -10 dBZ as threshold for drizzle fraction %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if i==1
        dum = find(cr(i).radar.data.dbz<-35); %%%% Masking for EPIC
        cr(i).radar.data.dbz(dum) = NaN;
        cr(i).radar.data.dbz_max = max(cr(i).radar.data.dbz(:,3:40)');
    elseif i==2
        ks = find(cr(i).radar.data.snr<-5);   %%%% Masking for PACS
        cr(i).radar.data.dbz(ks)= NaN;
        cr(i).radar.data.dbz_max = max(cr(i).radar.data.dbz(:,5:40)');
    end

    lam = 0:24;  		  % hours
    jdx = lam/24;		  % fraction of day

    jdd  = cr(i).radar.data.time;	   % date vector
    jdbz = cr(i).radar.data.dbz_max;    % Max dBz vector 

    jd0 = floor(jdd(1));   % first Julian Day of data series  % 340 for Stratus04
    jxm = floor(jdd(end)); % last Julian Day of data series   % 355 -//-
    jx=1;				   % counts number of hours of data

    while jd0<=jxm        % the current day is l.t.e. last day
       ii = find(jdd>jd0 & jdd<(jd0+1));  % find all data in the current day

       jj=1;		      % start first hour calc.
       while jj<25        % not the end of the day
            fpd = jdd(ii)-floor(jdd(ii));           % fractional part of day
            ij = find(fpd>jdx(jj) & fpd<jdx(jj+1)); % select hour
            
            alldata(jx) = length(fpd(ij));		   % all data in selected hour            
            isdrizzle   = find(jdbz(ii(ij))>-10);   % Define the threshold for drizzle
            drizzle(jx) = length(isdrizzle); 
            df_hr(jx)   = drizzle(jx)/alldata(jx);
            
            %isdrizzle2 = find(jdbz(ii(ij))>-17);   % Define the threshold for drizzle
            %drizzle2(jx) = length(isdrizzle2);    
            %drfr2(jx) = drizzle2(jx)/alldata(jx);

            jdh(jx) = jd0+jdx(jj);
            jj = jj+1;
            jx = jx+1;
        end;        %end while jj
        jd0 = jd0+1;
    end;            
    
    %cr(i).radar.hr.mday  = jdh-jdc(i);
    
    if i==1
        cr(i).radar.hr.drfr  = df_hr;
    elseif i==2
        cr(i).radar.hr.drfr  = [(1:48)*NaN df_hr];
    end

    cr(i).radar.diu.drfr = nanmean(reshape(df_hr,24,[]),2);
    
    keep cr i radar_path jdc 
end

keep cr

%%%% Addition: Mean & Std Values for Drizzle Fraction %%%%

ds = [16 16 11.2]; 
de = [22 21 16.2];

for i=1:3

    ind_df = find(cr(i).radar.hr.mday>=ds(i) & cr(i).radar.hr.mday<de(i));
     
    cr(i).radar.hr.mean.drfr = nanmean(cr(i).radar.hr.drfr(ind_df));
    cr(i).radar.hr.std.drfr  = nanstd(cr(i).radar.hr.drfr(ind_df));
   
    clear ind_df
end

keep cr


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Diurnal cycle of Cloud Base, Cloud Top and Cloud Thickness %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


for i=1:3
    cr(i).ceilo.hr.clba(find(cr(i).ceilo.hr.clba<300 | cr(i).ceilo.hr.clba>1650)) = NaN;
    cr(i).ceilo.diu.clba  = nanmean(reshape(cr(i).ceilo.hr.clba',24,[]),2);
    
    cr(i).radar.diu.cltop = nanmean(reshape(cr(i).radar.hr.cltop,24,[]),2);
    cr(i).radar.diu.clthi = nanmean(reshape(cr(i).radar.hr.clthi,24,[]),2);    
end


%%%% Using NaNs for the hours of the ceimometer malfunction %%%%

cr(3).ceilo.diu.clba_init = cr(3).ceilo.diu.clba;
cr(3).ceilo.diu.clba(14:23) = NaN;

cr(3).radar.diu.clthi_init = cr(3).radar.diu.clthi;
cr(3).radar.diu.clthi(14:23) = NaN;

%{
%%%% Correction for the hours of the ceilometer malfunction %%%%

ind=[];
i=5;
for m=1:16
    if (i==6 || i==7 || i==8)
        kl=find((cr(3).ceilo.hr.mday'>=(i+0.6666) & cr(3).ceilo.hr.mday'<=(i+0.75)));
    elseif (i==5 || i==9 || i==10)
        kl=[];    
    elseif (i==11 || i==12)
        kl=find((cr(3).ceilo.hr.mday'>=(i+0.6666) & cr(3).ceilo.hr.mday'<=(i+0.8334)));
    elseif (i==13 || i==14 || i==15 || i==16 || i==17 || i==18 || i==19 || i==20)
        kl=find((cr(3).ceilo.hr.mday'>=(i+0.5416) & cr(3).ceilo.hr.mday'<=(i+0.9167)));
    end
   ind=[ind kl];
   i=i+1;
   clear kl
end

cr(3).ceilo.diu.clba_init = cr(3).ceilo.diu.clba;

clba_dummy = cr(3).ceilo.hr.clba';
clba_dummy(ind) = NaN;
cr(3).ceilo.diu.clba = nanmean(reshape(clba_dummy,24,[]),2);

cr(3).radar.diu.clthi_init = cr(3).radar.diu.clthi;

clthi_dummy = cr(3).radar.hr.clthi;
clthi_dummy(ind) = NaN;
cr(3).radar.diu.clthi = nanmean(reshape(clthi_dummy,24,[]),2);

keep cr
%}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Estimate hourly average/diurnal cycle of flux data %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%clear all;
%warning off;

flux_path = {'C:\THESIS\DATA\EPIC\data\flux\5-min\flux5_epic2_nohd.txt';'C:\THESIS\DATA\PACS03\DATA\final\flux_5hf_weller03_nohd.txt';'C:\THESIS\DATA\Stratus04\flux\flux_5hf_stratus_04a_nohd.txt'};
jdc = [273 304 335];


for i=1:3    %%%% 1 corresponds to EPIC, 2 to Stratus03 and 3 to Stratus04
    
    fx = load(char(flux_path(i)));

    wspd = fx(:,2); 
    wdir = fx(:,3);
    sst  = fx(:,4);
    tair = fx(:,7);
    tdif = sst-tair;
    qsea = fx(:,8);
    qair = fx(:,9);
    sw   = fx(:,10);
    lw   = fx(:,11);
    sh   = fx(:,22);
    lh   = fx(:,23);
    vh   = cr(i).flux.min.vh;

    tlcl = cr(i).flux.min.tlcl;
    zlcl = cr(i).flux.min.zlcl;
    
    lam = 0:24;  		   % hours
    jdx = lam/24;		   % fraction of day

    jdd = fx(:,1);	       % date vector 
    jd0 = floor(jdd(1));   % first Julian Day of data series  % 340 for Stratus04
    jxm = floor(jdd(end)); % last Julian Day of data series   % 355 -//-
    jx=1;				   % counts number of hours of data

    while jd0<=jxm         % the current day is l.t.e. last day
       ii = find(jdd>jd0 & jdd<(jd0+1));  % find all data in the current day

       jj=1;		      % start first hour calc.
       while jj<25        % not the end of the day
            fpd = jdd(ii)-floor(jdd(ii));           % fractional part of day
            ij = find(fpd>jdx(jj) & fpd<jdx(jj+1)); % select hour

            cr(i).flux.hr.wspd(jx) = nanmean(wspd(ii(ij)));
            cr(i).flux.hr.wdir(jx) = nanmean(wdir(ii(ij)));
            cr(i).flux.hr.sst(jx)  = nanmean(sst(ii(ij)));
            cr(i).flux.hr.tair(jx) = nanmean(tair(ii(ij)));
            cr(i).flux.hr.tdif(jx) = nanmean(tdif(ii(ij)));
            cr(i).flux.hr.qsea(jx) = nanmean(qsea(ii(ij)));
            cr(i).flux.hr.qair(jx) = nanmean(qair(ii(ij)));
            cr(i).flux.hr.sw(jx)   = nanmean(sw(ii(ij)));
            cr(i).flux.hr.lw(jx)   = nanmean(lw(ii(ij)));
            cr(i).flux.hr.sh(jx)   = nanmean(sh(ii(ij)));
            cr(i).flux.hr.lh(jx)   = nanmean(lh(ii(ij)));
            cr(i).flux.hr.vh(jx)   = nanmean(vh(ii(ij)));
            
            cr(i).flux.hr.tlcl(jx) = nanmean(tlcl(ii(ij)));
            cr(i).flux.hr.zlcl(jx) = nanmean(zlcl(ii(ij)));
            
            jdh(jx) = jd0+jdx(jj);
            jj = jj+1;
            jx = jx+1;
        end;        %end while jj
        jd0 = jd0+1;
    end;
    
    %%%% Estimation of Potential Cloud Depth and Decoupling Index %%%%
    
    if i==3
        cr(i).radar.hr.pcd = cr(i).radar.hr.cltop - cr(i).flux.hr.zlcl(25:408);
    else
        cr(i).radar.hr.pcd = cr(i).radar.hr.cltop - cr(i).flux.hr.zlcl;
    end
    
    cr(i).radar.hr.di = cr(i).radar.hr.clthi./cr(i).radar.hr.pcd;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    cr(i).flux.hr.mday = jdh-jdc(i);   
 
    cr(i).flux.diu.wspd = nanmean(reshape(cr(i).flux.hr.wspd,24,[]),2);
    cr(i).flux.diu.wdir = nanmean(reshape(cr(i).flux.hr.wdir,24,[]),2);
    cr(i).flux.diu.sst  = nanmean(reshape(cr(i).flux.hr.sst,24,[]),2);
    cr(i).flux.diu.tair = nanmean(reshape(cr(i).flux.hr.tair,24,[]),2);
    cr(i).flux.diu.tdif = nanmean(reshape(cr(i).flux.hr.tdif,24,[]),2);
    cr(i).flux.diu.qsea = nanmean(reshape(cr(i).flux.hr.qsea,24,[]),2);
    cr(i).flux.diu.qair = nanmean(reshape(cr(i).flux.hr.qair,24,[]),2);
    cr(i).flux.diu.sw   = nanmean(reshape(cr(i).flux.hr.sw,24,[]),2);
    cr(i).flux.diu.lw   = nanmean(reshape(cr(i).flux.hr.lw,24,[]),2);
    cr(i).flux.diu.sh   = nanmean(reshape(cr(i).flux.hr.sh,24,[]),2);
    cr(i).flux.diu.lh   = nanmean(reshape(cr(i).flux.hr.lh,24,[]),2);
    cr(i).flux.diu.vh   = nanmean(reshape(cr(i).flux.hr.vh,24,[]),2);
    
    cr(i).flux.diu.tlcl = nanmean(reshape(cr(i).flux.hr.tlcl,24,[]),2);
    cr(i).flux.diu.zlcl = nanmean(reshape(cr(i).flux.hr.zlcl,24,[]),2);
    
    cr(i).radar.diu.pcd = nanmean(reshape(cr(i).radar.hr.pcd,24,[]),2);
    cr(i).radar.diu.di  = nanmean(reshape(cr(i).radar.hr.di,24,[]),2);

    keep cr i flux_path jdc   %%%% for 2004, values (25:408) should be considered to match with CF values
end

keep cr


%%%% Addition: Mean & Std Values for Potential Cloud Depth and Decoupling Index %%%%

ds = [16 16 11.2]; 
de = [22 21 16.2];

for i=1:3

    ind_pcd = find(cr(i).radar.hr.mday>=ds(i) & cr(i).radar.hr.mday<de(i));
     
    cr(i).radar.hr.mean.pcd = nanmean(cr(i).radar.hr.pcd(ind_pcd));
    cr(i).radar.hr.mean.di  = nanmean(cr(i).radar.hr.di(ind_pcd));
    
    cr(i).radar.hr.std.pcd  = nanstd(cr(i).radar.hr.pcd(ind_pcd));
    cr(i).radar.hr.std.di   = nanstd(cr(i).radar.hr.di(ind_pcd));
    
    clear ind_pcd
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

keep cr

    
    
    


