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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
2     !
3     !  subroutine name: thermal_mobility(s,mi,alpha,ni,mu,sigma,chi,zeta0,
4     !                                theta,Ti,DF,gammaij,omega,Lij)
5     !
6     !  author:  C. Hrenya, Feb 2009
7     !
8     !  Purpose: find thermal mobility coefficient according to GHD polydisperse KT
9     !
10     !  Literature/References:
11     !     C. Hrenya handwritten notes & Garzo, Hrenya, Dufty papers (PRE, 2007)
12     !
13     !  Modifications:
14     !     Sof: removed /ni in sum1, which yield zero based on Bill's mods.
15     !
16     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
17     
18           subroutine thermal_mobility(s,mi,alpha,ni,mu,sigma,chi,zeta0, &
19                                      theta,Ti,DF,gammaij,omega,Lij)
20           Implicit NONE
21     
22           integer s, indx(s)
23     
24           double precision mi(s),alpha(s,s),ni(s),mu(s,s),sigma(s,s), &
25                           chi(s,s),zeta0,theta(s),Ti(s),DF(s,s), &
26                           gammaij(s,s),omega(s,s),Lij(s,s)
27     
28           integer i,j,k,p,kk
29           double precision kronecker(s,s),sum1(s,s),l_bar(s,s),lkj(s,s), &
30                           Lkin(s,s),Lcol(s,s),Amat(s,s),bmat(s,s), &
31                           Amat0(s,s),bmat0(s)
32           double precision d
33     
34           integer NP
35           parameter (NP=15)     !max no. of linear equations to solve
36     
37           double precision pi
38           parameter (pi=3.14159265458979323846d0)
39     
40           do i=1,s
41              do j=1,s
42                 if (i.eq.j) then
43                    kronecker(i,j) = 1.d0
44                 else
45                    kronecker(i,j) = 0.d0
46                 endif
47              enddo
48           enddo
49     
50     !calculate summation used in l_bar (b vector) - p. 19 of CMH notes
51           do i=1,s
52              do j=1,s
53                 sum1(i,j) = 0.d0
54              enddo
55           enddo
56           do i=1,s
57              do j=1,s
58                 do k=1,s
59     ! in case ni(k) == 0d0, the contribution to sum1 can be neglected based on earlier code by Bill Holloway.
60                    if(ni(k) > 0d0) sum1(i,j) = sum1(i,j) + (omega(i,k)-zeta0* &
61                                             kronecker(i,k))/(ni(k)*Ti(k))*DF(k,j)
62                 enddo
63              enddo
64           enddo
65     
66     !calculate l_bar (p 19 of CMH notes)
67           do i=1,s
68              do j=1,s
69                 l_bar(i,j) = -2.5d0*ni(i)*Ti(i)**2/mi(i)*sum1(i,j)
70                 Amat(i,j) = gammaij(i,j) - 0.5d0*zeta0*kronecker(i,j)     !A matrix for solution of lkj (p 19 CMH notes)
71                 bmat(i,j) = l_bar(i,j)                                    !B matrix for solution of lkj (p 19 CMH notes)
72              enddo
73           enddo
74     
75     ! this extra kk loop and addition of Amat0 and bmat0 is necessary
76     ! since x & b in Ax=b are s by s matrices rather than vectors of length s,
77     ! whereas LUBSKB is specific to x & b vectors of length s
78           do kk=1,s
79              do i=1,s
80                 do j=1,s
81                     Amat0(i,j) = Amat(i,j)
82                 enddo
83              enddo
84              do i=1,s
85                 bmat0(i) = bmat(i,kk)
86              enddo
87     
88              CALL LUDCMP(Amat0, s, NP, indx, d, 'thermal_mobility') ! solve system of s linear equations using
89              CALL LUBKSB(Amat0, s, NP, indx, bmat0) ! LU decomposition
90     
91              do i=1,s
92                 lkj(i,kk) = bmat0(i)
93              enddo
94           enddo
95     
96     !kinetic contribution to thermal mobility (p 19 CMH notes)
97           do i = 1,s
98              do j=1,s
99                 Lkin(i,j) = lkj(i,j) + 2.5d0*Ti(i)/mi(i)*DF(i,j)
100              enddo
101           enddo
102     
103     !collisional contribution to thermal mobility (p 20 CMH notes)
104           do i=1,s
105              do j=1,s
106                 Lcol(i,j) = 0.d0
107              enddo
108           enddo
109           do i=1,s
110              do j=1,s
111                 do p=1,s
112                    Lcol(i,j) = Lcol(i,j) + (1.d0+alpha(i,p))/8.d0*mi(p)* &
113                      mu(i,p)*sigma(i,p)**3*chi(i,p)*(4.d0*pi/5.d0* &
114                      (1.d0-alpha(i,p))*(mu(i,p)-mu(p,i))*ni(i)  &
115                      *(2.d0/mi(p)*Lkin(p,j)+5.d0*Ti(i)/mi(i)/mi(p)* &
116                      DF(p,j))+48.d0*pi/15.d0*ni(i)*(2.d0*mu(p,i)/mi(p)* &
117                      Lkin(p,j)-5.d0*(2.d0*mu(i,p)-mu(p,i))*Ti(i)/mi(i)/ &
118                      mi(p)*DF(p,j)))
119                 enddo
120                 Lij(i,j) = Lkin(i,j) + Lcol(i,j)          !thermal mobility (p 19 CMH notes)
121              enddo
122           enddo
123     
124           return
125           end subroutine thermal_mobility
126