DOUAR.f90

      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
         call heap (threadinfo,'vo%node','main',size(vo%node),'int',-1)        ; deallocate (vo%node)
         call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',-1)        ; deallocate (vo%leaf)
         call heap (threadinfo,'vo%face','main',size(vo%face),'int',-1)        ; deallocate (vo%face)
         call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',-1)          ; deallocate (vo%rtf)
         call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',-1)          ; deallocate (vo%ftr)
         call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',-1) ; deallocate (vo%influid)
      end if

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

      iter=iter+1

   enddo 

   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
      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)         ; call heap (threadinfo,'iproc_col','main',size(iproc_col),'bool',+1)
      allocate(ja(nz_loc),stat=threadinfo%err)           ; call heap (threadinfo,'ja','main',size(ja),'int',+1)
      allocate(ia(n_iproc+1),stat=threadinfo%err)        ; call heap (threadinfo,'ia','main',size(ia),'int',+1)
      allocate(avals(nz_loc),stat=threadinfo%err)        ; call heap (threadinfo,'avals','main',size(avals),'dp',+1)
      allocate(b(ldb,nrhs),stat=threadinfo%err)          ; call heap (threadinfo,'b','main',size(b),'dp',+1)
      !allocate(weightel(osolve%nleaves),stat=threadinfo%err)   ; 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]----- 
      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,istep,iter,iter_nl,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$')

      call heap (threadinfo,'ia','main',size(ia),'int',-1)                ; deallocate(ia)
      call heap (threadinfo,'ja','main',size(ja),'int',-1)                ; deallocate(ja)
      call heap (threadinfo,'iproc_col','main',size(iproc_col),'bool',-1) ; deallocate(iproc_col)
      call heap (threadinfo,'avals','main',size(avals),'dp',-1)           ; deallocate(avals) 
      call heap (threadinfo,'b','main',size(b),'dp',-1)                   ; deallocate(b)
      !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)             ; call heap (threadinfo,'olsf%octree','main',size(olsf%octree),'int',+1)
   allocate (olsf%x(olsf%nnode),stat=threadinfo%err)                    ; call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',+1)
   allocate (olsf%y(olsf%nnode),stat=threadinfo%err)                    ; call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',+1)
   allocate (olsf%z(olsf%nnode),stat=threadinfo%err)                    ; call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',+1)
   allocate (olsf%lsf(olsf%nnode),stat=threadinfo%err)                  ; call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',+1)
   allocate (olsf%icon(params%mpe,olsf%nleaves),stat=threadinfo%err)    
   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,zi)
     if (params%isobc) zi%zisodisp=zi%zisodisp-zi%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
   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)

   call heap (threadinfo,'olsf%x','main',size(olsf%x),'dp',-1)            ; deallocate (olsf%x) 
   call heap (threadinfo,'olsf%y','main',size(olsf%y),'dp',-1)            ; deallocate (olsf%y) 
   call heap (threadinfo,'olsf%z','main',size(olsf%z),'dp',-1)            ; deallocate (olsf%z)
   call heap (threadinfo,'olsf%lsf','main',size(olsf%lsf),'dp',-1)        ; deallocate (olsf%lsf)
   call heap (threadinfo,'olsf%icon','main',size(olsf%icon),'int',-1)     ; deallocate (olsf%icon)
   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,zi,'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
  
   call heap (threadinfo,'osolve%octree','main',size(osolve%octree),'int',-1)           ; deallocate (osolve%octree)
   call heap (threadinfo,'osolve%icon','main',size(osolve%icon),'int',-1)               ; deallocate (osolve%icon)
   call heap (threadinfo,'osolve%x','main',size(osolve%x),'dp',-1)                      ; deallocate (osolve%x)
   call heap (threadinfo,'osolve%y','main',size(osolve%y),'dp',-1)                      ; deallocate (osolve%y)
   call heap (threadinfo,'osolve%z','main',size(osolve%z),'dp',-1)                      ; deallocate (osolve%z)
   call heap (threadinfo,'osolve%lsf','main',size(osolve%lsf),'dp',-1)                  ; deallocate (osolve%lsf)
   call heap (threadinfo,'osolve%kfix','main',size(osolve%kfix),'int',-1)               ; deallocate (osolve%kfix)
   call heap (threadinfo,'osolve%iface','main',size(osolve%iface),'int',-1)             ; deallocate (osolve%iface)
   call heap (threadinfo,'osolve%pressure','main',size(osolve%pressure),'dp',-1)        ; deallocate (osolve%pressure)
   call heap (threadinfo,'osolve%spressure','main',size(osolve%spressure),'dp',-1)      ; deallocate (osolve%spressure)
   call heap (threadinfo,'osolve%strain','main',size(osolve%strain),'dp',-1)            ; deallocate (osolve%strain)
   call heap (threadinfo,'osolve%kfixt','main',size(osolve%kfixt),'int',-1)             ; deallocate (osolve%kfixt)
   call heap (threadinfo,'osolve%u','main',size(osolve%u),'dp',-1)                      ; deallocate (osolve%u)
   call heap (threadinfo,'osolve%v','main',size(osolve%v),'dp',-1)                      ; deallocate (osolve%v)
   call heap (threadinfo,'osolve%w','main',size(osolve%w),'dp',-1)                      ; deallocate (osolve%w)
   call heap (threadinfo,'osolve%wiso','main',size(osolve%wiso),'dp',-1)                ; deallocate (osolve%wiso)
   call heap (threadinfo,'osolve%temp','main',size(osolve%temp),'dp',-1)                ; deallocate (osolve%temp)
   call heap (threadinfo,'osolve%crit','main',size(osolve%crit),'dp',-1)                ; deallocate (osolve%crit)
   call heap (threadinfo,'osolve%e2d','main',size(osolve%e2d),'dp',-1)                  ; deallocate (osolve%e2d)
   call heap (threadinfo,'osolve%e3d','main',size(osolve%e3d),'dp',-1)                  ; deallocate (osolve%e3d)
   call heap (threadinfo,'osolve%lode','main',size(osolve%lode),'dp',-1)                ; deallocate (osolve%lode)
   call heap (threadinfo,'osolve%eviscosity','main',size(osolve%eviscosity),'dp',-1)    ; deallocate (osolve%eviscosity)
   call heap (threadinfo,'osolve%is_plastic','main',size(osolve%is_plastic),'bool',-1)  ; deallocate (osolve%is_plastic)
   call heap (threadinfo,'vo%node','main',size(vo%node),'int',-1)                       ; deallocate (vo%node)
   call heap (threadinfo,'vo%leaf','main',size(vo%leaf),'int',-1)                       ; deallocate (vo%leaf)
   call heap (threadinfo,'vo%face','main',size(vo%face),'int',-1)                       ; deallocate (vo%face)
   call heap (threadinfo,'vo%rtf','main',size(vo%rtf),'int',-1)                         ; deallocate (vo%rtf)
   call heap (threadinfo,'vo%ftr','main',size(vo%ftr),'int',-1)                         ; deallocate (vo%ftr)
   call heap (threadinfo,'vo%influid','main',size(vo%influid),'bool',-1)                ; deallocate (vo%influid)
   !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
