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

%  % ---- difference field ----
%   1: routing_100m - no_routing
%   2: routing_100m - routing_250m
%   3: routing_100m - routing_500m
%   4: routing_100m - routing_1000m
%   5:  routing_100m - routing_250m
%   6: routing_250m - routing_500m
%   7: routing_500m - routing_1000m
soil_diff(:,:,1,:)  = soil(:,:,2,:) -  soil(:,:,1,:);
soil_diff(:,:,2,:)  = soil(:,:,2,:) -  soil(:,:,3,:);
soil_diff(:,:,3,:)  = soil(:,:,2,:) -  soil(:,:,4,:);
soil_diff(:,:,4,:)  = soil(:,:,2,:) -  soil(:,:,5,:);

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

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

end

