erosion.old.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              EROSION    Nov. 2006                                            |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine erosion (surface,olsf,is,zerosion)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine will need major improvement in the near future
! it is used to erode the surface "surface" according to a user supplied algorithm
! at the moment it simply "shaves" the top of the model at a given height

! surface are the surface/sheet object
! (the top free surface)
! olsf is a octreelsf object containig the geometry of the current velocity octree/object
! is is the number of the surface (0 playing a special role)
! zerosion is the level of erosion (between 0 and 1)

! Note that all surfaces are required because they can be affected by the erosion process

! for the top surface, we apply the erosion and on the octree
! we define an lsf that is simply the height of the eroded surface
! This allows us, for the other surfaces to find the height of the
! free/eroded surface at any point in space and check whether the
! vertical position of the point is above or below the free surface
! We then simply "erode" the following surfaces using this informatoin

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

use definitions
!use mpi

implicit none

include 'mpif.h'

type (sheet) surface
type (octreelsf) olsf
integer is
double precision zerosion

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer ierr,iproc,nproc,i,nt,nh,nohalt_hull,loc,nnn,nnlist,ntrilist
integer ntmax,nhmax,npmax,nmax,nmode,nvmax,nnpnmax,int_method
double precision,dimension(:),allocatable::field
double precision,dimension(:,:),allocatable::points,centres
integer,dimension(:,:),allocatable::vertices,neighbour
integer,dimension(:),allocatable::hulltriangles,lt_work,ln_work
integer,dimension(:),allocatable::vis_tlist,vis_elist,add_tlist
logical clockwise,out,nnext
double precision eps,zref
integer dmode,err
double precision x,y,f,df(2),ddf(3),v
double precision,dimension(:,:),allocatable::work_d1,work_d2,work_d3
double precision,dimension(:,:),allocatable::work_d4,work_d5,work_d6,work_d7
integer,dimension(:),allocatable::work_i1,work_i3
integer,dimension(:,:),allocatable::work_i2
real,dimension(:),allocatable::work_r2
real,dimension(:,:),allocatable::work_r1
double precision,dimension(:),allocatable::lsf

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

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

if (is.ne.1) then

   do i=1,surface%nsurface
      call octree_interpolate (olsf%octree,olsf%noctree,olsf%icon, &
                               olsf%nleaves,olsf%lsf,olsf%nnode, &
                               surface%x(i),surface%y(i),surface%z(i),f)

!      call octree_interpolate3 (olsf%octree,olsf%noctree,olsf%icon,olsf%nleaves,olsf%lsf,                  &
!                                olsf%x(1:olsf%nnode),olsf%y(1:olsf%nnode),olsf%z(1:olsf%nnode),olsf%nnode, &
!                                surface%x(i),surface%y(i),surface%z(i),f)

      surface%z(i)=min(f,surface%z(i))
   enddo

   return

endif

! first apply erosion to top surface

zref=zerosion

do i=1,surface%nsurface
   surface%z(i)=zref
enddo

! nn_setup (calculates delaunay triangles and other info)

ntmax=surface%nsurface*3
nhmax=surface%nsurface
npmax=surface%nsurface
nnpnmax=100
nmax=3*ntmax+npmax
nvmax=ntmax

allocate (points(2,surface%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc points in erode$')
allocate (field(surface%nsurface),stat=err)   ; if (err.ne.0) call stop_run ('Error alloc field in erode$')
allocate (vertices(3,ntmax),stat=err)         ; if (err.ne.0) call stop_run ('Error alloc vertices in erode$')
allocate (centres(3,ntmax),stat=err)          ; if (err.ne.0) call stop_run ('Error alloc centress in erode$')
allocate (neighbour(3,ntmax),stat=err)        ; if (err.ne.0) call stop_run ('Error alloc neighbour in erode$')
allocate (hulltriangles(nhmax),stat=err)      ; if (err.ne.0) call stop_run ('Error alloc hulltriangles in erode$')
allocate (vis_tlist(nvmax),stat=err)          ; if (err.ne.0) call stop_run ('Error alloc vis_tlist in erode$')
allocate (vis_elist(nvmax),stat=err)          ; if (err.ne.0) call stop_run ('Error alloc vis_elist in erode$')
allocate (add_tlist(nvmax),stat=err)          ; if (err.ne.0) call stop_run ('Error alloc add_tlist in erode$')
allocate (lt_work(ntmax),stat=err)            ; if (err.ne.0) call stop_run ('Error alloc lt_work in erode$')
allocate (ln_work(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ln_work in erode$')

points(1,:)=surface%x
points(2,:)=surface%y
field(:)=surface%z
dmode=-2
nmode=0
clockwise=.true.
eps=1.d-20
eps=0.d0
      
call nn2d_setup (surface%nsurface,ntmax,nhmax,npmax,nnpnmax,nmax, &
                 points,dmode,nmode,clockwise,field,nt,vertices, &
                 centres,neighbour,nh,hulltriangles,nohalt_hull, &
                 loc,nnn,nnlist,ntrilist, &
                 eps,nvmax,vis_tlist,vis_elist,add_tlist, &
                 lt_work,ln_work)

allocate (work_d1(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d1 in erode$')
allocate (work_d2(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d2 in erode$')
allocate (work_d3(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d3 in erode$')
allocate (work_d4(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d4 in erode$')
allocate (work_d5(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d5 in erode$')
allocate (work_d6(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d6 in erode$')
allocate (work_d7(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_d7 in erode$')
allocate (work_i1(nnpnmax),stat=err)   ; if (err.ne.0) call stop_run ('Error alloc work_i1 in erode$')
allocate (work_i2(2,nnpnmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_i2 in erode$')
allocate (work_i3(nnpnmax),stat=err)   ; if (err.ne.0) call stop_run ('Error alloc work_i3 in erode$')
allocate (work_r1(nnpnmax,2),stat=err) ; if (err.ne.0) call stop_run ('Error alloc work_r1 in erode$')
allocate (work_r2(nnpnmax),stat=err)   ; if (err.ne.0) call stop_run ('Error alloc work_r2 in erode$')
allocate (lsf(olsf%nnode),stat=err)    ; if (err.ne.0) call stop_run ('Error alloc lsf in erode$')

nnext=.true.
int_method=3
loc=1

lsf=0.d0

do i=iproc+1,olsf%nnode,nproc
   x=olsf%x(i)
   y=olsf%y(i)
   call nn2D (x,y,points,vertices,neighbour,nnext,int_method,centres, &
              hulltriangles,nh,loc,field,nnpnmax,work_r1,work_r2, &
              work_d1,work_i1,work_i2,work_i3,work_d2,work_d3,work_d4, &
              work_d5,work_d6,work_d7,lt_work,ln_work,out,f,df,ddf,v)
   lsf(i)=f
enddo

olsf%lsf(1:olsf%nnode)=0.d0
call mpi_allreduce (lsf,olsf%lsf,olsf%nnode, &
                    mpi_double_precision,mpi_sum,mpi_comm_world,ierr)

deallocate (points)
deallocate (field)
deallocate (vertices)
deallocate (centres)
deallocate (hulltriangles)
deallocate (vis_tlist)
deallocate (vis_elist)
deallocate (add_tlist)
deallocate (lt_work)
deallocate (ln_work)
deallocate (work_d1)
deallocate (work_d2)
deallocate (work_d3)
deallocate (work_d4)
deallocate (work_d5)
deallocate (work_d6)
deallocate (work_d7)
deallocate (work_r1)
deallocate (work_r2)
deallocate (work_i1)
deallocate (work_i2)
deallocate (work_i3)
deallocate (lsf)

return
end
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|