function plot_JAGUAR_fullres_3DST_20181216_2030_xy(lev)

% Download functions from Hindley et al. (2021) from
% https://doi.org/10.5281/zenodo.4721883 (Hindley, 2021)
addpath("../nhindley-acp-2020-465-2d5b3b8/AnalysisCode/", "../nhindley-acp-2020-465-2d5b3b8/Functions/")
datadir='../data/';
y=2018;
d=350;
g=196;
cdy='20181216';

fig=figure;
Position = [100 100 700 640];
fig.Position=Position;

t=tiledlayout(2,2);
t.Padding='compact';

f=gobjects(4);


cfile=[datadir '20181216-2030_Model_T_asAIRS196_300m.mat'];
cfile3DST=[datadir 'JAGUAR_fullres_3DST_20181216-2030_g196.mat'];

M = load(cfile);
load(cfile3DST,"MFx");
load(cfile3DST,"MFy");
load(cfile3DST,"vwn");
load(cfile3DST,"zwn");
load(cfile3DST,"mwn");
load(cfile3DST,"amp");
direction = vwn./abs(vwn);
for i=1:size(vwn,2)
    for j=1:size(vwn,3)
            zwn(:,i,j)=-zwn(:,i,j).*direction(:,i,j);
            mwn(:,i,j)=-mwn(:,i,j).*direction(:,i,j);
            vwn(:,i,j)= vwn(:,i,j).*direction(:,i,j);
            amp(:,i,j)= amp(:,i,j);
            MFx(:,i,j)=-MFx(:,i,j);
            MFy(:,i,j)=-MFy(:,i,j);
    end
end
lon  = squeeze(M.xg(:,:,1));
lonmin = min(lon,[],"all","omitnan");
lonmax = max(lon,[],"all","omitnan");
lon(lon>360) = lon(lon>360)-360; %0 - 360
if lonmax>360
    lonmax = lonmax-360;
end
lon(lon>180) = lon(lon>180)-360; %-180 - 180
if lonmin>180
    lonmin = lonmin-360;
end
if lonmax>180
    lonmax = lonmax-360;
end
lat  = squeeze(M.yg(:,:,1));
zg   = M.Model.Alt;
zq=knnsearch(zg,lev);
distx = zeros(89,135);
for j=1:135
    for i=1:89
        distx(i,j)=nph_haversine([lat(i,j) lon(i,j)],[lat(i+1,j) lon(i+1,j)]);
    end
end
disty = zeros(90,134);
for j=1:134
    for i=1:90
        disty(i,j)=nph_haversine([lat(i,j) lon(i,j)],[lat(i,j+1) lon(i,j+1)]);
    end
end

locx = zeros(90,135);
for j=1:135
    for i=44:-1:1
        locx(i,j)=locx(i+1,j)-distx(i,j);
    end
    for i=46:90
        locx(i,j)= locx(i-1,j)+distx(i-1,j);
    end
end
locy = zeros(90,135);
for j=2:135
    for i=1:90
        locy(i,j)= locy(i,j-1)+disty(i,j-1);
    end
end

load('../grid_granule.mat','xg','yg')
xmin=max(locx(1,:),[],'all');
xmax=min(locx(90,:),[],'all');
ymax=min(locy(:,135),[],'all');
[ygm,xgm]=meshgrid(yg(yg<=ymax),xg((xg>=xmin)&(xg<=xmax)));

lonq = griddata(double(locx),double(locy),double(lon),double(xgm),double(ygm));
latq = griddata(double(locx),double(locy),double(lat),double(xgm),double(ygm));

%Flon = scatteredInterpolant(locx,locy,lon,"linear","linear");
%lonq = Flon(xgm,ygm);
%Flat = scatteredInterpolant(locx,locy,lat,"linear","linear");
%latq = Flat(xgm,ygm);
lon1=reshape(lonq,[],1);
lat1=reshape(latq,[],1);
lonmax = max(lonq,[],"all");
lonmin = min(lonq,[],"all");
latmax = max(latq,[],"all");
latmin = min(latq,[],"all");

