function [] = func_plot_soil(lon2d, lat2d, hgt,soil_1, sm_label)
% PLOT the soil for different seasons

make_it_tight = true;
subplot = @(m,n,p) subtightplot (m, n, p, [0.07 0.07], [0.07 0.07], [0.1 0.1]);
if ~make_it_tight,  clear subplot;  end

figure
set(gcf,'PaperPositionMode', 'auto')
set(gcf,'Position', [100 100 1000 600])
subplot(2,2,1)
surf(lon2d, lat2d, hgt,soil_1(:,:,1,1)); 
shading interp; view(20,50); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
caxis([0.05 0.4])
colormap (gca, jet  )
% h1 = colorbar('southoutside');
zlim([0 800]);
title('DJF')

subplot(2,2,2)
surf(lon2d, lat2d, hgt,soil_1(:,:,1,2)); 
shading interp; view(20,50); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
caxis([0.05 0.4])
colormap (gca, jet  )
% h1 = colorbar;
zlim([0 800]);
title('MAM')

subplot(2,2,3)
surf(lon2d, lat2d, hgt,soil_1(:,:,1,3)); 
shading interp; view(20,50); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
caxis([0.05 0.4])
colormap (gca, jet  )
% hp4 = get(subplot(2,2,3),'Position')
% h = colorbar('southoutside')
% set(h, 'Position', [hp4(1)  hp4(2)-0.032  hp4(3)*2.1 0.02 ]);
% h.Label.String = [sm_label,' (m^3 m^{-3})'];
zlim([0 800]);
title('JJA')

subplot(2,2,4)
surf(lon2d, lat2d, hgt,soil_1(:,:,1,4)); 
shading interp; view(20,50); 
xlabel('Longitude')
ylabel('Latitude')
zlabel('Surface Elevation (m)')
caxis([0.05 0.4])
colormap (gca, jet  )
hp4 = get(subplot(2,2,4),'Position')
h = colorbar('Position', [hp4(1)+hp4(3)+0.03  hp4(2)  0.01  hp4(2)+hp4(3)*2.1]);
h.Label.String = [sm_label,' (m^3 m^{-3})'];

zlim([0 800]);
title('SON')

end

