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 ----
lh    = ncread('seavar_LHSH_06-13.nc', 'seavar_LH');
sh    = ncread('seavar_LHSH_06-13.nc', 'seavar_SH');

lh(:,:,:,5) = squeeze(mean(lh , 4));
sh(:,:,:,5) = squeeze(mean(sh , 4));

lh_diff = func_calc_diff(lh);
sh_diff = func_calc_diff(sh);

%% ----- 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
% ---------------------------
%  LH SH
% ---------------------------

% close all
% for season = 5:5
% func_plot_soil_vert_profile(lon2d, lat2d, hgt,lh_diff,...
%                                               sh_diff,...
%                                               season)
% 
% end

%% IF_LATERAL ALL DIFFERENCE CASES COMBINED
close all;
season = 5;
func_plot_if_lateral_soil_vert_profile_combined(lon2d, lat2d, hgt,lh_diff,...
                                              sh_diff,...
                                              season);

func_saveasnetcdf(lat2d,lon2d,lh_diff,'lh_diff');
func_saveasnetcdf(lat2d,lon2d,sh_diff,'sh_diff');
                                          
%% RELATED TO RESOLUTION IMPACT 
close all
season = 5;
func_plot_soil_res(lon2d, lat2d, hgt,lh_diff,...
                                     sh_diff,...
                                     season)

%% related to resolution impact, but in terms of percentage

[sm1_mean_pct_change] = func_cal_domain_mean(lh, lh_diff);
[sm2_mean_pct_change] = func_cal_domain_mean(sh, sh_diff);

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

%% PLOT PDF OF SOIL MOISTURE CONTENT
% 
% season = 5;
% func_plot_soil_pdf(lh_diff,lh,...
%                    sh_diff,sh,...
%                     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(lh_diff,sh_diff,slope);

%% Calculate and plot the bar plot for each slope bin
func_plot_soil_slope_bar(lh_diff,sh_diff,slope);
print('-depsc', '1km_slope_barplot')
%%
func_plot_soil_slope_bar_p1(lh_diff, sh_diff,slope);
print('-depsc', '1km_slope_barplot_p1')