clc 
clear 
close all 

%% PLOT the 3D soil moisture on surface. 
addpath('../ncfiles_1km/')
addpath('/Users/yang560/Dropbox/PNNL/Matlab_functions/utilities/')

soilt_gssurgo = ncread('geo_em.d04.nc.gssurgo', 'SCT_DOM');
lat2d     = ncread('geo_em.d04.nc.gssurgo', 'XLAT_M');
lon2d     = ncread('geo_em.d04.nc.gssurgo','XLONG_M');
hgt       = ncread('geo_em.d04.nc.gssurgo', 'HGT_M');

lat2d = double (lat2d);
lon2d = double (lon2d);
hgt   = double (hgt);

%---- add  soil moisture file ----
soil_1 = ncread('seavar_SM1_06-13.nc','seavar_SOIL_M1');
soil_2 = ncread('seavar_SM2_06-13.nc','seavar_SOIL_M2');
soil_3 = ncread('seavar_SM3_06-13.nc','seavar_SOIL_M3');
soil_4 = ncread('seavar_SM4_06-13.nc','seavar_SOIL_M4');

soil_1(:,:,:,5) = squeeze(mean(soil_1 , 4));
soil_2(:,:,:,5) = squeeze(mean(soil_2 , 4));
soil_3(:,:,:,5) = squeeze(mean(soil_3 , 4));
soil_4(:,:,:,5) = squeeze(mean(soil_4 , 4));

soil_1_diff = func_calc_diff(soil_1);
soil_2_diff = func_calc_diff(soil_2);
soil_3_diff = func_calc_diff(soil_3);
soil_4_diff = func_calc_diff(soil_4);

%% ----- PLOT SECTION -----
% close all
% func_plot_soil(lon2d, lat2d, hgt,soil_1,'SM1'); print('-depsc', 'SM1_4s') 
% func_plot_soil(lon2d, lat2d, hgt,soil_2,'SM2'); print('-depsc', 'SM2_4s')
% func_plot_soil(lon2d, lat2d, hgt,soil_3,'SM3'); print('-depsc', 'SM3_4s')
% func_plot_soil(lon2d, lat2d, hgt,soil_4,'SM4'); print('-depsc', 'SM4_4s')

%% ----- PLOT FOR DIFFERENT DIFFERENCE SCENARIOS -----
%    W/O surface/subsurface routing schemes. 

% close all
% season=1;
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_1_diff,'SM1',season);print('-depsc', 'SM1_diff_DJF')
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_2_diff,'SM2',season);print('-depsc', 'SM2_diff_DJF')
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_3_diff,'SM3',season);print('-depsc', 'SM3_diff_DJF')
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_4_diff,'SM4',season);print('-depsc', 'SM4_diff_DJF')
% season=3;
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_1_diff,'SM1',season);print('-depsc', 'SM1_diff_JJA')
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_2_diff,'SM2',season);print('-depsc', 'SM2_diff_JJA')
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_3_diff,'SM3',season);print('-depsc', 'SM3_diff_JJA')
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_4_diff,'SM4',season);print('-depsc', 'SM4_diff_JJA')
% 
% season=5;
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_1_diff,'SM1',season);print('-depsc', 'SM1_diff_ANN')
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_2_diff,'SM2',season);print('-depsc', 'SM2_diff_ANN')
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_3_diff,'SM3',season);print('-depsc', 'SM3_diff_ANN')
% func_plot_if_lateral(lon2d, lat2d, hgt,soil_4_diff,'SM4',season);print('-depsc', 'SM4_diff_ANN')

%% RELATED TO IF OR NOT LATERAL FLOW
% Individual plot for different scheme
% ---------
%  sm1 sm2 
%  sm3 sm4
% ---------

close all
for season = 5:5
func_plot_soil_vert_profile(lon2d, lat2d, hgt,soil_1_diff,...
                                              soil_2_diff,...
                                              soil_3_diff,...
                                              soil_4_diff,season)

end


%% SAVE AS SETCDF...

func_saveasnetcdf(lat2d,lon2d,soil_1_diff,'soil_1_diff');
func_saveasnetcdf(lat2d,lon2d,soil_2_diff,'soil_2_diff');
func_saveasnetcdf(lat2d,lon2d,soil_3_diff,'soil_3_diff');
func_saveasnetcdf(lat2d,lon2d,soil_4_diff,'soil_4_diff');

%% IF_LATERAL ALL DIFFERENCE CASES COMBINED
close all;
for season = 5:5
func_plot_if_lateral_soil_vert_profile_combined(lon2d, lat2d, hgt,soil_1_diff,...
                                              soil_2_diff,...
                                              soil_3_diff,...
                                              soil_4_diff,season)

