%    This Code is for L&OL article "Transport and environmental impact of ash 
%induced by the Hunga Tonga-Hunga Ha'apai volcanic eruption"
%    Authors:Yuntao Wang, Huaguo Zhang, Chang Liu, Mei Zheng, Huan-Huan Chen, 
%Wenting Cao, Shuangling Chen, Chen Wang, Xuehua Fan, Zhihua Mao, Fei Chai
%--------------------------------------------------------------------------
% Fig.4----Chlorophyll
clear;
load("Chlorophyll.mat")

xlim_Reg = [70 200];ylim_Reg = [-28 -8];
xtam=diff(xlim_Reg)*cos(mean(ylim_Reg)*pi/180);
ytam=diff(ylim_Reg);

dif = Chl-Chl_cli;
td = nan(721,3360);
tm = td;
for i = 1:721
    for j = 1:3360
        td(i,j) = mean(dif(i,j,14:41),'omitnan');
        tm(i,j) = mean(Chl_cli(i,j,14:41),'omitnan');
    end
end
subplot('position',[0.1 0.5 0.8 0.19]);
pcolor(lon_Chl,lat_Chl,td);shading flat;
hold on
contour(XX,YY,ZZ,[0 0],'k','linewidth',0.5);
contour(lon_Chl,lat_Chl,tm,[0.05,0.1],'color',[.4 .4 .4]);
text(90,-13.5,'0.05','fontsize',8,'fontweight','bold','color',[.4 .4 .4],'rotation',0);
text(190,-9.4,'0.05','fontsize',8,'fontweight','bold','color',[.4 .4 .4],'rotation',-10);
text(120,-11,'0.1','fontsize',8,'fontweight','bold','color',[.4 .4 .4],'rotation',-30);
text(155,-14,'0.1','fontsize',8,'fontweight','bold','color',[.4 .4 .4],'rotation',-30);
text(145,-11,'0.1','fontsize',8,'fontweight','bold','color',[.4 .4 .4],'rotation',-90);
shading flat;
caxis([-0.29 0.29]);
cmocean('balance');
xlim(xlim_Reg);ylim(ylim_Reg);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
set(gca,'fontsize',8,'fontweight','bold','xtick',80:20:200,'xticklabels',...
    {'80°E','100°E','120°E','140°E','160°E','180°','160°W'},'ytick',-25:5:-10,...
    'yticklabels',{'25°S','20°S','15°S','10°S'});
set(gca,'xaxislocation','top');
c = colorbar('southoutside');
set(c,'position',[0.12 0.542 0.21 0.01]);
text(104,-26,'mg/m^3','fontsize',8,'fontweight','bold','color','k');
plot([lon_Chl(1,2951) lon_Chl(1,3000)],[lat_Chl(465,1) lat_Chl(465,1)],'c','linewidth',1.2);
plot([lon_Chl(1,2951) lon_Chl(1,3000)],[lat_Chl(520,1) lat_Chl(520,1)],'c','linewidth',1.2);
plot([lon_Chl(1,2951) lon_Chl(1,2951)],[lat_Chl(465,1) lat_Chl(520,1)],'c','linewidth',1.2);
plot([lon_Chl(1,3000) lon_Chl(1,3000)],[lat_Chl(465,1) lat_Chl(520,1)],'c','linewidth',1.2);

plot([lon_Chl(1,1811) lon_Chl(1,1880)],[lat_Chl(331,1) lat_Chl(331,1)],'c','linewidth',1.2);
plot([lon_Chl(1,1811) lon_Chl(1,1880)],[lat_Chl(395,1) lat_Chl(395,1)],'c','linewidth',1.2);
plot([lon_Chl(1,1811) lon_Chl(1,1811)],[lat_Chl(331,1) lat_Chl(395,1)],'c','linewidth',1.2);
plot([lon_Chl(1,1880) lon_Chl(1,1880)],[lat_Chl(331,1) lat_Chl(395,1)],'c','linewidth',1.2);

plot([lon_Chl(1,1530) lon_Chl(1,1600)],[lat_Chl(270,1) lat_Chl(270,1)],'c','linewidth',1.2);
plot([lon_Chl(1,1530) lon_Chl(1,1600)],[lat_Chl(320,1) lat_Chl(320,1)],'c','linewidth',1.2);
plot([lon_Chl(1,1530) lon_Chl(1,1530)],[lat_Chl(270,1) lat_Chl(320,1)],'c','linewidth',1.2);
plot([lon_Chl(1,1600) lon_Chl(1,1600)],[lat_Chl(270,1) lat_Chl(320,1)],'c','linewidth',1.2);

plot(147.62,-19.3,'v','markersize',6,'color','b','MarkerFaceColor',[0.31,0.62,0.83]);

text(72,-11.4,'A','fontsize',8,'fontweight','bold','color','k','Backgroundcolor','w','edgecolor','k');
xlabel('Longitude','fontsize',8,'fontweight','bold');
ylab = ylabel('Latitude','fontsize',8,'fontweight','bold');
set(ylab,'Position',[63 -18 -1]);

%%
t = nan(108,6);
ttt = nan(108,3);
temp = nan(721,3360);

loc = [2951,465;3000,520;1811,331;1890,395;1530,270;1600,320];
    
for l = 1:3
    for i = 1:108
        temp = Chl(:,:,i);
        tt = temp(loc(2*l-1,2):loc(2*l,2),loc(2*l-1,1):loc(2*l,1));
        t(i,l) = mean(mean(tt),'omitnan');
        tt = tt(:); tt = tt(find(isnan(tt)==0));len = length(tt);
        ttt(i,l) = std(tt,'omitnan')/4;
        
        temp = Chl_cli(:,:,i);
        tt = temp(loc(2*l-1,2):loc(2*l,2),loc(2*l-1,1):loc(2*l,1));
        t(i,l+3) = mean(mean(tt),'omitnan'); 
    end
