       subroutine bdryi(n,tabr,tabz,itype,ri,zi,taub,normal)

c***********************************************************************
c      BDRYI.FOR
c      BounDaRY Interpolation
c      Faro, 19/08/2009 at 20:00
c      Written by Orlando Camargo Rodriguez
c****&******************************************************************

c      Include global (common) variables:

       include 'globli.for'
       
c      Define local variables:
       
       character*2 itype 
       
       integer*4 i,n
       
       real*8 tabr(n),tabz(n)
       real*8 xc(4),yc(4)
       real*8 xp(3),yp(3)
       real*8 xl(2),yl(2)
       real*8 taub(2),normal(2)
       real*8 ri,zi,zr,zrr

c***********************************************************************

c       All interpolation types are tested by readin.for, 
c       so the default option is a flat boundary:

       if (itype.eq.'3P') then
          
          if (ri.le.tabr(2)) then
           i = 1 
          elseif (ri.ge.tabr(n-1)) then
           i = n - 2 
          else 
           call brcket(n,tabr,ri,i)
          end if         
          
          xp(1) = tabr( i )
          xp(2) = tabr(i+1)
          xp(3) = tabr(i+2)
          yp(1) = tabz( i )
          yp(2) = tabz(i+1)
          yp(3) = tabz(i+2)
          
          call bpai1d(xp,yp,ri,zi,zr,zrr)
       
       elseif (itype.eq.'4P') then
          
          if (ri.le.tabr(2)) then
           i = 2 
          elseif (ri.ge.tabr(n-1)) then
           i = n - 2 
          else 
           call brcket(n,tabr,ri,i)
          end if         
          
          xc(1) = tabr(i-1)
          xc(2) = tabr( i )
          xc(3) = tabr(i+1)
          xc(4) = tabr(i+2)
          yc(1) = tabz(i-1)
          yc(2) = tabz( i )
          yc(3) = tabz(i+1)
          yc(4) = tabz(i+2)
          
          call bcui1d(xc,yc,ri,zi,zr,zrr)
       
       elseif (itype.eq.'2P') then

          call brcket(n,tabr,ri,i)
          
          xl(1) = tabr( i )
          xl(2) = tabr(i+1)
          yl(1) = tabz( i )
          yl(2) = tabz(i+1)
     
          call lini1d(xl,yl,ri,zi,zr)

       elseif (itype.eq.'SL') then
       
          xl(1) = tabr(1)
          xl(2) = tabr(2)
          yl(1) = tabz(1)
          yl(2) = tabz(2)
     
          call lini1d(xl,yl,ri,zi,zr)
       
       else 
          
          xl(1) = tabr(1)
          xl(2) = tabr(2)
          yl(1) = tabz(1)
          yl(2) = tabz(1)
     
          call lini1d(xl,yl,ri,zi,zr)
	  
       end if
       
       taub(1) =  cos( atan(zr) ) 
       taub(2) =  sin( atan(zr) )

c      Avoid rounding errors: 
         
       if ( abs(taub(1)).eq.1.0 ) then 
       taub(2) = 0.0 
       end if 
         
       if ( abs(taub(2)).eq.1.0 ) then 
       taub(1) = 0.0 
       end if         
               
       normal(1) = -taub(2) 
       normal(2) =  taub(1)
      
c***********************************************************************
c      Back to main:
c***********************************************************************
       
       return
       
       end

