File: /nfs/home/0/users/jenkins/mfix.git/model/GhdTheory/cooling_rate_tc.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16 subroutine cooling_rate_tc(s,mi,sigmai,alpha,ni,n,v0,mu,sigma, &
17 chi,T,zeta0,theta,Ti,zetau)
18 Implicit NONE
19
20 integer s, indx(s)
21
22 double precision mi(s),sigmai(s),alpha(s,s),ni(s),n,v0,mu(s,s), &
23 sigma(s,s),chi(s,s),T,zeta0,theta(s), &
24 Ti(s),zetau
25
26 integer i,j
27 double precision kronecker(s,s),sum1(s),psi(s,s),eid_bar(s), &
28 Amat(s,s),bmat(s),ejd(s),zeta10,zeta11
29 double precision d
30
31 integer NP
32 parameter (NP=15)
33
34 double precision pi
35 parameter (pi=3.14159265458979323846d0)
36
37 do i=1,s
38 do j=1,s
39 if (i.eq.j) then
40 kronecker(i,j) = 1.d0
41 else
42 kronecker(i,j) = 0.d0
43 endif
44 enddo
45 enddo
46
47 do i=1,s
48 sum1(i) = 0.d0
49 enddo
50 do i=1,s
51 do j=1,s
52 if (i.ne.j) then
53 sum1(i) = sum1(i)+dsqrt(pi)/2.d0*ni(j)*chi(i,j)* &
54 sigma(i,j)**2*v0*mu(j,i)*(1.d0+alpha(i,j))/ &
55 dsqrt(theta(i)*theta(j))/(theta(i)+theta(j))**2.5d0 &
56 *(-2.d0*(90.d0*theta(j)**3+231.d0*theta(i)* &
57 theta(j)**2+184.d0*theta(i)**2*theta(j)+40.d0* &
58 theta(i)**3)+3.d0*mu(j,i)*(1.d0+alpha(i,j))* &
59 (theta(i)+theta(j))*(70.d0*theta(j)**2+117.d0* &
60 theta(i)*theta(j)+44.d0*theta(i)**2)-24.d0* &
61 mu(j,i)**2*(1.d0+alpha(i,j))**2* &
62 (theta(i)+theta(j))**2*(5.d0*theta(j)+4.d0*theta(i)) &
63 +30.d0*mu(j,i)**3*(1.d0+alpha(i,j))**3*(theta(i)+ &
64 theta(j))**3+5.d0*theta(j)*(theta(i)+theta(j))* &
65 (2.d0*(4.d0*theta(i)+3.d0*theta(j))-3.d0*mu(j,i)* &
66 (1.d0+alpha(i,j))*(theta(i)+theta(j))))
67 endif
68 enddo
69 enddo
70
71
72 do i=1,s
73 eid_bar(i) = 0.d0
74 enddo
75 do i=1,s
76 do j=1,s
77 if (i.eq.j) then
78 psi(i,i) = -2.d0/15.d0*(sum1(i)+dsqrt(2.d0*pi) &
79 /32.d0*ni(i)*chi(i,i)*sigmai(i)**2*(1.d0+ &
80 alpha(i,i))/dsqrt(theta(i))*v0*(30.d0*alpha(i,i)**3 &
81 -126.d0*alpha(i,i)**2+177.d0*alpha(i,i)+48.d0* &
82 (3.d0*alpha(i,i)-7.d0)-137.d0)+3.d0/32.d0* &
83 dsqrt(2.d0*pi)*ni(i)*chi(i,i)*sigmai(i)**2* &
84 (1.d0+alpha(i,i))/dsqrt(theta(i))*v0 &
85 *(10.d0*alpha(i,i)**3+22.d0*alpha(i,i)**2+ &
86 11.d0*alpha(i,i)-3.d0))
87 else
88 psi(i,j) = -2.d0/15.d0*(dsqrt(pi)/2.d0* &
89 ni(j)*chi(i,j)*sigma(i,j)**2*v0*mu(j,i)*(1.d0+ &
90 alpha(i,j))*(theta(i)/theta(j))**1.5d0/(theta(i)+ &
91 theta(j))**2.5d0*((2.d0*theta(j)+5.d0*theta(i))* &
92 (2.d0*theta(j)+3.d0*mu(j,i)*(1.d0+alpha(i,j)) &
93 *(theta(i)+theta(j)))-24.d0*mu(j,i)**2*(1.d0+ &
94 alpha(i,j))**2*(theta(i)+theta(j))**2+30.d0* &
95 mu(j,i)**3*(1.d0+alpha(i,j))**3*(theta(i)+ &
96 theta(j))**3/theta(j)-5.d0*(theta(i)+theta(j))* &
97 (2.d0*theta(j)+3.d0*mu(j,i)*(1.d0+alpha(i,j))* &
98 (theta(i)+theta(j)))))
99 endif
100 eid_bar(i) = eid_bar(i) - pi/6.d0*ni(j)*chi(i,j)* &
101 sigma(i,j)**3*mu(j,i)*(1.d0+alpha(i,j)) &
102 *(40.d0*(mu(i,j)-1.d0)+4.d0*(19.d0+9.d0* &
103 alpha(i,j))*mu(j,i)-48.d0*mu(j,i)**2*(theta(i)+ &
104 theta(j))/theta(j)*(1.d0+alpha(i,j))**2 &
105 +15.d0*mu(j,i)**3*(theta(i)+theta(j))**2/ &
106 theta(j)**2*(1.d0+alpha(i,j))**3)
107 enddo
108 enddo
109
110 do i=1,s
111 do j=1,s
112 Amat(i,j) = psi(i,j) - 1.5d0*zeta0*kronecker(i,j)
113 enddo
114 bmat(i) = 2.d0*eid_bar(i)/15.d0
115 enddo
116
117 CALL LUDCMP(Amat, s, NP, indx, d, 'cooling_rate_tc')
118 CALL LUBKSB(Amat, s, NP, indx, bmat)
119
120 do i=1,s
121 ejd(i) = bmat(i)
122 enddo
123
124
125 = 0.d0
126 zeta11 = 0.d0
127 do i=1,s
128 do j=1,s
129 zeta10 = zeta10 - 2.d0*pi/(3.d0*n*T)*ni(i)*ni(j)*mu(j,i) &
130 *sigma(i,j)**3*chi(i,j)*(1.d0-alpha(i,j)**2)*Ti(i)
131 zeta11 = zeta11 + dsqrt(pi)*v0**3/(2.d0*n*T)*ejd(i)* &
132 ni(i)*ni(j)*sigma(i,j)**2*chi(i,j)*mi(j) &
133 *mu(i,j)*(1.d0-alpha(i,j)**2)*dsqrt(theta(j)/ &
134 (theta(i)**3*(theta(i)+theta(j))))
135 enddo
136 enddo
137
138 zetau = zeta10 + zeta11
139
140 return
141 end subroutine cooling_rate_tc
142
143