%    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.3----UVAI_Wind field figure
clear;
load("UVAI_Wind.mat");
load("HYSPLIT_simulation.mat")

tg_x = -175.3801033248204;%longitude of HTHH
tg_y = -20.568822097017062;%latitude of HTHH

pos = [0.06 0.6 0.43 0.09;0.5 0.6 0.43 0.09;...
       0.06 0.5 0.43 0.09;0.5 0.5 0.43 0.09;...
       0.06 0.4 0.43 0.09;0.5 0.4 0.43 0.09;...
       0.06 0.3 0.43 0.09;0.5 0.3 0.43 0.09;...
       0.94 0.45 0.01 0.2;0.94 0.31 0.01 0.07;];

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

yt = -25:5:-10;
xt = 80:20:200;
str1 = string(80:20:160)+'°E';
str2 = '180°';
str3 = string(160)+'°W';
labelStr = [str1,str2,str3];
color = [242 242 242;230 230 230;218 219 201;255 189 0;255 155 0;242 81 34;238 14 0;176 0 0;112 0 0;44 0 0]./256;
figure(1)
% a)
ax1 = subplot('position',pos(1,:));
temp = UVAI(:,:,1);
temp(lon_UVAI<=70.3)=nan;
temp(lon_UVAI>=199.7)=nan;
temp(lat_UVAI<=-27.7)=nan;
temp(lat_UVAI>=-8.3)=nan;
pcolor(lon_UVAI,lat_UVAI,temp);shading flat
hold on
quiver(lon_wind,lat_wind,u(:,:,1),v(:,:,1),0.4,'Color','#009051','LineWidth',0.5,'MaxHeadSize',1);
hold on
contour(lon_UVAI,lat_UVAI,Coast,[1,1],'-k')
colormap(color);
xlim(xlim_Reg);ylim(ylim_Reg);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
yticks(yt);
yticklabels(cellstr(string(25:-5:10)+'°S'));
caxis([0.5 3]);
hold on
plot(mod(360+tg_x,360),tg_y,'^k','LineWidth',0.5,'MarkerSize',10,'MarkerFaceColor','c');
hold on
text(73,-12,'                           ','Color','k','fontsize',6,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
quiver(73,-12,20,0,0.4,'Color','#009051','LineWidth',0.5,'MaxHeadSize',1);
quiver(73,-20,30,0,0.4,'Color','#009051','LineWidth',0.5,'MaxHeadSize',1);
hold on
text(83,-12,'20 m/s','Color','k','fontsize',10,'FontWeight','bold');
text(73,-23,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(73,-23,'A','Color','k','fontsize',11.5,'FontWeight','bold');
text(79,-23,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(79,-23,'15 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
set(gca,'ticklength',[0 0])

% b)
ax2 = subplot('position',pos(2,:));
temp = UVAI(:,:,2);
temp(lon_UVAI<=70.3)=nan;
temp(lon_UVAI>=199.7)=nan;
temp(lat_UVAI<=-27.7)=nan;
temp(lat_UVAI>=-8.3)=nan;
pcolor(lon_UVAI,lat_UVAI,temp);shading flat
hold on
quiver(lon_wind,lat_wind,u(:,:,2),v(:,:,2),0.4,'Color','#009051','LineWidth',0.5);
hold on
contour(lon_UVAI,lat_UVAI,Coast,[1,1],'-k')
colormap(ax1,color);
xlim(xlim_Reg);ylim(ylim_Reg);
colormap(ax2,color);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
caxis([0.5 3]);
hold on
plot(mod(360+tg_x,360),tg_y,'^k','LineWidth',0.5,'MarkerSize',10,'MarkerFaceColor','c');
text(73,-23,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(73,-23,'B','Color','k','fontsize',11.5,'FontWeight','bold');
text(79,-23,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(79,-23,'16 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
set(gca,'ticklength',[0 0])

% c)
ax3 = subplot('position',pos(3,:));
temp = UVAI(:,:,3);
temp(lon_UVAI<=70.3)=nan;
temp(lon_UVAI>=199.7)=nan;
temp(lat_UVAI<=-27.7)=nan;
temp(lat_UVAI>=-8.3)=nan;
pcolor(lon_UVAI,lat_UVAI,temp);shading flat
hold on
quiver(lon_wind,lat_wind,u(:,:,3),v(:,:,3),0.4,'Color','#009051','LineWidth',0.5);
hold on
contour(lon_UVAI,lat_UVAI,Coast,[1,1],'-k')
colormap(ax3,color);
xlim(xlim_Reg);ylim(ylim_Reg);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
yticks(yt);
yticklabels(cellstr(string(25:-5:10)+'°S'));
caxis([0.5 3]);
hold on
plot(mod(360+tg_x,360),tg_y,'^k','LineWidth',0.5,'MarkerSize',10,'MarkerFaceColor','c');
text(73,-23,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(73,-23,'C','Color','k','fontsize',11.5,'FontWeight','bold');
text(79,-23,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(79,-23,'17 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
text(55,-35,'Latitude','fontsize',11.5,'FontWeight','bold','rotation',90)
set(gca,'ticklength',[0 0])

% d)
ax4 = subplot('position',pos(4,:));
temp = UVAI(:,:,4);
temp(lon_UVAI<=70.3)=nan;
temp(lon_UVAI>=199.7)=nan;
temp(lat_UVAI<=-27.7)=nan;
temp(lat_UVAI>=-8.3)=nan;
pcolor(lon_UVAI,lat_UVAI,temp);shading flat
hold on
quiver(lon_wind,lat_wind,u(:,:,4),v(:,:,4),0.4,'Color','#009051','LineWidth',0.5);
hold on
contour(lon_UVAI,lat_UVAI,Coast,[1,1],'-k')
colormap(ax4,color);
xlim(xlim_Reg);ylim(ylim_Reg);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
caxis([0.5 3]);
hold on
plot(mod(360+tg_x,360),tg_y,'^k','LineWidth',0.5,'MarkerSize',10,'MarkerFaceColor','c');
text(73,-23,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(73,-23,'D','Color','k','fontsize',11.5,'FontWeight','bold'); 
text(79,-23,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(79,-23,'18 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
text(217,-5,'Aerosol Index','fontsize',11.5,'FontWeight','bold','rotation',-90)
set(gca,'ticklength',[0 0])

% e)
ax5 = subplot('position',pos(5,:));
temp = UVAI(:,:,5);
temp(lon_UVAI<=70.3)=nan;
temp(lon_UVAI>=199.7)=nan;
temp(lat_UVAI<=-27.7)=nan;
temp(lat_UVAI>=-8.3)=nan;
pcolor(lon_UVAI,lat_UVAI,temp);shading flat
hold on
quiver(lon_wind,lat_wind,u(:,:,5),v(:,:,5),0.4,'Color','#009051','LineWidth',0.5);
hold on
contour(lon_UVAI,lat_UVAI,Coast,[1,1],'-k')
colormap(ax5,color);
xlim(xlim_Reg);ylim(ylim_Reg);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
yticks(yt);
yticklabels(cellstr(string(25:-5:10)+'°S'));
caxis([0.5 3]);
hold on
plot(mod(360+tg_x,360),tg_y,'^k','LineWidth',0.5,'MarkerSize',10,'MarkerFaceColor','c');
text(73,-23,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(73,-23,'E','Color','k','fontsize',11.5,'FontWeight','bold');
text(79,-23,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(79,-23,'19 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
set(gca,'ticklength',[0 0])

% f)
ax6 = subplot('position',pos(6,:));
temp = UVAI(:,:,6);
temp(lon_UVAI<=70.3)=nan;
temp(lon_UVAI>=199.7)=nan;
temp(lat_UVAI<=-27.7)=nan;
temp(lat_UVAI>=-8.3)=nan;
pcolor(lon_UVAI,lat_UVAI,temp);shading flat
hold on
quiver(lon_wind,lat_wind,u(:,:,6),v(:,:,6),0.4,'Color','#009051','LineWidth',0.5);
hold on
contour(lon_UVAI,lat_UVAI,Coast,[1,1],'-k')
colormap(ax6,color);
xlim(xlim_Reg);ylim(ylim_Reg);
colormap(ax4,color);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
caxis([0.5 3]);
hold on
plot(mod(360+tg_x,360),tg_y,'^k','LineWidth',0.5,'MarkerSize',10,'MarkerFaceColor','c');
text(73,-23,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(73,-23,'F','Color','k','fontsize',11.5,'FontWeight','bold');
text(79,-23,'                       ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w');
text(79,-23,'20 Jan 2022','Color','k','fontsize',11.5,'FontWeight','bold');
set(gca,'ticklength',[0 0])
c1 = colorbar('eastoutside');
set(c1,'position',pos(9,:));
ct = 0.5:0.5:3;
set(c1,'Ticks',ct,'TickLabels',string(ct))

%g)
c = [1 86 153;250 192 15;243 118 74;95 198 201;79 89 109]/256;
ax7 = subplot('position',pos(7,:));
w = 1;
contour(lon_UVAI,lat_UVAI,UVAI(:,:,1),[1.46,1.46],'LineWidth',w,'LineColor',c(3,:))
hold on
contour(lon_UVAI,lat_UVAI,UVAI(:,:,2),[1.37,1.37],'LineWidth',w,'LineColor',c(3,:))
hold on
contour(lon_UVAI,lat_UVAI,UVAI(:,:,3),[1.23,1.23],'LineWidth',w,'LineColor',c(3,:))
hold on
contour(lon_UVAI,lat_UVAI,UVAI(:,:,4),[1.43,1.43],'LineWidth',w,'LineColor',c(3,:))
hold on
contour(lon_UVAI,lat_UVAI,UVAI(:,:,5),[1.44,1.44],'LineWidth',w,'LineColor',c(3,:))
hold on
contour(lon_UVAI,lat_UVAI,Coast,[1,1],'-k')
xlim(xlim_Reg);ylim(ylim_Reg);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
yticks(yt);
yticklabels(cellstr(string(25:-5:10)+'°S'));
xticks(xt(1:end-1));
xticklabels(cellstr(labelStr(1:end-1)));
hold on
plot(mod(360+tg_x,360),tg_y,'^k','LineWidth',0.5,'MarkerSize',10,'MarkerFaceColor','c');
text(73,-23,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(73,-23,'G','Color','k','fontsize',11.5,'FontWeight','bold');
text(171,-23,'15 Jan','Color','k','fontsize',10,'FontWeight','bold');
text(149,-23,'16 Jan','Color','k','fontsize',10,'FontWeight','bold');
text(128,-20,'17 Jan','Color','k','fontsize',10,'FontWeight','bold');
text(107,-17,'18 Jan','Color','k','fontsize',10,'FontWeight','bold');
text(86,-13,'19 Jan','Color','k','fontsize',10,'FontWeight','bold');
set(gca,'ticklength',[0 0])

%h)
ax8 = subplot('position',pos(8,:));
color_traj = [230 230 230;114 208 255;65 142 77;211 188 101;232 73 0;155 0 0]./256;
prbnan = HYSPLIT_simulation;
prbnan(lon_UVAI<=70.3) = nan;
pcolor(lon_UVAI,lat_UVAI,prbnan);shading flat
caxis([0 20])
colormap(ax8,color_traj);
hold on
contour(lon_UVAI,lat_UVAI,Coast,[1,1],'-k')
xlim(xlim_Reg);ylim(ylim_Reg);
hold on
cc = ["k","w","k","w","k"];
for i = 1:5
    temp_lon = lon_typical(:,i);
    temp_lat = lat_typical(:,i);
    temp_lon(find(temp_lon==0))=[];
    temp_lat(find(temp_lat==0))=[];
    l = length(temp_lon);
    plot(temp_lon,temp_lat,cc(i),'linewidth',2);
    plot(temp_lon(1),temp_lat(1),'ok','MarkerSize',5,'MarkerFaceColor','.7,.7,.7');
    hold on
end
plot(temp_lon(length(temp_lon)),temp_lat(length(temp_lat)),'ok','MarkerSize',5,'MarkerFaceColor','.7,.7,.7');
plot(mod(360+tg_x,360),tg_y,'^k','LineWidth',0.5,'MarkerSize',10,'MarkerFaceColor','c');
text(173,-24,'15 Jan','Color','k','fontsize',10,'FontWeight','bold','rotation',5);
text(158,-22,'16 Jan','Color','k','fontsize',10,'FontWeight','bold','rotation',-5);
text(142,-18.5,'17 Jan','Color','k','fontsize',10,'FontWeight','bold','rotation',-5);
text(122,-18.5,'18 Jan','Color','k','fontsize',10,'FontWeight','bold');
text(101.5,-16,'19 Jan','Color','k','fontsize',10,'FontWeight','bold','rotation',-5);
text(73,-23,'   ','Color','k','fontsize',10,'FontWeight','bold','BackgroundColor','w','edgecolor','k');
text(73,-23,'H','Color','k','fontsize',11.5,'FontWeight','bold');set(gca,'plotboxaspectratio',[1 ytam/xtam 1]);
set(gca,'Xtick',[],'YTick',[],'fontsize',11.5,'FontWeight','bold','linewidth',1);
xticks(xt)
xticklabels(cellstr(labelStr));
set(gca,'ticklength',[0 0])

c2 = colorbar('eastoutside');
set(c2,'position',pos(10,:));
text(58,-37,'Longitude','fontsize',11.5,'FontWeight','bold');
text(213,-18,'%','fontsize',11.5,'FontWeight','bold','rotation',0,'Interpreter','none')
set(gcf,'Position',[100 100 900 900])