function [] = func_plot_soil_res(lon2d, lat2d, hgt,soil_1,soil_2,soil_3,soil_4, season)
% Summary of this function goes here
%   
make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.045 0.09], [0.12 0.025], [0.095 0.05]);
if ~make_it_tight,  clear subplot;  end

season_list = {'DJF', 'MAM', 'JJA', 'SON','ANN'};
ann_list = char(97:122);

iplot = 1;

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

for ics = 4:5
%-----------------------------------------------------
subplot(4,2,iplot)
surf(lon2d, lat2d, hgt,soil_1(:,:,ics,season)); 
shading interp; view(20,50); 
if iplot ==1; zlabel('SM1'); end
caxis([-0.025 0.025])
colormap (gca, bluewhitered  )
% h = colorbar;
zlim([0 800]);
if iplot ==1; title('routing\_250m - routing\_500m'); end
if iplot ==2; title('routing\_500m - routing\_1000m'); end
ann = [ann_list(iplot), ')'];
text(-99, 35, 1200, ann);
%-----------------------------------------------------
subplot(4,2,iplot+2)
surf(lon2d, lat2d, hgt,soil_2(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
if iplot ==1; zlabel('SM2'); end
caxis([-0.025 0.025])
colormap (gca, bluewhitered  )
% h = colorbar;
zlim([0 800]);
ann = [ann_list(iplot+2), ')'];
text(-99, 35, 1200, ann);
%-----------------------------------------------------
subplot(4,2,iplot+4)
surf(lon2d, lat2d, hgt,soil_3(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
if iplot ==1; zlabel('SM3'); end
% zlabel('Surface Elevation (m)')
caxis([-0.025 0.025])
colormap (gca, bluewhitered  )
% h = colorbar;
zlim([0 800]);
ann = [ann_list(iplot+4), ')'];
text(-99, 35, 1200, ann);
%-----------------------------------------------------
subplot(4,2,iplot+6)
surf(lon2d, lat2d, hgt,soil_4(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
if iplot ==1; zlabel('SM4'); end
caxis([-0.025 0.025])
colormap (gca, bluewhitered  )

zlim([0 800]);
ann = [ann_list(iplot+6), ')'];
text(-99, 35, 1200, ann);

iplot = iplot + 1;


end

hp3 = get(subplot(4,2,7),'Position')
h = colorbar('southoutside')
set(h, 'Position', [hp3(1) + 0.05  hp3(2)-0.07  hp3(3)*1.8 0.01 ]);
h.Label.String = ['SM diff. (m^3 m^{-3})'];

print('-dpdf', 'soil_res_cumulative', '-r300');

%% 
iplot = 1;

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

for ics = 2:3
%-----------------------------------------------------
subplot(4,2,iplot)
surf(lon2d, lat2d, hgt,soil_1(:,:,ics,season)); 
shading interp; view(20,50); 
if iplot ==1; zlabel('SM1'); end
caxis([-0.025 0.025])
colormap (gca, bluewhitered  )
% h = colorbar;
zlim([0 800]);
if iplot ==1; title('routing\_250m - routing\_500m'); end
if iplot ==2; title('routing\_250m - routing\_1000m'); end
ann = [ann_list(iplot), ')'];
text(-99, 35, 1200, ann);
%-----------------------------------------------------
subplot(4,2,iplot+2)
surf(lon2d, lat2d, hgt,soil_2(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
if iplot ==1; zlabel('SM2'); end
caxis([-0.025 0.025])
colormap (gca, bluewhitered  )
% h = colorbar;
zlim([0 800]);
% title('SM2')
ann = [ann_list(iplot+2), ')'];
text(-99, 35, 1200, ann);
%-----------------------------------------------------
subplot(4,2,iplot+4)
surf(lon2d, lat2d, hgt,soil_3(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
if iplot ==1; zlabel('SM3'); end
% zlabel('Surface Elevation (m)')
caxis([-0.025 0.025])
colormap (gca, bluewhitered  )
% h = colorbar;
zlim([0 800]);
ann = [ann_list(iplot+4), ')'];
text(-99, 35, 1200, ann);
%-----------------------------------------------------
subplot(4,2,iplot+6)
surf(lon2d, lat2d, hgt,soil_4(:,:,ics,season)); 
shading interp; view(20,50); 
% xlabel('Longitude')
if iplot ==1; zlabel('SM4'); end
caxis([-0.025 0.025])
colormap (gca, bluewhitered  )

zlim([0 800]);
ann = [ann_list(iplot+6), ')'];
text(-99, 35, 1200, ann);

iplot = iplot + 1;


end

hp3 = get(subplot(4,2,7),'Position')
h = colorbar('southoutside')
set(h, 'Position', [hp3(1) + 0.05  hp3(2)-0.07  hp3(3)*1.8 0.01 ]);
h.Label.String = ['SM diff. (m^3 m^{-3})'];

% print('-dpdf', 'soil_res', '-r300');
print('-dpng', 'soil_res')

end

