Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
!***********************************************************************
!
!
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