fcloud.f

Go to the documentation of this file.

      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