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

program DOUAR

use threads
use definitions
!use mpi

implicit none

include 'mpif.h'

integer i,j,is,iter,istep,niter,nedge,nedgen 
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, external :: ioctree_number_of_elements

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

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

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=' '

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.'

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

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!     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
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
params%materialn(0)=material0
do i=1,params%ns
   params%materialn(i)=surface(i)%material
enddo

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!    compute plastic parameters 
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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)

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

   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   !     Are we in the main calculation phase or the spinup phase?
   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   if (istep.le.params%nstep_spinup) then
      params%dt=params%dt_spinup
      params%griditer=params%griditer_spinup
      params%nonlinear_iterations=params%nonlinear_iterations_spinup
   else
      params%dt=params%dt_main
      params%griditer=params%griditer_main
      params%nonlinear_iterations=params%nonlinear_iterations_main
   endif

   !---------------------------------------------------------------------------|
   !---------------------------------------------------------------------------|
   !     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
   if (istep.gt.1) params%init_e2d=.false.
   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

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

   iter=1

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

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

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     improve and update cloud fields 
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      call show_time (total,step,inc,1,'Change 3D cloud$')
      call update_cloud_structure (cl,osolve,params,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)

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      !     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)
      allocate (osolve%lode(osolve%nleaves),stat=threadinfo%err)
      if (params%debug.gt.1) call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',+1)
      allocate (osolve%is_plastic(osolve%nleaves),stat=threadinfo%err)  
      if (params%debug.gt.1) call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',+1)
      allocate (osolve%dilatr(osolve%nleaves),stat=threadinfo%err)         
      if (params%debug.gt.1) call heap (threadinfo,'osolve%dilatr','main',size(osolve%dilatr),'dp',+1)
      allocate (weightel(osolve%nleaves),stat=threadinfo%err)                               
      if (params%debug.gt.1) call heap (threadinfo,'weightel','main',size(weightel),'bool',+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%yield_ratio(osolve%nleaves),stat=threadinfo%err)  
      if (params%debug.gt.1) call heap (threadinfo,'osolve%yield_ratio','main',size(osolve%yield_ratio),'dp',+1)

      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      !     beginning of non linear iterations
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      do iter_nl=1,abs(params%nonlinear_iterations)                  
 
         if (iproc.eq.0) then
            write(*,*)  '-----------------------------------'
            write(8,*)  '-----------------------------------'
            call show_time(total,step,inc,-iter_nl,' start of non-linear iteration $')
            write(8,*)  'Doing inner iteration ',iter_nl
            write(*,*)  '-----------------------------------'
            write(8,*)  '-----------------------------------'
         end if
         if (params%debug.gt.1) call heap_hop3(threadinfo,iter_nl)
   
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         ! transfering previous solution from ov to osolve
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         if (iter_nl.gt.1) then    
            do i=1,ov%nnode             
               osolve%u(i)=ov%unode(i) 
               osolve%v(i)=ov%vnode(i)  
               osolve%w(i)=ov%wnode(i)
               osolve%wiso(i)=ov%wnodeiso(i)
            enddo  
         end if
