clear

casen=input('case No.: ');



ci=[-50:5:50];
ci2=[-20:2:20];
cc=[-50 50];
cc2=[-20 20];

lat=[-35.:.5:35.5]+0.25;
lon=[115.:.5:295.5]+0.25;

lev =  [5.0000   15.0000   25.0000   35.0000   45.0000   55.0000   65.0000   75.0000   85.0000   95.0000 ...
  105.0000  115.9789  129.7985  148.0428  171.8625  201.8625  238.0428  279.7985  325.9789  375.0000 ...
 425.0000  491.1513  604.1757  790.2066 1068.2310 1448.2310 1930.2066 2504.1757 3151.1513 3845.0000]';
 
lev=lev(1:20);

pdeni=[17:.2:28]';

figure

subplot(2,1,1)



filename = [sprintf('/export/t2d2/kelvin/iROAM/daily_output/case%d/vm.mat',casen)];
filename2 = [sprintf('/export/t2d2/kelvin/iROAM/daily_output/case%d/vrhom.mat',casen)];
filename3 = [sprintf('/export/t2d2/kelvin/iROAM/daily_output/case%d/hum.mat',casen)];


eval(['load ' filename])

%v=permute(v,[4,3,2,1]);

v=mean(vm,4);

[psi,latp,levp] = Bgrid_psi_ot(v,lat,lon,lev);

psi=psi/1000000;

contourf(latp,levp,psi',ci)
hold on
%contour(latp,levp,psi',[0 0],'k','LineWidth',2)
set(gca,'YDir','reverse')
axis([-30 30 0 250])
label_lat
caxis(cc)
nc=size(ci,2)-1;
%colormap(jet(nc))
cmap=b2rcolormap(ci);
colormap(cmap);
colorbar
xlabel('Latitude','FontSize',14)
ylabel('Depth (m)','FontSize',14)
title(['\Psi_E(z)   (Sv): Case ' int2str(casen)],'FontSize',16)

subplot(2,1,2)

eval(['load ' filename2])
eval(['load ' filename3])

vv=mean(vrhom,4);


% vv define on lat(1:end-1), lon(1:end-1)

lat=lat(1:end-1);

pden=av(pdeni,1);

[LATA,LONA,DZA]=ndgrid(lat,lon,pden);

deg2m=pi/180*6370000;

DX=diff(LONA.*cos(LATA*pi/180),1,2)*deg2m;

Vs=vv.*DX;

Vm=squeeze(nansum(Vs,2));

[na nb]=size(Vm);
psi=zeros(na,nb+1);
psi(:,2:end)=cumsum(Vm,2);

%psi=fliplr(psir);

psi=psi/1000000;

hm=squeeze(mean(hum,4));
hmm=squeeze(mean(hm,2));

% project onto z

z=[0:10:500];

[n0 nz]=size(z);
[na nb nc]=size(Vs);
psii=zeros(na,nb,nc+1);
psii(:,:,2:end)=cumsum(Vs,3);


psiiz=NaN*zeros(na,nb,nz);

for i=1:na
for k=1:nb

  hma=squeeze(hm(i,k,:))';
  jcount=0;
  for j=1:nc+1
   if hma(j) > -10
     jcount=jcount+1;
   else
     if jcount == 0
     hma(j)=-100+0.1*(j-1);
     else
     hma(j)=500+0.1*(j-1);
     end
   end
  end

if jcount > 1
psa=squeeze(psii(i,k,:));
psiiz(i,k,:)=interp1(hma,psa',z);
end

end
end

psiz=squeeze(nansum(psiiz,2))/1000000;


contourf(lat,z,psiz',ci)
hold on
%contour(lat,z,psiz',[0 0],'k','LineWidth',2)
set(gca,'YDir','reverse')
axis([-30 30 0 250])
label_lat
caxis(cc)
nc=size(ci,2)-1;
%colormap(jet(nc))
cmap=b2rcolormap(ci);
colormap(cmap);
colorbar
xlabel('Latitude','FontSize',14)
ylabel('Depth (m)','FontSize',14)
title(['\Psi_L(z)   (Sv): Case ' int2str(casen)],'FontSize',16)

figure
contourf(lat,pdeni,psi',ci2)
hold on
%contour(lat,pdeni,psi',[0 0],'k','LineWidth',2)
set(gca,'YDir','reverse')
axis([-30 30 21.5 26])
label_lat
caxis(cc2)
nc=size(ci2,2)-1;
%colormap(jet(nc))
cmap=b2rcolormap(ci2);
colormap(cmap);
colorbar
xlabel('Latitude','FontSize',14)
ylabel('\sigma','FontSize',14)
title(['\Psi_L(\sigma)   (Sv): Case ' int2str(casen)],'FontSize',16)

