%% Script Description 

clear all 
close all

format shortg

%% Display contents of netcdf file. 

%ncdisp('../ember/wrf-nlcd2011/WPS/met_em_new.d03.2010-12-31_00:00:00.nc')

%% Determine the values and ncread the files into matlab variables.

run_path = '../kingspeak/wrf-nlcd2011/';
domain = 'd03';

% Variables that won't be modified.
lat = ncread(strcat(run_path,'WPS/met_em_new.',domain,'.2010-12-31_00:00:00.nc'),'XLAT_M');
lon = ncread(strcat(run_path,'WPS/met_em_new.',domain,'.2010-12-31_00:00:00.nc'),'XLONG_M');
hgt = ncread(strcat(run_path,'WPS/met_em_new.',domain,'.2010-12-31_00:00:00.nc'),'HGT_M');
flg = ncread(strcat(run_path,'WPS/met_em_new.',domain,'.2010-12-31_00:00:00.nc'),'LANDMASK');

% Variables to modify.
snow  = ncread(strcat(run_path,'WPS/met_em_new.',domain,'.2010-12-31_00:00:00.nc'),'SNOW');
snowh = ncread(strcat(run_path,'WPS/met_em_new.',domain,'.2010-12-31_00:00:00.nc'),'SNOWH');

nam_snowh = snowh;

lon = double(lon(:,:,1));
lat = double(lat(:,:,1));
hgt = double(hgt);

latlim = [min(lat(:)) max(lat(:))];
lonlim = [min(lon(:)) max(lon(:))];

%% Set boundaries for Salt Lake Basin snow modification.

% Domain 3

x_min = 140;    % Boundaries for modification to snow.
x_max = 450;
y_min = 185;
y_max = 470;

% Domain 2

% x_min = 90;    
% x_max = 190;
% y_min = 100;
% y_max = 210;

% Domain 1

% x_min = 75;    
% x_max = 120;
% y_min = 85;
% y_max = 130;

% vx_min = 212;   % Valley domain boundaries (find min height).
% vx_max = 240;
% vy_min = 290;
% vy_max = 322;
% 
% lx_min = 155+13; 
% lx_max = 235;
% ly_min = 285;
% ly_max = 400;

new_snow = snow(x_min:x_max,y_min:y_max);
new_snowh = snowh(x_min:x_max,y_min:y_max);
new_hgt = hgt(x_min:x_max,y_min:y_max);

figure(1)
ax=usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)]);
surfm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),snowh(x_min:x_max,y_min:y_max));
title('NAM Initialization Snow Height (m)')
colorbar('YTick',[0:0.05:1])
contourcmap('jet',[0:0.025:1])
contourm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),new_hgt,'k','LevelStep',200,'LineWidth',0.5);

hgt_min = 1263;    % Valley (lake) floor.
hgt_bench = 1584;  % 5200 ft, location of benches.
hgt_max = 2800;    % Mountain tops.

snow_shore = (7)  * 0.0254;  %0.15;%(8) * 0.0254; % Snow depth in inches * m
snow_bench = (10.5) * 0.0254; %0.25;%(14) * 0.0254;
snow_peaks = (85) * 0.0254 %2.20;%(30) * 0.0254;

slope_valley = (snow_bench-snow_shore)/(hgt_bench-hgt_min);    % Snow increase w/ height from base to benches.
slope_mtn = (snow_peaks-snow_bench)/(hgt_max-hgt_bench);   % Snow increase w/ height from benches to mountaintops.

new_snowh(find(new_hgt >= hgt_min & new_hgt <= hgt_bench)) = (new_hgt(find(new_hgt >= hgt_min & new_hgt <= hgt_bench))-hgt_min) .* slope_valley + snow_shore;
new_snowh(find(new_hgt >= hgt_bench & new_hgt <= hgt_max)) = (new_hgt(find(new_hgt >= hgt_bench & new_hgt <= hgt_max))-hgt_bench) .* slope_mtn + snow_bench;
new_snowh(find(new_hgt >= hgt_max)) = snow_peaks;

