!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              REMOVE_POINT    May 2008                                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine remove_point (iii,surface,flag0,surface0)

use definitions
!use mpi

implicit none

include 'mpif.h'

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! This routine removes the node iii from the surface passed as argument.
! It especially ensures that the convex hull is kept intact, and repairs
! the triangulation (it will be Delaunay checked later) 
! Bear in mind that if flag0=0, then the surface0 arrays are deallocated, 
! so all operations on these arrays have to be toggled off. 

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

integer iii,flag0
type(sheet) surface,surface0

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

logical already_in_list,iii_found,jjj_found
double precision, dimension(:), allocatable :: new_x,new_y,new_z,new_r,new_s
double precision, dimension(:), allocatable :: new_xn,new_yn,new_zn
double precision, dimension(:), allocatable :: new_x0,new_y0,new_z0,new_r0,new_s0
double precision, dimension(:), allocatable :: new_xn0,new_yn0,new_zn0
double precision :: xiii,yiii,ziii,xjjj,yjjj,zjjj,dist,distt
integer, dimension(:,:), allocatable :: new_icon
integer,dimension(:),  allocatable :: list_test
integer new_nt, new_nsurface
integer list(20),nb_in_list,iptcl,nb_test,jjj
integer tr_counter,i,j,k,kk,counter,nside
integer ierr,iproc,nproc

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
      
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

!nside=2**surface%levelt+1
nside=2**surface%levelt+2   !opla

if (iproc.eq.0) write(8,*) '=================================================='
if (iproc.eq.0) write(8,*) 'removing point iii=',iii

!===============================================================================
!=====[find how many triangles have point iii as a vertex]======================
!===============================================================================

nb_in_list=1
list(nb_in_list)=iii

tr_counter=0
do i=1,surface%nt
   do k=1,3
      if (surface%icon(k,i)==iii) then 
         tr_counter=tr_counter+1
         do kk=1,3
            iptcl=surface%icon(kk,i)
            already_in_list=.false.
            do j=1,nb_in_list
               if (list(j)==iptcl) already_in_list=.true.
            end do
            if (.not.already_in_list) then
               nb_in_list=nb_in_list+1
               list(nb_in_list)=iptcl
            end if
         end do
      end if
   end do
end do

if(iproc.eq.0) write(8,*) 'nb of neighbouring triangles',tr_counter
if(iproc.eq.0) write(8,*) 'nb_in_list=',nb_in_list
if(iproc.eq.0) write(8,*) 'list of iii neighbours:',list(1:nb_in_list)

!===============================================================================
!=====[build list of neighbouring vertices that are not on convex hull]=========
!===============================================================================

nb_test=count(list(1:nb_in_list)>(4*nside-4))
allocate(list_test(nb_test))
list_test=-10
counter=0
do i=1,nb_in_list
   j=list(i)
   if(j>4*nside-4) then
       counter=counter+1
       list_test(counter)=j
   end if
end do

if(iproc.eq.0) write(8,*) 'nb_test=',nb_test
if(iproc.eq.0) write(8,*) 'list test=',list_test

if(counter/=nb_test) call stop_run('pb in test')

!===============================================================================
!=====[find vertex closest to iii]==============================================
!===============================================================================

dist=10.d0
do i=1,nb_test
   j=list_test(i)
   xiii=surface%x(iii)
   yiii=surface%y(iii)
   ziii=surface%z(iii)
   xjjj=surface%x(j)
   yjjj=surface%y(j)
   zjjj=surface%z(j)
   distt=sqrt((xiii-xjjj)**2+(yiii-yjjj)**2+(ziii-zjjj)**2) 
   if (iii/=j .and. distt<dist) then
      jjj=j
      dist=distt
   end if
end do

if (iproc.eq.0) write(8,*) 'closest point to iii is jjj=',jjj

!===============================================================================
!=====[carry out point removal]=================================================
!===============================================================================

new_nsurface=surface%nsurface-1
new_nt=surface%nt-2

allocate(new_icon(3,new_nt))
allocate(new_x(new_nsurface))
allocate(new_y(new_nsurface))
allocate(new_z(new_nsurface))
allocate(new_xn(new_nsurface))
allocate(new_yn(new_nsurface))
allocate(new_zn(new_nsurface))
allocate(new_r(new_nsurface))
allocate(new_s(new_nsurface))

