function plot_Td_granule_map(lev,y,d,g)

% 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/';

fig=figure;
Position = [100 100 1000 370];
fig.Position=Position;

t=tiledlayout(1,3);
t.Padding='compact';

f=gobjects(3);


f(1)=nexttile;
set(f(1),'Fontsize',10)

cfile=[datadir,'airs_',num2str(y),'_',num2str(d,'%03d'),'_',num2str(g,'%03d'),'.nc'];

data = ncread(cfile,'ret_temp');
lat  = ncread(cfile,'l1_lat');
lon  = ncread(cfile,'l1_lon');
time = ncread(cfile,'l1_time');
load('../grid_granule.mat','xg','yg','zg')
zq=knnsearch(zg,lev);
tst = int32(min(time,[],"all","omitnan"));
dy = idivide(int32(tst-596937600),int32(86400))+1;
hr = idivide(int32(tst-596937600-86400*(dy-1)),int32(3600));
% mi = idivide(int32(tst-596937600-86400*(dy-1)-3600*hr),int32(60));
if idivide(int32(dy),int32(32))==1
    mn = 1;
    dy = dy - 31;
    cmn = 'JAN';
else
    mn = 12;
    cmn = 'DEC';
end
starttimeindayssince2000 = double(tst)/60/60/24;
StartTime = datenum('01-Jan-2000') + starttimeindayssince2000;
hours = round((StartTime - floor(StartTime)) * 24);
mins = round((StartTime*24 - floor(StartTime*24)) * 60);
roundedtime = floor(StartTime) + (hours/24) + (mins/60/24);
str = datestr(roundedtime,'yyyymmdd-HHMM');
ctitle=append(str(10:11),":",str(12:13),"UTC ",str(7:8)," ",cmn," ",str(1:4));

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

xmin=max(locx(1,:),[],'all');
xmax=min(locx(90,:),[],'all');
ymin=max(locy(:,1),[],'all');
ymax=min(locy(:,135),[],'all');
[ygm,xgm]=meshgrid(yg((yg>ymin)&(yg<ymax)),xg((xg>xmin)&(xg<xmax)));
var = squeeze(data(zq,:,:));
vq  = griddata(double(locx),double(locy),double(var),double(xgm),double(ygm));
[td,tb] = nph_airs_4dp_detrend(vq,1,4);
td(1,:)=NaN;
td(end,:)=NaN;
td(:,1)=NaN;
td(:,end)=NaN;

lonq = griddata(double(locx),double(locy),double(lon),double(xgm),double(ygm));
latq = griddata(double(locx),double(locy),double(lat),double(xgm),double(ygm));
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);
lon1=reshape(lonq,[],1);
lat1=reshape(latq,[],1);
td1=reshape(td,[],1);
Vq=NaN(size(Xq));
F=scatteredInterpolant(lon1,lat1,td1,'linear','none');
Vq=F(Xq,Yq);    

axis([lonmin lonmax latmin latmax])
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' '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(1),Xq,Yq,Vq,'LineStyle','none','LevelStep',0.01);
colorbar(f(1),'eastoutside')
load('colormap_blue-red.mat')
colormap(f(1),CustomColormap1)
clim(f(1),[-4 4])
boundingbox = [lonmin latmin ; lonmax latmax];
nph_draw_coastline(boundingbox,0,20,'color','k');

title(f(1),'AIRS',ctitle,'FontSize',12,'FontWeight','normal','FontName','Helvetica')

%--- JAGUAR w/ obs filter ---%
f(2)=nexttile;
set(f(2),'Fontsize',10)

if g<=235
    if (d>=335) && (d<=365)
        nyr = y*10000+1200+d-334;
        cdy =int2str(nyr);
    elseif (d>=1) && (d<=31)
        nyr = y*10000+100+d;
        cdy =int2str(nyr);
    end
else
    if (d>=335) && (d<=364)
        nyr = y*10000+1200+d-333;
        cdy =int2str(nyr);
    elseif (d==365)
        cdy =[int2str(y+1) '0101'];
    elseif (d>=1) && (d<=30)
        nyr = y*10000+100+d+1;
        cdy =int2str(nyr);
    end
end
cfile=[datadir,cdy,'-*30_Model_T_asAIRS',num2str(g,'%0.3d'),'.mat'];
ctitle=append(F.name(10:11),":30","UTC ",str(7:8)," ",cmn," ",str(1:4));

M = load(cfile);
data = permute(M.Model.Ta_sm,[3 1 2]);
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
var = squeeze(data(zq,:,:));
vq  = griddata(double(locx),double(locy),double(var),double(xgm),double(ygm));
[td,tb] = nph_airs_4dp_detrend(vq,1,4);
td(1,:)=NaN;
td(end,:)=NaN;
td(:,1)=NaN;
td(:,end)=NaN;

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

axis([lonmin lonmax latmin latmax])
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' '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(2),Xq,Yq,Vq,'LineStyle','none','LevelStep',0.01);
colorbar(f(2),'eastoutside')
load('colormap_blue-red.mat')
colormap(f(2),CustomColormap1)
clim(f(2),[-4 4])
boundingbox = [lonmin latmin ; lonmax latmax];
nph_draw_coastline(boundingbox,0,20,'color','k');
title(f(2),'JAGUAR w/ obs. filter',ctitle,'FontSize',12,'FontWeight','normal','FontName','Helvetica')

%--- JAGUAR w/o obs. filter ---%
f(3)=nexttile;
set(f(3),'Fontsize',10)

cfile=[datadir,cdy,'-*30_Model_T_asAIRS',num2str(g,'%0.3d'),'_300m.mat'];
M = load(cfile);
data = permute(M.Model.Ta,[3 1 2]);
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
var = squeeze(data(zq,:,:));
vq  = griddata(double(locx),double(locy),double(var),double(xgm),double(ygm));
[td,tb] = nph_airs_4dp_detrend(vq,1,4);
td(1,:)=NaN;
td(end,:)=NaN;
td(:,1)=NaN;
td(:,end)=NaN;

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

axis([lonmin lonmax latmin latmax])
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' '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(3),Xq,Yq,Vq,'LineStyle','none','LevelStep',0.01);
colorbar(f(3),'eastoutside')
load('colormap_blue-red.mat')
colormap(f(3),CustomColormap1)
clim(f(3),[-6 6])
boundingbox = [lonmin latmin ; lonmax latmax];
nph_draw_coastline(boundingbox,0,20,'color','k');
title(f(3),'JAGUAR w/o obs. filter',ctitle,'FontSize',12,'FontWeight','normal','FontName','Helvetica')
end

