clear all;
ncid=netcdf.open('/data2/giant/OCEAN_OUTPUT/CLOGIN/yht_GTIAFa/TOPO.nc','nowrite');

ma=netcdf.inqVarID(ncid,'mask');kb=netcdf.inqVarID(ncid,'kb');
lo=netcdf.inqVarID(ncid,'lon');la=netcdf.inqVarID(ncid,'lat');
m=netcdf.getVar(ncid,la);n=netcdf.getVar(ncid,lo);kb=netcdf.getVar(ncid,kb);
lat=m(:,1);lon=n(:,1);
lev=netcdf.inqVarID(ncid,'lev_c');
zlev=netcdf.getVar(ncid,lev)./100;
lev2=netcdf.inqVarID(ncid,'lev_f');
zlev2=netcdf.getVar(ncid,lev2)./100;
dz=diff(zlev2);

mask=netcdf.getVar(ncid,ma);mask=repmat(double(mask(:,:,1)),1,1,12);
mask1=mask(:,:,1);
mask_R1=mask1-mask1;

mask_g=mask1;
mask_g(find(mask_g==0))=nan;

lat_s=-56;lat_n=65;lon_w=25;lon_e=290;
[cc,la0]=min(abs(lat-lat_s));[cc,lo0]=min(abs(lon-lon_w));[cc,la_e]=min(abs(lat-lat_n));[cc,lo_e]=min(abs(lon-lon_e));
mask_R1(lo0:lo_e,la0:la_e,:,:)=mask1(lo0:lo_e,la0:la_e,:,:);
lat_s=7;lat_n=51;lon_w=282;lon_e=290;
[cc,la0]=min(abs(lat-lat_s));[cc,lo0]=min(abs(lon-lon_w));[cc,la_e]=min(abs(lat-lat_n));[cc,lo_e]=min(abs(lon-lon_e));
mask_R1(lo0:lo_e,la0:la_e,:,:)=nan;
lat_s=16;lat_n=45;lon_w=260;lon_e=290;
[cc,la0]=min(abs(lat-lat_s));[cc,lo0]=min(abs(lon-lon_w));[cc,la_e]=min(abs(lat-lat_n));[cc,lo_e]=min(abs(lon-lon_e));
mask_R1(lo0:lo_e,la0:la_e,:,:)=nan;
lat_s=55;lat_n=70;lon_w=260;lon_e=300;
[cc,la0]=min(abs(lat-lat_s));[cc,lo0]=min(abs(lon-lon_w));[cc,la_e]=min(abs(lat-lat_n));[cc,lo_e]=min(abs(lon-lon_e));
mask_R1(lo0:lo_e,la0:la_e,:,:)=nan;
lat_s=10;lat_n=16;lon_w=277;lon_e=282;
[cc,la0]=min(abs(lat-lat_s));[cc,lo0]=min(abs(lon-lon_w));[cc,la_e]=min(abs(lat-lat_n));[cc,lo_e]=min(abs(lon-lon_e));
mask_R1(lo0:lo_e,la0:la_e,:,:)=nan;
lat_s=26;lat_n=42;lon_w=25;lon_e=40;
[cc,la0]=min(abs(lat-lat_s));[cc,lo0]=min(abs(lon-lon_w));[cc,la_e]=min(abs(lat-lat_n));[cc,lo_e]=min(abs(lon-lon_e));
mask_R1(lo0:lo_e,la0:la_e,:,:)=nan;


mask_R1(find(mask_R1==0))=nan;
mask_R2=mask1-mask1;

lat_s=-39;[cc,la0]=min(abs(lat-lat_s));mask_R1(:,1:la0,:,:)=nan;
lat_n=86;lat_s=-56;
[cc,la0]=min(abs(lat-lat_s));[cc,la_e]=min(abs(lat-lat_n));
mask_R2(:,la0:la_e,:,:)=mask1(:,la0:la_e,:,:);
mask_R2(find(mask_R1==1))=nan;
mask_R2(find(mask_R2==0))=nan;
lat_s=-35;[cc,la0]=min(abs(lat-lat_s));mask_R2(:,1:la0,:,:)=nan;
maskR1=repmat(mask_R1,1,1,45);
maskR2=repmat(mask_R2,1,1,45);