if (flag0==1) then
allocate(new_x0(new_nsurface))
allocate(new_y0(new_nsurface))
allocate(new_z0(new_nsurface))
allocate(new_xn0(new_nsurface))
allocate(new_yn0(new_nsurface))
allocate(new_zn0(new_nsurface))
allocate(new_r0(new_nsurface))
allocate(new_s0(new_nsurface))
end if

!=====[update icon]=====

tr_counter=0
do i=1,surface%nt
   iii_found=.false.
   jjj_found=.false.
   do k=1,3
      if (surface%icon(k,i)==iii) iii_found=.true.
      if (surface%icon(k,i)==jjj) jjj_found=.true.
   end do
   if (.not.(iii_found .and. jjj_found)) then
      tr_counter=tr_counter+1
      new_icon(1:3,tr_counter)=surface%icon(1:3,i)
      do k=1,3
         if (new_icon(k,tr_counter)==iii) new_icon(k,tr_counter)=jjj
         if (new_icon(k,tr_counter)>iii) new_icon(k,tr_counter)=new_icon(k,tr_counter)-1
      end do
   end if
end do

!=====[update coordinates]=====

counter=0
do i=1,surface%nsurface
   if (i/=iii) then
      counter=counter+1

      new_x  (counter)=surface%x(i)
      new_y  (counter)=surface%y(i)
      new_z  (counter)=surface%z(i)
      new_xn (counter)=surface%xn(i)
      new_yn (counter)=surface%yn(i)
      new_zn (counter)=surface%zn(i)
      new_r  (counter)=surface%r(i)
      new_s  (counter)=surface%s(i)

      if (flag0==1) then
      new_x0 (counter)=surface0%x(i)
      new_y0 (counter)=surface0%y(i)
      new_z0 (counter)=surface0%z(i)
      new_xn0(counter)=surface0%xn(i)
      new_yn0(counter)=surface0%yn(i)
      new_zn0(counter)=surface0%zn(i)
      new_r0 (counter)=surface0%r(i)
      new_s0 (counter)=surface0%s(i)
      end if

   end if
end do

deallocate(surface%x)
deallocate(surface%y)
deallocate(surface%z)
deallocate(surface%xn)
deallocate(surface%yn)
deallocate(surface%zn)
deallocate(surface%icon)
deallocate(surface%r)
deallocate(surface%s)

if (flag0==1) then
deallocate(surface0%x)
deallocate(surface0%y)
deallocate(surface0%z)
deallocate(surface0%xn)
deallocate(surface0%yn)
deallocate(surface0%zn)
deallocate(surface0%r)
end if

surface%nt=new_nt
surface%nsurface=new_nsurface

allocate(surface%x(surface%nsurface))
allocate(surface%y(surface%nsurface))
allocate(surface%z(surface%nsurface))
allocate(surface%xn(surface%nsurface))
allocate(surface%yn(surface%nsurface))
allocate(surface%zn(surface%nsurface))
allocate(surface%icon(3,surface%nt))
allocate(surface%r(surface%nsurface))
allocate(surface%s(surface%nsurface))

if (flag0==1) then
surface0%nsurface=new_nsurface
allocate(surface0%x(surface0%nsurface))
allocate(surface0%y(surface0%nsurface))
allocate(surface0%z(surface0%nsurface))
allocate(surface0%xn(surface0%nsurface))
allocate(surface0%yn(surface0%nsurface))
allocate(surface0%zn(surface0%nsurface))
allocate(surface0%r(surface0%nsurface))
allocate(surface0%s(surface0%nsurface))
end if

surface%x=new_x
surface%y=new_y
surface%z=new_z
surface%xn=new_xn
surface%yn=new_yn
surface%zn=new_zn
surface%r=new_r
surface%s=new_s
surface%icon=new_icon

if (flag0==1) then
surface0%x=new_x0
surface0%y=new_y0
surface0%z=new_z0
surface0%xn=new_xn0
surface0%yn=new_yn0
surface0%zn=new_zn0
surface0%r=new_r0
surface0%s=new_s0
end if

deallocate(new_x)
deallocate(new_y)
deallocate(new_z)
deallocate(new_icon)
deallocate(new_xn)
deallocate(new_yn)
deallocate(new_zn)
deallocate(new_r)
deallocate(new_s)

if (flag0==1) then
deallocate(new_x0)
deallocate(new_y0)
deallocate(new_z0)
deallocate(new_xn0)
deallocate(new_yn0)
deallocate(new_zn0)
deallocate(new_r0)
deallocate(new_s0)
end if

deallocate(list_test)

end subroutine

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------

