File: N:\mfix\model\GhdTheory\mass_mobility.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2     !
3     !  subroutine name: call mass_mobility(s,mi,ni,rho,zeta0,theta,nu,DF)
4     !
5     !  author:  C. Hrenya, Feb 2009
6     !
7     !  Purpose: find mass mobility coefficient according to GHD polydisperse KT
8     !
9     !  Literature/References:
10     !     C. Hrenya handwritten notes & Garzo, Hrenya, Dufty papers (PRE, 2007)
11     !
12     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
13     
14           subroutine mass_mobility(s,mi,ni,rho,zeta0,theta,nu,DF)
15           Implicit NONE
16     
17           integer s, indx(s)
18     
19           double precision mi(s),ni(s),rho, &
20                           zeta0,theta(s),nu(s,s),DF(s,s)
21     
22           integer i,j,kk
23           double precision kronecker(s,s),Amat(s,s),bmat(s,s), &
24                           Amat0(s,s),bmat0(s)
25           double precision d
26     
27           integer NP
28           parameter (NP=15)     !max no. of linear equations to solve
29     
30           do i=1,s
31              do j=1,s
32                 if (i.eq.j) then
33                    kronecker(i,j) = 1.d0
34                 else
35                    kronecker(i,j) = 0.d0
36                 endif
37              enddo
38           enddo
39     
40           do i=1,s
41              do j=1,s
42                  Amat(i,j) = (nu(i,j)+0.5d0*zeta0*kronecker(i,j))    !A matrix for solution of DF (p 8 CMH notes)
43                  bmat(i,j) = -ni(i)*mi(i)/mi(j)*(kronecker(i,j)-ni(j) &
44                     *mi(j)/rho)                                      !b matrix for solution of DF (p 8 CMH notes)
45              enddo
46           enddo
47     
48     ! this extra kk loop and addition of Amat0 and bmat0 is necessary
49     ! since x & b in Ax=b are s by s matrices rather than vectors of length s,
50     ! whereas LUBSKB is specific to x & b vectors of length s
51           do kk=1,s
52              do i=1,s
53                 do j=1,s
54                     Amat0(i,j) = Amat(i,j)
55                 enddo
56              enddo
57              do i=1,s
58                 bmat0(i) = bmat(i,kk)
59              enddo
60     
61              CALL LUDCMP(Amat0, s, NP, indx, d, 'mass_mobility') ! solve system of s linear equations using
62              CALL LUBKSB(Amat0, s, NP, indx, bmat0) ! LU decomposition
63     
64              do i=1,s
65                 DF(i,kk) = bmat0(i)
66              enddo
67     
68           enddo
69     
70           return
71           end subroutine mass_mobility
72