%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        %
%    This part of the program plot the   %
%    SR Histo as a function of altitude  %
%                                        %
%    Just modify the list with yours     %
%                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all;
close all;

%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/jja_list');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end



%%% initialization of the diag vector
mdiag(1:15,1:length(alt_mid),1:length(list))=0;

for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ;   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    diagSR=squeeze(nansum(nanval,4));  % average on the longitude column
    diagSR=squeeze(nansum(diagSR,3));  % average on the latitude column

    mdiag(:,:,i)=diagSR;             
 
clear nanval cfad_lidarsr532_Occ diagSR;


end

diagSR=squeeze(nansum(mdiag,3));
diagtot=nansum(diagSR(:));         % total of all SR occurence


z(1:length(alt),1:16)=0;
z(1:length(alt_mid),1:15)=diagSR(:,:)'; 

t=z/diagtot;      % conversion of occurence in fraction
t(t<0.0001)=nan;  % set the nan values

% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];


% plot of diagSR average

figure(3)
pcolor(1:16,alt,t)
hold on
plot([-90 90],[3.085 3.085],'r-')
hold on
plot([-90 90],[7.484 7.484],'r-')
set(gca,'xtick',[1:16])
set(gca,'xticklabel',{'0','0.01','1.2','3','5','7','10','15','20','25','30','40','50','60','80','100000'})
shading flat;
colorbar
colormap(a)
xlabel('SR')
ylabel('ALTITUDE (KM)')
axis([1 16 0 19.2])
caxis([0   0.001]);
title('SR Histo JJA GOCCP ')

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







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        %
%    This part of the program plot the   %
%    SR Histo as a function of different %
%    regions                             %
%                                        %
%    Just modify the list with yours     %
%                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



clear all;
%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/list_SRhist_JFM_CFMIP2_2.1');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end


%%% coordinates of the region
x=find(longitude>=70 & longitude<=150); % longitude
y=find(latitude>=-5 & latitude<=20);  % latitude

% fid=fopen('/home/gcesana/src/Figure_Helene/fig_papier/table/liste_jfm_night_cfmip_1.85.diag');

%%% initialization of the diag vector
mdiag(1:15,1:length(alt_mid),1:length(list))=0;

for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ(:,:,y,x);   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    diagSR=squeeze(nansum(nanval,4));  % average on the longitude column
    diagSR=squeeze(nansum(diagSR,3));  % average on the latitude column

    mdiag(:,:,i)=diagSR;             
 
clear nanval cfad_lidarsr532_Occ diagSR;


end

diagSR=squeeze(nansum(mdiag,3));
diagtot=nansum(diagSR(:));         % total of all SR occurence


z(1:length(alt),1:16)=0;
z(1:length(alt_mid),1:15)=diagSR(:,:)'; 

t=z/diagtot;      % conversion of occurence in fraction
t(t<0.0001)=nan;  % set the nan values

% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];


% plot of diagSR average

figure(1)

subplot(4,2,1)
pcolor(1:16,alt,t)
hold on
plot([-90 90],[3.085 3.085],'r-')
hold on
plot([-90 90],[7.484 7.484],'r-')
set(gca,'xtick',[1:16])
set(gca,'xticklabel',{'0','0.01','1.2','3','5','7','10','15','20','25','30','40','50','60','80','100000'})
shading flat;
xlabel('SR')
ylabel('ALTITUDE (KM)')
axis([1 16 0 19.2])
caxis([0   0.001]);
title('SR Histo Tropical Western Pacific JFM GOCCP ')




clear all;
%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/list_SRhist_JFM_CFMIP2_2.1');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end


x=find(longitude>=-180 & longitude<=-140);
y=find(latitude>=15 & latitude<=35);

% initialization of the diag vector
mdiag(1:15,1:length(alt_mid),1:length(list))=0;

for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ(:,:,y,x);   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    diagSR=squeeze(nansum(nanval,4));  % average on the longitude column
    diagSR=squeeze(nansum(diagSR,3));  % average on the latitude column

    mdiag(:,:,i)=diagSR;             
 
clear nanval cfad_lidarsr532_Occ diagSR;


end

diagSR=squeeze(nansum(mdiag,3));
diagtot=nansum(diagSR(:));         % total of all SR occurence


z(1:length(alt),1:16)=0;
z(1:length(alt_mid),1:15)=diagSR(:,:)'; 

t=z/diagtot;      % conversion of occurence in fraction
t(t<0.0001)=nan;  % set the nan values

% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];


% plot of diagSR average

figure(1)

