https://github.com/florentrenaud/nbody6tt
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
fcloud.f
SUBROUTINE FCLOUD(I,FREG,FDR,KCASE)
*
*
* Force due to interstellar clouds.
* ---------------------------------
*
INCLUDE 'common6.h'
COMMON/CLOUDS/ XCL(3,MCL),XDOTCL(3,MCL),BODYCL(MCL),RCL2(MCL),
& CLM(MCL),CLMDOT(MCL),CLDOT,VCL,SIGMA,RB2,PCL2,
& TCL,STEPCL,NCL,NEWCL
REAL*8 FREG(3),FDR(3),A(6)
*
*
* Distinguish between standard case and initialization.
IF (KCASE.EQ.1) THEN
RB3 = RB2*SQRT(RB2)
*
* Sum over all clouds.
DO 10 ICL = 1,NCL
DO 5 K = 1,3
A(K) = XCL(K,ICL) - X(K,I)
A(K+3) = XDOTCL(K,ICL) - XDOT(K,I)
5 CONTINUE
*
RIJ2 = A(1)**2 + A(2)**2 + A(3)**2 + RCL2(ICL)
A7 = BODYCL(ICL)/RB3
A8 = BODYCL(ICL)/(RIJ2*SQRT(RIJ2))
A9 = 3.0*(A(1)*A(4) + A(2)*A(5) + A(3)*A(6))/RIJ2
*
* Add cloud force & derivative but subtract the average field effect.
DO 8 K = 1,3
FREG(K) = FREG(K) + A(K)*A8 + (X(K,I) - RDENS(K))*A7
FDR(K) = FDR(K) + (A(K+3) - A(K)*A9)*A8 + XDOT(K,I)*A7
8 CONTINUE
10 CONTINUE
ELSE
*
* Obtain total cloud force & first derivative (for routine FPOLY1).
DO 20 ICL = 1,NCL
DO 15 K = 1,3
A(K) = XCL(K,ICL) - X(K,I)
A(K+3) = XDOTCL(K,ICL) - XDOT(K,I)
15 CONTINUE
*
RIJ2 = A(1)**2 + A(2)**2 + A(3)**2 + RCL2(ICL)
A8 = BODYCL(ICL)/(RIJ2*SQRT(RIJ2))
A9 = 3.0*(A(1)*A(4) + A(2)*A(5) + A(3)*A(6))/RIJ2
*
* Include direct cloud contributions but ignore the average field.
DO 18 K = 1,3
FR(K,I) = FR(K,I) + A(K)*A8
D1R(K,I) = D1R(K,I) + (A(K+3) - A(K)*A9)*A8
18 CONTINUE
20 CONTINUE
END IF
*
RETURN
*
END