latg=linspace(latmin,latmax,100);
long=linspace(lonmin,lonmax,100);
[Xq,Yq]=meshgrid(long,latg);

%------------------------------%
%--- Lon-Lat maps           ---%
%------------------------------%
i=1;
f(i)=nexttile;
set(f(i),'Fontsize',10)

td = squeeze(amp(zq,:,:));
td(1,:)=NaN;
td(end,:)=NaN;
td(:,1)=NaN;
td(:,end)=NaN;

td1=double(reshape(td,[],1));
F=scatteredInterpolant(lon1,lat1,td1,'linear','none');
Vq=F(Xq,Yq);    

%axis([lonmin lonmax latmin latmax])
axis([85 126 45 69])
xticks(-180:10:180)
xticklabels({'180W' '170W' '160W' '150W' '140W' '130W' '120W' '110W' '100W' '90W' '80W' '70W' '60W' '50W' '40W' '30W' '20W' '10W' '0' '10E' '20E' '30E' '40E' '50E' '60E' '70E' '80E' '90E' '100E' '110E' '120E' '130E' '140E' '150E' '160E' '170E' '180E'})
yticks(-90:5:90)
yticklabels({'90S' '85S' '80S' '75S' '70S' '65S' '60S' '55S' '50S' '45S' '40S' '35S' '30S' '25S' '20S' '15S' '10S' '5S' 'EQ' '5N' '10N' '15N' '20N' '25N' '30N' '35N' '40N' '45N' '50N' '55N' '60N' '65N' '70N' '75N' '80N' '85N' '90N'})
hold on
box on
contourf(f(i),Xq,Yq,Vq,'LineStyle','none','LevelStep',0.01);
colorbar(f(i),'eastoutside')
load('colormap_blue-red.mat')
colormap(f(i),'jet')
clim([0 8])
boundingbox = [lonmin latmin ; lonmax latmax];
nph_draw_coastline(boundingbox,0,20,'color','k');
title(f(i),'Amplitude (K)','FontSize',12,'FontWeight','normal','FontName','Helvetica')
%------------------------------%
i=2;
f(i)=nexttile;
set(f(i),'Fontsize',10)

td = squeeze(1/vwn(zq,:,:));
td(1,:)=NaN;
td(end,:)=NaN;
td(:,1)=NaN;
td(:,end)=NaN;

td1=double(reshape(td,[],1));
F=scatteredInterpolant(lon1,lat1,td1,'linear','none');
Vq=F(Xq,Yq);    

%axis([lonmin lonmax latmin latmax])
axis([85 126 45 69])
xticks(-180:10:180)
xticklabels({'180W' '170W' '160W' '150W' '140W' '130W' '120W' '110W' '100W' '90W' '80W' '70W' '60W' '50W' '40W' '30W' '20W' '10W' '0' '10E' '20E' '30E' '40E' '50E' '60E' '70E' '80E' '90E' '100E' '110E' '120E' '130E' '140E' '150E' '160E' '170E' '180E'})
yticks(-90:5:90)
yticklabels({'90S' '85S' '80S' '75S' '70S' '65S' '60S' '55S' '50S' '45S' '40S' '35S' '30S' '25S' '20S' '15S' '10S' '5S' 'EQ' '5N' '10N' '15N' '20N' '25N' '30N' '35N' '40N' '45N' '50N' '55N' '60N' '65N' '70N' '75N' '80N' '85N' '90N'})
hold on
box on
contourf(f(i),Xq,Yq,Vq,'LineStyle','none','LevelStep',0.5);
colorbar(f(i),'eastoutside')
colormap(f(i),'jet')
clim([5 25])
boundingbox = [lonmin latmin ; lonmax latmax];
nph_draw_coastline(boundingbox,0,20,'color','k');
title(f(i),'Vertical Wavelength (km)','FontSize',12,'FontWeight','normal','FontName','Helvetica')
%------------------------------%
i=3;
f(i)=nexttile;
set(f(i),'Fontsize',10)

