!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!          8888888b.   .d88888b.  888     888       d8888 8888888b.            |
!          888  "Y88b d88P" "Y88b 888     888      d88888 888   Y88b           |
!          888    888 888     888 888     888     d88P888 888    888           |
!          888    888 888     888 888     888    d88P 888 888   d88P           |
!          888    888 888     888 888     888   d88P  888 8888888P"            |
!          888    888 888     888 888     888  d88P   888 888 T88b             |
!          888  .d88P Y88b. .d88P Y88b. .d88P d8888888888 888  T88b            |
!          8888888P"   "Y88888P"   "Y88888P" d88P     888 888   T88b           |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              DEFINE_BC_NEST    May 2011                                      |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine define_bc_nest (nnode,kfix,kfixt,x,y,z,u,v,w,temp,vo,nest)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! this routine assigns the boundary condition for the boundaries of the nest in
! nested DOUAR simulations

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

use definitions

implicit none

integer :: nnode
integer :: kfix(nnode*3)
integer :: kfixt(nnode)
double precision :: x(nnode),y(nnode),z(nnode)
double precision :: u(nnode),v(nnode),w(nnode)
double precision :: temp(nnode)
type (void) vo
type (nest_info) nest

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

type (octreev) ovls
integer :: i,nface,nlsf,np,kfixls
logical :: is_plas
double precision :: eps,xls,yls,zls,s,crit,e2d,evisc,dilatr

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

eps=1.d-10

ovls%noctree=params%noctreemax
! Read in coarse model octree; needs to be stored in params @ end of step
open (9,file=trim(params%lsoutfile),status='old',form='unformatted')
  read (9) ovls%octree(3),&
           ovls%nnode,    &
           ovls%nleaves,  &
           nface,         &
           nlsf,          &
           np,            &
           current_time
  allocate (ovls%x(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%x in define_bc_nest$')
  allocate (ovls%y(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%y in define_bc_nest$')
  allocate (ovls%z(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%z in define_bc_nest$')
  allocate (ovls%icon(8,ovls%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%icon in define_bc_nest$')
  allocate (ovls%unode(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%unode in define_bc_nest$')
  allocate (ovls%vnode(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%vnode in define_bc_nest$')
  allocate (ovls%wnode(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%wnode in define_bc_nest$')
  allocate (ovls%wnodeiso(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%wnodeiso in define_bc_nest$')
  allocate (ovls%temp(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%temp in define_bc_nest$')
  allocate (ovls%pressure(ovls%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%pressure in define_bc_nest$') 
  allocate (ovls%spressure(ovls%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%spressure in define_bc_nest$')
  allocate (ovls%temporary_nodal_pressure(ovls%nnode),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%temp_nodal_pressure in define_bc_nest$')
  allocate (ovls%whole_leaf_in_fluid(ovls%nleaves),stat=err) ; if (err.ne.0) call stop_run ('Error alloc ovls%whole_leaf_in_fluid in define_bc_nest$')
  read (9) (ovls%x(i),                      &
           ovls%y(i),                       &
           ovls%z(i),                       &
           ovls%unode(i),                   &
           ovls%vnode(i),                   &
           ovls%wnode(i),                   &
           ovls%wnodeiso(i),                &
           (xlsf,j=1,nlsf),                 &
           ovls%temp(i),                    &
           ovls%temporary_nodal_pressure(i),&
           s,                               &
           kfixls,                          &
           kfixls,                          &
           kfixls,                          &
           kfixls,                          &
           i=1,ovls%nnode)
  read (9) ((ovls%icon(k,i),k=1,8),    &
           ovls%pressure(i),           & 
           ovls%spressure(i),          &
           crit,                       &
           e2d,                        &
           evisc,                      & 
           is_plas,                    &
           dilatr,                     &
           ovls%whole_leaf_in_fluid(i),&
           i=1,ovls%nleaves)
  read (9) (ovls%octree(i),i=1,ovls%octree(3))
close (9)

! Loop through all nodes in nest and interpolate B/C velocities from LS model
do i=1,nnode
  xls=x(i)*nest%sselemx+nest%xminls
  yls=y(i)*nest%sselemy+nest%yminls
  zls=z(i)*nest%sselemz+nest%zminls
  if (x(i).lt.eps) then
    kfix((i-1)*3+1)=1
    kfix((i-1)*3+2)=1
    kfix((i-1)*3+3)=1
    call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon,       &
                                  ovls%nleaves,ovls%nnode,xls,yls,zls,         &
                                  ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode,  &
                                  w(i))
  elseif (x(i).gt.1.d0-eps) then
    kfix((i-1)*3+1)=1
    kfix((i-1)*3+2)=1
    kfix((i-1)*3+3)=1
    call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon,       &
                                  ovls%nleaves,ovls%nnode,xls,yls,zls,         &
                                  ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode,  &
                                  w(i))
  elseif (y(i).lt.eps) then
    kfix((i-1)*3+1)=1
    kfix((i-1)*3+2)=1
    kfix((i-1)*3+3)=1
    call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon,       &
                                  ovls%nleaves,ovls%nnode,xls,yls,zls,         &
                                  ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode,  &
                                  w(i))
  elseif (y(i).gt.1.d0-eps) then
    kfix((i-1)*3+1)=1
    kfix((i-1)*3+2)=1
    kfix((i-1)*3+3)=1
    call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon,       &
                                  ovls%nleaves,ovls%nnode,xls,yls,zls,         &
                                  ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode,  &
                                  w(i))
  elseif (z(i).lt.eps) then
    kfix((i-1)*3+1)=1
    kfix((i-1)*3+2)=1
    kfix((i-1)*3+3)=1
    call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon,       &
                                  ovls%nleaves,ovls%nnode,xls,yls,zls,         &
                                  ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode,  &
                                  w(i))
    kfixt(i)=1
    temp(i)=1.d0
  elseif (z(i).gt.1.d0-eps) then
    kfix((i-1)*3+1)=1
    kfix((i-1)*3+2)=1
    kfix((i-1)*3+3)=1
    call octree_interpolate_three (3,ovls%octree,ovls%noctree,ovls%icon,       &
                                  ovls%nleaves,ovls%nnode,xls,yls,zls,         &
                                  ovls%unode,u(i),ovls%vnode,v(i),ovls%wnode,  &
                                  w(i))
    kfixt(i)=1
    temp(i)=0.d0
  endif
  if (.not.vo%influid(i)) then
    kfixt(i)=1
    temp(i)=0.d0
  endif
end do

end

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