function [] = func_plot_soil_vert_profile(lon2d, lat2d, hgt,lh, sh, 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'...
               };

for ics = 1:5
figure
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1000 600])
%-----------------------------------------------------
subplot(2,2,1)
surf(lon2d, lat2d, hgt,lh(:,:,ics,season)); 
shading interp; view(20,50); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
caxis([-4. 4.])
colormap (gca, bluewhitered  )
% h = colorbar;
zlim([0 800]);
title('LH')

%-----------------------------------------------------
subplot(2,2,2)
surf(lon2d, lat2d, hgt,sh(:,:,ics,season)); 
shading interp; view(20,50); 
xlabel('Longitude')
ylabel('Latitude')
caxis([-4. 4.])
colormap (gca, bluewhitered  )
% h = colorbar;
zlim([0 800]);
title('SH')

%-----------------------------------------------------


title_name = [char(case_list(ics)) , ' in ', char(season_list(season))];
suptitle(title_name)


print('-depsc', title_name);

end


end