fid=figure; FigHandle = fid;set(FigHandle, 'Position', [5, 31, 1900, 400]);hold on
subplot(1,3,1);pcolor(lon,lat,mask1');shading flat;title('Indo-Pacific ocean mask')
subplot(1,3,2);pcolor(lon,lat,mask_R1');shading flat;title('Indo-Pacific ocean mask')
subplot(1,3,3);pcolor(lon,lat,mask_R2');shading flat;title('Atlantic-Arctic ocean mask')
save basin_mask mask_R1 mask_R2 maskR1 maskR2 mask_g mask1


ncid=netcdf.open('/data2/giant/OCEAN_OUTPUT/CLOGIN/yht_GTJRAa/GTJRAa.cice.h.0001-11.nc','nowrite');
d=netcdf.inqVarID(ncid,'dxt');dx=netcdf.getVar(ncid,d)./1000;
d=netcdf.inqVarID(ncid,'dyt');dy=netcdf.getVar(ncid,d)./1000;
d=netcdf.inqVarID(ncid,'tarea');areaGd=netcdf.getVar(ncid,d);
d=netcdf.inqVarID(ncid,'lont_bounds');lon_bnd=netcdf.getVar(ncid,d);
d=netcdf.inqVarID(ncid,'latt_bounds');lat_bnd=netcdf.getVar(ncid,d);
kbmask=double(kb).*mask1;

% dx=zeros(size(lon,1),size(lat,1));dy=dx;
% for i=2:size(lon,1)-1
%     for j=2:size(lat,1)-1
% dx(i,j)=pos2dist(lat(j),(lon(i-1)+lon(i))./2,lat(j),(lon(i)+lon(i+1))./2,2);
% dy(i,j)=pos2dist((lat(j)+lat(j-1))./2,lon(i),(lat(j)+lat(j+1))./2,lon(i),2);
%     end
% end
% dx(1,1:size(lat,1))=dx(2,1:size(lat,1));
% dx(size(lon,1),1:size(lat,1))=dx(size(lon,1)-1,1:size(lat,1));
% dx(1:size(lon,1),1)=dx(1:size(lon,1),2);
% dx(1:size(lon,1),size(lat,1))=dx(1:size(lon,1),size(lat,1)-1);
% 
% dy(1,1:size(lat,1))=dy(2,1:size(lat,1));
% dy(size(lon,1),1:size(lat,1))=dy(size(lon,1)-1,1:size(lat,1));
% dy(1:size(lon,1),1)=dy(1:size(lon,1),2);
% dy(1:size(lon,1),size(lat,1))=dy(1:size(lon,1),size(lat,1)-1);

% subplot(1,2,1);
% pcolor(dx');shading flat;colorbar;title('dx')
% subplot(1,2,2)
% pcolor(dy');shading flat;colorbar;title('dy')
% print('-djpeg','-r300','dxdy.jpg');
%areaGd=dx.*dy.*1e6;
% start and end index  for MOC_basin1 
lat_s=-56;lat_n=65;lon_w=25;lon_e=250;
[cc,la0b1]=min(abs(lat-lat_s));[cc,lo0b1]=min(abs(lon-lon_w));
[cc,la_eb1]=min(abs(lat-lat_n));[cc,lo_eb1]=min(abs(lon-lon_e));

% start and end index  for MOC_basin2 
lat_s=-35;lat_n=80;lon_w=290;lon_e=360;
[cc,la0b2]=min(abs(lat-lat_s));[cc,lo0b2]=min(abs(lon-lon_w));
[cc,la_eb2]=min(abs(lat-lat_n));[cc,lo_eb2]=min(abs(lon-lon_e));

lat_s=-35;lat_n=80;lon_w=0;lon_e=25;
[cc,lo0b2a]=min(abs(lon-lon_w));[cc,lo_eb2a]=min(abs(lon-lon_e));


%earthellipsoid = referenceSphere('earth','km');
%for i=2:size(lon,1)-1
%    for j=2:size(lat,1)-1
%areaGd(i,j)=areaquad((lat(j)+lat(j-1))./2,(lon(i)+lon(i-1))./2,...
%    (lat(j)+lat(j+1))./2,(lon(i)+lon(i+1))./2,earthellipsoid).*1e6;
%    end
%end


dz1(1,1,:)=dz;
crosArea=repmat(dx,1,1,45).*repmat(dz1,size(dx,1),size(dx,2),1).*1000;
crosArea2=repmat(dy,1,1,45).*repmat(dz1,size(dx,1),size(dx,2),1).*1000;

for z=1:45
   A=find(kbmask<z);
   B=areaGd;
   B(A)=nan;
areaGd_z(:,:,z)=B;
end
sumarea_z=squeeze(nansum(nansum(areaGd_z,1),2));
dzz=permute(repmat(dz,1,size(areaGd_z,1),size(areaGd_z,2)),[2,3,1]);
vol_z=areaGd_z.*dzz;


% itt=1;
% lev=[1 20 25 30 41 42 43 44 45];
% fid=figure; FigHandle = fid;set(FigHandle, 'Position', [5, 31, 1900, 960]);hold on
% for z=lev(1:9)
% subplot(3,3,itt);
% pcolor(areaGd_z(:,:,z)');shading flat;title(strcat('depth=',num2str(zlev(z))));colorbar
% itt=itt+1;
% end
% print('-djpeg','-r300','area_z.jpg');

lat_s=-35;[cc,la_s35]=min(abs(lat-lat_s));lat_s=-39;[cc,la_s39]=min(abs(lat-lat_s));
dx0=dx;dy0=dy;

dep=zlev(kb);dep1=dep.*mask1;
[lon1 lat1]=meshgrid(lon, lat);
lev_bnd(:,1)=zlev2(1:45);lev_bnd(:,2)=zlev2(2:46);
save areaGd areaGd sumarea_z kb dx0 dy0 crosArea crosArea2 areaGd_z vol_z la_s35 la_s39


addpath('/home/yckuo/MATLAB/omip/giant/create_omip_nc/')
mkdir Ofx/areacello; cd Ofx/areacello ;create_omip_nc_area;cd ../../
mkdir Ofx/volcello; cd Ofx/volcello ;create_omip_nc_vol ;cd ../../
mkdir Ofx/deptho; cd Ofx/deptho ;create_omip_nc_deptho;cd ../../
mkdir Ofx/sftof; cd Ofx/sftof ;create_omip_nc_sftof;cd ../../

