
% Code used for Wills and DiNezio, submitted. Submitted for the purposes of
% peer review. Will be archived on Github upon final publication. 

% Note that the data files needed to run this code are too large to upload 
% (~500 GB, larger than 50 GB Zenodo upload limit). However, they are
% simply netcdf files for all included variables and models, in standard
% ESGF output format, but merged into 1 file per variable and model.
% Additionally, the data has been interpolated to 1° resolution for
% precipitation and 2.5° resolution for all other variables. The
% pre-processing workflow is documented in Brunner et al. 2020, https://doi.org/10.5281/zenodo.3734127

dataPath = '/Users/rjnglin/Data/ITCZ_SST_trend_biases/';

varnam = 'pr'; 

[precip_gpcp,lat,lon] = get_avg_field_nd('gpcp.precip.mon.ltm.1991-2020.g1.nc',{'precip','lat','lon'});
[precip_imerg,lat,lon] = get_avg_field_nd('IMERG_total_precipitation_mon_global_2001-2020.ltm.g1.nc',{'pr','lat','lon'});

%% Load and process CMIP and OBS precipitation data

dir = [dataPath,'cmip6/pr/'];
%dir = [dataPath,'cmip5/pr/'];
files = get_files(dir);

clear models
for i = 1:length(files)
    tmp = strsplit(files{i},'/');
    tmp2 = strsplit(tmp{end},'_');
    models{i} = tmp2{3};
end

precip_model = zeros([length(lon), length(lat),12,length(models)]);
for i = 1:length(models)
    field = get_avg_field_nd(files{i},varnam);
    if strcmp(varnam,'pr')
        units = 3600*24; % kg/m2/s to mm/day
        field = field.*units;
        field(field>1e10) = nan;
        if strcmp(models{i},'CIESM')
            field = field.*1000;
        end
    end
%     % the following 3 lines are only used to additionally compute 40-yr (not used in paper)
%     % precipitation trends (i.e., not used in the ITCZ bias analysis itself)
%     field_DJFMAM = reshape(field(:,:,end-480+1:end),[360 180 12 40]);
%     field_DJFMAM = squeeze(mean(field_DJFMAM(:,:,1:12,:),3));
%     [~,field_DJFMAM_trend_model(:,:,i)] = detrend(field_DJFMAM,3,1);
    try
        field = field(:,:,end-360+1:end); % 1985-2014
    catch 
        disp('Simulation has less than 30 years.')
        field = field(:,:,end-348+1:end);
    end
    % field = field(:,:,1:600); 50-yr trends; changes the ordering of models, but not which are best 6
    for j = 1:12
        precip_model(:,:,j,i) = nanmean(field(:,:,j:12:end),3);
    end
end

%% Calculation of Double-ITCZ bias index, including specification of the averaging region, season, and threshold

% was January-June, but now trying switching everything to December-May
%is_season = 1:6;
%is_season = [2 3 4];
is_season = [12 1 2 3 4 5];
%is_season = [3 4 5];

threshold = 4.2;
%rmse_threshold = 10;

lons = closest(lon,225.01):closest(lon,269.99);
%lons = closest(lon,225.01):closest(lon,274.99); earlier presented results were with this choice
lats = closest(lat,-14.99):closest(lat,14.99);
%latsSH = lats(1:15); % same top 6 models with this choice (also with JFMAMJ instead of DJFMAM); just alternates which ones are separated from the rest of the pack
latsSH = closest(lat,-14.5):closest(lat,-2.5); 

mse_gpcp = (precip_model-repmat(precip_gpcp,[1 1 1 length(models)])).^2;
msen_gpcp = ((precip_model-repmat(precip_gpcp,[1 1 1 length(models)]))./repmat(precip_gpcp,[1 1 1 length(models)])).^2;
precip_diff = precip_model-repmat(precip_gpcp,[1 1 1 length(models)]);
mse_imerg = (precip_model-repmat(precip_imerg,[1 1 1 length(models)])).^2;
msen_imerg = ((precip_model-repmat(precip_imerg,[1 1 1 length(models)]))./repmat(precip_gpcp,[1 1 1 length(models)])).^2;
precip_diff = precip_model-repmat(precip_imerg,[1 1 1 length(models)]);

precip_gpcp_smax = max(precip_gpcp,[],3);
precip_imerg_smax = max(precip_imerg,[],3);
precip_model_smax = max(precip_model,[],3);
precip_model_SHmax_smax_lmax = max(max(max(precip_model(lons,latsSH,:,:),[],2),[],3),[],1);
precip_gpcp_SHmax_smax_lmax = max(max(max(precip_gpcp(lons,latsSH,:),[],2),[],3),[],1);
precip_imerg_SHmax_smax_lmax = max(max(max(precip_imerg(lons,latsSH,:),[],2),[],3),[],1);

precip_gpcp_DJFMAM = mean(precip_gpcp(:,:,is_season),3);
precip_imerg_DJFMAM = mean(precip_imerg(:,:,is_season),3);
precip_model_DJFMAM = mean(precip_model(:,:,is_season,:),3);
precip_model_DJFMAM_SHmax_lmax = max(max(mean(precip_model(lons,latsSH,is_season,:),3),[],2),[],1);
precip_gpcp_DJFMAM_SHmax_lmax = max(max(mean(precip_gpcp(lons,latsSH,is_season),3),[],2),[],1);
precip_imerg_DJFMAM_SHmax_lmax = max(max(mean(precip_imerg(lons,latsSH,is_season),3),[],2),[],1);

precip_model_DJFMAM_SHmax_lmean = max(mean(mean(precip_model(lons,latsSH,is_season,:),3),1),[],2);
precip_gpcp_DJFMAM_SHmax_lmean = max(mean(mean(precip_gpcp(lons,latsSH,is_season),3),1),[],2);
precip_imerg_DJFMAM_SHmax_lmean = max(mean(mean(precip_imerg(lons,latsSH,is_season),3),1),[],2);

precip_model_DJFMAM_SHmean_lmean = mean(mean(mean(precip_model(lons,latsSH,is_season,:),3),2),1);
precip_gpcp_DJFMAM_SHmean_lmean = mean(mean(mean(precip_gpcp(lons,latsSH,is_season),3),2),1);
precip_imerg_DJFMAM_SHmean_lmean = mean(mean(mean(precip_imerg(lons,latsSH,is_season),3),2),1);

precip_model_DJFMAM_SHmax_120 = max(mean(precip_model(closest(lon,240),latsSH,is_season,:),3),[],2);
precip_gpcp_DJFMAM_SHmax_120 = max(mean(precip_gpcp(closest(lon,240),latsSH,is_season),3),[],2);
precip_imerg_DJFMAM_SHmax_120 = max(mean(precip_imerg(closest(lon,240),latsSH,is_season),3),[],2);

% To switch to other half year if needed
% precip_gpcp_JFMAMJ = mean(precip_gpcp(:,:,7:12),3);
% precip_imerg_JFMAMJ = mean(precip_imerg(:,:,7:12),3);
% precip_model_JFMAMJ = mean(precip_model(:,:,7:12,:),3);
% precip_model_JFMAMJ_SHmax_lmax = max(max(mean(precip_model(lons,latsSH,7:12,:),3),[],2),[],1);
% precip_gpcp_JFMAMJ_SHmax_lmax = max(max(mean(precip_gpcp(lons,latsSH,7:12),3),[],2),[],1);
% precip_imerg_JFMAMJ_SHmax_lmax = max(max(mean(precip_imerg(lons,latsSH,7:12),3),[],2),[],1);

precip_index = squeeze(precip_model_DJFMAM_SHmax_lmax);
rmse_imerg = mean_in_a_box(lon,lat,mean(mse_imerg,3),[-30 30],[0 360]);

