Subversion Repositories lagranto.ecmwf

function out = get_index (vert,position,grid);

% General index calculation for data grid
% Usage:   out = int_index(vert,coords,grid)   
% Example: grid.xmin = -180;
%          grid.dx   =    1;
%          grid.nx   =  361;
%          grid.ymin =    0;
%          grid.dy   =    1;
%          grid.ny   =   91;
%          out get_index(inp.P.data,[ 500 45 78; 250 56.7 78.4],grid)   

% Get the position arrays
xpos = position(:,3)';
ypos = position(:,2)';
zpos = position(:,1)';

% Get the grid dimension 
nx = size(vert,3);
ny = size(vert,2);
nz = size(vert,1);

% Get the horiozontal grid indices
rindx = (xpos-grid.xmin)/grid.dx+1;
rindy = (ypos-grid.ymin)/grid.dy+1;

% Special treatment if the the <vert> field is two-dimensional
if ( nz == 1 )
   out = [ ones(size(rindy))', rindy', rindx' ];
   return;
end 
    
% Direction of the vertical axis (Pressure: descencding, Height: ascending)
dir = vert(2:nz,:,:) - vert(1:nz-1,:,:);
if ( mean(dir(:)) > 0 ) 
     zpos = -zpos;
     vert = -vert;
end

% Get the weights for the individual grid points
fracx = 1-rindx+floor(rindx);
fracy = 1-rindy+floor(rindy);

% Combine the indices i and j to a common index ij
vert1 = reshape(vert,nz,nx*ny);

% The position array and the vertical index is needed in the same format
zpos1 = kron(ones(nz,1),zpos);
zind1 = kron([1:nz]',ones(size(zpos)));
shif1 = (0:length(zpos)-1) * nz;

% Get the vertical indices at position 00
indx           = floor(rindx);
indy           = floor(rindy);
indxy          = indy+(indx-1)*ny;
isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
indxy(isnok1)  = 1;  
vert           = vert1(:,indxy);
indz           = max((vert>=zpos1).*zind1);
isnok2         = (indz<1) | (indz>=nz);
indz(isnok2)   = 1;
indu           = indz     + shif1;
indo           = indz + 1 + shif1;
diff           = vert(indo)-vert(indu);
isnok3         = (diff==0);
diff(isnok3)   = 1;
rind00         = indz + (zpos-vert(indu)) ./ diff;
isnok4         = isnok1 | isnok2 | isnok3;
rind00(isnok4) = NaN;

% Get the vertical indices at positions 01, 10, 11 (if necessary)
if ( any(fracx ~= 1) | any(fracy ~=1) ) 
    
    indx           = floor(rindx+1);
    indy           = floor(rindy);
    indxy          = indy+(indx-1)*ny;
    isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
    indxy(isnok1)  = 1;  
    vert           = vert1(:,indxy);
    indz           = max((vert>zpos1).*zind1);
    isnok2         = (indz<1) | (indz>=nz);
    indz(isnok2)   = 1;
    indu           = indz     + shif1;  
    indo           = indz + 1 + shif1;
    diff           = vert(indo)-vert(indu);
    isnok3         = (diff==0);
    diff(isnok3)   = 1;
    rind01         = indz + (zpos-vert(indu)) ./ diff;
    isnok4         = isnok1 | isnok2 | isnok3;
    rind01(isnok4) = NaN;

    indx           = floor(rindx);
    indy           = floor(rindy+1);
    indxy          = indy+(indx-1)*ny;
    isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
    indxy(isnok1)  = 1;  
    vert           = vert1(:,indxy);
    indz           = max((vert>zpos1).*zind1);
    isnok2         = (indz<1) | (indz>=nz);
    indz(isnok2)   = 1;
    indu           = indz     + shif1;  
    indo           = indz + 1 + shif1;
    diff           = vert(indo)-vert(indu);
    isnok3         = (diff==0);
    diff(isnok3)   = 1;
    rind10         = indz + (zpos-vert(indu)) ./ diff;
    isnok4         = isnok1 | isnok2 | isnok3;
    rind10(isnok4) = NaN;

    indx           = floor(rindx+1);
    indy           = floor(rindy+1);
    indxy          = indy+(indx-1)*ny;
    isnok1         = (indx<1) | (indx>nx) | (indy<1) | (indy>ny);
    indxy(isnok1)  = 1;  
    vert           = vert1(:,indxy);
    indz           = max((vert>zpos1).*zind1);
    isnok2         = (indz<1) | (indz>=nz);
    indz(isnok2)   = 1;
    indu           = indz     + shif1;  
    indo           = indz + 1 + shif1;
    diff           = vert(indo)-vert(indu);
    isnok3         = (diff==0);
    diff(isnok3)   = 1;
    rind11         = indz + (zpos-vert(indu)) ./ diff;
    isnok4         = isnok1 | isnok2 | isnok3;
    rind11(isnok4) = NaN;

else

    rind01 = rind00;
    rind10 = rind00;
    rind11 = rind00;
    
end
    
% Vectorised interpolation
rindz = fracx     .* fracy     .* rind00 + ...
        (1-fracx) .* fracy     .* rind01 + ...
        fracx     .* (1-fracy) .* rind10 + ...
        (1-fracx) .* (1-fracy) .* rind11;
        
% Set the output of the function
nok = isnan(rindz') | isnan(rindy') | isnan(rindx');
out = [ rindz', rindy', rindx' ];