clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ttime=120; % for 2007 - 2016
coef=1; % for Sc Cu cloud fraction (use 0.7 for feedback computation)

% Reanalyse omega500
ncload5 /Volumes/T7_2TB/Reanalyses/ERA/omega500_ERA5-MERRA2-NCEP_200701-201612_2.5x2.5.nc

%%% GOCCP %%%
ncload5 ~/Downloads/MapLowMidHigh330m_200701-201612_avg_CFMIP2.5_sat_2.9.nc
% regrid 0 to 360
for t=1:ttime
  cllobs(:,:,t)=regridmap2(cllcalipso(:,:,t),144);  
end
clear cllcalipso 

%%% CERES %%%
ncload5('~/Desktop/Desktop_save/CERES/CERES_EBAF-TOA_Ed4.1_Subset_200701-201612_2.5x2.5.nc')
creswobs=double(toa_sw_clr_c_mon)-double(toa_sw_all_mon);
crelwobs=double(toa_lw_clr_c_mon)-double(toa_lw_all_mon);

ttime=120;
lat_obs=[-88.75:2.5:88.75];

% mask land/ocean
ncload5 /Volumes/T7_2TB/Observations/Landmask_200701-201612_2.5x2.5.nc 

% Tropics only
tata=find(lat_obs>35 | lat_obs<-35);
cllobs(:,tata,:)=nan;
wap500obs(:,tata,:)=nan;
creswobs(:,tata,:)=nan;

% ocean only
cllobs(landmask>0)=nan;
wap500obs(landmask>0)=nan;
creswobs(landmask>0)=nan;

% Subsidence only
cllobs(wap500obs<0)=nan;
creswobs(wap500obs<0)=nan;

% Sc and Cu regimes (cloud and cre)
[scmeanobsreg,cumeanobsreg,tmp1,tmp2,tmp3,tmp4,tmp5]=sccumask(cllobs,wap500obs,lat_obs,coef); % cloud
[scmeanobsreg2,cumeanobsreg2]=sccumask2(cllobs,creswobs,wap500obs,lat_obs,coef); % CRE

% Sc and Cu mean and frequency (cloud and cre)
[scmeanobs,cumeanobs,scfreqobs,cufreqobs,scfreqnormobs,cufreqnormobs,ratiomaskobs]=sccumasknoregimes(cllobs,wap500obs,lat_obs,coef); % cloud
[scmeanobs2,cumeanobs2]=sccumask2noregimes(cllobs,creswobs,wap500obs,lat_obs,coef); % CRE
clear tmp1 tmp2 tmp3 tmp4 tmp5 % dummy variables

% ensure nan values
scfreqobs(isnan(cllobs))=nan;
cufreqobs(isnan(cllobs))=nan;
scmeanobs(isnan(cllobs))=nan;
cumeanobs(isnan(cllobs))=nan;
scmeanobsreg(isnan(cllobs))=nan;
cumeanobsreg(isnan(cllobs))=nan;

% Load CASCCAD feedbacks and compute Sc Cu components using the mask method
ncload5 /Volumes/T7_2TB/Observations/Feedback_CALIPSO-CASCCAD_CERES_CFMIP2.5.nc

% 3 corresponds to p4K feedbacks
allobsfdk=squeeze(totfeedback(:,:,3)); % total
scobsfdk=squeeze(scfeedback(:,:,3));   % sc
cuobsfdk=squeeze(cufeedback(:,:,3));   % cu

% recompute the sc cu cloud fraction with coef=0.7
coef=0.7;
[scmeanobs,cumeanobs,scfreqobs,cufreqobs,scfreqnormobs,cufreqnormobs,ratiomaskobs]=sccumasknoregimes(cllobs,wap500obs,lat_obs,coef);

% weight the total low cloud feedback with sc cu cloud fraction
scobsfdkmask=squeeze(nanmean(scmeanobs./(cllobs),3)).*allobsfdk;
cuobsfdkmask=squeeze(nanmean(cumeanobs./(cllobs),3)).*allobsfdk;