end

ttt(47,1) = (ttt(46,1)+ttt(48,1))/2;
subplot('position',[0.1 0.37 0.55 0.12]);
plot(t(:,1),'r','linewidth',1.4);
hold on
lowt = [1:108,108:-1:1]'; uppt = [t(:,1)-ttt(:,1);t(108:-1:1,1)+ttt(108:-1:1,1)];
uppt(find(uppt < 0)) = 0;
patch(lowt,uppt,[0.85 0.1 0.15],'FaceAlpha',0.2,'linestyle','none');
plot([14 14],[-0.2 1.2],'-.k','linewidth',0.5);
plot(t(:,4),'--','linewidth',1.5,'color',[.7 .7 .7]);
plot(t(:,1),'r','linewidth',1.4);
xlim([1 59]);
ylim([0.02 0.72]);
set(gca,'fontsize',8,'fontweight','bold');
set(gca,'xtick',[1,11,21,32,42,52],'xticklabels','');
text(2.5,0.53,'B','fontsize',8,'fontweight','bold','color','k','Backgroundcolor','w','edgecolor','k');

subplot('position',[0.1 0.23 0.55 0.12]);
plot(t(:,2),'r','linewidth',1.4);
hold on
lowt = [1:108,108:-1:1]'; uppt = [t(:,2)-ttt(:,2);t(108:-1:1,2)+ttt(108:-1:1,2)];
uppt(find(uppt < 0)) = 0;
patch(lowt,uppt,[0.85 0.1 0.15],'FaceAlpha',0.2,'linestyle','none');
plot([17 17],[0 1.2],'-.k','linewidth',0.5);
plot(t(:,5),'--','linewidth',1.5,'color',[.7 .7 .7]);
plot(t(:,2),'r','linewidth',1.4);
xlim([1 59]);
ylim([0.12 1.15]);
set(gca,'fontsize',8,'fontweight','bold');
set(gca,'xtick',[1,11,21,32,42,52],'xticklabels',{''},'ytick',0.2:0.4:1);
text(12.2,-0.123,'January','fontsize',8,'fontweight','bold');
text(42.4,-0.123,'Feburary','fontsize',8,'fontweight','bold');
ylab = ylabel('mg/m^3','fontsize',8,'fontweight','bold');
set(ylab,'Position',[-2.6 0.63 -1]);
set(gca,'xtick',[1,11,21,32,42,52],'xticklabels','');
text(2.5,0.865,'C','fontsize',8,'fontweight','bold','color','k','Backgroundcolor','w','edgecolor','k');

subplot('position',[0.1 0.09 0.55 0.12]);
plot(t(:,3),'r','linewidth',1.4);
hold on
lowt = [1:108,108:-1:1]'; uppt = [t(:,3)-ttt(:,3);t(108:-1:1,3)+ttt(108:-1:1,3)];
uppt(find(uppt < 0)) = 0;
patch(lowt,uppt,[0.85 0.1 0.15],'FaceAlpha',0.2,'linestyle','none');
plot([18 18],[0 1.2],'-.k','linewidth',0.5);
plot(t(:,6),'--','linewidth',1.5,'color',[.7 .7 .7]);
plot(t(:,3),'r','linewidth',1.4);
xlim([1 59]);
ylim([0.05 1]);
set(gca,'fontsize',8,'fontweight','bold');
set(gca,'xtick',[1,11,21,32,42,52],'xticklabels',{'1','10','20','1','10','20'},'ytick',0.2:0.4:1);
text(12.8,-0.25,'January','fontsize',8,'fontweight','bold');
text(43.8,-0.25,'Feburary','fontsize',8,'fontweight','bold');
text(2.5,0.75,'D','fontsize',8,'fontweight','bold','color','k','Backgroundcolor','w','edgecolor','k');

%%
subplot('position',[0.67 0.09 0.23 0.4]);
xr= ((chl_mean+chl_mean_std)');  yr=dep_int(7:76,1)';
xL=((chl_mean-chl_mean_std)');
H=fill([xL(1,1) xr(1,1) xr xr(1,end) xL(1,end) fliplr(xL) xL(1,1)],[yr(1,1) ...
    yr(1,1) yr yr(1,end) yr(1,end) fliplr(yr) yr(1,1)],[0.6510 0.8078 0.8902]);
set(H,'linestyle','none','faceAlpha',0.4)
hold on
plot(chl_mean,dep_mean,'linewid',1.5,'color',[0.31,0.62,0.83])
plot(smoothdata(chl_int(:,13)),depth_int(:,13),'linewid',1.5,'color',[0.99,0.56,0.55])
set(gca,'ydir','reverse','xtick',[0.1:0.1:0.5],'ytick',[3:3:15])
set(gca,'yAxisLocation','right');
grid on
set(gca,'TickDir','in','layer','top','fontsize',8,'fontweight','bold','GridLineStyle','--')
xlab = xlabel('mg/m^3','fontsize',8,'fontweight','bold');
set(xlab,'Position',[0.25 15.8 -1]);
ylab = ylabel('Depth (m)','fontsize',8,'fontweight','bold');
set(ylab,'Position',[0.45 9.1 -1],'rotation',270);
ylim([2.9 15]);
xlim([0.08 0.4])
text(0.1,4,'E','fontsize',8,'fontweight','bold','color','k','Backgroundcolor','w','edgecolor','k');