ncid = netcdf.open('heat_retrieval.nc','NC_NOWRITE');
xax = netcdf.getVar(ncid,0);
yax = netcdf.getVar(ncid,1);
lhr = netcdf.getVar(ncid,4);

lhr = reshape(lhr,60*60*10*9,1);
lhr(lhr == 0.0) = NaN;
mean1 = nanmean(lhr)

ind = lhr > -0.25 & lhr < 0.25 | lhr > 300.0 | lhr < -150.0;
%ind = find(lhr == 0.0 | lhr > 300.0 | lhr < -150.0);
lhr(ind) = NaN;

lhrcold = lhr(lhr < 0.0);   %latent cooling
lhrwarm = lhr(lhr > 0.0);   %latent heating

ncid = netcdf.open('heat_retrieval_saturated.nc','NC_NOWRITE');
lhs = netcdf.getVar(ncid,4);

lhs = reshape(lhs,60*60*10*9,1);
lhs(lhs == 0.0) = NaN;
mean3 = nanmean(lhs)

return;

%ind = find(lhs >= -0.5 & lhs <= 0.5 | abs(lhs) > 300.0);
%lhs(ind) = NaN;

ncid = netcdf.open('model_heat_steps_2-10.nc','NC_NOWRITE');
lhm = netcdf.getVar(ncid,4);

%avgf = fspecial('average');
%lhm2 = imfilter(lhm,avgf,'replicate');
%ind = find(isnan(lhm2));
%lhm2(ind) = lhm(ind);   %this step fixes boundary problem when filtering

lhm = reshape(lhm,60*60*10*9,1);
lhm(lhm == 0.0) = NaN;
mean2 = nanmean(lhm);
%lhm2 = reshape(lhm2,60*60*10*9,1);

ind = lhm > -0.25 & lhm < 0.25 | lhm > 300.0 | lhm < -150.0;
%ind = find(lhm == 0.0 | lhm > 300.0 | lhm < -150.0);
lhm(ind) = NaN;

%ind = find(lhm2 >= -0.5 & lhm2 <= 0.5 | abs(lhm2) > 300.0);
%lhm2(ind) = NaN;

lhmcold = lhm(lhm < 0.0);   %latent cooling
lhmwarm = lhm(lhm > 0.0);   %latent heating

xx1 = -150:10:0;
xx2 = 0:10:300;

subplot(2,2,1)
hist(lhrcold,xx1);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','b','EdgeColor','w')
%bar(xx1,n1,'blue')
xlim([-150 5])
ylim([0 6.5*10^4])
set(gca,'XTick',[-150:25:0])
xlabel('Latent Cooling Rate (K/h)')
ylabel('Count')


subplot(2,2,2)
hist(lhrwarm,xx2);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','r','EdgeColor','w')
%%bar(xx2,n2,'red')
xlim([-5 300])
ylim([0 3*10^4])
set(gca,'XTick',[0:25:300])
xlabel('Latent Heating Rate (K/h)')
ylabel('Count')


text(-100,3.5*10^4,['Mean = ',num2str(mean1),' K/h'])

subplot(2,2,3)
hist(lhmcold,xx1);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','b','EdgeColor','w')
xlim([-150 5])
ylim([0 6.5*10^4])
set(gca,'XTick',[-150:25:0])
xlabel('Latent Cooling Rate (K/h)')
ylabel('Count')

subplot(2,2,4)
hist(lhmwarm,xx2);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','r','EdgeColor','w')
xlim([-5 300])
ylim([0 3*10^4])
set(gca,'XTick',[0:25:300])
xlabel('Latent Heating Rate (K/h)')
ylabel('Count')

text(-100,3.5*10^4,['Mean = ',num2str(mean2),' K/h'])
