clear;

load -ASCII coastline;

directory='/mnt/lfs4/HFIP/hfip-fiu/zhup/output/HAFS_B/track_base/';

file_name='natl00l.2018100718.trak.hafs.atcfunix.all';
url_1(1:length(directory))=directory;
url_1(1+length(directory):length(directory)+length(file_name))=file_name;


fid=fopen(url_1,'r');
tmp=fgets(fid);

nm=1;

while length(tmp) > 100
data(1,nm)=str2double(tmp(5:6));
data(2,nm)=str2double(tmp(31:33));
data(3,nm)=str2double(tmp(36:38));
data(4,nm)=str2double(tmp(43:45));
data(5,nm)=str2double(tmp(49:51));
data(6,nm)=str2double(tmp(54:57));
nm=nm+1;
tmp=fgets(fid);
end % while

nm=1;
var=data(1,:);
while isempty(var) == 0
sno(nm)=var(1);
var=var(find(var~=sno(nm)));
nm=nm+1;
end

n1=1;
for n=1:length(data)
if data(1,n) == sno(2)
track(:,n1)=data(:,n);
n1=n1+1;
end
end

n1=1;
track_1(:,n1)=track(:,1);
ch_time=track_1(2,n1);
for n=2:length(track)
if track(2,n) ~= ch_time
n1=n1+1;
track_1(:,n1)=track(:,n);
ch_time=track_1(:,n1);
end
end

nfile=57;

for n=1:nfile
file_o1(n,:)='/mnt/lfs4/HFIP/hfip-fiu/zhup/output/HAFS_B/base/base_2018100718/';
end
fid=fopen('infile');
file2=fscanf(fid,'%s',[50 inf]);
fclose(fid);
file2=file2';
url(:,1:length(file_o1))=file_o1;
url(:,1+length(file_o1):length(file_o1)+50)=file2(1:nfile,:);

%nn=25;  %---Oct. 10, 18:00;
%nn=24;  %---Oct. 10, 15:00;
%nn=23;  %---Oct. 10, 12:00;
nn=22;  %---Oct. 10, 09:00;
%nn=21;  %---Oct. 10, 06:00;
%nn=20;  %---Oct. 10, 03:00;


c_lat=track_1(3,nn)/10;
c_lon=-track_1(4,nn)/10;


grib=ncgeodataset(url(nn,:));
variables = grib.variables;

lon1 = grib.geovariable('lon');
lat1 = grib.geovariable('lat');
[LON0,LAT0]=meshgrid(lon1,lat1);

lon0(:)=lon1(:,:)';
lat0(:)=lat1(:,:)';

[aaa,nx]=min(abs(lon0-c_lon));
[aaa,ny]=min(abs(lat0-c_lat));

nd=150;
%nd=75;

b_LON=LON0(ny-nd:ny+nd,nx-nd:nx+nd);
b_LAT=LAT0(ny-nd:ny+nd,nx-nd:nx+nd);

p1=grib.geovariable('isobaric');
p0=p1(:);

t1 = grib.geovariable('Temperature_isobaric');
t0(:,:,:)=t1(1,:,:,:);
b_t=t0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qv1 = grib.geovariable('Specific_humidity_isobaric');
qv0(:,:,:)=qv1(1,:,:,:);
b_qv=qv0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qc1 = grib.geovariable('Cloud_mixing_ratio_isobaric');
qc0(:,:,:)=qc1(1,:,:,:);
b_qc=qc0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qi1 = grib.geovariable('Ice_water_mixing_ratio_isobaric');
qi0(:,:,:)=qi1(1,:,:,:);
b_qi=qi0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qr1 = grib.geovariable('Rain_mixing_ratio_isobaric');
qr0(:,:,:)=qr1(1,:,:,:);
b_qr=qr0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qs1 = grib.geovariable('Snow_mixing_ratio_isobaric');
qs0(:,:,:)=qs1(1,:,:,:);
b_qs=qs0(:,ny-nd:ny+nd,nx-nd:nx+nd);


refl1=grib.geovariable('Reflectivity_isobaric');
refl0(:,:)=mean(refl1(1,20:29,:,:),2);

%btemp=grib.geovariable('Brightness_temperature_atmosphere_top');
%btemp=grib.geovariable('Total_Column-Integrated_Condensate_entire_atmosphere');
%btemp0(:,:)=btemp(1,:,:);

cref1 = grib.geovariable('Composite_reflectivity_entire_atmosphere');
cref0(:,:)=cref1(1,:,:);
%b_cref=cref0(ny-nd:ny+nd,nx-nd:nx+nd);
b_cref=refl0(ny-nd:ny+nd,nx-nd:nx+nd);


directory='/mnt/lfs4/HFIP/hfip-fiu/zhup/output/HAFS_B/track_zhup/';

file_name='natl00l.2018100718.trak.hafs.atcfunix.all';
url_1(1:length(directory))=directory;
url_1(1+length(directory):length(directory)+length(file_name))=file_name;


