function [] = func_plot_soil_pdf(soil_1,s1, soil_2,s2, season)
% Summary of this function goes here
%   
make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.055 0.075], [0.12 0.025], [0.05 0.05]);
if ~make_it_tight,  clear subplot;  end

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(11)
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1000 600])
x = [-0.25:0.0001:0.25];

%-----------------------------------------------------
subplot(2,2,1)
for ics = 1:5
sdata = soil_1(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1.); hold on; grid on; 
end
legend('both routing - no routing','surface routing - no routing', 'subsurface routing - no routing', 'both routing - subsurface routing','both routing - surface routing');
xlim([-0.025 0.025])
ylim([0 0.025])

title_left('a) SM1')

%-----------------------------------------------------
subplot(2,2,2)
for ics = 1:5
sdata = soil_2(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1); hold on; grid on; 
end
xlim([-4 4])
ylim([0 0.25])
title_left('a) SM2')

%--------------------
% print('-depsc', '1km_if_lateral_sm_diff_distribution')

%%
figure
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1000 600])
x = [-0.25:0.0001:0.25];

%-----------------------------------------------------
subplot(2,2,1)
for ics = 6:8
sdata = soil_1(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1.); hold on; grid on; 
end
legend('routing 100m - routing 250m','routing 100m - routing 500m', 'routing 100m - routing 1000m');
xlim([-4 4])
ylim([0 0.12])

title_left('a) SM1')

%-----------------------------------------------------
subplot(2,2,2)
for ics = 6:8
sdata = soil_2(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1); hold on; grid on; 
end
xlim([-4 4])
ylim([0 0.12])
title_left('a) SM2')


% print('-depsc', '1km_routing_res_sm_diff_distribution')


%%
figure
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1000 600])
x = [0:0.001:0.6];

%-----------------------------------------------------
subplot(2,2,1)
for ics = 1:4
sdata = s1(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1.); hold on; grid on; 
end
legend('no routing','surface routing only', 'subsurface routing only', 'both routing ');
xlim([0 40])
ylim([0 0.03])

title_left('a) SM1')

%-----------------------------------------------------
subplot(2,2,2)
for ics = 1:4
sdata = s2(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1); hold on; grid on; 
end
xlim([0 40])
ylim([0 0.03])
title_left('a) SM2')

% print('-depsc', '1km_if_lateral_sm_distribution')


%%
figure
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1000 600])
x = [0:0.001:0.6];

%-----------------------------------------------------
subplot(2,2,1)
for ics = 5:8
sdata = s1(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1.); hold on; grid on; 
end
legend('routing 100m','routing 250m', 'routing 500m', 'routing 1000');
xlim([0 40])
ylim([0 0.03])

title_left('a) SM1')

%-----------------------------------------------------
subplot(2,2,2)
for ics = 5:8
sdata = s2(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1); hold on; grid on; 
end
xlim([0 40])
ylim([0 0.03])
title_left('a) SM2')


% print('-depsc', '1km_routing_res_sm_distribution')

end