!--------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------
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)

call heap (threadinfo,'ov%octree','main',size(ov%octree),'int',-1)   ; deallocate (ov%octree)
call heap (threadinfo,'ov%x','main',size(ov%x),'dp',-1)              ; deallocate (ov%x)
call heap (threadinfo,'ov%y','main',size(ov%y),'dp',-1)              ; deallocate (ov%y)
call heap (threadinfo,'ov%z','main',size(ov%z),'dp',-1)              ; deallocate (ov%z)
call heap (threadinfo,'ov%unode','main',size(ov%unode),'dp',-1)      ; deallocate (ov%unode)
call heap (threadinfo,'ov%vnode','main',size(ov%vnode),'dp',-1)      ; deallocate (ov%vnode)
call heap (threadinfo,'ov%wnode','main',size(ov%wnode),'dp',-1)      ; deallocate (ov%wnode)
call heap (threadinfo,'ov%wnodeiso','main',size(ov%wnodeiso),'dp',-1); deallocate (ov%wnodeiso)
call heap (threadinfo,'ov%icon','main',size(ov%icon),'int',-1)       ; deallocate (ov%icon)
call heap (threadinfo,'ov%temp','main',size(ov%temp),'dp',-1)        ; deallocate (ov%temp)
call heap (threadinfo,'ov%temporary_nodal...','main',size(ov%temporary_nodal_pressure),'dp',-1) 
 deallocate (ov%temporary_nodal_pressure)
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
   call heap (threadinfo,'surface(is)%x','main',size(surface(is)%x),'dp',-1)         ; deallocate (surface(is)%x)
   call heap (threadinfo,'surface(is)%y','main',size(surface(is)%y),'dp',-1)         ; deallocate (surface(is)%y)
   call heap (threadinfo,'surface(is)%z','main',size(surface(is)%z),'dp',-1)         ; deallocate (surface(is)%z)
   call heap (threadinfo,'surface(is)%xn','main',size(surface(is)%xn),'dp',-1)       ; deallocate (surface(is)%xn)
   call heap (threadinfo,'surface(is)%yn','main',size(surface(is)%yn),'dp',-1)       ; deallocate (surface(is)%yn)
   call heap (threadinfo,'surface(is)%zn','main',size(surface(is)%zn),'dp',-1)       ; deallocate (surface(is)%zn)
   call heap (threadinfo,'surface(is)%r','main',size(surface(is)%r),'dp',-1)         ; deallocate (surface(is)%r)
   call heap (threadinfo,'surface(is)%s','main',size(surface(is)%s),'dp',-1)         ; deallocate (surface(is)%s)
   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) deallocate (zisodisp)
if (params%isobc) then
  deallocate (zi%zisodisp)
  deallocate (zi%zisoinc)
endif
deallocate (params%materialn)

call heap_final (threadinfo)

close (threadinfo%Logunit)
close (threadinfo%mem_heap_unit)

call mpi_finalize (ierr)

end