!      call define_bc (params,osolve,vo)

         !if (iproc.eq.0) then
         !write(*,*) minval(osolve%u),maxval(osolve%u)
         !write(*,*) minval(osolve%v),maxval(osolve%v)
         !write(*,*) minval(osolve%w),maxval(osolve%w)
         !end if

         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         ! building the FEM matrix and rhs 
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         call show_time (total,step,inc,1,'build system$')

         call build_system_wsmp (n,n_iproc,n_iproc_st,nz_loc,iproc_col,ia,ja,  &
                                ldb,nrhs,avals,b,params,osolve,ndof,mat,vo,    &
                                threadinfo,weightel)

         ! Set reference strain rate flag to false after first nl iteration
         params%init_e2d=.false.

         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         ! if this is the first iteration of the first time step, write output file
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------

         if (istep.eq.1 .and. iter.eq.1 .and. iter_nl.eq.1) then
           call show_time (total,step,inc,1,'Write output initial step$')
           do is=1,osolve%nlsf
             allocate(surface(is)%u(surface(is)%nsurface))
             allocate(surface(is)%v(surface(is)%nsurface))
             allocate(surface(is)%w(surface(is)%nsurface))
             surface(is)%u=0.d0;surface(is)%v=0.d0;surface(is)%w=0.d0
           enddo
           call write_global_output (params,0,0,0.d0,osolve,ov,vo,surface,cl,bcdef,nest,'final')
           do is=1,osolve%nlsf
             deallocate (surface(is)%u,surface(is)%v,surface(is)%w)
           enddo
           call mpi_barrier (mpi_comm_world,ierr)
         endif

         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         ! solve system with wsmp solver 
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         call show_time (total,step,inc,1,'wsmp solve$')

         call solve_with_pwssmp (n,nz_loc,n_iproc,n_iproc_st,n_iproc_end,ia,ja,  &
                                 ldb,nrhs,avals,b,params,osolve,ov,vo,threadinfo,&
                                 ndof,istep,iter,iter_nl)

         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         !     calculate pressures
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         call show_time (total,step,inc,1,'Calculate pressures$')

         call compute_pressure (params,osolve,ov,vo,mat)

         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         !     smoothing pressures
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         call show_time (total,step,inc,1,'Smoothing pressures$')

         call smooth_pressures (osolve,ov,params)

         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         !     leaf measurements
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         call show_time (total,step,inc,1,'do leaf measurements$')

         call do_leaf_measurements (osolve,ov,params,istep,iter,iter_nl) 
 
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         !     write all informations in a text file
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         if (params%debug>=2) then
            call show_time (total,step,inc,1,'write global output$')
            call write_global_output(params,istep,iter,current_time,osolve,ov,vo,surface,cl,bcdef,nest,'debug')
            call mpi_barrier (mpi_comm_world,ierr)
!            call slices (params,osolve,ov,vo,sections,istep,iter,iter_nl)
         end if

         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         !     convergence criterion computation
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         call show_time (total,step,inc,1,'compute convergence criterion$')

         call compute_convergence_criterion(osolve,ov,vo,params,istep,iter,    &
                                           iter_nl,velocity_converged)

         cltemp=current_level  !!++!!

         if (velocity_converged .and. increase_current_level) current_level=min(current_level+1,params%levelmax_oct)

         if (iproc.eq.0 .and. params%debug>=1) then
            write(*,'(a,L1)') shift//'increase_current_level ->',increase_current_level
            write(*,'(a,i2)') shift//'current_level ->',current_level
         end if

         if (cltemp+1 == params%levelmax_oct) increase_current_level=.false. !!++!!

         if (params%nonlinear_iterations.lt.0 .and. velocity_converged) exit
      
      end do

      if (params%debug.gt.1) call heap_hop3_f (threadinfo)

      call DoRuRe_nonlin_stats(params%doDoRuRe,istep,iter,iter_nl)
 
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      !    end of non linear iterations
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------

      call show_time (total,step,inc,1,'slicing the cube$')
