function [] = func_plot_soil_pdf(soil_1,s1, soil_2,s2, soil_3,s3, soil_4,s4, 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'};


%%
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 = 2:4
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([-0.012 0.012])
ylim([0 0.12])

title_left('a) SM1')

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

%-----------------------------------------------------
subplot(2,2,3)
for ics = 2:4
sdata = soil_3(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1); hold on; grid on; 
end
xlim([-0.012 0.012])
ylim([0 0.12])
xlabel('SM diff (m^3 m^{-3})')
title_left('c) SM3')

%-----------------------------------------------------
subplot(2,2,4)
for ics = 2:4
sdata = soil_4(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1); hold on; grid on; 
end
xlim([-0.012 0.012])
ylim([0 0.12])
xlabel('SM diff (m^3 m^{-3})')
title_left('d) SM4')

print('-depsc', '4km_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 = 2:5
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 1000m');
xlim([0 0.5])
ylim([0 0.03])

title_left('a) SM1')

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

%-----------------------------------------------------
subplot(2,2,3)
for ics = 2:5
sdata = s3(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1); hold on; grid on; 
end
xlim([0 0.5])
ylim([0 0.03])
xlabel('SM (m^3 m^{-3})')
title_left('c) SM3')

%-----------------------------------------------------
subplot(2,2,4)
for ics = 2:5
sdata = s4(:,:,ics,season);
[p,xc] = hist(sdata(:),x); plot(xc,p/sum(p), 'LineWidth', 1); hold on; grid on; 
end
xlim([0 0.5])
ylim([0 0.03])
xlabel('SM (m^3 m^{-3})')
title_left('d) SM4')
print('-depsc', '4km_routing_res_sm_distribution')

end

