DOUAR.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

program DOUAR

use threads
use definitions
use version
!use mpi

implicit none

include 'mpif.h'

integer i,j,is,iter,istep,niter,nedge,nedgen,ileaves,matel
integer ntriangle,ndof,ndoft,material0,nsurfacen
integer numarg,iargc
integer ierr,nproc,iproc,k,nrem
integer current_level
integer ie_ov,ie_loc,ie_level,iter_nl
integer err,iunit,cltemp
integer nremove,ninject,nnmax,nadd,naddp,ref_count
integer nleaves_new_mem1
integer nleaves_new_mem2
integer nleaves_old_mem1
integer nleaves_old_mem2
integer ematnump,surfcount
integer size_str(2)
integer, external :: ioctree_number_of_elements
integer,dimension(:,:),allocatable :: cloud_elem_mat_bins
integer,dimension(:),allocatable :: leaf_mat_bin

double precision x,y,z,z0,u,v,w
double precision exx,eyy,ezz,exy,eyz,ezx,e2dmax
double precision duvw,uvw,dtcourant,current_time
double precision umax,xxx,yyy,zzz,x0_leaf,y0_leaf,z0_leaf
double precision total,step,inc
double precision height_max
real  time1,time2,time3 
      
character     :: ic*7
character*3   :: ciproc
character*4   :: cs,cs2
character*40  :: string
character*72  :: shift

logical converge,increase_current_level,velocity_converged,usecourant
logical all_surf_fixed_spinup,all_surf_fixed
logical surf_fixed_spinup,surf_fixed

type (sheet),dimension(:),allocatable::surface,surface0
type (octreev) ov
type (octreelsf) olsf
type (octreesolve) osolve
type (material),dimension(:),allocatable::mat
type (void) vo
type (cloud) cl
type (topology),dimension(:),allocatable::tpl
type (box),dimension(:),allocatable::boxes
type (cross_section),dimension(:),allocatable::sections
type (face),dimension(6) :: cube_faces
type (edge),dimension(:),allocatable::ed,edswap
type (parameters) params
type (thread) threadinfo
type (bc_definition) bcdef
type (nest_info) nest
type (string),dimension(:,:),allocatable::density_str

integer n_iproc_st,n_iproc_end,n_iproc
integer ldb,nrhs,n,nz_loc
integer, dimension(:), allocatable :: ia,ja
logical, dimension(:), allocatable :: iproc_col
double precision, dimension(:), allocatable :: avals
double precision, dimension(:,:), allocatable :: b
double precision, dimension(:), allocatable :: weightel
double precision, dimension(:), allocatable :: xyz_t

shift=' '

! DOUAR version number
! vermajor is incremented for major rewrites
! verminor is incremented for significant revisions or changes to the input or
!   output files
! verstat is used to designated the current code status: 0=alpha, 1=beta, 2=rc
!   3=release
! verrev is incremented for minor changes and bugfixes
params%vermajor=0
params%verminor=2
params%verstat=0
params%verrev=1

ndof=3
ndoft=1
params%mpe=8
nrhs = 1

nleaves_old_mem1=0
nleaves_old_mem2=0

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

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

!call write_streepjes(6,2)

