function [v] = ll_regrid(loni,lati,vi,lon,lat)

% Function [v] = ll_regrid(loni,lati,vi,lon,lat)
%
% INPUTS
% loni   matrix of input lon coords on the interval [0,360)
% lati   matrix of input lat coords on the interval [0,180)
% vi     matrix of values on the lati,loni points
% lon    matrix of output lon coords on the interval [0,360)
% lat    matrix of output lat coords on the interval [0,180)
%
% OUTPUTS
% v      matrix of values on the lat,lon points
%
% This function implements a very simple regridding method for
% projecting scalars defined on general high-res grids
% (eg. cubesphere) onto lower-res Lat--Lon grids.

v = [];

if nargin ~= 5
  disp('Error: there must be exactly 5 arguments!');
  return
end

%  Mandatory arguments
if (prod(size(lati)) ~= prod(size(loni))) ...
    || (prod(size(lati)) ~= prod(size(vi))) 
  disp(['Error: the lati, loni, and vi must all ' ...
        'be the same size']);
  return
end

%  lat = 0:10:350
%  lon = -80:10:80

% slon = sort(lon(:));
% dlon = diff(slon);
% lons = slon(1)-dlon(1)/2;
% lone = slon(length(slon))+dlon(length(dlon))/2;
% mlon = [ lons ; lons+cumsum(dlon) ; lone ];
% klon = [ 1 ; [1:length(dlon)]' ; length(lon) ];

% slat = sort(lat(:));
% dlat = diff(slat);
% lats = slat(1)-dlat(1)/2;
% late = slat(length(slat))+dlat(length(dlat))/2;
% mlat = [ lats ; lats+cumsum(dlat) ; late ];
% klat = [ 1 ; [1:length(dlat)]' ; length(lat) ];

%  slon = sort(lon);
%  slat = sort(lat);
%  klon = [1:length(slon)];
%  klat = [1:length(slat)];

alon = min(lon);
blon = max(lon);
dlon = blon - alon;
nlon = length(lon);

alat = min(lat);
blat = max(lat);
dlat = blat - alat;
nlat = length(lat);

vsum = zeros(length(lon),length(lat));
vcnt = zeros(length(lon),length(lat));

for kk = 1:(prod(size(vi)))

  %  ii = interp1(slon,klon,loni(kk),'nearest');
  %  jj = interp1(slat,klat,lati(kk),'nearest');
  
  flon = (loni(kk) - alon) / dlon;
  flat = (lati(kk) - alat) / dlat;

  if  (flon >= 0) && (flon <= 1) && (flat >= 0) && (flat <= 1)
    ii = 1 + round(flon*(dlon-1));
    jj = 1 + round(flat*(dlat-1));
    vcnt(ii,jj) = vcnt(ii,jj) + 1;
    vsum(ii,jj) = vsum(ii,jj) + vi(kk);
  end
  
end

inz = find(vcnt ~= 0);
vsum(inz) = vsum(inz) ./ vcnt(inz);

iz = find(vcnt == 0);
vsum(iz) = NaN;

v = vsum;