function [soil_diff] = func_calc_diff(soil)
% ========== Summary of this function goes here ==========
% ----  Origional data description soil ----
%   1: no routing 
%   2: surface routing 
%   3: subsurface routing 
%   4: both routing 
%   5: both routing @ 100m  + channel + lake
%   6: both routing @ 250m  + channel + lake
%   7: both routing @ 500m  + channel + lake
%   8: both routing @ 1000m + channel + lake


% ----  Difference data description soil ----
%   1: both routing    - no routing -> total routing impact 
%   2: surface routing - no routing -> merely surface routing impact 
%   3: subsurface routing - no routing -> merely subsurface routing impact 
%   4: both routing - subsurface routing -> surface routing impact 
%   5: both routing - surface routing    -> subsurface routing impact
%   6: routing_100m - routing_250m
%   7: routing_100m - routing_500m
%   8: routing_100m - routing_1000m
%   9:  routing_100m - routing_250m
%   10: routing_250m - routing_500m
%   11: routing_500m - routing_1000m


soil_diff(:,:,1,:) = soil(:,:,4,:) -  soil(:,:,1,:);
soil_diff(:,:,2,:) = soil(:,:,2,:) -  soil(:,:,1,:);
soil_diff(:,:,3,:) = soil(:,:,3,:) -  soil(:,:,1,:);
soil_diff(:,:,4,:) = soil(:,:,4,:) -  soil(:,:,3,:);
soil_diff(:,:,5,:) = soil(:,:,4,:) -  soil(:,:,2,:);
soil_diff(:,:,6,:) = soil(:,:,5,:) -  soil(:,:,6,:);
soil_diff(:,:,7,:) = soil(:,:,5,:) -  soil(:,:,7,:);
soil_diff(:,:,8,:) = soil(:,:,5,:) -  soil(:,:,8,:);

soil_diff(:,:,9,:)  = soil(:,:,5,:) -  soil(:,:,6,:);
soil_diff(:,:,10,:) = soil(:,:,6,:) -  soil(:,:,7,:);
soil_diff(:,:,11,:) = soil(:,:,7,:) -  soil(:,:,8,:);

% ---- Calculate the surface/subsurface interactive effect
soil_diff(:,:,12,:)  = soil_diff(:,:,4,:)  -  soil_diff(:,:,2,:);
soil_diff(:,:,13,:)  = soil_diff(:,:,5,:)  -  soil_diff(:,:,3,:);
soil_diff(:,:,14,:)  = soil_diff(:,:,13,:) -  soil_diff(:,:,12,:);



soil_diff(:,:,:,5) = squeeze(mean(soil_diff , 4)); % obtain the annual value

end