iselect = precip_index<threshold;
iselect_biased = precip_index>threshold;


%% Plotting precipitaiton climatology in IMERG, high-bias models, and low-bias models (Fig. 1)

plot_field(lon,lat,precip_imerg_DJFMAM,linspace(0,10.5,22),magma(21)); set(gca,'xlim',[100 300]); set(gca,'ylim',[-20 20])
set(gca,'ytick',-20:10:20); set(gca,'yticklabel',{'20°S','10°S','EQ','10°N','20°N'})
set(gca,'xtick',60:60:300); set(gca,'xticklabel',{'60°E','120°E','180°','120°W','60°W'})
hold on; plot([225 270],[-2 -2],'w','linewidth',1.5); plot([270 270],[-15 -2],'w','linewidth',1.5); plot([225 225],[-15 -2],'w','linewidth',1.5);  hold on; plot([225 270],[-15 -15],'w','linewidth',1.5);
coast_vector
pretty_figure(600,200,'none','none','none','none',16);

plot_field(lon,lat,nanmean(precip_model_DJFMAM(:,:,:,iselect_biased),4),linspace(0,10.5,22),magma(21)); set(gca,'xlim',[100 300]); set(gca,'ylim',[-20 20])
set(gca,'ytick',-20:10:20); set(gca,'yticklabel',{'20°S','10°S','EQ','10°N','20°N'})
set(gca,'xtick',60:60:300); set(gca,'xticklabel',{'60°E','120°E','180°','120°W','60°W'})
hold on; plot([225 270],[-2 -2],'w','linewidth',1.5); plot([270 270],[-15 -2],'w','linewidth',1.5); plot([225 225],[-15 -2],'w','linewidth',1.5);  hold on; plot([225 270],[-15 -15],'w','linewidth',1.5);
coast_vector
pretty_figure(600,200,'none','none','none','none',16);

plot_field(lon,lat,nanmean(precip_model_DJFMAM(:,:,:,iselect),4),linspace(0,10.5,22),magma(21)); set(gca,'xlim',[100 300]); set(gca,'ylim',[-20 20])
set(gca,'ytick',-20:10:20); set(gca,'yticklabel',{'20°S','10°S','EQ','10°N','20°N'})
set(gca,'xtick',60:60:300); set(gca,'xticklabel',{'60°E','120°E','180°','120°W','60°W'})
hold on; plot([225 270],[-2 -2],'w','linewidth',1.5); plot([270 270],[-15 -2],'w','linewidth',1.5); plot([225 225],[-15 -2],'w','linewidth',1.5);  hold on; plot([225 270],[-15 -15],'w','linewidth',1.5);
coast_vector
pretty_figure(600,200,'none','none','none','none',16);

%% Plotting Double-ITCZ bias index in models and observations (Fig. 1d)

figure; plot(1:length(models),sort(precip_index),'ko','markersize',4.25,'markerfacecolor','k','linewidth',0.5)
hold on; plot([1 70],[precip_imerg_DJFMAM_SHmax_lmax precip_imerg_DJFMAM_SHmax_lmax],'k--')
hold on; plot([1 70],[precip_gpcp_DJFMAM_SHmax_lmax precip_gpcp_DJFMAM_SHmax_lmax],'k--')
pretty_figure(400,300,'Model','Max DJFMAM Precip. (mm day^{-1})','none','none',16);
set(gca,'ylim',[0 13]); set(gca,'ytick',0:2:14)
set(gca,'xlim',[0 70]); set(gca,'xtick',0:10:70)
tmp = sort(precip_index(iselect));
tmp = tmp(tmp<threshold);
colors = lines(2);
hold on; plot(1:length(tmp),tmp,'ko','markersize',4.25,'markerfacecolor',colors(1,:),'linewidth',0.5)

%% Other Figures Illustrating Double-ITCZ Bias (Extended Data Figure 1)

plot_field(lon,lat,sqrt(nanmean(nanmean(mse_imerg,4),3))./mean(precip_imerg,3),linspace(0,16,17),magma(16)); set(gca,'ylim',[-20 20]); set(gca,'xlim',[90 360+90])%set(gca,'xlim',[100 300]);
set(gca,'ytick',-20:10:20); set(gca,'yticklabel',{'20°S','10°S','EQ','10°N','20°N'})
set(gca,'xtick',120:60:420); set(gca,'xticklabel',{'120°E','180°','120°W','60°W','0°','60°E'})
%hold on; plot([225 275],[-2 -2],'w','linewidth',1.5); plot([275 275],[-15 -2],'w','linewidth',1.5); plot([225 225],[-15 -2],'w','linewidth',1.5);  hold on; plot([225 275],[-15 -15],'w','linewidth',1.5);
hold on; plot([225 270],[-2 -2],'w','linewidth',1.5); plot([270 270],[-15 -2],'w','linewidth',1.5); plot([225 225],[-15 -2],'w','linewidth',1.5);  hold on; plot([225 270],[-15 -15],'w','linewidth',1.5);
coast_vector
pretty_figure(800,200,'none','none','none','none',16);

figure; plot(1:12,sqrt(mean_in_a_box(lon,lat,nanmean(mse_imerg,4),[-14.5 -2.5],[225 275])),'k','linewidth',2); set(gca,'xtick',1:12); set(gca,'xticklabel',{'J','F','M','A','M','J','J','A','S','O','N','D'})
pretty_figure(300,150,'none','none','none','none',16); set(gca,'ytick',0:2:10); set(gca,'xlim',[1 12]); set(gca,'xgrid','on')
ylabel('RMSE (mm day^{-1})','fontsize',16)
title('SE Pacific','fontsize',14)
set(gca,'xgrid','on')

%% Observational SST pattern and mechanism analysis (Data Handling and initialization)

season = 'Amean'; 
trendlength = 40;

% this is for a pattern correlation that is not currently used in the paper
boxlat = [-45 15];
boxlon = [140 280];

[sst,lat25,lon25] = get_avg_field_nd('ersst5.sst.mnmean.185401-202405.g025.nc',{'sst','lat','lon'});
sst(abs(sst)>1e10) = nan;
sst = sst(:,:,1:end-5);

[u10,v10] = get_avg_field_nd('era5_u10_v10_1940-2023_g025.nc',{'u10','v10'});
slp = get_avg_field_nd('era5.msl.194001-202312.g025.nc','msl');
u10(abs(u10)>1000) = nan;
v10(abs(v10)>1000) = nan;
slp(slp<0) = nan;

% pr = get_avg_field_nd('gpcp.mon.mean.197901-202403.g1.nc','precip');
% pr = pr(:,:,1:end-3);
pr = 1000.*get_avg_field_nd('era5.tp.194001-202312.g1.nc','tp');
pr(abs(pr)>1000) = nan;

sstcut = sst(:,:,553:end);
sstcut2 = sst(:,:,end-1008+1:end);
s = size(sstcut2);
sst_clim = mean(reshape(sstcut2(:,:,1:50*12),[s(1) s(2) 12 50]),4);
sst_onepoint = mean_in_a_box(lon25,lat25,sstcut2,[-20 -15],[258 262]);

