clear;

url='/home/n04d3/zhup/Patricia/p_2106_ntke_1/wrfout_d03_2015-10-23_08:00:00';

lon=hdf5read(url,'HLON');
lat=hdf5read(url,'HLAT');
t=hdf5read(url,'T');
var=hdf5read(url,'Z');
z=0.5*(var(:,:,1:end-1)+var(:,:,2:end));
cwm=hdf5read(url,'CWM');
var=hdf5read(url,'PINT');
pres=0.5*(var(:,:,1:end-1)+var(:,:,2:end));
qv=hdf5read(url,'QVAPOR');

qt=qv+cwm;

xlv=2500000;
cp=1004;
rd=287;

th=t.*(1e5./pres).^(rd/cp);

th_e=th.*(1+xlv./(cp*t).*qv);


subplot(2,2,1);
cwm_h=squeeze(cwm(:,:,30))*1000;
contourf(lon,lat,cwm_h,[0:0.5:5],'edgecolor','none');
caxis([0 5]);
colorbar;
hold on;
%plot(lon(110:150,268),lat(110:150,268),'w');
plot([-105.5,-103],[16.43 16.43],'w');
axis([-104.55,-104.05,16.2,16.7]);
plot(lon(120,268),lat(120,268),'g*');
plot(lon(121,268),lat(121,268),'r*');
plot(lon(122,268),lat(122,268),'k*');
hold off;

z_m=squeeze(z(120,268,1:36))/1000;
th_e1=squeeze(th_e(120,268,1:36));
th_e2=squeeze(th_e(121,268,1:36));
th_e3=squeeze(th_e(122,268,1:36));
qt_1=squeeze(qt(120,268,1:36));
qt_2=squeeze(qt(121,268,1:36));
qt_3=squeeze(qt(122,268,1:36));
th1=squeeze(th(120,268,1:36));
th2=squeeze(th(121,268,1:36));
th3=squeeze(th(122,268,1:36));

subplot(2,2,2);
plot(th_e1,z_m,'g');
hold on;
plot(th_e2,z_m,'r');
plot(th_e3,z_m,'k');
hold off;

subplot(2,2,3);
plot(qt_1*1000,z_m,'g');
hold on;
plot(qt_2*1000,z_m,'r');
plot(qt_3*1000,z_m,'k');
hold off;

subplot(2,2,4);
dqt=(qt_1-qt_2);
m_th=0.5*(th1+th2);
cri=m_th.*dqt./0.5;
dth_e=(th_e1-th_e2);

dqt1=(qt_2-qt_3);
m_th1=0.5*(th2+th3);
cri1=m_th1.*dqt1./0.5;
dth_e1=(th_e1-th_e2);

dqt2=(qt_1-qt_3);
m_th2=0.5*(th1+th3);
cri2=m_th2.*dqt2./0.5;
dth_e2=(th_e1-th_e3);


plot(dth_e-cri,z_m);
hold on;
plot(dth_e1-cri1,z_m,'r');
plot(dth_e2-cri2,z_m,'g');
plot([0 0],[0 6],'--','color',[0.6 0.6 0.6]);
hold off;

save data/en_instability_dat.mat lon lat cwm_h z_m th_e1 th_e2 th_e3 qt_1 qt_2 qt_3 th1 th2 th3 cri dth_e cri1 dth_e1 cri2 dth_e2;
 