new_snowh = new_snowh .* flg(x_min:x_max,y_min:y_max); % Get rid of snow where there's lake.
new_snow = new_snowh * 100 ; % kg/m-3

figure(2)
ax=usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)]);
surfm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),new_snowh);
title('Modified Snow Height Based on Observations (m)')
colorbar('YTick',[0:0.05:1])
contourcmap('jet',[0:0.025:1])
contourm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),new_hgt,'k','LevelStep',200,'LineWidth',0.5);

%% Replace modified values in whole domain.

snowh(x_min:x_max,y_min:y_max) = new_snowh;
snow(x_min:x_max,y_min:y_max) = new_snow;

% figure(3)
% ax=usamap(latlim,lonlim);
% surfm(lat,lon,snowh);
% title('Modified Snow Height (m) Embedded within NAM Domain 3 Intialization Field')
% colorbar('YTick',[0:0.1:2])
% contourcmap('cool',[0:0.05:2])
% contourm(lat,lon,hgt,'k','LevelStep',200,'LineWidth',0.5);
% 
% diff_snowh = snowh - nam_snowh;
% 
% figure(4)
% ax=usamap([lat(x_min,y_min),lat(x_max,y_max)],[lon(x_min,y_max),lon(x_max,y_min)]);
% surfm(lat(x_min:x_max,y_min:y_max),lon(x_min:x_max,y_min:y_max),diff_snowh(x_min:x_max,y_min:y_max));
% title('Difference Between Modified Snow Cover and NAM Initialization Snow Cover')
% colorbar('YTick',[-0.5:0.1:1.5])
% contourcmap('cool',[-0.5:0.05:1.5])
% contourm(lat,lon,hgt,'k','LevelStep',200,'LineWidth',0.5);
% 
% %%
% 
% figure('Position', [100, 100, 900, 475])
% 
% modis = imread('modis_lake.png');
% 
% [X,map] = rgb2ind(modis,1000);
% 
% subplot(1,2,1)
% subimage(X,map)
% set(gca,'XTickLabel','')
% set(gca,'YTickLabel','')
% hold on
% patch([40 5 5 40],[5 5 40 40],'white')
% text(12.5,22.5,'a','FontSize',14)
% 
% subplot(1,2,2), 
% 
% ax=usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
% contourfm(lat,lon,double(diff_snowh),[-0.5:0.05:0.5],'LineStyle','none');
% c = colorbar;
% contourcmap('jet',[-0.5:0.05:0.5])%colormap(flipud(gray))
% caxis([-0.5 0.5])
% set(c,'fontsize',10);
% cpos = c.Position;
% cpos(3) = 0.5*cpos(3);
% c.Position = cpos;
% set(c,'position',[0.85 0.195 0.02 0.655])
% hL = title(c,'m','FontSize',12); 
% 
% contourm(lat,lon,hgt,'k','LevelStep',125,'LineWidth',0.5);
% contourm(lat,lon,double(flg),[1],'k','LineWidth',1);
% 
% fl = figm_label('b','ul',14,0.0425,0.04,12000,12000);
% set(fl,'color','k');
% 
% fig_nudge(ax, -0.05, 0)
% 
% %%
% 
% figure(6)
% ax=usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
% surfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),snowh(lx_min:lx_max,ly_min:ly_max));
% title('Modified Snow Height Based on Observations (m)')
% colorbar('YTick',[0:0.01:0.2])
% contourcmap('cool',[0:0.01:0.2])
% contourm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),hgt(lx_min:lx_max,ly_min:ly_max),'k','LevelStep',100,'LineWidth',0.5);


%% Save modified files in met_em netcdf file.

ncwrite(strcat(run_path,'WPS/met_em.',domain,'.2010-12-31_00:00:00.nc'),'SNOWH',snowh);
ncwrite(strcat(run_path,'WPS/met_em.',domain,'.2010-12-31_00:00:00.nc'),'SNOW',snow);