%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Find mean and Std values for All properties %%%%%%%%%%%%%%
%%%% during Buoy Periods (Use permanent mat files) %%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all;
close all;
warning off;


%%%%%%%%%%%%%%%%%%%
%%%% soundings %%%%
%%%%%%%%%%%%%%%%%%%


sound_path = {'C:\THESIS\Processing_Results\EPIC2001\sondes\epic_sondes_reader'; 'C:\THESIS\Processing_Results\Stratus03\sondes\stratus03_sondes_reader';
              'C:\THESIS\Processing_Results\Stratus04\sondes\stratus04_sondes_reader'};

cr_name  = {'EPIC 2001'; 'Stratus 2003'; 'Stratus 2004'};
cr_month = {'October'; 'November'; 'December'};
jdc = [273 304 335];    %%%% correction coefficent to get day of month from Julian day for each cruise

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

    cr(i).name  = char(cr_name(i));
    cr(i).month = char(cr_month(i));
    
    if i==1
        cr(i).sound.data.mday = mday+0.0417; %%%% Correction for the soundings time for EPIC - all soundings are apparently launched 1 hour ahead
    else
       cr(i).sound.data.mday = mday;
    end
    
    t   = t_int;       clear t_int;  
    rh  = rh_int;      clear rh_int;
    p   = press_int;   clear press_int;
    h   = h_int;       clear h_int;
    spd = spd_int;     clear spd_int;
    dir = dir_int;     clear dir_int;

    es        = 6.112*exp((17.67*t)./(t+243.5)); %%% mb
    e         = (rh.*es)/100;                    %%% mb
    r         = 622*(e./(p-e));                  %%% g/kg
    rs        = 622*(es./(p-es));                %%% g/kg   
    tk        = t+273.15;                        %%% K
    a         = 1./(tk-55);
    b         = log(rh/100)/2480;
    tl        = (1./(a-b))+55;                                              %%% K
    theta     = tk.*((1000./p).^(0.2854*(1-0.28*10^(-3)*r)));               %%% K
    thetav    = theta.*(1+0.608*(r/1000));                                  %%% K
    thetae    = theta.*exp(((3.376./tl)-0.00254).*r.*(1+0.81*10^(-3)*r));   %%% K
    thetaes   = theta.*exp(2.675*rs./tk);                                   %%% K
    
    clear U_int V_int a b dew_int e es h_all lat lon rs td_int time_all tl;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% Incorporation of the sounding data to the "CR" structure %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    cr(i).sound.data.h       = h;     
    cr(i).sound.data.p       = p;
    cr(i).sound.data.t       = t;
    cr(i).sound.data.theta   = theta;
    cr(i).sound.data.thetav  = thetav;
    cr(i).sound.data.thetae  = thetae;
    cr(i).sound.data.thetaes = thetaes;
    cr(i).sound.data.r       = r;
    cr(i).sound.data.rh      = rh;
    cr(i).sound.data.wspd    = spd;
    cr(i).sound.data.wdir    = dir;
        
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
        %%%% Estimation of Surface and 700 mb values %%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        cr(i).sound.level.mday = cr(i).sound.data.mday;
        
        %%%% 1000 mb values %%%%
        
        for w = 1:length(mday)
            
            pr_1000_dummy = find(abs(p(:,w)-1000)==min(abs(p(:,w)-1000)));
            
            if length(pr_1000_dummy)==0  %%%% Correction for NaN sondes
       
                cr(i).sound.level.mb1000.t(w)       = nan;
                cr(i).sound.level.mb1000.theta(w)   = nan;
                cr(i).sound.level.mb1000.thetav(w)  = nan;
                cr(i).sound.level.mb1000.thetae(w)  = nan;
                cr(i).sound.level.mb1000.thetaes(w) = nan;
                cr(i).sound.level.mb1000.r(w)       = nan;
                cr(i).sound.level.mb1000.rh(w)      = nan;
                cr(i).sound.level.mb1000.wspd(w)    = nan;
                cr(i).sound.level.mb1000.wdir(w)    = nan;
                continue    
            elseif length(pr_1000_dummy)==2     %%%% Correction if the difference of the values right above and below 700 mb are exactly the same 
                pr_1000_ind = pr_1000_dummy(1);
            else
                pr_1000_ind = pr_1000_dummy;
            end
            
            cr(i).sound.level.mb1000.t(w)       = tk(pr_1000_ind,w);
            cr(i).sound.level.mb1000.theta(w)   = theta(pr_1000_ind,w);
            cr(i).sound.level.mb1000.thetav(w)  = thetav(pr_1000_ind,w);
            cr(i).sound.level.mb1000.thetae(w)  = thetae(pr_1000_ind,w);
            cr(i).sound.level.mb1000.thetaes(w) = thetaes(pr_1000_ind,w);
            cr(i).sound.level.mb1000.r(w)       = r(pr_1000_ind,w);
            cr(i).sound.level.mb1000.rh(w)      = rh(pr_1000_ind,w);
            cr(i).sound.level.mb1000.wspd(w)    = spd(pr_1000_ind,w);
            cr(i).sound.level.mb1000.wdir(w)    = dir(pr_1000_ind,w);
            
        end
        clear w;
        
        %%%% Values at 700 mb %%%%
        
        for w = 1:length(mday) 
            
            pr_700_dummy = find(abs(p(:,w)-700)==min(abs(p(:,w)-700)));
            
            if length(pr_700_dummy)==0 || (i==1 && (w==4 || w==14 || w==15 || w==81)) || (i==2 && w==37) %%%% EPIC Sondes # 4,14,15,81 + PACS Sonde # 37 don't reach 700 mb - Also correction for NaN sondes
       
                cr(i).sound.level.mb700.t(w)       = nan;
                cr(i).sound.level.mb700.theta(w)   = nan;
                cr(i).sound.level.mb700.thetav(w)  = nan;
                cr(i).sound.level.mb700.thetae(w)  = nan;
                cr(i).sound.level.mb700.thetaes(w) = nan;
                cr(i).sound.level.mb700.r(w)       = nan;
                cr(i).sound.level.mb700.rh(w)      = nan;
                cr(i).sound.level.mb700.wspd(w)    = nan;
                cr(i).sound.level.mb700.wdir(w)    = nan;     
                continue    
            elseif length(pr_700_dummy)==2     %%%% Correction if the difference of the values right above and below 700 mb are exactly the same 
                pr_700_ind = pr_700_dummy(1);
            else
                pr_700_ind = pr_700_dummy;
            end
            
            cr(i).sound.level.mb700.t(w)       = tk(pr_700_ind,w);
            cr(i).sound.level.mb700.theta(w)   = theta(pr_700_ind,w);
            cr(i).sound.level.mb700.thetav(w)  = thetav(pr_700_ind,w);
            cr(i).sound.level.mb700.thetae(w)  = thetae(pr_700_ind,w);
            cr(i).sound.level.mb700.thetaes(w) = thetaes(pr_700_ind,w);
            cr(i).sound.level.mb700.r(w)       = r(pr_700_ind,w);
            cr(i).sound.level.mb700.rh(w)      = rh(pr_700_ind,w);
            cr(i).sound.level.mb700.wspd(w)    = spd(pr_700_ind,w);
            cr(i).sound.level.mb700.wdir(w)    = dir(pr_700_ind,w);
        end
        clear w;
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% Estimation of Inversion Values %%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        for w = 1:length(mday)
        
        i_1005 = 10;

            for k = i_1005:200;
                beta(k,w)      = 0.608*theta(k,w)/(1+0.608*r(k,w)/1000);
                dtheta_dp(k,w) = (theta(k+1,w)-theta(k,w))/(p(k+1,w)-p(k,w));
                dr_dp(k,w)     = (10^(-3)*(r(k+1,w)-r(k,w)))/(p(k+1,w)-p(k,w));
                mu(k,w)        = -(dtheta_dp(k,w) - beta(k,w)*dr_dp(k,w));
            end

            if isnan(mu(41:100,w)) | (i==1 && w==14) %%%% Exclude all problematic soundings with NaN values (limits - 410 to 1000 m - work for all cruises)  (w==14 is only for EPIC)
    
                cr(i).sound.level.inv.top_h(w)    = nan;
                cr(i).sound.level.inv.top_th(w)   = nan;
                cr(i).sound.level.inv.top_r(w)    = nan;
                cr(i).sound.level.inv.base_h(w)   = nan;
                cr(i).sound.level.inv.base_th(w)  = nan;
                cr(i).sound.level.inv.base_r(w)   = nan;
                cr(i).sound.level.inv.thick(w)    = nan;
                cr(i).sound.level.inv.dtheta(w)   = nan;
                cr(i).sound.level.inv.dr(w)       = nan;
                cr(i).sound.level.inv.shear(w)    = nan;
                cr(i).sound.level.inv.strength(w) = nan;
            else
                peak  = find(mu(1:200,w)==max(mu(1:200,w)));  
                dummy = find(mu(1:200,w)<0.3);                  %%%% Maybe need to change 0.3
                inv_b = dummy(max(find(peak-dummy>0)));
                inv_t = dummy(min(find(peak-dummy<0)));
                cr(i).sound.level.inv.top_h(w)   = h(inv_t);
                cr(i).sound.level.inv.top_th(w)  = theta(inv_t,w);
                cr(i).sound.level.inv.top_r(w)   = r(inv_t,w);
                cr(i).sound.level.inv.base_h(w)  = h(inv_b);
                cr(i).sound.level.inv.base_th(w) = theta(inv_b,w);
                cr(i).sound.level.inv.base_r(w)  = r(inv_b,w);
  
                cr(i).sound.level.inv.thick(w)  = cr(i).sound.level.inv.top_h(w)-cr(i).sound.level.inv.base_h(w);
                cr(i).sound.level.inv.dtheta(w) = cr(i).sound.level.inv.top_th(w)-cr(i).sound.level.inv.base_th(w);
                cr(i).sound.level.inv.dr(w)     = cr(i).sound.level.inv.top_r(w)-cr(i).sound.level.inv.base_r(w);
                
                wind_invb  = spd(inv_b,w); 
                wind_invt  = spd(inv_t,w);
                cr(i).sound.level.inv.shear(w) = wind_invt-wind_invb;
                
                cr(i).sound.level.inv.strength(w) = nanmean(mu(inv_b:inv_t,w));  %%%% Inversion strength
                
                %inv_max = find(mu(:,w)==max(mu(:,w)));          
                %cr(i).sound.inv.strength2(w) = mu(inv_max,w);
                
                %inv_top_pr  = p(inv_t,w);
                %inv_base_pr = p(inv_b,w);
                %inv_dth_dp  = (cr(i).sound.inv_top_th(w)-cr(i).sound.inv_base_th(w))/(inv_top_pr-inv_base_pr);
                %inv_dr_dp   = (cr(i).sound.inv_top_r(w)-cr(i).sound.inv_base_r(w))/(inv_top_pr-inv_base_pr);
                %inv_beta    = 0.608*theta(inv_b,w)/(1+0.608*r(inv_b,w)/1000);
                %cr(i).sound.inv.strength3(w) = -(inv_dth_dp - inv_beta*inv_dr_dp);
                
                %cr(i).sound.inv.rich(w)  = (9.81*(thetav(inv_t,w)-thetav(inv_b,w))*(h(inv_t)-h(inv_b)))/(thetav(inv_b,w)*((wind_invt-wind_invb)^2));
                
                               
            end
            
        end
        
        clear th* pr* wi* h t p r* rh spd dir mday w k
        clear inv* peak dummy beta dtheta_dp dr_dp mu* i_1005
        
        if i==1     %%%% Define starting (ws) and ending (we) sondes for buoy period
             ws=45; we=92;
        elseif i==2
             ws=10; we=29;
        else ws=21; we=50;
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% Estimation of Mean and Std Values %%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        cr(i).sound.level.mb700.mean.t       = nanmean(cr(i).sound.level.mb700.t(ws:we));
        cr(i).sound.level.mb700.mean.theta   = nanmean(cr(i).sound.level.mb700.theta(ws:we));
        cr(i).sound.level.mb700.mean.thetav  = nanmean(cr(i).sound.level.mb700.thetav(ws:we));
        cr(i).sound.level.mb700.mean.thetae  = nanmean(cr(i).sound.level.mb700.thetae(ws:we));
        cr(i).sound.level.mb700.mean.thetaes = nanmean(cr(i).sound.level.mb700.thetaes(ws:we));
        cr(i).sound.level.mb700.mean.r       = nanmean(cr(i).sound.level.mb700.r(ws:we));
        cr(i).sound.level.mb700.mean.rh      = nanmean(cr(i).sound.level.mb700.rh(ws:we));
        cr(i).sound.level.mb700.mean.wspd    = nanmean(cr(i).sound.level.mb700.wspd(ws:we));
        cr(i).sound.level.mb700.mean.wdir    = nanmean(cr(i).sound.level.mb700.wdir(ws:we));
        
        cr(i).sound.level.mb700.std.t       = nanstd(cr(i).sound.level.mb700.t(ws:we));
        cr(i).sound.level.mb700.std.theta   = nanstd(cr(i).sound.level.mb700.theta(ws:we));
        cr(i).sound.level.mb700.std.thetav  = nanstd(cr(i).sound.level.mb700.thetav(ws:we));
        cr(i).sound.level.mb700.std.thetae  = nanstd(cr(i).sound.level.mb700.thetae(ws:we));
        cr(i).sound.level.mb700.std.thetaes = nanstd(cr(i).sound.level.mb700.thetaes(ws:we));
        cr(i).sound.level.mb700.std.r       = nanstd(cr(i).sound.level.mb700.r(ws:we));
        cr(i).sound.level.mb700.std.rh      = nanstd(cr(i).sound.level.mb700.rh(ws:we));
        cr(i).sound.level.mb700.std.wspd    = nanstd(cr(i).sound.level.mb700.wspd(ws:we));
        cr(i).sound.level.mb700.std.wdir    = nanstd(cr(i).sound.level.mb700.wdir(ws:we));
        
        cr(i).sound.level.mb1000.mean.t       = nanmean(cr(i).sound.level.mb1000.t(ws:we));
        cr(i).sound.level.mb1000.mean.theta   = nanmean(cr(i).sound.level.mb1000.theta(ws:we));
        cr(i).sound.level.mb1000.mean.thetav  = nanmean(cr(i).sound.level.mb1000.thetav(ws:we));
        cr(i).sound.level.mb1000.mean.thetae  = nanmean(cr(i).sound.level.mb1000.thetae(ws:we));
        cr(i).sound.level.mb1000.mean.thetaes = nanmean(cr(i).sound.level.mb1000.thetaes(ws:we));
        cr(i).sound.level.mb1000.mean.r       = nanmean(cr(i).sound.level.mb1000.r(ws:we));
        cr(i).sound.level.mb1000.mean.rh      = nanmean(cr(i).sound.level.mb1000.rh(ws:we));
        cr(i).sound.level.mb1000.mean.wspd    = nanmean(cr(i).sound.level.mb1000.wspd(ws:we));
        cr(i).sound.level.mb1000.mean.wdir    = nanmean(cr(i).sound.level.mb1000.wdir(ws:we));
        
        cr(i).sound.level.mb1000.std.t       = nanstd(cr(i).sound.level.mb1000.t(ws:we));
        cr(i).sound.level.mb1000.std.theta   = nanstd(cr(i).sound.level.mb1000.theta(ws:we));
        cr(i).sound.level.mb1000.std.thetav  = nanstd(cr(i).sound.level.mb1000.thetav(ws:we));
        cr(i).sound.level.mb1000.std.thetae  = nanstd(cr(i).sound.level.mb1000.thetae(ws:we));
        cr(i).sound.level.mb1000.std.thetaes = nanstd(cr(i).sound.level.mb1000.thetaes(ws:we));
        cr(i).sound.level.mb1000.std.r       = nanstd(cr(i).sound.level.mb1000.r(ws:we));
        cr(i).sound.level.mb1000.std.rh      = nanstd(cr(i).sound.level.mb1000.rh(ws:we));
        cr(i).sound.level.mb1000.std.wspd    = nanstd(cr(i).sound.level.mb1000.wspd(ws:we));
        cr(i).sound.level.mb1000.std.wdir    = nanstd(cr(i).sound.level.mb1000.wdir(ws:we));
        
        cr(i).sound.level.inv.mean.top_h    = nanmean(cr(i).sound.level.inv.top_h(ws:we));
        cr(i).sound.level.inv.mean.top_th   = nanmean(cr(i).sound.level.inv.top_th(ws:we));
        cr(i).sound.level.inv.mean.top_r    = nanmean(cr(i).sound.level.inv.top_r(ws:we));
        cr(i).sound.level.inv.mean.base_h   = nanmean(cr(i).sound.level.inv.base_h(ws:we));
        cr(i).sound.level.inv.mean.base_th  = nanmean(cr(i).sound.level.inv.base_th(ws:we));
        cr(i).sound.level.inv.mean.base_r   = nanmean(cr(i).sound.level.inv.base_r(ws:we));
        cr(i).sound.level.inv.mean.thick    = nanmean(cr(i).sound.level.inv.thick(ws:we));
        cr(i).sound.level.inv.mean.dtheta   = nanmean(cr(i).sound.level.inv.dtheta(ws:we));
        cr(i).sound.level.inv.mean.dr       = nanmean(cr(i).sound.level.inv.dr(ws:we));
        cr(i).sound.level.inv.mean.strength = nanmean(cr(i).sound.level.inv.strength(ws:we));
        cr(i).sound.level.inv.mean.shear    = nanmean(cr(i).sound.level.inv.shear(ws:we));
        
        cr(i).sound.level.inv.std.top_h    = nanstd(cr(i).sound.level.inv.top_h(ws:we));
        cr(i).sound.level.inv.std.top_th   = nanstd(cr(i).sound.level.inv.top_th(ws:we));
        cr(i).sound.level.inv.std.top_r    = nanstd(cr(i).sound.level.inv.top_r(ws:we));
        cr(i).sound.level.inv.std.base_h   = nanstd(cr(i).sound.level.inv.base_h(ws:we));
        cr(i).sound.level.inv.std.base_th  = nanstd(cr(i).sound.level.inv.base_th(ws:we));
        cr(i).sound.level.inv.std.base_r   = nanstd(cr(i).sound.level.inv.base_r(ws:we));
        cr(i).sound.level.inv.std.thick    = nanstd(cr(i).sound.level.inv.thick(ws:we));
        cr(i).sound.level.inv.std.dtheta   = nanstd(cr(i).sound.level.inv.dtheta(ws:we));
        cr(i).sound.level.inv.std.dr       = nanstd(cr(i).sound.level.inv.dr(ws:we));
        cr(i).sound.level.inv.std.strength = nanstd(cr(i).sound.level.inv.strength(ws:we));
        cr(i).sound.level.inv.std.shear    = nanstd(cr(i).sound.level.inv.shear(ws:we));
              
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Ceilometer & Radar %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%clear all
keep cr

