clear
path(path,'/fluid/disk2/yhtseng/research/MITGCM/matlab')
path(path,'/fluid/disk2/yhtseng/research/plotting')
%resol:=80,160,240
resol=240;
if resol==80
RESOL='80_62_9/';
TIME=1440;
elseif resol==160
RESOL='160_122_9/';
TIME=14400;
elseif resol==240
RESOL='240_182_9/';
TIME=86400;
end
%1.full_cell_data
%2.part_cell_data
%3.IBM_data
item=2
if item==1
DIRS='full_cell_data/';
elseif item==2
DIRS='part_cell_data/';
else
DIRS='data/';%data/,full_cell_data/,part_cell_data/
end
depth=rdmds([RESOL 'part_cell_data/' 'Depth']);
xg=rdmds([RESOL DIRS 'XG']);yg=rdmds([RESOL DIRS 'YG']);
xc=rdmds([RESOL DIRS 'XC']);yc=rdmds([RESOL DIRS 'YC']);
dxg=rdmds([RESOL DIRS 'DXG']);dyg=rdmds([RESOL DIRS 'DYG']);
dxc=rdmds([RESOL DIRS 'DXC']);dyc=rdmds([RESOL DIRS 'DYC']);
%hc=rdmds([RESOL DIRS 'hFacC']);hs=rdmds([RESOL DIRS 'hFacS']);hw=rdmds([RESOL DIRS 'hFacW']);
hc=rdmds([RESOL 'part_cell_data/' 'hFacC']);hs=rdmds([RESOL 'part_cell_data/' 'hFacS']);hw=rdmds([RESOL 'part_cell_data/' 'hFacW']);

[u,itu]=rdmds([RESOL DIRS 'U'],TIME);
[v,itv]=rdmds([RESOL DIRS 'V'],TIME);
[w,itw]=rdmds([RESOL DIRS 'W'],TIME);
[t,itt]=rdmds([RESOL DIRS 'T'],TIME);
[s,its]=rdmds([RESOL DIRS 'S'],TIME);
if itu~=itv
 error('Different time levels in U and V')
end
[nx ny nz]=size(u);
dx=dxg(1,1);dy=dx;
H=4500;f=1e-4;dz=H/(nz-1);z=-dz/2:-dz:-H+dz/2;
L=25e3;dh=0.9*H;x_center=xg(nx/2+1,1);y_center=yg(1,ny/2+1);
if item==3
U_temp=u.*hw*dz;
V_temp=v.*hs*dz;
else
U_temp=u.*hw*dz;
V_temp=v.*hs*dz;
end
for k=1:nz-1
VOR(:,:,k)=(V_temp(:,:,k)'-V_temp([end 1:end-1],:,k)')/dx-(U_temp(:,:,k)'-U_temp(:,[end 1:end-1],k)')/dy;
end
U=sum(U_temp(:,:,1:nz-1),3);
V=sum(V_temp(:,:,1:nz-1),3);
figure(item)
display=[-950 -750 -650 -550 -450 -350 -250 -150 -50 150 250 350 450 550 650 750 850];
%display=[-0.85:.1:1.05]*1000;
clf
Z_temp=(V-V([end 1:end-1],:))/dx-(U-U(:,[end 1:end-1]))/dy;
Z_TEMP=Z_temp(1:round(nx/4),:);
Z_TEMP(Z_temp(1:round(nx/4),:)<100*f)=0;
Z_temp(1:round(nx/4),:)=Z_TEMP;
Z_temp(:,1:2)=0;Z_temp(:,ny)=0;
contourf(xg(:,1:ny-2)',yg(:,1:ny-2)',Z_temp(:,2:ny-1)'/f,display);
%contourf(xg(:,1:ny-2)',yg(:,1:ny-2)',VOR(2:ny-1,:,1)/f,display1);

colorbar('h')
figure(item+3)
clf
ind=ny/2;
contourfv(xc(:,ind),z,depth(:,ind),t(:,ind,1:nz-1));
[X Y Z]=meshgrid(xc(:,1),yc(1,:),z);
figure(item+6)
ind=7;quiver(X(:,:,ind)',Y(:,:,ind)',U_temp(:,:,ind),V_temp(:,:,ind))
%Output the max of the vortice
vor_is=round(nx/4*3);vor_ie=round(nx-1);
vor_js=round(ny/6);vor_je=round(ny/2);
max(max(Z_temp(vor_is:vor_ie,vor_js:vor_je)))/f
if nz>10
figure(10)
clf
sr=35*f;
j_s=round(ny/8);j_e=ny-j_s+1;
i_s=round(nx/4);i_e=nx;
k_s=1;k_e=nz-1;
isosurface(X(j_s:j_e,i_s:i_e,k_s:k_e),Y(j_s:j_e,i_s:i_e,k_s:k_e),Z(j_s:j_e,i_s:i_e,k_s:k_e),VOR(j_s:j_e,i_s:i_e,k_s:k_e),sr);
%lighting phong
camlight right
view([-20,20]);
axis([0 4e5 0 3e5 -H 0])
hold on
h=-H+dh*exp( -((X(:,:,1)-x_center).^2+(Y(:,:,1)-y_center).^2)/(2*(L^2)));
surf(X(:,:,1),Y(:,:,1),h)
lighting phong
camlight right
camlight left
end