!      call slices(params,osolve,ov,vo,sections,istep,iter,iter_nl)

      !if (params%debug.gt.1) call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1)
      !deallocate (osolve%eviscosity)
      if (params%debug.gt.1) call heap (threadinfo,'osolve%q','main',size(osolve%q),'dp',-1)
      deallocate (osolve%q)                                                    

      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      !     estimating the divergence (to check if incompressibility has been respected)
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      call show_time (total,step,inc,1,'Compute divergence$')

      call compute_divergence (params,osolve,ov,vo,istep,iter)
      
      !-------------------------------------------------------------------------!
      !-------------------------------------------------------------------------!
      ! deallocate memory used by the solver and terminates the solver's job
      !-------------------------------------------------------------------------!
      !-------------------------------------------------------------------------!
      call show_time (total,step,inc,1,'wsmp_cleanup$')

      if (params%debug.gt.1) call heap (threadinfo,'ia','main',size(ia),'int',-1)
      deallocate(ia)
      if (params%debug.gt.1) call heap (threadinfo,'ja','main',size(ja),'int',-1)
      deallocate(ja)
      if (params%debug.gt.1) call heap (threadinfo,'iproc_col','main',size(iproc_col),'bool',-1)
      deallocate(iproc_col)
      if (params%debug.gt.1) call heap (threadinfo,'avals','main',size(avals),'dp',-1)
      deallocate(avals) 
      if (params%debug.gt.1) call heap (threadinfo,'b','main',size(b),'dp',-1)
      deallocate(b)
      !if (params%debug.gt.1) call heap (threadinfo,'weightel','main',size(weightel),'dp',-1)
      !deallocate(weightel)

      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      ! check for convergence
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------

      if (iter.eq.abs(params%griditer)) converge=.true.

      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      ! deallocate osolve
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      ! note that osolve will be needed in the temperature calculations
      ! so if this is the last iteration (converge is true), we will not deallocate osolve

      if (.not.converge) then
         if (params%debug.gt.1) call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',-1)
         deallocate (osolve%octree)
         if (params%debug.gt.1) call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',-1)
         deallocate (osolve%icon)                                     
         if (params%debug.gt.1) call heap (threadinfo,'osolve%x','main',size(osolve%x),'dp',-1)
         deallocate (osolve%x)    
         if (params%debug.gt.1) call heap (threadinfo,'osolve%y','main',size(osolve%y),'dp',-1)
         deallocate (osolve%y) 
         if (params%debug.gt.1) call heap (threadinfo,'osolve%z','main',size(osolve%z),'dp',-1)
         deallocate (osolve%z)
         if (params%debug.gt.1) call heap (threadinfo,'osolve%lsf','main',size(osolve%lsf),'dp',-1)
         deallocate (osolve%lsf) 
         if (params%debug.gt.1) call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',-1)
         deallocate (osolve%kfix)  
         if (params%debug.gt.1) call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',-1)
         deallocate (osolve%iface) 
         if (params%debug.gt.1) call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',-1)
         deallocate (osolve%strain)  
         if (params%debug.gt.1) call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',-1)
         deallocate (osolve%pressure)  
         if (params%debug.gt.1) call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',-1)
         deallocate (osolve%spressure)
         if (params%debug.gt.1) call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',-1)
         deallocate (osolve%kfixt) 
         if (params%debug.gt.1) call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',-1)
         deallocate (osolve%u) 
         if (params%debug.gt.1) call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',-1)
         deallocate (osolve%v) 
         if (params%debug.gt.1) call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',-1)
         deallocate (osolve%w)
         if (params%debug.gt.1) call heap (threadinfo,'osolve%wiso','main',size(osolve%wiso),'dp',-1)
         deallocate (osolve%wiso) 
         if (params%debug.gt.1) call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',-1)
         deallocate (osolve%temp) 
         if (params%debug.gt.1) call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',-1)
         deallocate (osolve%crit) 
         if (params%debug.gt.1) call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1)
         deallocate (osolve%e2d)  
         if (params%debug.gt.1) call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',-1)
         deallocate (osolve%e3d)  
         if (params%debug.gt.1) call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',-1)
         deallocate (osolve%lode)
         if (params%debug.gt.1) call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1)
         deallocate (osolve%eviscosity)
         if (params%debug.gt.1) call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',-1)
         deallocate (osolve%is_plastic)
         if (params%debug.gt.1) call heap (threadinfo,'osolve%dilatr','main',size(osolve%dilatr),'dp',-1)
         deallocate (osolve%dilatr)
         if (params%debug.gt.1) call heap (threadinfo,'osolve%matnum','main',size(osolve%matnum),'int',-1)
         deallocate (osolve%matnum)
         !if (params%debug.gt.1) call heap (threadinfo,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1)
         !deallocate (ov%wpreiso)
         if (params%debug.gt.1) call heap (threadinfo,'osolve%yield_ratio','main',size(osolve%yield_ratio),'dp',-1)
         deallocate (osolve%yield_ratio)
      end if

      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      ! finds courant condition time step in case this is the first iteration
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      if (iter.eq.1) then
         umax=0.d0
         do i=1,ov%nnode
            if (vo%influid(i)) umax=max(umax,sqrt(ov%unode(i)**2+ov%vnode(i)**2+ov%wnode(i)**2))
         enddo
         dtcourant=.5d0**params%levelmax_oct/umax*params%courant
         if (usecourant) then
            params%dt=dtcourant
            if(iproc.eq.0 .and. params%debug .ge. 1) write(*,'(a,es15.6)') shift//'dt (CFL) =',params%dt 
         end if

      endif

      if (converge) iter=abs(params%griditer)

      call DoRuRe_dt_stats (params%doDoRuRe,istep,params%dt)

      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      ! reset surface geometry to original geometry
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      call show_time(total,step,inc,1,'Reset surface geometry$')
      do is=1,params%ns
         surface(is)%r=surface0(is)%r
         surface(is)%s=surface0(is)%s
         surface(is)%x=surface0(is)%x
         surface(is)%y=surface0(is)%y
         surface(is)%z=surface0(is)%z
         surface(is)%xn=surface0(is)%xn
         surface(is)%yn=surface0(is)%yn
         surface(is)%zn=surface0(is)%zn
      enddo

      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      ! update surface geometry by midpoint rule if not converge
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------

      if (.not.converge) then

         olsf%nnode=ov%nnode
         olsf%nleaves=ov%nleaves
         olsf%noctree=ov%noctree
   
         allocate (olsf%octree(olsf%noctree),stat=threadinfo%err)             
         if (params%debug.gt.1) call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
         allocate (olsf%x(olsf%nnode),stat=threadinfo%err)                    
         if (params%debug.gt.1) call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
         allocate (olsf%y(olsf%nnode),stat=threadinfo%err)                    
         if (params%debug.gt.1) call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
         allocate (olsf%z(olsf%nnode),stat=threadinfo%err)                    
         if (params%debug.gt.1) call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
         allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)                  
         if (params%debug.gt.1) call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
         allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)    
         if (params%debug.gt.1) call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)

         olsf%octree(1:olsf%noctree)=ov%octree(1:ov%noctree)
         olsf%x(1:olsf%nnode)=ov%x(1:ov%nnode)
         olsf%y(1:olsf%nnode)=ov%y(1:ov%nnode)
         olsf%z(1:olsf%nnode)=ov%z(1:ov%nnode)
         olsf%icon(:,1:olsf%nleaves)=ov%icon(:,1:ov%nleaves)

         do is=1,params%ns
            if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
               call move_surface (surface(is),surface0(is),1,0,ov,params%dt/2.d0,params,istep,is)
            else
               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
            endif

            ! erosion added by Jean on Dec 12 2007
            if (material0.eq.0.and.params%erosion) then
               call int_to_char(ic,2,is)
               call show_time (total,step,inc,1,'Erode surface '//ic(1:2)//'$')
               if (current_time+tiny(current_time).ge. surface(is)%activation_time) then    
                  call erosion (surface(is),olsf,is,params)
               else
                  if (iproc.eq.0) write(8,*) 'Surface ',is,' is attached to Surface0'
                  surface(is)%z=surface(1)%z   
               endif
            end if
         enddo

         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         !     compute isostasy and move surfaces again
         !--------------------------------------------------------------------------------------
         !--------------------------------------------------------------------------------------
         if (params%isostasy) then
           call show_time (total,step,inc,1,'Compute isostasy and adjust vertical velocity$')
           !allocate(ov%wpreiso(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wpreiso in main$')
           call isostasy (params,weightel,ov,surface,mat,1,bcdef)

           call show_time (total,step,inc,1,'Apply isostasy to surfaces$')
           do is=1,params%ns
             if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
               call move_surface (surface(is),surface0(is),1,1,ov,params%dt/2.d0,params,istep,is)
             else
               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
             endif
           enddo
         endif

         if (params%debug.gt.1) call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',-1)
         deallocate (olsf%x)
         if (params%debug.gt.1) call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',-1)
         deallocate (olsf%y)
         if (params%debug.gt.1) call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',-1)
         deallocate (olsf%z)
         if (params%debug.gt.1) call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',-1)
         deallocate (olsf%lsf)
         if (params%debug.gt.1) call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',-1)
         deallocate (olsf%icon)
         if (params%debug.gt.1) call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',-1)
         deallocate (olsf%octree)
      endif

      if (.not.converge) then
        if (params%debug.gt.1) call heap (threadinfo,'weightel','main',size(weightel),'dp',-1)
        deallocate(weightel)
      endif

      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      ! deallocate void
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      ! note that vo will be needed in the temperature calculations
      ! so if this is the last iteration (converge is true), we will not deallocate vo

      if (.not.converge) then
         if (params%debug.gt.1) call heap (threadinfo,'vo%node','main',size(vo%node),'int',-1)
         deallocate (vo%node)
         if (params%debug.gt.1) call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',-1)
         deallocate (vo%leaf)
         if (params%debug.gt.1) call heap (threadinfo,'vo%face','main',size(vo%face),'int',-1)
         deallocate (vo%face)
         if (params%debug.gt.1) call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',-1)
         deallocate (vo%rtf)
         if (params%debug.gt.1) call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',-1)
         deallocate (vo%ftr)
         if (params%debug.gt.1) call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',-1)
         deallocate (vo%influid)
      end if

      if (converge) then
         do is=1,params%ns
            if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%r','main',size(surface0(is)%r),'dp',-1)
            deallocate (surface0(is)%r)
            if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%s','main',size(surface0(is)%s),'dp',-1)
            deallocate (surface0(is)%s)
            if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%x','main',size(surface0(is)%x),'dp',-1)
            deallocate (surface0(is)%x)
            if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%y','main',size(surface0(is)%y),'dp',-1)
            deallocate (surface0(is)%y)
            if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%z','main',size(surface0(is)%z),'dp',-1)
            deallocate (surface0(is)%z)
            if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%xn','main',size(surface0(is)%xn),'dp',-1)
            deallocate (surface0(is)%xn)
            if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%yn','main',size(surface0(is)%yn),'dp',-1)
            deallocate (surface0(is)%yn)
            if (params%debug.gt.1) call heap (threadinfo,'surface0(is)%zn','main',size(surface0(is)%zn),'dp',-1)
            deallocate (surface0(is)%zn)
         enddo
      endif

      iter=iter+1

   enddo 

   if (params%debug.gt.1) call heap_hop2_f (threadinfo)

   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   !     end of grid iterations
   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------

   ! we update current_time and can only do it now because this is where the time step is known
   current_time=current_time+params%dt
   if (iproc.eq.0) write(8,*) 'current time =', current_time

   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   !     solve for temperature
   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------

   if (iproc.eq.0) then
      write(*,*)  '-----------------------------------------------------------------------'
      write(8,*)  '-----------------------------------------------------------------------'
      call show_time (total,step,inc,1,'Temperature calculations $')
      write (8,*) 'Start of temperature calculations'
      write(*,*)  '-----------------------------------------------------------------------'
      write(8,*)  '-----------------------------------------------------------------------'
   end if

   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   ! transfers velocity and temperature solution from ov to osolve
   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------

   do i=1,ov%nnode
      osolve%u(i)=ov%unode(i)
      osolve%v(i)=ov%vnode(i)
      osolve%w(i)=ov%wnode(i)
      !osolve%wpreiso(i)=ov%wnodepreiso(i)
   enddo

   if (params%calculate_temp) then

      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|
      ! build matrix arrays, allocate memory for wsmp
      !------------------------------------------------------------------------|
      !------------------------------------------------------------------------|

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

      n = vo%nnode * ndoft
     
      !---[topology]----- 
      allocate (tpl(n)) 
      tpl%nheightmax=27*ndoft
      do i=1,n
         allocate (tpl(i)%icol(tpl(i)%nheightmax),stat=err) 
      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,ndoft,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)
      !allocate(weightel(osolve%nleaves),stat=threadinfo%err)
      !if (params%debug.gt.1) call heap (threadinfo,'weightel','main',size(weightel),'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,ndoft,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]----- 

      !----------------------------------------------------------------------------
      !----------------------------------------------------------------------------
      ! building the FEM matrix and rhs 
      !----------------------------------------------------------------------------
      !----------------------------------------------------------------------------

      call show_time (total,step,inc,1,'build system$')

      call build_system_wsmp (n,n_iproc,n_iproc_st,nz_loc,iproc_col,ia,ja,ldb, &
                             nrhs,avals,b,params,osolve,ndoft,mat,vo,          &
                             threadinfo,weightel)

      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------
      ! solve system with wsmp solver 
      !--------------------------------------------------------------------------------------
      !--------------------------------------------------------------------------------------

      call solve_with_pwgsmp (n,nz_loc,n_iproc,n_iproc_st,n_iproc_end,ia,ja,ldb,&
                              nrhs,avals,b,params,osolve,ov,vo,threadinfo,ndoft,&
                              istep,iter,iter_nl)

      !-------------------------------------------------------------------------!
      !-------------------------------------------------------------------------!
      ! deallocate memory used by the solver and terminates the solver's job
      !-------------------------------------------------------------------------!
      !-------------------------------------------------------------------------!
      call show_time (total,step,inc,1,'wsmp_cleanup$')

      if (params%debug.gt.1) call heap (threadinfo,'ia','main',size(ia),'int',-1)
      deallocate(ia)
      if (params%debug.gt.1) call heap (threadinfo,'ja','main',size(ja),'int',-1)
      deallocate(ja)
      if (params%debug.gt.1) call heap (threadinfo,'iproc_col','main',size(iproc_col),'bool',-1)
      deallocate(iproc_col)
      if (params%debug.gt.1) call heap (threadinfo,'avals','main',size(avals),'dp',-1)
      deallocate(avals) 
      if (params%debug.gt.1) call heap (threadinfo,'b','main',size(b),'dp',-1)
      deallocate(b)
      !if (params%debug.gt.1) call heap (threadinfo,'weightel','main',size(weightel),'dp',-1)
      !deallocate(weightel)

   else

      if(iproc.eq.0) then
         write(*,'(a,f15.7,a)') shift//'================================='
         write(*,'(a)')         shift//'skip temperature calculation'
         write(*,'(a,f15.7,a)') shift//'================================='
      end if

   end if

   !-------------------------------------------------------------------------!
   !-------------------------------------------------------------------------!

   call DoRuRe_temp_stats(params%doDoRuRe,istep,ov%nnode,ov%temp)

   olsf%nnode=ov%nnode
   olsf%nleaves=ov%nleaves
   olsf%noctree=ov%noctree
   
   allocate (olsf%octree(olsf%noctree),stat=threadinfo%err)
   if (params%debug.gt.1) call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
   allocate (olsf%x(olsf%nnode),stat=threadinfo%err)
   if (params%debug.gt.1) call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
   allocate (olsf%y(olsf%nnode),stat=threadinfo%err)
   if (params%debug.gt.1) call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
   allocate (olsf%z(olsf%nnode),stat=threadinfo%err)
   if (params%debug.gt.1) call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
   allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)
   if (params%debug.gt.1) call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
   allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)    
   if (params%debug.gt.1) call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',+1)

   olsf%octree(1:olsf%noctree)=ov%octree(1:ov%noctree)
   olsf%x(1:olsf%nnode)=ov%x(1:ov%nnode)
   olsf%y(1:olsf%nnode)=ov%y(1:ov%nnode)
   olsf%z(1:olsf%nnode)=ov%z(1:ov%nnode)
   olsf%icon(:,1:olsf%nleaves)=ov%icon(:,1:ov%nleaves)

   ! write time step to Log file
   if (iproc.eq.0) then
      write (8,*) 'Current time step ',params%dt
      write (8,*) 'Courant time step ',dtcourant
   end if

   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   !     move surfaces and apply erosion
   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   call show_time (total,step,inc,1,'Advect the surfaces$')

   do is=1,params%ns
      write (ic(1:2),'(i2)') is
      call show_time (total,step,inc,1,'Advect surface '//ic(1:2)//'$')
      if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
         call move_surface (surface(is),surface0(is),0,0,ov,params%dt,params,istep,is)
      else
         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
      endif

      ! We apply erosion
      if (material0.eq.0.and.params%erosion) then
         call show_time (total,step,inc,1,'Erode surface '//ic(1:2)//'$')
         if (current_time+tiny(current_time).ge. surface(is)%activation_time) then    
            call erosion (surface(is),olsf,is,params)
         else
            if (iproc.eq.0) write(8,*) 'Surface ',is,' is attached to Surface0'
            surface(is)%z=surface(1)%z   
         endif
      end if
      if (iproc.eq.0) write (8,*) 'Min-Max z surf ',is,':',minval(surface(is)%z),maxval(surface(is)%z)

      allocate(surface(is)%u(surface(is)%nsurface))
      allocate(surface(is)%v(surface(is)%nsurface))
      allocate(surface(is)%w(surface(is)%nsurface))

   enddo

   !call interpolate_velocity_on_surface(params,surface,ov)

   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   !     compute isostasy, move surfaces again, update Moho displacement
   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------

   if (params%isostasy) then
     call show_time (total,step,inc,1,'Compute isostasy and adjust vertical velocity$')
     !allocate(ov%wpreiso(ov%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ov%wpreiso in main$')
     call isostasy (params,weightel,ov,surface,mat,0,bcdef)
     if (params%isobc) bcdef%zisodisp=bcdef%zisodisp-bcdef%zisoinc

     do is=1,params%ns
       if (current_time+tiny(current_time).ge. surface(is)%activation_time) then
         call move_surface (surface(is),surface0(is),0,1,ov,params%dt,params,istep,is)
       else
         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
       endif
     enddo
   endif

   ! Done with isostasy, deallocate weightel
   if (params%debug.gt.1) call heap (threadinfo,'weightel','main',size(weightel),'dp',-1)
   deallocate(weightel)

   call interpolate_velocity_on_surface(params,surface,ov)

   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   !     update strain, pressure, temperature on 3D cloud
   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   call show_time (total,step,inc,1,'Update cloud fields$')
   call update_cloud_fields (cl,ov,osolve,vo,params)

   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   !     advect cloud
   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   call show_time (total,step,inc,1,'Advect the cloud$')
!   call move_cloud (params,cl,ov%octree,ov%noctree,ov%unode,ov%vnode,ov%wnode,ov%nnode,ov%icon,ov%nleaves)

   call move_cloud (params,cl,ov)

   if (params%debug.gt.1) call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',-1)
   deallocate (olsf%x) 
   if (params%debug.gt.1) call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',-1)
   deallocate (olsf%y) 
   if (params%debug.gt.1) call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',-1)
   deallocate (olsf%z)
   if (params%debug.gt.1) call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',-1)
   deallocate (olsf%lsf)
   if (params%debug.gt.1) call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',-1)
   deallocate (olsf%icon)
   if (params%debug.gt.1) call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',-1)
   deallocate (olsf%octree)

   osolve%u=ov%unode
   osolve%v=ov%vnode
   osolve%w=ov%wnode
   osolve%wiso=ov%wnodeiso
   osolve%temp=ov%temp

   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   !     write global output before deallocating arrays
   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   call show_time (total,step,inc,1,'Write output$')
   call write_global_output (params,istep,iter-1,current_time,osolve,ov,vo,surface,cl,bcdef,nest,'final')
   call mpi_barrier (mpi_comm_world,ierr)

   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   ! deallocate osolve and vo
   !--------------------------------------------------------------------------------------
   !--------------------------------------------------------------------------------------
   do is=1,params%ns
      deallocate(surface(is)%u)
      deallocate(surface(is)%v)
      deallocate(surface(is)%w)
   end do
  
   if (params%debug.gt.1) call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',-1)
   deallocate (osolve%octree)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',-1)
   deallocate (osolve%icon)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%x','main',size(osolve%x),'dp',-1)
   deallocate (osolve%x)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%y','main',size(osolve%y),'dp',-1)
   deallocate (osolve%y)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%z','main',size(osolve%z),'dp',-1)
   deallocate (osolve%z)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%lsf','main',size(osolve%lsf),'dp',-1)
   deallocate (osolve%lsf)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',-1)
   deallocate (osolve%kfix)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',-1)
   deallocate (osolve%iface)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',-1)
   deallocate (osolve%pressure)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',-1)
   deallocate (osolve%spressure)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',-1)
   deallocate (osolve%strain)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',-1)
   deallocate (osolve%kfixt)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',-1)
   deallocate (osolve%u)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',-1)
   deallocate (osolve%v)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',-1)
   deallocate (osolve%w)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%wiso','main',size(osolve%wiso),'dp',-1)
   deallocate (osolve%wiso)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',-1)
   deallocate (osolve%temp)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',-1)
   deallocate (osolve%crit)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1)
   deallocate (osolve%e2d)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',-1)
   deallocate (osolve%e3d)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',-1)
   deallocate (osolve%lode)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1)
   deallocate (osolve%eviscosity)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',-1)
   deallocate (osolve%is_plastic)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%dilatr','main',size(osolve%dilatr),'dp',-1)
   deallocate (osolve%dilatr)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%matnum','main',size(osolve%matnum),'int',-1)
   deallocate (osolve%matnum)
   if (params%debug.gt.1) call heap (threadinfo,'osolve%yield_ratio','main',size(osolve%yield_ratio),'dp',-1)
   deallocate (osolve%yield_ratio)
   if (params%debug.gt.1) call heap (threadinfo,'vo%node','main',size(vo%node),'int',-1)
   deallocate (vo%node)
   if (params%debug.gt.1) call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',-1)
   deallocate (vo%leaf)
   if (params%debug.gt.1) call heap (threadinfo,'vo%face','main',size(vo%face),'int',-1)
   deallocate (vo%face)
   if (params%debug.gt.1) call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',-1)
   deallocate (vo%rtf)
   if (params%debug.gt.1) call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',-1)
   deallocate (vo%ftr)
   if (params%debug.gt.1) call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',-1)
   deallocate (vo%influid)
   !if (params%debug.gt.1) call heap (threadinfo,'ov%wpreiso','main',size(ov%wpreiso),'dp',-1)
   !deallocate (ov%wpreiso)

   call show_time (total,step,inc,1,'End of time step$')

   ! end of big time loop

   istep=istep+1

enddo
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
!     end of time stepping
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
if (params%debug.gt.1) call heap_hop1_f (threadinfo)

call show_time (total,step,inc,1,'End of run$')

call io_DoRuRe2 (params,'close')
   
deallocate (cl%x)
deallocate (cl%y)
deallocate (cl%z)
deallocate (cl%x0)
deallocate (cl%y0)
deallocate (cl%z0)
deallocate (cl%strain)
deallocate (cl%lsf0)
deallocate (cl%temp)
deallocate (cl%press)
deallocate (cl%tag)

if (params%debug.gt.1) call heap (threadinfo,'ov%octree','main',size(ov%octree),'int',-1)
deallocate (ov%octree)
if (params%debug.gt.1) call heap (threadinfo,'ov%x','main',size(ov%x),'dp',-1)
deallocate (ov%x)
if (params%debug.gt.1) call heap (threadinfo,'ov%y','main',size(ov%y),'dp',-1)
deallocate (ov%y)
if (params%debug.gt.1) call heap (threadinfo,'ov%z','main',size(ov%z),'dp',-1)
deallocate (ov%z)
if (params%debug.gt.1) call heap (threadinfo,'ov%unode','main',size(ov%unode),'dp',-1)
deallocate (ov%unode)
if (params%debug.gt.1) call heap (threadinfo,'ov%vnode','main',size(ov%vnode),'dp',-1)
deallocate (ov%vnode)
if (params%debug.gt.1) call heap (threadinfo,'ov%wnode','main',size(ov%wnode),'dp',-1)
deallocate (ov%wnode)
if (params%debug.gt.1) call heap (threadinfo,'ov%wnodeiso','main',size(ov%wnodeiso),'dp',-1)
deallocate (ov%wnodeiso)
if (params%debug.gt.1) call heap (threadinfo,'ov%icon','main',size(ov%icon),'int',-1)
deallocate (ov%icon)
if (params%debug.gt.1) call heap (threadinfo,'ov%temp','main',size(ov%temp),'dp',-1)
deallocate (ov%temp)
if (params%debug.gt.1) call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pressure),'dp',-1) 
deallocate (ov%temporary_nodal_pressure)
if (params%debug.gt.1) call heap (threadinfo,'ov%whole_leaf...','main',size(ov%whole_leaf_in_fluid),'bool',-1)         
deallocate (ov%whole_leaf_in_fluid)

do is=1,params%ns
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%x','main',size(surface(is)%x),'dp',-1)
   deallocate (surface(is)%x)
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%y','main',size(surface(is)%y),'dp',-1)
   deallocate (surface(is)%y)
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%z','main',size(surface(is)%z),'dp',-1)
   deallocate (surface(is)%z)
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%xn','main',size(surface(is)%xn),'dp',-1)
   deallocate (surface(is)%xn)
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%yn','main',size(surface(is)%yn),'dp',-1)
   deallocate (surface(is)%yn)
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%zn','main',size(surface(is)%zn),'dp',-1)
   deallocate (surface(is)%zn)
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%r','main',size(surface(is)%r),'dp',-1)
   deallocate (surface(is)%r)
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%s','main',size(surface(is)%s),'dp',-1)
   deallocate (surface(is)%s)
   if (params%debug.gt.1) call heap (threadinfo,'surface(is)%icon','main',size(surface(is)%icon),'int',-1)
   deallocate (surface(is)%icon)
end do

deallocate (surface)
deallocate (surface0)
deallocate (mat)

if (params%isobc) then
  deallocate (bcdef%zisodisp)
  deallocate (bcdef%zisoinc)
endif
deallocate (params%materialn)

if (params%debug.gt.1) then
  call heap_final (threadinfo)
  close (threadinfo%Logunit)
  close (threadinfo%mem_heap_unit)
endif

call mpi_finalize (ierr)

end
