Skip to content
Snippets Groups Projects
Select Git revision
  • feature/mumps
  • develop
  • feature/parabolic_subduction
  • feature/hdf5out
  • feature/optimisedMoveCloud
  • feature/thrust_isobc
  • feature/particleAdvectionFast
  • feature/orography_sinclair
  • feature/mixed_cloudmat
  • feature/improved_cloud_material_numbering
  • feature/multibox2
  • feature/EGU2020
  • feature/mumps_final
  • feature/oroprecip_tuebingen
  • MfN_Berlin
  • feature/cloudmat-sphere
  • feature/baselevel
  • feature/LK_manuscript_2018
  • feature/sb04_oroprecip
  • feature/oroprecip2
  • v0.5.1
  • v0.5
  • v0.4.3.0
  • v0.4.2.0
  • v0.4.1.0
  • v0.4.0.2
  • v0.4.0.1
  • v0.4
  • v0.4rc1
  • DOUAR-WSMP_v0.1.1
  • DOUAR-WSMP_v0.1
  • 2009-09-02
  • fromDavidWhipp
33 results
Blame Permalink
erosion.f90 16.46 KiB 
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              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

implicit none
      
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,lt_work,ln_work
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 alpha
double precision zerosion
double precision length_scale,velocity_scale,dt
double precision fluvial_erosion,diffusion_erosion
integer baselevelx0,baselevelx1,baselevely0,baselevely1

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

INCLUDE 'mpif.h'

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
   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+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
      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 (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)

! 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)
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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|