clear all;
%web.corral.tacc.utexas.edu/OceanProjects/ASTE/Release1/nctiles_daily/
dirroot='/scratch/projects/ecco/atnguyen/aste_270x450x180/ASTE_daily_from_pfe/';
dirgrid=[dirroot 'nctiles_grid/'];
dirdata=[dirroot 'output_tiles_nc/UVELMASS/'];
dirout=['/work/03901/atnguyen/BE_for_Rebecca/'];if(~exist(dirout));mkdir(dirout);end;

%from this figure web.corral.tacc.utexas.edu/OceanProjects/ASTE/Release1/ASTE_tiles.png
%we want tile 19
tilestr=sprintf('%4.4i',19);

%uvelmass = uvel*hfacw = mass-weighted velocity [m/s]
%web.corral.tacc.utexas.edu/OceanProjects/ASTE/Release1/nctiles_daily/time_list_daily.txt
tt=load([dirroot 'time_list_daily.txt']);
ttp=datenum(tt(:,1:3));ttp=ttp-ttp(1);ttp=ttp./365.25+2002;	%yr since 2002

%web.corral.tacc.utexas.edu/OceanProjects/ASTE/Release1/nctiles_daily/UVELMASS.0019.nc.gz
%download file .gz and unzip using gunzip:
tic;ncdisp([dirdata 'UVELMASS.' tilestr '.nc']);
tic;u=ncread([dirdata 'UVELMASS.' tilestr '.nc'],'UVELMASS');toc;%90x90x50x5844days

%web.corral.tacc.utexas.edu/OceanProjects/ASTE/Release1/nctiles_grid/
%load in grid, 
%hfacw is just for cosmetic location of Bering strait 
hfw=ncread([dirgrid 'GRID.' tilestr '.nc'],'hFacW');

%latitude at the left bottom corner of cell
yg=ncread([dirgrid 'GRID.0019.nc'],'YG');
%centering yg along the left edge to correspond to u point on c-grid
ygc=nan(size(yg));ygc(:,1:end-1)=(yg(:,1:end-1)+yg(:,2:end))./2;

%longitude at the left bottom corner of cell
xg=ncread([dirgrid 'GRID.0019.nc'],'XG');
%centering xg along the left edge to correspond to u point on c-grid
xgc=nan(size(xg));xgc(:,1:end-1)=(xg(:,1:end-1)+xg(:,2:end))./2;

%depth cell center:
rc=ncread([dirgrid 'GRID.0019.nc'],'RC'); %[m], negative

%depth dz:
drf=ncread([dirgrid 'GRID.0019.nc'],'DRF'); %[m]

%dy along left edge , for use with velocity u:
dyg=ncread([dirgrid 'GRID.0019.nc'],'DYG'); %[m]

