!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              COMPUTE_CONVERGENCE_CRITERION    Apr. 2007                      |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

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

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! subroutine that computes a convergence criterion based on the difference 
! between the velocity field obtained at this iteration (osolve) and the 
! previous velocity field (ov)

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

use definitions
   
implicit none

type (octreesolve) osolve
type (octreev) ov
type (void) vo
type (parameters) params
integer istep,iter, iter_nl
double precision maxu,maxv,maxw
integer refine_level
logical velocity_converged

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|
  
integer err,ierr,iproc,nproc,i
double precision dduvw,duvw,uvw,time1,time2
double precision velocity_diff_norm,maxdiff,oldmax
character*72 shift

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

INCLUDE 'mpif.h'
call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

shift=' '

duvw=0.d0
uvw=0.d0
maxu=0.d0
maxv=0.d0
maxw=0.d0
oldmax= 0.d0
maxdiff=0.d0

do i=1,osolve%nnode
   if (vo%influid(i)) then
      dduvw=(osolve%u(i)-ov%unode(i))**2&
           +(osolve%v(i)-ov%vnode(i))**2&
           +(osolve%w(i)-ov%wnode(i))**2
      !duvw=duvw+dduvw
      !uvw=uvw+(ov%unode(i))**2+(ov%vnode(i))**2+(ov%wnode(i))**2
      uvw=(ov%unode(i))**2+(ov%vnode(i))**2+(ov%wnode(i))**2
      !maxu=max(maxu,abs(osolve%u(i)-ov%unode(i)))
      !maxv=max(maxv,abs(osolve%v(i)-ov%vnode(i)))
      !maxw=max(maxw,abs(osolve%w(i)-ov%wnode(i)))
      oldmax=max(oldmax,uvw)
      maxdiff=max(maxdiff,dduvw)
   endif
enddo

! Uncommented for testing - dwhipp 12/10
!velocity_diff_norm=sqrt(duvw)/sqrt(uvw)

! New convergence criterion - dwhipp 12/10
velocity_diff_norm=sqrt(maxdiff)/sqrt(oldmax)

! Commented out for testing - dwhipp 12/10
!velocity_diff_norm=max(maxu,maxv,maxw)/2.d0

!velocity_converged=(velocity_diff_norm<params%tol)

velocity_converged=(velocity_diff_norm<params%tol .and. iter_nl > params%nb_iter_nl_min) 

!velocity_converged=(velocity_diff_norm<params%tol &
!                    .and. maxu<params%tol &
!                    .and. maxv<params%tol &
!                    .and. maxw<params%tol &
!                    .and. iter_nl > 5) 


if (iproc.eq.0) then
   write(8,*) 'in compute_convergence_criterion, velocity_diff_norm=', velocity_diff_norm
   write(8,*) 'minval(osolve%u)=',minval(osolve%u),'maxval(osolve%u)',maxval(osolve%u)
   write(8,*) 'minval(osolve%v)=',minval(osolve%v),'maxval(osolve%v)',maxval(osolve%v)
   write(8,*) 'minval(osolve%w)=',minval(osolve%w),'maxval(osolve%w)',maxval(osolve%w)

   if (params%debug>=1) then
      write(*,'(a,f10.7,a,f10.7,a)') shift//'velocity_diff_norm =',velocity_diff_norm,' (tol =',params%tol,')'
      if (velocity_converged) then
         write(*,'(a)')             shift//' '
         write(*,'(a)')             shift//'convergence on this octree reached'
      end if
   end if
end if   

call DoRuRe_conv_stats (params%doDoRuRe,istep,iter,iter_nl,velocity_diff_norm,params%tol,maxu,maxv,maxw)

if (params%adaptive_tol .and. velocity_converged) params%tol=velocity_diff_norm

return
end

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