
% LWCRETotal = 11(COpaqueZT + CThin?ThinZT )

clear
% goccpuse = 0 ==> CERES 2007-2018
% goccpuse = 1 ==> GOCCP-OPAQ 2007-2018
% goccpuse = 2 ==> FLXHR-LIDAR 2007-2010

goccpuse=1;
if goccpuse ==1
 ncload5 /Volumes/T7_2TB/Observations/GOCCP/Map_OPAQ330m_200701-202112_avg_CFMIP2_sat_3.1.4.nc

lat=[-89:2:89];
lon=[-179:2:179];

% Creating the weighting array (airemat) for global mean computations
clear aire
for i=1:length(lat)/2
clear cotes tribase x1

x1=lat(length(lat)/2+i);
x1=x1*2*pi/360;
angle=2*2*pi/360;
hauteur=sin(angle);

cotes=2*sqrt(1-cos(x1)^2);
tribase=2*sqrt(cos(x1)^2-(cos(x1)*cos(angle/2))^2);
aire(i)=hauteur*tribase;
end

airemat(1:length(lon),1:length(lat)/2)=0;
airemat(1,1:length(lat)/2)=aire(length(lat)/2:-1:1);
airemat(1,length(lat)/2+1:length(lat))=aire(1:length(lat)/2);

for i=2:length(lon)
airemat(i,:)=airemat(1,:);
end

calcre{1}=cltcalipso_opaque;
calcre{2}=cltcalipso_opaque_z_se;
calcre{3}=cltcalipso_thin;
calcre{4}=cltcalipso_thin_emis;
calcre{5}=cltcalipso_thin_z_se;

for i=1:5
topcalipso=calcre{i};
topcalipso2(:,:,:)=topcalipso(:,:,1:144);   % 2007-2018

calcre2{i}=topcalipso2;
clear topcalipso topcalipso2
end

% From Vaillant de Guelis et al 2018
LWCRE=11.*(calcre2{1}.*calcre2{2}+calcre2{3}.*calcre2{4}.*calcre2{5});
sizearray=size(LWCRE);
ttime=sizearray(3);

clear cltcalipso_opaque cltcalipso_opaque_z_se cltcalipso_thin cltcalipso_thin_emis cltcalipso_thin_z_se
cltopaq(:,:,:,1)=calcre2{1};  %clt opaq
cltopaq(:,:,:,2)=calcre2{2};  %clt opaq z
cltopaq(:,:,:,3)=calcre2{3};  %clt thin
cltopaq(:,:,:,5)=calcre2{5};  %clt thin emis
cltopaq(:,:,:,4)=calcre2{4}; %clt thin z

for ii=1:5
for t=1:ttime
    cltopaqmean(t,ii)=score2dtot(squeeze(cltopaq(:,:,t,ii)),airemat);
    
end
end

%%% creating climatologies
for ii=1:5
for t=1:12
 cltclim2(:,:,t,ii)=squeeze(nanmean(cltopaq(:,:,t:12:end,ii),3));  
end
end

for ii=1:5
for t=1:sizearray(3)/12
    kk=(t-1)*12;
    cltclim(:,:,kk+1:kk+12,ii)=cltclim2(:,:,:,ii);
end  
end
 clear   cltclim2

for ii=1:5
for t=1:ttime
    cltclimmean(t,ii)=score2dtot(squeeze(cltclim(:,:,t,ii)),airemat);   
end
end
 
%%% Contributions to the feedback page 6 Vaillant de Guelis et al 2018
for t=1:ttime
% Partial derivative terms: replace terms in eq 6 and 7 from those in eq 1
dopaq(t,1)=11.*( cltopaqmean(t,2)+ cltclimmean(t,2))./2;
dopaq(t,2)=11.*( cltopaqmean(t,1)+ cltclimmean(t,1))./2;
dopaq(t,3)=11.*( (cltopaqmean(t,4)+ cltclimmean(t,4))./2 ).* ( (cltopaqmean(t,5)+ cltclimmean(t,5))./2 );
dopaq(t,4)=11.*( (cltopaqmean(t,3)+ cltclimmean(t,3))./2 ).* ( (cltopaqmean(t,5)+ cltclimmean(t,5))./2 );
dopaq(t,5)=11.*( (cltopaqmean(t,3)+ cltclimmean(t,3))./2 ).* ( (cltopaqmean(t,4)+ cltclimmean(t,4))./2 );

