%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        %
%    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];
jetc=[1 1 1; 0 0 1.; 0 0.3333 1.; 0 0.6667 1.; 0 1. 1.; 0.3333 1. 0.6667; 0.6667 1. 0.3333; 1. 1. 0; 1. 0.6667 0; 1. 0.3333 0; 1. 0 0; 0.666 0 0]

%%% Create the JJA list files

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

ncload([list{1}])

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

%%% 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

%%% Initialize tot low mid high vectors
jjacol(1:length(latitude),1:length(longitude),1:length(list))=0;
jjalow(1:length(latitude),1:length(longitude),1:length(list))=0;
jjamid(1:length(latitude),1:length(longitude),1:length(list))=0;
jjahigh(1:length(latitude),1:length(longitude),1:length(list))=0;


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

  cltcalipso(cltcalipso<-9998)=nan;
  cllcalipso(cllcalipso<-9998)=NaN ;
  clmcalipso(clmcalipso<-9998)=NaN ;
  clhcalipso(clhcalipso<-9998)=NaN ;

  jjacol(:,:,i)=cltcalipso;
  jjalow(:,:,i)=cllcalipso ;
  jjamid(:,:,i)=clmcalipso;
  jjahigh(:,:,i)=clhcalipso;

  clear cltcalipso cllcalipso clmcalipso clhcalipso;
end


meanjjacol=squeeze(nanmean(jjacol,3));   % total cloud vector
meanjjalow=squeeze(nanmean(jjalow,3));   % low cloud
meanjjamid=squeeze(nanmean(jjamid,3));   % mid cloud
meanjjahigh=squeeze(nanmean(jjahigh,3)); % high cloud



figure(1)
%%% Map (hammer graph)
subplot(4,2,1)
m_proj('hammer-aitoff','clongitude',0);
m_pcolor(longitude,latitude,meanjjacol)

% 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
title('TOTAL CLOUD')


subplot(4,2,3)
m_proj('hammer-aitoff','clongitude',0);
m_pcolor(longitude,latitude,meanjjahigh)
 
 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
title('HIGH CLOUD')



subplot(4,2,5)
m_proj('hammer-aitoff','clongitude',0);
m_pcolor(longitude,latitude,meanjjamid)

 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 1])
shading flat;   
title('MID CLOUD')


subplot(4,2,7)
m_proj('hammer-aitoff','clongitude',0);
m_pcolor(longitude,latitude,meanjjalow)

 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 1])
shading flat;   
title('LOW CLOUD')




     % Planisphere (miller graph)
subplot(4,2,2)
m_proj('miller','clongitude',0);
m_pcolor(longitude,latitude,meanjjacol)


 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('TOTAL CLOUD')


subplot(4,2,4)
m_proj('miller','clongitude',0);
m_pcolor(longitude,latitude,meanjjahigh)
 
 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('HIGH CLOUD')



subplot(4,2,6)
m_proj('miller','clongitude',0);
m_pcolor(longitude,latitude,meanjjamid)

 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 1])
shading flat; 
colorbar;  
title('MID CLOUD')


subplot(4,2,8)
m_proj('miller','clongitude',0);
m_pcolor(longitude,latitude,meanjjalow)

 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 1])
colorbar;
shading flat;   
title('LOW CLOUD')

suplabel('Map GOCCP2 AVG JJA 2006-2009','t')
