%function interpolate_obs_to_model_grid

clear

% load obs original grid
origin_lon=double(ncread('/work/uo0122/u241089/DATA/INSITU/TS/EN4/DATA/file.000','lon')); 
origin_lat=double(ncread('/work/uo0122/u241089/DATA/INSITU/TS/EN4/DATA/file.000','lat')); 
origin_dep=double(ncread('/work/uo0122/u241089/DATA/INSITU/TS/EN4/DATA/file.000','depth'));

new_lon_EN=origin_lon;
new_lon_EN(new_lon_EN>180)=new_lon_EN(new_lon_EN>180)-360;
[y,x]=meshgrid(origin_lat,new_lon_EN);
%x_new=x;
%x_new(361:360+50,:)=x(1:50,:);
%x_new(1:50,:)=[];


% build obs original grid 3d matrices
orig_lon=repmat(x,[1 1 42]);
orig_lat=repmat(y,[1 1 42]);
orig_dep=repmat(origin_dep,[1 173 360]);orig_dep=permute(orig_dep,[3 2 1]);


% load obs.(model) destination grid
destin_lon=double(ncread('/work/uo0122/u241089/ATL/ATL24/grid.cdf','xc')); 
destin_lat=double(ncread('/work/uo0122/u241089/ATL/ATL24/grid.cdf','yc'));
destin_dep=double(ncread('/work/uo0122/u241089/ATL/ATL24/grid_vert.cdf','ZT'));

% build obs.(model) destination grid 3d matrices
%[dest_lat,dest_lon,dest_dep]=meshgrid(destin_lat,destin_lon,destin_dep);

dest_lon=repmat(destin_lon,[1 1 100]);
dest_lat=repmat(destin_lat,[1 1 100]);
dest_dep=repmat(destin_dep,[1 3456 3072]);dest_dep=permute(dest_dep,[3 2 1]);


% first find grid non-NaN values
good_values=find(~isnan(orig_lon));

% load model (obs) field to be interpolated
load /work/uo0122/u241272/Compare_EN4_GECCO2/EN4_salt.mat

 
saltEN_in_ATL4km= zeros(3072,3456,100);

% loop over time
for n=637:744
tic
n

% select time step
orig_field=saltEN(:,:,:,n);

% build interpolant
F=TriScatteredInterp(orig_lon(good_values),orig_lat(good_values),orig_dep(good_values),orig_field(good_values),'linear');

% evaluate field in the destination grid
dest_field=F(dest_lon,dest_lat,dest_dep);

% store in matrix
saltEN_in_ATL4km=dest_field;

% save in matlab format
save(['/work/uo0122/u241272/Compare_EN4_GECCO2/EN4_in_ATL4grid/EN4_salt_in_ATL4km_grid_',int2str(n-636),'.mat'],'saltEN_in_ATL4km', '-v7.3'); 


% clear temporary variables
clear orig_field F dest_field saltEN_in_ATL4km

toc
end %n

%figure(1);pcolor(squeeze(saltEN_on_ATL4km(:,:,20))'); shading flat;
