!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              IMPROVE_OSOLVE    Feb. 2007                                     |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine improve_osolve (osolve,ov,params,threadinfo,istep,iter,total,step,  &
                           inc,refine_level,boxes,cube_faces,                  &
                           increase_current_level)
      
!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine calculates the second invariant of the strain rate, and a value 
! for the criterion used by the improve_octree subroutine, inside each element
! by using the velocity information at the nodes. Then, the routine improves the
! osolve on which the solution is calculated. Also, it refines the octree if the
! user has defined a box (or several boxes) in which the octree is artificially
! refined to a desired level. boxes and nboxes are read from the input file. 
! The routine then refines user supplied rectangular areas of each of the six
! faces of the cubic simulation domain. The information are stored in cube_faces
! read from the input file. Finally, the octree is smoothed and super-smoothed
! if the user has set ismooth to 1
!
! ov is the velocity octree containing the velocity solution
! osolve is the octree to be improved
! refine_ratio is a parameter read in input.txt. When multiplied by the 
! maximum of the criterion previously computed, one obtains the threshold 
! used to determine whether a leaf is to be refined or not.  
! refine_level is the level at which the osolve octree is to be refined 

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|



!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use definitions
use gauss
use invariants
!use mpi
use threads

implicit none

include 'mpif.h'

type (octreesolve) osolve
type (octreev) ov
type (parameters) params
type (thread) threadinfo
integer istep,iter
double precision total,step,inc
integer refine_level
type (box) boxes(params%nboxes)
type (face), dimension(6) ::  cube_faces
logical increase_current_level

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

double precision jcb(3,3),jcbi(3,3)
double precision,dimension(:),allocatable :: x,y,z,vx,vy,vz
double precision,dimension(:),allocatable :: dhdx,dhdy,dhdz
double precision,dimension(:),allocatable :: crit,crittemp,xc,yc,zc
double precision :: eps,r,s,t,volume,w
double precision :: exx,eyy,ezz,exy,eyz,ezx
double precision critmax,xxx,yyy,zzz,b,l
double precision :: x0,y0,z0,x1,y1,z1,E2d
integer err,ierr,nproc,iproc,k,i,iint,nint,mpe
integer :: j,nn,nb_of_boxes,level_box,nx,ny,nz,level,nv,nh
character*72 shift

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

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

eps=tiny(eps)

nint=8
mpe=8
shift=' '

!------------------------------------------------------------------
! first compute the refinement criterion measured on the ov octree
!------------------------------------------------------------------
      
call show_time (total,step,inc,1,'compute ref. criterion on ov$')

allocate (x(mpe),y(mpe),z(mpe),stat=threadinfo%err)            
if (params%debug.gt.1) call heap (threadinfo,'x/y/z','improve_osolve',size(x)+size(y)+size(z),'dp',+1)
allocate (vx(mpe),vy(mpe),vz(mpe),stat=threadinfo%err)         
if (params%debug.gt.1) call heap (threadinfo,'vx/y/z','improve_osolve',size(vx)+size(vy)+size(vz),'dp',+1)
allocate (dhdx(mpe),dhdy(mpe),dhdz(mpe),stat=threadinfo%err)   
if (params%debug.gt.1) call heap (threadinfo,'dhdx/y/z','improve_osolve',size(dhdx)+size(dhdy)+size(dhdz),'dp',+1)
allocate (crit(ov%nleaves),stat=threadinfo%err)                
if (params%debug.gt.1) call heap (threadinfo,'crit','improve_osolve',size(crit),'dp',+1)
allocate (crittemp(ov%nleaves),stat=threadinfo%err)            
if (params%debug.gt.1) call heap (threadinfo,'crittemp','improve_osolve',size(crittemp),'dp',+1)

crit=0.d0
crittemp=0.d0

r=0.d0
s=0.d0
t=0.d0