clear tmptmp
tmptmp(:,:,1)=cltopaq(:,:,t,2);
tmptmp(:,:,2)=cltclim(:,:,t,2);
tmptmp=squeeze(nansum(tmptmp,3));
masktmp=squeeze(cltopaq(:,:,t,2));
tmptmp(isnan(masktmp))=nan;
dopaqmap(:,:,t,1)=11.*tmptmp./2;

clear tmptmp
tmptmp(:,:,1)=cltopaq(:,:,t,1);
tmptmp(:,:,2)=cltclim(:,:,t,1);
tmptmp=squeeze(nansum(tmptmp,3));
masktmp=squeeze(cltopaq(:,:,t,1));
tmptmp(isnan(masktmp))=nan;
dopaqmap(:,:,t,2)=11.*tmptmp./2;

clear tmptmp
tmptmp(:,:,1)=cltopaq(:,:,t,4);
tmptmp(:,:,2)=cltclim(:,:,t,4);
tmptmp=squeeze(nansum(tmptmp,3));
masktmp=squeeze(cltopaq(:,:,t,4));
tmptmp(isnan(masktmp))=nan;
clear tmptmp2
tmptmp2(:,:,1)=cltopaq(:,:,t,5);
tmptmp2(:,:,2)=cltclim(:,:,t,5);
tmptmp2=squeeze(nansum(tmptmp2,3));
masktmp=squeeze(cltopaq(:,:,t,5));
tmptmp2(isnan(masktmp))=nan;
dopaqmap(:,:,t,3)=11.*tmptmp./2.*tmptmp2./2;

clear tmptmp3
tmptmp3(:,:,1)=cltopaq(:,:,t,3);
tmptmp3(:,:,2)=cltclim(:,:,t,3);
tmptmp3=squeeze(nansum(tmptmp3,3));
masktmp=squeeze(cltopaq(:,:,t,3));
tmptmp3(isnan(masktmp))=nan;
dopaqmap(:,:,t,4)=11.*tmptmp2./2.*tmptmp3./2;
dopaqmap(:,:,t,5)=11.*tmptmp./2.*tmptmp3./2;

end



for ii=1:5
for t=1:ttime
% Delta             
Deltaopaq(t,ii)= cltopaqmean(t,ii)-squeeze(nanmean(cltclimmean(:,ii),1));
Deltaopaqmap(:,:,t,ii)= cltopaq(:,:,t,ii)-squeeze(nanmean(cltclim(:,:,:,ii),3));

% Delta Contributions to LW CRE
DeltaLWopaq(t,ii)=dopaq(t,ii).*Deltaopaq(t,ii);
DeltaLWopaqmap(:,:,t,ii)=dopaqmap(:,:,t,ii).*Deltaopaqmap(:,:,t,ii);
end
end



elseif goccpuse ==0
ncload5('~/Desktop/Desktop_save/CERES/CERES_EBAF-TOA_Ed4.1_Subset_200101-201912.nc')

tt=0;
for t=73:72+144 % 2007 2018
    tt=tt+1;
LWCRE(:,:,tt)=regridmap2(toa_lw_clr_c_mon(:,:,t)-toa_lw_all_mon(:,:,t),360);
end
sizearray=size(LWCRE);
ttime=1:sizearray(3);


elseif goccpuse ==2
lat=[-88.75:2.5:88.75];
lon=[-178.75:2.5:178.75];
ncload5 /Volumes/T7_2TB/Observations/CloudSat-CALIPSO/FLXHR_R05_TOA_LWCRE_200607-201012.nc
 LWCRE(1:144,1:72,1:48)=nan;
 LWCRE(1:144,4:69,:)=Month_LWCRE(:,:,7:end);
  
sizearray=size(LWCRE);
ttime=1:48;
end


ttime=sizearray(3);
% Creating the weighting array for GMST (airemat)
clear aire
for i=1:length(lat)/2
clear cotes tribase x1

x1=lat(length(lat)/2+i);
x1=x1*2*pi/360;
angle=2*2*pi/360;
hauteur=sin(angle);

cotes=2*sqrt(1-cos(x1)^2);
tribase=2*sqrt(cos(x1)^2-(cos(x1)*cos(angle/2))^2);
aire(i)=hauteur*tribase;
end

airemat(1:length(lon),1:length(lat)/2)=0;
airemat(1,1:length(lat)/2)=aire(length(lat)/2:-1:1);
airemat(1,length(lat)/2+1:length(lat))=aire(1:length(lat)/2);