months = repmat(1:12,[1 length(sst_onepoint)/12]); years = floor(1:1/12:length(sst_onepoint)/12+0.99);
sst_onepoint = monthly_to_seasonal(months,years,sst_onepoint,[9:12 1:2],2); % was [12 1 2], was [9:12 1 2]
sst_onepoint = sst_onepoint - mean(sst_onepoint);
sst_onepoint = sst_onepoint./std(sst_onepoint);
ist = 13;
sst_centered = monthly_to_seasonal(months(ist:end),years(ist:end),sstcut2(:,:,ist:end),3:8,3); % was [3 4 5], was 3:8
sst_centered = sst_centered - mean(sst_centered,3);
u10_centered = monthly_to_seasonal(months(ist:end),years(ist:end),u10(:,:,ist:end),3:8,3); % was [3 4 5]
u10_centered = u10_centered - mean(u10_centered,3);
v10_centered = monthly_to_seasonal(months(ist:end),years(ist:end),v10(:,:,ist:end),3:8,3); % was [3 4 5]
v10_centered = v10_centered - mean(v10_centered,3);
slp_centered = monthly_to_seasonal(months(ist:end),years(ist:end),slp(:,:,ist:end),3:8,3); % was [3 4 5]
slp_centered = slp_centered - mean(slp_centered,3);
pr_centered = monthly_to_seasonal(months(ist:end),years(ist:end),pr(:,:,ist:end),3:8,3); % was [3 4 5] %13:end
pr_centered = pr_centered - mean(pr_centered,3);
sst_onepoint_original = sst_onepoint;
sst_onepoint_map = mean(sst_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3)./sqrt(mean(sst_centered.^2,3));
u10_onepoint_map = mean(u10_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3)./sqrt(mean(u10_centered.^2,3));
v10_onepoint_map = mean(v10_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3)./sqrt(mean(v10_centered.^2,3));
slp_onepoint_map = mean(slp_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3)./sqrt(mean(slp_centered.^2,3));
pr_onepoint_map = mean(pr_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[360 180 1]),3)./sqrt(mean(pr_centered.^2,3));
sst_onepointr_map = mean(sst_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3);
u10_onepointr_map = mean(u10_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3);
v10_onepointr_map = mean(v10_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3);
slp_onepointr_map = mean(slp_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3);
pr_onepointr_map = mean(pr_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[360 180 1]),3);

%% Observational SST pattern and mechanism analysis 

sst = sst(:,:,1:end-120);
sst = sst(:,:,end-12*trendlength+1:end);
s = size(sst);
switch season
    case 'Amean'
        sst = squeeze(mean(reshape(sst,[s(1) s(2) 12 trendlength]),3));
    case 'JFMAMJ'
        tmp = reshape(sst,[s(1) s(2) 12 trendlength]);
        sst = squeeze(mean(tmp(:,:,1:6,:),3));
    case 'JASOND'
        tmp = reshape(sst,[s(1) s(2) 12 trendlength]);
        sst = squeeze(mean(tmp(:,:,7:12,:),3));
end
[~,sst_trend] = detrend(sst,3,1);

obs_trend_centered = sst_trend - mean_in_a_box(lon25,lat25,sst_trend,boxlat,boxlon);
%obs_centered = sst - mean(sst,3);
obs_norm = sqrt(mean_in_a_box(lon25,lat25,obs_trend_centered.^2,boxlat,boxlon));
sst_trend_pattern = sst_trend./global_mean(lon25,lat25,sst_trend);

[LFCs, LFPs, r, pvar, weightsf, PCs, EOFs, N, pvar_slow, pvar_lfc, r_eofs, pvar_slow_eofs] = run_lfca(lon25,lat25,sstcut,30,120);