do i=1+iproc,ov%nleaves,nproc

   do k=1,mpe
      x(k)=ov%x(ov%icon(k,i))
      y(k)=ov%y(ov%icon(k,i))
      z(k)=ov%z(ov%icon(k,i))
      vx(k)=ov%unode(ov%icon(k,i))
      vy(k)=ov%vnode(ov%icon(k,i))
      vz(k)=ov%wnode(ov%icon(k,i))
   enddo

   call compute_dhdx_dhdy_dhdz (mpe,r,s,t,x,y,z,dhdx,dhdy,dhdz,volume)

   exx=0.d0
   eyy=0.d0
   ezz=0.d0
   exy=0.d0
   eyz=0.d0
   ezx=0.d0
   do k=1,mpe
      exx=exx+dhdx(k)*vx(k)
      eyy=eyy+dhdy(k)*vy(k)
      ezz=ezz+dhdz(k)*vz(k)
      exy=exy+(dhdx(k)*vy(k)+dhdy(k)*vx(k))/2.d0
      eyz=eyz+(dhdy(k)*vz(k)+dhdz(k)*vy(k))/2.d0
      ezx=ezx+(dhdz(k)*vx(k)+dhdx(k)*vz(k))/2.d0
   enddo

   call deviatoric (params%invariants_2d,exx,eyy,ezz,exy,eyz,ezx)

   E2d=second_invariant(exx,eyy,ezz,exy,eyz,ezx)

   select case (params%refine_criterion)
      case(1) 
         crittemp(i)=sqrt(E2d)
!         crittemp(i)=log(1.d0+sqrt(E2d))
      case(2) 
         crittemp(i)=sqrt(exx**2+eyy**2+ezz**2)
      case(3)
         crittemp(i)=sqrt(E2d)*abs(x(1)-x(2)) 
      case default
         crittemp(i)=0.d0
   end select

end do

call mpi_allreduce (crittemp,crit,ov%nleaves,mpi_double_precision,mpi_sum,mpi_comm_world,ierr)

if (params%debug.gt.1) call heap (threadinfo,'x/y/z','improve_osolve',size(x)+size(y)+size(z),'dp',-1)
deallocate (x,y,z)            
if (params%debug.gt.1) call heap (threadinfo,'vx/vy/vz','improve_osolve',size(vx)+size(vy)+size(vz),'dp',-1)
deallocate (vx,vy,vz)         
if (params%debug.gt.1) call heap (threadinfo,'dhdx/y/z','improve_osolve',size(dhdx)+size(dhdy)+size(dhdz),'dp',-1)
deallocate (dhdx,dhdy,dhdz)   
if (params%debug.gt.1) call heap (threadinfo,'crittemp','improve_osolve',size(crittemp),'dp',-1)
deallocate (crittemp)         

!---------------------------------------------------
! we can now use crit to refine the osolve octree
!---------------------------------------------------

if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)

call show_time (total,step,inc,1,'improve osolve for criterion$')

if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,es15.6)') shift//'max(crit)        =',maxval(crit)
if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,es15.6)') shift//'max(crit,w_l_i_f)=',maxval(crit,ov%whole_leaf_in_fluid)

if(istep.eq.1 .and. iter.eq.1) then
   critmax=maxval(crit)
else
   critmax=maxval(crit,ov%whole_leaf_in_fluid)
end if

if (critmax.gt.tiny(critmax)) then
   do i=1,ov%nleaves


!if(increase_current_level) then
!      if (crit(i).gt. critmax*2.d0*params%refine_ratio  ) then
!         xxx=0.5d0*(ov%x(ov%icon(1,i))+ov%x(ov%icon(8,i)))
!         yyy=0.5d0*(ov%y(ov%icon(1,i))+ov%y(ov%icon(8,i)))
!         zzz=0.5d0*(ov%z(ov%icon(1,i))+ov%z(ov%icon(8,i)))
!         if (ov%whole_leaf_in_fluid(i)) &
!         call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,refine_level)
!      end if
!else
      if (crit(i).gt. critmax*params%refine_ratio  ) then
         xxx=0.5d0*(ov%x(ov%icon(1,i))+ov%x(ov%icon(8,i)))
         yyy=0.5d0*(ov%y(ov%icon(1,i))+ov%y(ov%icon(8,i)))
         zzz=0.5d0*(ov%z(ov%icon(1,i))+ov%z(ov%icon(8,i)))
         if (ov%whole_leaf_in_fluid(i)) &
         call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,refine_level)
      end if
!end if


   enddo
end if

if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)

if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after crit refinement   ',osolve%octree(2)


if (params%debug.gt.1) call heap (threadinfo,'crit','improve_osolve',size(crit),'dp',-1)
deallocate (crit)  
    
!-----------------------------------------------
! osolve has yet to be refined in boxes, if any
!-----------------------------------------------

