fname='NVELMASS.nc';
readncfile

% You need the llgrid.mat in the directory if you are using lonlat data
% You need the intgrid.mat in the directory if you are using interp data
if exist('llgrid.mat')==2 
  load llgrid.mat  
  DRF=llgrid.DRF;
  DXV=llgrid.DXV;
  RC=llgrid.RC;
  RF=llgrid.RF;
  basins=llgrid.basins;
elseif exist('intgrid.mat')==2 
  load intgrid.mat
  DRF=intgrid.DRF;
  RC=intgrid.RC;
  RF=intgrid.RF;
  DXV=0.5.*111000.*cosd(lat);  %This is to reconstruct the size of 0.5 interpolated grid
  basins=intgrid.basins;
end

% Next we will calculate the meridional transport streamfunction
% It is v(y,z,t)=\int_West^East \int_z^0 v dz dx

ntrans=zeros(size(squeeze(NVELMASS(1,:,:))));
gmoc=zeros(size(squeeze(NVELMASS(1,:,:))));

% zonal transport integral
for zz=1:50
 ntrans(:,zz)=squeeze(nansum(NVELMASS(:,:,zz).*DXV,1)); 
end

figure(1)
contourf(lat(1,:),RC,ntrans');shading('flat');colorbar
xlabel('latitude (degN)')
ylabel('depth (m)')
title('Zonally-Integrated Northward Transport (m^2/s)')

% vertical integral
for zz=2:50
 gmoc(:,zz)=gmoc(:,zz-1)+ntrans(:,zz-1).*DRF(zz-1);
end

%gmoc=gmoc/1e6;  % Convert to Sverdrups

figure(2)
contourf(lat(1,:),RF((end-1):-1:1),gmoc(:,end:-1:1)'./1e6);shading('flat');caxis([-20 20]);colorbar
xlabel('latitude (degN)')
ylabel('depth (m)')
title('Global Meridonal Overturning Streamfunction (Sv)')


figure(3)
contourf(lon,lat,llgrid.basins)
title('Basin Code')
xlabel('longitude (degE)')
xlabel('latitude (degN)')

basmoc=zeros([size(gmoc) 17]);

for bas=1:17
 ntransb=0.*ntrans;
 gmocb=0.*gmoc;
for zz=1:50
 ntransb(:,zz)=squeeze(nansum(NVELMASS(:,:,zz).*DXV.*(basins==bas),1)); 
end

% vertical integral
for zz=2:50
 gmocb(:,zz)=gmocb(:,zz-1)+ntransb(:,zz-1).*DRF(zz-1);
end
   
 basmoc(:,:,bas)=gmocb;
end

figure(4)
num=0;
for bas=[1:3 5:17]
num=num+1;
subplot(4,4,num)
contourf(lat(1,:),RF((end-1):-1:1),squeeze(basmoc(:,end:-1:1,bas))'./1e6);shading('flat');caxis([-20 20]);colorbar
xlabel('latitude (degN)')
ylabel('depth (m)')
title(['MOC (Sv) in Basin ',num2str(bas)])
end