clear all
addpath('./ASTE_mlibrary/');			%<-- modify absolute path as needed
addpath('./ASTE_mlibrary/gcmfaces_mod/');	%<-- modify absolute path as needed 

%the most important structure needed is the tile, located at
dirdiags_nc='./';				%<-- modify absolute path as needed
dirgrid_nc='./nctiles_grid/';			%<-- modify absolute path as needed
load(['./ASTE_mlibrary/tiles_params.mat']);
%tiles = 
%  struct with fields:
%      tot: 45
%      wet: 29
%     list: [45×2 double]
%    tsize: [90 90]
%    fsize: [270 1350 1]
%      nfx: [270 0 270 180 450]
%      nfy: [450 0 270 270 270]
%     mask: {[270×450 double]  []  [270×270 double]  [180×270 double]  [450×270 double]}
%      lat: 13.7500
%    field: []
%    index: {1×29 cell}

%Now, we read in the nc tiles for the ENTIRE ASTE domain, output is compact format [270 1350 length(iz) length(it)];
%Note this will be VERY BIG array, watch out for memory.
%It is recommended, if one want to work with full ASTE domain, to read 1 timestep or 2 time-steps at a time instead of the full time-series.
%e.g., set itime=1:2 or itime = a number between 1 and 192;

%=============================================================================
%=============================================================================
%======================FULL ASTE  DOMAIN =====================================
%=============================================================================
%
%example with tracer field:
iz=[1,5];		%0-10m (center at 5m, dz=10m) , 40-50m (center at 45m, dz = 10m)
itime=[167:168];	%entire time series 192 months cover Jan/2002 to Dec/2017, 168 = Dec/2016;
itile=1:29;		%use all 29 tiles to get entire ASTE domain
flag_aste=1;		%1 to put output into full ASTE domain compact size, 0 to leave as individual tiles (to save memory)
[field_full_aste time_step]=nctiles2aste_v2('THETA',[dirdiags_nc '/nctiles_monthly/'],itile,iz,itime,tiles,1);

size(field_full_aste)	%[270 1350 2-depth-levels 2-monthly-means] , takes some time to generate if itime=1:192;

%==========================================
%put the full ASTE compact domain into the 5 llc faces:
field_faces=get_aste_faces(field_full_aste,tiles.nfx,tiles.nfy);	%5 faces

%field_faces =
%  1×5 cell array
%    {270×450×2×2 double}    {0×0 double}    {270×270×2×2 double}    {180×270×2×2 double}    {450×270×2×2 double}

%plot first level, month 167:
for i=[1,3:5];tmp{i}=field_faces{i}(:,:,1,1);end;
fignum=11;crange=[];plot_aste_faces(tmp,fignum,crange);

%==========================================
%now showing handeling with vector fields
[u_full_aste time_step]=nctiles2aste_v2('UVELMASS',[dirdiags_nc '/nctiles_monthly/'],itile,iz,itime(1),tiles,1);	%[270 1350 2]
[v_full_aste time_step]=nctiles2aste_v2('VVELMASS',[dirdiags_nc '/nctiles_monthly/'],itile,iz,itime(1),tiles,1);	%[270 1350 2]

u_faces_Cgrid=get_aste_faces(u_full_aste,tiles.nfx,tiles.nfy);	%5 faces, at C-grid locations
v_faces_Cgrid=get_aste_faces(v_full_aste,tiles.nfx,tiles.nfy);	%5 faces, at C-grid locations

%==========================================
% Some example calculations:
%put into 5 faces, add 1 extra row/column padding to enable centering an/or calculation of divergence and curl
%note that the padding gives 1 extra row/column, in contrast to the simple get_aste_faces function
[ufaces,vfaces]=llc_uv_addpadding(u_full_aste,v_full_aste,tiles.nfx,tiles.nfy);

%ufaces =
%  1×5 cell array
%    {271×450×2 double}    {0×0 double}    {271×270×2 double}    {181×270×2 double}    {451×270×2 double}
%vfaces =
%  1×5 cell array
%    {270×451×2 double}    {0×0 double}    {270×271×2 double}    {180×271×2 double}    {450×271×2 double}

%==========================================
% Get centered vel (which requires/calls padding function) and calc speed 
[uc,vc]=aste_uv_centering(u_full_aste,v_full_aste,tiles.nfx,tiles.nfy);

%calculate speed
sp=[];
for i=[1,3:5];
  sp{i}=sqrt(uc{i}.^2+vc{i}.^2);
