!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              MOVE_SURFACE    Nov. 2007                                       |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine move_surface (surface,surface0,flag0,ov,dt,params,istep,is)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! subroutine to move a set of particles on a surface from a velocity
! field known at the nodes of an octree

! surface is the surface on which the particles need to be moved
! surface0 is the surface in its initial geometry
! flag0 is 1 if in midpoint configuration and 0 if at the end of the timestep
! if flag0=0, then the surface0 arrays are deallocated!!
! octree is the velocity octree
! noctree is its size
! unode,vnode and wnode are the components of the velocity
! nnode is the number of nodes on the octree
! icon_octree is the connectivity matrix on the octree
! nleaves is the number of leaves in the octree
! dt is the time step
! niter_move is the number of iterations used to move the particles
! it is read in input.txt

! here we use a simple method that performs a given (niter_move)
! number of iterations in a mid_point algorithm

! we also advect the normals by using the velocity gradient computed
! from the finite element shape function derivatives

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

use definitions

implicit none

type (sheet) surface
type (sheet) surface0
type (octreev) ov
integer flag0
double precision dt
type(parameters) params
integer istep,is
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer i,ip,iter,nremove,nedge,nrem
integer i1,i2,i3,ij,k,err,nside
integer iproc,nproc,ierr
integer,dimension(:),allocatable::new,old
type (sheet) surfacep,surfacep0
double precision u,v,w,rat,eps
double precision x0,y0,z0,xi,yi,zi,x,y,z
double precision,dimension(:),allocatable::xbuf,ybuf,zbuf
double precision,dimension(:),allocatable::xnbuf,ynbuf,znbuf
double precision,dimension(:),allocatable::pos,zzz
double precision ux,uy,uz,vx,vy,vz,wx,wy,wz,xyzn
double precision xn0,yn0,zn0,xyn,theta,phi
double precision xn1,yn1,zn1,xn2,yn2,zn2
double precision xln1,yln1,zln1,xln2,yln2,zln2
double precision x1,x2,x3,y1,y2,y3,z1,z2,z3,xne,yne,zne
logical remove
character*72  :: shift

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

INCLUDE 'mpif.h'
      
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

shift=' '