for i=2:length(lon)
airemat(i,:)=airemat(1,:);
end


for t=1:ttime
LWCREmean(t)=score2dtot(squeeze(LWCRE(:,:,t)),airemat);
end

for t=1:12
 LWCREclim2(:,:,t)=squeeze(nanmean(LWCRE(:,:,t:12:end),3));  
 LWCREmeanclim2(t)=squeeze(nanmean(LWCREmean(t:12:end)));
end

for t=1:sizearray(3)/12
    kk=(t-1)*12;
    LWCREclim(:,:,kk+1:kk+12)=LWCREclim2(:,:,:);
    LWCREmeanclim(kk+1:kk+12)=LWCREmeanclim2(:);
end  
clear LWCREclim2

LWCREmapdiff=LWCRE-LWCREclim;


for t=1:sizearray(3)
   LWCRE2(t)=score2dtot(squeeze(LWCRE(:,:,t)),airemat);  
  % CLT2(t)=score2dtot(squeeze(CLT(:,:,t)),airemat);  
end

  LWCREdiff=LWCREmean-LWCREmeanclim;


% Looping on surface temperature datasets to compute GMST anomalies
ttime=sizearray(3);  
for n=1:3
instru=[{'Berkeley'} {'gistemp1200'} {'HadCRUT4'} ];
sat=instru{n}
%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%% CHANGE %%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%

%% BERKELEY
if isequal(sat,'Berkeley')
ncload5 /Volumes/T7_2TB/for_Assia_Juliet/Berkeley.185001-201812.nc
if goccpuse==2
    timed=8*12; % ensures the same time period as CloudSat
else
    timed=0;
end
lat=latitude;
lon=longitude;

% Creating the weighting array for GMST (airemat)
 clear aire
for i=1:length(lat)/2
clear cotes tribase x1

% x1 = angle i
x1=lat(length(lat)/2+i);
x1=x1*2*pi/360;
angle=2*2*pi/360;
hauteur=sin(angle);

cotes=2*sqrt(1-cos(x1)^2);
tribase=2*sqrt(cos(x1)^2-(cos(x1)*cos(angle/2))^2);
aire(i)=hauteur*tribase;
end

airemat(1:length(lon),1:length(lat)/2)=0;
airemat(1,1:length(lat)/2)=aire(length(lat)/2:-1:1);
airemat(1,length(lat)/2+1:length(lat))=aire(1:length(lat)/2);

for i=2:length(lon)
airemat(i,:)=airemat(1,:);
end

% Climatology of each month in an array
tempanomaly=temperature(:,:,length(time)-timed-ttime+1:end-timed); % 2007 2018
for t=1:12
 tempanomalyclim2(:,:,t)=squeeze(nanmean(tempanomaly(:,:,t:12:end),3));  
end

% Creating a climatology array with the full time dimension 
for t=1:sizearray(3)/12
    kk=(t-1)*12;
    tempanomalyclim(:,:,kk+1:kk+12)=tempanomalyclim2;
end   
tempanomaly=tempanomaly-tempanomalyclim;

%%% Weighted global mean surface temperature anomaly
time=1:sizearray(3);
for t=1:length(time) % April 2020
temp(t)=score2dtot(squeeze(tempanomaly(:,:,t)),airemat);
end
clear tempanomaly tempanomalyclim tempanomalyclim2

%% GISTEMP
elseif isequal(sat,'gistemp1200')
ncload5 /Volumes/T7_2TB/for_Assia_Juliet/gistemp1200_GHCNv4_ERSSTv5_188001-202012_2x2.nc
if goccpuse==2
    timed=8*12+24; % ensures the same time period as CloudSat
else
    timed=24; % Going back to December 2018 
end

tempanomaly2=tempanomaly(:,:,length(time)-timed-ttime+1:end-timed); % 2007 2018
clear tempanomaly
tempanomaly=tempanomaly2;
clear tempanomaly2

for t=1:12
 tempanomalyclim2(:,:,t)=squeeze(nanmean(tempanomaly(:,:,t:12:end),3));  
end

for t=1:sizearray(3)/12
    kk=(t-1)*12;
    tempanomalyclim(:,:,kk+1:kk+12)=tempanomalyclim2;
end   
tempanomaly=tempanomaly-tempanomalyclim;

% Creating the weighting array for GMST (airemat)
 clear aire
for i=1:length(lat)/2
clear cotes tribase x1

