%function interpolate_model_to_obs_grid

clear

% load model original grid
origin_lon=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/Atl_1deg/monthly_temp_01_2003.nc','longitude');origin_lon(origin_lon==0)=NaN;origin_lon(origin_lon<0)=origin_lon(origin_lon<0)+360;
origin_lat=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/Atl_1deg/monthly_temp_01_2003.nc','latitude');origin_lat(origin_lat==0)=NaN;
origin_dep=ncread('/work/uo0122/u241272/Compare_EN4_GECCO2/Atl_1deg/monthly_temp_01_2003.nc','depth');

% build model original grid 3d matrices
orig_lon=repmat(origin_lon,[1 1 100]);
orig_lat=repmat(origin_lat,[1 1 100]);
orig_dep=repmat(origin_dep,[1 138 122]);orig_dep=permute(orig_dep,[3 2 1]);

% load obs. destination grid
destin_lon=double(ncread('/work/uo0122/u241089/DATA/INSITU/TS/EN4/DATA/file.000','lon'));
destin_lat=double(ncread('/work/uo0122/u241089/DATA/INSITU/TS/EN4/DATA/file.000','lat'));
destin_dep=double(ncread('/work/uo0122/u241089/DATA/INSITU/TS/EN4/DATA/file.000','depth'));

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

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

year=[2003:2011];  %#ok<NBRAK>
tempATL4_in_EN4 = zeros(360,173,42);

% loop over time

for i=1:9
tic
i

for j=1:12
tic
j

if j<10
mnth=['0',int2str(j)];
else
mnth=int2str(j);
end %if


% load model field to be interpolated
tempATL4=double(ncread(['/work/uo0122/u241272/Compare_EN4_GECCO2/Atl_1deg/monthly_temp_',mnth,'_',int2str(year(i)),'.nc'],'tempdeg'));


% select time step
orig_field=tempATL4;

% 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
tempATL4_in_EN4=dest_field;

%save in nc format
nccreate(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.nc'],'longitude','dimensions',{'longitude' 360 'latitude' 173});
nccreate(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.nc'],'latitude','dimensions',{'longitude' 360 'latitude' 173});
nccreate(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.nc'],'depth','dimensions',{'depth' 42});

ncwrite(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.nc'],'longitude',destin_lon);
ncwrite(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.nc'],'latitude',destin_lat);
ncwrite(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.nc'],'depth',destin_dep);


nccreate(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.nc'],'tempATL4_in_EN4','dimensions',{'longitude' 360 'latitude' 173 'depth' 42});
ncwriteatt(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.nc'],'tempATL4_in_EN4','units','psu');

ncwrite(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.nc'],'tempATL4_in_EN4',tempATL4_in_EN4);


% save in matlab format
%save(['/work/uo0122/u241272/Compare_EN4_GECCO2/ATL4km_in_EN4grid/ATL4km_temp_in_EN4grid_',mnth,'_',int2str(year(i)),'.mat'],'tempATL4_in_EN4','-v7.3');


% clear temporary variables
clear orig_field F dest_field tempATL4_in_EN4


toc
end %j

toc
end %i








