clear all
close all

set(0,'defaultfigureposition',[500 500 500 420])

% Get static data of interest.

terrain = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/usgs-dec31/met_em.d03.2010-12-31_00:00:00.nc','HGT_M');
mask = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/usgs-dec31/met_em.d03.2010-12-31_00:00:00.nc','LANDMASK');
lon = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/usgs-dec31/met_em.d03.2010-12-31_00:00:00.nc','XLONG_M');
lat = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/usgs-dec31/met_em.d03.2010-12-31_00:00:00.nc','XLAT_M');
flg = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/usgs-dec31/met_em.d03.2010-12-31_00:00:00.nc','LANDMASK');

snow = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/usgs-dec31/met_em.d03.2010-12-31_00:00:00.nc','SNOW');
snowh_31 = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/usgs-dec31/wrfout_d03_2010-12-31_00:00:00','SNOWH');
snowh_01 = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/usgs-dec31/wrfout_d03_2011-01-01_00:00:00','SNOWH');


lat = double(lat(:,:,1));
lon = double(lon(:,:,1));
terrain = double(terrain(:,:,1));
flg = double(flg(:,:,1));

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

states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);

% lx_min = 155; 
% lx_max = 235;
% ly_min = 285;
% ly_max = 400;

lx_min = 150;   
lx_max = 425;
ly_min = 180;
ly_max = 425;

snow_diff = snowh_01(:,:,1) - snowh_31(:,:,1);

figure(1)
title('NAM Model Snow Height , Jan 01 - Dec 31 00 Z (m)')
ax1 = usamap(latlim,lonlim);
temp_ax = surfm(lat,lon,snow_diff);
geoshow(ax1, states, 'FaceColor', 'none');
colorbar('YTick',[0:0.025:0.5])
contourcmap('jet',[0:0.05:0.5])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',1.5)
contourm(lat,lon,terrain,'k','LevelStep',300,'LineWidth',0.5)

figure(2)

title('NAM Initialize Snow Height (m)')
ax1 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
temp_ax = contourfm(double(lat),double(lon),double(snowh_31(:,:,1)),[0:0.05:1],'LineStyle','none');
colorbar
contourcmap('cool',[0:0.05:1])
caxis([0 1])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',1.5)
contourm(lat,lon,terrain,'k','LevelStep',300,'LineWidth',0.5)