select_slice Subroutine

  subroutine select_slice(dir,xslice,writeonly,file_handle,normconv)
    use mod_forest, only: iglevel1_sfc, sfc_iglevel1, tree_node_ptr,&
        tree_root, Morton_sub_start, Morton_sub_stop
    use mod_global_parameters
    integer, intent(in) :: dir
    double precision, intent(in) :: xslice
    integer, intent(in) :: file_handle
    logical, intent(in) :: writeonly
    double precision,dimension(0:nw+nwauxio),intent(out)       :: normconv
    ! .. local ..
    integer :: ig1,ig2,ig3, jgrid, slice_fh, ipe, mylevmax
    integer, dimension(nlevelshi) :: igslice
    integer, allocatable :: iglevel1_sub_sfc(:)
    integer :: igmin1, igmin2, igmin3, igmax1, igmax2, igmax3, ngsub1,&
       jglevel1, igsorted, kgrid, ngmax

    kgrid = 0
    jgrid = 0
    mylevmax = 0
    jglevel1 = 0

    ! Find the global slice index for every level:
    call get_igslice(dir,xslice,igslice)
    igmin1=1; igmin2=1; igmin3=1;
    igmax1=ng1(1); igmax2=ng2(1); igmax3=ng3(1);

    ! Traverse forest to find grids indicating slice:

    select case(dir)
    case (1)
       igmin1=igslice(1); igmax1=igslice(1);
    case (2)
       igmin2=igslice(1); igmax2=igslice(1);
    case (3)
       igmin3=igslice(1); igmax3=igslice(1);
    case default
       call mpistop("slice direction not clear in select_slice")
    end select

    ngsub1=((igmax1-igmin1+1)*(igmax2-igmin2+1)*(igmax3-igmin3+1))
    allocate(iglevel1_sub_sfc(ngsub1))
    if (.not. allocated(sfc_sub)) then
       ngmax = ngsub1*(4**refine_max_level-1)
       allocate(sfc_sub(0:ngmax,1:3))
       sfc_sub = 0
    end if

    do ig3=igmin3,igmax3
       do ig2=igmin2,igmax2
          do ig1=igmin1,igmax1
             jglevel1=jglevel1+1
             iglevel1_sub_sfc(jglevel1)=iglevel1_sfc(ig1,ig2,ig3)
          end do
       end do
    end do

    call quicksort_integer(iglevel1_sub_sfc, 1, ngsub1)

    do igsorted=1,ngsub1
       ig1 = sfc_iglevel1(1,iglevel1_sub_sfc(igsorted))
       ig2 = sfc_iglevel1(2,iglevel1_sub_sfc(igsorted))
       ig3 = sfc_iglevel1(3,iglevel1_sub_sfc(igsorted))
       call traverse_slice(tree_root(ig1,ig2,ig3))
    end do
    deallocate(iglevel1_sub_sfc)




    if (.not.writeonly) then
       ! Synchronize the levmax_sub for output (only rank 0 needs it):
       levmax_sub = mylevmax
       call MPI_ALLREDUCE(MPI_IN_PLACE,levmax_sub,1,MPI_INTEGER,MPI_MAX,icomm,&
          ierrmpi)

       ! Communicate the subgrid indices according to new Morton sub-sfc:
       Morton_sub_start(:) = 1
       do ipe=0,npe-1
          call MPI_GATHER(jgrid,1,MPI_INTEGER,Morton_sub_stop,1,MPI_INTEGER,&
             ipe,icomm,ierrmpi)
       end do

       do ipe = 0, npe-2
          Morton_sub_start(ipe+1) = Morton_sub_stop(ipe)+Morton_sub_start(ipe+&
             1)
          Morton_sub_stop(ipe+1)  = Morton_sub_stop(ipe)+Morton_sub_stop(ipe+&
             1)
       end do
    end if

  contains

    recursive subroutine traverse_slice(tree)
      implicit none
      type(tree_node_ptr) :: tree
      integer :: ic1,ic2,ic3

      if (mype==0) then
         kgrid = kgrid + 1
         sfc_sub(kgrid,1) = tree%node%igrid
         sfc_sub(kgrid,2) = tree%node%ipe
         if (tree%node%leaf) then
            sfc_sub(0,1) = sfc_sub(0,1) + 1
            sfc_sub(kgrid,3) = 1
         else
            sfc_sub(0,2) = sfc_sub(0,2) + 1
            sfc_sub(kgrid,3) = 0
         end if
      end if

      if (tree%node%leaf) then
         if (tree%node%ipe == mype.and..not.writeonly) then
            mylevmax = max(mylevmax,tree%node%level)
            call fill_subnode(tree%node%igrid,tree%node%active,jgrid,dir,&
               xslice,normconv)
         end if
         return
      end if
      ! We are out for leaves now, continue for branches

      ! Get the correct child:
      select case (dir)
         case (1)
         ic1 = igslice(tree%node%level+1) - 2 * tree%node%ig1 + 2

         case (2)
         ic2 = igslice(tree%node%level+1) - 2 * tree%node%ig2 + 2

         case (3)
         ic3 = igslice(tree%node%level+1) - 2 * tree%node%ig3 + 2

      case default
         call mpistop("slice direction not clear in traverse_slice")
      end select

      ! Recursively descend into the correct child branch:

      select case(dir)
      case (1)
         do ic3=1,2
            do ic2=1,2
               call traverse_slice(tree%node%child(ic1,ic2,ic3))
            end do
         end do
      case (2)
         do ic3=1,2
            do ic1=1,2
               call traverse_slice(tree%node%child(ic1,ic2,ic3))
            end do
         end do
      case (3)
         do ic2=1,2
            do ic1=1,2
               call traverse_slice(tree%node%child(ic1,ic2,ic3))
            end do
         end do
      case default
         call mpistop("slice direction not clear in traverse_slice")
      end select




    end subroutine traverse_slice

  end subroutine select_slice