subplot(4,2,3)
pcolor(1:16,alt,t)
hold on
plot([-90 90],[3.085 3.085],'r-')
hold on
plot([-90 90],[7.484 7.484],'r-')
set(gca,'xtick',[1:16])
set(gca,'xticklabel',{'0','0.01','1.2','3','5','7','10','15','20','25','30','40','50','60','80','100000'})
shading flat;
xlabel('SR')
ylabel('ALTITUDE (KM)')

axis([1 16 0 19.2])
caxis([0   0.001]);
title('SR Histo Hawai Trade Cumulus JFM GOCCP ')





clear all;

%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/list_SRhist_JFM_CFMIP2_2.1');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end



x=find(longitude>=-140 & longitude<=-110);
y=find(latitude>=15 & latitude<=35);

% initialization of the diag vector
mdiag(1:15,1:length(alt_mid),1:length(list))=0;

for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ(:,:,y,x);   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    diagSR=squeeze(nansum(nanval,4));  % average on the longitude column
    diagSR=squeeze(nansum(diagSR,3));  % average on the latitude column

    mdiag(:,:,i)=diagSR;             
 
clear nanval cfad_lidarsr532_Occ diagSR;


end

diagSR=squeeze(nansum(mdiag,3));
diagtot=nansum(diagSR(:));         % total of all SR occurence


z(1:length(alt),1:16)=0;
z(1:length(alt_mid),1:15)=diagSR(:,:)'; 

t=z/diagtot;      % conversion of occurence in fraction
t(t<0.0001)=nan;  % set the nan values

% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];


% plot of diagSR average

figure(1)

subplot(4,2,5)
pcolor(1:16,alt,t)
hold on
plot([-90 90],[3.085 3.085],'r-')
hold on
plot([-90 90],[7.484 7.484],'r-')
set(gca,'xtick',[1:16])
set(gca,'xticklabel',{'0','0.01','1.2','3','5','7','10','15','20','25','30','40','50','60','80','100000'})
shading flat;
xlabel('SR')
ylabel('ALTITUDE (KM)')
axis([1 16 0 19.2])
caxis([0   0.001]);
title('SR Histo California StratoCumulus JFM GOCCP ')



clear all;
%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/list_SRhist_JFM_CFMIP2_2.1');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end


xx=find(longitude>=-180 & longitude<=-125);
xxx=find(longitude>=160);
toto=length(xx)+length(xxx);
x(1:length(xx))=xx;
x(length(xx)+1:toto)=xxx;

y=find(latitude>=40 & latitude<=60);


% initialization of the diag vector
mdiag(1:15,1:length(alt_mid),1:length(list))=0;

for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ(:,:,y,x);   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    diagSR=squeeze(nansum(nanval,4));  % average on the longitude column
    diagSR=squeeze(nansum(diagSR,3));  % average on the latitude column

    mdiag(:,:,i)=diagSR;             
 
clear nanval cfad_lidarsr532_Occ diagSR;


end

diagSR=squeeze(nansum(mdiag,3));
diagtot=nansum(diagSR(:));         % total of all SR occurence


z(1:length(alt),1:16)=0;
z(1:length(alt_mid),1:15)=diagSR(:,:)'; 

t=z/diagtot;      % conversion of occurence in fraction
t(t<0.0001)=nan;  % set the nan values

% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];


% plot of diagSR average

figure(1)

subplot(4,2,7)
pcolor(1:16,alt,t)
hold on
plot([-90 90],[3.085 3.085],'r-')
hold on
plot([-90 90],[7.484 7.484],'r-')
set(gca,'xtick',[1:16])
set(gca,'xticklabel',{'0','0.01','1.2','3','5','7','10','15','20','25','30','40','50','60','80','100000'})
shading flat;
xlabel('SR')
xlabel('SR')
ylabel('ALTITUDE (KM)')
axis([1 16 0 19.2])
caxis([0   0.001]);
colormap(a)
title('SR Histo Mid Latitude North Pacific JFM GOCCP ')




clear all;
%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/list_SRhist_JJA_CFMIP2_2.1');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end


x=find(longitude>=70 & longitude<=150);
y=find(latitude>=-5 & latitude<=20);


% initialization of the diag vector
mdiag(1:15,1:length(alt_mid),1:length(list))=0;

for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ(:,:,y,x);   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    diagSR=squeeze(nansum(nanval,4));  % average on the longitude column
    diagSR=squeeze(nansum(diagSR,3));  % average on the latitude column

    mdiag(:,:,i)=diagSR;            
 
clear nanval cfad_lidarsr532_Occ diagSR;

end

diagSR=squeeze(nansum(mdiag,3));
diagtot=nansum(diagSR(:));         % total of all SR occurence