%visually locating Bering STrait by plotting hfw:
figure(1);clf;colormap(jet(21));
    subplot(321);imagesc(1:90,1:90,hfw(:,:,1)');colorbar;
    set(gca,'layer','top','ydir','normal');grid;
    hold on;plot(ix.*ones(size(iy)),iy,'k.-');hold off;

%now pick the lcoation, and plot cross section:
ix=42;iy=52:60;iz=1:6;%(after picking, plot that line above)
figure(1);subplot(322);
   imagesc(squeeze(xgc(ix,iy)),rc(iz),squeeze(hfw(ix,iy,iz))');colorbar;shading flat;
   set(gca,'layer','top','ydir','normal');grid;xlabel('longitude');ylabel('depth [m]');

%now compute transport as m/s*dyg*drf , and note from subplot 221 that NEGATIVE u is northward:
tmpu=squeeze(u(ix,iy,iz,:));	%9 x 6 x 5844

%cross-sectional trsp, note the minus for northward, size 9 x 6 x 5844
trspu=-tmpu.*repmat(dyg(ix,iy)',[1 length(iz) length(tt)]).*repmat(drf(iz)',[length(ix) 1 length(tt)]);

trspu_xsec_mean=mean(trspu,3);	   %9x6
trsputot=squeeze(sum(sum(trspu))); %1x5844

%create a structure:
be=[];
be.tilestr=tilestr;
be.ix=ix;
be.jy=iy;
be.iz=iz;
be.dyg=squeeze(dyg(ix,iy));	%distance in model y-dir in [m] (not along lat or lon, just numerical dir)
be.drf=squeeze(drf(iz));	%distance in model z-dir in [m] (depth dir)
be.rc=squeeze(rc(iz));	        %depth at cell center [m] (negative)
be.xgc=squeeze(xgc(ix,iy));	%LONGITUDE across Bering Strait (deg E)
be.ygc=squeeze(ygc(ix,iy));	%LATITUDE across Bering Strait (deg N)
be.trsp = trspu;		%volume flux in cross-sectional and with time, m3/s
be.trsptot = trsputot;		%net volume flux (total across whole section), m3/s
be.tt = tt;			%time array
be.ttp=ttp;			%time using datenum, minus first datenum(2002,1,1)
%plot time-series:
figure(1);
   subplot(323);imagesc(xgc(ix,iy),rc(iz),trspu_xsec_mean');colorbar;shading flat
   set(gca,'layer','top','ydir','normal');grid;xlabel('longitude');ylabel('depth [m]');

   subplot(324);plot(ttp,trsputot./1e6,'-');grid;xlabel('yr');ylabel('northward trsp Sv');
   title(['time mean transp: ' num2str(mean(trsputot)/1e6) ' Sv']);

%in the case you'd like to compute salt transp and heat transport (ref to 0degC)
%the field to download is SALT/SALT.0019.nc , THETA/THETA.0019.nc
%here i show an example to compute heatflux:
%read in theta, pot temp, at cell center 
clear T
tic;T=ncread([dirroot 'output_tiles_nc/THETA/THETA.' tilestr '.nc'],'THETA');toc;%90x90x50x5844days

%note u is at left edge, so we can't just yet mult u*T, but rather we just move the T onto the u point:
T_at_u=nan(size(T));
T_at_u(2:90,:,:,:)=(T(1:89,:,:,:)+T(2:90,:,:,:))./2;
rhoconst=1029;	%kg/m3
cp=3994;	%specific heat J/kg/degC
T_ref=-1.9;	%i just show for T_ref=-1.9 here;

%
tmpu=squeeze(u(ix,iy,iz,:));	%9 x 6 x 5844
tmpT=squeeze(T_at_u(ix,iy,iz,:));%9 x 6 x 5844
heatflux_adv=-tmpu.*(tmpT-T_ref)...	%note minus sign, so that now positive u is northward
             .*repmat(dyg(ix,iy)',[1 length(iz) length(tt)])...
             .*repmat(drf(iz)',[length(ix) 1 length(tt)]) ...
             .*rhoconst.*cp;		%m/s*degC*m*m* kg/m3 * J/kg/degC = J/s, 9 x 6 x 5844

heatflux_adv_tot=squeeze(sum(sum(heatflux_adv)));	%J/s, net flux , 1 x 5844 

be.heatflux_adv = heatflux_adv;	        %eulerian advective part, based on daily velocity and T
be.heatflux_adv_tot = heatflux_adv_tot;	%eulerian advective part, based on daily velocity and T

figure(1);
   subplot(325);imagesc(xgc(ix,iy),rc(iz),mean(heatflux_adv,3)');colorbar;shading flat
   set(gca,'layer','top','ydir','normal');grid;xlabel('longitude');ylabel('depth [m]');

   subplot(326);plot(ttp,heatflux_adv_tot./1e15,'-');grid;
   xlabel('yr');ylabel('northward heat trsp PW ref to -1.9');

figure(1);set(gcf,'paperunits','inches','paperposition',[0 0 16 12]);
fpr=[dirout 'BeringStrait_volflux_ASTER1_daily.png'];print(fpr,'-dpng');

fsave=[dirout 'BeringStrait_volflux_ASTER1_daily.mat'];
save(fsave,'be');
