function [] = func_plot_if_lateral_soil_vert_profile_combined(lon2d, lat2d, hgt,soil_1,soil_2, season)
% Summary of this function goes here
%   
make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.042 0.045], [0.14 0.04], [0.05 0.05]);
if ~make_it_tight,  clear subplot;  end

ann_list = char(97:122);

season_list = {'DJF', 'MAM', 'JJA', 'SON','ANN'};
case_list   = {'both routing - no routing',...
               'surface routing - no routing',...
               'subsurface routing - no routing', ...
               'both routing - subsurface routing', ...
               'both routing - surface routing'...
               };

figure
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1200 400])

iplot = 1;
for ics = 1:5

%-----------------------------------------------------
subplot(2,5,iplot)
surf(lon2d, lat2d, hgt,soil_1(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
% ylabel('Latitude')
% zlabel('Surface Elevation (m)')
caxis([-8 8])
colormap (gca, bluewhitered(20)  )
% h = colorbar;
zlim([0 800]);
if iplot ==1; zlabel('LH '); end
title_name = [char(case_list(ics))];
title(title_name)

ann = [ann_list(iplot), ')'];
text(-100, 35, 1200, ann);
%-----------------------------------------------------
subplot(2,5,iplot+5)
surf(lon2d, lat2d, hgt,soil_2(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
% ylabel('Latitude')
caxis([-8 8])
colormap (gca, bluewhitered(20)  )
% h = colorbar;
zlim([0 800]);
if iplot ==1; zlabel('SH '); end
ann = [ann_list(iplot+5), ')'];
text(-100, 35, 1200, ann);


% suptitle(title_name)

iplot = iplot +1;

end

% hp3 = get(subplot(2,5,2),'Position')
% h = colorbar('southoutside');
% set(h, 'Position', [hp3(1) + 0.1  hp3(2)-0.06  hp3(3)*1.8 0.01 ]);
% h.Label.String = ['LH diff. (W m^{-2})'];

hp3 = get(subplot(2,5,7),'Position')
h = colorbar('southoutside');
set(h, 'Position', [hp3(1) + 0.1  hp3(2)-0.06  hp3(3)*1.8 0.01 ]);
h.Label.String = ['SH or LH diff. (W m^{-2})'];

title_name = ['1km_if_lateral_all'];
print('-depsc', title_name);

%%

figure
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1200 400])

iplot = 1;
for ics = 1:3

%-----------------------------------------------------
subplot(2,3,iplot)
surf(lon2d, lat2d, hgt,soil_1(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
% ylabel('Latitude')
% zlabel('Surface Elevation (m)')
caxis([-8 8])
colormap (gca, bluewhitered(20)  )
% h = colorbar;
zlim([0 800]);
if iplot ==1; zlabel('LH '); end
title_name = [char(case_list(ics))];
title(title_name)

ann = [ann_list(iplot), ')'];
text(-100, 35, 1200, ann);
%-----------------------------------------------------
subplot(2,3,iplot+3)
surf(lon2d, lat2d, hgt,soil_2(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
% ylabel('Latitude')
caxis([-8 8])
colormap (gca, bluewhitered(20)  )
% h = colorbar;
zlim([0 800]);
if iplot ==1; zlabel('SH '); end
ann = [ann_list(iplot+3), ')'];
text(-100, 35, 1200, ann);


% suptitle(title_name)

iplot = iplot +1;

end

% hp3 = get(subplot(2,5,2),'Position')
% h = colorbar('southoutside');
% set(h, 'Position', [hp3(1) + 0.1  hp3(2)-0.06  hp3(3)*1.8 0.01 ]);
% h.Label.String = ['LH diff. (W m^{-2})'];

hp3 = get(subplot(2,3,4),'Position')
h = colorbar('southoutside');
set(h, 'Position', [hp3(1) + 0.16  hp3(2)-0.06  hp3(3)*1.8 0.01 ]);
h.Label.String = ['SH or LH diff. (W m^{-2})'];

title_name = ['1km_if_lateral'];
print('-depsc', title_name);



%% Calculate the interactive effect... 

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.003 0.065], [0.16 0.045], [0.08 0.05]);
if ~make_it_tight,  clear subplot;  end


int_list = {'Interaction between surface and subsurface flow', 'Interaction from subsurface lateral flow'};
figure
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 800 400])

iplot = 1;
for ics = 12:12

%-----------------------------------------------------
subplot(1,2,iplot)
surf(lon2d, lat2d, hgt,soil_1(:,:,ics,season)); 
shading interp; view(20,50); 
caxis([-8 8])
colormap (gca, bluewhitered(20)  )
% h = colorbar;
zlim([0 800]);
% if iplot ==1; zlabel('SM1'); end

title_name = [char(int_list(iplot))];
% title(title_name)

ann = [ann_list(iplot), ') ','LH'];
text(-100, 35, 1200, ann);

%-----------------------------------------------------
subplot(1,2,iplot+1)
surf(lon2d, lat2d, hgt,soil_2(:,:,ics,season)); 
shading interp; view(20,50); 

caxis([-8 8])
colormap (gca, bluewhitered(20)  )

zlim([0 800]);
% if iplot ==1; zlabel('SM2'); end
ann = [ann_list(iplot+1), ') ','SH'];
text(-100, 35, 1200, ann);



iplot = iplot +1;

end

hp3 = get(subplot(1,2,1),'Position')
h = colorbar('southoutside');
set(h, 'Position', [hp3(1) + 0.125  hp3(2)-0.07  hp3(3)*1.25 0.01 ]);
h.Label.String = ['SH or LH diff. (W m^{-2})'];
suptitle(title_name)


title_name = ['1km_if_lateral_interactive'];
print('-depsc', title_name);


end

