% Katja - 
%  WARNING: the rnt_oa3d_bdc.m code uses a max file that needs to be compiled. 
% To compile it go in the rnt toolbox within matlab and type 
% >> mex rnt_oa3d_bruce_version_mex.f

%ssh -p 9822 edl@memg.ocean.dal.ca

% 
% 3D Mapping Tutorial using the RNT toolbox
%
%   TUTORIAL 4    3D Objective Mapping (or Optimal Interpolation)
%   TUTORIAL 5    Mapping 3D volumes between s-grid (e.g. nesting)
%
% RNT - E. Di Lorenzo (edl@gatetch.edu)



%   TUTORIAL 4    3D Objective Mapping (or Optimal Interpolation)
% ------------------------------------------------------------------

% Load the outout from a North Pacific Grid
grdin=rnt_gridload('nena');
ctl=rnt_timectl( {'his_nena_856_0306.nc'},'ocean_time');
tempin=rnt_loadvar(ctl,1,'temp');
zeta=rnt_loadvar(ctl,1,'zeta');
grdin.zr=rnt_setdepth(zeta,grdin);
mask3din=repmat(grdin.maskr, [1 1 grdin.N]);

% Load destination grid NEPD
grd=rnt_gridload('mabgom4');
grd.zr=rnt_setdepth(0,grd);
mask3d=repmat(grd.maskr, [1 1 grd.N]);

% Build or Load the PMAP
if ~exist('pmap_nena_2_mabgom4.mat')
  [pmap]=rnt_oapmap(grdin.lonr, grdin.latr, grdin.maskr, grd.lonr, grd.latr, 15);
  save pmap_nena_2_mabgom4 pmap
else
  load pmap_nena_2_mabgom4.mat
end

% Do the 3D OA Mapping
a=10; b=10;   % decorrelation length scale in lon, lat. Keep large in this case.

% If the volume of the CHILD is NOT fully contained in the PARENT use this:
tic
[temp2,temp2_err]=rnt_oa3d_bdc(grdin.lonr, grdin.latr, grdin.zr, tempin.*mask3din, ...
                                  grd.lonr, grd.latr, grd.zr ,a,b,pmap);
toc
temp2=temp2.*mask3d;
temp2(isnan(temp2))=0;

% NOTE: if you do velocities you need to rotate them before or after.

%   TUTORIAL 5    Mapping 3D volumes between s-grid (e.g. nesting)
% ------------------------------------------------------------------

% Let us check how good we expect the 3D Mapping 
% bring topography to from source to destination grid
% Study the geometry of child grid.
h=rnt_griddata(grdin.lonr, grdin.latr, grdin.h.*grdin.maskr, grd.lonr, grd.latr, 'linear');

% check were the destination volume is deeper than the source
hdiff=h-grd.h;
in=find(isnan(grd.maskr));
hdiff(in)=0;

% make a map of the regions where 3D mapping is not accurate
mfig;
pcolor(grd.lonr, grd.latr, hdiff); shading flat; caxis([-300 0]); colorbar
hold on



% Check the result

OPT.interp = 'cubic';
OPT.res    = 8;
x=[-74.1 -76.3]; y=[33 35.15];

hold on;
plot(x,y,'color','k','linewidth',3); 

mfig;  colormap(getpmap(5));
subplot(2,1,1)
rnt_section(grdin.lonr,grdin.latr,grdin.zr,tempin.*mask3din,x,y,OPT);
title 'PARENT GRID'
ax=caxis;
subplot(2,1,2)
rnt_section(grd.lonr,grd.latr,grd.zr,temp2.*mask3d,x,y,OPT);
title 'CHILD GRID'
%caxis(ax); colorbar

mfig;
subplot(2,1,1);
rnt_plc2(tempin(:,:,end),grdin, 0,7,0,0);
rnt_gridbox(grd,'b');
ax=caxis;
subplot(2,1,2);
rnt_plc2(temp2(:,:,end),grd, 0,7,0,0);
caxis(ax); colorbar