MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
eos_polar.f90
1 MODULE eos_polar
2 
3  IMPLICIT NONE
4 
5  PRIVATE
6  PUBLIC :: f_polar, p_polar, phi_polar
7 
8 CONTAINS
9 
10 
11 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
12 !
13 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
14 
15 SUBROUTINE f_polar ( fdd, fqq, fdq )
16 
17  USE eos_variables, ONLY: ncomp, parame, dd_term, qq_term, dq_term
18 
19  !-----------------------------------------------------------------------------
20  REAL, INTENT(OUT) :: fdd, fqq, fdq
21 
22  !-----------------------------------------------------------------------------
23  INTEGER :: dipole
24  INTEGER :: quadrupole
25  INTEGER :: dipole_quad
26  !-----------------------------------------------------------------------------
27 
28  fdd = 0.0
29  fqq = 0.0
30  fdq = 0.0
31 
32  dipole = 0
33  quadrupole = 0
34  dipole_quad = 0
35  IF ( sum( parame(1:ncomp,6) ) /= 0.0 ) dipole = 1
36  IF ( sum( parame(1:ncomp,7) ) /= 0.0 ) quadrupole = 1
37  IF ( dipole == 1 .AND. quadrupole == 1 ) dipole_quad = 1
38 
39  !-----------------------------------------------------------------------------
40  ! dipole-dipole term
41  !-----------------------------------------------------------------------------
42  IF (dipole == 1) THEN
43 
44  IF (dd_term == 'GV') CALL f_dd_gross_vrabec( fdd )
45  ! IF (dd_term == 'SF') CALL F_DD_SAAGER_FISCHER( k )
46  ! IF (dd_term /= 'GV' .AND. dd_term /= 'SF') write (*,*) 'specify dipole term !'
47 
48  ENDIF
49 
50  !-----------------------------------------------------------------------------
51  ! quadrupole-quadrupole term
52  !-----------------------------------------------------------------------------
53  IF (quadrupole == 1) THEN
54 
55  IF (qq_term == 'JG') CALL f_qq_gross( fqq )
56  ! IF (qq_term == 'SF') CALL F_QQ_SAAGER_FISCHER( k )
57  ! IF (qq_term /= 'JG' .AND. qq_term /= 'SF') write (*,*) 'specify quadrupole term !'
58 
59  ENDIF
60 
61  !-----------------------------------------------------------------------------
62  ! dipole-quadrupole cross term
63  !-----------------------------------------------------------------------------
64  IF (dipole_quad == 1) THEN
65 
66  IF (dq_term == 'VG') CALL f_dq_vrabec_gross( fdq )
67  ! IF (dq_term /= 'VG' ) write (*,*) 'specify DQ-cross term !'
68 
69  ENDIF
70 
71 END SUBROUTINE f_polar
72 
73 
74 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
75 !
76 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
77 
78 SUBROUTINE f_dd_gross_vrabec( fdd )
79 
80  USE parameters, ONLY: pi, kbol
81  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
82  USE eos_constants, ONLY: ddp2, ddp3, ddp4
83 
84  !-----------------------------------------------------------------------------
85  REAL, INTENT(IN OUT) :: fdd
86  !-----------------------------------------------------------------------------
87  INTEGER :: i, j, k, m
88  INTEGER :: ddit, ddmax
89  REAL :: factor2, factor3
90  REAL :: xijfa, xijkfa, xijf_j, xijkf_j, eij
91  REAL :: fdd2, fdd3
92  REAL, DIMENSION(nc) :: my2dd, my0, alph_tst, z1dd, z2dd, dderror
93  REAL, DIMENSION(nc) :: fdd2m, fdd3m, fdd2m2, fdd3m2, fddm, fddm2
94  REAL, DIMENSION(nc,nc) :: idd2, idd4
95  REAL, DIMENSION(nc,nc,nc) :: idd3
96  !-----------------------------------------------------------------------------
97 
98  fdd = 0.0
99  ddit = 0
100  ddmax = 0 ! value assigned, if polarizable compound is present
101  fddm(:) = 0.0
102  DO i = 1, ncomp
103  IF ( uij(i,i) == 0.0 ) write (*,*) 'F_DD_GROSS_VRABEC: do not use dimensionless units'
104  IF ( uij(i,i) == 0.0 ) stop
105  my2dd(i) = (parame(i,6))**2 *1.e-49 /(uij(i,i)*kbol* mseg(i)*sig_ij(i,i)**3 *1.e-30)
106  alph_tst(i) = parame(i,11) / (mseg(i)*sig_ij(i,i)**3 ) * t/parame(i,3)
107  IF ( alph_tst(i) /= 0.0 ) ddmax = 25 ! set maximum number of polarizable RGT-iterations
108  z1dd(i) = my2dd(i) + 3.0*alph_tst(i)
109  z2dd(i) = 3.0*alph_tst(i)
110  my0(i) = my2dd(i)
111  END DO
112 
113  DO i = 1, ncomp
114  DO j = 1, ncomp
115  idd2(i,j) = 0.0
116  idd4(i,j) = 0.0
117  ! IF (parame(i,6).NE.0.0 .AND. parame(j,6).NE.0.0) THEN
118  DO m = 0, 4
119  idd2(i,j) = idd2(i,j) + ddp2(i,j,m)*eta**m
120  idd4(i,j) = idd4(i,j) + ddp4(i,j,m)*eta**m
121  END DO
122  DO k = 1, ncomp
123  idd3(i,j,k) = 0.0
124  ! IF (parame(k,6).NE.0.0) THEN
125  DO m = 0, 4
126  idd3(i,j,k) = idd3(i,j,k) + ddp3(i,j,k,m)*eta**m
127  END DO
128  ! ENDIF
129  END DO
130  ! ENDIF
131  END DO
132  END DO
133 
134  factor2 = -pi *rho
135  factor3 = -4.0/3.0*pi**2 * rho**2
136 
137 9 CONTINUE
138 
139  fdd2m(:) = 0.0
140  fdd2m2(:) = 0.0
141  fdd3m(:) = 0.0
142  fdd3m2(:) = 0.0
143  fdd2 = 0.0
144  fdd3 = 0.0
145  DO i = 1, ncomp
146  DO j = 1, ncomp
147  ! IF (parame(i,6).NE.0.0 .AND. parame(j,6).NE.0.0) THEN
148  xijfa =x(i)*parame(i,3)/t*parame(i,2)**3 * x(j)*parame(j,3)/t*parame(j,2)**3 &
149  /((parame(i,2)+parame(j,2))/2.0)**3 * (z1dd(i)*z1dd(j)-z2dd(i)*z2dd(j)) ! * (1.0-lij(i,j))
150  eij = (parame(i,3)*parame(j,3))**0.5
151  fdd2= fdd2 + factor2 * xijfa * ( idd2(i,j) + eij/t*idd4(i,j) )
152  xijf_j = parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
153  /((parame(i,2)+parame(j,2))/2.0)**3 ! * (1.0-lij(i,j))
154  fdd2m(i)=fdd2m(i)+4.0*sqrt(my2dd(i))*z1dd(j)*factor2* xijf_j *(idd2(i,j)+eij/t*idd4(i,j))
155  fdd2m2(i)=fdd2m2(i) + 4.0*z1dd(j)*factor2* xijf_j *(idd2(i,j)+eij/t*idd4(i,j))
156  IF (j == i) fdd2m2(i) =fdd2m2(i) +8.0*factor2* xijf_j*my2dd(i) *(idd2(i,j)+eij/t*idd4(i,j))
157  DO k = 1, ncomp
158  ! IF (parame(k,6).NE.0.0) THEN
159  xijkfa = x(i)*parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
160  *x(k)*parame(k,3)/t*parame(k,2)**3 / ((parame(i,2)+parame(j,2))/2.0) &
161  /((parame(i,2)+parame(k,2))/2.0) / ((parame(j,2)+parame(k,2))/2.0) &
162  *(z1dd(i)*z1dd(j)*z1dd(k)-z2dd(i)*z2dd(j)*z2dd(k))
163  ! *(1.0-lij(i,j))*(1.0-lij(i,k))*(1.0-lij(j,k))
164  fdd3 = fdd3 + factor3 * xijkfa * idd3(i,j,k)
165  xijkf_j = parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
166  *x(k)*parame(k,3)/t*parame(k,2)**3 /((parame(i,2)+parame(j,2))/2.0) &
167  /((parame(i,2)+parame(k,2))/2.0) /((parame(j,2)+parame(k,2))/2.0)
168  ! *(1.0-lij(i,j))*(1.0-lij(i,k))*(1.0-lij(j,k))
169  fdd3m(i)=fdd3m(i)+6.0*factor3*sqrt(my2dd(i))*z1dd(j)*z1dd(k) *xijkf_j*idd3(i,j,k)
170  fdd3m2(i)=fdd3m2(i)+6.0*factor3*z1dd(j)*z1dd(k) *xijkf_j*idd3(i,j,k)
171  IF(j == i) fdd3m2(i) =fdd3m2(i)+24.0*factor3*my2dd(i)*z1dd(k) *xijkf_j*idd3(i,j,k)
172  ! ENDIF
173  END DO
174  ! ENDIF
175  END DO
176  END DO
177 
178  IF (fdd2 < -1.e-50 .AND. fdd3 /= 0.0) THEN
179  fdd = fdd2 / ( 1.0 - fdd3/fdd2 )
180  IF ( ddmax /= 0 ) THEN
181  DO i = 1, ncomp
182  ddit = ddit + 1
183  fddm(i) =fdd2*(fdd2*fdd2m(i) -2.0*fdd3*fdd2m(i)+fdd2*fdd3m(i)) /(fdd2-fdd3)**2
184  fddm2(i) = fdd2m(i) * (fdd2*fdd2m(i)-2.0*fdd3*fdd2m(i) +fdd2*fdd3m(i)) / (fdd2-fdd3)**2 &
185  + fdd2*(fdd2*fdd2m2(i) -2.0*fdd3*fdd2m2(i)+fdd2m(i)**2 &
186  -fdd2m(i)*fdd3m(i) +fdd2*fdd3m2(i)) / (fdd2-fdd3)**2 &
187  - 2.0*fdd2*(fdd2*fdd2m(i) -2.0*fdd3*fdd2m(i) +fdd2*fdd3m(i)) /(fdd2-fdd3)**3 &
188  *(fdd2m(i)-fdd3m(i))
189  dderror(i)= sqrt( my2dd(i) ) - sqrt( my0(i) ) + alph_tst(i)*fddm(i)
190  my2dd(i) = ( sqrt( my2dd(i) ) - dderror(i) / (1.0+alph_tst(i)*fddm2(i)) )**2
191  z1dd(i) = my2dd(i) + 3.0 * alph_tst(i)
192  ENDDO
193  DO i = 1, ncomp
194  IF (abs(dderror(i)) > 1.e-11 .AND. ddit < ddmax) GOTO 9
195  ENDDO
196  fdd = fdd + sum( 0.5*x(1:ncomp)*alph_tst(1:ncomp)*fddm(1:ncomp)**2 )
197  ENDIF
198  END IF
199 
200 
201 END SUBROUTINE f_dd_gross_vrabec
202 
203 
204 
205 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
206 !
207 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
208 
209 SUBROUTINE f_qq_gross( fqq )
210 
211  USE parameters, ONLY: pi, kbol
212  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
213  USE eos_constants, ONLY: qqp2, qqp3, qqp4
214 
215  !-----------------------------------------------------------------------------
216  REAL, INTENT(IN OUT) :: fqq
217  !-----------------------------------------------------------------------------
218  INTEGER :: i, j, k, m
219  REAL :: factor2, factor3
220  REAL :: xijfa, xijkfa, eij
221  REAL :: fqq2, fqq3
222  REAL, DIMENSION(nc) :: qq2
223  REAL, DIMENSION(nc,nc) :: iqq2, iqq4
224  REAL, DIMENSION(nc,nc,nc) :: iqq3
225  !-----------------------------------------------------------------------------
226 
227 
228  fqq = 0.0
229  DO i = 1, ncomp
230  qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
231  END DO
232 
233  DO i = 1, ncomp
234  DO j = 1, ncomp
235  iqq2(i,j) = 0.0
236  iqq4(i,j) = 0.0
237  IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) THEN
238  DO m = 0, 4
239  iqq2(i,j) = iqq2(i,j) + qqp2(i,j,m)*eta**m
240  iqq4(i,j) = iqq4(i,j) + qqp4(i,j,m)*eta**m
241  END DO
242  DO k = 1, ncomp
243  iqq3(i,j,k) = 0.0
244  IF (parame(k,7) /= 0.0) THEN
245  DO m = 0, 4
246  iqq3(i,j,k) = iqq3(i,j,k) + qqp3(i,j,k,m)*eta**m
247  END DO
248  END IF
249  END DO
250  END IF
251  END DO
252  END DO
253 
254  factor2 = -9.0/16.0*pi *rho
255  factor3 = 9.0/16.0*pi**2 * rho**2
256 
257  fqq2 = 0.0
258  fqq3 = 0.0
259  DO i = 1, ncomp
260  DO j = 1, ncomp
261  IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) THEN
262  xijfa=x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t &
263  *x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,j)**7
264  eij = (parame(i,3)*parame(j,3))**0.5
265  fqq2= fqq2 +factor2* xijfa * (iqq2(i,j)+eij/t*iqq4(i,j))
266  DO k = 1, ncomp
267  IF (parame(k,7) /= 0.0) THEN
268  xijkfa=x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t/sig_ij(i,j)**3 &
269  *x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,k)**3 &
270  *x(k)*uij(k,k)*qq2(k)*sig_ij(k,k)**5 /t/sig_ij(j,k)**3
271  fqq3 = fqq3 + factor3 * xijkfa * iqq3(i,j,k)
272  END IF
273  END DO
274  END IF
275  END DO
276  END DO
277 
278  IF ( fqq2 < -1.e-50 .AND. fqq3 /= 0.0 ) THEN
279  fqq = fqq2 / ( 1.0 - fqq3/fqq2 )
280  END IF
281 
282 
283 
284 END SUBROUTINE f_qq_gross
285 
286 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
287 !
288 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
289 
290 SUBROUTINE f_dq_vrabec_gross( fdq )
291 
292  USE parameters, ONLY: pi, kbol
293  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
294  USE eos_constants, ONLY: dqp2, dqp3, dqp4
295 
296  !-----------------------------------------------------------------------------
297  REAL, INTENT(IN OUT) :: fdq
298  !-----------------------------------------------------------------------------
299  INTEGER :: i, j, k, m
300  REAL :: factor2, factor3
301  REAL :: xijfa, xijkfa, eij
302  REAL :: fdq2, fdq3
303  REAL, DIMENSION(nc) :: my2dd, myfac, qq2, q_fac
304  REAL, DIMENSION(nc,nc) :: idq2, idq4
305  REAL, DIMENSION(nc,nc,nc) :: idq3
306  !-----------------------------------------------------------------------------
307 
308 
309  fdq = 0.0
310  DO i = 1, ncomp
311  my2dd(i) = (parame(i,6))**2 *1.e-49 /(uij(i,i)*kbol* mseg(i)*sig_ij(i,i)**3 *1.e-30)
312  myfac(i) = parame(i,3)/t*parame(i,2)**4 *my2dd(i)
313  ! myfac(i)=parame(i,3)/T*parame(i,2)**4 *my2dd_renormalized(i)
314  qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
315  q_fac(i) = parame(i,3)/t*parame(i,2)**4 *qq2(i)
316  END DO
317 
318  DO i = 1, ncomp
319  DO j = 1, ncomp
320  idq2(i,j) = 0.0
321  idq4(i,j) = 0.0
322  IF (myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0) THEN
323  DO m = 0, 4
324  idq2(i,j) = idq2(i,j) + dqp2(i,j,m)*eta**m
325  idq4(i,j) = idq4(i,j) + dqp4(i,j,m)*eta**m
326  END DO
327  DO k = 1, ncomp
328  idq3(i,j,k) = 0.0
329  IF (myfac(k) /= 0.0 .OR. q_fac(k) /= 0.0) THEN
330  DO m = 0, 4
331  idq3(i,j,k) = idq3(i,j,k) + dqp3(i,j,k,m)*eta**m
332  END DO
333  END IF
334  END DO
335  END IF
336  END DO
337  END DO
338 
339  factor2 = -9.0/4.0 * pi *rho
340  factor3 = pi**2 * rho**2
341 
342  fdq2 = 0.0
343  fdq3 = 0.0
344  DO i = 1, ncomp
345  DO j = 1, ncomp
346  IF (myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0) THEN
347  xijfa = x(i)*myfac(i) * x(j)*q_fac(j) /sig_ij(i,j)**5
348  eij = (parame(i,3)*parame(j,3))**0.5
349  fdq2 = fdq2 +factor2* xijfa*(idq2(i,j)+eij/t*idq4(i,j))
350  DO k = 1, ncomp
351  IF (myfac(k) /= 0.0 .OR. q_fac(k) /= 0.0) THEN
352  xijkfa=x(i)*x(j)*x(k)/(sig_ij(i,j)*sig_ij(i,k)*sig_ij(j,k))**2 &
353  *( myfac(i)*q_fac(j)*myfac(k) + myfac(i)*q_fac(j)*q_fac(k)*1.1937350 )
354  fdq3 = fdq3 + factor3*xijkfa*idq3(i,j,k)
355  END IF
356  END DO
357  END IF
358  END DO
359  END DO
360 
361  IF (fdq2 < -1.e-50 .AND. fdq3 /= 0.0) THEN
362  fdq = fdq2 / ( 1.0 - fdq3/fdq2 )
363  END IF
364 
365 END SUBROUTINE f_dq_vrabec_gross
366 
367 
368 
369 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
370 !
371 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
372 
373 SUBROUTINE p_polar ( zdd, zddz, zddz2, zddz3, zqq, zqqz, zqqz2, zqqz3, zdq, zdqz, zdqz2, zdqz3 )
374 
375  USE eos_variables, ONLY: ncomp, parame, dd_term, qq_term, dq_term
376 
377  !-----------------------------------------------------------------------------
378  REAL, INTENT(OUT) :: zdd, zddz, zddz2, zddz3
379  REAL, INTENT(OUT) :: zqq, zqqz, zqqz2, zqqz3
380  REAL, INTENT(OUT) :: zdq, zdqz, zdqz2, zdqz3
381 
382  !-----------------------------------------------------------------------------
383  INTEGER :: dipole
384  INTEGER :: quadrupole
385  INTEGER :: dipole_quad
386  !-----------------------------------------------------------------------------
387 
388  zdd = 0.0
389  zddz = 0.0
390  zddz2 = 0.0
391  zddz3 = 0.0
392  zqq = 0.0
393  zqqz = 0.0
394  zqqz2 = 0.0
395  zqqz3 = 0.0
396  zdq = 0.0
397  zdqz = 0.0
398  zdqz2 = 0.0
399  zdqz3 = 0.0
400 
401  dipole = 0
402  quadrupole = 0
403  dipole_quad = 0
404  IF ( sum( parame(1:ncomp,6) ) /= 0.0 ) dipole = 1
405  IF ( sum( parame(1:ncomp,7) ) /= 0.0 ) quadrupole = 1
406  IF ( dipole == 1 .AND. quadrupole == 1 ) dipole_quad = 1
407 
408  !-----------------------------------------------------------------------------
409  ! dipole-dipole term
410  !-----------------------------------------------------------------------------
411  IF (dipole == 1) THEN
412 
413  IF (dd_term == 'GV') CALL p_dd_gross_vrabec( zdd, zddz, zddz2, zddz3 )
414  ! IF (dd_term == 'SF') CALL F_DD_SAAGER_FISCHER( k )
415  IF (dd_term /= 'GV' .AND. dd_term /= 'SF') write (*,*) 'specify dipole term !'
416 
417  ENDIF
418 
419  !-----------------------------------------------------------------------------
420  ! quadrupole-quadrupole term
421  !-----------------------------------------------------------------------------
422  IF (quadrupole == 1) THEN
423 
424  !IF (qq_term == 'SF') CALL F_QQ_SAAGER_FISCHER( k )
425  IF (qq_term == 'JG') CALL p_qq_gross( zqq, zqqz, zqqz2, zqqz3 )
426  IF (qq_term /= 'JG' .AND. qq_term /= 'SF') write (*,*) 'specify quadrupole term !'
427 
428  ENDIF
429 
430  !-----------------------------------------------------------------------------
431  ! dipole-quadrupole cross term
432  !-----------------------------------------------------------------------------
433  IF (dipole_quad == 1) THEN
434 
435  IF (dq_term == 'VG') CALL p_dq_vrabec_gross( zdq, zdqz, zdqz2, zdqz3 )
436  IF (dq_term /= 'VG' ) write (*,*) 'specify DQ-cross term !'
437 
438  ENDIF
439 
440 END SUBROUTINE p_polar
441 
442 
443 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
444 !
445 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
446 
447 SUBROUTINE p_dd_gross_vrabec( zdd, zddz, zddz2, zddz3 )
448 
449  USE parameters, ONLY: pi, kbol
450  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
451  USE eos_constants, ONLY: ddp2, ddp3, ddp4
452 
453  !-----------------------------------------------------------------------------
454  REAL, INTENT(IN OUT) :: zdd, zddz, zddz2, zddz3
455  !-----------------------------------------------------------------------------
456  INTEGER :: i, j, k, m
457  REAL :: factor2, factor3, z3
458  REAL :: xijfa, xijkfa, eij
459  REAL :: fdddr, fddd2, fddd3, fddd4
460  REAL :: fdd2, fdd2z, fdd2z2, fdd2z3, fdd2z4
461  REAL :: fdd3, fdd3z, fdd3z2, fdd3z3, fdd3z4
462  REAL, DIMENSION(nc) :: my2dd
463  REAL, DIMENSION(nc,nc) :: idd2, idd2z, idd2z2, idd2z3, idd2z4
464  REAL, DIMENSION(nc,nc) :: idd4, idd4z, idd4z2, idd4z3, idd4z4
465  REAL, DIMENSION(nc,nc,nc) :: idd3, idd3z, idd3z2, idd3z3, idd3z4
466  !-----------------------------------------------------------------------------
467 
468 
469  zdd = 0.0
470  zddz = 0.0
471  zddz2 = 0.0
472  zddz3 = 0.0
473  z3 = eta
474  DO i = 1, ncomp
475  my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
476  END DO
477 
478  DO i = 1, ncomp
479  DO j = 1, ncomp
480  idd2(i,j) = 0.0
481  idd4(i,j) = 0.0
482  idd2z(i,j) = 0.0
483  idd4z(i,j) = 0.0
484  idd2z2(i,j) = 0.0
485  idd4z2(i,j) = 0.0
486  idd2z3(i,j) = 0.0
487  idd4z3(i,j) = 0.0
488  idd2z4(i,j) = 0.0
489  idd4z4(i,j) = 0.0
490  ! IF (paramei,6).NE.0.0 .AND. parame(j,6).NE.0.0) THEN
491  DO m = 0, 4
492  idd2(i,j) =idd2(i,j) + ddp2(i,j,m) *z3**(m+1)
493  idd4(i,j) =idd4(i,j) + ddp4(i,j,m) *z3**(m+1)
494  idd2z(i,j) =idd2z(i,j) +ddp2(i,j,m)*REAL(m+1) *z3**m
495  idd4z(i,j) =idd4z(i,j) +ddp4(i,j,m)*REAL(m+1) *z3**m
496  idd2z2(i,j)=idd2z2(i,j)+ddp2(i,j,m)*REAL((m+1)*m) *z3**(m-1)
497  idd4z2(i,j)=idd4z2(i,j)+ddp4(i,j,m)*REAL((m+1)*m) *z3**(m-1)
498  idd2z3(i,j)=idd2z3(i,j)+ddp2(i,j,m)*REAL((m+1)*m*(m-1)) *z3**(m-2)
499  idd4z3(i,j)=idd4z3(i,j)+ddp4(i,j,m)*REAL((m+1)*m*(m-1)) *z3**(m-2)
500  idd2z4(i,j)=idd2z4(i,j)+ddp2(i,j,m)*REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
501  idd4z4(i,j)=idd4z4(i,j)+ddp4(i,j,m)*REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
502  END DO
503  DO k = 1, ncomp
504  idd3(i,j,k) = 0.0
505  idd3z(i,j,k) = 0.0
506  idd3z2(i,j,k) = 0.0
507  idd3z3(i,j,k) = 0.0
508  idd3z4(i,j,k) = 0.0
509  ! IF (parame(k,6).NE.0.0) THEN
510  DO m = 0, 4
511  idd3(i,j,k) =idd3(i,j,k) +ddp3(i,j,k,m)*z3**(m+2)
512  idd3z(i,j,k) =idd3z(i,j,k) +ddp3(i,j,k,m)*REAL(m+2)*z3**(m+1)
513  idd3z2(i,j,k)=idd3z2(i,j,k)+ddp3(i,j,k,m)*REAL((m+2)*(m+1))*z3**m
514  idd3z3(i,j,k)=idd3z3(i,j,k)+ddp3(i,j,k,m)*REAL((m+2)*(m+1)*m)*z3**(m-1)
515  idd3z4(i,j,k)=idd3z4(i,j,k)+ddp3(i,j,k,m)*REAL((m+2)*(m+1)*m*(m-1)) *z3**(m-2)
516  END DO
517  ! ENDIF
518  END DO
519  ! ENDIF
520  END DO
521  END DO
522 
523  factor2= -pi *rho/z3
524  factor3= -4.0/3.0*pi**2 * (rho/z3)**2
525 
526  fdd2 = 0.0
527  fdd2z = 0.0
528  fdd2z2 = 0.0
529  fdd2z3 = 0.0
530  fdd2z4 = 0.0
531  fdd3 = 0.0
532  fdd3z = 0.0
533  fdd3z2 = 0.0
534  fdd3z3 = 0.0
535  fdd3z4 = 0.0
536  DO i = 1, ncomp
537  DO j = 1, ncomp
538  ! IF (parame(i,6).NE.0.0 .AND. parame(j,6).NE.0.0) THEN
539  xijfa = x(i)*parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
540  / ((parame(i,2)+parame(j,2))/2.0)**3 *my2dd(i)*my2dd(j)
541  eij = (parame(i,3)*parame(j,3))**0.5
542  fdd2 = fdd2 +factor2*xijfa*(idd2(i,j) +eij/t*idd4(i,j))
543  fdd2z = fdd2z +factor2*xijfa*(idd2z(i,j) +eij/t*idd4z(i,j))
544  fdd2z2 = fdd2z2+factor2*xijfa*(idd2z2(i,j)+eij/t*idd4z2(i,j))
545  fdd2z3 = fdd2z3+factor2*xijfa*(idd2z3(i,j)+eij/t*idd4z3(i,j))
546  fdd2z4 = fdd2z4+factor2*xijfa*(idd2z4(i,j)+eij/t*idd4z4(i,j))
547  DO k = 1, ncomp
548  ! IF (parame(k,6).NE.0.0) THEN
549  xijkfa= x(i)*parame(i,3)/t*parame(i,2)**3 *x(j)*parame(j,3)/t*parame(j,2)**3 &
550  *x(k)*parame(k,3)/t*parame(k,2)**3 /((parame(i,2)+parame(j,2))/2.0) &
551  /((parame(i,2)+parame(k,2))/2.0) /((parame(j,2)+parame(k,2))/2.0) &
552  *my2dd(i)*my2dd(j)*my2dd(k)
553  fdd3 = fdd3 + factor3 * xijkfa*idd3(i,j,k)
554  fdd3z = fdd3z + factor3 * xijkfa*idd3z(i,j,k)
555  fdd3z2 = fdd3z2 + factor3 * xijkfa*idd3z2(i,j,k)
556  fdd3z3 = fdd3z3 + factor3 * xijkfa*idd3z3(i,j,k)
557  fdd3z4 = fdd3z4 + factor3 * xijkfa*idd3z4(i,j,k)
558  ! ENDIF
559  END DO
560  ! ENDIF
561  END DO
562  END DO
563  IF (fdd2 < -1.e-50 .AND. fdd3 /= 0.0 .AND. fdd2z /= 0.0 .AND. fdd3z /= 0.0) THEN
564 
565  fdddr= fdd2* (fdd2*fdd2z - 2.0*fdd3*fdd2z+fdd2*fdd3z) / (fdd2-fdd3)**2
566  fddd2=(2.0*fdd2*fdd2z*fdd2z +fdd2*fdd2*fdd2z2 &
567  -2.0*fdd2z**2 *fdd3-2.0*fdd2*fdd2z2*fdd3+fdd2*fdd2*fdd3z2) &
568  /(fdd2-fdd3)**2 + fdddr * 2.0*(fdd3z-fdd2z)/(fdd2-fdd3)
569  fddd3=(2.0*fdd2z**3 +6.0*fdd2*fdd2z*fdd2z2+fdd2*fdd2*fdd2z3 &
570  -6.0*fdd2z*fdd2z2*fdd3-2.0*fdd2z**2 *fdd3z &
571  -2.0*fdd2*fdd2z3*fdd3 -2.0*fdd2*fdd2z2*fdd3z &
572  +2.0*fdd2*fdd2z*fdd3z2+fdd2*fdd2*fdd3z3) /(fdd2-fdd3)**2 &
573  + 2.0/(fdd2-fdd3)* ( 2.0*fddd2*(fdd3z-fdd2z) &
574  + fdddr*(fdd3z2-fdd2z2) &
575  - fdddr/(fdd2-fdd3)*(fdd3z-fdd2z)**2 )
576  fddd4=( 12.0*fdd2z**2 *fdd2z2+6.0*fdd2*fdd2z2**2 &
577  +8.0*fdd2*fdd2z*fdd2z3+fdd2*fdd2*fdd2z4-6.0*fdd2z2**2 *fdd3 &
578  -12.0*fdd2z*fdd2z2*fdd3z -8.0*fdd2z*fdd2z3*fdd3 &
579  -2.0*fdd2*fdd2z4*fdd3-4.0*fdd2*fdd2z3*fdd3z &
580  +4.0*fdd2*fdd2z*fdd3z3+fdd2**2 *fdd3z4 ) /(fdd2-fdd3)**2 &
581  + 6.0/(fdd2-fdd3)* ( fddd3*(fdd3z-fdd2z) &
582  -fddd2/(fdd2-fdd3)*(fdd3z-fdd2z)**2 &
583  - fdddr/(fdd2-fdd3)*(fdd3z-fdd2z)*(fdd3z2-fdd2z2) &
584  + fddd2*(fdd3z2-fdd2z2) +1.0/3.0*fdddr*(fdd3z3-fdd2z3) )
585  zdd = fdddr*eta
586  zddz = fddd2*eta + fdddr
587  zddz2 = fddd3*eta + 2.0* fddd2
588  zddz3 = fddd4*eta + 3.0* fddd3
589 
590  END IF
591 
592 
593 END SUBROUTINE p_dd_gross_vrabec
594 
595 
596 
597 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
598 !
599 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
600 
601 SUBROUTINE p_qq_gross( zqq, zqqz, zqqz2, zqqz3 )
602 
603  USE parameters, ONLY: pi, kbol
604  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
605  USE eos_constants, ONLY: qqp2, qqp3, qqp4
606 
607  !-----------------------------------------------------------------------------
608  REAL, INTENT(IN OUT) :: zqq, zqqz, zqqz2, zqqz3
609  !-----------------------------------------------------------------------------
610  INTEGER :: i, j, k, m
611  REAL :: factor2, factor3, z3
612  REAL :: xijfa, xijkfa, eij
613  REAL :: fqqdr, fqqd2, fqqd3, fqqd4
614  REAL :: fqq2, fqq2z, fqq2z2, fqq2z3, fqq2z4
615  REAL :: fqq3, fqq3z, fqq3z2, fqq3z3, fqq3z4
616  REAL, DIMENSION(nc) :: qq2
617  REAL, DIMENSION(nc,nc) :: iqq2, iqq2z, iqq2z2, iqq2z3, iqq2z4
618  REAL, DIMENSION(nc,nc) :: iqq4, iqq4z, iqq4z2, iqq4z3, iqq4z4
619  REAL, DIMENSION(nc,nc,nc) :: iqq3, iqq3z, iqq3z2, iqq3z3, iqq3z4
620  !-----------------------------------------------------------------------------
621 
622  zqq = 0.0
623  zqqz = 0.0
624  zqqz2 = 0.0
625  zqqz3 = 0.0
626  z3 = eta
627  DO i=1,ncomp
628  qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
629  END DO
630 
631  DO i = 1, ncomp
632  DO j = 1, ncomp
633  iqq2(i,j) = 0.0
634  iqq4(i,j) = 0.0
635  iqq2z(i,j) = 0.0
636  iqq4z(i,j) = 0.0
637  iqq2z2(i,j) = 0.0
638  iqq4z2(i,j) = 0.0
639  iqq2z3(i,j) = 0.0
640  iqq4z3(i,j) = 0.0
641  iqq2z4(i,j) = 0.0
642  iqq4z4(i,j) = 0.0
643  IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) THEN
644  DO m = 0, 4
645  iqq2(i,j) =iqq2(i,j) + qqp2(i,j,m)*z3**(m+1)
646  iqq4(i,j) =iqq4(i,j) + qqp4(i,j,m)*z3**(m+1)
647  iqq2z(i,j) =iqq2z(i,j) +qqp2(i,j,m)*REAL(m+1)*z3**m
648  iqq4z(i,j) =iqq4z(i,j) +qqp4(i,j,m)*REAL(m+1)*z3**m
649  iqq2z2(i,j)=iqq2z2(i,j)+qqp2(i,j,m)*REAL((m+1)*m)*z3**(m-1)
650  iqq4z2(i,j)=iqq4z2(i,j)+qqp4(i,j,m)*REAL((m+1)*m)*z3**(m-1)
651  iqq2z3(i,j)=iqq2z3(i,j)+qqp2(i,j,m)*REAL((m+1)*m*(m-1)) *z3**(m-2)
652  iqq4z3(i,j)=iqq4z3(i,j)+qqp4(i,j,m)*REAL((m+1)*m*(m-1)) *z3**(m-2)
653  iqq2z4(i,j)=iqq2z4(i,j)+qqp2(i,j,m)*REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
654  iqq4z4(i,j)=iqq4z4(i,j)+qqp4(i,j,m)*REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
655  END DO
656  DO k=1,ncomp
657  iqq3(i,j,k) = 0.0
658  iqq3z(i,j,k) = 0.0
659  iqq3z2(i,j,k) = 0.0
660  iqq3z3(i,j,k) = 0.0
661  iqq3z4(i,j,k) = 0.0
662  IF (parame(k,7) /= 0.0) THEN
663  DO m=0,4
664  iqq3(i,j,k) =iqq3(i,j,k) + qqp3(i,j,k,m)*z3**(m+2)
665  iqq3z(i,j,k)=iqq3z(i,j,k)+qqp3(i,j,k,m)*REAL(m+2)*z3**(m+1)
666  iqq3z2(i,j,k)=iqq3z2(i,j,k)+qqp3(i,j,k,m)*REAL((m+2)*(m+1)) *z3**m
667  iqq3z3(i,j,k)=iqq3z3(i,j,k)+qqp3(i,j,k,m)*REAL((m+2)*(m+1)*m) *z3**(m-1)
668  iqq3z4(i,j,k)=iqq3z4(i,j,k)+qqp3(i,j,k,m) *REAL((m+2)*(m+1)*m*(m-1)) *z3**(m-2)
669  END DO
670  END IF
671  END DO
672 
673  END IF
674  END DO
675  END DO
676 
677  factor2= -9.0/16.0*pi *rho/z3
678  factor3= 9.0/16.0*pi**2 * (rho/z3)**2
679 
680  fqq2 = 0.0
681  fqq2z = 0.0
682  fqq2z2 = 0.0
683  fqq2z3 = 0.0
684  fqq2z4 = 0.0
685  fqq3 = 0.0
686  fqq3z = 0.0
687  fqq3z2 = 0.0
688  fqq3z3 = 0.0
689  fqq3z4 = 0.0
690  DO i = 1, ncomp
691  DO j = 1, ncomp
692  IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) THEN
693  xijfa =x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t &
694  *x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,j)**7
695  eij = (parame(i,3)*parame(j,3))**0.5
696  fqq2= fqq2 +factor2*xijfa*(iqq2(i,j) +eij/t*iqq4(i,j) )
697  fqq2z =fqq2z +factor2*xijfa*(iqq2z(i,j) +eij/t*iqq4z(i,j) )
698  fqq2z2=fqq2z2+factor2*xijfa*(iqq2z2(i,j)+eij/t*iqq4z2(i,j))
699  fqq2z3=fqq2z3+factor2*xijfa*(iqq2z3(i,j)+eij/t*iqq4z3(i,j))
700  fqq2z4=fqq2z4+factor2*xijfa*(iqq2z4(i,j)+eij/t*iqq4z4(i,j))
701  DO k = 1, ncomp
702  IF (parame(k,7) /= 0.0) THEN
703  xijkfa=x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t/sig_ij(i,j)**3 &
704  *x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,k)**3 &
705  *x(k)*uij(k,k)*qq2(k)*sig_ij(k,k)**5 /t/sig_ij(j,k)**3
706  fqq3 = fqq3 + factor3 * xijkfa*iqq3(i,j,k)
707  fqq3z = fqq3z + factor3 * xijkfa*iqq3z(i,j,k)
708  fqq3z2 = fqq3z2 + factor3 * xijkfa*iqq3z2(i,j,k)
709  fqq3z3 = fqq3z3 + factor3 * xijkfa*iqq3z3(i,j,k)
710  fqq3z4 = fqq3z4 + factor3 * xijkfa*iqq3z4(i,j,k)
711  END IF
712  END DO
713  END IF
714  END DO
715  END DO
716  IF (fqq2 < -1.e-50 .AND. fqq3 /= 0.0 .AND. fqq2z /= 0.0 .AND. fqq3z /= 0.0) THEN
717  fqqdr = fqq2* (fqq2*fqq2z - 2.0*fqq3*fqq2z+fqq2*fqq3z) /(fqq2-fqq3)**2
718  fqqd2= (2.0*fqq2*fqq2z*fqq2z +fqq2*fqq2*fqq2z2 &
719  -2.0*fqq2z**2 *fqq3-2.0*fqq2*fqq2z2*fqq3+fqq2*fqq2*fqq3z2) &
720  /(fqq2-fqq3)**2 + fqqdr * 2.0*(fqq3z-fqq2z)/(fqq2-fqq3)
721  fqqd3=(2.0*fqq2z**3 +6.0*fqq2*fqq2z*fqq2z2+fqq2*fqq2*fqq2z3 &
722  -6.0*fqq2z*fqq2z2*fqq3-2.0*fqq2z**2 *fqq3z &
723  -2.0*fqq2*fqq2z3*fqq3 -2.0*fqq2*fqq2z2*fqq3z &
724  +2.0*fqq2*fqq2z*fqq3z2+fqq2*fqq2*fqq3z3) /(fqq2-fqq3)**2 &
725  + 2.0/(fqq2-fqq3)* ( 2.0*fqqd2*(fqq3z-fqq2z) &
726  + fqqdr*(fqq3z2-fqq2z2) - fqqdr/(fqq2-fqq3)*(fqq3z-fqq2z)**2 )
727  fqqd4=( 12.0*fqq2z**2 *fqq2z2+6.0*fqq2*fqq2z2**2 &
728  +8.0*fqq2*fqq2z*fqq2z3+fqq2*fqq2*fqq2z4-6.0*fqq2z2**2 *fqq3 &
729  -12.0*fqq2z*fqq2z2*fqq3z -8.0*fqq2z*fqq2z3*fqq3 &
730  -2.0*fqq2*fqq2z4*fqq3-4.0*fqq2*fqq2z3*fqq3z &
731  +4.0*fqq2*fqq2z*fqq3z3+fqq2**2 *fqq3z4 ) /(fqq2-fqq3)**2 &
732  + 6.0/(fqq2-fqq3)* ( fqqd3*(fqq3z-fqq2z) &
733  -fqqd2/(fqq2-fqq3)*(fqq3z-fqq2z)**2 &
734  - fqqdr/(fqq2-fqq3)*(fqq3z-fqq2z)*(fqq3z2-fqq2z2) &
735  + fqqd2*(fqq3z2-fqq2z2) +1.0/3.0*fqqdr*(fqq3z3-fqq2z3) )
736  zqq = fqqdr*eta
737  zqqz = fqqd2*eta + fqqdr
738  zqqz2 = fqqd3*eta + 2.0* fqqd2
739  zqqz3 = fqqd4*eta + 3.0* fqqd3
740  END IF
741 
742 
743 END SUBROUTINE p_qq_gross
744 
745 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
746 !
747 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
748 
749 SUBROUTINE p_dq_vrabec_gross( zdq, zdqz, zdqz2, zdqz3 )
750 
751  USE parameters, ONLY: pi, kbol
752  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
753  USE eos_constants, ONLY: dqp2, dqp3, dqp4
754 
755  !-----------------------------------------------------------------------------
756  REAL, INTENT(IN OUT) :: zdq, zdqz, zdqz2, zdqz3
757  !-----------------------------------------------------------------------------
758  INTEGER :: i, j, k, m
759  REAL :: factor2, factor3, z3
760  REAL :: xijfa, xijkfa, eij
761  REAL :: fdqdr, fdqd2, fdqd3, fdqd4
762  REAL :: fdq2, fdq2z, fdq2z2, fdq2z3, fdq2z4
763  REAL :: fdq3, fdq3z, fdq3z2, fdq3z3, fdq3z4
764  REAL, DIMENSION(nc) :: my2dd, myfac, qq2, q_fac
765  REAL, DIMENSION(nc,nc) :: idq2, idq2z, idq2z2, idq2z3, idq2z4
766  REAL, DIMENSION(nc,nc) :: idq4, idq4z, idq4z2, idq4z3, idq4z4
767  REAL, DIMENSION(nc,nc,nc) :: idq3, idq3z, idq3z2, idq3z3, idq3z4
768  !-----------------------------------------------------------------------------
769 
770  zdq = 0.0
771  zdqz = 0.0
772  zdqz2 = 0.0
773  zdqz3 = 0.0
774  z3 = eta
775  DO i = 1, ncomp
776  my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
777  myfac(i) = parame(i,3)/t*parame(i,2)**4 *my2dd(i)
778  qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
779  q_fac(i) = parame(i,3)/t*parame(i,2)**4 *qq2(i)
780  END DO
781 
782 
783  DO i = 1, ncomp
784  DO j = 1, ncomp
785  idq2(i,j) = 0.0
786  idq4(i,j) = 0.0
787  idq2z(i,j) = 0.0
788  idq4z(i,j) = 0.0
789  idq2z2(i,j) = 0.0
790  idq4z2(i,j) = 0.0
791  idq2z3(i,j) = 0.0
792  idq4z3(i,j) = 0.0
793  idq2z4(i,j) = 0.0
794  idq4z4(i,j) = 0.0
795  IF (myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0) THEN
796  DO m = 0, 4
797  idq2(i,j) =idq2(i,j) + dqp2(i,j,m)*z3**(m+1)
798  idq4(i,j) =idq4(i,j) + dqp4(i,j,m)*z3**(m+1)
799  idq2z(i,j) =idq2z(i,j) +dqp2(i,j,m)*REAL(m+1)*z3**m
800  idq4z(i,j) =idq4z(i,j) +dqp4(i,j,m)*REAL(m+1)*z3**m
801  idq2z2(i,j)=idq2z2(i,j)+dqp2(i,j,m)*REAL((m+1)*m)*z3**(m-1)
802  idq4z2(i,j)=idq4z2(i,j)+dqp4(i,j,m)*REAL((m+1)*m)*z3**(m-1)
803  idq2z3(i,j)=idq2z3(i,j)+dqp2(i,j,m)*REAL((m+1)*m*(m-1)) *z3**(m-2)
804  idq4z3(i,j)=idq4z3(i,j)+dqp4(i,j,m)*REAL((m+1)*m*(m-1)) *z3**(m-2)
805  idq2z4(i,j)=idq2z4(i,j)+dqp2(i,j,m)*REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
806  idq4z4(i,j)=idq4z4(i,j)+dqp4(i,j,m)*REAL((m+1)*m*(m-1)*(m-2)) *z3**(m-3)
807  END DO
808  DO k = 1, ncomp
809  idq3(i,j,k) = 0.0
810  idq3z(i,j,k) = 0.0
811  idq3z2(i,j,k) = 0.0
812  idq3z3(i,j,k) = 0.0
813  idq3z4(i,j,k) = 0.0
814  IF (myfac(k) /= 0.0.OR.q_fac(k) /= 0.0) THEN
815  DO m = 0, 4
816  idq3(i,j,k) =idq3(i,j,k) + dqp3(i,j,k,m)*z3**(m+2)
817  idq3z(i,j,k)=idq3z(i,j,k)+dqp3(i,j,k,m)*REAL(m+2)*z3**(m+1)
818  idq3z2(i,j,k)=idq3z2(i,j,k)+dqp3(i,j,k,m)*REAL((m+2)*(m+1)) *z3**m
819  idq3z3(i,j,k)=idq3z3(i,j,k)+dqp3(i,j,k,m)*REAL((m+2)*(m+1)*m) *z3**(m-1)
820  idq3z4(i,j,k)=idq3z4(i,j,k)+dqp3(i,j,k,m) &
821  *REAL((m+2)*(m+1)*m*(m-1)) *z3**(m-2)
822  END DO
823  END IF
824  END DO
825 
826  END IF
827  END DO
828  END DO
829 
830  factor2= -9.0/4.0*pi *rho/z3
831  factor3= pi**2 * (rho/z3)**2
832 
833  fdq2 = 0.0
834  fdq2z = 0.0
835  fdq2z2 = 0.0
836  fdq2z3 = 0.0
837  fdq2z4 = 0.0
838  fdq3 = 0.0
839  fdq3z = 0.0
840  fdq3z2 = 0.0
841  fdq3z3 = 0.0
842  fdq3z4 = 0.0
843  DO i = 1, ncomp
844  DO j = 1, ncomp
845  IF (myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0) THEN
846  xijfa =x(i)*myfac(i) * x(j)*q_fac(j) /sig_ij(i,j)**5
847  eij = (parame(i,3)*parame(j,3))**0.5
848  fdq2= fdq2 +factor2*xijfa*(idq2(i,j) +eij/t*idq4(i,j) )
849  fdq2z =fdq2z +factor2*xijfa*(idq2z(i,j) +eij/t*idq4z(i,j) )
850  fdq2z2=fdq2z2+factor2*xijfa*(idq2z2(i,j)+eij/t*idq4z2(i,j))
851  fdq2z3=fdq2z3+factor2*xijfa*(idq2z3(i,j)+eij/t*idq4z3(i,j))
852  fdq2z4=fdq2z4+factor2*xijfa*(idq2z4(i,j)+eij/t*idq4z4(i,j))
853  DO k = 1, ncomp
854  IF (myfac(k) /= 0.0.OR.q_fac(k) /= 0.0) THEN
855  xijkfa=x(i)*x(j)*x(k)/(sig_ij(i,j)*sig_ij(i,k)*sig_ij(j,k))**2 &
856  *( myfac(i)*q_fac(j)*myfac(k) + myfac(i)*q_fac(j)*q_fac(k)*1.193735 )
857  fdq3 =fdq3 + factor3 * xijkfa*idq3(i,j,k)
858  fdq3z =fdq3z + factor3 * xijkfa*idq3z(i,j,k)
859  fdq3z2=fdq3z2 + factor3 * xijkfa*idq3z2(i,j,k)
860  fdq3z3=fdq3z3 + factor3 * xijkfa*idq3z3(i,j,k)
861  fdq3z4=fdq3z4 + factor3 * xijkfa*idq3z4(i,j,k)
862  END IF
863  END DO
864  END IF
865  END DO
866  END DO
867  IF (fdq2 < -1.e-50 .AND. fdq3 /= 0.0 .AND. fdq2z /= 0.0 .AND. fdq3z /= 0.0) THEN
868  fdqdr = fdq2* (fdq2*fdq2z - 2.0*fdq3*fdq2z+fdq2*fdq3z) /(fdq2-fdq3)**2
869  fdqd2= (2.0*fdq2*fdq2z*fdq2z +fdq2*fdq2*fdq2z2 &
870  -2.0*fdq2z**2 *fdq3-2.0*fdq2*fdq2z2*fdq3+fdq2*fdq2*fdq3z2) &
871  /(fdq2-fdq3)**2 + fdqdr * 2.0*(fdq3z-fdq2z)/(fdq2-fdq3)
872  fdqd3=(2.0*fdq2z**3 +6.0*fdq2*fdq2z*fdq2z2+fdq2*fdq2*fdq2z3 &
873  -6.0*fdq2z*fdq2z2*fdq3-2.0*fdq2z**2 *fdq3z &
874  -2.0*fdq2*fdq2z3*fdq3 -2.0*fdq2*fdq2z2*fdq3z &
875  +2.0*fdq2*fdq2z*fdq3z2+fdq2*fdq2*fdq3z3) /(fdq2-fdq3)**2 &
876  + 2.0/(fdq2-fdq3)* ( 2.0*fdqd2*(fdq3z-fdq2z) &
877  + fdqdr*(fdq3z2-fdq2z2) - fdqdr/(fdq2-fdq3)*(fdq3z-fdq2z)**2 )
878  fdqd4=( 12.0*fdq2z**2 *fdq2z2+6.0*fdq2*fdq2z2**2 &
879  +8.0*fdq2*fdq2z*fdq2z3+fdq2*fdq2*fdq2z4-6.0*fdq2z2**2 *fdq3 &
880  -12.0*fdq2z*fdq2z2*fdq3z -8.0*fdq2z*fdq2z3*fdq3 &
881  -2.0*fdq2*fdq2z4*fdq3-4.0*fdq2*fdq2z3*fdq3z &
882  +4.0*fdq2*fdq2z*fdq3z3+fdq2**2 *fdq3z4 ) /(fdq2-fdq3)**2 &
883  + 6.0/(fdq2-fdq3)* ( fdqd3*(fdq3z-fdq2z) &
884  -fdqd2/(fdq2-fdq3)*(fdq3z-fdq2z)**2 &
885  - fdqdr/(fdq2-fdq3)*(fdq3z-fdq2z)*(fdq3z2-fdq2z2) &
886  + fdqd2*(fdq3z2-fdq2z2) +1.0/3.0*fdqdr*(fdq3z3-fdq2z3) )
887  zdq = fdqdr*eta
888  zdqz = fdqd2*eta + fdqdr
889  zdqz2 = fdqd3*eta + 2.0* fdqd2
890  zdqz3 = fdqd4*eta + 3.0* fdqd3
891  END IF
892 
893 
894 END SUBROUTINE p_dq_vrabec_gross
895 
896 
897 
898 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
899 !
900 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
901 
902 SUBROUTINE phi_polar ( k, z3_rk, fdd_rk, fqq_rk, fdq_rk )
903 
904  USE eos_variables, ONLY: ncomp, parame, dd_term, qq_term, dq_term
905 
906  !-----------------------------------------------------------------------------
907  INTEGER, INTENT(IN) :: k
908  REAL, INTENT(IN) :: z3_rk
909  REAL, INTENT(OUT) :: fdd_rk, fqq_rk, fdq_rk
910 
911  !-----------------------------------------------------------------------------
912  INTEGER :: dipole
913  INTEGER :: quadrupole
914  INTEGER :: dipole_quad
915  !-----------------------------------------------------------------------------
916 
917  fdd_rk = 0.0
918  fqq_rk = 0.0
919  fdq_rk = 0.0
920 
921  dipole = 0
922  quadrupole = 0
923  dipole_quad = 0
924  IF ( sum( parame(1:ncomp,6) ) /= 0.0 ) dipole = 1
925  IF ( sum( parame(1:ncomp,7) ) /= 0.0 ) quadrupole = 1
926  IF ( dipole == 1 .AND. quadrupole == 1 ) dipole_quad = 1
927 
928  !-----------------------------------------------------------------------------
929  ! dipole-dipole term
930  !-----------------------------------------------------------------------------
931  IF (dipole == 1) THEN
932 
933  IF (dd_term == 'GV') CALL phi_dd_gross_vrabec( k, z3_rk, fdd_rk )
934  ! IF (dd_term == 'SF') CALL PHI_DD_SAAGER_FISCHER( k )
935 
936  IF (dd_term /= 'GV' .AND. dd_term /= 'SF') write (*,*) 'specify dipole term !'
937 
938  ENDIF
939 
940  !-----------------------------------------------------------------------------
941  ! quadrupole-quadrupole term
942  !-----------------------------------------------------------------------------
943  IF (quadrupole == 1) THEN
944 
945  !IF (qq_term == 'SF') CALL PHI_QQ_SAAGER_FISCHER( k )
946  IF (qq_term == 'JG') CALL phi_qq_gross( k, z3_rk, fqq_rk )
947 
948  IF (qq_term /= 'JG' .AND. qq_term /= 'SF') write (*,*) 'specify quadrupole term !'
949 
950  ENDIF
951 
952  !-----------------------------------------------------------------------------
953  ! dipole-quadrupole cross term
954  !-----------------------------------------------------------------------------
955  IF (dipole_quad == 1) THEN
956 
957  IF (dq_term == 'VG') CALL phi_dq_vrabec_gross( k, z3_rk, fdq_rk )
958 
959  IF (dq_term /= 'VG' ) write (*,*) 'specify DQ-cross term !'
960 
961  ENDIF
962 
963 END SUBROUTINE phi_polar
964 
965 
966 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
967 !
968 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
969 
970 SUBROUTINE phi_dd_gross_vrabec( k, z3_rk, fdd_rk )
971 
972  USE parameters, ONLY: pi, kbol
973  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
974  USE eos_constants, ONLY: ddp2, ddp3, ddp4
975 
976  !-----------------------------------------------------------------------------
977  INTEGER, INTENT(IN) :: k
978  REAL, INTENT(IN) :: z3_rk
979  REAL, INTENT(IN OUT) :: fdd_rk
980 
981  !-----------------------------------------------------------------------------
982  INTEGER :: i, j, l, m
983 
984  REAL :: factor2, factor3, z3
985  REAL :: xijfa, xijkfa, xijfa_x, xijkf_x, eij
986  REAL :: fdd2, fdd3, fdd2x, fdd3x
987  REAL, DIMENSION(nc) :: my2dd
988  REAL, DIMENSION(nc,nc) :: idd2, idd4, idd2x, idd4x
989  REAL, DIMENSION(nc,nc,nc) :: idd3, idd3x
990  !-----------------------------------------------------------------------------
991 
992 
993  fdd_rk = 0.0
994  z3 = eta
995  DO i = 1, ncomp
996  IF ( uij(i,i) == 0.0 ) write (*,*) 'PHI_DD_GROSS_VRABEC: do not use dimensionless units'
997  IF ( uij(i,i) == 0.0 ) stop
998  my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
999  END DO
1000 
1001  DO i = 1, ncomp
1002  DO j = 1, ncomp
1003  idd2(i,j) = 0.0
1004  idd4(i,j) = 0.0
1005  idd2x(i,j) = 0.0
1006  idd4x(i,j) = 0.0
1007  IF (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0) THEN
1008  DO m=0,4
1009  idd2(i,j) =idd2(i,j) + ddp2(i,j,m)*z3**m
1010  idd4(i,j) =idd4(i,j) + ddp4(i,j,m)*z3**m
1011  idd2x(i,j) =idd2x(i,j)+ ddp2(i,j,m)*REAL(m)*z3**(m-1)*z3_rk
1012  idd4x(i,j) =idd4x(i,j)+ ddp4(i,j,m)*REAL(m)*z3**(m-1)*z3_rk
1013  END DO
1014  DO l = 1, ncomp
1015  idd3(i,j,l) = 0.0
1016  idd3x(i,j,l) = 0.0
1017  IF (parame(l,6) /= 0.0) THEN
1018  DO m=0,4
1019  idd3(i,j,l) =idd3(i,j,l) +ddp3(i,j,l,m)*z3**m
1020  idd3x(i,j,l)=idd3x(i,j,l)+ddp3(i,j,l,m)*REAL(m)*z3**(m-1)*z3_rk
1021  END DO
1022  END IF
1023  END DO
1024  END IF
1025  END DO
1026  END DO
1027 
1028  factor2= -pi
1029  factor3= -4.0/3.0*pi**2
1030 
1031  fdd2 = 0.0
1032  fdd3 = 0.0
1033  fdd2x = 0.0
1034  fdd3x = 0.0
1035  DO i = 1, ncomp
1036  xijfa_x = 2.0*x(i)*rho*uij(i,i)*my2dd(i)*sig_ij(i,i)**3 /t &
1037  *uij(k,k)*my2dd(k)*sig_ij(k,k)**3 /t/sig_ij(i,k)**3
1038  eij = (parame(i,3)*parame(k,3))**0.5
1039  fdd2x = fdd2x + factor2*xijfa_x*( idd2(i,k) + eij/t*idd4(i,k) )
1040  DO j = 1, ncomp
1041  IF (parame(i,6) /= 0.0 .AND. parame(j,6) /= 0.0) THEN
1042  xijfa =x(i)*rho*uij(i,i)*my2dd(i)*sig_ij(i,i)**3 /t &
1043  *x(j)*rho*uij(j,j)*my2dd(j)*sig_ij(j,j)**3 /t/sig_ij(i,j)**3
1044  eij = (parame(i,3)*parame(j,3))**0.5
1045  fdd2= fdd2 +factor2*xijfa*(idd2(i,j) +eij/t*idd4(i,j) )
1046  fdd2x =fdd2x +factor2*xijfa*(idd2x(i,j)+eij/t*idd4x(i,j))
1047  !---------------------
1048  xijkf_x=x(i)*rho*uij(i,i)*my2dd(i)*sig_ij(i,i)**3 /t/sig_ij(i,j) &
1049  *x(j)*rho*uij(j,j)*my2dd(j)*sig_ij(j,j)**3 /t/sig_ij(i,k) &
1050  *3.0* uij(k,k)*my2dd(k)*sig_ij(k,k)**3 /t/sig_ij(j,k)
1051  fdd3x=fdd3x+factor3*xijkf_x*idd3(i,j,k)
1052  DO l=1,ncomp
1053  IF (parame(l,6) /= 0.0) THEN
1054  xijkfa= x(i)*rho*uij(i,i)/t*my2dd(i)*sig_ij(i,i)**3 &
1055  *x(j)*rho*uij(j,j)/t*my2dd(j)*sig_ij(j,j)**3 &
1056  *x(l)*rho*uij(l,l)/t*my2dd(l)*sig_ij(l,l)**3 &
1057  /sig_ij(i,j)/sig_ij(i,l)/sig_ij(j,l)
1058  fdd3 =fdd3 + factor3 * xijkfa *idd3(i,j,l)
1059  fdd3x =fdd3x + factor3 * xijkfa *idd3x(i,j,l)
1060  END IF
1061  END DO
1062  END IF
1063  END DO
1064  END DO
1065 
1066  IF (fdd2 < -1.e-50 .AND. fdd3 /= 0.0 .AND. fdd2x /= 0.0 .AND. fdd3x /= 0.0)THEN
1067 
1068  fdd_rk = fdd2* (fdd2*fdd2x - 2.0*fdd3*fdd2x+fdd2*fdd3x) / (fdd2-fdd3)**2
1069 
1070  END IF
1071 
1072 END SUBROUTINE phi_dd_gross_vrabec
1073 
1074 
1075 
1076 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1077 !
1078 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1079 
1080 SUBROUTINE phi_qq_gross( k, z3_rk, fqq_rk )
1081 
1082  USE parameters, ONLY: pi, kbol
1083  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
1084  USE eos_constants, ONLY: qqp2, qqp3, qqp4
1085 
1086  !-----------------------------------------------------------------------------
1087  INTEGER, INTENT(IN) :: k
1088  REAL, INTENT(IN) :: z3_rk
1089  REAL, INTENT(IN OUT) :: fqq_rk
1090 
1091  !-----------------------------------------------------------------------------
1092  INTEGER :: i, j, l, m
1093 
1094  REAL :: factor2, factor3, z3
1095  REAL :: xijfa, xijkfa, xijfa_x, xijkf_x, eij
1096  REAL :: fqq2, fqq3, fqq2x, fqq3x
1097  REAL, DIMENSION(nc) :: qq2
1098  REAL, DIMENSION(nc,nc) :: iqq2, iqq4, iqq2x, iqq4x
1099  REAL, DIMENSION(nc,nc,nc) :: iqq3, iqq3x
1100  !-----------------------------------------------------------------------------
1101 
1102 
1103  fqq_rk = 0.0
1104  z3 = eta
1105  DO i = 1, ncomp
1106  IF ( uij(i,i) == 0.0 ) write (*,*) 'PHI_QQ_GROSS: do not use dimensionless units'
1107  qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
1108  END DO
1109 
1110  DO i = 1, ncomp
1111  DO j = 1, ncomp
1112  iqq2(i,j) = 0.0
1113  iqq4(i,j) = 0.0
1114  iqq2x(i,j) = 0.0
1115  iqq4x(i,j) = 0.0
1116  IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) THEN
1117  DO m = 0, 4
1118  iqq2(i,j) = iqq2(i,j) + qqp2(i,j,m) * z3**m
1119  iqq4(i,j) = iqq4(i,j) + qqp4(i,j,m) * z3**m
1120  iqq2x(i,j) = iqq2x(i,j) + qqp2(i,j,m) * REAL(m)*z3**(m-1)*z3_rk
1121  iqq4x(i,j) = iqq4x(i,j) + qqp4(i,j,m) * REAL(m)*z3**(m-1)*z3_rk
1122  END DO
1123  DO l = 1, ncomp
1124  iqq3(i,j,l) = 0.0
1125  iqq3x(i,j,l) = 0.0
1126  IF (parame(l,7) /= 0.0) THEN
1127  DO m = 0, 4
1128  iqq3(i,j,l) = iqq3(i,j,l) + qqp3(i,j,l,m)*z3**m
1129  iqq3x(i,j,l) = iqq3x(i,j,l) + qqp3(i,j,l,m)*REAL(m) *z3**(m-1)*z3_rk
1130  END DO
1131  END IF
1132  END DO
1133  END IF
1134  END DO
1135  END DO
1136 
1137  factor2= -9.0/16.0*pi
1138  factor3= 9.0/16.0*pi**2
1139 
1140  fqq2 = 0.0
1141  fqq3 = 0.0
1142  fqq2x = 0.0
1143  fqq3x = 0.0
1144  DO i = 1, ncomp
1145  xijfa_x = 2.0*x(i)*rho*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t &
1146  *uij(k,k)*qq2(k)*sig_ij(k,k)**5 /t/sig_ij(i,k)**7
1147  eij = (parame(i,3)*parame(k,3))**0.5
1148  fqq2x =fqq2x +factor2*xijfa_x*(iqq2(i,k)+eij/t*iqq4(i,k))
1149  DO j=1,ncomp
1150  IF (parame(i,7) /= 0.0 .AND. parame(j,7) /= 0.0) THEN
1151  xijfa =x(i)*rho*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t &
1152  *x(j)*rho*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,j)**7
1153  eij = (parame(i,3)*parame(j,3))**0.5
1154  fqq2= fqq2 +factor2*xijfa*(iqq2(i,j) +eij/t*iqq4(i,j) )
1155  fqq2x =fqq2x +factor2*xijfa*(iqq2x(i,j)+eij/t*iqq4x(i,j))
1156  ! ------------------
1157  xijkf_x=x(i)*rho*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t/sig_ij(i,j)**3 &
1158  *x(j)*rho*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,k)**3 &
1159  *3.0* uij(k,k)*qq2(k)*sig_ij(k,k)**5 /t/sig_ij(j,k)**3
1160  fqq3x = fqq3x + factor3*xijkf_x*iqq3(i,j,k)
1161  DO l = 1, ncomp
1162  IF (parame(l,7) /= 0.0) THEN
1163  xijkfa=x(i)*rho*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /t/sig_ij(i,j)**3 &
1164  *x(j)*rho*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /t/sig_ij(i,l)**3 &
1165  *x(l)*rho*uij(l,l)*qq2(l)*sig_ij(l,l)**5 /t/sig_ij(j,l)**3
1166  fqq3 =fqq3 + factor3 * xijkfa *iqq3(i,j,l)
1167  fqq3x =fqq3x + factor3 * xijkfa *iqq3x(i,j,l)
1168  END IF
1169  END DO
1170  END IF
1171  END DO
1172  END DO
1173 
1174  IF (fqq2 < -1.e-50 .AND. fqq3 /= 0.0 .AND. fqq2x /= 0.0 .AND. fqq3x /= 0.0) THEN
1175  fqq_rk = fqq2* (fqq2*fqq2x - 2.0*fqq3*fqq2x+fqq2*fqq3x) / (fqq2-fqq3)**2
1176  END IF
1177 
1178 END SUBROUTINE phi_qq_gross
1179 
1180 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1181 !
1182 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1183 
1184 SUBROUTINE phi_dq_vrabec_gross( k, z3_rk, fdq_rk )
1185 
1186  USE parameters, ONLY: pi, kbol
1187  USE eos_variables, ONLY: nc, ncomp, uij, parame, mseg, sig_ij, rho, eta, x, t
1188  USE eos_constants, ONLY: dqp2, dqp3, dqp4
1189 
1190  !-----------------------------------------------------------------------------
1191  INTEGER, INTENT(IN) :: k
1192  REAL, INTENT(IN) :: z3_rk
1193  REAL, INTENT(IN OUT) :: fdq_rk
1194 
1195  !-----------------------------------------------------------------------------
1196  INTEGER :: i, j, l, m
1197 
1198  REAL :: factor2, factor3, z3
1199  REAL :: xijfa, xijkfa, xijfa_x, xijkf_x, eij
1200  REAL :: fdq2, fdq3, fdq2x, fdq3x
1201  REAL, DIMENSION(nc) :: my2dd, myfac, qq2, q_fac
1202  REAL, DIMENSION(nc,nc) :: idq2, idq4, idq2x, idq4x
1203  REAL, DIMENSION(nc,nc,nc) :: idq3, idq3x
1204  !-----------------------------------------------------------------------------
1205 
1206  fdq_rk = 0.0
1207  z3 = eta
1208  DO i=1,ncomp
1209  my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
1210  myfac(i) = parame(i,3)/t*parame(i,2)**4 *my2dd(i)
1211  qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
1212  q_fac(i) = parame(i,3)/t*parame(i,2)**4 *qq2(i)
1213  END DO
1214 
1215  DO i = 1, ncomp
1216  DO j = 1, ncomp
1217  idq2(i,j) = 0.0
1218  idq4(i,j) = 0.0
1219  idq2x(i,j) = 0.0
1220  idq4x(i,j) = 0.0
1221  DO m = 0, 4
1222  idq2(i,j) = idq2(i,j) + dqp2(i,j,m)*z3**m
1223  idq4(i,j) = idq4(i,j) + dqp4(i,j,m)*z3**m
1224  idq2x(i,j) = idq2x(i,j) + dqp2(i,j,m)*REAL(m)*z3**(m-1) *z3_rk
1225  idq4x(i,j) = idq4x(i,j) + dqp4(i,j,m)*REAL(m)*z3**(m-1) *z3_rk
1226  END DO
1227  DO l = 1, ncomp
1228  idq3(i,j,l) = 0.0
1229  idq3x(i,j,l) = 0.0
1230  DO m = 0, 4
1231  idq3(i,j,l) =idq3(i,j,l) +dqp3(i,j,l,m)*z3**m
1232  idq3x(i,j,l)=idq3x(i,j,l)+dqp3(i,j,l,m)*REAL(m)*z3**(m-1)*z3_rk
1233  END DO
1234  END DO
1235  END DO
1236  END DO
1237 
1238  factor2= -9.0/4.0*pi
1239  factor3= pi**2
1240 
1241  fdq2 = 0.0
1242  fdq3 = 0.0
1243  fdq2x = 0.0
1244  fdq3x = 0.0
1245  DO i = 1, ncomp
1246  xijfa_x = x(i)*rho*( myfac(i)*q_fac(k) + myfac(k)*q_fac(i) ) / sig_ij(i,k)**5
1247  eij = (parame(i,3)*parame(k,3))**0.5
1248  fdq2x =fdq2x +factor2*xijfa_x*(idq2(i,k)+eij/t*idq4(i,k))
1249  DO j=1,ncomp
1250  xijfa =x(i)*rho*myfac(i) * x(j)*rho*q_fac(j) /sig_ij(i,j)**5
1251  eij = (parame(i,3)*parame(j,3))**0.5
1252  fdq2= fdq2 +factor2*xijfa*(idq2(i,j) +eij/t*idq4(i,j) )
1253  fdq2x =fdq2x +factor2*xijfa*(idq2x(i,j) +eij/t*idq4x(i,j))
1254  !---------------------
1255  xijkf_x=x(i)*rho*x(j)*rho/(sig_ij(i,j)*sig_ij(i,k)*sig_ij(j,k))**2 &
1256  *( myfac(i)*q_fac(j)*myfac(k) + myfac(i)*q_fac(k)*myfac(j) &
1257  + myfac(k)*q_fac(i)*myfac(j) +myfac(i)*q_fac(j)*q_fac(k)*1.1937350 &
1258  +myfac(i)*q_fac(k)*q_fac(j)*1.193735 &
1259  +myfac(k)*q_fac(i)*q_fac(j)*1.193735 )
1260  fdq3x = fdq3x + factor3*xijkf_x*idq3(i,j,k)
1261  DO l = 1, ncomp
1262  xijkfa=x(i)*rho*x(j)*rho*x(l)*rho/(sig_ij(i,j)*sig_ij(i,l)*sig_ij(j,l))**2 &
1263  *( myfac(i)*q_fac(j)*myfac(l) &
1264  +myfac(i)*q_fac(j)*q_fac(l)*1.193735 )
1265  fdq3 =fdq3 + factor3 * xijkfa *idq3(i,j,l)
1266  fdq3x =fdq3x + factor3 * xijkfa *idq3x(i,j,l)
1267  END DO
1268  END DO
1269  END DO
1270 
1271  IF (fdq2 < -1.e-50 .AND. fdq3 /= 0.0 .AND. fdq2x /= 0.0 .AND. fdq3x /= 0.0)THEN
1272 
1273  fdq_rk = fdq2* (fdq2*fdq2x - 2.0*fdq3*fdq2x+fdq2*fdq3x) / (fdq2-fdq3)**2
1274 
1275  END IF
1276 
1277 END SUBROUTINE phi_dq_vrabec_gross
1278 
1279 END MODULE eos_polar
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
eos_constants
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:200
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120