fid=fopen(url_1,'r');
tmp=fgets(fid);

nm=1;

while length(tmp) > 100
data(1,nm)=str2double(tmp(5:6));
data(2,nm)=str2double(tmp(31:33));
data(3,nm)=str2double(tmp(36:38));
data(4,nm)=str2double(tmp(43:45));
data(5,nm)=str2double(tmp(49:51));
data(6,nm)=str2double(tmp(54:57));
nm=nm+1;
tmp=fgets(fid);
end % while

nm=1;
var=data(1,:);
while isempty(var) == 0
sno(nm)=var(1);
var=var(find(var~=sno(nm)));
nm=nm+1;
end

n1=1;
for n=1:length(data)
if data(1,n) == sno(2)
track(:,n1)=data(:,n);
n1=n1+1;
end
end

n1=1;
track_1(:,n1)=track(:,1);
ch_time=track_1(2,n1);
for n=2:length(track)
if track(2,n) ~= ch_time
n1=n1+1;
track_1(:,n1)=track(:,n);
ch_time=track_1(:,n1);
end
end

nfile=57;

for n=1:nfile
file_o1(n,:)='/mnt/lfs4/HFIP/hfip-fiu/zhup/output/HAFS_B/zhup/zhup_2018100718/';
end
fid=fopen('infile');
file2=fscanf(fid,'%s',[50 inf]);
fclose(fid);
file2=file2';
url(:,1:length(file_o1))=file_o1;
url(:,1+length(file_o1):length(file_o1)+50)=file2(1:nfile,:);

c_lat=track_1(3,nn)/10;
c_lon=-track_1(4,nn)/10;

grib=ncgeodataset(url(nn,:));
variables = grib.variables;

lon1 = grib.geovariable('lon');
lat1 = grib.geovariable('lat');
[LON0,LAT0]=meshgrid(lon1,lat1);

lon0(:)=lon1(:,:)';
lat0(:)=lat1(:,:)';

[aaa,nx]=min(abs(lon0-c_lon));
[aaa,ny]=min(abs(lat0-c_lat));

s_LON=LON0(ny-nd:ny+nd,nx-nd:nx+nd);
s_LAT=LAT0(ny-nd:ny+nd,nx-nd:nx+nd);

t1 = grib.geovariable('Temperature_isobaric');
t0(:,:,:)=t1(1,:,:,:);
s_t=t0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qv1 = grib.geovariable('Specific_humidity_isobaric');
qv0(:,:,:)=qv1(1,:,:,:);
s_qv=qv0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qc1 = grib.geovariable('Cloud_mixing_ratio_isobaric');
qc0(:,:,:)=qc1(1,:,:,:);
s_qc=qc0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qi1 = grib.geovariable('Ice_water_mixing_ratio_isobaric');
qi0(:,:,:)=qi1(1,:,:,:);
s_qi=qi0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qr1 = grib.geovariable('Rain_mixing_ratio_isobaric');
qr0(:,:,:)=qr1(1,:,:,:);
s_qr=qr0(:,ny-nd:ny+nd,nx-nd:nx+nd);

qs1 = grib.geovariable('Snow_mixing_ratio_isobaric');
qs0(:,:,:)=qs1(1,:,:,:);
s_qs=qs0(:,ny-nd:ny+nd,nx-nd:nx+nd);

refl1=grib.geovariable('Reflectivity_isobaric');
refl0(:,:)=mean(refl1(1,20:29,:,:),2);

%btemp=grib.geovariable('Brightness_temperature_atmosphere_top');
%btemp=grib.geovariable('Total_Column-Integrated_Condensate_entire_atmosphere');
%btemp0(:,:)=btemp(1,:,:);

cref1 = grib.geovariable('Composite_reflectivity_entire_atmosphere');
cref0(:,:)=cref1(1,:,:);
%s_cref=cref0(ny-nd:ny+nd,nx-nd:nx+nd);
s_cref=refl0(ny-nd:ny+nd,nx-nd:nx+nd);

b_qct=b_qc+b_qi+b_qr+b_qs;
s_qct=s_qc+s_qi+s_qr+s_qs;

b_qt=b_qct+b_qv;
s_qt=s_qct+s_qv;

load colorscale;

%colormap(colorscale);
colormap('jet');
colormap(1,:)=[1 1 1];

x1=max(min(b_LON(:)),min(s_LON(:)));
x2=min(max(b_LON(:)),max(s_LON(:)));
y1=max(min(b_LAT(:)),min(s_LAT(:)));
y2=min(max(b_LAT(:)),max(s_LAT(:)));