z(1:length(alt),1:16)=0;
z(1:length(alt_mid),1:15)=diagSR(:,:)'; 

t=z/diagtot;      % conversion of occurence in fraction
t(t<0.0001)=nan;  % set the nan values

% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];


% plot of diagSR average
figure(1)

subplot(4,2,2)
pcolor(1:16,alt,t)
hold on
plot([-90 90],[3.085 3.085],'r-')
hold on
plot([-90 90],[7.484 7.484],'r-')
set(gca,'xtick',[1:16])
set(gca,'xticklabel',{'0','0.01','1.2','3','5','7','10','15','20','25','30','40','50','60','80','100000'})
shading flat;
xlabel('SR')
ylabel('ALTITUDE (KM)')
axis([1 16 0 19.2])
caxis([0   0.001]);
title('SR Histo Tropical Western Pacific JJA GOCCP ')


clear all;
%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/list_SRhist_JJA_CFMIP2_2.1');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end


x=find(longitude>=-180 & longitude<=-140);
y=find(latitude>=15 & latitude<=35);


% initialization of the diag vector
mdiag(1:15,1:length(alt_mid),1:length(list))=0;

for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ(:,:,y,x);   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    diagSR=squeeze(nansum(nanval,4));  % average on the longitude column
    diagSR=squeeze(nansum(diagSR,3));  % average on the latitude column

    mdiag(:,:,i)=diagSR;             

end

diagSR=squeeze(nansum(mdiag,3));
diagtot=nansum(diagSR(:));         % total of all SR occurence


z(1:length(alt),1:16)=0;
z(1:length(alt_mid),1:15)=diagSR(:,:)'; 

t=z/diagtot;      % conversion of occurence in fraction
t(t<0.0001)=nan;  % set the nan values

% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];


% plot of diagSR average
figure(1)

subplot(4,2,4)
pcolor(1:16,alt,t)
hold on
plot([-90 90],[3.085 3.085],'r-')
hold on
plot([-90 90],[7.484 7.484],'r-')
set(gca,'xtick',[1:16])
set(gca,'xticklabel',{'0','0.01','1.2','3','5','7','10','15','20','25','30','40','50','60','80','100000'})
shading flat;
xlabel('SR')
ylabel('ALTITUDE (KM)')
axis([1 16 0 19.2])
caxis([0   0.001]);
title('SR Histo Hawai Trade Cumulus JJA GOCCP ')


clear all;
%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/list_SRhist_JJA_CFMIP2_2.1');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end


x=find(longitude>=-140 & longitude<=-110);
y=find(latitude>=15 & latitude<=35);

% initialization of the diag vector
mdiag(1:15,1:length(alt_mid),1:length(list))=0;

for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ(:,:,y,x);   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    diagSR=squeeze(nansum(nanval,4));  % average on the longitude column
    diagSR=squeeze(nansum(diagSR,3));  % average on the latitude column

    mdiag(:,:,i)=diagSR;             

end

diagSR=squeeze(nansum(mdiag,3));
diagtot=nansum(diagSR(:));         % total of all SR occurence


z(1:length(alt),1:16)=0;
z(1:length(alt_mid),1:15)=diagSR(:,:)'; 

t=z/diagtot;      % conversion of occurence in fraction
t(t<0.0001)=nan;  % set the nan values

% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];


% plot of diagSR average

figure(1)

subplot(4,2,6)
pcolor(1:16,alt,t)
hold on
plot([-90 90],[3.085 3.085],'r-')
hold on
plot([-90 90],[7.484 7.484],'r-')
set(gca,'xtick',[1:16])
set(gca,'xticklabel',{'0','0.01','1.2','3','5','7','10','15','20','25','30','40','50','60','80','100000'})
shading flat;
xlabel('SR')
ylabel('ALTITUDE (KM)')
axis([1 16 0 19.2])
caxis([0   0.001]);
title('SR Histo California StratoCumulus JJA GOCCP ')



clear all;

%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/list_SRhist_JJA_CFMIP2_2.1');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end


xx=find(longitude>=-180 & longitude<=-125);
xxx=find(longitude>=160);
toto=length(xx)+length(xxx);
x(1:length(xx))=xx;
x(length(xx)+1:toto)=xxx;

y=find(latitude>=40 & latitude<=60);


% initialization of the diag vector
mdiag(1:15,1:length(alt_mid),1:length(list))=0;

for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ(:,:,y,x);   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    diagSR=squeeze(nansum(nanval,4));  % average on the longitude column
    diagSR=squeeze(nansum(diagSR,3));  % average on the latitude column

    mdiag(:,:,i)=diagSR;             