sst_lfp = reshape((LFCs(:,1:2)*LFPs(1:2,:))',[length(lon25) length(lat25) 124*12]);
sst_lfp_trendlength = sst_lfp(:,:,1:end-120);
sst_lfp_trendlength = sst_lfp_trendlength(:,:,end-12*trendlength+1:end);
sst_lfp_trendlength = squeeze(mean(reshape(sst_lfp_trendlength,[s(1) s(2) 12 trendlength]),3));
[~,sst_lfp_trend] = detrend(sst_lfp_trendlength,3,1);
sst_lfp_trend_pattern = sst_lfp_trend./global_mean(lon25,lat25,sst_lfp_trend);

% for observational esimtate of LFvar (after removing 3 LFCs)

ssta = sstcut-sst_lfp;
switch season
    case 'Amean'
        ssta = squeeze(mean(reshape(ssta,[s(1) s(2) 12 124]),3));
    case 'JFMAMJ'
        tmp = reshape(ssta,[s(1) s(2) 12 124]);
        ssta = squeeze(mean(tmp(:,:,1:6,:),3));
    case 'JASOND'
        tmp = reshape(ssta,[s(1) s(2) 12 124]);
        ssta = squeeze(mean(tmp(:,:,7:12,:),3));
end
ssta = ssta - mean(ssta,3);
for l = 1:length(lon25)
    for k = 1:length(lat25)
        sst_lf_obs(l,k,:) = lanczos(squeeze(ssta(l,k,:)),1,10);
    end
end
sst_lf_mse_obs = mean(sst_lf_obs.^2,3);

%% CMIP6 SST pattern and mechanism analysis 

dir = [dataPath,'cmip6/tos/'];

varnam = 'tos'; 

pcorr_all = nan.*ones(length(models),68);
corr_all = nan.*ones(length(models),68);

for i = 1:length(models)
    i
    % find relevant files for each model
    files = get_files(dir);
    files = select_files(files,[models{i},'_']);
    files_tauu = get_files([dataPath,'cmip6/tauu/']);
    files_tauu = select_files(files_tauu,[models{i},'_']);
    files_tauv = get_files([dataPath,'cmip6/tauv/']);
    files_tauv = select_files(files_tauv,[models{i},'_']);
    files_slp = get_files([dataPath,'cmip6/psl/']);
    files_slp = select_files(files_slp,[models{i},'_']);
    files_pr = get_files([dataPath,'cmip6/pr_50to50_seasonal_all/']);
    files_pr = select_files(files_pr,[models{i},'_']);

    % determine member indices
    clear members
    for j = 1:length(files)
        tmp = strsplit(files{j},'/');
        tmp2 = strsplit(tmp{end},'_');
        members{j} = tmp2{5};
    end
    members_all{i} = members;
    sst_all = zeros([length(lon25), length(lat25), 165, length(members)]);
    sst_lf = zeros([length(lon25), length(lat25), 165, 3]);

    for j = 1:length(members)
        field = get_avg_field_nd(files{j},varnam);
        field(abs(field)>1e10) = nan;
        tauu = get_avg_field_nd(files_tauu{j},'tauu');
        tauu(abs(tauu)>1e10) = nan;
        tauv = get_avg_field_nd(files_tauv{j},'tauv');
        tauv(abs(tauv)>1e10) = nan;
        slp = get_avg_field_nd(files_slp{j},'psl');
        slp(abs(slp)>1e10) = nan;
        [pr,lat_cut,lon_cut] = get_avg_field_nd(files_pr{j},{'pr','lat','lon'}); % loading in seasonal pr (versus monthly for other fields)
        pr(pr>1e10) = nan;
        pr = pr.*units;
        if strcmp(models{i},'CIESM')
            pr = pr.*1000;
        elseif strcmp(models{i},'FGOALS-g3')
            pr = pr(:,:,1:end-8); % FGOALS-g3 goes through 2016...
        elseif strcmp(models{i},'NorCPM1')
            pr = pr(:,:,1:end-60);
        end
        %pr_centered = (pr(:,:,6:4:end)+pr(:,:,5:4:end-1))/2; % DJFMAM = (MAM + DJF)/2, skipping first JFMAM and last D
        pr_centered = (pr(:,:,6:4:end)+pr(:,:,7:4:end))/2; % MAMJJA = (MAM + JJA)/2, skipping first year
        pr_centered = pr_centered - mean(pr_centered,3);

        sst_onepoint = mean_in_a_box(lon25,lat25,field,[-20 -15],[258 262]);

        months = repmat(1:12,[1 length(sst_onepoint)/12]); years = floor(1:1/12:length(sst_onepoint)/12+0.99);
        sst_onepoint = monthly_to_seasonal(months,years,sst_onepoint,[9:12 1 2],2);  % was [12 1 2], [9:12 1 2]
        sst_onepoint = sst_onepoint - mean(sst_onepoint);
        sst_onepoint = sst_onepoint./std(sst_onepoint);
        s = size(field);
        %sst_centered = monthly_anomalies(field(:,:,13:end));
        ist = 13;
        sst_centered = monthly_to_seasonal(months(ist:end),years(ist:end),field(:,:,ist:end),3:8,3); % was [3 4 5], 3:8
        sst_centered = sst_centered - mean(sst_centered,3);
        if size(field,3)<1980
            if size(field,3)==1380 % IITM-ESM exception
                tauu = tauu(:,:,end-1380+1:end);
                tauv = tauv(:,:,end-1380+1:end);
                slp = slp(:,:,end-1380+1:end);
                pr_centered = pr_centered(:,:,end-114+1:end); % pr is already seasonal, instead of monthly
                pr_centered = pr_centered - mean(pr_centered,3);
            else
                disp('Warning: unexpected size of field')
            end
        end
        tauu_centered = monthly_to_seasonal(months(ist:end),years(ist:end),tauu(:,:,ist:end),3:8,3); % was [3 4 5], 3:8
        tauu_centered = tauu_centered - mean(tauu_centered,3);
        tauv_centered = monthly_to_seasonal(months(ist:end),years(ist:end),tauv(:,:,ist:end),3:8,3); % was [3 4 5], 3:8
        tauv_centered = tauv_centered - mean(tauv_centered,3);
        slp_centered = monthly_to_seasonal(months(ist:end),years(ist:end),slp(:,:,ist:end),3:8,3); % was [3 4 5], 3:8
        slp_centered = slp_centered - mean(slp_centered,3);
        %sst_onepoint = sst_onepoint(1:end-1);
        sst_onepoint_map_model(:,:,i,j) = mean(sst_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3)./sqrt(mean(sst_centered.^2,3));
        tauu_onepoint_map_model(:,:,i,j) = mean(tauu_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3)./sqrt(mean(tauu_centered.^2,3));
        tauv_onepoint_map_model(:,:,i,j) = mean(tauv_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3)./sqrt(mean(tauv_centered.^2,3));
        slp_onepoint_map_model(:,:,i,j) = mean(slp_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3)./sqrt(mean(slp_centered.^2,3));
        pr_onepoint_map_model(:,:,i,j) = mean(pr_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[size(pr,1) size(pr,2) 1]),3)./sqrt(mean(pr_centered.^2,3));
        sst_onepointr_map_model(:,:,i,j) = mean(sst_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3);
        tauu_onepointr_map_model(:,:,i,j) = mean(tauu_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3);
        tauv_onepointr_map_model(:,:,i,j) = mean(tauv_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3);
        slp_onepointr_map_model(:,:,i,j) = mean(slp_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[s(1) s(2) 1]),3);
        pr_onepointr_map_model(:,:,i,j) = mean(pr_centered.*repmat(reshape(sst_onepoint,[1 1 length(sst_onepoint)]),[size(pr,1) size(pr,2) 1]),3);

        tos_clim_model(:,:,:,i,j) = mean(reshape(field(:,:,1:50*12),[s(1) s(2) 12 50]),4);

        % Computing SST trend in chosen season
        s = size(field);
        switch season
            case 'Amean'
                field = squeeze(mean(reshape(field,[s(1) s(2) 12 s(3)/12]),3));
            case 'JFMAMJ'
                tmp = reshape(field,[s(1) s(2) 12 s(3)/12]);
                field = squeeze(mean(tmp(:,:,1:6,:),3));
            case 'JASOND'
                tmp = reshape(field,[s(1) s(2) 12 s(3)/12]);
                field = squeeze(mean(tmp(:,:,7:12,:),3));
        end
        field_trendlength = field(:,:,end-trendlength+1:end);
        [~,tos_trend] = detrend(field_trendlength,3,1);
        tos_trend_model(:,:,i,j) = tos_trend;
        model_trend_centered = tos_trend - mean_in_a_box(lon25,lat25,tos_trend,boxlat,boxlon);
        model_norm = sqrt(mean_in_a_box(lon25,lat25,model_trend_centered.^2,boxlat,boxlon));
        pcorr_all(i,j) = mean_in_a_box(lon25,lat25,obs_trend_centered.*model_trend_centered,boxlat,boxlon)./(obs_norm.*model_norm);
        %model_centered = field_trendlength - mean(field_trendlength,3);
        %corr_all(i,j) = mean_in_a_box(lon25,lat25,mean(obs_centered.*model_centered,3),boxlat,boxlon)./sqrt(mean_in_a_box(lon25,lat25,mean(obs_centered.^2,3),boxlat,boxlon).*mean_in_a_box(lon25,lat25,mean(model_centered.^2,3),boxlat,boxlon));
        tos_trend_pattern_model(:,:,i,j) = tos_trend_model(:,:,i,j)./global_mean(lon25,lat25,tos_trend_model(:,:,i,j));

        % outputing SST for all models and members for use later in code
        try
            sst_all(:,:,:,j) = field;
        catch
            sst_all(:,:,end-115+1:end,j) = field;
        end

    end
    
    % low-frequency SST variance is computed from 3 members in all models,
    % and models with less than 3 members are excluded
    if length(members)>2

        clear ssta

        for j = 1:3
            ssta = sst_all(:,:,:,j) - mean(sst_all,4);
            for l = 1:length(lon25)
                for k = 1:length(lat25)
                    sst_lf(l,k,:,j) = lanczos(squeeze(ssta(l,k,:)),1,10);
                end
            end
        end
    end
    sst_lf_mse(:,:,i) = mean(mean(sst_lf.^2,3),4);

    % fill in output data with NaNs for member index higher than the number
    % of members in each model
    tos_clim_model(:,:,:,i,length(members)+1:68) = nan;
    tos_trend_model(:,:,i,length(members)+1:68) = nan;
    tos_trend_pattern_model(:,:,i,length(members)+1:68) = nan;
    sst_onepoint_map_model(:,:,i,length(members)+1:68) = nan;
    tauu_onepoint_map_model(:,:,i,length(members)+1:68) = nan;
    tauv_onepoint_map_model(:,:,i,length(members)+1:68) = nan;
    slp_onepoint_map_model(:,:,i,length(members)+1:68) = nan;
    pr_onepoint_map_model(:,:,i,length(members)+1:68) = nan;
    sst_onepointr_map_model(:,:,i,length(members)+1:68) = nan;
    tauu_onepointr_map_model(:,:,i,length(members)+1:68) = nan;
    tauv_onepointr_map_model(:,:,i,length(members)+1:68) = nan;
    slp_onepointr_map_model(:,:,i,length(members)+1:68) = nan;
    pr_onepointr_map_model(:,:,i,length(members)+1:68) = nan;
end

precip_index_plot = repmat(precip_index,[1 68]);

%%

% % SST-sharp calculation (legacy code, not used)
% 
% tos_amean_clim_model = squeeze(mean(nanmean(tos_clim_model,5),3));
% [Y,X] = meshgrid(lat25,lon25);
% 
% for i = 1
%     mask = nan.*ones([length(lon25), length(lat25)]);
%     
%     tmp = squeeze(tos_amean_clim_model(:,31:42,j));
%     tclim_prctile = prctile(tmp(:),30);
%     tmp = squeeze(tos_amean_clim_model(:,:,j));
%     mask(abs(Y)<15 & tmp<tclim_prctile) = 1;
% 
%     SSTantisharp(i,:) = global_mean(lon25,lat25,squeeze(tos_trend_pattern_model(:,:,i,:)).*repmat(mask,[1 1 68]));
% end