end;

%calculate divergence, note need to read in the grid dxg and dyg
grid_nc=nctiles_grid({'DXG','DYG'},dirgrid_nc,itile,tiles,1);
%grid_nc = 
%  struct with fields:
%    DXG: [270×1350 double]
%    DYG: [270×1350 double]

hconv=aste_uv_hconv(u_full_aste,v_full_aste,grid_nc.DXG,grid_nc.DYG,tiles.nfx,tiles.nfy);
%hconv =
%  1×5 cell array
%    {270×450×2 double}    {0×0 double}    {270×270×2 double}    {270×270×2 double}    {450×270×2 double}

%=============================================================================
%=============================================================================
%======================INDIVIDUAL TILE =======================================
%=============================================================================

%Read in 1 tile instead of expanding out to entire ASTE compact size (to save memory)

iz=1:50;		%for this example, choose entire water column
itile=[11,27];		%choose 11 = Fram Strait region, 27=Davis Strait region
flag_aste=0;		%1 to put output into full ASTE domain compact size, 0 to leave as individual tiles (to save memory)
itime=1:192;		%read entire time-series, 192 months
[field_tile time_step]=nctiles2aste_v2('THETA',[dirdiags_nc '/nctiles_monthly/'],itile,iz,itime,tiles,0);
[u_tile time_step]=nctiles2aste_v2('UVELMASS',[dirdiags_nc '/nctiles_monthly/'],itile,iz,itime,tiles,0);
%field_tile =
%  1×2 cell array
%    {90×90×50×192 double}    {90×90×50×192 double}

%read in also grid [XC,YC,hFacC]
grid_nc=nctiles_grid({'XC','YC','hFacC'},dirgrid_nc,itile,tiles,0);

%>> grid_nc =
%  1×2 cell array
%    {1×1 struct}    {1×1 struct}
%
%>> grid_nc{1}
%  struct with fields:
%       XC: [90×90 double]
%       YC: [90×90 double]
%    hFacC: [90×90×50 double]

%look at Fram Strait:
tmp=squeeze(field_tile{1});	%[90 x 90 x 50 x 192]
tmpu=squeeze(u_tile{1});	%[90 x 90 x 50 x 192], for this tile, +u is flowing from Nordic Seas into the Arctic/Barents Sea
sz=size(tmp);

ix=55;iy=89:-1:39;		%eye-balling indices of Fram Strait
pt_i=47;pt_j=43;		%a point within the West Spitsbergen Current

figure(12);clf;colormap(parula(21));
  subplot(221);imagesc(1:sz(1),1:sz(2),tmp(:,:,1,192)');colorbar;grid;set(gca,'ydir','normal');title('surface temperature for Dec/2017');
    hold on;plot([ix ix],[iy(1) iy(end)],'-','color',.9.*[1 1 1],'linewidth',2);plot(pt_i,pt_j,'ms');hold off;
  subplot(212);plot(1:sz(4),squeeze(tmp(pt_i,pt_j,1,:)));grid;
        title(['Point [' num2str(grid_nc{1}.XC(pt_i,pt_j)) 'E,' num2str(grid_nc{1}.YC(pt_i,pt_j)) 'N], within West Spitsbergen Current']);
        ylabel('surface temperature (degC)');xlabel('month since Jan/2002');
  subplot(222);imagesc(iy,-[1:sz(3)],squeeze(tmp(ix,iy,:,192))');set(gca,'ydir','normal','xdir','reverse');colorbar;grid;
	title('Section across Fram Strait for Dec/2017');xlabel('point west to east');ylabel('depth level');

%put contour of velocity, positive is into the Arctic
  hold on;[aa,bb]=contour(iy,-[1:sz(3)],squeeze(tmpu(ix,iy,:,192))',[-.01 -.005]);hold off;set(bb,'color',.9.*[1 1 1],'linestyle','--','linewidth',2);
  hold on;[aa,bb]=contour(iy,-[1:sz(3)],squeeze(tmpu(ix,iy,:,192))',[ .005  .01]);hold off;set(bb,'color',.9.*[1 1 1],'linestyle','-','linewidth',2);

figure(11); %this figure is plotted from the function plot_aste_faces.m
fpr=['./figure_plot_aste_faces.png'];print(fpr,'-dpng');

figure(12);set(gcf,'paperunits','inches','paperposition',[0 0 10 8]);
fpr=['./figure_example_lookat_nc_ASTEr1.png'];print(fpr,'-dpng');