end
%% RELATED TO RESOLUTION IMPACT 
close all
for season = 5:5
func_plot_soil_res(lon2d, lat2d, hgt,soil_1_diff,...
                                     soil_2_diff,...
                                     soil_3_diff,...
                                     soil_4_diff,season)
end
%% related to resolution impact, but in terms of percentage

[sm1_mean_pct_change] = func_cal_domain_mean(soil_1, soil_1_diff);
[sm2_mean_pct_change] = func_cal_domain_mean(soil_2, soil_2_diff);
[sm3_mean_pct_change] = func_cal_domain_mean(soil_3, soil_3_diff);
[sm4_mean_pct_change] = func_cal_domain_mean(soil_4, soil_4_diff);

func_plot_soil_pct_change_res(lon2d, lat2d, hgt,sm1_mean_pct_change,sm2_mean_pct_change,sm3_mean_pct_change,sm4_mean_pct_change);


%% PLOT PDF OF SOIL MOISTURE CONTENT
% 
% season = 5;
% func_plot_soil_pdf(soil_1_diff,soil_1,...
%                    soil_2_diff,soil_2,...
%                    soil_3_diff,soil_3,... 
%                    soil_4_diff,soil_4,season)
% 

%% Calculate and plot the slope/gradient for each slope bins
close all
[aspect,slope,gradN,gradE] = gradientm(lat2d,lon2d,hgt);

func_plot_soil_scatter(soil_1_diff,soil_2_diff,soil_3_diff,soil_4_diff,slope);

gradT = (gradN.^2 + gradE.^2).^.5;

close all
figure
surf(lon2d, lat2d, slope); view(2); shading interp;colorbar
slope90 = prctile(slope, 95); 
slope_view = slope;
slope_view(slope_view<slope90)=nan;

figure
surf(lon2d, lat2d, slope_view); view(2); shading interp;colorbar



%% Check where the vertical profile of changes in subsurface impact is...
% variables needed, soil_1_diff soil_4_diff

s1_diff = soil_1_diff(:,:,3,5);   % subsurface impact  + annual 
s4_diff = soil_4_diff(:,:,3,5); 

var1_temp = aspect; 
var1_temp(s4_diff>-0.01  ) = nan; 

var2_temp = slope; 
var2_temp(s4_diff>-0.01  ) = nan; 

var3_temp = hgt; 
var3_temp(s4_diff>-0.01  ) = nan; 

figure
subplot(331)
surf(lon2d, lat2d, var1_temp); view(2); shading interp;colorbar;caxis([1 360])
title ('aspect')
subplot(334)
surf(lon2d, lat2d, var2_temp); view(2); shading interp;colorbar;caxis([0 3])
title('slope')
subplot(337)
surf(lon2d, lat2d, var3_temp); view(2); shading interp;colorbar;caxis([200 650])
title('elevation')

subplot(332)
surf(lon2d, lat2d, aspect); view(2); shading interp;colorbar;caxis([1 360])
title ('aspect')
subplot(335)
surf(lon2d, lat2d, slope); view(2); shading interp;colorbar;caxis([0 3])
title('slope')
subplot(338)
surf(lon2d, lat2d, hgt); view(2); shading interp;colorbar;caxis([200 650])
title('elevation')


var1_temp = aspect; 
var1_temp(s4_diff<=0  ) = nan; 

var2_temp = slope; 
var2_temp(s4_diff<=0  ) = nan; 

var3_temp = hgt; 
var3_temp(s4_diff<=0  ) = nan; 

subplot(333)
surf(lon2d, lat2d, var1_temp); view(2); shading interp;colorbar;caxis([1 360])
title ('aspect')
subplot(336)
surf(lon2d, lat2d, var2_temp); view(2); shading interp;colorbar;caxis([0 3])
title('slope')
subplot(339)
surf(lon2d, lat2d, var3_temp); view(2); shading interp;colorbar;caxis([200 650])
title('elevation')




% figure
% histogram(slope_temp)
% figure;
% histogram(aspect)

%% Calculate and plot the bar plot for each slope bin
func_plot_soil_slope_bar(soil_1_diff,soil_2_diff,soil_3_diff,soil_4_diff,slope);
print('-depsc', '1km_slope_barplot')
%%
func_plot_soil_slope_bar_p1(soil_1_diff,soil_2_diff,soil_3_diff,soil_4_diff,slope);
print('-depsc', '1km_slope_barplot_p1')