clear;

load data/midterm_isabel_dat;
load data/xland_isabel;
u1=squeeze(u(1,:,:));
v1=squeeze(v(1,:,:));
wind1=sqrt(u1.^2+v1.^2);

s_pres=squeeze(pres(1,:,:));
s_pres(XLAND<2)=NaN;

[a1, ny]=min(min(s_pres,[],2));
[a1, nx]=min(min(s_pres,[],1));

xlat1=XLAT;
xlong1=XLONG;


[d3 d2 d1]=size(Hydrometer);

lat_center1=xlat1(ny,nx);
lon_center1=xlong1(ny,nx);


for i=1:d1
for j=1:d2
x1=xlong1(j,i)-lon_center1;
y1=xlat1(j,i)-lat_center1;
cour(j,i)=course(x1,y1);
dist=sqrt((i-nx)^2+(j-ny)^2);
if dist > 1.0
 alpha=atan(v1(j,i)/u1(j,i))*180/pi;
 if cour(j,i)<=180
  if alpha > 0
  phi1(j,i)=alpha+cour(j,i)-180;
   if phi1(j,i) < -90
   phi1(j,i)=phi1(j,i)+180;
   end
  else
  phi1(j,i)=cour(j,i)+alpha;
   if phi1(j,i) > 90
   phi1(j,i)=phi1(j,i)-180;
   end
  end
 else
  if alpha >0
  phi1(j,i)=alpha+cour(j,i)-360;
   if phi1(j,i) < -90
   phi1(j,i)=phi1(j,i)+180;
   end
  else
  phi1(j,i)=alpha+cour(j,i)-180;
   if phi1(j,i) > 90
   phi1(j,i)=phi1(j,i)-180;
   end
  end
 end
else
 phi1(j,i)=0;
end
end
end
spd_t=wind1.*cos(phi1*pi/180);
spd_r=wind1.*sin(phi1*pi/180);


%---method 1
lwp(1:d2,1:d1)=0;
for k=1:d3
    lwp(:,:)=lwp(:,:)+squeeze(Hydrometer(k,:,:).*(z(k+1,:,:)-z(k,:,:)));
end

%---method 2
dz=z(2:end,:,:)-z(1:end-1,:,:);
var=Hydrometer.*dz;
lwp1=squeeze(sum(var,1));

figure(1);
subplot(2,2,1);
contourf(lwp,'edgecolor','none');
colorbar;

subplot(2,2,2);
contourf(lwp1,'edgecolor','none');
colorbar;

figure(2);
contourf(lwp,[0:0.5:50],'edgecolor','none');
colorbar;