!!=====[allocate threadinfo and open mpi log and memory heap files]=========
!params%nmpi=nproc
!
!call int_to_char(ciproc,3,iproc)
!
!threadinfo%Logunit=1000+iproc
!open (unit=threadinfo%Logunit,file='./DEBUG/mpilogs/Log_'//ciproc//'.dat',status='replace')
!write(threadinfo%Logunit,*) 'This is DOUAR-WSMP v0.1 (2011-03-23)'
!if (iproc.eq.0) write(*,*)  'This is DOUAR-WSMP v0.1 (2011-03-23), using ', nproc, ' processors.'
!write(threadinfo%Logunit,'(a,i3)') 'Log file of mpi process',iproc
!
!call heap_init(threadinfo,2000+iproc,'./DEBUG/mpilogs/mem_heap_'//ciproc//'.dat')

params%nmpi=nproc
!if (iproc.eq.0) write(*,'(a,i3,a)')  'This is DOUAR-WSMP v0.1 (2011-03-23), using ', nproc, ' processors.'

if (iproc.eq.0) then
  write (*,'(a)') '!----------------------------------------------------------------------|'
  write (*,'(a)') '!----------------------------------------------------------------------|'
  write (*,'(a)') '!                                                                      |'
  write (*,'(a)') '!      8888888b.   .d88888b.  888     888       d8888 8888888b.        |'
  write (*,'(a)') '!      888  "Y88b d88P" "Y88b 888     888      d88888 888   Y88b       |'
  write (*,'(a)') '!      888    888 888     888 888     888     d88P888 888    888       |'
  write (*,'(a)') '!      888    888 888     888 888     888    d88P 888 888   d88P       |'
  write (*,'(a)') '!      888    888 888     888 888     888   d88P  888 8888888P"        |'
  write (*,'(a)') '!      888    888 888     888 888     888  d88P   888 888 T88b         |'
  write (*,'(a)') '!      888  .d88P Y88b. .d88P Y88b. .d88P d8888888888 888  T88b        |'
  write (*,'(a)') '!      8888888P"   "Y88888P"   "Y88888P" d88P     888 888   T88b       |'
  write (*,'(a)') '!                                                                      |'
  write(*,'(a,i1,a,i1,a,i1,a,i1,a,a8,a)') '!     DOUAR-WSMP version ',params%vermajor,'.',params%verminor,'.',params%verstat,'.',params%verrev,' - r', svnrev,'                           |'
  write(*,'(a,i3,a)') '!     Running on ', nproc, ' processor cores                                   |'
  write (*,'(a)') '!----------------------------------------------------------------------|'
  write (*,'(a)') '!----------------------------------------------------------------------|'
endif

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     is there any input file to douar ?
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

numarg=iargc()
! 2009-09-03: Douglas Guptill
! mpich doesn't allow us to read the command line
! so we simulate an empty one
numarg = 0

if (numarg==0) then
   params%infile='input.txt'
   if (iproc.eq.0) write(*,'(a)') 'Program called with no argument'
else
   call getarg (1,params%infile)
   if (iproc.eq.0) then
      write(*,'(a,a)') 'program called with input file ',params%infile
   end if
end if

!call write_streepjes(6,2)
call write_streepjes(6,1)
call show_time (total,step,inc,0,'Start of Computations$')
call show_time (total,step,inc,1,'Reading Input$')

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     read debug level
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call read_controlling_parameters (params,threadinfo,'debug')

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     allocate threadinfo, open MPI log and memory heap files
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call int_to_char(ciproc,3,iproc)
if (params%debug.gt.1) then
  threadinfo%Logunit=1000+iproc
  open (unit=threadinfo%Logunit,file='./DEBUG/mpilogs/Log_'//ciproc//'.dat',status='replace')
  write(threadinfo%Logunit,*) 'This is DOUAR-WSMP v0.1 (2011-03-23)'
  write(threadinfo%Logunit,'(a,i3)') 'Log file of mpi process',iproc
  write(threadinfo%Logunit,'(a16,i3)') 'debug',params%debug
  call heap_init(threadinfo,2000+iproc,'./DEBUG/mpilogs/mem_heap_'//ciproc//'.dat')
endif

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     read controlling parameters
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call read_controlling_parameters (params,threadinfo,'main')

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     open and prepare output files
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call io_DoRuRe2 (params,'init')

allocate (surface         (  params%ns  ),stat=err); if (err.ne.0) call stop_run ('Error alloc surface in main$')
allocate (surface0        (  params%ns  ),stat=err); if (err.ne.0) call stop_run ('Error alloc surface0 in main$')
allocate (mat             (0:params%nmat),stat=err); if (err.ne.0) call stop_run ('Error alloc mat in main$')
allocate (params%materialn(0:params%ns  ),stat=err); if (err.ne.0) call stop_run ('Error alloc materialn in main$')
allocate (params%surface_for_mat_props(params%ns),stat=err); if (err.ne.0) call stop_run ('Error alloc surface_for_mat_props in main$')
      
if (params%nboxes.gt.0) then
    allocate (boxes(params%nboxes),stat=err) ; if (err.ne.0) call stop_run ('Error alloc boxes arrays$')
else
    allocate (boxes(1),stat=err) ! necessary to avoid nil size argument
endif
if (params%nsections.gt.0) then
    allocate (sections(params%nsections),stat=err) ; if (err.ne.0) call stop_run ('Error alloc cross_section arrays$')
else
    allocate (sections(1),stat=err) ! necessary to avoid nil size argument
end if











!Allocate density_str, %density, %densityp if using compaction
if (params%compaction) then
  allocate (density_str(2**params%levelmax_oct,2**params%levelmax_oct),stat=err); if (err.ne.0) call stop_run ('Error alloc density_str in main$')
  do j=1,2**params%levelmax_oct
    do i=1,2**params%levelmax_oct
      allocate (density_str(i,j)%density(1),stat=threadinfo%err)       
      allocate (density_str(i,j)%densityp(1),stat=threadinfo%err) 
    enddo 
  enddo
endif








!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     read input file for all parameters of the run
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call read_input_file (params,threadinfo,material0,mat,surface,boxes,sections,  &
                     cube_faces,nest,bcdef)

if (params%irestart==0) then
   current_level=params%initial_refine_level
else
   current_level=params%levelmax_oct
end if

call show_time (total,step,inc,1,'Problem Setup$')

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     either reads or creates the surface(s)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call show_time (total,step,inc,1,'define surfaces$')
call define_surface(params,surface,threadinfo,total,step,inc,current_time)

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     either reads or creates the cloud  
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call show_time (total,step,inc,1,'define cloud$')
call define_cloud(cl,params,bcdef)

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     extract material information to be passed to solver (moved down below)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!params%materialn(0)=material0
!do i=1,params%ns
!   params%materialn(i)=surface(i)%material
!enddo

! Pass information about which surfaces should be used for material property
! assignment when using the cloud to define material properties
if (params%materials_on_cloud) then
  do i=1,params%ns
    params%surface_for_mat_props(i)=surface(i)%surf_for_mat_props
  enddo
endif

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!    compute plastic parameters (obsolete: calculations done in vrm.f90)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!call compute_plastic_params (params,mat)

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     either reads or creates the velocity octree
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
call show_time (total,step,inc,1,'define velocity octree$')
call define_ov (ov,params,threadinfo)

istep=1+params%irestart

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     is the CFL condition used for the timestep ?
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

!usecourant = (params%dt .lt. 0.d0)

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     Will all or some surfaces have fixed positions?
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
all_surf_fixed_spinup=.true.
all_surf_fixed=.true.
surf_fixed_spinup=.false.
surf_fixed=.false.
do is=1,params%ns
  surface(is)%fixed_surf_spinup=(surface(is)%fixed_surf_spinup .or. params%all_surf_fixed_spinup)
  all_surf_fixed_spinup=(all_surf_fixed_spinup .and. surface(is)%fixed_surf_spinup)
  all_surf_fixed=(all_surf_fixed .and. surface(is)%fixed_surf)
  surf_fixed_spinup=(surf_fixed_spinup .or. surface(is)%fixed_surf_spinup)
  surf_fixed=(surf_fixed .or. surface(is)%fixed_surf)
enddo


!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     beginning of time stepping
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
increase_current_level=.false.
do while (istep.le.params%nstep) 

   if (istep==1) params%first_step = .true.

   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   !     Are we in the main calculation phase or the spinup phase?
   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   if (istep.le.params%nstep_spinup) then
     params%in_spinup=.true.
   else
     params%in_spinup=.false.
   endif
   if (bcdef%bctype == 'iso_only') params%in_spinup=.true.
   if (params%in_spinup) then
     params%dt=params%dt_spinup
     params%griditer=params%griditer_spinup
     params%nonlinear_iterations=params%nonlinear_iterations_spinup
     params%tol=params%tol_spinup
     params%materialn(0)=material0
     do i=1,params%ns
       surface(i)%material=surface(i)%material_spinup
     enddo        
   else     
     params%dt=params%dt_main
     params%griditer=params%griditer_main
     params%nonlinear_iterations=params%nonlinear_iterations_main
     params%tol=params%tol_main
     params%materialn(0)=material0
     do i=1,params%ns
       surface(i)%material=surface(i)%material_main
     enddo        
   endif

   ! Update material order for LSF assignment
   do i=1,params%ns
     params%materialn(i)=surface(i)%material
   enddo        

   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   !     is the CFL condition used for the timestep ?
   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   usecourant = (params%dt .lt. 0.d0)

   call write_streepjes(6,2)
   call write_streepjes(8,2)
   call show_time (total,step,inc,-istep,'Start of Step $')
   if (iproc.eq.0) write(8,*)  'Doing step ',istep
   call write_streepjes(6,2)
   call write_streepjes(8,2)
   if (params%debug.gt.1) call heap_hop1(threadinfo,istep)

   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   !     refining surfaces
   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   
   call show_time (total,step,inc,1,'Loop over surfaces$')
   
   do is=1,params%ns
      call int_to_char(ic,2,is)
      call show_time (total,step,inc,1,'surface '//ic(1:2)//'$')

      if (iproc.eq.0 ) then
         write (8,*) 'Surface ',is,' has ',surface(is)%nsurface,'particles before refinement'
         if (params%debug>=1) write(*,'(a,i2,a,i8,a)') shift//'S.',is,': ',surface(is)%nsurface,'points'
      end if

      if (surface(is)%activation_time.ge.current_time+tiny(current_time)) then
         !------------------------------------------------------------------------!
         ! if surface is slave to top surface, no refinement
         !------------------------------------------------------------------------!
         deallocate(surface(is)%x,surface(is)%y,surface(is)%z)
         deallocate(surface(is)%xn,surface(is)%yn,surface(is)%zn)
         deallocate(surface(is)%r,surface(is)%s)
         deallocate(surface(is)%icon)
         surface(is)%nsurface=surface(1)%nsurface
         surface(is)%nt=surface(1)%nt
         allocate (surface(is)%x(surface(is)%nsurface),surface(is)%y(surface(is)%nsurface),surface(is)%z(surface(is)%nsurface))
         allocate (surface(is)%xn(surface(is)%nsurface),surface(is)%yn(surface(is)%nsurface),surface(is)%zn(surface(is)%nsurface))
         allocate (surface(is)%r(surface(is)%nsurface),surface(is)%s(surface(is)%nsurface))
         allocate (surface(is)%icon(3,surface(is)%nt))
         surface(is)%x=surface(1)%x
         surface(is)%y=surface(1)%y
         surface(is)%z=surface(1)%z
         surface(is)%xn=surface(1)%xn
         surface(is)%yn=surface(1)%yn
         surface(is)%zn=surface(1)%zn
         surface(is)%r=surface(1)%r
         surface(is)%s=surface(1)%s
         surface(is)%icon=surface(1)%icon(:,1:surface(is)%nt)
      else
         !------------------------------------------------------------------------!
         ! if surface is no slave to top surface, do refine
         !------------------------------------------------------------------------!
         ref_count=1   
         nadd=1
         nrem=1
         do while (nadd+nrem>0)
            call show_time (total,step,inc,1,'refine surface$')
            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+nrem)
            ref_count=ref_count+1
         end do
      endif

      if (iproc.eq.0 ) then
         write (8,*) 'Surface ',is,' has ',surface(is)%nsurface,'particles after refinement'
         if (params%debug>=1) write(*,'(a,i2,a,i8,a)') shift//'S.',is,': ',surface(is)%nsurface,'points'
      end if

   end do

   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   !     Stores the geometry of the surfaces for the mid-point iterative scheme
   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   do is=1,params%ns
      surface0(is)%nsurface=surface(is)%nsurface
      allocate (surface0(is)%r(surface0(is)%nsurface),stat=threadinfo%err)       
      if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%r','main',size(surface0(is)%r),'dp',+1)
      allocate (surface0(is)%s(surface0(is)%nsurface),stat=threadinfo%err)       
      if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%s','main',size(surface0(is)%s),'dp',+1)
      allocate (surface0(is)%x(surface0(is)%nsurface),stat=threadinfo%err)       
      if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%x','main',size(surface0(is)%x),'dp',+1)
      allocate (surface0(is)%y(surface0(is)%nsurface),stat=threadinfo%err)       
      if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%y','main',size(surface0(is)%y),'dp',+1)
      allocate (surface0(is)%z(surface0(is)%nsurface),stat=threadinfo%err)       
      if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%z','main',size(surface0(is)%z),'dp',+1)
      allocate (surface0(is)%xn(surface0(is)%nsurface),stat=threadinfo%err)      
      if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%xn','main',size(surface0(is)%xn),'dp',+1)
      allocate (surface0(is)%yn(surface0(is)%nsurface),stat=threadinfo%err)      
      if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%yn','main',size(surface0(is)%yn),'dp',+1)
      allocate (surface0(is)%zn(surface0(is)%nsurface),stat=threadinfo%err)      
      if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%zn','main',size(surface0(is)%zn),'dp',+1)
      surface0(is)%r=surface(is)%r
      surface0(is)%s=surface(is)%s
      surface0(is)%x=surface(is)%x
      surface0(is)%y=surface(is)%y
      surface0(is)%z=surface(is)%z
      surface0(is)%xn=surface(is)%xn
      surface0(is)%yn=surface(is)%yn
      surface0(is)%zn=surface(is)%zn
   enddo




   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   !     Store the cloud particle positiions for the mid-point iterative scheme|
   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   do i=1,cl%np
     cl%x0mp(i)=cl%x(i)
     cl%y0mp(i)=cl%y(i)
     cl%z0mp(i)=cl%z(i)
   enddo


   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   !     beginning of grid iterations
   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|

   iter=1

   do while (iter.le.abs(params%griditer).or.iter.eq.1)

      if (iter==1 .and. params%first_step) params%first_iter = .true.

      call write_streepjes(6,1)
      call write_streepjes(8,1)
      call show_time (total,step,inc,-iter,' Start of iteration $')
      if(iproc.eq.0) write(8,*)  'Doing iteration ',iter
      call write_streepjes(6,1)
      call write_streepjes(8,1)
      if (params%debug.gt.1) call heap_hop2(threadinfo,iter)
      call show_time (total,step,inc,1,'Create octree solve$')

      converge=.false.

!      if (params%griditer .lt. 0 .and. current_level==params%levelmax_oct) converge=.true.

      if (iproc.eq.0 .and. params%debug>=1) then
         write(*,'(a,L1)') shift//'(1) params%griditer < 0               ',(params%griditer.lt.0) 
         write(*,'(a,L1)') shift//'(2) current_level==params%levelmax_oct',(current_level==params%levelmax_oct) 
         write(*,'(a,L1)') shift//'(3) increase_current_level            ',increase_current_level
         write(*,'(a,L1)') shift//'converge = (1) & (2) & (3) -> ',(params%griditer .lt. 0 .and. current_level==params%levelmax_oct .and. increase_current_level) 
      end if

      if (params%griditer .lt. 0 .and. current_level==params%levelmax_oct .and. increase_current_level) converge=.true.

      if (iproc.eq.0 .and. converge .and. params%debug>=1) then
         write(*,'(a)') shift//'+++++++++++++++++'
         write(*,'(a)') shift//'grid conv reached'
         write(*,'(a)') shift//'+++++++++++++++++'
      end if

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     creation of osolve
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      ! we create the large "solve" octree, ie the octree on which the finite 
      ! elements will be built.

      osolve%noctree=params%noctreemax

      allocate (osolve%octree(osolve%noctree),stat=threadinfo%err)      
      if (params%debug.gt.1) call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',+1)

      call octree_init (osolve%octree,osolve%noctree)
      call octree_create_uniform (osolve%octree,osolve%noctree,params%leveluniform_oct)

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     refinement/improvement of osolve
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call improve_osolve (osolve,ov,params,threadinfo,istep,iter,total,step,  &
                          inc,current_level,boxes,cube_faces)
      !call improve_osolve (osolve,ov,params,threadinfo,istep,iter,total,step,inc, &
      !                      current_level,boxes,cube_faces,increase_current_level )

      nleaves_new_mem1=osolve%octree(2)

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     find connectivity for osolve 
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'Build osolve icon$')
      
      osolve%nleaves = ioctree_number_of_elements (osolve%octree,osolve%noctree)
      osolve%nnode   = osolve%nleaves*3
      osolve%nlsf    = params%ns

      if (iproc.eq.0) write (8,*) osolve%nleaves,' leaves in solve octree '

      allocate (osolve%icon(8,osolve%nleaves),stat=threadinfo%err)        
      if (params%debug.gt.1) call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',+1)
      allocate (osolve%x(osolve%nnode),stat=threadinfo%err)               
      if (params%debug.gt.1) call heap (threadinfo,'osolve%x',   'main',size(osolve%x),'dp',+1)
      allocate (osolve%y(osolve%nnode),stat=threadinfo%err)               
      if (params%debug.gt.1) call heap (threadinfo,'osolve%y',   'main',size(osolve%y),'dp',+1)
      allocate (osolve%z(osolve%nnode),stat=threadinfo%err)               
      if (params%debug.gt.1) call heap (threadinfo,'osolve%z',   'main',size(osolve%z),'dp',+1)

      call octree_find_node_connectivity (osolve%octree,osolve%noctree, &
                                          osolve%icon,osolve%nleaves,osolve%x,osolve%y,osolve%z,osolve%nnode)

      ! osolve%nnode has been changed by octree_find_node_connectivity, so  re-size x,y,z
      call octreesolve_shrink_xyz(osolve, threadinfo, params)

      ! now that osolve%nnode is known we can allocate osolve%lsf
      allocate (osolve%lsf(osolve%nnode,osolve%nlsf),stat=threadinfo%err) 
      if (params%debug.gt.1) call heap (threadinfo,'osolve%lsf', 'main',size(osolve%lsf),'dp',+1)
      osolve%lsf=0.d0

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     embed the surfaces in osolve
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'Embed surface in osolve$')
      do is=1,params%ns
         write (ic(1:2),'(i2)') is
         call show_time (total,step,inc,1,'embedding surface '//ic(1:2)//'$')
         call embed_surface_in_octree (osolve,params,surface(is),is,istep,iter,threadinfo)
      enddo

      nleaves_new_mem2=osolve%octree(2)
 
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     assess grid convergence 
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'Assessing octree stability$')

      increase_current_level=  (nleaves_old_mem2 .ge. nleaves_new_mem2 .and.  & 
                               (dble(nleaves_old_mem2-nleaves_new_mem2)) / (dble(nleaves_old_mem2+nleaves_new_mem2))*2.d0 &
                               .lt. params%octree_refine_ratio )

!      increase_current_level=  ((dble(nleaves_old_mem2-nleaves_new_mem2)) / (dble(nleaves_old_mem2+nleaves_new_mem2))*2.d0 &
!                               .lt. params%octree_refine_ratio )
 
      if(iproc.eq.0 .and. params%debug>=1) then 
        write(*,'(a,i5,a)') shift//'current refine level ',current_level
        write(*,'(a)')      shift//'After criterion based refinement:'
        write(*,'(a,i8,a)') shift//'nleaves_old= ',nleaves_old_mem1,' leaves'
        write(*,'(a,i8,a)') shift//'nleaves_new= ',nleaves_new_mem1,' leaves'
        write(*,'(a)')      shift//'After surfaces embedding:'
        write(*,'(a,i8,a)') shift//'nleaves_old= ',nleaves_old_mem2,' leaves'
        write(*,'(a,i8,a)') shift//'nleaves_new= ',nleaves_new_mem2,' leaves'
        write(*,'(a,l1)')   shift//'C2: authorise increase of refine level = ',increase_current_level
      end if

      nleaves_old_mem1=nleaves_new_mem1
      nleaves_old_mem2=nleaves_new_mem2

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     allocate fields of osolve 
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      allocate (osolve%u(osolve%nnode),stat=threadinfo%err)             
      if (params%debug.gt.1) call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',+1)
      allocate (osolve%v(osolve%nnode),stat=threadinfo%err)             
      if (params%debug.gt.1) call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',+1)
      allocate (osolve%w(osolve%nnode),stat=threadinfo%err)             
      if (params%debug.gt.1) call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',+1)
      allocate (osolve%wiso(osolve%nnode),stat=threadinfo%err)             
      if (params%debug.gt.1) call heap (threadinfo,'osolve%wiso','main',size(osolve%wiso),'dp',+1)
      allocate (osolve%temp(osolve%nnode),stat=threadinfo%err)          
      if (params%debug.gt.1) call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',+1)
      allocate (osolve%pressure(osolve%nleaves),stat=threadinfo%err)    
      if (params%debug.gt.1) call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',+1)
      allocate (osolve%spressure(osolve%nleaves),stat=threadinfo%err)   
      if (params%debug.gt.1) call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',+1)
      allocate (osolve%strain(osolve%nnode),stat=threadinfo%err)        
      if (params%debug.gt.1) call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',+1)
      allocate (osolve%e2dp(osolve%nnode),stat=threadinfo%err)        
      if (params%debug.gt.1) call heap (threadinfo,'osolve%e2dp','main',size(osolve%e2dp),'dp',+1)



      allocate (osolve%matnum(osolve%nleaves),stat=threadinfo%err)                               
      if (params%debug.gt.1) call heap (threadinfo,'osolve%matnum','main',size(osolve%matnum),'int',+1)
      allocate (osolve%matnump(osolve%nleaves),stat=threadinfo%err)                               
      if (params%debug.gt.1) call heap (threadinfo,'osolve%matnump','main',size(osolve%matnump),'int',+1)




      if (params%compaction) then
        allocate (osolve%compaction_density(osolve%nleaves),stat=threadinfo%err)   
        if (params%debug.gt.1) call heap (threadinfo,'osolve%compaction_density','main',size(osolve%compaction_density),'dp',+1)
        allocate (osolve%wcompact(osolve%nnode),stat=threadinfo%err)             
        if (params%debug.gt.1) call heap (threadinfo,'osolve%wcompact','main',size(osolve%wcompact),'dp',+1)
        osolve%compaction_density = 0.d0
        osolve%wcompact = 0.d0
      endif



      osolve%u               = 0.d0
      osolve%v               = 0.d0
      osolve%w               = 0.d0
      osolve%wiso            = 0.d0
      osolve%temp            = 0.d0
      osolve%pressure        = 0.d0
      osolve%spressure       = 0.d0
      osolve%strain          = 0.d0
      osolve%e2dp            = 0.d0




      osolve%matnum          = 0
      osolve%matnump         = 0



      ! At this point, should we 'reconstruct' the elemental material numbers
      ! based on the update cloud particle positions?
      !
      ! For the case of materials_on_cloud and first_iter, this should not be
      ! done as the cloud particles do not yet have their material numbers
      
      ! If using the cloud and not in the first iteration, recalculate the
      ! elemental material numbers from the cloud particles
      !
      ! CAN WE PARALLELIZE THIS LOOP TO BE A BIT QUICKER???
      !
      ! PUT THIS IN A SEPARATE SUBROUTINE???
      !
      ! WHAT DO WE DO WITH NEW ELEMENTS? WHAT IS THE MATNUM FOR PARTICLES OUTSIDE THE
      ! MODEL? WHAT HAPPENS IF FIND_LEAF CAN'T FIND A LEAF?
      if (params%materials_on_cloud .and. .not.params%first_iter) then
        call show_time (total,step,inc,1,'Assign osolve mat nums from cloud$')
        allocate (cloud_elem_mat_bins(osolve%nleaves,0:params%nmat),stat=err)
        if (err.ne.0) call stop_run ('Error alloc cloud_elem_mat_bins in build_system$')
        ! Calculate number of cloud particles of each material in all elements
        call find_mat_numbers_from_cloud(params,cl,osolve,cloud_elem_mat_bins)
        ! Allocate small array to store cloud particle material types for individual
        ! elements
        allocate (leaf_mat_bin(0:params%nmat),stat=err)
        if (err.ne.0) call stop_run ('Error alloc leaf_mat_bin in build_system$')
        do ileaves=1,osolve%nleaves
          ! Grab material bins for current element
          leaf_mat_bin=cloud_elem_mat_bins(ileaves,:)
          ematnump=osolve%matnump(ileaves)
          call get_elem_mat_number_from_cloud(params,ileaves,leaf_mat_bin,ematnump,matel)
          osolve%matnum(ileaves) = matel
        enddo
        ! Deallocate material bin arrays
        deallocate(cloud_elem_mat_bins,leaf_mat_bin)
      endif


      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     improve and update cloud fields 
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'Update 3D cloud structure$')
      call update_cloud_structure (cl,osolve,params,mat,ninject,nremove,istep)

!      if (istep==1) call strain_history (params,osolve,cl)
!      call strain_history (params,osolve,cl)

      call DoRuRe_cloud_stats (params%doDoRuRe,istep,iter,cl%np,nremove,ninject,cl%strain,cl%temp,cl%press)

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      ! interpolate velocity from ov to solve, prepare ov receive the new solution
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'Interpolate velo onto osolve$')

      call interpolate_ov_on_osolve (osolve,ov,params,threadinfo)

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     find bad faces of solve octree
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'Find bad faces$')

      osolve%nface=osolve%nleaves
      allocate (osolve%iface(9,osolve%nface),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',+1)

      call octree_find_bad_faces (osolve%octree,osolve%noctree, &
                                  osolve%iface,osolve%nface, &
                                  osolve%icon,osolve%nleaves)

      if (osolve%nface.gt.osolve%nleaves) call stop_run ('nface larger than nleaves$')
      if (iproc.eq.0) write (8,*) osolve%nface,' bad faces in solve octree'
      if (iproc.eq.0) write (8,'(a,i9)') shift//'nb of bad faces',osolve%nface

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      ! find void (which nodes are in the void) and which elements are entirely in 
      ! the void and should not be included in the fe calculations
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'Find void$')

      allocate (vo%node(osolve%nnode),stat=threadinfo%err)              
      if (params%debug.gt.1) call heap (threadinfo,'vo%node','main',size(vo%node),'int',+1)
      allocate (vo%leaf(osolve%nnode),stat=threadinfo%err)              
      if (params%debug.gt.1) call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',+1)
      allocate (vo%face(osolve%nface),stat=threadinfo%err)              
      if (params%debug.gt.1) call heap (threadinfo,'vo%face','main',size(vo%face),'int',+1)
      allocate (vo%rtf(osolve%nnode),stat=threadinfo%err)               
      if (params%debug.gt.1) call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',+1)
      allocate (vo%ftr(osolve%nnode),stat=threadinfo%err)               
      if (params%debug.gt.1) call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',+1)
      allocate (vo%influid(osolve%nnode),stat=threadinfo%err)           
      if (params%debug.gt.1) call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',+1)

      call find_void_nodes (params,osolve,vo,ov,threadinfo)
      !call find_void_nodes (params,osolve,vo,istep,ov,threadinfo)



      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     Assign material numbers to cloud (if using cloud for material types)
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      if (params%materials_on_cloud .and. params%first_iter) then
        call assign_cloud_mat_number(params,osolve,cl,mat)
      endif






      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !                  Resize density_str instances                          |
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|


	  if (params%compaction) then
	    height_max=maxval(surface(1)%z)
        call resize_density_profiles(params,density_str,height_max)
      endif



      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     Assign elemental compaction densities                              |
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      if (params%compaction) then
        if (params%first_iter) call initialize_compaction(params,density_str,osolve,mat)
        call put_str_density_on_osolve(params,density_str,osolve)
      endif










      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     Remove surfaces flagged for removal (if desired)
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      if (params%materials_on_cloud .and. params%first_iter) then
        ! Modify params%ns and remove surfaces
        call show_time (total,step,inc,1,'Remove flagged surfaces$')
        !This is not going to work on restart as coded
        params%nsorig=params%ns
        ! Create array to store info about surfaces that are removed or kept
        allocate(params%surfremindex(params%nsorig))
        surfcount=0
        do is=1,params%nsorig
          if (surface(is)%remove_after_mat_def) then
            params%surfremindex(is)=-1
            if (iproc.eq.0 .and. params%debug.ge.1) write (*,'(a,i2,a)') shift//'Surface ',is,' will be removed'
          else
            surfcount=surfcount+1
            params%surfremindex(is)=surfcount
          endif
        enddo
        params%nsmod = surfcount
        params%ns = params%nsmod
        if (params%nsorig /= params%ns) then
          osolve%nlsf=params%ns

          ! Remove flagged surfaces
          call remove_surfaces(params,surface,0,threadinfo)
          call remove_surfaces(params,surface0,1,threadinfo)

          ! Remove unneeded LSFs
          call remove_lsfs(params,osolve,threadinfo)

          ! Adjust materialn for reduced number of surfaces
          call adjust_materialn(params,threadinfo)
        else
          if (iproc.eq.0 .and. params%debug.ge.1) write (*,'(a)') shift//'No surfaces/LSFs to remove'
        endif
      endif



      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     definition of the bc
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'Define boundary conditions$')
      allocate (osolve%kfix(osolve%nnode*3),stat=threadinfo%err)        
      if (params%debug.gt.1) call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',+1)
      allocate (osolve%kfixt(osolve%nnode),stat=threadinfo%err)         
      if (params%debug.gt.1) call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',+1)

      call define_bc (params,osolve,vo,bcdef,nest)

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      ! build matrix arrays, allocate memory for wsmp
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'wsmp setup1$')

      n = vo%nnode * ndof
     
      !---[topology]----- 
      allocate (tpl(n)) 
      tpl(:)%nheightmax=17*ndof
      do i=1,n
         allocate (tpl(i)%icol(tpl(i)%nheightmax)) 
      enddo
      if (params%debug.gt.1) call heap (threadinfo,'tpl(:)%icol(:)','main',n*tpl(1)%nheightmax,'int',+1)
      !---[topology]----- 

      call wsmp_setup1 (n,n_iproc_st,n_iproc_end,n_iproc,nz_loc,ldb,ndof,vo,osolve,tpl,params,threadinfo,istep,iter)

      allocate(iproc_col(n),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'iproc_col','main',size(iproc_col),'bool',+1)
      allocate(ja(nz_loc),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'ja','main',size(ja),'int',+1)
      allocate(ia(n_iproc+1),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'ia','main',size(ia),'int',+1)
      allocate(avals(nz_loc),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'avals','main',size(avals),'dp',+1)
      allocate(b(ldb,nrhs),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'b','main',size(b),'dp',+1)

      call show_time (total,step,inc,1,'wsmp setup2$')

      call wsmp_setup2 (n,n_iproc,n_iproc_st,n_iproc_end,nz_loc,iproc_col,     &
                        ia,ja,ndof,vo,osolve,tpl,params,threadinfo,istep,iter)

      !---[topology]----- 
      if (params%debug.gt.1) call heap (threadinfo,'tpl(:)%icol(:)','main',n*tpl(1)%nheightmax,'int',-1)
      do i=1,n
         deallocate (tpl(i)%icol)
      enddo
      deallocate (tpl)
      !---[topology]----- 

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      ! allocate measurement arrays 
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|

      allocate (osolve%eviscosity(osolve%nleaves),stat=threadinfo%err)  
      if (params%debug.gt.1) call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',+1)
      allocate (osolve%q(osolve%nleaves),stat=threadinfo%err)           
      if (params%debug.gt.1) call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',+1)
      allocate (osolve%crit(osolve%nleaves),stat=threadinfo%err)        
      if (params%debug.gt.1) call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',+1)
      allocate (osolve%e2d(osolve%nleaves),stat=threadinfo%err)         
      if (params%debug.gt.1) call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',+1)
      allocate (osolve%e3d(osolve%nleaves),stat=threadinfo%err)         
      if (params%debug.gt.1) call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',+1)