%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                        %
%    This program plot zonal average of  %
%    cloud / clear / undef fraction of   %
%    GOCCP files.                        %
%                                        %
%    Just modify the list with yours     %
%                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



clear all;
close ('all');

%%% Create the JJA list files

 fid=fopen('/bdd/CFMIP/GOCCP/3D_CloudFraction/list_3DCF_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


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


zonalmeanclear(1:length(alt),1:length(latitude),1:length(list))=0;
zonalmeancloud(1:length(alt),1:length(latitude),1:length(list))=0;
zonalmeanundef(1:length(alt),1:length(latitude),1:length(list))=0;


%%% Average on JJA files
for i = 1:length(list)
  list{i}
  ncload([list{i}]);

%%% Clear Fraction
nanval=clrcalipso;
nanval(nanval<0)=nan;  % met ttes les valeurs de nanval <0 egales à nan
zonalmeanclear(1:length(alt_mid),1:length(latitude),i)=squeeze(nanmean(nanval(:,:,:),3));  % moyenne sur la 2e colonne de tpr2obs puis squeeze
clear nanval

%%% Cloud Fraction
nanval=clcalipso ;
nanval(nanval<0)=nan;  % met ttes les valeurs de nanval <0 egales à nan
zonalmeancloud(1:length(alt_mid),1:length(latitude),i)=squeeze(nanmean(nanval(:,:,:),3));  % moyenne sur la 2e colonne de tpr2obs puis squeeze
clear nanval

%%% Undefined Fraction
nanval=uncalipso;
nanval(nanval<0)=nan;  % met ttes les valeurs de nanval <0 egales à nan
zonalmeanundef(1:length(alt_mid),1:length(latitude),i)=squeeze(nanmean(nanval(:,:,:),3));  % moyenne sur la 2e colonne de tpr2obs puis squeeze
clear nanval

clear uncalipso clcalipso clrcalipso


end


zclear=squeeze(nanmean(zonalmeanclear,3));
zcloud=squeeze(nanmean(zonalmeancloud,3));
zundef=squeeze(nanmean(zonalmeanundef,3));





%%% plot
figure(1)

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

subplot(2,2,1)
pcolor(latitude,alt,zcloud)
hold on
w(1:length(latitude))=3.36;
y(1:length(latitude))=7.2;
plot(latitude,w,'r',latitude,y,'r')
shading flat;
xlabel('LATITUDE')
ylabel('ALTITUDE (KM)')
caxis([0 0.33])
axis([-80 80 0 16])
colormap(a)
colorbar
title('CLOUDY FRACTION JJA')


subplot(2,2,3)
pcolor(latitude,alt,zclear)
shading flat;
xlabel('LATITUDE')
ylabel('ALTITUDE (KM)')
caxis([0 1])
axis([-80 80 0 16])
colormap(a)
colorbar
title('CLEAR FRACTION JJA')


subplot(2,2,4)
pcolor(latitude,alt,zundef)
shading flat;
caxis([0 0.33])
axis([-80 80 0 16])
xlabel('LATITUDE')
ylabel('ALTITUDE (KM)')
colorbar
title('UNDEFINED FRACTION JJA')




%%%%%%%%%%%% SAME PROCESS FOR JFM DATA %%%%%%%%%%%%%%%%%%%

clear all;



%%% Create the JFM list files

 fid=fopen('/bdd/CFMIP/GOCCP/3D_CloudFraction/list_3DCF_JFM_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


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


zonalmeanclear(1:length(alt),1:length(latitude),1:length(list))=0;
zonalmeancloud(1:length(alt),1:length(latitude),1:length(list))=0;
zonalmeanundef(1:length(alt),1:length(latitude),1:length(list))=0;


%%% Average on JFM files
for i = 1:length(list)
  list{i}
  ncload([list{i}]);

%%% Clear Fraction
nanval=clrcalipso;
nanval(nanval<0)=nan;  % met ttes les valeurs de nanval <0 egales à nan
zonalmeanclear(1:length(alt_mid),1:length(latitude),i)=squeeze(nanmean(nanval(:,:,:),3));  % moyenne sur la 2e colonne de tpr2obs puis squeeze
clear nanval

%%% Cloud Fraction
nanval=clcalipso ;
nanval(nanval<0)=nan;  % met ttes les valeurs de nanval <0 egales à nan
zonalmeancloud(1:length(alt_mid),1:length(latitude),i)=squeeze(nanmean(nanval(:,:,:),3));  % moyenne sur la 2e colonne de tpr2obs puis squeeze
clear nanval

%%% Undefined Fraction
nanval=uncalipso;
nanval(nanval<0)=nan;  % met ttes les valeurs de nanval <0 egales à nan
zonalmeanundef(1:length(alt_mid),1:length(latitude),i)=squeeze(nanmean(nanval(:,:,:),3));  % moyenne sur la 2e colonne de tpr2obs puis squeeze
clear nanval

clear uncalipso clcalipso clrcalipso


end


zclear=squeeze(nanmean(zonalmeanclear,3));
zcloud=squeeze(nanmean(zonalmeancloud,3));
zundef=squeeze(nanmean(zonalmeanundef,3));



figure(1)
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];

subplot(2,2,2)
pcolor(latitude,alt,zcloud)
hold on
w(1:length(latitude))=3.36;
y(1:length(latitude))=7.2;
plot(latitude,w,'r',latitude,y,'r')
shading flat;
xlabel('LATITUDE')
ylabel('ALTITUDE (KM)')
caxis([0 0.33])
axis([-80 80 0 16])
colormap(a)
colorbar
title('CLOUD FRACTION JFM')

