clear
ncload5 /Volumes/T7_2TB/Observations/PHACT/PHACT2D_200701-202012_avg_2.5x5_3.1.2_v1.5.nc
titl=[{'Ice only'} {'Ice seeding'} {'Ice and liq'} {'Liq single'} {'Liq multi'} {'Liq contiguous'} ...
      {'Mix contiguous'} {'Mix multi'} {'Likely mix contiguous'} {'Likely mix multi'} ];

% Bug that needs to be fix in next version some NaN are = 9999
ice_O_S(ice_O_S==9999)=nan;
ice_O_ML(ice_O_ML==9999)=nan;

% Ice only total cloud cover
ice2(:,:,:,1)=ice_O_S;
ice2(:,:,:,2)=ice_O_ML;
iceall=squeeze(nansum(ice2,4));
iceall(isnan(ice_O_S)|isnan(ice_O_ML))=nan;
[min(iceall(:)) max(iceall(:))]

% total cloud cover = ice only + liquid-bearing
liq(liq==9999)=nan; % liquid-bearing cloud cover
tot(:,:,:,1)=liq;
tot(:,:,:,2)=iceall;
tot=squeeze(nansum(tot,4));
tot(isnan(liq)|isnan(iceall))=nan;

clear ice2 liq

ncload5 /Volumes/T7_2TB/Observations/PHACT/PHACT2Dthin_200701-202012_avg_2.5x5_3.1.2_v1.5.nc
ice_O_S_thin(ice_O_S_thin==9999)=nan;
ice_O_ML_thin(ice_O_ML_thin==9999)=nan;
ice2(:,:,:,1)=ice_O_S_thin;
ice2(:,:,:,2)=ice_O_ML_thin;
iceallthin=squeeze(nansum(ice2,4));
iceallthin(isnan(ice_O_S_thin)|isnan(ice_O_ML_thin))=nan;
[min(iceallthin(:)) max(iceallthin(:))]


liq_thin(liq_thin==9999)=nan;
totthin(:,:,:,1)=liq_thin;
totthin(:,:,:,2)=iceallthin;
totthin=squeeze(nansum(totthin,4));
totthin(isnan(liq_thin)|isnan(iceallthin))=nan;

% opaque cloud cover = total - thin
totopaq(:,:,:,1)=-totthin;
totopaq(:,:,:,2)=tot;
totopaq=squeeze(nansum(totopaq,4));
totopaq(isnan(tot)|isnan(totthin))=nan;


% Visualization
load coast
panoply=colormap_cpt('~/Desktop/Desktop_save/YOTC/panoply.cpt');

figure
subplot(2,2,1)
pcolor(lon,lat,squeeze(nanmean(totopaq(:,:,13:108),3))');shading flat
hold on
plot(long,lati,'black-','linewidth',2); plot(long+360,lati,'black-','linewidth',2);caxis([-1.8 1.8])
caxis([0 1])
set(gca,'layer','top')
grid on
daspect([1 1 1])
title('PHACT2D','Interpreter','none','Fontname','times','fontsize',14)
ylabel('Opaque','fontname','times','color','black','fontsize',14)

subplot(2,2,3)
pcolor(lon,lat,squeeze(nanmean(totthin(:,:,13:108),3))');shading flat
hold on
plot(long,lati,'black-','linewidth',2); plot(long+360,lati,'black-','linewidth',2);caxis([-1.8 1.8])
caxis([0 1])
set(gca,'layer','top')
grid on
daspect([1 1 1])
ylabel('Thin','fontname','times','color','black','fontsize',14)

set(gcf,'color','white')
colormap(panoply)
tutu=colorbar('east','ytick',[-0:1/16:1],'fontname','times','fontsize',12);
ylabel(tutu,'Cloud Fraction','fontname','times','color','black','fontsize',14)



% Now looking at the OPAQ product
clear

ncload5 /Volumes/T7_2TB/Observations/GOCCP/Map_OPAQ330m_200606-201612_avg_CFMIP2_sat_3.1.2.nc
load coast

% 20 is January 2008 and 108+7=115 December 2015
subplot(2,2,2)
pcolor(longitude,latitude,squeeze(nanmean(cltcalipso_opaque(:,:,20:108+7),3))');shading flat
hold on
plot(long,lati,'black-','linewidth',2); plot(long+360,lati,'black-','linewidth',2);caxis([-1.8 1.8])
caxis([0 1])
set(gca,'layer','top')
grid on
daspect([1 1 1])
title('OPAQ','Interpreter','none','Fontname','times','fontsize',14)

subplot(2,2,4)
pcolor(longitude,latitude,squeeze(nanmean(cltcalipso_thin(:,:,20:108+7),3))');shading flat
hold on
plot(long,lati,'black-','linewidth',2); plot(long+360,lati,'black-','linewidth',2);caxis([-1.8 1.8])
caxis([0 1])
set(gca,'layer','top')
grid on
daspect([1 1 1])

