function plot_impala_heatwaves_month_with_most(outdir)

% % Plot month with most HWs %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

cmap = cptcmap('/nfs/see-fs-02_users/earceb/matlab/cptfiles/Cat12.cpt');
load /nfs/see-fs-02_users/earceb/gshhs/shoreline.mat

% Load data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 97p 98.9p
load /nfs/a240/earceb/IMPALA/MetUM/25km/update_26012019/impala_25km_wbmax_HWMId_97p_98pt9p_metrics_3dayHWs.mat
load /nfs/a240/earceb/IMPALA/MetUM/25km/fc/impala_25kmfc_PDBL_97p_98pt9p_wbmax_HWMId_metrics_3dayHWs.mat
load /nfs/a240/earceb/IMPALA/MetUM/4km/update_25012018/impala_4km_wbmax_HWMId_metrics_interp_97p_98pt9p_3dayHWs.mat
load /nfs/a240/earceb/IMPALA/MetUM/4km/fc/impala_4kmfc_PDBL_wbmax_HWMId_metrics_interp_97p_98pt9p_3dayHWs.mat
load /nfs/a240/earceb/IMPALA/ERA5/impala_ERA5_wbmax_HWMId_97p_98pt9p_metrics_3dayHWs.mat

% Read in lats and lons of P25&CP4
lat_25km = ncread('/nfs/a37/IMPALA/data/25km/d03236/d03236_Aday_mean_ay488_25km_20010701-20010730.nc','latitude');
lon_25km = ncread('/nfs/a37/IMPALA/data/25km/d03236/d03236_Aday_mean_ay488_25km_20010701-20010730.nc','longitude')-360;

% Field names for all datasets 
field_names = fieldnames(HWMId_metrics_25km);

% Land-sea mask %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lsm_25km = ncread('/nfs/a277/IMPALA/data/25km/ANCILS/landseamask_ancil_25km.nc','lsm');
lsm_25km = permute(lsm_25km,[2,1]);

% Interpolate lsm onto ERA5 grid and set values<1 to zero
load /nfs/a240/earceb/IMPALA/ERA5/impala_ERA5_wbmax_HWMId_97p_98pt9p_1997-2006.mat
heatwaves_ERA5_1997 = heatwaves_ERA5;
lsm_ERA5 = interp2(lon_25km',lat_25km,lsm_25km,heatwaves_ERA5_1997.lon,heatwaves_ERA5_1997.lat');
[m,n] = size(lsm_ERA5);
for mm=1:m
  ii=find(lsm_ERA5(mm,:)<1);
  lsm_ERA5(mm,ii) = 0;
end  
lsm_ERA5(1:10,:) = 0;
lsm_ERA5(335:345,:) = 0;
lsm_ERA5(:,1:10) = 0;
lsm_ERA5(:,327:337) = 0;

% Set points over ocean to NaN --------------------------------------------
for ff=3:length(field_names)
  HWMId_metrics_25km.(field_names{ff})(lsm_25km==0)=NaN;
  HWMId_metrics_25kmfc_PDBL.(field_names{ff})(lsm_25km==0)=NaN;  
  HWMId_metrics_4km.(field_names{ff})(lsm_25km==0)=NaN;
  HWMId_metrics_4kmfc_PDBL.(field_names{ff})(lsm_25km==0)=NaN; 
  HWMId_metrics_ERA5.(field_names{ff})(lsm_ERA5==0)=NaN;
end

% Plotting %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ERA5
figure
contourf(HWMId_metrics_ERA5.lon,HWMId_metrics_ERA5.lat,HWMId_metrics_ERA5.month_with_most,(1:1:12),'LineColor','None')
caxis([1 12])
colormap(cmap)
colorbar
hold on
plot(shoreline(:,1),shoreline(:,2),'k-','Linewidth',1.5)
set(gca,'FontSize',16)
xlim([-25 55])
ylim([-45 40])
title('Month with most HWs: ERA5')

