get_Morton_range_active Subroutine

  subroutine get_Morton_range_active
    use mod_forest
    use mod_global_parameters

    ! Cut the sfc based on weighted decision.  
    ! Oliver Porth, 02.02.2012

    !!!Here you choose the weithts:!!
    integer, parameter :: wa=3, wp=1
    ! wp : Weight for passive block
    ! wa : Weight for active block
    ! wp=0 : balance load (active blocks) exactly and 
    ! don't care about memory imbalance
    ! wp=wa : balance memory exactly and don't care about load
    ! Maximum possible memory imbalance is X=wa/wp.

    ! If you run into memory issues, decrease this ratio.
    ! Scaling should be better if you allow for higher ratio.  
    ! Note also that passive cells still do regridding and boundary swap.  
    ! I have best results with a ratio 2:1, but it is problem dependent.
    ! It can't do magic though...
    ! Best to make sure that the sfc is properly aligned with your problem. 

    integer :: ipe, Morton_no
    integer :: Mtot, Mstop, Mcurr
    ! For debugging: 
    integer :: nactive(0:npe-1),npassive(0:npe-1)
    !double precision, save :: ptasum=0

    if (allocated(sfc_to_igrid)) deallocate(sfc_to_igrid)
    

    Mtot  = nleafs_active*wa+(nleafs-nleafs_active)*wp
    ipe = 0 
    Mcurr = 0

    nactive=0
    npassive=0

    Morton_start(0) = 1
    do Morton_no=1,nleafs
       ! This is where we ideally would like to make the cuts:
       Mstop  = (ipe+1)*int(Mtot/npe)+min(ipe+1,mod(Mtot,npe))
       ! Build up mass:
       Mcurr = Mcurr + (wa*sfc(3,Morton_no)+wp*(1-sfc(3,Morton_no)))

       if (sfc(3,Morton_no)==1) then 
          nactive(ipe) = nactive(ipe) +1
          else
             npassive(ipe) = npassive(ipe) +1
          end if

       if (Mcurr >= Mstop) then 
          Morton_stop(ipe) = Morton_no
          ipe = ipe +1
          if (ipe>=npe) exit
          Morton_start(ipe) = Morton_no + 1
       end if
    end do

    Xmemory=dble(maxval(npassive+nactive))/dble(minval(npassive+nactive))
    Xload=dble(maxval(nactive))/dble(minval(nactive))

    !ptasum = ptasum +dble(nleafs-nleafs_active)/dble(nleafs_active)

    !if (mype == 0) print*, 'nleafs_passive:',nleafs-nleafs_active, 'nleafs_active:',nleafs_active,'ratio:',dble(nleafs-nleafs_active)/dble(nleafs_active),'mean ratio:',ptasum/it

    if (Morton_stop(mype)>=Morton_start(mype)) then
       allocate(sfc_to_igrid(Morton_start(mype):Morton_stop(mype)))
    
    end if

  end subroutine get_Morton_range_active