%%

% % Alternative SST trend indeces (legacy code, not used)

%figure; plot(precip_index_plot(:),pcorr_all(:),'k.')

% mask = ones([length(lon25), length(lat25)]);
% mask(lon25<170,:) = nan;
% mask(lon25>240,:) = nan;
% [Y,X] = meshgrid(lat25,lon25);
% mask(Y>(X-170).*(3/7)) = nan;
% mask(-Y>(X-170).*(3/7)) = nan;
% sst_triangle = global_mean(lon25,lat25,tos_trend_model.*mask);

%figure; plot(precip_index_plot(:),sst_triangle(:),'k.')

% mask = ones([length(lon25), length(lat25)]);
% mask(lon25<160,:) = nan;
% mask(lon25>270,:) = nan;
% [Y,X] = meshgrid(lat25,lon25);
% mask(Y>0 & X<200) = nan;
% mask(-Y>(X-160).*(3/11)) = nan;
% mask(-Y<(X-200).*(3/11)) = nan;
% sst_diagonal = global_mean(lon25,lat25,tos_trend_model.*mask);

% mask = ones([length(lon25), length(lat25)]);
% mask(lon25<200,:) = nan;
% mask(lon25>270,:) = nan;
% [Y,X] = meshgrid(lat25,lon25);
% mask(Y>0 & X<250) = nan;
% mask(-Y>(X-200).*(3/7)) = nan;
% mask(-Y<(X-250).*(3/7)) = nan;
% sst_diagonal_obs = global_mean(lon25,lat25,sst_trend.*mask);

% figure; plot(precip_index_plot(:),sst_diagonal(:),'k.')
% hold on; plot([3 13],sst_diagonal_obs.*[1 1],'k-')

%% Equatorial Warming Minimum (Fig. 1e)

for i = 1:68
    for j = 1:68
        eq_warming_min(i,j) = latmean_lonmin_in_a_box(lon25,lat25,tos_trend_pattern_model(:,:,i,j),[-2.5 2.5],[150 240]);
    end
end
for i = 1:68
    if length(members_all{i})>2
        eq_warming_min_ens(i) = nanmedian(eq_warming_min(i,:));
    else
        eq_warming_min_ens(i) = nan;
    end
end

eq_warming_min_obs = latmean_lonmin_in_a_box(lon25,lat25,sst_trend_pattern,[-3 3],[150 240]);

figure; plot(eq_warming_min(:),precip_index_plot(:),'.','color',[0.5 0.5 0.5],'markersize',9)
hold on; plot(eq_warming_min_ens,precip_index,'ko','markerfacecolor','k')
hold on; plot(eq_warming_min_obs.*[1 1],[0 14],'k-')
set(gca,'ylim',[0 13]); set(gca,'ytick',0:2:14)
set(gca,'xlim',[-2.5 2]); set(gca,'xtick',-2:1:2)
hold on; plot([-3 3],[precip_imerg_DJFMAM_SHmax_lmax precip_imerg_DJFMAM_SHmax_lmax],'k--')
hold on; plot([-3 3],[precip_gpcp_DJFMAM_SHmax_lmax precip_gpcp_DJFMAM_SHmax_lmax],'k--')
pretty_figure(400,300,'Equatorial Pacific Warming Minimum (K K^{-1})','Max DJFMAM Precip. (mm day^{-1})','none','none',16);

%% SST Trend Maps (Fig. 2)

