MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
eos_enthalpy.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 !
3 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
4 
5 subroutine enthalpy_ig_qspr ( cpig, hig, sig )
6 
7  use parameters, only: rgas
8  use basic_variables, only: nc, np, ncomp, nphas, t, p, xi, parame
9  implicit none
10 
11  !-----------------------------------------------------------------------------
12  real, dimension(np), intent(out) :: cpig
13  real, dimension(np), intent(out) :: hig
14  real, dimension(np), intent(out) :: sig
15 
16  !-----------------------------------------------------------------------------
17  integer :: i, ph
18  real :: p0
19  real :: t0
20  real, dimension(nc) :: cp_ig, h_ig, s_ig
21  real :: m_par, sig_par, eps_par
22  real :: dip_par, quad_par, eps_ab_par, kap_ab_par
23  real :: p1, p2, p3, p4, p5, p6
24  real :: c1a, c2a, c3a, c4a, c5a, c6a
25  real :: c1b, c2b, c3b, c4b, c5b, c6b
26  real :: a_coeff, b_coeff
27  real :: icpc300, icpc400
28  real :: t1, t2
29  !-----------------------------------------------------------------------------
30 
31  t0 = 298.0
32 
33  do i = 1, ncomp
34 
35  m_par = parame(i, 1)
36  sig_par = parame(i, 2)
37  eps_par = parame(i, 3)
38 
39  eps_ab_par = parame(i,14)
40  kap_ab_par = parame(i,13)
41  dip_par = parame(i, 6)
42  quad_par = parame(i, 7)
43 
44  p1 = m_par * eps_par / t
45  p2 = m_par * sig_par**3
46  p3 = m_par * sig_par**3 * eps_par / t
47  p4 = m_par * sig_par**6 * kap_ab_par *( exp( eps_ab_par/t ) - 1.0 )
48  p5 = m_par * sig_par**3 * quad_par
49  p6 = 1.0
50 
51  t1 = 300.0
52  t2 = 400.0
53 
54  !--------------------------------------------------------------------------
55  ! nAnP
56  !--------------------------------------------------------------------------
57 
58  ! ICPC @ 300 K
59 
60  c1a = -5763.04893
61  c2a = 1232.30607
62  c3a = -239.3513996
63  c4a = 0.0
64  c5a = 0.0
65  c6a = -15174.28321
66 
67  ! ICPC @ 400 K
68 
69  c1b = -8171.26676935062
70  c2b = 1498.01217504596
71  c3b = -315.515836223387
72  c4b = 0.0
73  c5b = 0.0
74  c6b = -19389.5468655708
75 
76  !--------------------------------------------------------------------------
77  ! nAP
78  !--------------------------------------------------------------------------
79 
80  !
81  !! ICPC @ 300K:
82  ! c1a = 5177.19095226181
83  ! c2a = 919.565206504576
84  ! c3a = -108.829105648889
85  ! c4a = 0
86  ! c5a = -3.93917830677682
87  ! c6a = -13504.5671858292
88  !
89  !! ICPC @ 400K:
90  ! c1b = 10656.1018362315
91  ! c2b = 1146.10782703748
92  ! c3b = -131.023645998081
93  ! c4b = 0
94  ! c5b = -9.93789225413177
95  ! c6b = -24430.12952497
96 
97  !--------------------------------------------------------------------------
98  ! AP
99  !--------------------------------------------------------------------------
100 
101  !! ICPC @ 300K:
102  ! c1a = 3600.32322462175
103  ! c2a = 1006.20461224949
104  ! c3a = -151.688378113974
105  ! c4a = 7.81876773647109d-07
106  ! c5a = 8.01001754473385
107  ! c6a = -8959.37140957179
108 
109  !! ICPC @ 400K:
110  ! c1b = 7248.0697641199
111  ! c2b = 1267.44346171358
112  ! c3b = -208.738557800023
113  ! c4b = 0.000170238690157906
114  ! c5b = -6.7841792685616
115  ! c6b = -12669.4196622924
116 
117  icpc300 = c1a*p1/300.0 + c2a*p2 + c3a*p3/300.0 + c4a*p4/300.0 + c5a*p5 + c6a*p6
118 
119  icpc300 = icpc300 /1000.0
120 
121  icpc400 = c1b*p1/400.0 + c2b*p2 + c3b*p3/400.0 + c4b*p4/400.0 + c5b*p5 + c6b*p6
122 
123  icpc400 = icpc400 / 1000.0
124 
125  !--------------------------------------------------------------------------
126  ! linear approximation
127  !--------------------------------------------------------------------------
128 
129  b_coeff = ( icpc400 - icpc300 ) / ( t2 - t1 )
130  a_coeff = icpc300 - b_coeff * t1
131 
132  t0 = 298.0
133  p0 = 1.e5
134 
135  cp_ig(i)= a_coeff + b_coeff * t
136  h_ig(i) = a_coeff * (t-t0) + 0.5 * b_coeff * ( t**2 - t0**2 )
137  s_ig(i) = a_coeff * log( t/t0 ) + b_coeff * (t-t0) - rgas * log( p/p0 )
138 
139  end do
140 
141  !-----------------------------------------------------------------------------
142  ! Ideal gas properties
143  !-----------------------------------------------------------------------------
144  do ph = 1, nphas
145 
146  hig(ph) = sum( xi( ph, 1:ncomp ) * h_ig( 1:ncomp ) )
147  cpig(ph) = sum( xi( ph, 1:ncomp ) *cp_ig( 1:ncomp ) )
148 
149  sig(ph) = 0.0
150  do i = 1, ncomp
151  if ( xi(ph,i) > 0.0 ) sig(ph) = sig(ph) + xi(ph,i) * ( s_ig(i) - rgas * log( xi(ph,i) ) )
152  end do
153 
154  end do
155 
156 end subroutine enthalpy_ig_qspr
157 
158 
159 
160 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
161 !
162 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
163 
164 subroutine enthalpy_ig ( coef_ig, cpig, hig, sig )
165 
166  use parameters, only: rgas
167  use basic_variables, only: nc, np, ncomp, nphas, t, p, xi
168  implicit none
169 
170  !-----------------------------------------------------------------------------
171  real, dimension(nc,5), intent(in) :: coef_ig
172  real, dimension(np), intent(out) :: cpig
173  real, dimension(np), intent(out) :: hig
174  real, dimension(np), intent(out) :: sig
175 
176  !-----------------------------------------------------------------------------
177  integer :: i, ph
178  real :: p0, t0, t02, t03, t04
179  real, dimension(nc) :: cp_ig, h_ig, s_ig
180  !-----------------------------------------------------------------------------
181 
182  p0 = 1.e5
183 
184  t0 = 298.15
185  t02 = t0 * t0
186  t03 = t02 * t0
187  t04 = t02 * t02
188  do i = 1, ncomp
189 
190  cp_ig(i)= coef_ig(i,1) + coef_ig(i,2)*t + coef_ig(i,3)*t*t + coef_ig(i,4)*t**3
191  h_ig(i) = coef_ig(i,1)*(t-t0) + coef_ig(i,2)/2.0*(t*t-t02) &
192  + coef_ig(i,3)/3.0*(t**3-t03) + coef_ig(i,4)/4.0*(t**4 - t04)
193  s_ig(i) = coef_ig(i,1) * log( t/t0 ) + coef_ig(i,2) * (t-t0) &
194  + coef_ig(i,3)/2.0*(t*t-t02) + coef_ig(i,4)/3.0*(t**3-t03) - rgas * log( p/p0 )
195 
196  end do
197 
198  do ph = 1, nphas
199 
200  hig(ph) = sum( xi( ph, 1:ncomp ) * h_ig( 1:ncomp ) )
201  cpig(ph) = sum( xi( ph, 1:ncomp ) *cp_ig( 1:ncomp ) )
202  sig(ph) = 0.0
203  do i = 1, ncomp
204  if ( xi(ph,i) > 0.0 ) sig(ph) = sig(ph) + xi(ph,i) * ( s_ig(i) - rgas * log( xi(ph,i) ) )
205  end do
206 
207  end do
208 
209 END SUBROUTINE enthalpy_ig
210 
211 
212 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
213 !
214 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
215 
216 SUBROUTINE h_eos
217 
218  USE parameters, ONLY: rgas
219  ! USE BASIC_VARIABLES, ONLY: mm ! only for speed of sound
220  USE eos_constants
221  USE eos_variables
222  IMPLICIT NONE
223 
224  !-----------------------------------------------------------------------------
225  REAL :: zges, df_dt, dfdr, ddfdrdr
226  REAL :: cv_res, df_dt2, df_drdt
227  REAL :: fact, dist, t_tmp, rho_0
228  REAL :: fdr1, fdr2, fdr3, fdr4
229  ! REAL :: cp_ig_mix, w_square, mm_mean ! only for speed of sound
230 
231 
232  INTEGER :: i, m
233  REAL :: dhsdt(nc), dhsdt2(nc)
234  REAL :: z0, z1, z2, z3, z1tdt, z2tdt, z3tdt
235  REAL :: z1dt, z2dt, z3dt, zms, gii
236  REAL :: fhsdt, fhsdt2
237  REAL :: fchdt, fchdt2
238  REAL :: fdspdt, fdspdt2
239  REAL :: fhbdt, fhbdt2
240  REAL :: sumseg, i1, i2, i1dt, i2dt, i1dt2, i2dt2
241  REAL :: c1_con, c2_con, c3_con, c1_dt, c1_dt2
242  REAL :: z1tdt2, z2tdt2, z3tdt2
243  REAL :: z1dt2, z2dt2, z3dt2
244 
245  INTEGER :: j, k, l, no, ass_cnt, max_eval
246  LOGICAL :: assoc
247  REAL :: dij, dijdt, dijdt2
248  REAL :: gij1dt, gij2dt, gij3dt
249  REAL :: gij1dt2, gij2dt2, gij3dt2
250  REAL, DIMENSION(nc,nc) :: gijdt, gijdt2, kap_hb
251  REAL, DIMENSION(nc,nc,nsite,nsite) :: ass_d_dt, ass_d_dt2, eps_hb, delta, deltadt, deltadt2
252  REAL, DIMENSION(nc,nsite) :: mxdt, mxdt2, mx_itr, mx_itrdt, mx_itrdt2
253  REAL :: attenu, tol, suma, sumdt, sumdt2, err_sum
254 
255  INTEGER :: dipole
256  REAL :: fdddt, fdddt2
257  REAL, DIMENSION(nc) :: my2dd, my0
258  REAL, DIMENSION(nc,nc) :: idd2, idd2dt, idd2dt2, idd4, idd4dt, idd4dt2
259  REAL, DIMENSION(nc,nc,nc) :: idd3, idd3dt, idd3dt2
260  REAL :: factor2, factor3
261  REAL :: fdd2, fdd3, fdd2dt, fdd3dt, fdd2dt2, fdd3dt2
262  REAL :: eij, xijmt, xijkmt
263 
264  INTEGER :: qudpole
265  REAL :: fqqdt, fqqdt2
266  REAL, DIMENSION(nc) :: qq2
267  REAL, DIMENSION(nc,nc) :: iqq2, iqq2dt, iqq2dt2, iqq4, iqq4dt, iqq4dt2
268  REAL, DIMENSION(nc,nc,nc) :: iqq3, iqq3dt, iqq3dt2
269  REAL :: fqq2, fqq2dt, fqq2dt2, fqq3, fqq3dt, fqq3dt2
270 
271  INTEGER :: dip_quad
272  REAL :: fdqdt, fdqdt2
273  REAL, DIMENSION(nc) :: myfac, q_fac
274  REAL, DIMENSION(nc,nc) :: idq2, idq2dt, idq2dt2, idq4, idq4dt, idq4dt2
275  REAL, DIMENSION(nc,nc,nc) :: idq3, idq3dt, idq3dt2
276  REAL :: fdq2, fdq2dt, fdq2dt2, fdq3, fdq3dt, fdq3dt2
277  !-----------------------------------------------------------------------------
278 
279 
280  !-----------------------------------------------------------------------------
281  ! Initializing
282  !-----------------------------------------------------------------------------
283 
284  CALL perturbation_parameter
285 
286  rho = eta/z3t
287  z0 = z0t*rho
288  z1 = z1t*rho
289  z2 = z2t*rho
290  z3 = z3t*rho
291 
292  sumseg = z0t / (pi/6.0)
293  zms = 1.0 - z3
294 
295 
296  !-----------------------------------------------------------------------------
297  ! first and second derivative of f to density (dfdr,ddfdrdr)
298  !-----------------------------------------------------------------------------
299 
300  CALL p_eos
301 
302  zges = (pges * 1.e-30)/(kbol*t*rho)
303 
304  dfdr = pges/(eta*rho*(kbol*t)/1.e-30)
305  ddfdrdr = pgesdz/(eta*rho*(kbol*t)/1.e-30) - dfdr*2.0/eta - 1.0/eta**2
306 
307 
308  !-----------------------------------------------------------------------------
309  ! Helmholtz Energy f/kT = fres
310  !-----------------------------------------------------------------------------
311 
312  CALL f_eos
313 
314 
315  !-----------------------------------------------------------------------------
316  ! derivative of some auxilliary properties to temperature
317  !-----------------------------------------------------------------------------
318 
319  DO i = 1, ncomp
320  dhsdt(i)=parame(i,2) *(-3.0*parame(i,3)/t/t)*0.12*exp(-3.0*parame(i,3)/t)
321  dhsdt2(i) = dhsdt(i)*3.0*parame(i,3)/t/t &
322  + 6.0*parame(i,2)*parame(i,3)/t**3 *0.12*exp(-3.0*parame(i,3)/t)
323  END DO
324 
325  z1tdt = 0.0
326  z2tdt = 0.0
327  z3tdt = 0.0
328  DO i = 1, ncomp
329  z1tdt = z1tdt + x(i) * mseg(i) * dhsdt(i)
330  z2tdt = z2tdt + x(i) * mseg(i) * 2.0*dhs(i)*dhsdt(i)
331  z3tdt = z3tdt + x(i) * mseg(i) * 3.0*dhs(i)*dhs(i)*dhsdt(i)
332  END DO
333  z1dt = pi / 6.0*z1tdt *rho
334  z2dt = pi / 6.0*z2tdt *rho
335  z3dt = pi / 6.0*z3tdt *rho
336 
337 
338  z1tdt2 = 0.0
339  z2tdt2 = 0.0
340  z3tdt2 = 0.0
341  DO i = 1, ncomp
342  z1tdt2 = z1tdt2 + x(i)*mseg(i)*dhsdt2(i)
343  z2tdt2 = z2tdt2 + x(i)*mseg(i)*2.0 *( dhsdt(i)*dhsdt(i) +dhs(i)*dhsdt2(i) )
344  z3tdt2 = z3tdt2 + x(i)*mseg(i)*3.0 *( 2.0*dhs(i)*dhsdt(i)* &
345  dhsdt(i) +dhs(i)*dhs(i)*dhsdt2(i) )
346  END DO
347  z1dt2 = pi / 6.0*z1tdt2 *rho
348  z2dt2 = pi / 6.0*z2tdt2 *rho
349  z3dt2 = pi / 6.0*z3tdt2 *rho
350 
351 
352  !-----------------------------------------------------------------------------
353  ! 1st & 2nd derivative of f/kT hard spheres to temp. (fhsdt)
354  !-----------------------------------------------------------------------------
355 
356  fhsdt = 6.0/pi/rho*( 3.0*(z1dt*z2+z1*z2dt)/zms + 3.0*z1*z2*z3dt/zms/zms &
357  + 3.0*z2*z2*z2dt/z3/zms/zms &
358  + z2**3 *(2.0*z3*z3dt-z3dt*zms)/(z3*z3*zms**3 ) &
359  + (3.0*z2*z2*z2dt*z3-2.0*z2**3 *z3dt)/z3**3 *log(zms) &
360  + (z0-z2**3 /z3/z3)*z3dt/zms )
361 
362  fhsdt2= 6.0/pi/rho*( 3.0*(z1dt2*z2+2.0*z1dt*z2dt+z1*z2dt2)/zms &
363  + 6.0*(z1dt*z2+z1*z2dt)*z3dt/zms/zms &
364  + 3.0*z1*z2*z3dt2/zms/zms + 6.0*z1*z2*z3dt*z3dt/zms**3 &
365  + 3.0*z2*(2.0*z2dt*z2dt+z2*z2dt2)/z3/zms/zms &
366  - z2*z2*(6.0*z2dt*z3dt+z2*z3dt2)/(z3*z3*zms*zms) &
367  + 2.0*z2**3 *z3dt*z3dt/(z3**3 *zms*zms) &
368  - 4.0*z2**3 *z3dt*z3dt/(z3*z3 *zms**3 ) &
369  + (12.0*z2*z2*z2dt*z3dt+2.0*z2**3 *z3dt2)/(z3*zms**3 ) &
370  + 6.0*z2**3 *z3dt*z3dt/(z3*zms**4 ) &
371  - 2.0*(3.0*z2*z2*z2dt/z3/z3-2.0*z2**3 *z3dt/z3**3 ) *z3dt/zms &
372  -(z2**3 /z3/z3-z0)*(z3dt2/zms+z3dt*z3dt/zms/zms) &
373  + ( (6.0*z2*z2dt*z2dt+3.0*z2*z2*z2dt2)/z3/z3 &
374  - (12.0*z2*z2*z2dt*z3dt+2.0*z2**3 *z3dt2)/z3**3 &
375  + 6.0*z2**3 *z3dt*z3dt/z3**4 )* log(zms) )
376 
377 
378  !-----------------------------------------------------------------------------
379  ! 1st & 2nd derivative of f/kT of chain term to T (fchdt)
380  !-----------------------------------------------------------------------------
381 
382  fchdt = 0.0
383  fchdt2 = 0.0
384  DO i = 1, ncomp
385  DO j = 1, ncomp
386  dij=dhs(i)*dhs(j)/(dhs(i)+dhs(j))
387  dijdt =(dhsdt(i)*dhs(j) + dhs(i)*dhsdt(j)) / (dhs(i)+dhs(j)) &
388  - dhs(i)*dhs(j)/(dhs(i)+dhs(j))**2 *(dhsdt(i)+dhsdt(j))
389  dijdt2=(dhsdt2(i)*dhs(j) + 2.0*dhsdt(i)*dhsdt(j) &
390  + dhs(i)*dhsdt2(j)) / (dhs(i)+dhs(j)) &
391  - 2.0*(dhsdt(i)*dhs(j) + dhs(i)*dhsdt(j)) &
392  / (dhs(i)+dhs(j))**2 *(dhsdt(i)+dhsdt(j)) &
393  + 2.0* dhs(i)*dhs(j) / (dhs(i)+dhs(j))**3 &
394  * (dhsdt(i)+dhsdt(j))**2 &
395  - dhs(i)*dhs(j)/(dhs(i)+dhs(j))**2 *(dhsdt2(i)+dhsdt2(j))
396  gij1dt = z3dt/zms/zms
397  gij2dt = 3.0*( z2dt*dij+z2*dijdt )/zms/zms +6.0*z2*dij*z3dt/zms**3
398  gij3dt = 4.0*(dij*z2)* (dijdt*z2 + dij*z2dt) /zms**3 &
399  + 6.0*(dij*z2)**2 * z3dt /zms**4
400  gij1dt2 = z3dt2/zms/zms +2.0*z3dt*z3dt/zms**3
401  gij2dt2 = 3.0*( z2dt2*dij+2.0*z2dt*dijdt+z2*dijdt2 )/zms/zms &
402  + 6.0*( z2dt*dij+z2*dijdt )/zms**3 * z3dt &
403  + 6.0*(z2dt*dij*z3dt+z2*dijdt*z3dt+z2*dij*z3dt2) /zms**3 &
404  + 18.0*z2*dij*z3dt*z3dt/zms**4
405  gij3dt2 = 4.0*(dijdt*z2+dij*z2dt)**2 /zms**3 &
406  + 4.0*(dij*z2)* (dijdt2*z2+2.0*dijdt*z2dt+dij*z2dt2) /zms**3 &
407  + 24.0*(dij*z2) *(dijdt*z2+dij*z2dt)/zms**4 *z3dt &
408  + 6.0*(dij*z2)**2 * z3dt2 /zms**4 &
409  + 24.0*(dij*z2)**2 * z3dt*z3dt /zms**5
410  gijdt(i,j) = gij1dt + gij2dt + gij3dt
411  gijdt2(i,j) = gij1dt2 + gij2dt2 + gij3dt2
412  END DO
413  END DO
414 
415  DO i = 1, ncomp
416  gii = 1.0/zms + 3.0*dhs(i)/2.0*z2/zms/zms + 2.0*dhs(i)*dhs(i)/4.0*z2*z2/zms**3
417  fchdt = fchdt + x(i) * (1.0-mseg(i)) * gijdt(i,i) / gii
418  fchdt2= fchdt2+ x(i) * (1.0-mseg(i)) &
419  * (gijdt2(i,i) / gii - (gijdt(i,i)/gii)**2 )
420  END DO
421 
422 
423  !-----------------------------------------------------------------------------
424  ! 1st & 2nd derivative of f/kT dispersion term to T (fdspdt)
425  !-----------------------------------------------------------------------------
426 
427  i1 = 0.0
428  i2 = 0.0
429  i1dt = 0.0
430  i2dt = 0.0
431  i1dt2= 0.0
432  i2dt2= 0.0
433  DO m = 0, 6
434  i1 = i1 + apar(m)*z3**m
435  i2 = i2 + bpar(m)*z3**m
436  i1dt = i1dt + apar(m)*z3dt*REAL(m)*z3**(m-1)
437  i2dt = i2dt + bpar(m)*z3dt*REAL(m)*z3**(m-1)
438  i1dt2= i1dt2+ apar(m)*z3dt2*REAL(m)*z3**(m-1) &
439  + apar(m)*z3dt*z3dt *REAL(m)*REAL(m-1)*z3**(m-2)
440  i2dt2= i2dt2+ bpar(m)*z3dt2*REAL(m)*z3**(m-1) &
441  + bpar(m)*z3dt*z3dt *REAL(m)*REAL(m-1)*z3**(m-2)
442  END DO
443 
444  c1_con= 1.0/ ( 1.0 + sumseg*(8.0*z3-2.0*z3**2 )/zms**4 &
445  + (1.0 - sumseg)*(20.0*z3-27.0*z3**2 +12.0*z3**3 -2.0*z3**4 ) &
446  /(zms*(2.0-z3))**2 )
447  c2_con= - c1_con*c1_con *(sumseg*(-4.0*z3**2 +20.0*z3+8.0)/zms**5 &
448  + (1.0 - sumseg) *(2.0*z3**3 +12.0*z3**2 -48.0*z3+40.0) &
449  /(zms*(2.0-z3))**3 )
450  c3_con= 2.0 * c2_con*c2_con/c1_con - c1_con*c1_con &
451  *( sumseg*(-12.0*z3**2 +72.0*z3+60.0)/zms**6 + (1.0 - sumseg) &
452  *(-6.0*z3**4 -48.0*z3**3 +288.0*z3**2 -480.0*z3+264.0) &
453  /(zms*(2.0-z3))**4 )
454  c1_dt = c2_con*z3dt
455  c1_dt2 = c3_con*z3dt*z3dt + c2_con*z3dt2
456 
457  fdspdt = - 2.0*pi*rho*(i1dt-i1/t)*order1 &
458  - pi*rho*sumseg*(c1_dt*i2+c1_con*i2dt-2.0*c1_con*i2/t)*order2
459 
460  fdspdt2 = - 2.0*pi*rho*(i1dt2-2.0*i1dt/t+2.0*i1/t/t)*order1 &
461  - pi*rho*sumseg*order2*( c1_dt2*i2 +2.0*c1_dt*i2dt -4.0*c1_dt*i2/t &
462  + 6.0*c1_con*i2/t/t -4.0*c1_con*i2dt/t +c1_con*i2dt2)
463 
464 
465  !-----------------------------------------------------------------------------
466  ! 1st & 2nd derivative of f/kT association term to T (fhbdt)
467  !-----------------------------------------------------------------------------
468 
469  fhbdt = 0.0
470  fhbdt2 = 0.0
471  assoc = .false.
472  DO i = 1, ncomp
473  IF ( nhb_typ(i) /= 0 ) assoc = .true.
474  END DO
475  IF (assoc) THEN
476 
477  DO i = 1, ncomp
478  IF ( nhb_typ(i) /= 0 ) THEN
479  kap_hb(i,i) = parame(i,13)
480  no = 0
481  DO j = 1,nhb_typ(i)
482  DO k = 1,nhb_typ(i)
483  eps_hb(i,i,j,k) = parame(i,(14+no))
484  no = no + 1
485  END DO
486  END DO
487  DO j = 1,nhb_typ(i)
488  no = no + 1
489  END DO
490  ELSE
491  kap_hb(i,i) = 0.0
492  DO k = 1,nsite
493  DO l = 1,nsite
494  eps_hb(i,i,k,l) = 0.0
495  END DO
496  END DO
497  END IF
498  END DO
499 
500  DO i = 1, ncomp
501  DO j = 1, ncomp
502  IF ( i /= j .AND. (nhb_typ(i) /= 0 .AND. nhb_typ(j) /= 0) ) THEN
503  kap_hb(i,j)= (kap_hb(i,i)*kap_hb(j,j))**0.5 &
504  *((parame(i,2)*parame(j,2))**3 )**0.5 &
505  /(0.5*(parame(i,2)+parame(j,2)))**3
506  ! kap_hb(i,j)= kap_hb(i,j)*(1.0-k_kij)
507  DO k = 1,nhb_typ(i)
508  DO l = 1,nhb_typ(j)
509  IF (k /= l) THEN
510  eps_hb(i,j,k,l)=(eps_hb(i,i,k,l)+eps_hb(j,j,l,k))/2.0
511  ! eps_hb(i,j,k,l)=eps_hb(i,j,k,l)*(1.0-eps_kij)
512  END IF
513  END DO
514  END DO
515  END IF
516  END DO
517  END DO
518  IF (nhb_typ(1) == 3) THEN
519  eps_hb(1,2,3,1)=0.5*(eps_hb(1,1,3,1)+eps_hb(2,2,1,2))
520  eps_hb(2,1,1,3) = eps_hb(1,2,3,1)
521  END IF
522  IF (nhb_typ(2) == 3) THEN
523  eps_hb(2,1,3,1)=0.5*(eps_hb(2,2,3,1)+eps_hb(1,1,1,2))
524  eps_hb(1,2,1,3) = eps_hb(2,1,3,1)
525  END IF
526 
527  DO i = 1, ncomp
528  DO k = 1, nhb_typ(i)
529  DO j = 1, ncomp
530  DO l = 1, nhb_typ(j)
531  ! ass_d(i,j,k,l)=kap_hb(i,j) *sig_ij(i,j)**3 *(EXP(eps_hb(i,j,k,l)/t)-1.0)
532  ass_d_dt(i,j,k,l)= - kap_hb(i,j) *sig_ij(i,j)**3 * &
533  eps_hb(i,j,k,l)/t/t*exp(eps_hb(i,j,k,l)/t)
534  ass_d_dt2(i,j,k,l)= - kap_hb(i,j) *sig_ij(i,j)**3 &
535  * eps_hb(i,j,k,l)/t/t*exp(eps_hb(i,j,k,l)/t) &
536  * (-2.0/t - eps_hb(i,j,k,l)/t/t)
537  END DO
538  END DO
539  END DO
540  END DO
541 
542  DO i = 1, ncomp
543  DO k = 1, nhb_typ(i)
544  DO j = 1, ncomp
545  DO l = 1, nhb_typ(j)
546  delta(i,j,k,l)=gij(i,j)*ass_d(i,j,k,l)
547  deltadt(i,j,k,l) = gijdt(i,j)*ass_d(i,j,k,l) + gij(i,j)*ass_d_dt(i,j,k,l)
548  deltadt2(i,j,k,l)= gijdt2(i,j)*ass_d(i,j,k,l) &
549  + 2.0*gijdt(i,j)*ass_d_dt(i,j,k,l) +gij(i,j)*ass_d_dt2(i,j,k,l)
550  END DO
551  END DO
552  END DO
553  END DO
554 
555 
556  !--- constants for iteration ----------------------------------------------
557  attenu = 0.7
558  tol = 1.e-10
559  IF (eta < 0.2) tol = 1.e-11
560  IF (eta < 0.01) tol = 1.e-12
561  max_eval = 200
562 
563  !--- initialize mxdt(i,j) -------------------------------------------------
564  DO i = 1, ncomp
565  DO k = 1, nhb_typ(i)
566  mxdt(i,k) = 0.0
567  mxdt2(i,k) = 0.0
568  END DO
569  END DO
570 
571 
572  !--- iterate over all components and all sites ----------------------------
573  DO ass_cnt = 1, max_eval
574 
575  DO i = 1, ncomp
576  DO k = 1, nhb_typ(i)
577  suma = 0.0
578  sumdt = 0.0
579  sumdt2= 0.0
580  DO j = 1, ncomp
581  DO l = 1, nhb_typ(j)
582  suma = suma + x(j)*nhb_no(j,l)* mx(j,l) *delta(i,j,k,l)
583  sumdt = sumdt + x(j)*nhb_no(j,l)*( mx(j,l) *deltadt(i,j,k,l) &
584  + mxdt(j,l)*delta(i,j,k,l) )
585  sumdt2 = sumdt2 + x(j)*nhb_no(j,l)*( mx(j,l)*deltadt2(i,j,k,l) &
586  + 2.0*mxdt(j,l)*deltadt(i,j,k,l) + mxdt2(j,l)*delta(i,j,k,l) )
587  END DO
588  END DO
589  mx_itr(i,k) = 1.0 / (1.0 + suma * rho)
590  mx_itrdt(i,k)= - mx_itr(i,k)**2 * sumdt*rho
591  mx_itrdt2(i,k)= +2.0*mx_itr(i,k)**3 * (sumdt*rho)**2 - mx_itr(i,k)**2 *sumdt2*rho
592  END DO
593  END DO
594 
595  err_sum = 0.0
596  DO i = 1, ncomp
597  DO k = 1, nhb_typ(i)
598  err_sum = err_sum + abs(mx_itr(i,k) - mx(i,k)) &
599  + abs(mx_itrdt(i,k) - mxdt(i,k)) + abs(mx_itrdt2(i,k) - mxdt2(i,k))
600  mx(i,k) = mx_itr(i,k) * attenu + mx(i,k) * (1.0 - attenu)
601  mxdt(i,k)=mx_itrdt(i,k)*attenu +mxdt(i,k)* (1.0 - attenu)
602  mxdt2(i,k)=mx_itrdt2(i,k)*attenu +mxdt2(i,k)* (1.0 - attenu)
603  END DO
604  END DO
605  IF(err_sum <= tol) GO TO 10
606 
607  END DO
608  WRITE(6,*) 'CAL_PCSAFT: max_eval violated err_sum = ',err_sum,tol
609  stop
610 10 CONTINUE
611 
612  DO i = 1, ncomp
613  DO k = 1, nhb_typ(i)
614  ! fhb = fhb + x(i)* nhb_no(i,k)* ( 0.5 * ( 1.0 - mx(i,k) ) + LOG(mx(i,k)) )
615  fhbdt = fhbdt + x(i)*nhb_no(i,k) *mxdt(i,k)*(1.0/mx(i,k)-0.5)
616  fhbdt2= fhbdt2 + x(i)*nhb_no(i,k) *(mxdt2(i,k)*(1.0/mx(i,k)-0.5) &
617  -(mxdt(i,k)/mx(i,k))**2 )
618  END DO
619  END DO
620 
621  END IF
622 
623 
624  !-----------------------------------------------------------------------------
625  ! derivatives of f/kT of dipole-dipole term to temp. (fdddt)
626  !-----------------------------------------------------------------------------
627 
628  fdddt = 0.0
629  fdddt2 = 0.0
630  dipole = 0
631  DO i = 1, ncomp
632  my2dd(i) = 0.0
633  IF ( parame(i,6) /= 0.0 .AND. uij(i,i) /= 0.0 ) THEN
634  dipole = 1
635  my2dd(i) = (parame(i,6))**2 *1.e-49 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**3 *1.e-30)
636  END IF
637  my0(i) = my2dd(i) ! needed for dipole-quadrupole-term
638  END DO
639 
640  IF ( dipole == 1 ) THEN
641  DO i = 1, ncomp
642  DO j = 1, ncomp
643  idd2(i,j) =0.0
644  idd4(i,j) =0.0
645  idd2dt(i,j) =0.0
646  idd4dt(i,j) =0.0
647  idd2dt2(i,j)=0.0
648  idd4dt2(i,j)=0.0
649  DO m = 0, 4
650  idd2(i,j) = idd2(i,j) +ddp2(i,j,m)*z3**m
651  idd4(i,j) = idd4(i,j) +ddp4(i,j,m)*z3**m
652  idd2dt(i,j)= idd2dt(i,j) +ddp2(i,j,m)*z3dt*REAL(m)*z3**(m-1)
653  idd4dt(i,j)= idd4dt(i,j) +ddp4(i,j,m)*z3dt*REAL(m)*z3**(m-1)
654  idd2dt2(i,j)=idd2dt2(i,j)+ddp2(i,j,m)*z3dt2*REAL(m)*z3**(m-1) &
655  + ddp2(i,j,m)*z3dt**2 *REAL(m)*REAL(m-1)*z3**(m-2)
656  idd4dt2(i,j)=idd4dt2(i,j)+ddp4(i,j,m)*z3dt2*REAL(m)*z3**(m-1) &
657  + ddp4(i,j,m)*z3dt**2 *REAL(m)*REAL(m-1)*z3**(m-2)
658  END DO
659  DO k = 1, ncomp
660  idd3(i,j,k) =0.0
661  idd3dt(i,j,k) =0.0
662  idd3dt2(i,j,k)=0.0
663  DO m = 0, 4
664  idd3(i,j,k) = idd3(i,j,k) +ddp3(i,j,k,m)*z3**m
665  idd3dt(i,j,k) = idd3dt(i,j,k)+ddp3(i,j,k,m)*z3dt*REAL(m) *z3**(m-1)
666  idd3dt2(i,j,k)= idd3dt2(i,j,k)+ddp3(i,j,k,m)*z3dt2*REAL(m) &
667  *z3**(m-1) +ddp3(i,j,k,m)*z3dt**2 *REAL(m)*REAL(m-1) *z3**(m-2)
668  END DO
669  END DO
670  END DO
671  END DO
672 
673 
674  factor2= - pi *rho
675  factor3= - 4.0 / 3.0 * pi**2 * rho**2
676 
677  fdd2 = 0.0
678  fdd3 = 0.0
679  fdd2dt= 0.0
680  fdd3dt= 0.0
681  fdd2dt2= 0.0
682  fdd3dt2= 0.0
683  DO i = 1, ncomp
684  DO j = 1, ncomp
685  xijmt = x(i)*parame(i,3)*parame(i,2)**3 *x(j)*parame(j,3)*parame(j,2)**3 &
686  / ((parame(i,2)+parame(j,2))/2.0)**3 *my2dd(i)*my2dd(j)
687  eij = (parame(i,3)*parame(j,3))**0.5
688  fdd2 = fdd2 +factor2* xijmt/t/t*(idd2(i,j)+eij/t*idd4(i,j))
689  fdd2dt= fdd2dt+ factor2* xijmt/t/t*(idd2dt(i,j)-2.0*idd2(i,j)/t &
690  +eij/t*idd4dt(i,j)-3.0*eij/t/t*idd4(i,j))
691  fdd2dt2=fdd2dt2+factor2*xijmt/t/t*(idd2dt2(i,j)-4.0*idd2dt(i,j)/t &
692  +6.0*idd2(i,j)/t/t+eij/t*idd4dt2(i,j) &
693  -6.0*eij/t/t*idd4dt(i,j)+12.0*eij/t**3 *idd4(i,j))
694  DO k = 1, ncomp
695  xijkmt=x(i)*parame(i,3)*parame(i,2)**3 &
696  *x(j)*parame(j,3)*parame(j,2)**3 &
697  *x(k)*parame(k,3)*parame(k,2)**3 &
698  /((parame(i,2)+parame(j,2))/2.0) /((parame(i,2)+parame(k,2))/2.0) &
699  /((parame(j,2)+parame(k,2))/2.0) *my2dd(i)*my2dd(j)*my2dd(k)
700  fdd3 =fdd3 +factor3*xijkmt/t**3 *idd3(i,j,k)
701  fdd3dt =fdd3dt +factor3*xijkmt/t**3 * (idd3dt(i,j,k)-3.0*idd3(i,j,k)/t)
702  fdd3dt2=fdd3dt2+factor3*xijkmt/t**3 &
703  *( idd3dt2(i,j,k)-6.0*idd3dt(i,j,k)/t+12.0*idd3(i,j,k)/t/t )
704  END DO
705  END DO
706  END DO
707 
708  IF ( fdd2 < -1.e-100 .AND. fdd3 /= 0.0 ) THEN
709  fdddt = fdd2* (fdd2*fdd2dt - 2.0*fdd3*fdd2dt+fdd2*fdd3dt) / (fdd2-fdd3)**2
710  fdddt2 = ( 2.0*fdd2*fdd2dt*fdd2dt +fdd2*fdd2*fdd2dt2 &
711  -2.0*fdd2dt**2 *fdd3 -2.0*fdd2*fdd2dt2*fdd3 +fdd2*fdd2*fdd3dt2 ) &
712  /(fdd2-fdd3)**2 + fdddt * 2.0*(fdd3dt-fdd2dt)/(fdd2-fdd3)
713  END IF
714  END IF
715 
716 
717  !-----------------------------------------------------------------------------
718  ! derivatives f/kT of quadrupole-quadrup. term to T (fqqdt)
719  !-----------------------------------------------------------------------------
720 
721  fqqdt = 0.0
722  fqqdt2 = 0.0
723  qudpole = 0
724  DO i = 1, ncomp
725  qq2(i) = (parame(i,7))**2 *1.e-69 / (uij(i,i)*kbol*mseg(i)*sig_ij(i,i)**5 *1.e-50)
726  IF (qq2(i) /= 0.0) qudpole = 1
727  END DO
728 
729  IF (qudpole == 1) THEN
730 
731  DO i = 1, ncomp
732  DO j = 1, ncomp
733  iqq2(i,j) = 0.0
734  iqq4(i,j) = 0.0
735  iqq2dt(i,j) = 0.0
736  iqq4dt(i,j) = 0.0
737  iqq2dt2(i,j)= 0.0
738  iqq4dt2(i,j)= 0.0
739  DO m = 0, 4
740  iqq2(i,j) = iqq2(i,j) + qqp2(i,j,m)*z3**m
741  iqq4(i,j) = iqq4(i,j) + qqp4(i,j,m)*z3**m
742  iqq2dt(i,j) = iqq2dt(i,j)+ qqp2(i,j,m)*z3dt*REAL(m)*z3**(m-1)
743  iqq4dt(i,j) = iqq4dt(i,j)+ qqp4(i,j,m)*z3dt*REAL(m)*z3**(m-1)
744  iqq2dt2(i,j)= iqq2dt2(i,j)+qqp2(i,j,m)*z3dt2*REAL(m)*z3**(m-1) &
745  + qqp2(i,j,m)*z3dt**2 *REAL(m)*REAL(m-1)*z3**(m-2)
746  iqq4dt2(i,j)= iqq4dt2(i,j)+qqp4(i,j,m)*z3dt2*REAL(m)*z3**(m-1) &
747  + qqp4(i,j,m)*z3dt**2 *REAL(m)*REAL(m-1)*z3**(m-2)
748  END DO
749  DO k = 1, ncomp
750  iqq3(i,j,k) =0.0
751  iqq3dt(i,j,k) =0.0
752  iqq3dt2(i,j,k)=0.0
753  DO m = 0, 4
754  iqq3(i,j,k) = iqq3(i,j,k) + qqp3(i,j,k,m)*z3**m
755  iqq3dt(i,j,k) = iqq3dt(i,j,k)+ qqp3(i,j,k,m)*z3dt*REAL(m) * z3**(m-1)
756  iqq3dt2(i,j,k)= iqq3dt2(i,j,k)+qqp3(i,j,k,m)*z3dt2*REAL(m) * z3**(m-1) &
757  + qqp3(i,j,k,m)*z3dt**2 *REAL(m)*REAL(m-1)*z3**(m-2)
758  END DO
759  END DO
760  END DO
761  END DO
762 
763  factor2 = -9.0/16.0 * pi *rho
764  factor3 = 9.0/16.0 * pi**2 * rho**2
765 
766  fqq2 = 0.0
767  fqq3 = 0.0
768  fqq2dt = 0.0
769  fqq3dt = 0.0
770  fqq2dt2= 0.0
771  fqq3dt2= 0.0
772  DO i = 1, ncomp
773  DO j = 1, ncomp
774  xijmt = x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 &
775  * x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /sig_ij(i,j)**7.0
776  eij = (parame(i,3)*parame(j,3))**0.5
777  fqq2 = fqq2 +factor2* xijmt/t/t*(iqq2(i,j)+eij/t*iqq4(i,j))
778  fqq2dt= fqq2dt +factor2* xijmt/t/t*(iqq2dt(i,j)-2.0*iqq2(i,j)/t &
779  + eij/t*iqq4dt(i,j)-3.0*eij/t/t*iqq4(i,j))
780  fqq2dt2=fqq2dt2+factor2*xijmt/t/t*(iqq2dt2(i,j)-4.0*iqq2dt(i,j)/t &
781  + 6.0*iqq2(i,j)/t/t+eij/t*iqq4dt2(i,j) &
782  - 6.0*eij/t/t*iqq4dt(i,j)+12.0*eij/t**3 *iqq4(i,j))
783  DO k = 1, ncomp
784  xijkmt = x(i)*uij(i,i)*qq2(i)*sig_ij(i,i)**5 /sig_ij(i,j)**3 &
785  * x(j)*uij(j,j)*qq2(j)*sig_ij(j,j)**5 /sig_ij(i,k)**3 &
786  * x(k)*uij(k,k)*qq2(k)*sig_ij(k,k)**5 /sig_ij(j,k)**3
787  fqq3 = fqq3 +factor3*xijkmt/t**3 *iqq3(i,j,k)
788  fqq3dt = fqq3dt +factor3*xijkmt/t**3 *(iqq3dt(i,j,k)-3.0*iqq3(i,j,k)/t)
789  fqq3dt2= fqq3dt2+factor3*xijkmt/t**3 &
790  * ( iqq3dt2(i,j,k)-6.0*iqq3dt(i,j,k)/t+12.0*iqq3(i,j,k)/t/t )
791  END DO
792  END DO
793  END DO
794 
795  IF ( fqq2 /= 0.0 .AND. fqq3 /= 0.0 ) THEN
796  fqqdt = fqq2* (fqq2*fqq2dt - 2.0*fqq3*fqq2dt+fqq2*fqq3dt) / (fqq2-fqq3)**2
797  fqqdt2 = ( 2.0*fqq2*fqq2dt*fqq2dt +fqq2*fqq2*fqq2dt2 &
798  - 2.0*fqq2dt**2 *fqq3 -2.0*fqq2*fqq2dt2*fqq3 +fqq2*fqq2*fqq3dt2 ) &
799  / (fqq2-fqq3)**2 + fqqdt * 2.0*(fqq3dt-fqq2dt)/(fqq2-fqq3)
800  END IF
801 
802  END IF
803 
804 
805  !-----------------------------------------------------------------------------
806  ! derivatives f/kT of dipole-quadruppole term to T (fdqdt)
807  !-----------------------------------------------------------------------------
808 
809  fdqdt = 0.0
810  fdqdt2= 0.0
811  dip_quad = 0
812  DO i = 1, ncomp
813  DO j = 1, ncomp
814  IF (parame(i,6) /= 0.0 .AND. parame(j,7) /= 0.0) dip_quad = 1
815  END DO
816  myfac(i) = parame(i,3) * parame(i,2)**4 * my0(i)
817  q_fac(i) = parame(i,3) * parame(i,2)**4 * qq2(i)
818  END DO
819 
820  IF (dip_quad == 1) THEN
821 
822  DO i = 1, ncomp
823  DO j = 1, ncomp
824  idq2(i,j) = 0.0
825  idq4(i,j) = 0.0
826  idq2dt(i,j) = 0.0
827  idq4dt(i,j) = 0.0
828  idq2dt2(i,j)= 0.0
829  idq4dt2(i,j)= 0.0
830  IF ( myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0 ) THEN
831  DO m = 0, 4
832  idq2(i,j) = idq2(i,j) + dqp2(i,j,m)*z3**m
833  idq4(i,j) = idq4(i,j) + dqp4(i,j,m)*z3**m
834  idq2dt(i,j) = idq2dt(i,j)+ dqp2(i,j,m)*z3dt*REAL(m)*z3**(m-1)
835  idq4dt(i,j) = idq4dt(i,j)+ dqp4(i,j,m)*z3dt*REAL(m)*z3**(m-1)
836  idq2dt2(i,j)= idq2dt2(i,j)+dqp2(i,j,m)*z3dt2*REAL(m)*z3**(m-1) &
837  + dqp2(i,j,m)*z3dt**2 *REAL(m)*REAL(m-1)*z3**(m-2)
838  idq4dt2(i,j)= idq4dt2(i,j)+dqp4(i,j,m)*z3dt2*REAL(m)*z3**(m-1) &
839  + dqp4(i,j,m)*z3dt**2 *REAL(m)*REAL(m-1)*z3**(m-2)
840  END DO
841 
842  DO k = 1, ncomp
843  idq3(i,j,k) = 0.0
844  idq3dt(i,j,k) = 0.0
845  idq3dt2(i,j,k)= 0.0
846  IF ( myfac(k) /= 0.0 .OR. q_fac(k) /= 0.0 ) THEN
847  DO m = 0, 4
848  idq3(i,j,k) = idq3(i,j,k) + dqp3(i,j,k,m)*z3**m
849  idq3dt(i,j,k)= idq3dt(i,j,k)+ dqp3(i,j,k,m)*z3dt*REAL(m) *z3**(m-1)
850  idq3dt2(i,j,k)= idq3dt2(i,j,k)+dqp3(i,j,k,m)*z3dt2*REAL(m) *z3**(m-1) &
851  + dqp3(i,j,k,m)*z3dt**2 *REAL(m)*REAL(m-1)*z3**(m-2)
852  END DO
853  END IF
854  END DO
855  END IF
856  END DO
857  END DO
858 
859  factor2= -9.0/4.0 * pi * rho
860  factor3= pi**2 * rho**2
861 
862  fdq2 = 0.0
863  fdq3 = 0.0
864  fdq2dt= 0.0
865  fdq3dt= 0.0
866  fdq2dt2=0.0
867  fdq3dt2=0.0
868  DO i = 1, ncomp
869  DO j = 1, ncomp
870  IF ( myfac(i) /= 0.0 .AND. q_fac(j) /= 0.0 ) THEN
871  xijmt = x(i)*myfac(i) * x(j)*q_fac(j) /sig_ij(i,j)**5
872  eij = (parame(i,3)*parame(j,3))**0.5
873  fdq2 = fdq2 + factor2* xijmt/t/t*(idq2(i,j)+eij/t*idq4(i,j))
874  fdq2dt= fdq2dt+ factor2* xijmt/t/t*(idq2dt(i,j)-2.0*idq2(i,j)/t &
875  + eij/t*idq4dt(i,j)-3.0*eij/t/t*idq4(i,j))
876  fdq2dt2 = fdq2dt2+factor2*xijmt/t/t*(idq2dt2(i,j)-4.0*idq2dt(i,j)/t &
877  + 6.0*idq2(i,j)/t/t+eij/t*idq4dt2(i,j) &
878  - 6.0*eij/t/t*idq4dt(i,j)+12.0*eij/t**3 *idq4(i,j))
879  DO k = 1, ncomp
880  IF ( myfac(k) /= 0.0 .OR. q_fac(k) /= 0.0 ) THEN
881  xijkmt= x(i)*x(j)*x(k)/(sig_ij(i,j)*sig_ij(i,k)*sig_ij(j,k))**2 &
882  * ( myfac(i)*q_fac(j)*myfac(k) &
883  + myfac(i)*q_fac(j)*q_fac(k)*1.193735 )
884 
885  fdq3 =fdq3 + factor3*xijkmt/t**3 *idq3(i,j,k)
886  fdq3dt=fdq3dt+ factor3*xijkmt/t**3 * (idq3dt(i,j,k)-3.0*idq3(i,j,k)/t)
887  fdq3dt2=fdq3dt2+factor3*xijkmt/t**3 &
888  *( idq3dt2(i,j,k)-6.0*idq3dt(i,j,k)/t+12.0*idq3(i,j,k)/t/t )
889  END IF
890  END DO
891  END IF
892  END DO
893  END DO
894 
895  IF (fdq2 /= 0.0 .AND. fdq3 /= 0.0) THEN
896  fdqdt = fdq2* (fdq2*fdq2dt - 2.0*fdq3*fdq2dt+fdq2*fdq3dt) / (fdq2-fdq3)**2
897  fdqdt2 = ( 2.0*fdq2*fdq2dt*fdq2dt +fdq2*fdq2*fdq2dt2 &
898  - 2.0*fdq2dt**2 *fdq3 -2.0*fdq2*fdq2dt2*fdq3 +fdq2*fdq2*fdq3dt2 ) &
899  / (fdq2-fdq3)**2 + fdqdt * 2.0*(fdq3dt-fdq2dt)/(fdq2-fdq3)
900  END IF
901 
902  END IF
903  !-----------------------------------------------------------------------------
904 
905 
906 
907 
908  !-----------------------------------------------------------------------------
909  ! total derivative of fres/kT to temperature
910  !-----------------------------------------------------------------------------
911 
912  df_dt = fhsdt + fchdt + fdspdt + fhbdt + fdddt + fqqdt + fdqdt
913 
914 
915 
916  !-----------------------------------------------------------------------------
917  ! second derivative of fres/kT to T
918  !-----------------------------------------------------------------------------
919 
920  df_dt2 = fhsdt2 +fchdt2 +fdspdt2 +fhbdt2 +fdddt2 +fqqdt2 +fdqdt2
921 
922 
923 
924  !-----------------------------------------------------------------------------
925  ! derivatives of fres/kt to density and to T
926  !-----------------------------------------------------------------------------
927 
928  !-----------------------------------------------------------------------------
929  ! the analytic derivative of fres/kT to (density and T) (df_drdt)
930  ! is still missing. A numerical differentiation is implemented.
931  !-----------------------------------------------------------------------------
932 
933  fact = 1.0
934  dist = t * 100.e-5 * fact
935  t_tmp = t
936  rho_0 = rho
937 
938 
939  t = t_tmp - 2.0*dist
940  CALL perturbation_parameter
941  eta = z3t*rho_0
942  CALL p_eos
943  fdr1 = pges / (eta*rho_0*(kbol*t)/1.e-30)
944  t = t_tmp - dist
945  CALL perturbation_parameter
946  eta = z3t*rho_0
947  CALL p_eos
948  fdr2 = pges / (eta*rho_0*(kbol*t)/1.e-30)
949 
950  t = t_tmp + dist
951  CALL perturbation_parameter
952  eta = z3t*rho_0
953  CALL p_eos
954  fdr3 = pges / (eta*rho_0*(kbol*t)/1.e-30)
955 
956  t = t_tmp + 2.0*dist
957  CALL perturbation_parameter
958  eta = z3t*rho_0
959  CALL p_eos
960  fdr4 = pges / (eta*rho_0*(kbol*t)/1.e-30)
961 
962  t = t_tmp
963  CALL perturbation_parameter
964  eta = z3t*rho_0
965  CALL p_eos
966 
967 
968  df_drdt = (-fdr4+8.0*fdr3-8.0*fdr2+fdr1)/(12.0*dist)
969 
970 
971 
972 
973  !-----------------------------------------------------------------------------
974  ! thermodynamic properties
975  !-----------------------------------------------------------------------------
976 
977  s_res = ( - df_dt *t - fres )*rgas + rgas * log(zges)
978  h_res = ( - t*df_dt + zges-1.0 ) * rgas *t
979  cv_res = - ( t*df_dt2 + 2.0*df_dt ) * rgas*t
980  cp_res = cv_res - rgas + rgas*(zges +eta*t*df_drdt)**2 &
981  / (1.0 + 2.0*eta*dfdr +eta**2 *ddfdrdr)
982 
983  !-----------------------------------------------------------------------------
984  ! speed of sound
985  !-----------------------------------------------------------------------------
986 
987 !!! cp_ig_mix = 1.57797619E-04 * t * t + 1.76863690E-01 * t + 2.98855179E+01 ! in units J/(mol*K)
988  ! mm_mean = SUM( x(1:ncomp) * mm(1:ncomp) )
989 
990  ! w_square = 1.0 + 2.0*eta*dfdr + eta**2 *ddfdrdr &
991  ! + 1.0/( cp_ig_mix/RGAS - 1.0 + cv_res/RGAS ) &
992  ! * ( zges + eta*t *df_drdt )**2
993  ! speed_sound = SQRT( w_square * RGAS * t / (mm_mean/1000.0) )
994 
995 
996  !-----------------------------------------------------------------------------
997  ! write (*,*) 'df_... ', df_dt,df_dt2
998  ! write (*,*) 'kreuz ', zges,eta*t*df_drdt,eta*dfdr, eta**2 *ddfdrdr
999  ! write (*,*) 'h,cv,cp', h_res,cv_res,cp_res
1000 
1001 
1002 END SUBROUTINE h_eos
1003 
1004 
1005 
1006 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1007 ! SUBROUTINE H_EOS_num
1008 !
1009 ! This subroutine calculates enthalpies and heat capacities (cp) by
1010 ! taking numerical derivatieves.
1011 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1012 
1013 SUBROUTINE h_eos_num
1014 
1015  USE parameters, ONLY: rgas
1016  USE eos_variables
1017  IMPLICIT NONE
1018 
1019  !-----------------------------------------------------------------------------
1020  REAL :: dist, fact, rho_0
1021  REAL :: fres1, fres2, fres3, fres4, fres5
1022  REAL :: f_1, f_2, f_3, f_4
1023  REAL :: cv_res, t_tmp, zges
1024  REAL :: df_dt, df_dtdt, df_drdt, dfdr, ddfdrdr
1025  !-----------------------------------------------------------------------------
1026 
1027 
1028  CALL perturbation_parameter
1029  rho_0 = eta/z3t
1030 
1031 
1032  fact = 1.0
1033  dist = t * 100.e-5 * fact
1034 
1035  t_tmp = t
1036 
1037  t = t_tmp - 2.0*dist
1038  CALL perturbation_parameter
1039  eta = z3t * rho_0
1040  CALL f_eos
1041  fres1 = fres
1042  t = t_tmp - dist
1043  CALL perturbation_parameter
1044  eta = z3t * rho_0
1045  CALL f_eos
1046  fres2 = fres
1047  t = t_tmp + dist
1048  CALL perturbation_parameter
1049  eta = z3t * rho_0
1050  CALL f_eos
1051  fres3 = fres
1052  t = t_tmp + 2.0*dist
1053  CALL perturbation_parameter
1054  eta = z3t * rho_0
1055  CALL f_eos
1056  fres4 = fres
1057  t = t_tmp
1058  CALL perturbation_parameter
1059  eta = z3t * rho_0
1060  CALL f_eos
1061  fres5 = fres
1062  ! *(KBOL*T)/1.E-30
1063 
1064  zges = (p * 1.e-30)/(kbol*t*rho_0)
1065 
1066 
1067  df_dt = (-fres4+8.0*fres3-8.0*fres2+fres1)/(12.0*dist)
1068  df_dtdt = (-fres4+16.0*fres3-3.d1*fres5+16.0*fres2-fres1) &
1069  /(12.0*(dist**2 ))
1070 
1071 
1072  s_res = (- df_dt -fres/t)*rgas*t + rgas * log(zges)
1073  h_res = ( - t*df_dt + zges-1.0 ) * rgas*t
1074  cv_res = -( t*df_dtdt + 2.0*df_dt ) * rgas*t
1075 
1076 
1077 
1078  t = t_tmp - 2.0*dist
1079  CALL perturbation_parameter
1080  eta = z3t * rho_0
1081  CALL p_eos
1082  f_1 = pges/(eta*rho_0*(kbol*t)/1.e-30)
1083 
1084  t = t_tmp - dist
1085  CALL perturbation_parameter
1086  eta = z3t * rho_0
1087  CALL p_eos
1088  f_2 = pges/(eta*rho_0*(kbol*t)/1.e-30)
1089 
1090  t = t_tmp + dist
1091  CALL perturbation_parameter
1092  eta = z3t * rho_0
1093  CALL p_eos
1094  f_3 = pges/(eta*rho_0*(kbol*t)/1.e-30)
1095 
1096  t = t_tmp + 2.0*dist
1097  CALL perturbation_parameter
1098  eta = z3t * rho_0
1099  CALL p_eos
1100  f_4 = pges/(eta*rho_0*(kbol*t)/1.e-30)
1101 
1102  t = t_tmp
1103  CALL perturbation_parameter
1104  eta = z3t * rho_0
1105  CALL p_eos
1106 
1107  dfdr = pges / ( eta * rho_0 * (kbol*t) / 1.e-30 )
1108  ddfdrdr = pgesdz / ( eta * rho_0*(kbol*t)/1.e-30 ) - dfdr * 2.0 / eta - 1.0 / eta**2
1109 
1110  df_drdt = ( -f_4 + 8.0*f_3 - 8.0*f_2 + f_1) / ( 12.0 * dist )
1111 
1112  cp_res = cv_res - rgas + rgas *(zges + eta*t*df_drdt)**2 / (1.0 +2.0*eta*dfdr +eta**2 *ddfdrdr)
1113 
1114  ! write (*,*) 'n',df_dt,df_dtdt
1115  ! write (*,*) 'kreuz ', zges,eta*t*df_drdt,eta*dfdr, eta**2 *ddfdrdr
1116  ! write (*,*) 'h, cv', h_res, cv_res
1117  ! write (*,*) h_res - t*s_res
1118  ! write (*,*) cv_res,cp_res,eta
1119  ! pause
1120 
1121 END SUBROUTINE h_eos_num
1122 
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
sig
double sig
correlated to variance of initial distribution
Definition: experimentalInputs.h:138
eos_constants
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:200
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29