printname = ['impala_wbmax_HWMId_ERA5_m_with_most_3dayHWs.png'];
eval(['print ''',outdir,'/',printname,''' -dpng -r200']);
close

% 25km
figure
contourf(HWMId_metrics_25km.lon,HWMId_metrics_25km.lat,HWMId_metrics_25km.month_with_most,(1:1:12),'LineColor','None')
caxis([1 12])
colormap(cmap)
colorbar
hold on
plot(shoreline(:,1),shoreline(:,2),'k-','Linewidth',1.5)
set(gca,'FontSize',16)
xlim([-25 55])
ylim([-45 40])
title('Month with most HWs: P25')

% for rr=1:length(regions)
%   plot([left_bnds(rr),left_bnds(rr)],[bot_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)   
%   plot([right_bnds(rr),right_bnds(rr)],[bot_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)
%   plot([left_bnds(rr),right_bnds(rr)],[bot_bnds(rr),bot_bnds(rr)],'k-','LineWidth',1.2)
%   plot([left_bnds(rr),right_bnds(rr)],[top_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)  
% end 

printname = ['impala_wbmax_HWMId_25km_m_with_most_3dayHWs.png'];
eval(['print ''',outdir,'/',printname,''' -dpng -r200']);
close

% 25kmfcPDBL
figure
contourf(HWMId_metrics_25km.lon,HWMId_metrics_25km.lat,HWMId_metrics_25kmfc_PDBL.month_with_most,[1:1:12],'LineColor','None')
caxis([1 12])
colormap(cmap)
colorbar
hold on
plot(shoreline(:,1),shoreline(:,2),'k-','Linewidth',1.5)
set(gca,'FontSize',16)
xlim([-25 55])
ylim([-45 40])
title('Month with most HWs: P25FUT')

% for rr=1:length(regions)
%   plot([left_bnds(rr),left_bnds(rr)],[bot_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)   
%   plot([right_bnds(rr),right_bnds(rr)],[bot_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)
%   plot([left_bnds(rr),right_bnds(rr)],[bot_bnds(rr),bot_bnds(rr)],'k-','LineWidth',1.2)
%   plot([left_bnds(rr),right_bnds(rr)],[top_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)  
% end 

printname = ['impala_wbmax_HWMId_25kmfc_PDBL_m_with_most_3dayHWs.png'];
eval(['print ''',outdir,'/',printname,''' -dpng -r200']);
close

% 4km
figure
contourf(HWMId_metrics_4km.lon,HWMId_metrics_4km.lat,HWMId_metrics_4km.month_with_most,[1:1:12],'LineColor','None')
caxis([1 12])
colormap(cmap)
colorbar
hold on
plot(shoreline(:,1),shoreline(:,2),'k-','Linewidth',1.5)
set(gca,'FontSize',16)
xlim([-25 55])
ylim([-45 40])
title('Month with most HWs: CP4')

% for rr=1:length(regions)
%   plot([left_bnds(rr),left_bnds(rr)],[bot_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)   
%   plot([right_bnds(rr),right_bnds(rr)],[bot_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)
%   plot([left_bnds(rr),right_bnds(rr)],[bot_bnds(rr),bot_bnds(rr)],'k-','LineWidth',1.2)
%   plot([left_bnds(rr),right_bnds(rr)],[top_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)  
% end 

printname = ['impala_wbmax_HWMId_4km_m_with_most_3dayHWs.png'];
eval(['print ''',outdir,'/',printname,''' -dpng -r200']);
close

% 4kmfcPDBL
figure
contourf(HWMId_metrics_4km.lon,HWMId_metrics_4km.lat,HWMId_metrics_4kmfc_PDBL.month_with_most,[1:1:12],'LineColor','None')
caxis([1 12])
colormap(cmap)
colorbar
hold on
plot(shoreline(:,1),shoreline(:,2),'k-','Linewidth',1.5)
set(gca,'FontSize',16)
xlim([-25 55])
ylim([-45 40])
title('Month with most HWs: CP4FUT')

% for rr=1:length(regions)
%   plot([left_bnds(rr),left_bnds(rr)],[bot_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)   
%   plot([right_bnds(rr),right_bnds(rr)],[bot_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)
%   plot([left_bnds(rr),right_bnds(rr)],[bot_bnds(rr),bot_bnds(rr)],'k-','LineWidth',1.2)
%   plot([left_bnds(rr),right_bnds(rr)],[top_bnds(rr),top_bnds(rr)],'k-','LineWidth',1.2)  
% end 

printname = ['impala_wbmax_HWMId_4kmfc_PDBL_m_with_most_3dayHWs.png'];
eval(['print ''',outdir,'/',printname,''' -dpng -r200']);
close