td = squeeze(MFx(zq,:,:))*1000;
td(1,:)=NaN;
td(end,:)=NaN;
td(:,1)=NaN;
td(:,end)=NaN;

td1=double(reshape(td,[],1));
F=scatteredInterpolant(lon1,lat1,td1,'linear','none');
Vq=F(Xq,Yq);    

%axis([lonmin lonmax latmin latmax])
axis([85 126 45 69])
xticks(-180:10:180)
xticklabels({'180W' '170W' '160W' '150W' '140W' '130W' '120W' '110W' '100W' '90W' '80W' '70W' '60W' '50W' '40W' '30W' '20W' '10W' '0' '10E' '20E' '30E' '40E' '50E' '60E' '70E' '80E' '90E' '100E' '110E' '120E' '130E' '140E' '150E' '160E' '170E' '180E'})
yticks(-90:5:90)
yticklabels({'90S' '85S' '80S' '75S' '70S' '65S' '60S' '55S' '50S' '45S' '40S' '35S' '30S' '25S' '20S' '15S' '10S' '5S' 'EQ' '5N' '10N' '15N' '20N' '25N' '30N' '35N' '40N' '45N' '50N' '55N' '60N' '65N' '70N' '75N' '80N' '85N' '90N'})
hold on
box on
contourf(f(i),Xq,Yq,Vq,'LineStyle','none','LevelStep',0.1);
colorbar(f(i),'eastoutside')
colormap(f(i),CustomColormap1)
clim([-20 20])
boundingbox = [lonmin latmin ; lonmax latmax];
nph_draw_coastline(boundingbox,0,20,'color','k');
title(f(i),'Zonal Momentum Flux (mPa)','FontSize',12,'FontWeight','normal','FontName','Helvetica')
%------------------------------%
i=4;
f(i)=nexttile;
set(f(i),'Fontsize',10)

td = squeeze(MFy(zq,:,:))*1000;
td(1,:)=NaN;
td(end,:)=NaN;
td(:,1)=NaN;
td(:,end)=NaN;

td1=double(reshape(td,[],1));
F=scatteredInterpolant(lon1,lat1,td1,'linear','none');
Vq=F(Xq,Yq);    

%axis([lonmin lonmax latmin latmax])
axis([85 126 45 69])
xticks(-180:10:180)
xticklabels({'180W' '170W' '160W' '150W' '140W' '130W' '120W' '110W' '100W' '90W' '80W' '70W' '60W' '50W' '40W' '30W' '20W' '10W' '0' '10E' '20E' '30E' '40E' '50E' '60E' '70E' '80E' '90E' '100E' '110E' '120E' '130E' '140E' '150E' '160E' '170E' '180E'})
yticks(-90:5:90)
yticklabels({'90S' '85S' '80S' '75S' '70S' '65S' '60S' '55S' '50S' '45S' '40S' '35S' '30S' '25S' '20S' '15S' '10S' '5S' 'EQ' '5N' '10N' '15N' '20N' '25N' '30N' '35N' '40N' '45N' '50N' '55N' '60N' '65N' '70N' '75N' '80N' '85N' '90N'})
hold on
box on
contourf(f(i),Xq,Yq,Vq,'LineStyle','none','LevelStep',0.1);
colorbar(f(i),'eastoutside')
colormap(f(i),CustomColormap1)
clim([-20 20])
boundingbox = [lonmin latmin ; lonmax latmax];
nph_draw_coastline(boundingbox,0,20,'color','k');
title(f(i),'Meridional Momentum Flux (mPa)','FontSize',12,'FontWeight','normal','FontName','Helvetica')


%saveas(fig,['Figures/granules/',num2str(y),'_',num2str(d,'%03d'),'_',num2str(g,'%03d'),'_map.png'])
end

