!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! EROSION Nov. 2006 |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
subroutine erosion (surface,olsf,is,params)
!------------------------------------------------------------------------------|
!(((((((((((((((( 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
! length_scale is the length scale in km (corresponding to the unit length of the box)
! velocity_scale is the scale for the velocity in km/Myr
! dt is the length of the dimensionless time step
! fluvial_erosion is the fluvial erosion constant
! diffusion_erosion is the diffusion erosion constant
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine arguments ))))))))))))))))))))
!------------------------------------------------------------------------------|
use definitions
!use mpi
implicit none
include 'mpif.h'
type (sheet) surface
type (octreelsf) olsf
type (parameters) params
integer is
!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
integer ierr,iproc,nproc,i,nt,nh,nohalt_hull,loc,nnn,nnlist,ntrilist,k,kp,ic,icp
integer ntmax,nhmax,npmax,nmax,nmode,nvmax,nnpnmax,int_method
double precision,dimension(:),allocatable::field,wpoints
double precision,dimension(:,:),allocatable::points,centres,pointsmin,pointsmax,ppoints
integer,dimension(:,:),allocatable::vertices,neighbour
integer,dimension(:),allocatable::hulltriangles,ln_work
logical*1,dimension(:),allocatable::lt_work_l,ln_work_l
integer,dimension(:),allocatable::vis_tlist,vis_elist,add_tlist
logical clockwise,out,nnext
double precision eps
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
double precision,dimension(:),allocatable::xdouar,ydouar,hdouar,hdouarp
integer nnodedouar
double precision dtdouar,href
double precision dhminfluv,dhmaxfluv,dhmindiff,dhmaxdiff
double precision zerosion
double precision length_scale,velocity_scale,dt
double precision fluvial_erosion,diffusion_erosion
integer baselevelx0,baselevelx1,baselevely0,baselevely1
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)
zerosion=params%zerosion
length_scale=params%length_scale
velocity_scale=params%velocity_scale
dt=params%dt
fluvial_erosion=params%fluvial_erosion
diffusion_erosion=params%diffusion_erosion
baselevelx0=params%baselevelx0
baselevelx1=params%baselevelx1
baselevely0=params%baselevely0
baselevely1=params%baselevely1
if (is.ne.1) then
if (surface%spread_surface_points.gt.0) then
allocate (points(2,surface%nsurface),stat=err); if (err.ne.0) call stop_run ('Error alloc points in erode$')
allocate (pointsmin(2,surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc pointsmin in erode$')
allocate (pointsmax(2,surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc pointsmax in erode$')
allocate (ln_work(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ln_work in erode$')
points=0.d0
pointsmin=2.d0
pointsmax=-1.d0
ln_work=0
do i=1,surface%nt
do k=1,3
kp=mod(k,3)+1
ic=surface%icon(k,i)
icp=surface%icon(kp,i)
points(1,ic)=points(1,ic)+surface%x(icp)
points(2,ic)=points(2,ic)+surface%y(icp)
pointsmin(1,ic)=min(pointsmin(1,ic),surface%x(icp))
pointsmax(1,ic)=max(pointsmax(1,ic),surface%x(icp))
pointsmin(2,ic)=min(pointsmin(2,ic),surface%y(icp))
pointsmax(2,ic)=max(pointsmax(2,ic),surface%y(icp))
ln_work(ic)=ln_work(ic)+1
enddo
enddo
do i=1,surface%nsurface
if (surface%x(i).lt.1.d-10) ln_work(i)=0
if (surface%x(i).gt.1.d0-1.d-10) ln_work(i)=0
if (surface%y(i).lt.1.d-10) ln_work(i)=0
if (surface%y(i).gt.1.d0-1.d-10) ln_work(i)=0
enddo
do i=1,surface%nsurface
if (ln_work(i).ne.0) then
if (surface%spread_surface_points.eq.1) then
surface%x(i)=points(1,i)/ln_work(i)
surface%y(i)=points(2,i)/ln_work(i)
elseif (surface%spread_surface_points.eq.2) then
surface%x(i)=(pointsmin(1,i)+pointsmax(1,i))/2.d0
surface%y(i)=(pointsmin(2,i)+pointsmax(2,i))/2.d0
elseif (surface%spread_surface_points.eq.3) then
points(1,i)=points(1,i)/ln_work(i)
points(2,i)=points(2,i)/ln_work(i)
endif
endif
enddo
if (surface%spread_surface_points.eq.3) then
allocate (ppoints(2,surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ppoints in erode$')
allocate (wpoints(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc wpoints in erode$')
wpoints=0.d0
do i=1,surface%nt
do k=1,3
kp=mod(k,3)+1
ic=surface%icon(k,i)
icp=surface%icon(kp,i)
wpoints(ic)=wpoints(ic)+sqrt((surface%x(icp)-points(1,ic))**2+(surface%y(icp)-points(2,ic))**2)*0.1d0
enddo
enddo
ppoints=0.d0
do i=1,surface%nt
do k=1,3
kp=mod(k,3)+1
ic=surface%icon(k,i)
icp=surface%icon(kp,i)
ppoints(1,ic)=ppoints(1,ic)+sqrt((surface%x(icp)-points(1,ic))**2+(surface%y(icp)-points(2,ic))**2)*0.1d0*surface%x(icp)
ppoints(2,ic)=ppoints(2,ic)+sqrt((surface%x(icp)-points(1,ic))**2+(surface%y(icp)-points(2,ic))**2)*0.1d0*surface%y(icp)
enddo
enddo
do i=1,surface%nsurface
if (ln_work(i).ne.0) then
surface%x(i)=ppoints(1,i)/wpoints(i)
surface%y(i)=ppoints(2,i)/wpoints(i)
endif
enddo
deallocate (ppoints,wpoints)
endif
deallocate (points)
deallocate (pointsmin)
deallocate (pointsmax)
deallocate (ln_work)
endif
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)
surface%z(i)=min(f,surface%z(i))
enddo
return
endif
! first apply erosion to top surface
if (length_scale.gt.0.d0) then
nnodedouar=surface%nsurface
allocate (xdouar(nnodedouar))
allocate (ydouar(nnodedouar))
allocate (hdouar(nnodedouar))
allocate (hdouarp(nnodedouar))
xdouar=surface%x*length_scale*1.d3
ydouar=surface%y*length_scale*1.d3
if (iproc.eq.0) call random_number (hdouar)
call mpi_bcast (hdouar,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
hdouar=(hdouar+1.0)/2.d0
hdouar=hdouar+(surface%z-zerosion)*length_scale*1.d3*params%vex
href=minval(hdouar)
hdouar=hdouar-href
dtdouar=dt*length_scale/velocity_scale*1.d6
if (iproc.eq.0) write (8,*) 'Cascade time scale ',dtdouar
if (iproc.eq.0) write (8,*) 'Topo in ',minval(hdouar),maxval(hdouar)
if (iproc.eq.0) write (8,*) 'X in ',minval(xdouar),maxval(xdouar)
if (iproc.eq.0) write (8,*) 'Y in ',minval(ydouar),maxval(ydouar)
hdouarp=hdouar
call cascade (xdouar,ydouar,hdouar,nnodedouar,dtdouar, &
dhminfluv,dhmaxfluv, &
dhmindiff,dhmaxdiff, &
fluvial_erosion,diffusion_erosion, &
baselevelx0,baselevelx1,baselevely0,baselevely1)
if (iproc.eq.0) write (8,*) 'Topo out ',minval(hdouar),maxval(hdouar)
if (iproc.eq.0) write (8,*) 'Delta Topo ',minval(hdouar-hdouarp),maxval(hdouar-hdouarp)
do i=1,surface%nsurface
surface%z(i)=(hdouar(i)+href)/1.d3/length_scale/params%vex+zerosion
surface%r(i)=(hdouar(i)-hdouarp(i))/1.d3/length_scale/dt
enddo
deallocate (xdouar)
deallocate (ydouar)
deallocate (hdouar)
deallocate (hdouarp)
else ! perfect erosion
do i=1,surface%nsurface
if (surface%z(i).gt.zerosion) surface%z(i)=zerosion
enddo
endif
! 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 (ln_work(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ln_work in erode$')
allocate (lt_work_l(ntmax),stat=err) ; if (err.ne.0) call stop_run ('Error alloc lt_work_l in erode$')
allocate (ln_work_l(surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ln_work_l 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_l,ln_work_l)
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_l,ln_work_l,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)
! here spreads the points on the top surface to avoid collapse of the surface due to erosion
! for that it uses a barrycenter algorithm and a reinterpolation of the surface to avoid undesired diffusion effects
if (surface%spread_surface_points.gt.0) then
allocate (pointsmin(2,surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc pointsmin in erode$')
allocate (pointsmax(2,surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc pointsmax in erode$')
points=0.d0
pointsmin=2.d0
pointsmax=-1.d0
ln_work=0
do i=1,surface%nt
do k=1,3
kp=mod(k,3)+1
ic=surface%icon(k,i)
icp=surface%icon(kp,i)
points(1,ic)=points(1,ic)+surface%x(icp)
points(2,ic)=points(2,ic)+surface%y(icp)
pointsmin(1,ic)=min(pointsmin(1,ic),surface%x(icp))
pointsmax(1,ic)=max(pointsmax(1,ic),surface%x(icp))
pointsmin(2,ic)=min(pointsmin(2,ic),surface%y(icp))
pointsmax(2,ic)=max(pointsmax(2,ic),surface%y(icp))
ln_work(ic)=ln_work(ic)+1
enddo
enddo
do i=1,surface%nsurface
if (surface%x(i).lt.1.d-10) ln_work(i)=0
if (surface%x(i).gt.1.d0-1.d-10) ln_work(i)=0
if (surface%y(i).lt.1.d-10) ln_work(i)=0
if (surface%y(i).gt.1.d0-1.d-10) ln_work(i)=0
enddo
do i=1,surface%nsurface
if (ln_work(i).ne.0) then
if (surface%spread_surface_points.eq.1) then
surface%x(i)=points(1,i)/ln_work(i)
surface%y(i)=points(2,i)/ln_work(i)
else
surface%x(i)=(pointsmin(1,i)+pointsmax(1,i))/2.d0
surface%y(i)=(pointsmin(2,i)+pointsmax(2,i))/2.d0
endif
endif
enddo
deallocate (pointsmin,pointsmax)
endif
deallocate (points)
deallocate (field)
deallocate (vertices)
deallocate (centres)
deallocate (hulltriangles)
deallocate (vis_tlist)
deallocate (vis_elist)
deallocate (add_tlist)
deallocate (lt_work_l)
deallocate (ln_work)
deallocate (ln_work_l)
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|