x1=lat(length(lat)/2+i);
x1=x1*2*pi/360;
angle=2*2*pi/360;
hauteur=sin(angle);

cotes=2*sqrt(1-cos(x1)^2);
tribase=2*sqrt(cos(x1)^2-(cos(x1)*cos(angle/2))^2);
aire(i)=hauteur*tribase;
end

airemat(1:length(lon),1:length(lat)/2)=0;
airemat(1,1:length(lat)/2)=aire(length(lat)/2:-1:1);
airemat(1,length(lat)/2+1:length(lat))=aire(1:length(lat)/2);

for i=2:length(lon)
airemat(i,:)=airemat(1,:);
end

%%% Weighted global mean surface temperature (GMST) anomaly
time=1:sizearray(3);
for t=1:length(time) 
temp(t)=score2dtot(squeeze(tempanomaly(:,:,t)),airemat);
end
clear tempanomaly tempanomalyclim tempanomalyclim2

%% HADCRUT
elseif isequal(sat,'HadCRUT4')      
ncload5 /Volumes/T7_2TB/for_Assia_Juliet/HadCRUT4.6.0.0.185001-201812.nc
if goccpuse==2
    timed=8*12; % ensures the same time period as CloudSat
else
    timed=0; 
end

lat=latitude;
lon=longitude;
tempanomaly=temperature_anomaly(:,:,length(time)-timed-ttime+1:end-timed); % 2007 2018

% Climatology of each month in an array
for t=1:12
 tempanomalyclim2(:,:,t)=squeeze(nanmean(tempanomaly(:,:,t:12:end),3));  
end

% Creating a climatology array with the full time dimension 
for t=1:sizearray(3)/12
    kk=(t-1)*12;
    tempanomalyclim(:,:,kk+1:kk+12)=tempanomalyclim2;
end   
% Anomalies over the 2007-2018 time period
tempanomaly=tempanomaly-tempanomalyclim;

% Creating the weighting array for GMST (airemat)
 clear aire
for i=1:length(lat)/2
clear cotes tribase x1

x1=lat(length(lat)/2+i);
x1=x1*2*pi/360;
angle=2*2*pi/360;
hauteur=sin(angle);

cotes=2*sqrt(1-cos(x1)^2);
tribase=2*sqrt(cos(x1)^2-(cos(x1)*cos(angle/2))^2);
aire(i)=hauteur*tribase;
end

airemat(1:length(lon),1:length(lat)/2)=0;
airemat(1,1:length(lat)/2)=aire(length(lat)/2:-1:1);
airemat(1,length(lat)/2+1:length(lat))=aire(1:length(lat)/2);

for i=2:length(lon)
airemat(i,:)=airemat(1,:);
end

%%% global surface temperature mean anomaly
for t=1:ttime % April 2020
temp(t)=score2dtot(squeeze(tempanomaly(:,:,t)),airemat);
end
clear tempanomaly tempanomalyclim tempanomalyclim2
%%
end % case to select temperature dataset
 

%% Now computing LW Feedback
% Reshaping LWCRE array for regress
LWCREmapdiff=reshape(LWCREmapdiff,sizearray(1)*sizearray(2),sizearray(3));
   
 % Regress against GMST anomaly at each grid box
 for i=1:sizearray(1)*sizearray(2)
   [coeflw(n,i) coeflwint]=regress(squeeze(LWCREmapdiff(i,:))',temp(:)); 
 end


% Contributions of each parameter to the feedback for OPAQ only 
if goccpuse == 1
  for ii=1:5
  % Delta Contributions to the total LW feedback
   tmptmp=squeeze(DeltaLWopaqmap(:,:,:,ii)); % temporary array for each parameter
   tmptmp=reshape(tmptmp,sizearray(1)*sizearray(2),sizearray(3)); % reshaping to make it easy to loop

    % Regression against GMST (interannual feedback)
    for i=1:sizearray(1)*sizearray(2)
     [coeflwopaq(n,i,ii) coeflwint]=regress(squeeze(tmptmp(i,:))',temp(:)); 
    end
 
  end 
end

 clear tmp tempanomaly sst
 
   
end %% n loop

% Reshaping the total LW feedback array in a map array
coeflw=reshape(coeflw,3,sizearray(1),sizearray(2));

% Reshaping the LW feedback array of each contribution in a map array
if goccpuse ==1
coeflwopaq=reshape(coeflwopaq,3,sizearray(1),sizearray(2),5);
end

