clear all
close all

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

run = 'nlcd2011-jan01';

% Get static data of interest.

terrain = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'HGT');
mask = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'LANDMASK');
lon = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'XLONG');
lat = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'XLAT');
flg = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_00:00:00'),'LANDMASK');

temp = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_12:00:00'),'T2')-273;
cloud = ncread(strcat('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/',run,'/wrfout_d03_2011-01-01_12:00:00'),'QCLOUD');
snowh_NLCD = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/output/nlcd2011-dec31/wrfout_d03_2010-12-31_00:00:00','SNOWH');
snowh_USGS = ncread('/uufs/chpc.utah.edu/common/home/horel-group4/model/cfoster/kingspeak/wrf-usgs/WPS/met_em.d03.2010-12-31_00:00:00.nc','SNOWH');


cloud = reshape(cloud(:,:,1:41,1),[474,498,41]);
cloud = mean(cloud,3);

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 = 160; 
lx_max = 240;
ly_min = 290;
ly_max = 405;

figure(1)
title('NLCD Jan 01 Run Temperature (C)')
ax1 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
temp_ax = surfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),temp(lx_min:lx_max,ly_min:ly_max,1));
colorbar
contourcmap('jet',[0:-1:-35])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',1.5)
contourm(lat,lon,terrain,'k','LevelStep',300,'LineWidth',0.5)

figure(2)
title('NLCD Jan 01 Run QCLOUD')
ax1 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
temp_ax = surfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),cloud(lx_min:lx_max,ly_min:ly_max,1));
colorbar
contourcmap('jet',[0:1e-6:10e-5])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',1.5)
contourm(lat,lon,terrain,'k','LevelStep',300,'LineWidth',0.5)

figure(3)
title('NLCD Dec 31 Run SNOWH')
ax1 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
temp_ax = surfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),snowh_NLCD(lx_min:lx_max,ly_min:ly_max,1));
colorbar
contourcmap('cool',[0:0.1:2])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',1.5)
contourm(lat,lon,terrain,'k','LevelStep',300,'LineWidth',0.5)

figure(4)
title('NLCD Jan 01 NAM Snow Initial Fields')
ax1 = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
temp_ax = surfm(lat(lx_min:lx_max,ly_min:ly_max),lon(lx_min:lx_max,ly_min:ly_max),snowh_USGS(lx_min:lx_max,ly_min:ly_max,1));
colorbar
contourcmap('cool',[0:0.1:2])
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',1.5)
contourm(lat,lon,terrain,'k','LevelStep',300,'LineWidth',0.5)

snow_diff = snowh_NLCD(:,:,1)-snowh_USGS;

figure(8)
subplot(1,2,2)
title('Landmask Difference, NLCD 2011 - USGS')
ax = usamap([lat(lx_min,ly_min),lat(lx_max,ly_max)],[lon(lx_min,ly_max),lon(lx_max,ly_min)]);
contourfm(double(lat),double(lon),double(snow_diff),[-0.5:0.05:0.5],'LineStyle','none')
contourcmap('cool',[-0.5:0.05:0.5])
caxis([-0.5 0.5])
colorbar('YTick',[-0.5:0.1:0.5],'FontSize',12)
hold on
contourm(lat,lon,flg,'k','LevelStep',1,'LineWidth',1.5)
contourm(lat,lon,terrain,'k','LevelStep',125,'LineWidth',0.5)