end

diagSR=squeeze(nansum(mdiag,3));
diagtot=nansum(diagSR(:));         % total of all SR occurence


z(1:length(alt),1:16)=0;
z(1:length(alt_mid),1:15)=diagSR(:,:)'; 

t=z/diagtot;      % conversion of occurence in fraction
t(t<0.0001)=nan;  % set the nan values

% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];


% plot of diagSR average

figure(1)

subplot(4,2,8)
pcolor(1:16,alt,t)
hold on
plot([-90 90],[3.085 3.085],'r-')
hold on
plot([-90 90],[7.484 7.484],'r-')
set(gca,'xtick',[1:16])
set(gca,'xticklabel',{'0','0.01','1.2','3','5','7','10','15','20','25','30','40','50','60','80','100000'})
shading flat;
axis([1 16 0 19.2])
caxis([0   0.001]);
colormap(a)
colorbar
title('SR Histo Mid Latitude North Pacific JJA GOCCP ')


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








%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        %
%    This part of the program plot the   %
%    zonal average of cloud fraction     %
%    from the SR histo occurence         %
%                                        %
%    Just modify the list with yours     %
%                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all;
%%% Set the Colormap
a=[0.75 0 1; 0.5 0 1; 0 0 1; 0 0.5 1; 0 1 1; 0 1 0.5; 0 1 0; 0.5 1 0; 1 1 0; 1 0.5 0; 1 0 0; 1 0 0.5];

%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/SR_histo/list_SRhist_JJA_CFMIP2_2.1');
 c=textscan(fid,'%s','delimiter','');      
 list=c{1};

ncload([list{1}])

%%% Initialize lat/lon/alt/diagbox
lat(1:length(latitude)+1)=90;
lon(1:length(longitude)+1)=-180;
alt(1:41)=0;
vecdiag=(1:16);   %%% diagbox

%%% Create top latitude / top longitude
for i=1:length(latitude)
  lat(i+1)=lat(i)-180/length(latitude);
end

for i=1:length(longitude)
  lon(i+1)=lon(i)+360/length(longitude);
end


%%% Create top level
for i=1:40
  alt(i+1)=alt(i)+0.48;
end


%%% initialization of the SR Histo vector
mdiag(1:15,1:length(alt_mid),1:length(latitude),1:length(longitude),1:length(list))=0;

%%% Create the SR histo matrix with all values
for i = 1:length(list)
  list{i}
  ncload([list{i}]);

    nanval=cfad_lidarsr532_Occ;   % select the diagSR data
    nanval(nanval<0)=nan;              % set nan values 
    mdiag(:,:,:,:,i)=nanval;             

end

%%% Sum of SR histo occurence
diagSR=squeeze(nansum(mdiag,5));
diagSR=squeeze(nansum(diagSR,4));


%%% initialization of zonalmean vector
x(1:length(alt),1:length(latitude),1:6)=0;

%%% Vector of different cloud detection treshold
tt(1:6)=0;
tt(1)=1.2;
tt(2)=3;
tt(3)=5;
tt(4)=7;
tt(5)=10;
tt(6)=15;


for t=1:6
frac(1:length(alt_mid),1:length(latitude))=0;

%%% Compute the SR histogram occurence in SR histogram fraction
for i=1:length(alt_mid)
  for j=1:length(latitude)
     if nansum(diagSR(1:15,i,j)) ~= 0
       frac(i,j)=nansum(diagSR(2+t:15,i,j))/nansum(diagSR(1:15,i,j));   % computation of the fraction
     end
  end
end


x(1:length(alt_mid),:,t)=frac;



%%% plot 
figure(2)
subplot(3,2,t)
pcolor(latitude',alt,x(:,:,t))
hold on
plot([-90 90],[3.0 3.0],'r-')
hold on
plot([-90 90],[6.3 6.3],'r-')
%set(gca,'YDir','Reverse')
shading flat;
xlabel('LATITUDE')
ylabel('ALTITUDE (KM)')
caxis([0 0.33])
axis([-80 70 0 19])
     colormap(a)
colorbar
if t == 1
title('CF FROM SR Histo TRESHOLD SR=1.2')
     elseif t == 2
title('CF FROM SR Histo TRESHOLD SR=3')
     elseif t == 3
title('CF FROM SR Histo TRESHOLD SR=5')
     elseif t == 4
title('CF FROM SR Histo TRESHOLD SR=7')
     elseif t == 5
title('CF FROM SR Histo TRESHOLD SR=10')
     else t == 6
title('CF FROM SR Histo TRESHOLD SR=15')
end


end
