%  This program plot zonalmean of        %
%  of CloudSat Obs.                      %
%                                        %
%  Just modify the list with yours       %
%                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


close ('all');
clear all;

%%% colormap vector
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];


ncload('/bdd/CFMIP/GOCCP/CloudSat/CloudSat_Reflectivity/CloudSat_TotalCloudCoverage_200606-201011_ver0.cdf');


%%% Create the matrix
clt_cloudsat(clt_cloudsat==-999)=nan;  % Set NaN values
clt=squeeze(nanmean(clt_cloudsat,3));  % Average over the time


%%% Create the latitude vector
lat(1:90)=0;
lat(1)=-89;

for i =1:89 
lat(i+1)=lat(i)+2;
end


%%% New matrix to match to m_proj graphic
clt2(1:180,1:90)=0;
clt2(:,5:86)=clt;



%%% Create the plot
figure;

%%% Map (hammer graph)
subplot(2,1,1)
m_proj('hammer-aitoff','clongitude',0);
m_pcolor(longitude,lat',clt2') 
% some latitude line to make easier the reading of the plot
 y(1:length(longitude))=0;
m_line(longitude,y','color','black')
 y(1:length(longitude))=30;
m_line(longitude,y','color','black')
 y(1:length(longitude))=60;
m_line(longitude,y','color','black')
 y(1:length(longitude))=-60;
m_line(longitude,y','color','black')
 y(1:length(longitude))=-30;
m_line(longitude,y','color','black')

m_coast('color','black') ; % pour les cotes continents
caxis([0 1])
m_grid('xtick', [], 'ytick', []);
shading flat;     % pour ne pas mettre de grille
colorbar
title('CloudSat Total Cloud Coverage')


%%% Load the file
ncload('/bdd/CFMIP/GOCCP/CloudSat/CloudSat_Reflectivity/CloudSat_RadarReflectivity_CFAD_200606-201011_ver0.cdf');


%%% Create the matrix
cfadDbze94(cfadDbze94==-999)=nan;  % NaN values
 
%%% average to make the 2dim matrix 
cloud1=squeeze(nanmean(cfadDbze94,5)); 
cloud2=squeeze(nansum(cloud1,4));
cloud3=squeeze(nanmean(cloud2,1));


%%% Create the plot
subplot(2,1,2)
imagesc(latitude,alt40/1000,cloud3');
set(gca,'ydir','normal')
caxis([0 0.33])
axis([-80 80 0 16])
colorbar
colormap(a)
ylabel('ALTITUDE (KM)')
xlabel('LATITUDE')
title('CloudSat Cloud Fraction 3D')