allocate (xbuf(surface%nsurface),stat=err)  ; if (err.ne.0) call stop_run ('Error alloc xbuf  in move_surface$')
allocate (ybuf(surface%nsurface),stat=err)  ; if (err.ne.0) call stop_run ('Error alloc ybuf  in move_surface$')
allocate (zbuf(surface%nsurface),stat=err)  ; if (err.ne.0) call stop_run ('Error alloc zbuf  in move_surface$')
allocate (xnbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc xnbuf in move_surface$')
allocate (ynbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ynbuf in move_surface$')
allocate (znbuf(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc znbuf in move_surface$')

xbuf=0.d0
ybuf=0.d0
zbuf=0.d0
xnbuf=0.d0
ynbuf=0.d0
znbuf=0.d0

rat=1.d0-1.d-10
eps=1.d-10

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

do i=iproc+1,surface%nsurface,nproc
   u=0.d0
   v=0.d0
   w=0.d0
   x0=surface%x(i)
   y0=surface%y(i)
   z0=surface%z(i)
   x=x0
   y=y0
   z=z0
   do iter=1,params%niter_move
      xi=(x0+x)/2.d0
      yi=(y0+y)/2.d0
      zi=(z0+z)/2.d0
      call octree_interpolate_three (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
                                    ov%nnode,xi,yi,zi,ov%unode,u,ov%vnode,v,ov%wnode,w)

      if (i==1 .or. i==nside .or. i==2*nside-1 .or. i==3*nside-2) then !corner points
         z=z0+dt*w
      elseif (i>1 .and. i<nside) then                                  !y=0 edge
         x=x0+dt*u
         z=z0+dt*w
      elseif (i>nside .and. i<2*nside-1) then                          !x=1 edge 
         y=y0+dt*v
         z=z0+dt*w
      elseif (i>2*nside-1 .and. i<3*nside-2) then                      !y=1 edge
         x=x0+dt*u
         z=z0+dt*w
      elseif (i>3*nside-2 .and. i<4*nside-3) then                      !x=0 edge
         y=y0+dt*v
         z=z0+dt*w
      else 
         x=x0+dt*u
         y=y0+dt*v
         z=z0+dt*w
      end if
   enddo
   xbuf(i)=x
   ybuf(i)=y
   zbuf(i)=z
   if (params%normaladvect) then
      call octree_interpolate_three_derivative (3,ov%octree,ov%noctree,ov%icon,ov%nleaves, &
                                               ov%nnode,xi,yi,zi, &
                                               ov%unode,u,ux,uy,uz, &
                                               ov%vnode,v,vx,vy,vz, &
                                               ov%wnode,w,wx,wy,wz)
      xn0=surface%xn(i)
      yn0=surface%yn(i)
      zn0=surface%zn(i)
      xyn=xn0**2+yn0**2
      if (xyn.gt.0.d0) xyn=sqrt(xyn)
      theta=atan2(xyn,zn0)
      phi=0.d0
      if (sin(theta).ne.0.d0) phi=atan2(yn0,xn0)
      xn1=cos(theta)*cos(phi)
      yn1=cos(theta)*sin(phi)
      zn1=-sin(theta)
      xn2=-sin(phi)
      yn2=cos(phi)
      zn2=0.d0
      xln1=ux*xn1+uy*yn1+uz*zn1
      yln1=vx*xn1+vy*yn1+vz*zn1
      zln1=wx*xn1+wy*yn1+wz*zn1
      xln2=ux*xn2+uy*yn2+uz*zn2
      yln2=vx*xn2+vy*yn2+vz*zn2
      zln2=wx*xn2+wy*yn2+wz*zn2
      xnbuf(i)=xn0+(yn1*zln2-zn1*yln2+yln1*zn2-zln1*yn2)*dt+(yln1*zln2-zln1*yln2)*dt**2
      ynbuf(i)=yn0+(zn1*xln2-xn1*zln2+zln1*xn2-xln1*zn2)*dt+(zln1*xln2-xln1*zln2)*dt**2
      znbuf(i)=zn0+(xn1*yln2-yn1*xln2+xln1*yn2-yln1*xn2)*dt+(xln1*yln2-yln1*xln2)*dt**2
      if (xnbuf(i).gt.1.d50) write(*,*) 'on proc ', iproc, ', xnbuf is large: ', xnbuf(i)
      if (ynbuf(i).gt.1.d50) write(*,*) 'on proc ', iproc, ', ynbuf is large: ', ynbuf(i)
      if (znbuf(i).gt.1.d50) write(*,*) 'on proc ', iproc, ', znbuf is large: ', znbuf(i)
      xyzn=sqrt(xnbuf(i)**2+ynbuf(i)**2+znbuf(i)**2)
      xnbuf(i)=xnbuf(i)/xyzn
      ynbuf(i)=ynbuf(i)/xyzn
      znbuf(i)=znbuf(i)/xyzn
   endif
enddo

surface%x=0.d0
surface%y=0.d0
surface%z=0.d0
call mpi_allreduce (xbuf,surface%x,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (ybuf,surface%y,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
call mpi_allreduce (zbuf,surface%z,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)

if (params%normaladvect) then
   surface%xn=0.d0
   surface%yn=0.d0
   surface%zn=0.d0
   call mpi_allreduce (xnbuf,surface%xn,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
   call mpi_allreduce (ynbuf,surface%yn,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
   call mpi_allreduce (znbuf,surface%zn,surface%nsurface,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)
else
   call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
                         surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
endif

deallocate (xbuf)
deallocate (ybuf)
deallocate (zbuf)
deallocate (xnbuf)
deallocate (ynbuf)
deallocate (znbuf)

! -------------------------------------------------------------------------------
! At this point, all points of the surface (and possibly the normal they carry) 
! have been advected. A little bit of extra work is needed for the points on the 
! convex hull (those at the intersection of the surface with the cube). On each 
! edge of the square, there are nside points of coordinates (x_i,z(x_i)). 
! Because of advection, the x_i are no more equidistant, which after several time 
! steps can lead to one these points to be advected outside the cube and removed, 
! which in turn leads to other problems. In order to avoid this problem, and also 
! to keep the triangulation somewhat regular near the edges, the curve (x_i,z(x_i)) 
! is reinterpolated at equidistant X_i by means of the 'resample' routine in 
! /RESAMPLE/libresample.a, based on linear interpolation (k=1) or spline 
! interpolation (k=3,5). This procedure is carried out only at the end of a time 
! step, i.e. when flag0=0.

allocate(pos(nside))
allocate(zzz(nside))

pos=surface%x(1:nside)
zzz=surface%z(1:nside)
call resample (pos,zzz,nside,1)
surface%x(1:nside)=pos
surface%z(1:nside)=zzz

pos=surface%y(nside:2*nside-1)
zzz=surface%z(nside:2*nside-1)
call resample (pos,zzz,nside,1)
surface%y(nside:2*nside-1)=pos
surface%z(nside:2*nside-1)=zzz

pos=surface%x(3*nside-2:2*nside-1:-1)
zzz=surface%z(3*nside-2:2*nside-1:-1)
call resample (pos,zzz,nside,1)
surface%x(3*nside-2:2*nside-1:-1)=pos
surface%z(3*nside-2:2*nside-1:-1)=zzz

pos(1)=surface%y(1)
zzz(1)=surface%z(1)
pos(2:nside)=surface%y(4*nside-4:3*nside-2:-1)
zzz(2:nside)=surface%z(4*nside-4:3*nside-2:-1)
call resample (pos,zzz,nside,1)
surface%y(4*nside-4:3*nside-2:-1)=pos(2:nside)
surface%z(4*nside-4:3*nside-2:-1)=zzz(2:nside)

deallocate(pos)
deallocate(zzz)


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

if (surface%itype.eq.0) then
   surface%r=surface%x
   surface%s=surface%y
endif

!-------------------------------------------------------------------------------
! We now check whether points have been advected outside the simulation domain 
! (unit cube), and if so, we need to remove them. However, because the triangulation 
! is only computed at startup and further updated, the removal of a point requires 
! the triangulation to be sown back by hand. This will most probably lead to a 
! non-Delaunay triangulation (locally), but this will be taken care of in the 
! check_delaunay subroutine. The first 4*nside-4 points in the surface 
! array correspond to those on the convex hull so that we need not test them.


nrem=0
remove=.true.
do while (remove)
   remove=.false.
   do i=4*nside-3,surface%nsurface
      if (surface%x(i).lt.-eps .or. surface%x(i).gt.1.d0+eps .or. &
          surface%y(i).lt.-eps .or. surface%y(i).gt.1.d0+eps .or. &
          surface%z(i).lt.-eps .or. surface%z(i).gt.1.d0+eps) then
         call remove_point(i,surface,flag0,surface0)
         nrem=nrem+1
         remove=.true.
         exit
      endif
   end do
end do

if (iproc.eq.0 .and. nrem/=0) write(*,'(a,i4,a)') shift//'removing',nrem,' particle from surface'

call check_delaunay (params,surface,is,istep,'delaunay_aft_remove',0)

if (.not.params%normaladvect) then
   call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
                         surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)
end if

if (flag0==0) then
   nedge=surface%nsurface+surface%nt-1
   call DoRuRe_surf_stats (params%doDoRuRe,istep,50,is,surface%nt,surface%nsurface,nedge,-nrem)
end if

if (iproc.eq.0) write (8,*) 'Removal complete'

end subroutine move_surface
