%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        %
%    This program plot Map low mid high  %
%    and total cloud of GOCCP files.     %
%                                        %
%    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];

%%% Load files
ncload clisccp_198307-200806.nc
ncload ctpisccp_198307-200806.nc
ncload cttisccp_198307-200806.nc 
ncload albisccp_198307-200806.nc
ncload cltisccp_198307-200806.nc

%%% Parameter vectors
cl(1:7,1:8)=0;
tau1(1:7)=0;
plev7(1:8)=0;

tau1(1:6)=tau_bnds(1,:);
tau1(7)=tau_bnds(2,6);
plev(1:7)=plev7_bnds(1,:);
plev(8)=plev7_bnds(2,7);


%%% Nan value initialization
clisccp(clisccp==-1000)=nan;
albisccp(albisccp==-1000)=nan;
cttisccp(cttisccp==-1000)=nan;
ctpisccp(ctpisccp==-1000)=nan;
cltisccp(cltisccp==-1000)=nan;


%%% Average over the 198307-200806 period
cl(1:6,1:7)=squeeze(nanmean(squeeze(nanmean(squeeze(nanmean(clisccp,1)),3)),3));
alb=squeeze(nanmean(albisccp,1));
ctt=squeeze(nanmean(cttisccp,1));
ctp=squeeze(nanmean(ctpisccp,1));
clt=squeeze(nanmean(cltisccp,1));



%%% Plot the results
figure(1)

subplot(2,2,1)
m_proj('miller',0);
m_pcolor(longitude,latitude,alb)
 
 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('Albedo')



subplot(2,2,2)
m_proj('miller',0);
m_pcolor(longitude,latitude,ctt)

 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
m_grid('xtick', [], 'ytick', []);
caxis([0 330])
shading flat; 
colorbar;  
title('Cloud Top Temperature')


subplot(2,2,3)
m_proj('miller',0);
m_pcolor(longitude,latitude,ctp)

 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
m_grid('xtick', [], 'ytick', []);
colormap(a)
caxis([0 110000])
colorbar;
shading flat;   
title('Cloud Top Pressure')


subplot(2,2,4)
m_proj('miller',0);
m_pcolor(longitude,latitude,clt)

 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
m_grid('xtick', [], 'ytick', []);
colormap(a)
caxis([0 100])
colorbar;
shading flat;   
title('Total Cloud Fraction')
 suplabel('GCM simulator-oriented ISCCP Data for 198307-200806','t')

