!***********************************************************************
!
!
                          SUBROUTINE cell_prop
!
!
!***********************************************************************
!
!
!     Subroutine to compute masses, center of mass coordinates,
!     and optional quadrupole moments of cells, processing cells
!     in order of increasing size.  
!
!
!=======================================================================

	 USE fly_h
	 implicit none
	INCLUDE 'mpif.h'

!   Declaration of local variables.
!-----------------------------------
	
        INTEGER		:: i,j,k,l,m,n,ind_l1,ind1_th
        INTEGER(KIND=4), DIMENSION(:), ALLOCATABLE	:: isub 
        INTEGER		:: subpv,i_l
        REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE	:: pl
	REAL(KIND=8), DIMENSION(:), ALLOCATABLE	::pcl
	REAL(KIND=8), DIMENSION(:), ALLOCATABLE	::mc

	REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE:: qadd1
			INTEGER:: status

!-----------------------------------------------------------------------
!   Generate permutation of cells, according to cellsize.
!
!   Initialize properties of cells.
!-----------------------------------------------------------------------

           mass_cell=0.
           pos_cell=0.

        IF(usequad) THEN
	
	   quad=0.

        ENDIF

!----------------------------------------------------------------------- 
!   Process cells in order of increasing size.
!-----------------------------------------------------------------------
	
 	DO 40 i_l=lmax,1,-1

!-----------------------------------------------------------------------
!   Compute properties of the selected cells, looping over subcells locally residing.
!-----------------------------------------------------------------------
!
!loop only on the local residing cells:
!check my  starting address on the cell list
           
!$OMP PARALLEL PRIVATE(ind1_th,isub,subpv,ind_l1,mc,pl,pcl,qadd1,i,l,j,m,n,k) 
     ALLOCATE(isub(nsubcell), STAT=status) 
     ALLOCATE(mc(nsubcell), STAT=status) 
     ALLOCATE(pl(ndim,nsubcell), STAT=status) 
     ALLOCATE(pcl(ndim), STAT=status) 
     ALLOCATE(qadd1(2*ndim-1,nsubcell), STAT=status) 
      
!$OMP DO
	   DO  i=cell_ss(i_l,1),cell_ss(i_l,2)
	       ind1_th=i-nb_loc     
           isub(1:8)=subp(1:8,ind1_th)
		
             DO j=1,nsubcell  
                subpv=isub(j)
				IF(subpv.gt.nb_loc) THEN
					ind_l1=subpv-nb_loc
					mc(j)=mass_cell(ind_l1)
					pl(1:ndim,j)=pos_cell(1:ndim,ind_l1)
					IF(usequad) qadd1(1:5,j)=quad(1:5,ind_l1)
        		ENDIF
        		IF(subpv.GT.0.and.subpv.le.nb_loc) THEN
                   mc(j) = mass_read
		           qadd1(1:5,j)=0.
					pl(1:ndim,j)=pos(1:ndim,subpv)
        		ENDIF !if (subpv.GT.0.and.subpv.le.nb_loc)
              ENDDO ! do j
              
	      	  DO j=1,nsubcell
                   IF (isub(j) .gt. 0) THEN
                     mass_cell(ind1_th)=mass_cell(ind1_th)+mc(j)
                     pos_cell(1,ind1_th)=pos_cell(1,ind1_th)+mc(j)*pl(1,j)
                     pos_cell(2,ind1_th)=pos_cell(2,ind1_th)+mc(j)*pl(2,j)
                     pos_cell(3,ind1_th)=pos_cell(3,ind1_th)+mc(j)*pl(3,j)                 
                   ENDIF
             ENDDO ! do j

           pos_cell(1,ind1_th)=pos_cell(1,ind1_th)/mass_cell(ind1_th)
           pos_cell(2,ind1_th)=pos_cell(2,ind1_th)/mass_cell(ind1_th)
           pos_cell(3,ind1_th)=pos_cell(3,ind1_th)/mass_cell(ind1_th)

!-----------------------------------------------------------------------
!  Compute optional quadrupole moments.
!-----------------------------------------------------------------------

           IF(usequad .AND. (mass_cell(ind1_th).GT.mass_read)) THEN

             DO j =1,nsubcell
 				pcl(1:ndim)=pos_cell(1:ndim,ind1_th)            
                If (subp(j,ind1_th) .eq. 0) CYCLE
                DO 200 m=1,2
                   DO 190 n = m,ndim
                     l = (m-1)*(ndim-1)+n
                     quad(l,ind1_th) = quad(l,ind1_th)+                 &
                      mc(j)*(3.*(pl(m,j)-pcl(m))*(pl(n,j)-pcl(n)))+      &
                       qadd1(l,j)
                     IF ( m .eq. n) THEN
                        DO k=1,ndim
                           quad(l,ind1_th) = quad(l,ind1_th) -  & 
                           		mc(j)*(pl(k,j)-pcl(k))**2
                        ENDDO
                     ENDIF   
190                continue
200              continue
               ENDDO ! do j=...
           
	   		ENDIF   !IF usequad
	   
	   ENDDO  ! do i=
!$OMP END DO
!----------------------------------------------------------------------- 
!   Start new lower level
!-----------------------------------------------------------------------
     DEALLOCATE(isub) 
     DEALLOCATE(mc) 
     DEALLOCATE(pl) 
     DEALLOCATE(pcl) 
     DEALLOCATE(qadd1) 

!$OMP END PARALLEL  !**QUI sarebbe meglio mettere una barrier e non deallocare e fare la sessione parallela includendo anche il loop 40 CONTINUE
	
 40 	CONTINUE

!monitor check

	IF(me.EQ.0) THEN
	write(uterm,*) 'CELL_PROP: mass root', mass_cell(1),                   &
        mass_read*nbodies, 'pos. root', pos_cell(1:3,1),                 &
        'quad. root', quad(1:5,1)
	call flush(uterm)
	
	ENDIF
        
	RETURN
        END