s_qct_500=squeeze(mean(s_qct(21:21,:,:),1))*1e3;
p_ti=squeeze(s_t(:,141,146));
p_tm=squeeze(s_t(:,140,146));
p_to=squeeze(s_t(:,139,146));
p_qvi=squeeze(s_qv(:,141,146));
p_qvm=squeeze(s_qv(:,140,146));
p_qvo=squeeze(s_qv(:,139,146));
p_thi=p_ti.*(1e5./p0).^(287/1004);
p_thei=p_thi.*(1+2500000/1004./p_ti.*p_qvi);
p_thm=p_tm.*(1e5./p0).^(287/1004);
p_them=p_thm.*(1+2500000/1004./p_tm.*p_qvm);
p_tho=p_to.*(1e5./p0).^(287/1004);
p_theo=p_tho.*(1+2500000/1004./p_to.*p_qvo);

save data/en_inst_michael_dat.mat p0 s_LON s_LAT s_qct_500 p_ti p_tm p_to p_qvi p_qvm p_qvo p_thi p_thei p_thm p_them p_tho p_theo;



k1=subplot(2,2,1);
contourf(s_LON,s_LAT,s_qct_500,[0:0.25:3.5],'edgecolor','none');
caxis([0 3.5]);
%contourf(b_LON,b_LAT,b_cref,[20:2.5:55],'edgecolor','none');
%caxis([20 55]);
colorbar;
hold on;
plot(coastline(:,1),coastline(:,2),'color',[0.7 0.7 0.7]);
set(gca,'xtick',[-90:1:-82],'ytick',[24:1:34]);
axis([x1,x2,y1,y2]);
h=title('(a) HAFB');
set(h,'fontweight','normal');
xlabel('Longitude');
ylabel('Latitude');
%h=text(-88.5,33.2,'Composite reflectivity (dBZ) of the entire atmosphere at 12 UTC Oct. 10, 2018');
%set(h,'fontsize',9); 
hold off;

k2=subplot(2,2,2);
contourf(s_LON,s_LAT,squeeze(mean(s_qct(21:21,:,:),1))*1e3,[0:0.25:3.5],'edgecolor','none');
caxis([0 3.5]);
%contourf(b_LON,b_LAT,b_cref,[20:2.5:55],'edgecolor','none');
%caxis([20 55]);
hold on;
plot(s_LON(142,146),s_LAT(142,146),'k*');
plot(s_LON(141,146),s_LAT(141,146),'g*');
plot(s_LON(140,146),s_LAT(140,146),'r*');
plot(coastline(:,1),coastline(:,2),'color',[0.7 0.7 0.7]);
set(gca,'xtick',[-90:1:-82],'ytick',[24:1:34]);
xlabel('Longitude');
ylabel('Latitude');
axis([-86.8,-85,27.7,29.1]);
hold off;

k4=subplot(2,3,4);
h=plot(p_thei,p0/100,'k');
hold on;
h=plot(p_them,p0/100,'g');
h=plot(p_theo,p0/100,'r');
set(gca,'ydir','reverse');
axis([345,360,400,1000]);
hold off;

k5=subplot(2,3,5);
h=plot(p_qvi*1e3,p0/100,'k');
hold on;
h=plot(p_qvm*1e3,p0/100,'g');
h=plot(p_qvo*1e3,p0/100,'r');
set(gca,'ydir','reverse');
axis([5,25,400,1000]);
hold off;

k6=subplot(2,3,6);
den1=(p_them-p_thei)-(p_qvm-p_qvi).*(p_thi+p_thm);
den2=(p_theo-p_them)-(p_qvo-p_qvm).*(p_tho+p_thm);
den3=(p_theo-p_thei)-(p_qvo-p_qvi).*(p_tho+p_thi);

h=plot(den1,p0/100,'k');
hold on;
h=plot(den2,p0/100,'g');
h=plot(den3,p0/100,'r');
plot([0 0],[400 1000],'--');
set(gca,'ydir','reverse');
axis([-10,1,400,1000]);
hold off;


---can

k2=subplot(2,2,2);
%contourf(s_LON,s_LAT,s_cref,[0:1:20],'edgecolor','none');
%caxis([0 20]);
contourf(s_LON,s_LAT,s_cref,[20:2.5:55],'edgecolor','none');
caxis([20 55]);
c2=colorbar('hori');
hold on;
plot(coastline(:,1),coastline(:,2),'color',[0.8 0.8 0.8]);
axis([-89.09,-82.49,25.99,32.41]);
axis([x1,x2,y1,y2]);
h=title('(b) HAFN');
set(h,'fontweight','normal');
set(gca,'xtick',[-90:1:-82],'ytick',[24:1:34]);
xlabel('Longitude');
hold off;

set(k1,'position',[0.065 0.30 0.43 0.57],'fontsize',9);
set(k2,'position',[0.550 0.30 0.43 0.57],'fontsize',9);
%set(c2,'position',[0.225 0.18 0.6 0.03],'xtick',[-20:10:55],'fontsize',9);
set(c2,'position',[0.225 0.18 0.6 0.03],'fontsize',9);

---can

print -dpng Fig9_ref_hori_com