ceilo_path_sec  = {'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'};
ceilo_path_min  = {'C:\THESIS\DATA\EPIC\data\ceilo\ETL_web\epic2_ceilo_10min.txt';'C:\THESIS\DATA\PACS03\DATA\ceilo\STRATUS03_ceilo_10min_a.txt';'C:\THESIS\DATA\Stratus04\ceilo\STRATUS04_ceilo_10min_a.txt'};
ceilo_path_hr = {'C:\THESIS\DATA\EPIC\data\ceilo\ETL_web\epic2_ceilo_h.txt';'C:\THESIS\DATA\PACS03\DATA\ceilo\STRATUS03_ceilo_h_a.txt';'C:\THESIS\DATA\Stratus04\ceilo\STRATUS04_ceilo_h_a.txt'};

radar_path_hr = {'C:\THESIS\Processing_Results\EPIC2001\MMCR\EPIC_boundaries_hr';'C:\THESIS\Processing_Results\Stratus03\MMCR\PACS_boundaries_hr';'C:\THESIS\Processing_Results\Stratus04\windprof\stratus04_boundaries_hr'};
    
jdc = [273 304 335];    %%%% correction coefficent to get day of month from Julian day for each cruise



%i=1        
for i=1:3    %%%% It can be reduced to one loop without the "ifs" using cr(i).ind etc.
    
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% Mean and Std values for Cloud Boundaries, Fraction & Thickness %%%%
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        cr(i).ceilo.sec.data = load(char(ceilo_path_sec(i)));
        cr(i).ceilo.min.data = load(char(ceilo_path_min(i)));
        cr(i).ceilo.hr.data  = load(char(ceilo_path_hr(i)));
        
        cr(i).ceilo.sec.mday = cr(i).ceilo.sec.data(:,1)-jdc(i);
        cr(i).ceilo.min.mday = cr(i).ceilo.min.data(:,1)-jdc(i);
        cr(i).ceilo.hr.mday = cr(i).ceilo.hr.data(:,1)-jdc(i);
        
        cr(i).ceilo.sec.clba = cr(i).ceilo.sec.data(:,3);
        cr(i).ceilo.min.clba = cr(i).ceilo.min.data(:,11);
        cr(i).ceilo.hr.clba = cr(i).ceilo.hr.data(:,11);
        
        cr(i).ceilo.sec.clla = cr(i).ceilo.sec.data(:,2);  
        cr(i).ceilo.min.clfr = cr(i).ceilo.min.data(:,10);   %%%% We include obscure returns (9 for not including obscure)
        %cr(i).ceilo.hr.clfr  = cr(i).ceilo.hr.data(:,10);    %%%% We include obscure returns (9 for not including obscure)
        
        cr(i).radar.hr = load(char(radar_path_hr(i)));
                
        if i==1     %%%% Define starting (ds) and ending (de) time for buoy period
            ds1=16; de1=22;
            ind1_sec_cf = find(cr(i).ceilo.sec.mday>=ds1 & cr(i).ceilo.sec.mday<de1);
            ind1_sec_cb = find(cr(i).ceilo.sec.mday>=ds1 & cr(i).ceilo.sec.mday<de1 & cr(i).ceilo.sec.clba>50);
            ind1_min = find(cr(i).ceilo.min.mday>=ds1 & cr(i).ceilo.min.mday<de1);
            ind1_hr  = find(cr(i).ceilo.hr.mday>=ds1 & cr(i).ceilo.hr.mday<de1);
   
            alldata1 = length(ind1_sec_cf); 
            iscloud1 = find(cr(i).ceilo.sec.clla(ind1_sec_cf)>0 & cr(i).ceilo.sec.clla(ind1_sec_cf)<5);
            cloud1 = length(iscloud1);
            cr(i).ceilo.sec.mean.clfr = cloud1/alldata1;
            cr(i).ceilo.sec.mean.clba = nanmean(cr(i).ceilo.sec.clba(ind1_sec_cb));
            cr(i).ceilo.sec.std.clba  = nanstd(cr(i).ceilo.sec.clba(ind1_sec_cb));
            
            cr(i).ceilo.min.mean.clfr = nanmean(cr(i).ceilo.min.clfr(ind1_min));
            cr(i).ceilo.min.mean.clba = nanmean(cr(i).ceilo.min.clba(ind1_min));
            cr(i).ceilo.min.std.clfr  = nanstd(cr(i).ceilo.min.clfr(ind1_min));
            cr(i).ceilo.min.std.clba  = nanstd(cr(i).ceilo.min.clba(ind1_min));
            
            %cr(i).ceilo.hr.mean.clfr = nanmean(cr(i).ceilo.hr.clfr(ind1_hr));
            cr(i).ceilo.hr.mean.clba = nanmean(cr(i).ceilo.hr.clba(ind1_hr));
            %cr(i).ceilo.hr.std.clfr  = nanstd(cr(i).ceilo.hr.clfr(ind1_hr));
            cr(i).ceilo.hr.std.clba  = nanstd(cr(i).ceilo.hr.clba(ind1_hr));
            
            
            ind1_ra  = find(cr(i).radar.hr.mday>=ds1 & cr(i).radar.hr.mday<de1);
            
            cr(i).radar.hr.mean.cltop = nanmean(cr(i).radar.hr.cltop(ind1_ra));
            cr(i).radar.hr.mean.clthi = nanmean(cr(i).radar.hr.clthi(ind1_ra));
            cr(i).radar.hr.std.cltop  = nanstd(cr(i).radar.hr.cltop(ind1_ra));
            cr(i).radar.hr.std.clthi  = nanstd(cr(i).radar.hr.clthi(ind1_ra));
                        
        elseif i==2
            ds2=16; de2=21;
            ind2_sec_cf = find(cr(i).ceilo.sec.mday>=ds2 & cr(i).ceilo.sec.mday<de2);
            ind2_sec_cb = find(cr(i).ceilo.sec.mday>=ds2 & cr(i).ceilo.sec.mday<de2 & cr(i).ceilo.sec.clba>50 & cr(i).ceilo.sec.clba<1600);
            ind2_min_cf = find(cr(i).ceilo.min.mday>=ds2 & cr(i).ceilo.min.mday<de2);
            ind2_min_cb = find(cr(i).ceilo.min.mday>=ds2 & cr(i).ceilo.min.mday<de2 & cr(i).ceilo.min.clba<1600);
            ind2_hr_cf  = find(cr(i).ceilo.hr.mday>=ds2 & cr(i).ceilo.hr.mday<de2);
            ind2_hr_cb  = find(cr(i).ceilo.hr.mday>=ds2 & cr(i).ceilo.hr.mday<de2 & cr(i).ceilo.hr.clba<1600);
            
            alldata2 = length(ind2_sec_cf); 
            iscloud2 = find(cr(i).ceilo.sec.clla(ind2_sec_cf)>0 & cr(i).ceilo.sec.clla(ind2_sec_cf)<5);
            cloud2 = length(iscloud2);
            cr(i).ceilo.sec.mean.clfr = cloud2/alldata2;
            cr(i).ceilo.sec.mean.clba = nanmean(cr(i).ceilo.sec.clba(ind2_sec_cb));
            cr(i).ceilo.sec.std.clba  = nanstd(cr(i).ceilo.sec.clba(ind2_sec_cb));
            
            cr(i).ceilo.min.mean.clfr = nanmean(cr(i).ceilo.min.clfr(ind2_min_cf));
            cr(i).ceilo.min.mean.clba = nanmean(cr(i).ceilo.min.clba(ind2_min_cb));
            cr(i).ceilo.min.std.clfr  = nanstd(cr(i).ceilo.min.clfr(ind2_min_cf));
            cr(i).ceilo.min.std.clba  = nanstd(cr(i).ceilo.min.clba(ind2_min_cb));
            
            %cr(i).ceilo.hr.mean.clfr = nanmean(cr(i).ceilo.hr.clfr(ind2_hr_cf));
            cr(i).ceilo.hr.mean.clba = nanmean(cr(i).ceilo.hr.clba(ind2_hr_cb));
            %cr(i).ceilo.hr.std.clfr  = nanstd(cr(i).ceilo.hr.clfr(ind2_hr_cf));
            cr(i).ceilo.hr.std.clba  = nanstd(cr(i).ceilo.hr.clba(ind2_hr_cb));
            
            
            ind2_ra  = find(cr(i).radar.hr.mday>=ds2 & cr(i).radar.hr.mday<de2);
            
            cr(i).radar.hr.mean.cltop = nanmean(cr(i).radar.hr.cltop(ind2_ra));
            cr(i).radar.hr.mean.clthi = nanmean(cr(i).radar.hr.clthi(ind2_ra));
            cr(i).radar.hr.std.cltop  = nanstd(cr(i).radar.hr.cltop(ind2_ra));
            cr(i).radar.hr.std.clthi  = nanstd(cr(i).radar.hr.clthi(ind2_ra));
            
        else
            ds3=11.2; de3=16.2;
            ind3_sec_cf = find(cr(i).ceilo.sec.mday>=ds3 & cr(i).ceilo.sec.mday<de3);
            ind3_sec_cb = find(cr(i).ceilo.sec.mday>=ds3 & cr(i).ceilo.sec.mday<de3 & cr(i).ceilo.sec.clba>200 & cr(i).ceilo.sec.clba<1600);
            ind3_min_cf = find(cr(i).ceilo.min.mday>=ds3 & cr(i).ceilo.min.mday<de3);
            ind3_min_cb = find(cr(i).ceilo.min.mday>=ds3 & cr(i).ceilo.min.mday<de3 & cr(i).ceilo.min.clba>200 & cr(i).ceilo.min.clba<1600);
            ind3_hr_cf  = find(cr(i).ceilo.hr.mday>=ds3 & cr(i).ceilo.hr.mday<de3);
            ind3_hr_cb  = find(cr(i).ceilo.hr.mday>=ds3 & cr(i).ceilo.hr.mday<de3 & cr(i).ceilo.hr.clba>200 & cr(i).ceilo.hr.clba<1600);
            
            alldata3 = length(ind3_sec_cf); 
            iscloud3 = find(cr(i).ceilo.sec.clla(ind3_sec_cf)>0 & cr(i).ceilo.sec.clla(ind3_sec_cf)<5);
            cloud3 = length(iscloud3);
            cr(i).ceilo.sec.mean.clfr = cloud3/alldata3;
            cr(i).ceilo.sec.mean.clba = nanmean(cr(i).ceilo.sec.clba(ind3_sec_cb));
            cr(i).ceilo.sec.std.clba  = nanstd(cr(i).ceilo.sec.clba(ind3_sec_cb));
            
            cr(i).ceilo.min.mean.clfr = nanmean(cr(i).ceilo.min.clfr(ind3_min_cf));
            cr(i).ceilo.min.mean.clba = nanmean(cr(i).ceilo.min.clba(ind3_min_cb));
            cr(i).ceilo.min.std.clfr  = nanstd(cr(i).ceilo.min.clfr(ind3_min_cf));
            cr(i).ceilo.min.std.clba  = nanstd(cr(i).ceilo.min.clba(ind3_min_cb));
            
            %cr(i).ceilo.hr.mean.clfr = nanmean(cr(i).ceilo.hr.clfr(ind3_hr_cf));
            cr(i).ceilo.hr.mean.clba = nanmean(cr(i).ceilo.hr.clba(ind3_hr_cb));
            %cr(i).ceilo.hr.std.clfr  = nanstd(cr(i).ceilo.hr.clfr(ind3_hr_cf));
            cr(i).ceilo.hr.std.clba  = nanstd(cr(i).ceilo.hr.clba(ind3_hr_cb));
            
            
            ind3_ra  = find(cr(i).radar.hr.mday>=ds3 & cr(i).radar.hr.mday<de3);
            
            cr(i).radar.hr.mean.cltop = nanmean(cr(i).radar.hr.cltop(ind3_ra));
            cr(i).radar.hr.mean.clthi = nanmean(cr(i).radar.hr.clthi(ind3_ra));
            cr(i).radar.hr.std.cltop  = nanstd(cr(i).radar.hr.cltop(ind3_ra));
            cr(i).radar.hr.std.clthi  = nanstd(cr(i).radar.hr.clthi(ind3_ra));
        end       
