!----------------------------------------------------------------------- ! ! SUBROUTINE force_group(p,nterms_gr,iterms_gr,drdotdr_gr,,dx_gr,dy_gr,dz_gr,pmass_gr,acc_g,pquad_gr,option) ! ! !----------------------------------------------------------------------- ! ! ! Subroutine to compute the force F_far componete in a cell-group ! ! !----------------------------------------------------------------------- USE fly_h implicit none ! Declaration of local variables. ! ------------------------------- INTEGER(KIND=4) :: p,i, nterms_gr,nqterms,j,n,it,iqin INTEGER(KIND=4), DIMENSION (maxnterm) :: qindex, smindex INTEGER(KIND=4), DIMENSION (ndim) ::i_mesh,ip,csign INTEGER :: il_sh REAL(KIND=8) :: drdeldrg,phsm,drsm,accsm,mass_s REAL(KIND=8) :: sdrdotdr, epsp,epsi,rpinveff REAL(KIND=8) :: xx,yy,zz,wt1,wt2,wt3,wt4,wt5,wt6,wt7,wt8 REAL(KIND=8) :: wi1,wj1,wk1,wi2,wj2,wk2 REAL(KIND=8), DIMENSION (maxnterm) ::r3inveff,rinveff,qr5inv,phiquad REAL(KIND=8), DIMENSION (maxnterm) :: r2inveff,acci REAL(KIND=8), DIMENSION (1:2*ndim-1) ::quad_l REAL(KIND=8), DIMENSION (ndim) ::pos_l,r_abs,f INTEGER(KIND=4), INTENT(inout) :: iterms_gr(:) REAL(KIND=8), INTENT(inout) ::drdotdr_gr(:),dx_gr(:),dy_gr(:),dz_gr(:) REAL(KIND=8), INTENT(inout) :: pmass_gr(:) REAL(KIND=8), INTENT(inout) :: acc_g(:) REAL(KIND=8), INTENT(inout) :: pquad_gr(:,:) CHARACTER(LEN=4) :: option !----------------------------------------------------------------------- il_sh=p-nbodsmax pos_l(1:ndim)=pos_cell(1:ndim,il_sh) !----------------------------------------------------------------------- ! Compute monopole contribution; temporarily set mass of body p to ! zero to avoid possible self-interaction contribution. !----------------------------------------------------------------------- rpinveff=1./(0.+tiny) epsp=eps !----------------------------------------------------------------------- ! Loop over interaction list. !----------------------------------------------------------------------- DO 30 i=1,nterms_gr it=iterms_gr(i) sdrdotdr=SQRT(drdotdr_gr(i)) rinveff(i)=1./(sdrdotdr+tiny) r3inveff(i)=rinveff(i)/(drdotdr_gr(i)+tiny) epsi=eps drdeldrg=sdrdotdr*ninterp/(epsp+epsi) smindex(i)=drdeldrg smindex(i)=MIN(ninterp,smindex(i)) drsm=MIN(one,drdeldrg-smindex(i)) phsm=(1.-drsm)*phsmooth(smindex(i))+ & drsm*phsmooth(1+smindex(i)) accsm=(1.-drsm)*acsmooth(smindex(i))+ & drsm*acsmooth(1+smindex(i)) rinveff(i)=phsm*rinveff(i) r3inveff(i)=accsm*r3inveff(i) acci(i)=pmass_gr(i)*r3inveff(i) !----------------------------------------------------------------------- ! Ewald boundary periodic conditions section !----------------------------------------------------------------------- xx=dx_gr(i) yy=dy_gr(i) zz=dz_gr(i) r_abs(1)=abs(xx) r_abs(2)=abs(yy) r_abs(3)=abs(zz) if(r_abs(1).eq.0.0) then csign(1)=0 else csign(1)=xx/r_abs(1) endif if(r_abs(2).eq.0.0) then csign(2)=0. else csign(2)=yy/r_abs(2) endif if(r_abs(3).eq.0.0) then csign(3)=0. else csign(3)=zz/r_abs(3) endif i_mesh=r_abs*linv ip=i_mesh+1 wi1=1.-linv*(r_abs(1)-rk(i_mesh(1),1)) wi2=1.-linv*(rk(ip(1),1)-r_abs(1)) wj1=1.-linv*(r_abs(2)-rk(i_mesh(2),2)) wj2=1.-linv*(rk(ip(2),2)-r_abs(2)) wk1=1.-linv*(r_abs(3)-rk(i_mesh(3),3)) wk2=1.-linv*(rk(ip(3),3)-r_abs(3)) wt1=wi1*wj1*wk1 wt2=wi2*wj1*wk1 wt3=wi1*wj2*wk1 wt4=wi1*wj1*wk2 wt5=wi2*wj2*wk1 wt6=wi2*wj1*wk2 wt7=wi1*wj2*wk2 wt8=wi2*wj2*wk2 f(1:3)=wt1*fc(i_mesh(1),i_mesh(2),i_mesh(3),1:3) f(1:3)=f(1:3)+wt2*fc(ip(1),i_mesh(2),i_mesh(3),1:3) f(1:3)=f(1:3)+wt3*fc(i_mesh(1),ip(2),i_mesh(3),1:3) f(1:3)=f(1:3)+wt4*fc(i_mesh(1),i_mesh(2),ip(3),1:3) f(1:3)=f(1:3)+wt5*fc(ip(1),ip(2),i_mesh(3),1:3) f(1:3)=f(1:3)+wt6*fc(ip(1),i_mesh(2),ip(3),1:3) f(1:3)=f(1:3)+wt7*fc(i_mesh(1),ip(2),ip(3),1:3) f(1:3)=f(1:3)+wt8*fc(ip(1),ip(2),ip(3),1:3) f=csign*f*pmass_gr(i) acc_g(1)=acc_g(1)+f(1) acc_g(2)=acc_g(2)+f(2) acc_g(3)=acc_g(3)+f(3) !----------------------------------------------------------------------- ! END Ewald section !----------------------------------------------------------------------- 30 CONTINUE IF(option.NE.'pot ') THEN DO 50 i=1,nterms_gr acc_g(1)=acc_g(1)-dx_gr(i)*acci(i) acc_g(2)=acc_g(2)-dy_gr(i)*acci(i) acc_g(3)=acc_g(3)-dz_gr(i)*acci(i) 50 CONTINUE ENDIF !----------------------------------------------------------------------- ! If required, compute quadrupole contribution. !----------------------------------------------------------------------- IF(usequad) THEN !----------------------------------------------------------------------- ! Filter out bodies. !----------------------------------------------------------------------- nqterms=0 DO i=1,nterms_gr IF(iterms_gr(i).GT.nbodsmax) THEN nqterms=nqterms+1 qindex(nqterms)=i ENDIF ENDDO IF(option.NE.'pot ') THEN DO i=1,nqterms iqin=qindex(i) quad_l(1:2*ndim-1)=pquad_gr(1:2*ndim-1,iqin) r2inveff(i)=rinveff(iqin)*rinveff(iqin) qr5inv(i)=r3inveff(iqin)*r2inveff(i) phiquad(i)=((-.5*((dx_gr(iqin)**2-dz_gr(iqin)**2)* & quad_l(1)+(dy_gr(iqin)**2- & dz_gr(iqin)**2)*quad_l(4))- & (dx_gr(iqin)*dy_gr(iqin)*quad_l(2)+ & dx_gr(iqin)*dz_gr(iqin)*quad_l(3)+ & dy_gr(iqin)*dz_gr(iqin)*quad_l(5)))* & qr5inv(i))*5.*r2inveff(i) acc_g(1)=acc_g(1)+dx_gr(iqin)*phiquad(i)+ & (dx_gr(iqin)*quad_l(1)+ & dy_gr(iqin)*quad_l(2)+ & dz_gr(iqin)*quad_l(3))*qr5inv(i) acc_g(2)=acc_g(2)+dy_gr(iqin)*phiquad(i)+ & (dy_gr(iqin)*quad_l(4)+ & dx_gr(iqin)*quad_l(2)+ & dz_gr(iqin)*quad_l(5))*qr5inv(i) acc_g(3)=acc_g(3)+dz_gr(iqin)*phiquad(i)+ & (dz_gr(iqin)*(-quad_l(1)- & quad_l(4))+dx_gr(iqin)* & quad_l(3)+dy_gr(iqin)* & quad_l(5))*qr5inv(i) ENDDO ENDIF !IF(option.NE.'pot ') ENDIF ! IF(usequad) RETURN END