if (params%nboxes.gt.0) then
   call show_time (total,step,inc,1,'Improve osolve for imposed boxes$')

   do nn=1,params%nboxes
      x0=boxes(nn)%x0
      x1=boxes(nn)%x1
      y0=boxes(nn)%y0
      y1=boxes(nn)%y1
      z0=boxes(nn)%z0
      z1=boxes(nn)%z1
      level_box=boxes(nn)%level
      nx=int((x1-x0)/2.d0**(-level_box))! nb of box particles to insert on the x axis 
      ny=int((y1-y0)/2.d0**(-level_box))! nb of box particles to insert on the y axis
      nz=int((z1-z0)/2.d0**(-level_box)) ! nb of box particles to insert on the z axis
      do i=1,nx
         do j=1,ny
            do k=1,nz
               xxx=x0+(real(i)-.5d0)*2.d0**(-level_box)
               yyy=y0+(real(j)-.5d0)*2.d0**(-level_box)
               zzz=z0+(real(k)-.5d0)*2.d0**(-level_box)
               if (.not.((xxx>=x0).and.(xxx<=x1) .and. &
                         (yyy>=y0).and.(yyy<=y1) .and. &
                         (zzz>=z0).and.(zzz<=z1))) call stop_run ('pb with box cloud$')
               ! if (iproc.eq.0) print*,i,j,k
               call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level_box)
            enddo
         enddo
      enddo
   enddo
   if (params%debug>=1 .and. iproc.eq.0) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
   if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after box refinement   ',osolve%octree(2)
endif


!-----------------------------------------------
! refine the octree on the cube faces
!----------------------------------------------
if (params%ref_on_faces) then

   call show_time (total,step,inc,1,'Improve osolve on cube faces$')

   do k=1,6

!   if (k==5) then
!      level=refine_level
!   else
     level=cube_faces(k)%level
!   end if

   l=cube_faces(k)%l
   r=cube_faces(k)%r
   t=cube_faces(k)%t
   b=cube_faces(k)%b

   nv=int((t-b)/2.d0**(-(level+1))) ! nb of box particles to insert on the vertical axis 
   nh=int((r-l)/2.d0**(-(level+1))) ! nb of box particles to insert on the horizontal axis
   
   if (l > r .or. b > t) call stop_run ('pb in improve_octree_on_faces$')
   
   select case (k)
      case (1)  ! x=0 face
        xxx=0.00001
        do i=1,nh
           do j=1,nv
              yyy=l+(real(i)-.5d0)*2.d0**(-(level+1))
              zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
              call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
           enddo
        enddo

      case (2)  ! x=1 face
        xxx=0.99999
        do i=1,nh
           do j=1,nv
              yyy=l+(real(i)-.5d0)*2.d0**(-(level+1))
              zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
              call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
           enddo
        enddo

      case (3)  ! y=0 face
        yyy=0.00001
        do i=1,nh
           do j=1,nv
              xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
              zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
              call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
           enddo
        enddo

      case (4)  ! y=1 face
        yyy=0.99999
        do i=1,nh
           do j=1,nv
              xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
              zzz=b+(real(j)-.5d0)*2.d0**(-(level+1))
              call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
           enddo
        enddo

      case (5)  ! z=0 face
        zzz=0.00001
        do i=1,nh
           do j=1,nv
              xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
              yyy=b+(real(j)-.5d0)*2.d0**(-(level+1))
              call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
           enddo
        enddo
      case (6)  ! z=1 face
        zzz=0.99999
        do i=1,nh
           do j=1,nv
              xxx=l+(real(i)-.5d0)*2.d0**(-(level+1))
              yyy=b+(real(j)-.5d0)*2.d0**(-(level+1))
              call octree_create_from_particle(osolve%octree,osolve%noctree,xxx,yyy,zzz,1,level)
           enddo
        enddo
   end select

   enddo

   if (params%debug>=1 .and. iproc.eq.0) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
   if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after faces refinement ',osolve%octree(2)
   
end if


!-----------------------------------------------
! finally smooth the octree
!-----------------------------------------------
      
call show_time (total,step,inc,1,'Smoothen octree after refinement$')

call octree_smoothen (osolve%octree,osolve%noctree)

if (params%debug>=1 .and. iproc.eq.0) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after smoothing ',osolve%octree(2)

!-----------------------------------------------
! super smooth the octree if needed
!-----------------------------------------------
      
if (params%ismooth) then
   call show_time (total,step,inc,1,'Super smoothing of osolve$')

   call octree_super_smoothen (osolve%octree,osolve%noctree)

   if (params%debug>=1 .and. iproc.eq.0) write(*,'(a,i8,a)') shift//'osolve%nleaves=',osolve%octree(2)
   if (params%debug.gt.1) write(threadinfo%Logunit,*) '[same] improve_osolve: osolve%octree(2) after super smoothing ',osolve%octree(2)
   
end if

if (params%debug.gt.1) call flush(threadinfo%Logunit) 

end

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