end

%%%%%%%%%%%%%%%%%%%
%%%% Flux data %%%%
%%%%%%%%%%%%%%%%%%%

%clear all;
keep cr
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];

%jdc2 = [9 10 4];
%i=3;

for i=1:3    %%%% 1 corresponds to EPIC, 2 to Stratus03 and 3 to Stratus04 (all are 5-min resolution)
    
    if i==1     %%%% Define starting (ds) and ending (de) time for buoy period
        ds=16; de=22;
    elseif i==2
        ds=16; de=21;
    else
        ds=11.2; de=16.2;
    end
     
    cr(i).flux.min.data = load(char(flux_path(i)));
    cr(i).flux.min.data(:,25:end) = [];
   
    cr(i).flux.min.mday = cr(i).flux.min.data(:,1)-jdc(i);
    
    cr(i).flux.min.wspd = cr(i).flux.min.data(:,2);
    cr(i).flux.min.wdir = cr(i).flux.min.data(:,3);
    cr(i).flux.min.sst  = cr(i).flux.min.data(:,4);
    cr(i).flux.min.tair = cr(i).flux.min.data(:,7);
    cr(i).flux.min.tdif = cr(i).flux.min.sst-cr(i).flux.min.tair;
    cr(i).flux.min.qsea = cr(i).flux.min.data(:,8);
    cr(i).flux.min.qair = cr(i).flux.min.data(:,9);
    cr(i).flux.min.sw   = cr(i).flux.min.data(:,10);
    cr(i).flux.min.lw   = cr(i).flux.min.data(:,11);
    cr(i).flux.min.sh   = cr(i).flux.min.data(:,22);
    cr(i).flux.min.lh   = cr(i).flux.min.data(:,23);
    cr(i).flux.min.vh   = cr(i).flux.min.sh+0.608*292*0.0004*cr(i).flux.min.lh;
    
    %%%% Addition of LCL data %%%%
    
    z0 = 15.5;
    p_lcl = interp1(cr(i).sound.data.mday, cr(i).sound.data.p(2,:), cr(i).flux.min.mday);            %%%% Interpolation of surface pressure data from soundings to the flux data
    r_lcl = ((cr(i).flux.min.qair/1000)./(1-(cr(i).flux.min.qair/1000)))*1000;                       %%%% Surface mixing ratio in g/kg
    e_lcl = (p_lcl.*r_lcl)./(622+r_lcl);                                                             %%%% Surface vapor pressure in mb
    cr(i).flux.min.tlcl = (2840./(3.5*log(cr(i).flux.min.tair+273.15)-log(e_lcl)-4.805))+55-273.15;  %%%% Temperature of LCL in Celsius
    cr(i).flux.min.zlcl = z0 + (cr(i).flux.min.tair-cr(i).flux.min.tlcl)/0.0098;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    cr(i).ind = find(cr(i).flux.min.mday>=ds & cr(i).flux.min.mday<de);
    
    cr(i).flux.min.mean.wspd = nanmean(cr(i).flux.min.wspd(cr(i).ind));
    cr(i).flux.min.mean.wdir = nanmean(cr(i).flux.min.wdir(cr(i).ind));
    cr(i).flux.min.mean.sst  = nanmean(cr(i).flux.min.sst(cr(i).ind));
    cr(i).flux.min.mean.tair = nanmean(cr(i).flux.min.tair(cr(i).ind));
    cr(i).flux.min.mean.tdif = nanmean(cr(i).flux.min.tdif(cr(i).ind));
    cr(i).flux.min.mean.qsea = nanmean(cr(i).flux.min.qsea(cr(i).ind));
    cr(i).flux.min.mean.qair = nanmean(cr(i).flux.min.qair(cr(i).ind));
    cr(i).flux.min.mean.sw   = nanmean(cr(i).flux.min.sw(cr(i).ind));
    cr(i).flux.min.mean.lw   = nanmean(cr(i).flux.min.lw(cr(i).ind));
    cr(i).flux.min.mean.sh   = nanmean(cr(i).flux.min.sh(cr(i).ind));
    cr(i).flux.min.mean.lh   = nanmean(cr(i).flux.min.lh(cr(i).ind));
    cr(i).flux.min.mean.vh   = nanmean(cr(i).flux.min.vh(cr(i).ind));
    
    cr(i).flux.min.mean.tlcl = nanmean(cr(i).flux.min.tlcl(cr(i).ind));
    cr(i).flux.min.mean.zlcl = nanmean(cr(i).flux.min.zlcl(cr(i).ind));
    
    cr(i).flux.min.std.wspd = nanstd(cr(i).flux.min.wspd(cr(i).ind));
    cr(i).flux.min.std.wdir = nanstd(cr(i).flux.min.wdir(cr(i).ind));
    cr(i).flux.min.std.sst  = nanstd(cr(i).flux.min.sst(cr(i).ind));
    cr(i).flux.min.std.tair = nanstd(cr(i).flux.min.tair(cr(i).ind));
    cr(i).flux.min.std.tdif = nanstd(cr(i).flux.min.tdif(cr(i).ind));
    cr(i).flux.min.std.qsea = nanstd(cr(i).flux.min.qsea(cr(i).ind));
    cr(i).flux.min.std.qair = nanstd(cr(i).flux.min.qair(cr(i).ind));
    cr(i).flux.min.std.sw   = nanstd(cr(i).flux.min.sw(cr(i).ind));
    cr(i).flux.min.std.lw   = nanstd(cr(i).flux.min.lw(cr(i).ind));
    cr(i).flux.min.std.sh   = nanstd(cr(i).flux.min.sh(cr(i).ind));
    cr(i).flux.min.std.lh   = nanstd(cr(i).flux.min.lh(cr(i).ind));
    cr(i).flux.min.std.vh   = nanstd(cr(i).flux.min.vh(cr(i).ind));
    
    cr(i).flux.min.std.tlcl = nanstd(cr(i).flux.min.tlcl(cr(i).ind));
    cr(i).flux.min.std.zlcl = nanstd(cr(i).flux.min.zlcl(cr(i).ind));
        
    clear p_lcl r_lcl e_lcl
end


keep cr
cr = rmfield(cr, 'ind');



