function oct_bry_interp(zbryname,lon,lat,seas_datafile,ann_datafile,...
                    dataname,vname,obcndx,Roa);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%  function oct_bry_interp(zbryname,lon,lat,seas_datafile,ann_datafile,...
%                      dataname,vname,obcndx,Roa);
% 
%  Interpole data for the lateral boundaries (bry_file) along 
%  horizontal z levels.
% 
%  Further Information:  
%  http://www.croco-ocean.org
%  
%  This file is part of CROCOTOOLS
%
%  CROCOTOOLS is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published
%  by the Free Software Foundation; either version 2 of the License,
%  or (at your option) any later version.
%
%  CROCOTOOLS is distributed in the hope that it will be useful, but
%  WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
%  MA  02111-1307  USA
%
%  Copyright (c) 2001-2006 by Pierrick Penven 
%  e-mail:Pierrick.Penven@ird.fr  
%
%  Updated    5-Oct-2006 by Pierrick Penven (test for negative salinity)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Infos
% dl=1.6;  %-> good for WOA2009
% dl=2.0;   %-> good for CARS2009  
%
[M,L]=size(lon);
%
% set the default value if no data
%
default=NaN;
%
% Read in the datafile
%
ncidseas = netcdf.open(seas_datafile, 'NC_NOWRITE');
X=netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, 'X'));
Y=netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, 'Y'));
Zseas=-netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, 'Z'));
T=netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, 'T'));
tlen=length(T);
Nzseas=length(Zseas);
%
% get the boundary position
%
if obcndx==1
%
% Southern boundary
% 
  icroco=(1:L);
  jcroco=1;
elseif obcndx==2
%
% Eastern boundary
% 
  icroco=L;
  jcroco=(1:M);
elseif obcndx==3
%
% Northern boundary
% 
  icroco=(1:L);
  jcroco=M;
elseif obcndx==4
%
% Western boundary
% 
  icroco=1;
  jcroco=(1:M);
end
%
lon=lon(jcroco,icroco);
lat=lat(jcroco,icroco);
%
% get a data subgrid (dependant of the OBC used)
%
dl=1.6;  %-> good for WOA2009
%dl=2.0;   %-> good for CARS2009  

lonmin=min(min(lon))-dl;
lonmax=max(max(lon))+dl;
latmin=min(min(lat))-dl;
latmax=max(max(lat))+dl;
%
j=find(Y>=latmin & Y<=latmax);
i1=find(X-360>=lonmin & X-360<=lonmax);
i2=find(X>=lonmin & X<=lonmax);
i3=find(X+360>=lonmin & X+360<=lonmax);
x=cat(1,X(i1)-360,X(i2),X(i3)+360);
y=Y(j);
%
% Open the Z-boundary file
%
ncid = netcdf.open(zbryname, 'NC_WRITE');
Z=-netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'Z'));
Nz=length(Z);
%
% Check the time
tbry=netcdf.getVar(ncid, netcdf.inqVarID(ncid, 'bry_time')); 
T=T*30; % if time in month in the dataset !!!
%
%% Comment => check only if woa_time=(15:30:345) and woa_cycle=360
% if tbry~=T ;    
%     error(['time mismatch  tbry = ',num2str(tbry),...
%            '  t = ',num2str(T)])
% end
%
% Read the annual dataset
%
if Nz > Nzseas
    ncidann = netcdf.open(ann_datafile, 'NC_NOWRITE');
    zann=-netcdf.getVar(ncidann, netcdf.inqVarID(ncidann, 'Z'));
    if (Z~=zann)
        error('Vertical levels mismatch')
    end
%
% Interpole the annual dataset on the horizontal CROCO grid
%
    disp(' Ext tracers: horizontal interpolation of the annual data')
    if Zseas~=zann(1:length(Zseas)) 
        error('vertical levels dont match')
    end
    dims=size(lon);
    datazgrid=zeros(Nz,length(lon));
    missval=netcdf.getAtt(ncidann, netcdf.inqVarID(ncidann, dataname), 'missing_value');
    for k=Nzseas+1:Nz
        if ~isempty(i2)
            data=squeeze(netcdf.getVar(ncidann, netcdf.inqVarID(ncidann, dataname)));
        else
            data=[];
        end
        if ~isempty(i1)
            data=cat(2,squeeze(netcdf.getVar(ncidann, netcdf.inqVarID(ncidann, dataname))),data);
        end
        if ~isempty(i3)
            data=cat(2,data,squeeze(netcdf.getVar(ncidann, netcdf.inqVarID(ncidann, dataname))));
        end
        data=oct_get_missing_val(x,y,data,missval,Roa,default);
        if dims(1)==1
            datazgrid(k,:)=interp2(x,y,data,lon,lat,'cubic');
        else
            datazgrid(k,:)=(interp2(x,y,data,lon,lat,'cubic'))';    
        end
    end
    netcdf.close(ncidann);
end

%Else read seasonal datafile 
%
% interpole the seasonal dataset on the horizontal croco grid
%
disp([' Ext tracers: horizontal interpolation of the seasonal data'])
%
% loop on time
%
missval=netcdf.getAtt(ncidseas, netcdf.inqVarID(ncidseas, dataname), 'missing_value');
for l=1:tlen
    %for l=1:1
    disp(['time index: ',num2str(l),' of total: ',num2str(tlen)])
    dims=size(lon);
    if Nz <= Nzseas
        %    datazgrid=zeros(Nz,M,L);
        datazgrid=zeros(Nz,length(lon));
    end
    for k=1:min([Nz Nzseas])
        if ~isempty(i2)
            data=squeeze(netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, dataname)));
        else
            data=[];
        end
        if ~isempty(i1)
            data=cat(2,squeeze(netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, dataname))),data);
        end
        if ~isempty(i3)
            data=cat(2,data,squeeze(netcdf.getVar(ncidseas, netcdf.inqVarID(ncidseas, dataname))));
        end
        data=oct_get_missing_val(x,y,data,missval,Roa,default);
        if dims(1)==1
            datazgrid(k,:)=interp2(x,y,data,lon,lat,'cubic');
        else
            datazgrid(k,:)=(interp2(x,y,data,lon,lat,'cubic'))';    
        end
    end
%
% Test for salinity (no negative salinity !)
%
    if strcmp(vname,'salt_south') | strcmp(vname,'salt_north') | ...
            strcmp(vname,'salt_east') | strcmp(vname,'salt_west')
        disp('salinity test')
        datazgrid(datazgrid<2)=2;
    end
%
    netcdf.putVar(ncid, netcdf.inqVarID(ncid, vname), l,:,:,:-1, 1, datazgrid);  % [conv] 0-based
end
netcdf.close(ncid);
netcdf.close(ncidseas);
return