Modified to use surface0

src/DOUAR.f90

@@ -295,12 +295,13 @@ do while (istep.le.params%nstep)
          !------------------------------------------------------------------------!
          ref_count=1   
          nadd=1
-         do while (nadd>0)
+         nrem=1
+         do while (nadd+nrem>0)
             call show_time (total,step,inc,1,'refine surface$')
-            call refine_surface(params,surface(is),threadinfo,nadd,is,istep,ref_count) 
+            call refine_surface(params,surface(is),surface0(is),threadinfo,nadd,nrem,is,istep,ref_count) 
             call show_time (total,step,inc,1,'check delaunay$')
             call check_delaunay (params,surface(is),is,istep,'delaunay_after_refine',ref_count)
-            call DoRuRe_surf_stats (params%doDoRuRe,istep,ref_count,is,surface(is)%nt,surface(is)%nsurface,nedge,nadd)
+            call DoRuRe_surf_stats (params%doDoRuRe,istep,ref_count,is,surface(is)%nt,surface(is)%nsurface,nedge,nadd+nrem)
             ref_count=ref_count+1
          end do
       endif

src/refine_surface.f90

@@ -16,7 +16,8 @@
 !------------------------------------------------------------------------------|
 !------------------------------------------------------------------------------|
       
-subroutine refine_surface (params,surface,threadinfo,nadd,is,istep,ref_count) 
+subroutine refine_surface(params,surface,surface0,threadinfo,nadd,nrem,is,&
+           istep,ref_count) 
 
 !------------------------------------------------------------------------------|
 !(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
@@ -167,7 +168,7 @@ distmin=1.d0/(2.d0**surface%levelt+4)*sqrt(2.d0)
 allocate (refine(nedge),stat=err)
 allocate (refinep(nedge),stat=err)
 allocate (remove(nedge),stat=err)
-allocate (removep(ndege),stat=err)
+allocate (removep(nedge),stat=err)
 
 refine=.false.
 refinep=.false.
@@ -335,7 +336,7 @@ do i=1,nrem
    iedge=remove_list(i)
    i1=ed(iedge)%n1
    i2=ed(iedge)%n2
-   call remove_point(i1,surface,1,surface0)
+   call remove_point(i1,surface,0,surface0)
 end do
 
 if (.not.params%normaladvect) then