plot_field_div(lon25,lat25,sst_trend_pattern,linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45]) % originally [-40 40]
hold on; contour(lon,lat,precip_imerg_DJFMAM',[4 4],'k','linewidth',1)
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

% & rmse_imerg<rmse_threshold
plot_field_div(lon25,lat25,mean(nanmean(tos_trend_pattern_model(:,:,precip_index>threshold,:),4),3),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
hold on; contour(lon,lat,squeeze(nanmean(precip_model_DJFMAM(:,:,:,precip_index>threshold),4))',[4 4],'k','linewidth',1)
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

plot_field_div(lon25,lat25,mean(nanmean(tos_trend_pattern_model(:,:,precip_index<threshold,:),4),3),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
hold on; contour(lon,lat,squeeze(nanmean(precip_model_DJFMAM(:,:,:,precip_index<threshold),4))',[4 4],'k','linewidth',1)
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
% t = (mean(nanmean(tos_trend_pattern_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(tos_trend_pattern_model(:,:,precip_index>threshold,:),4),3))./sqrt(var(nanmean(tos_trend_pattern_model(:,:,precip_index<threshold,:),4),0,3)./sum(precip_index<threshold) + var(nanmean(tos_trend_pattern_model(:,:,precip_index>threshold,:),4),0,3)./sum(precip_index>threshold));
% plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.96,7)
% plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.645,3)

plot_field_div(lon25,lat25,mean(nanmean(tos_trend_pattern_model(:,:,precip_index<threshold,:),4),3)-nanmean(nanmean(tos_trend_pattern_model(:,:,precip_index>threshold,:),3),4),linspace(-0.8,0.8,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
t = (mean(nanmean(tos_trend_pattern_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(tos_trend_pattern_model(:,:,precip_index>threshold,:),4),3))./sqrt(var(nanmean(tos_trend_pattern_model(:,:,precip_index<threshold,:),4),0,3)./sum(precip_index<threshold) + var(nanmean(tos_trend_pattern_model(:,:,precip_index>threshold,:),4),0,3)./sum(precip_index>threshold));
plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.96,7)
plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.645,3)

%% Decadal Variance Maps (Fig. 3)

plot_field_div(lon25,lat25,sqrt(sst_lf_mse_obs),linspace(0,0.5,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

plot_field_div(lon25,lat25,sqrt(mean(sst_lf_mse(:,:,precip_index>threshold),3)),linspace(0,0.5,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

plot_field_div(lon25,lat25,sqrt(mean(sst_lf_mse(:,:,precip_index<threshold),3)),linspace(0,0.5,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
% t = (mean(sst_lf_mse(:,:,precip_index<threshold),3)-mean(sst_lf_mse(:,:,precip_index>threshold),3))./sqrt(var(sst_lf_mse(:,:,precip_index<threshold),0,3)./sum(precip_index<threshold) + var(sst_lf_mse(:,:,precip_index>threshold),0,3)./sum(precip_index>threshold));
% plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.96,7)

plot_field_div(lon25,lat25,sqrt(mean(sst_lf_mse(:,:,precip_index<threshold),3))-sqrt(mean(sst_lf_mse(:,:,precip_index>threshold),3)),linspace(-0.2,0.2,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
t = (mean(sst_lf_mse(:,:,precip_index<threshold),3)-mean(sst_lf_mse(:,:,precip_index>threshold),3))./sqrt(var(sst_lf_mse(:,:,precip_index<threshold),0,3)./sum(precip_index<threshold) + var(sst_lf_mse(:,:,precip_index>threshold),0,3)./sum(precip_index>threshold));
plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.96,7)

%% One-point Correlation Maps (Fig. 4)

qs = 3; % quiver spacing

plot_field_div(lon25,lat25,sst_onepoint_map,linspace(-1,1,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
[X,Y] = meshgrid(lon25(qs:qs:end),lat25(qs:qs:end));
hold on; ncquiverref(X,Y,u10_onepoint_map(qs:qs:end,qs:qs:end)',v10_onepoint_map(qs:qs:end,qs:qs:end)','',0.5,'true','k');
hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
% gfilt = fspecial('gaussian',6,1);
% field = conv2(pr_onepoint_map,gfilt,'same');
% hold on; contour(lon,lat,field',0.25:0.25:1,'color',[6 64 43]/256,'linewidth',1)
% hold on; contour(lon,lat,-field',0.25:0.25:1,'color',[82 32 129]/256,'linewidth',1)

plot_field_div(lon25,lat25,mean(nanmean(sst_onepoint_map_model(:,:,precip_index>threshold,:),4),3),linspace(-1,1,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
[X,Y] = meshgrid(lon25(qs:qs:end),lat25(qs:qs:end));
hold on; ncquiverref(X,Y,mean(nanmean(tauu_onepoint_map_model(qs:qs:end,qs:qs:end,precip_index>threshold,:),4),3)',mean(nanmean(tauv_onepoint_map_model(qs:qs:end,qs:qs:end,precip_index>threshold,:),4),3)','',0.5,'true','k');
hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
% hold on; contour(lon,lat,mean(nanmean(pr_onepoint_map_model(:,:,precip_index>threshold,:),4),3)',0.25:0.25:1,'color',[6 64 43]/256,'linewidth',1)
% hold on; contour(lon,lat,-mean(nanmean(pr_onepoint_map_model(:,:,precip_index>threshold,:),4),3)',0.25:0.25:1,'color',[82 32 129]/256,'linewidth',1)

plot_field_div(lon25,lat25,mean(nanmean(sst_onepoint_map_model(:,:,precip_index<threshold,:),4),3),linspace(-1,1,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
[X,Y] = meshgrid(lon25(qs:qs:end),lat25(qs:qs:end));
hold on; ncquiverref(X,Y,mean(nanmean(tauu_onepoint_map_model(qs:qs:end,qs:qs:end,precip_index<threshold,:),4),3)',mean(nanmean(tauv_onepoint_map_model(qs:qs:end,qs:qs:end,precip_index<threshold,:),4),3)','',0.5,'true','k');
hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
% hold on; contour(lon,lat,mean(nanmean(pr_onepoint_map_model(:,:,precip_index<threshold,:),4),3)',0.25:0.25:1,'color',[6 64 43]/256,'linewidth',1)
% hold on; contour(lon,lat,-mean(nanmean(pr_onepoint_map_model(:,:,precip_index<threshold,:),4),3)',0.25:0.25:1,'color',[82 32 129]/256,'linewidth',1)

plot_field_div(lon25,lat25,mean(nanmean(sst_onepoint_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(sst_onepoint_map_model(:,:,precip_index>threshold,:),4),3),linspace(-0.4,0.4,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
t = (mean(nanmean(sst_onepoint_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(sst_onepoint_map_model(:,:,precip_index>threshold,:),4),3))./sqrt(var(nanmean(sst_onepoint_map_model(:,:,precip_index<threshold,:),4),0,3)./sum(precip_index<threshold) + var(nanmean(sst_onepoint_map_model(:,:,precip_index<threshold,:),4),0,3)./sum(precip_index>threshold));
plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.96,7)
plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.645,3)
hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
% hold on; contour(lon,lat,(mean(nanmean(pr_onepoint_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(pr_onepoint_map_model(:,:,precip_index>threshold,:),4),3))',0.1:0.1:1,'color',[6 64 43]/256,'linewidth',1)
% hold on; contour(lon,lat,-(mean(nanmean(pr_onepoint_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(pr_onepoint_map_model(:,:,precip_index>threshold,:),4),3))',0.1:0.1:1,'color',[82 32 129]/256,'linewidth',1)

% %% One-point Regression Maps

% qs = 3; % quiver spacing
% 
% plot_field_div(lon25,lat25,sst_onepointr_map,linspace(-0.6,0.6,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
% pretty_figure(800,175,'none','none','none','none',16);
% set(gca,'color',[0.95 0.95 0.95])
% [X,Y] = meshgrid(lon25(qs:qs:end),lat25(qs:qs:end));
% hold on; ncquiverref(X,Y,u10_onepoint_map(qs:qs:end,qs:qs:end)',v10_onepoint_map(qs:qs:end,qs:qs:end)','',0.5,'true','k');
% hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
% hold on; contour(lon,lat,pr_onepointr_map',0.4:0.4:2,'color',[6 64 43]/256,'linewidth',1)
% hold on; contour(lon,lat,-pr_onepointr_map',0.4:0.4:2,'color',[82 32 129]/256,'linewidth',1)
% 
% plot_field_div(lon25,lat25,mean(nanmean(sst_onepointr_map_model(:,:,precip_index>threshold,:),4),3),linspace(-0.6,0.6,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
% pretty_figure(800,175,'none','none','none','none',16);
% set(gca,'color',[0.95 0.95 0.95])
% [X,Y] = meshgrid(lon25(qs:qs:end),lat25(qs:qs:end));
% hold on; ncquiverref(X,Y,mean(nanmean(tauu_onepoint_map_model(qs:qs:end,qs:qs:end,precip_index>threshold,:),4),3)',mean(nanmean(tauv_onepoint_map_model(qs:qs:end,qs:qs:end,precip_index>threshold,:),4),3)','',0.5,'true','k');
% hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
% hold on; contour(lon,lat,mean(nanmean(pr_onepointr_map_model(:,:,precip_index>threshold,:),4),3)',0.4:0.4:2,'color',[6 64 43]/256,'linewidth',1)
% hold on; contour(lon,lat,-mean(nanmean(pr_onepointr_map_model(:,:,precip_index>threshold,:),4),3)',0.4:0.4:2,'color',[82 32 129]/256,'linewidth',1)
% 
% plot_field_div(lon25,lat25,mean(nanmean(sst_onepointr_map_model(:,:,precip_index<threshold,:),4),3),linspace(-0.6,0.6,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
% pretty_figure(800,175,'none','none','none','none',16);
% set(gca,'color',[0.95 0.95 0.95])
% [X,Y] = meshgrid(lon25(qs:qs:end),lat25(qs:qs:end));
% hold on; ncquiverref(X,Y,mean(nanmean(tauu_onepoint_map_model(qs:qs:end,qs:qs:end,precip_index<threshold,:),4),3)',mean(nanmean(tauv_onepoint_map_model(qs:qs:end,qs:qs:end,precip_index<threshold,:),4),3)','',0.5,'true','k');
% hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
% hold on; contour(lon,lat,mean(nanmean(pr_onepointr_map_model(:,:,precip_index<threshold,:),4),3)',0.4:0.4:2,'color',[6 64 43]/256,'linewidth',1)
% hold on; contour(lon,lat,-mean(nanmean(pr_onepointr_map_model(:,:,precip_index<threshold,:),4),3)',0.4:0.4:2,'color',[82 32 129]/256,'linewidth',1)
% 
% plot_field_div(lon25,lat25,mean(nanmean(sst_onepointr_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(sst_onepointr_map_model(:,:,precip_index>threshold,:),4),3),linspace(-0.3,0.3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
% pretty_figure(800,175,'none','none','none','none',16);
% set(gca,'color',[0.95 0.95 0.95])
% t = (mean(nanmean(sst_onepointr_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(sst_onepointr_map_model(:,:,precip_index>threshold,:),4),3))./sqrt(var(nanmean(sst_onepointr_map_model(:,:,precip_index<threshold,:),4),0,3)./sum(precip_index<threshold) + var(nanmean(sst_onepointr_map_model(:,:,precip_index<threshold,:),4),0,3)./sum(precip_index>threshold));
% plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.96,7)
% plot_stipple(lon25(2:2:end),lat25(2:2:end),abs(t(2:2:end,2:2:end))>1.645,3)
% hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
% hold on; contour(lon,lat,(mean(nanmean(pr_onepointr_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(pr_onepointr_map_model(:,:,precip_index>threshold,:),4),3))',0.2:0.2:2,'color',[6 64 43]/256,'linewidth',1)
% hold on; contour(lon,lat,-(mean(nanmean(pr_onepointr_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(pr_onepointr_map_model(:,:,precip_index>threshold,:),4),3))',0.2:0.2:2,'color',[82 32 129]/256,'linewidth',1)

%% One-point Regression Maps PR/SLP (Extended Data Figure 4)

plot_field_robinson(lon25,lat25,slp_onepointr_map,linspace(-120,120,25));
pretty_figure(800,350,'none','none','none','none',16);
gfilt = fspecial('gaussian',6,0.5);
field = conv2(pr_onepointr_map,gfilt,'same');
hold on; contourm(lat,lon,field',0.1667:0.333:4,'k','linewidth',1)
hold on; contourm(lat,lon,-field',0.1667:0.333:4,'k--','linewidth',1)
colorbar('southoutside')

plot_field_robinson(lon25,lat25,mean(nanmean(slp_onepointr_map_model(:,:,precip_index>threshold,:),4),3),linspace(-120,120,25));
pretty_figure(800,350,'none','none','none','none',16);
hold on; contourm(lat,lon,mean(nanmean(pr_onepointr_map_model(:,:,precip_index>threshold,:),4),3)',0.1667:0.333:4,'k','linewidth',1)
hold on; contourm(lat,lon,-mean(nanmean(pr_onepointr_map_model(:,:,precip_index>threshold,:),4),3)',0.1667:0.333:4,'k--','linewidth',1)
colorbar('southoutside')

plot_field_robinson(lon25,lat25,mean(nanmean(slp_onepointr_map_model(:,:,precip_index<threshold,:),4),3),linspace(-120,120,25)); 
pretty_figure(800,350,'none','none','none','none',16);
hold on; contourm(lat,lon,mean(nanmean(pr_onepointr_map_model(:,:,precip_index<threshold,:),4),3)',0.1667:0.333:4,'k','linewidth',1)
hold on; contourm(lat,lon,-mean(nanmean(pr_onepointr_map_model(:,:,precip_index<threshold,:),4),3)',0.1667:0.333:4,'k--','linewidth',1)
colorbar('southoutside')

plot_field_robinson(lon25,lat25,mean(nanmean(slp_onepointr_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(slp_onepointr_map_model(:,:,precip_index>threshold,:),4),3),linspace(-40,40,25)); 
pretty_figure(800,350,'none','none','none','none',16);
hold on; contourm(lat,lon,(mean(nanmean(pr_onepointr_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(pr_onepointr_map_model(:,:,precip_index>threshold,:),4),3))',0.0833:0.1667:2,'k','linewidth',1)
hold on; contourm(lat,lon,-(mean(nanmean(pr_onepointr_map_model(:,:,precip_index<threshold,:),4),3)-mean(nanmean(pr_onepointr_map_model(:,:,precip_index>threshold,:),4),3))',0.0833:0.1667:2,'k--','linewidth',1)
colorbar('southoutside')

%% SST Trends for specific models (Extended Data Figures 2, 6)

is = find(precip_index<threshold);

for i = 1:length(is)
    if is(i) == 45 || is(i) == 3
        ICON = mean(nanmean(tos_trend_pattern_model(:,:,is(i),:),4),3); % ICON or AWI
        CESM2 = mean(nanmean(tos_trend_pattern_model(:,:,12,:),4),3);
        ICON(isnan(CESM2)) = nan;
        plot_field_div(lon25,lat25,ICON,linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
        hold on; contour(lon,lat,squeeze(precip_model_DJFMAM(:,:,:,is(i)))',[4 4],'k','linewidth',1)
        pretty_figure(800,175,'none','none','none','none',16);
        set(gca,'color',[0.95 0.95 0.95])
    else
        plot_field_div(lon25,lat25,mean(nanmean(tos_trend_pattern_model(:,:,is(i),:),4),3),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
        hold on; contour(lon,lat,squeeze(precip_model_DJFMAM(:,:,:,is(i)))',[4 4],'k','linewidth',1)
        pretty_figure(800,175,'none','none','none','none',16);
        set(gca,'color',[0.95 0.95 0.95])
    end
end

plot_field_div(lon25,lat25,mean(nanmean(tos_trend_pattern_model(:,:,12,:),4),3),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
hold on; contour(lon,lat,squeeze(precip_model_DJFMAM(:,:,:,12))',[4 4],'k','linewidth',1)
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

plot_field_div(lon25,lat25,mean(nanmean(tos_trend_pattern_model(:,:,9,:),4),3),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
hold on; contour(lon,lat,squeeze(precip_model_DJFMAM(:,:,:,9))',[4 4],'k','linewidth',1)
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

plot_field_div(lon25,lat25,mean(nanmean(tos_trend_pattern_model(:,:,10,:),4),3),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
hold on; contour(lon,lat,squeeze(precip_model_DJFMAM(:,:,:,10))',[4 4],'k','linewidth',1)
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

plot_field_div(lon25,lat25,mean(nanmean(tos_trend_pattern_model(:,:,63,:),4),3),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
hold on; contour(lon,lat,squeeze(precip_model_DJFMAM(:,:,:,10))',[4 4],'k','linewidth',1)
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

%% Individual ensemble members for CESM2-WACCM and NorESM2-MM (Extended Data Figure 3)

for i = 1:3
    plot_field_div(lon25,lat25,tos_trend_pattern_model(:,:,11,i),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
    hold on; contour(lon,lat,squeeze(precip_model_DJFMAM(:,:,:,11))',[4 4],'k','linewidth',1)
    pretty_figure(800,175,'none','none','none','none',16);
    set(gca,'color',[0.95 0.95 0.95])
end

for i = 1:3
    plot_field_div(lon25,lat25,tos_trend_pattern_model(:,:,64,i),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
    hold on; contour(lon,lat,squeeze(precip_model_DJFMAM(:,:,:,64))',[4 4],'k','linewidth',1)
    pretty_figure(800,175,'none','none','none','none',16);
    set(gca,'color',[0.95 0.95 0.95])
end

%% Decadal variance for specific models (Extended Data Figures 5, 7)

is = find(precip_index<threshold);

for i = 1:length(is)
    if is(i) == 45 || is(i) == 3
        ICON = sqrt(mean(sst_lf_mse(:,:,is(i)),3)); % ICON or AWI
        CESM2 = sqrt(mean(sst_lf_mse(:,:,12),3));
        ICON(isnan(CESM2)) = nan;
        plot_field_div(lon25,lat25,ICON,linspace(0,0.5,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
        pretty_figure(800,175,'none','none','none','none',16);
        set(gca,'color',[0.95 0.95 0.95])
        set(gcf,'renderer','painters')
    else
        plot_field_div(lon25,lat25,sqrt(mean(sst_lf_mse(:,:,is(i)),3)),linspace(0,0.5,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
        pretty_figure(800,175,'none','none','none','none',16);
        set(gca,'color',[0.95 0.95 0.95])
        set(gcf,'renderer','painters')
    end
end

plot_field_div(lon25,lat25,sqrt(mean(sst_lf_mse(:,:,12),3)),linspace(0,0.5,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
set(gcf,'renderer','painters')

plot_field_div(lon25,lat25,sqrt(mean(sst_lf_mse(:,:,9),3)),linspace(0,0.5,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
set(gcf,'renderer','painters')

plot_field_div(lon25,lat25,sqrt(mean(sst_lf_mse(:,:,10),3)),linspace(0,0.5,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
set(gcf,'renderer','painters')

plot_field_div(lon25,lat25,sqrt(mean(sst_lf_mse(:,:,63),3)),linspace(0,0.5,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
set(gcf,'renderer','painters')

%% One-point Correlation maps for specific models

is = find(precip_index<threshold);

for i = 1:length(is)
    if is(i) == 45 || is(i) == 3
        ICON = nanmean(sst_onepoint_map_model(:,:,is(i),:),4); % ICON or AWI
        CESM2 = nanmean(sst_onepoint_map_model(:,:,12,:),4);
        ICON(isnan(CESM2)) = nan;
        plot_field_div(lon25,lat25,ICON,linspace(-1,1,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
        pretty_figure(800,175,'none','none','none','none',16);
        set(gca,'color',[0.95 0.95 0.95])
        [X,Y] = meshgrid(lon25(2:2:end),lat25(2:2:end));
        hold on; ncquiverref(X,Y,nanmean(tauu_onepoint_map_model(2:2:end,2:2:end,is(i),:),4)',nanmean(tauv_onepoint_map_model(2:2:end,2:2:end,is(i),:),4)','',0.5,'true','k');
        hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
    else
        plot_field_div(lon25,lat25,nanmean(sst_onepoint_map_model(:,:,is(i),:),4),linspace(-1,1,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
        pretty_figure(800,175,'none','none','none','none',16);
        set(gca,'color',[0.95 0.95 0.95])
        [X,Y] = meshgrid(lon25(2:2:end),lat25(2:2:end));
        hold on; ncquiverref(X,Y,nanmean(tauu_onepoint_map_model(2:2:end,2:2:end,is(i),:),4)',nanmean(tauv_onepoint_map_model(2:2:end,2:2:end,is(i),:),4)','',0.5,'true','k');
        hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
    end
end

for i = [12 9 10 63]
    plot_field_div(lon25,lat25,nanmean(sst_onepoint_map_model(:,:,i,:),4),linspace(-1,1,21)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45])
    pretty_figure(800,175,'none','none','none','none',16);
    set(gca,'color',[0.95 0.95 0.95])
    [X,Y] = meshgrid(lon25(2:2:end),lat25(2:2:end));
    hold on; ncquiverref(X,Y,nanmean(tauu_onepoint_map_model(2:2:end,2:2:end,i,:),4)',nanmean(tauv_onepoint_map_model(2:2:end,2:2:end,i,:),4)','',0.5,'true','k');
    hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
end

%% Panels for schematic (Fig. 5)

plot_field_div(lon25,lat25,mean(nanmean(sst_onepointr_map_model(:,:,precip_index>threshold,:),4),3),linspace(-0.6,0.6,25)); coast_vector; set(gca,'xlim',[100 300]); set(gca,'ylim',[-60 45])
%camproj('perspective')
%view(12,65)
set(gca,'box','off')
set(gca,'color',[0.95 0.95 0.95])
hold on; contour(lon,lat,nanmean(precip_model_DJFMAM(:,:,:,precip_index>threshold),4)',[4 4],'k','linewidth',1)
hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
set(gca,'xtick','')
set(gca,'ytick','')

plot_field_div(lon25,lat25,mean(nanmean(sst_onepointr_map_model(:,:,precip_index<threshold,:),4),3),linspace(-0.6,0.6,25)); coast_vector; set(gca,'xlim',[100 300]); set(gca,'ylim',[-60 45])
%camproj('perspective')
%view(12,65)
set(gca,'box','off')
set(gca,'color',[0.95 0.95 0.95])
hold on; contour(lon,lat,nanmean(precip_model_DJFMAM(:,:,:,precip_index<threshold),4)',[4 4],'k','linewidth',1)
hold on; plot([260 260],[-17.5 -17.5],'ko','markersize',8','markerfacecolor','k')
set(gca,'xtick','')
set(gca,'ytick','')

%% Extended Data Figure 8 (separate calculation)

ddir = '/Users/rjnglin/Data/ITCZ_SST_trend_biases/cmip6/tos';

trendlength = 40;
tlmo = 12*trendlength;

files = get_files(ddir);
files_M = select_files(files,'MIROC6');
files_E = select_files(files,'E3SM-2-0');
tos_M = get_avg_field_nd(files_M,'tos');
tos_E = get_avg_field_nd(files_E,'tos');
tos_M(abs(tos_M)>1000) = nan;
tos_E(abs(tos_E)>1000) = nan;
tos_M_amean = squeeze(mean(reshape(tos_M(:,:,end-tlmo+1:end),[144 72 12 trendlength]),3));
tos_E_amean = squeeze(mean(reshape(tos_E(:,:,end-tlmo+1:end),[144 72 12 trendlength]),3));
[~,tos_M_trend_hist] = detrend(tos_M_amean,3,1);
[~,tos_E_trend_hist] = detrend(tos_E_amean,3,1);

ddir = '/Users/rjnglin/Data/ITCZ_SST_trend_biases/cmip6/tos_hist-ghg';

files = get_files(ddir);
files_M = select_files(files,'MIROC6');
files_E = select_files(files,'E3SM-2-0');
tos_M = get_avg_field_nd(files_M,'tos');
tos_E = get_avg_field_nd(files_E,'tos');
tos_M(abs(tos_M)>1000) = nan;
tos_E(abs(tos_E)>1000) = nan;
tos_M_amean = squeeze(mean(reshape(tos_M(:,:,end-tlmo+1:end),[144 72 12 trendlength]),3));
tos_E_amean = squeeze(mean(reshape(tos_E(:,:,end-tlmo+1:end),[144 72 12 trendlength]),3));
[~,tos_M_trend_hist_ghg] = detrend(tos_M_amean,3,1);
[~,tos_E_trend_hist_ghg] = detrend(tos_E_amean,3,1);

ddir = '/Users/rjnglin/Data/ITCZ_SST_trend_biases/cmip6/tos_hist-aer';

files = get_files(ddir);
files_M = select_files(files,'MIROC6');
files_E = select_files(files,'E3SM-2-0');
tos_M = get_avg_field_nd(files_M,'tos');
tos_E = get_avg_field_nd(files_E,'tos');
tos_M(abs(tos_M)>1000) = nan;
tos_E(abs(tos_E)>1000) = nan;
tos_M_amean = squeeze(mean(reshape(tos_M(:,:,end-tlmo+1:end),[144 72 12 trendlength]),3));
tos_E_amean = squeeze(mean(reshape(tos_E(:,:,end-tlmo+1:end),[144 72 12 trendlength]),3));
[~,tos_M_trend_hist_aer] = detrend(tos_M_amean,3,1);
[~,tos_E_trend_hist_aer] = detrend(tos_E_amean,3,1);

%%

plot_field_div(lon25,lat25,tos_M_trend_hist./global_mean(lon25,lat25,tos_M_trend_hist),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45]) % originally [-40 40]
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
plot_field_div(lon25,lat25,tos_E_trend_hist./global_mean(lon25,lat25,tos_E_trend_hist),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45]) % originally [-40 40]
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

plot_field_div(lon25,lat25,tos_M_trend_hist_aer./global_mean(lon25,lat25,tos_M_trend_hist),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45]) % originally [-40 40]
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
plot_field_div(lon25,lat25,tos_E_trend_hist_aer./global_mean(lon25,lat25,tos_E_trend_hist),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45]) % originally [-40 40]
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])

plot_field_div(lon25,lat25,(tos_M_trend_hist-tos_M_trend_hist_aer)./global_mean(lon25,lat25,tos_M_trend_hist),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45]) % originally [-40 40]
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
plot_field_div(lon25,lat25,(tos_E_trend_hist-tos_E_trend_hist_aer)./global_mean(lon25,lat25,tos_E_trend_hist),linspace(-3,3,25)); coast_vector; set(gca,'xlim',[0 360]); set(gca,'ylim',[-45 45]) % originally [-40 40]
pretty_figure(800,175,'none','none','none','none',16);
set(gca,'color',[0.95 0.95 0.95])
