interp_1_to_2=interpolate_data(nc1, nc2, data1, data2)

%interpolate_data_xlongitude_ylatitude.m:
%interpolates the data from the modified netcdf file onto a common grid
%This file features a fragment of code to get rid of the NaNs as a result of the 
%difference in grids during the interpolation.
%The KEY IDEA is that all variables with '1' get interpolated into '2's grid.
%This particular code also prepares the data such that the interpolation does not 
%introduce any NaNs.
%written by Max Kaznady from Professor Paul Kushner's hints, May 2005

function interp_1_to_2=interpolate_data(nc1, nc2, data1, data2)

[r1,c1]=size(data1);
[r2,c2]=size(data2);

%get variable vectors from files
lon1=nc1{'lon'}(:);
lat1=nc1{'lat'}(:);
lon2=nc2{'lon'}(:);
lat2=nc2{'lat'}(:);

%preparing data for interpolation
%for each row, add some columns for wraparound

%NOTE: because of absolute values in the code, the interpolation
%is symmetric for large to small AND small to large grids... but, 
%this is still in the testing stage

%if interpolating onto the smaller grid
%if (r1>r2 | c1>c2) 
    
dx=abs(abs(nc1{'lon'}(2))-abs(nc1{'lon'}(1)));
%find Nshift - smallest integer possible such that Nshift * dx > a
true=1;
Nshift=0;
a=abs(abs(nc1{'lon'}(1))-abs(nc2{'lon'}(1)));
while (true)
    Nshift=Nshift+1;
    if (Nshift*dx>a & ~(dx==0))
        true=0;
    end            
end
%Nshift=c1;
%now form a new extended variable
extended_data=ones(r1, Nshift+c1);
extended_lon=ones(1,Nshift+c1);
for row=1:r1,
    for shift=1:Nshift,
        extended_data(row, shift)=data1(row, shift)-(Nshift-shift+1)*dx;
    end           
    extended_data(row, Nshift+1:end)=data1(row, :);
end
for shift=1:Nshift,
    extended_lon(1,shift)=nc1{'lon'}(shift)-(Nshift-shift+1)*dx;
end 
extended_lon(1,(Nshift+1):end)=nc1{'lon'}(:);          
data1=extended_data;
lon1=extended_lon;

%Renew the size.
[r1,c1]=size(data1);            

dx=abs(abs(nc1{'lat'}(2))-abs(nc1{'lat'}(1)));
%find Nshift - smallest integer possible such that Nshift * dx > a
true=1;
Nshift=0;
a=abs(abs(nc1{'lat'}(1))-abs(nc2{'lat'}(1)));
while (true)
    Nshift=Nshift+1;
    if (Nshift*dx>a & ~(dx==0))
        true=0;
    end            
end
%Nshift=r1;
%now form a new extended variable
extended_data=ones(r1+Nshift, c1);
extended_lat=ones(1,r1+Nshift);
for col=1:c1,
    for shift=1:Nshift,
        extended_data(shift, col)=data1(shift, col)-(Nshift-shift+1)*dx;
    end           
    extended_data(Nshift+1:end, col)=data1(:, col);
end
for shift=1:Nshift,
    extended_lat(1,shift)=nc1{'lat'}(shift)-(Nshift-shift+1)*dx;
end 
extended_lat(1,(Nshift+1):end)=nc1{'lat'}(:);          
data1=extended_data;
lat1=extended_lat;

%  size(data1)
%  size(lat1)
%  size(lon1)

%if interpolating onto the larger grid
%else

%end    

%perform interpolation
data1=transpose(data1);
%dupclicate data points result in warnings
warning off
[lat1_mg, lon1_mg]=meshgrid(lat1, lon1);
[lat2_mg, lon2_mg]=meshgrid(lat2, lon2);

%final step
interp_1_to_2=transpose(griddata(lat1_mg, lon1_mg, data1, lat2_mg, lon2_mg));

warning on
return;