MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
eos.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! SUBROUTINE FUGACITY
3 !
4 ! This subroutine serves as an interface to the eos-subroutines.
5 ! (1) case 1, when ensemble_flag = 'tp'
6 ! The subroutine gives the residual chemical potential:
7 ! mu_i^res(T,p,x)/kT = ln( phi_i )
8 ! and in addition, the densities that satisfy the specified p
9 ! (2) case 2, when ensemble_flag = 'tv'
10 ! The subroutine gives the residual chemical potential:
11 ! --> mu_i^res(T,rho,x)/kT
12 ! and in addition the resulting pressure for the given density.
13 ! The term "residual" means difference of the property and the same
14 ! property for an ideal gas mixture.
15 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
16 
17 SUBROUTINE fugacity (ln_phi)
18 
19  USE basic_variables, only: nc, np, eos, nphas, ncomp, xi, dense, densta, p_cal, &
20  z_cal, my_cal, rhoi_cal, f_res, gibbs, ensemble_flag
21  USE eos_variables, ONLY: phas, x, t, eta, eta_start, lnphi, fres, rho, pges, kbol
22  IMPLICIT NONE
23 
24  !-----------------------------------------------------------------------------
25  REAL, INTENT(OUT) :: ln_phi(np,nc)
26 
27  !-----------------------------------------------------------------------------
28  INTEGER :: i, ph
29  LOGICAL :: trivial_result
30  !-----------------------------------------------------------------------------
31 
32  IF (eos < 2) THEN
33 
34  DO ph = 1,nphas
35 
36  phas = ph
37  eta_start = densta(ph)
38  x(1:ncomp) = xi(ph,1:ncomp)
39 
40  CALL check_eos_variables ( trivial_result )
41  IF ( trivial_result ) return
42 
43  !-----------------------------------------------------------------------
44  ! calculate chemical potential and other quantities
45  !-----------------------------------------------------------------------
46 
47  CALL phi_eos_interface
48 
49  !-----------------------------------------------------------------------
50  ! densities, pressure and Z
51  !-----------------------------------------------------------------------
52  dense(ph) = eta
53  rhoi_cal(ph,1:ncomp) = rho * x(1:ncomp)
54 
55  p_cal(ph) = pges
56  z_cal(ph) = pges / ( kbol*1.e30 * t * rho )
57 
58  !-----------------------------------------------------------------------
59  ! chemical potential / kT and ln( fugacity coefficient )
60  !-----------------------------------------------------------------------
61  ln_phi(ph,1:ncomp) = lnphi(1:ncomp)
62 
63  do i = 1,ncomp
64  if ( ( x(i) * rho ) > 1.e-200 ) then
65  my_cal( ph, i ) = lnphi( i ) + log( x(i) * rho )
66  else
67  my_cal( ph, i ) = - 1.e200
68  end if
69  end do
70  if ( ensemble_flag == 'tp' ) my_cal(ph,1:ncomp) = my_cal(ph,1:ncomp) + log(z_cal(ph))
71 
72  !-----------------------------------------------------------------------
73  ! Gibbs energy / kT and Helmholtz energy density / kT ( not Helmholtz energy!)
74  !-----------------------------------------------------------------------
75  gibbs(ph) = sum( x(1:ncomp) * lnphi(1:ncomp) )
76  f_res(ph) = fres * rho
77  do i = 1,ncomp
78  if ( x(i) > 1.e-200 ) gibbs(ph) = gibbs(ph) + x(i) * log( x(i)*rho )
79  if ( x(i) > 1.e-200 ) f_res(ph) = f_res(ph) + x(i) * rho * ( log( x(i)*rho) - 1.0 )
80  end do
81  if ( ensemble_flag == 'tp' ) gibbs(ph) = gibbs(ph) + log( z_cal(ph) )
82 
83  ! write (*,'(i3,4G20.11)') ph,eta,lnphi(1),lnphi(2),x(1)
84  ! DO i = 1,ncomp
85  ! DO j=1,NINT(parame(i,12))
86  ! mxx(ph,i,j) = mx(i,j)
87  ! END DO
88  ! END DO
89 
90  END DO
91 
92  ELSE
93 
94  ! IF (eos == 2) CALL srk_eos (ln_phi)
95  ! IF (eos == 3) CALL pr_eos (ln_phi)
96  ! dense(1) = 0.01
97  ! dense(2) = 0.3
98  ! IF (eos == 4.OR.eos == 5.OR.eos == 6.OR.eos == 8) CALL lj_fugacity(ln_phi)
99  ! IF (eos == 7) CALL sw_fugacity(ln_phi)
100  ! IF (eos == 9) CALL lj_bh_fug(ln_phi)
101 
102  END IF
103 
104 END SUBROUTINE fugacity
105 
106 
107 
108 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
109 ! SUBROUTINE CHECK_EOS_VARIABLES
110 !
111 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
112 
113 SUBROUTINE check_eos_variables ( trivial_result )
114 
115  USE basic_variables
116  USE eos_variables, ONLY: x, eta, eta_start, lnphi, fres, rho, pges, kbol, z3t
117  USE utilities
118  IMPLICIT NONE
119 
120  !-----------------------------------------------------------------------------
121  LOGICAL, INTENT(OUT) :: trivial_result
122 
123  !-----------------------------------------------------------------------------
124  INTEGER :: i
125  REAL :: sum_x
126  !-----------------------------------------------------------------------------
127 
128  trivial_result = .false.
129 
130  !-----------------------------------------------------------------------------
131  ! verify specification: either (T,p,x) or (T,rho,x)-variables
132  !-----------------------------------------------------------------------------
133 
134  IF ( ensemble_flag /= 'tp' .AND. ensemble_flag /= 'tv' ) THEN
135  WRITE(*,*) ' FUGACITY: variable ensemble_flag is undefined ',ensemble_flag
136  stop
137  END IF
138 
139  !-----------------------------------------------------------------------------
140  ! check for NaN
141  !-----------------------------------------------------------------------------
142 
143  IF ( eta_start /= eta_start ) THEN
144  WRITE(*,*) ' FUGACITY: density input is "not a number" !'
145  trivial_result = .true.
146  END IF
147 
148  DO i = 1, ncomp
149  IF ( x(i) /= x(i) ) THEN
150  WRITE(*,*) ' FUGACITY: composition input x is "not a number" ! Species:',i
151  trivial_result = .true.
152  END IF
153  END DO
154 
155  !-----------------------------------------------------------------------------
156  ! verify proper specification of mole fractions x
157  !-----------------------------------------------------------------------------
158 
159  sum_x = sum( x(1:ncomp) )
160  IF ( sum_x /= 1.0 ) THEN
161  IF ( (sum_x - 1.0 ) < 1.e-4 .AND. (sum_x - 1.0 ) > -1.e-4 ) THEN
162  x( 1:ncomp ) = x( 1:ncomp ) / sum_x
163  ! WRITE(*,*) ' FUGACITY: rescale composition',sum_x, x(1:ncomp)
164  END IF
165  END IF
166 
167  !-----------------------------------------------------------------------------
168  ! if sum is not unity within 10.E-8, then at this point, the sum has to be
169  ! outside of 10.E-4 bandwidth
170  !-----------------------------------------------------------------------------
171  sum_x = sum( x(1:ncomp) )
172  IF ( (sum_x - 1.0 ) > 1.e-8 .OR. (sum_x - 1.0 ) < -1.e-8 ) THEN
173  ! WRITE(*,*) ' FUGACITY: composition not properly defined',sum( x(1:ncomp) )
174  if ( sum( x(1:ncomp) ) < 0.2 ) x(1:ncomp) = x(1:ncomp) / sum( x(1:ncomp) )
175  ! call pause
176  END IF
177 
178  IF ( minval( x(1:ncomp) ) < 0.0 ) THEN
179  WRITE(*,*) ' FUGACITY: mole fraction x is negative, of species', minloc(x(1:ncomp))
180  stop
181  END IF
182  DO i = 1, ncomp
183  IF ( x(i) < 1.e-50 ) x(i) = 0.0
184  END DO
185 
186  !-----------------------------------------------------------------------------
187  ! verify proper specification of either pressure (p), or of density (eta)
188  !-----------------------------------------------------------------------------
189 
190  IF ( ensemble_flag == 'tp' .AND. p < 1.e-100 ) THEN
191  WRITE(*,*) ' FUGACITY: PRESSURE TOO LOW', p
192  p = 1.e-6
193  END IF
194 
195  IF ( ensemble_flag == 'tv' .AND. eta_start < 1.e-100 ) THEN
196  WRITE(*,*) ' FUGACITY: DENSITY TOO LOW', eta_start
197  eta = 1.e-100
198  trivial_result = .true.
199  END IF
200 
201  !-----------------------------------------------------------------------------
202  ! for too low density, don't execute the subroutine PHI_EOS
203  !-----------------------------------------------------------------------------
204 
205  IF ( trivial_result ) THEN
206  CALL perturbation_parameter
207  rho = eta / z3t
208  pges = kbol*1.e30 * t * rho
209  lnphi(1:ncomp) = 1.0
210  fres = 0.0
211  END IF
212 
213 END SUBROUTINE check_eos_variables
214 
215 
216 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
217 ! subroutine p_eos_interface
218 !
219 ! This subroutine serves as interface to the different versions of the PC-SAFT eos.
220 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
221 
222 subroutine p_eos_interface
223 
224  use basic_variables, only: num
225  use eos_variables
226  use eos_numerical
227 
228  IF (num == 0) THEN
229  CALL p_eos
230  ELSE IF (num == 1) THEN
231  CALL p_numerical
232  ELSE IF (num == 2) THEN
233  CALL f_eos_rn
234  ELSE
235  write (*,*) 'p_eos_interface: define calculation option (num)'
236  END IF
237 
238 end subroutine p_eos_interface
239 
240 
241 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
242 ! subroutine f_eos_interface
243 !
244 ! This subroutine serves as interface to the different versions of the PC-SAFT eos.
245 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
246 
247 subroutine f_eos_interface
248 
249  use basic_variables, only: num
250  use eos_variables
251  use eos_numerical
252 
253  IF (num == 0) THEN
254  CALL f_eos
255  ELSE IF (num == 1) THEN
256  CALL f_numerical
257  ELSE IF (num == 2) THEN
258  CALL f_eos_rn
259  ELSE
260  write (*,*) 'f_eos_interface: define calculation option (num)'
261  END IF
262 
263 
264 end subroutine f_eos_interface
265 
266 
267 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
268 ! subroutine phi_eos_interface
269 !
270 ! This subroutine serves as interface to the different versions of the PC-SAFT eos.
271 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
272 
273 subroutine phi_eos_interface
274 
275  use basic_variables, only: num
276  use eos_variables
277  use eos_numerical
278 
279  IF (num == 0) THEN
280  CALL phi_eos
281  ELSE IF (num == 1) THEN
282  CALL phi_numerical
283  ELSE IF (num == 2) THEN
284  CALL phi_critical_renorm
285  ELSE
286  write (*,*) 'phi_eos_interface: define calculation option (num)'
287  END IF
288 
289 end subroutine phi_eos_interface
290 
291 
292 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
293 ! subroutine h_eos_interface
294 !
295 ! This subroutine serves as interface to the different versions of the PC-SAFT eos.
296 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
297 
298 subroutine h_eos_interface
299 
300  use basic_variables, only: num
301  use eos_variables
302  use eos_numerical
303 
304  IF (num == 0) THEN
305  CALL h_eos
306  ELSE IF (num == 1) THEN
307  CALL h_numerical
308  ELSE IF (num == 2) THEN
309  write (*,*) 'enthalpy_etc: incorporate H_EOS_RN'
310  stop
311  ! CALL H_EOS_rn
312  ELSE
313  write (*,*) 'phi_eos_interface: define calculation option (num)'
314  END IF
315 
316 end subroutine h_eos_interface
317 
318 
319 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
320 ! SUBROUTINE PHI_EOS
321 !
322 ! This subroutine gives the residual chemical potential:
323 ! --> mu_i^res(T,p,x)/kT = ln( phi_i ) when ensemble_flag = 'tp'
324 ! The required input for this case (T, p, x(nc)) and as a starting value
325 ! eta_start
326 !
327 ! or it gives
328 !
329 ! --> mu_i^res(T,rho,x)/kT when ensemble_flag = 'tv'
330 ! The required input for this case (T, eta_start, x(nc)). Note that
331 ! eta_start is the specified density (packing fraction) in this case.
332 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
333 
334 SUBROUTINE phi_eos
335 
336  USE parameters
337  USE eos_variables
338  USE eos_constants
339  USE eos_polar, only: phi_polar
340  IMPLICIT NONE
341 
342  !-----------------------------------------------------------------------------
343  INTEGER :: i, j, k, ii, jj, kk, m
344  REAL :: z0, z1, z2, z3, z0_rk, z1_rk, z2_rk, z3_rk
345  REAL :: ome, ome2, ome3, m_mean
346  REAL, DIMENSION(nc) :: mhs, mdsp, mhc, myres
347  REAL :: z3_m
348  REAL :: m_rk
349  REAL :: gij_rk(nc,nc)
350  REAL :: zres, zges
351  REAL :: dpdz, dpdz2
352 
353  REAL :: i1, i2, i1_rk, i2_rk
354  REAL :: ord1_rk, ord2_rk
355  REAL :: c1_con, c2_con, c1_rk
356  REAL, DIMENSION(nc,0:6) :: ap_rk, bp_rk
357 
358  LOGICAL :: assoc
359  REAL :: ass_s2, m_hbon(nc)
360 
361  REAL :: fdd_rk, fqq_rk, fdq_rk
362  REAL, DIMENSION(nc) :: my_dd, my_qq, my_dq
363  !-----------------------------------------------------------------------------
364 
365  !-----------------------------------------------------------------------------
366  ! obtain parameters and density independent expressions
367  !-----------------------------------------------------------------------------
368  CALL perturbation_parameter
369 
370 
371  !-----------------------------------------------------------------------------
372  ! density iteration: (pTx)-ensemble OR p calc.: (pvx)-ensemble
373  !-----------------------------------------------------------------------------
374  IF ( ensemble_flag == 'tp' ) THEN
375  CALL density_iteration
376  ELSEIF ( ensemble_flag == 'tv' ) THEN
377  eta = eta_start
378  CALL p_eos
379  END IF
380 
381  ! --- Eq.(A.8) ---------------------------------------------------------------
382  rho = eta / z3t
383  IF ( rho /= rho ) write (*,*) 'PHI_EOS: error in density',eta, z3t
384  IF ( rho /= rho ) stop
385  z0 = z0t * rho
386  z1 = z1t * rho
387  z2 = z2t * rho
388  z3 = z3t * rho
389 
390  m_mean = z0t / (pi/6.0)
391  ome = 1.0 - eta
392  ome2 = ome * ome
393  ome3 = ome * ome2
394 
395  !-----------------------------------------------------------------------------
396  ! compressibility factor z = p/(kT*rho)
397  !-----------------------------------------------------------------------------
398  zges = (p * 1.e-30) / (kbol*t*rho)
399  IF ( ensemble_flag == 'tv' ) zges = (pges * 1.e-30) / (kbol*t*rho)
400  zres = zges - 1.0
401 
402 
403 
404  !=============================================================================
405  ! calculate the derivatives of f to mole fraction x ( d(f)/d(rho_k) )
406  !=============================================================================
407 
408  DO k = 1, ncomp
409 
410  z0_rk = pi/6.0 * mseg(k)
411  z1_rk = z0_rk * dhs(k)
412  z2_rk = z1_rk * dhs(k)
413  z3_rk = z2_rk * dhs(k)
414 
415  ! --- derivative d(m_mean)/d(rho_k) ---------------------------------------
416  m_rk = ( mseg(k) - m_mean ) / rho
417  ! lij(1,2)= -0.050
418  ! lij(2,1)=lij(1,2)
419  ! r_m2dx(k)=0.0
420  ! m_mean2=0.0
421  ! DO i =1,ncomp
422  ! r_m2dx(k)=r_m2dx(k)+2.0*x(i)*(mseg(i)+mseg(k))/2.0*(1.0-lij(i,k))
423  ! DO j =1,ncomp
424  ! m_mean2=m_mean2+x(i)*x(j)*(mseg(i)+mseg(j))/2.0*(1.0-lij(i,j))
425  ! ENDDO
426  ! ENDDO
427 
428  !--------------------------------------------------------------------------
429  ! d(f)/d(rho_k) : hard sphere contribution
430  !--------------------------------------------------------------------------
431  if ( z3**3 > 0.0 ) then
432  mhs(k) = 6.0/pi* ( 3.0*(z1_rk*z2+z1*z2_rk)/ome + 3.0*z1*z2*z3_rk/ome2 &
433  + 3.0*z2*z2*z2_rk/z3/ome2 + z2**3 *z3_rk*(3.0*z3-1.0)/z3/z3/ome3 &
434  + ((3.0*z2*z2*z2_rk*z3-2.0*z2**3 *z3_rk)/z3**3 -z0_rk) *log(ome) &
435  + (z0-z2**3 /z3/z3)*z3_rk/ome )
436  end if
437 
438  !--------------------------------------------------------------------------
439  ! d(f)/d(rho_k) : chain term
440  !--------------------------------------------------------------------------
441  DO i = 1, ncomp
442  DO j = 1, ncomp
443  gij(i,j) = 1.0/ome + 3.0*dij_ab(i,j)*z2/ome2 + 2.0*(dij_ab(i,j)*z2)**2 /ome3
444  gij_rk(i,j) = z3_rk/ome2 &
445  + 3.0*dij_ab(i,j)*(z2_rk+2.0*z2*z3_rk/ome)/ome2 &
446  + dij_ab(i,j)**2 *z2/ome3 *(4.0*z2_rk+6.0*z2*z3_rk/ome)
447  END DO
448  END DO
449 
450  mhc(k) = 0.0
451  DO i = 1, ncomp
452  mhc(k) = mhc(k) + x(i)*rho * (1.0-mseg(i)) / gij(i,i) * gij_rk(i,i)
453  END DO
454  mhc(k) = mhc(k) + ( 1.0-mseg(k)) * log( gij(k,k) )
455 
456 
457  !--------------------------------------------------------------------------
458  ! PC-SAFT: d(f)/d(rho_k) : dispersion contribution
459  !--------------------------------------------------------------------------
460 
461  ! --- derivatives of apar, bpar to rho_k -------------------------------
462  DO m = 0, 6
463  ap_rk(k,m) = m_rk/m_mean**2 * ( ap(m,2) + (3.0 -4.0/m_mean) *ap(m,3) )
464  bp_rk(k,m) = m_rk/m_mean**2 * ( bp(m,2) + (3.0 -4.0/m_mean) *bp(m,3) )
465  END DO
466 
467  i1 = 0.0
468  i2 = 0.0
469  i1_rk = 0.0
470  i2_rk = 0.0
471  DO m = 0, 6
472  z3_m = eta**m
473  i1 = i1 + apar(m) * z3_m
474  i2 = i2 + bpar(m) * z3_m
475  i1_rk = i1_rk + apar(m) * REAL(m) * eta**(m-1) * z3_rk + ap_rk(k,m) * z3_m
476  i2_rk = i2_rk + bpar(m) * REAL(m) * eta**(m-1) * z3_rk + bp_rk(k,m) * z3_m
477  END DO
478 
479  ord1_rk = 0.0
480  ord2_rk = 0.0
481  DO i = 1,ncomp
482  ord1_rk = ord1_rk + 2.0*mseg(k)*rho*x(i)*mseg(i)*sig_ij(i,k)**3 *uij(i,k)/t
483  ord2_rk = ord2_rk + 2.0*mseg(k)*rho*x(i)*mseg(i)*sig_ij(i,k)**3 *(uij(i,k)/t)**2
484  END DO
485 
486  c1_con= 1.0/ ( 1.0 + m_mean*(8.0*z3-2.0*z3*z3)/ome2/ome2 &
487  + (1.0 - m_mean)*(20.0*z3-27.0*z3*z3 +12.0*z3**3 -2.0*z3**4 ) &
488  /(ome*(2.0-z3))**2 )
489  c2_con= - c1_con*c1_con *( m_mean*(-4.0*z3*z3+20.0*z3+8.0)/ome2/ome3 &
490  + (1.0 - m_mean) *(2.0*z3**3 +12.0*z3*z3-48.0*z3+40.0) &
491  /(ome*(2.0-z3))**3 )
492  c1_rk= c2_con*z3_rk - c1_con*c1_con*m_rk * ( (8.0*z3-2.0*z3*z3)/ome2/ome2 &
493  - (-2.0*z3**4 +12.0*z3**3 -27.0*z3*z3+20.0*z3) / (ome*(2.0-z3))**2 )
494 
495  mdsp(k) = -2.0*pi* ( order1*rho*rho*i1_rk + ord1_rk*i1 ) &
496  - pi* c1_con*m_mean * ( order2*rho*rho*i2_rk + ord2_rk*i2 ) &
497  - pi* ( c1_con*m_rk + c1_rk*m_mean ) * order2*rho*rho*i2
498 
499 
500  !--------------------------------------------------------------------------
501  ! TPT-1-association according to Chapman et al.
502  !--------------------------------------------------------------------------
503  m_hbon(k) = 0.0
504  assoc = .false.
505  DO i = 1,ncomp
506  IF (nhb_typ(i) /= 0) assoc = .true.
507  END DO
508  IF (assoc) THEN
509 
510  ass_s2 = 0.0
511  DO kk = 1, nhb_typ(k)
512  ass_s2 = ass_s2 + nhb_no(k,kk) * log(mx(k,kk))
513  END DO
514 
515  m_hbon(k)=ass_s2
516  DO i = 1, ncomp
517  DO ii = 1, nhb_typ(i)
518  DO j = 1, ncomp
519  DO jj = 1, nhb_typ(j)
520  m_hbon(k) = m_hbon(k) - rho * rho / 2.0 * x(i) * x(j) * mx(i,ii) * mx(j,jj) &
521  * nhb_no(i,ii)*nhb_no(j,jj) *gij_rk(i,j) *ass_d(i,j,ii,jj)
522  END DO
523  END DO
524  END DO
525  END DO
526 
527  END IF
528 
529 
530  !--------------------------------------------------------------------------
531  ! polar terms
532  !--------------------------------------------------------------------------
533  CALL phi_polar ( k, z3_rk, fdd_rk, fqq_rk, fdq_rk )
534  my_dd(k) = fdd_rk
535  my_qq(k) = fqq_rk
536  my_dq(k) = fdq_rk
537 
538 
539  !--------------------------------------------------------------------------
540  ! d(f)/d(rho_k) : summation of all contributions
541  !--------------------------------------------------------------------------
542  myres(k) = mhs(k) +mhc(k) +mdsp(k) +m_hbon(k) +my_dd(k) +my_qq(k) +my_dq(k)
543 
544  END DO
545 
546 
547  !-----------------------------------------------------------------------------
548  ! finally calculate
549  ! mu_i^res(T,p,x)/kT = ln( phi_i ) when ensemble_flag = 'tp'
550  ! mu_i^res(T,rho,x)/kT when ensemble_flag = 'tv'
551  !-----------------------------------------------------------------------------
552 
553  DO k = 1, ncomp
554  ! write (*,*) k,myres(k) +LOG(rho*x(k)),rho*32000.0
555  IF (ensemble_flag == 'tp' ) lnphi(k) = myres(k) - log(zges)
556  IF (ensemble_flag == 'tv' ) lnphi(k) = myres(k)
557  ! write (*,*) 'in',k,lnphi(k),LOG(zges),eta
558  END DO
559  ! write (*,'(a,5G18.10)') 'fug.coeff.lnphi 1,2',lnphi(1),lnphi(2),rho
560 
561  dpdz = pgesdz
562  dpdz2 = pgesd2
563 
564  tfr= mhs(1)
565 
566 
567 END SUBROUTINE phi_eos
568 
569 
570 
571 
572 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
573 ! SUBROUTINE ddA_dhoi_drhoi_EOS
574 !
575 ! This subroutine gives the second derivatives
576 ! dd( F/VkT ) / d(rhoi)d(rhoi)
577 ! The variables are (T, rhoi), corrsponding to a case: ensemble_flag = 'tv'
578 ! The required input is (T, rho, x(nc)).
579 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
580 
581 SUBROUTINE dda_drhoi_drhoj_eos ( n_comp, rhoi, A_rr, Aig_rr )
582 
583  USE parameters
584  USE eos_variables
585  USE eos_constants
586  IMPLICIT NONE
587 
588  !-----------------------------------------------------------------------------
589  integer, intent(in) :: n_comp
590  real, dimension(n_comp), intent(in) :: rhoi
591  real, dimension(n_comp,n_comp), intent(out) :: a_rr
592  real, dimension(n_comp,n_comp), intent(out) :: aig_rr
593 
594  !-----------------------------------------------------------------------------
595  integer :: i, j, k, l, m
596  integer :: n_dim, m_dim
597  integer :: ii, jj, ll, iii, jjj, kk, lll
598  real :: z0, z1, z2, z3
599  real, allocatable, dimension(:) :: z0_r, z1_r, z2_r, z3_r
600  real :: m_mean
601  real :: ome, ome2, ome3, pi_6
602  real :: ahs_rkrl, ahc_rkrl, adsp_rkrl
603  real, allocatable, dimension(:) :: m_r
604  real :: a_term, b_term, m_rkrl
605  real, allocatable, dimension(:,:,:) :: gij_r
606  real, allocatable, dimension(:,:,:,:) :: gij_rr
607 
608  real, allocatable, dimension(:,:) :: q_xx, q_xr, q_xr_transpose, q_rr
609  real, allocatable, dimension(:,:) :: ahb_rr, a_polar_rr
610  real :: determinant
611 
612  real, allocatable, dimension(:,:) :: ap_r, bp_r
613  real, dimension(0:6) :: ap_rkrl, bp_rkrl
614  real :: i1, i2
615  real, allocatable, dimension(:) :: i1_r, i2_r
616  real :: i1_rkrl, i2_rkrl
617  real, allocatable, dimension(:) :: ord1_r, ord2_r
618  real :: ord1_rkrl, ord2_rkrl
619  real :: eta_m, aux_term_k
620  real :: c1_con, c2_con, c3_con, c2_dm
621  real, allocatable, dimension(:) :: c1_r
622  real :: c1_rkrl
623  real :: chi_dm, chi_dmdeta
624 
625  logical :: assoc
626  !-----------------------------------------------------------------------------
627 
628  allocate( z0_r(ncomp), z1_r(ncomp), z2_r(ncomp), z3_r(ncomp) )
629  allocate( gij_r(ncomp,ncomp,ncomp) )
630  allocate( gij_rr(ncomp,ncomp,ncomp,ncomp) )
631  allocate( m_r(ncomp) )
632  allocate( ap_r(ncomp,0:6), bp_r(ncomp,0:6) )
633  allocate( i1_r(ncomp), i2_r(ncomp) )
634  allocate( ord1_r(ncomp), ord2_r(ncomp) )
635  allocate( c1_r(ncomp) )
636  allocate( a_polar_rr( ncomp, ncomp ) )
637 
638  rho = sum( rhoi(1:ncomp) )
639  x(1:ncomp) = rhoi(1:ncomp) / rho
640 
641  !-----------------------------------------------------------------------------
642  ! obtain parameters and density independent expressions
643  !-----------------------------------------------------------------------------
644  CALL perturbation_parameter
645 
646  IF ( rho /= rho ) write (*,*) 'ddA_dhoi_drhoi_EOS: error in density',eta, z3t
647  IF ( rho /= rho ) stop
648  ! rhoi( 1:ncomp ) = x(1:ncomp ) * rho
649  z0 = z0t * rho
650  z1 = z1t * rho
651  z2 = z2t * rho
652  z3 = z3t * rho
653  eta = z3
654 
655  m_mean = z0t / (pi/6.0)
656  ome = 1.0 - z3
657  ome2 = ome * ome
658  ome3 = ome * ome2
659 
660  pi_6 = pi/6.0
661 
662  !=============================================================================
663  ! calculate the derivatives of f to rho_k ( d(f)/d(rho_k) )
664  !=============================================================================
665 
666  c1_con= 1.0/ ( 1.0 + m_mean*(8.0*z3-2.0*z3*z3)/ome2/ome2 &
667  + (1.0 - m_mean)*(20.0*z3-27.0*z3*z3 +12.0*z3**3 -2.0*z3**4 ) &
668  /(ome*(2.0-z3))**2 )
669  c2_con= - c1_con*c1_con *( m_mean*(-4.0*z3*z3+20.0*z3+8.0)/ome**5 &
670  + (1.0 - m_mean) *(2.0*z3**3 +12.0*z3*z3-48.0*z3+40.0) &
671  /(ome*(2.0-z3))**3 )
672  c3_con= 2.0 * c2_con*c2_con/c1_con - c1_con*c1_con &
673  *( m_mean*(-12.0*eta**2 +72.0*eta+60.0)/ome**6 &
674  + (1.0 - m_mean) &
675  *(-6.0*eta**4 -48.0*eta**3 +288.0*eta**2 &
676  -480.0*eta+264.0) /(ome*(2.0-eta))**4 )
677  chi_dm = (8.0*z3-2.0*z3*z3)/ome**4 &
678  - (-2.0*z3**4 +12.0*z3**3 -27.0*z3*z3+20.0*z3) / (ome*(2.0-z3))**2
679  chi_dmdeta = (-4.0*z3*z3+20.0*z3+8.0)/ome**5 &
680  - (2.0*z3**3 +12.0*z3*z3-48.0*z3+40.0) / (ome*(2.0-z3))**3
681  c2_dm = - 2.0*c2_con*c1_con*chi_dm - c1_con*c1_con*chi_dmdeta
682 
683  DO k = 1, ncomp
684 
685  z0_r(k) = pi_6 * mseg(k)
686  z1_r(k) = z0_r(k) * dhs(k)
687  z2_r(k) = z1_r(k) * dhs(k)
688  z3_r(k) = z2_r(k) * dhs(k)
689 
690  DO i = 1, ncomp
691  DO j = 1, ncomp
692  gij(i,j) = 1.0/ome + 3.0*dij_ab(i,j)*z2/ome2 + 2.0*(dij_ab(i,j)*z2)**2 /ome3
693  gij_r(k,i,j) = z3_r(k)/ome2 &
694  + 3.0*dij_ab(i,j)*(z2_r(k)+2.0*z2*z3_r(k)/ome)/ome2 &
695  + dij_ab(i,j)**2 *z2/ome3 *(4.0*z2_r(k)+6.0*z2*z3_r(k)/ome)
696  END DO
697  END DO
698 
699  m_r(k) = ( mseg(k) - m_mean ) / rho
700 
701  DO m = 0, 6
702  a_term = ( ap(m,2) + (3.0 -4.0/m_mean) *ap(m,3) ) / m_mean/m_mean
703  b_term = ( bp(m,2) + (3.0 -4.0/m_mean) *bp(m,3) ) / m_mean/m_mean
704  ap_r(k,m) = m_r(k) * a_term
705  bp_r(k,m) = m_r(k) * b_term
706  END DO
707 
708  i1 = 0.0
709  i2 = 0.0
710  i1_r(k) = 0.0
711  i2_r(k) = 0.0
712  DO m = 0, 6
713  eta_m = eta**m
714  i1 = i1 + apar(m) * eta_m
715  i2 = i2 + bpar(m) * eta_m
716  i1_r(k) = i1_r(k) + apar(m) * real(m) * eta**(m-1) * z3_r(k) + ap_r(k,m) * eta_m
717  i2_r(k) = i2_r(k) + bpar(m) * real(m) * eta**(m-1) * z3_r(k) + bp_r(k,m) * eta_m
718  END DO
719 
720  ord1_r(k) = 0.0
721  ord2_r(k) = 0.0
722  DO i = 1,ncomp
723  aux_term_k = 2.0*mseg(k)*mseg(i)*sig_ij(i,k)**3 *uij(i,k)/t
724  ord1_r(k) = ord1_r(k) + rhoi(i) * aux_term_k
725  ord2_r(k) = ord2_r(k) + rhoi(i) * aux_term_k * uij(i,k)/t
726  END DO
727 
728  c1_r(k)= c2_con*z3_r(k) - c1_con*c1_con * chi_dm * m_r(k)
729 
730  end do
731 
732  do k = 1, ncomp
733  do l = 1, ncomp
734  DO i = 1, ncomp
735  DO j = 1, ncomp
736  gij_rr(k,l,i,j) = 2.0*z3_r(k)*z3_r(l)/ome3 &
737  + 6.0*dij_ab(i,j)*( z2_r(k)*z3_r(l)+z2_r(l)*z3_r(k)+3.0*z2*z3_r(k)*z3_r(l)/ome ) /ome3 &
738  +dij_ab(i,j)**2/ome3 *( 4.0*z2_r(k)*z2_r(l) + 12.0*z2*(z2_r(k)*z3_r(l)+z2_r(l)*z3_r(k))/ome &
739  + 24.0*z2*z2*z3_r(k)*z3_r(l)/ome2 )
740  END DO
741  END DO
742  end do
743  end do
744 
745 
746  !=============================================================================
747  ! calculate the derivatives of f to rho_k and rho_l ( dd(f)/d(rho_k)d(rho_l) )
748  !=============================================================================
749 
750  DO k = 1, ncomp
751 
752  DO l = 1, ncomp
753 
754  !--------------------------------------------------------------------------
755  ! d(f)/d(rho_k) : hard sphere contribution
756  !--------------------------------------------------------------------------
757  if ( z3**4 > 0.0 ) then
758  ahs_rkrl = ( 3.0*(z1_r(k)*z2_r(l)+z1_r(l)*z2_r(k))/ome + 3.0*(z1_r(k)*z2+z1*z2_r(k))*z3_r(l)/ome2 &
759  + 3.0*(z1_r(l)*z2+z1*z2_r(l))*z3_r(k)/ome2 + 6.0*z1*z2*z3_r(k)*z3_r(l)/ome3 &
760  + 6.0*z2*z2_r(l)*z2_r(k)/z3/ome2 + 3.0*z2*z2*(z2_r(k)*z3_r(l)+z2_r(l)*z3_r(k))*(3.0*z3-1.0)/z3/z3/ome3 &
761  + 3.0*z2**3*z3_r(k)*z3_r(l) /z3/z3/ome3 &
762  + z2**3*z3_r(k)*z3_r(l)*(3.0*z3-1.0)*(5.0*z3-2.0)/z3**3/ome2/ome2 &
763  + ( 6.0*z2*z2*(z2_r(l)/z2*z2_r(k)*z3-z2_r(k)*z3_r(l)-z2_r(l)*z3_r(k)+z2*z3_r(k)*z3_r(l)/z3 )/z3**3 )*log(ome) &
764  + ( (2.0*z2*z3_r(k)-3.0*z2_r(k)*z3)*z2*z2/z3**3 + z0_r(k) )*z3_r(l)/ome &
765  + (z0_r(l)-z2*z2*(3.0*z2_r(l)*z3-2.0*z2*z3_r(l))/z3**3)*z3_r(k)/ome &
766  + (z0-z2**3/z3/z3)*z3_r(k)*z3_r(l)/ome2 ) / pi_6
767  else
768  ahs_rkrl = 0.0
769  end if
770 
771 
772  !--------------------------------------------------------------------------
773  ! d(f)/d(rho_k) : chain term
774  !--------------------------------------------------------------------------
775 
776  ahc_rkrl = 0.0
777  DO i = 1, ncomp
778  ahc_rkrl = ahc_rkrl + rhoi(i) * (1.0-mseg(i)) / gij(i,i) &
779  * ( gij_rr(k,l,i,i) - gij_r(k,i,i)*gij_r(l,i,i)/gij(i,i) )
780  END DO
781  ahc_rkrl = ahc_rkrl + (1.0-mseg(k)) / gij(k,k) * gij_r(l,k,k) + (1.0-mseg(l)) / gij(l,l) * gij_r(k,l,l)
782 
783 
784  !--------------------------------------------------------------------------
785  ! PC-SAFT: d(f)/d(rho_k) : dispersion contribution
786  !--------------------------------------------------------------------------
787 
788  m_rkrl = ( 2.0*m_mean - mseg(k) - mseg(l) ) / rho/rho
789 
790  ! --- derivatives of apar, bpar to rho_k -------------------------------
791  DO m = 0, 6
792  a_term = ( ap(m,2) + (3.0 -4.0/m_mean) *ap(m,3) ) / m_mean/m_mean
793  b_term = ( bp(m,2) + (3.0 -4.0/m_mean) *bp(m,3) ) / m_mean/m_mean
794  ap_rkrl(m) = m_rkrl * a_term + 2.0/m_mean*m_r(k)*m_r(l)*( 2.0*ap(m,3)/m_mean**3 - a_term )
795  bp_rkrl(m) = m_rkrl * b_term + 2.0/m_mean*m_r(k)*m_r(l)*( 2.0*bp(m,3)/m_mean**3 - b_term )
796  END DO
797 
798  i1_rkrl = 0.0
799  i2_rkrl = 0.0
800  DO m = 0, 6
801  eta_m = eta**m
802  i1_rkrl = i1_rkrl + ap_rkrl(m) * eta_m + real(m) *(z3_r(k)*ap_r(l,m)+z3_r(l)*ap_r(k,m)) *eta**(m-1) &
803  + real(m)*real(m-1)*apar(m)*z3_r(k)*z3_r(l)*eta**(m-2)
804  i2_rkrl = i2_rkrl + bp_rkrl(m) * eta_m + real(m) *(z3_r(k)*bp_r(l,m)+z3_r(l)*bp_r(k,m)) *eta**(m-1) &
805  + real(m)*real(m-1)*bpar(m)*z3_r(k)*z3_r(l)*eta**(m-2)
806  END DO
807 
808  ord1_rkrl = 2.0*mseg(k)*mseg(l)*sig_ij(k,l)**3 *uij(k,l)/t
809  ord2_rkrl = 2.0*mseg(k)*mseg(l)*sig_ij(k,l)**3 * (uij(k,l)/t)**2
810 
811  c1_rkrl = c3_con*z3_r(l)*z3_r(k) + c2_dm*m_r(l)*z3_r(k) &
812  - c1_con* ( 2.0*c1_r(l)*chi_dm*m_r(k) + c1_con*chi_dmdeta*z3_r(l)*m_r(k)+ c1_con*chi_dm*m_rkrl)
813 
814  adsp_rkrl = -2.0*pi* ( ord1_rkrl*i1 +ord1_r(k)*i1_r(l)+ord1_r(l)*i1_r(k)+order1*rho*rho*i1_rkrl ) &
815  - pi* ( m_rkrl*c1_con + m_r(k)*c1_r(l)+m_r(l)*c1_r(k)+ m_mean*c1_rkrl ) * order2*rho*rho*i2 &
816  - pi* (m_r(k)*c1_con+m_mean*c1_r(k)) * ( ord2_r(l)*i2 + order2*rho*rho*i2_r(l) ) &
817  - pi* (m_r(l)*c1_con+m_mean*c1_r(l)) * ( ord2_r(k)*i2 + order2*rho*rho*i2_r(k) ) &
818  - pi* m_mean*c1_con * ( ord2_rkrl*i2 +ord2_r(k)*i2_r(l)+ord2_r(l)*i2_r(k)+order2*rho*rho*i2_rkrl)
819 
820 
821  !--------------------------------------------------------------------------
822  ! second derivative of Helmholtz energy to rho_k and rho_l
823  !--------------------------------------------------------------------------
824 
825  a_rr(k,l) = ahs_rkrl + ahc_rkrl + adsp_rkrl
826  ! write (*,*) k, l, A_rr(k,l)
827 
828  end do
829  end do
830 
831  !-----------------------------------------------------------------------------
832  ! TPT-1-association according to Chapman et al.
833  !-----------------------------------------------------------------------------
834  assoc = .false.
835  DO i = 1,ncomp
836  IF (nhb_typ(i) /= 0) assoc = .true.
837  END DO
838  IF (assoc) THEN
839 
840  n_dim = 0
841  DO i = 1, ncomp
842  DO ii = 1, nhb_typ(i)
843  n_dim = n_dim + 1
844  END DO
845  END DO
846 
847  m_dim = ncomp
848 
849  allocate( q_xx( n_dim, n_dim ) )
850  allocate( q_xr( n_dim, m_dim ) )
851  allocate( q_xr_transpose( m_dim, n_dim ) )
852  allocate( q_rr( m_dim, m_dim ) )
853  allocate( ahb_rr( m_dim, m_dim ) )
854 
855 
856  iii = 0
857  DO i = 1, ncomp
858  DO ii = 1, nhb_typ(i)
859  iii = iii + 1
860  jjj = 0
861  DO j = 1, ncomp
862  DO jj = 1, nhb_typ(j)
863  jjj = jjj + 1
864  q_xx(iii,jjj) = - rhoi(i) * rhoi(j) * gij(i,j) *ass_d(i,j,ii,jj)
865  if ( iii == jjj ) q_xx(iii,jjj) = q_xx(iii,jjj) - rhoi(i) / ( mx(i,ii) * nhb_no(i,ii) )**2 * nhb_no(i,ii)
866  END DO
867  END DO
868  END DO
869  END DO
870 
871  do k = 1, ncomp
872  lll = 0
873  do l = 1, ncomp
874  DO ll = 1, nhb_typ(l)
875  lll = lll + 1
876  q_xr(lll,k) = 0.0
877  DO i = 1, ncomp
878  do ii = 1, nhb_typ(i)
879  q_xr(lll,k) = q_xr(lll,k) - rhoi(l)*rhoi(i)*mx(i,ii)*nhb_no(i,ii) *gij_r(k,i,l) *ass_d(i,l,ii,ll)
880  end do
881  END DO
882  do kk = 1, nhb_typ(k)
883  q_xr(lll,k) = q_xr(lll,k) - rhoi(l)* mx(k,kk)*nhb_no(k,kk) *gij(k,l) *ass_d(k,l,kk,ll)
884  end do
885  q_xr_transpose(k,lll) = q_xr(lll,k)
886  END DO
887  end do
888  end do
889 
890  q_rr(:,:) = 0.0
891  do k = 1, ncomp
892  do l = 1, ncomp
893  DO i = 1, ncomp
894  DO ii = 1, nhb_typ(i)
895  DO ll = 1, nhb_typ(l)
896  q_rr(k,l) = q_rr(k,l) - rhoi(i) * mx(i,ii)*nhb_no(i,ii) &
897  * mx(l,ll)*nhb_no(l,ll) * gij_r(k,i,l) *ass_d(i,l,ii,ll)
898  END DO
899  DO kk = 1, nhb_typ(k)
900  q_rr(k,l) = q_rr(k,l) - rhoi(i) * mx(i,ii)*nhb_no(i,ii) &
901  * mx(k,kk)*nhb_no(k,kk) * gij_r(l,i,k) *ass_d(i,k,ii,kk)
902  END DO
903  END DO
904  END DO
905  end do
906  end do
907 
908  do k = 1, ncomp
909  do l = 1, ncomp
910  DO kk = 1, nhb_typ(k)
911  DO ll = 1, nhb_typ(l)
912  q_rr(k,l) = q_rr(k,l) - mx(k,kk)*nhb_no(k,kk) * mx(l,ll)*nhb_no(l,ll) * gij(k,l) *ass_d(k,l,kk,ll)
913  end do
914  end do
915  DO i = 1, ncomp
916  DO ii = 1, nhb_typ(i)
917  DO j = 1, ncomp
918  DO jj = 1, nhb_typ(j)
919  q_rr(k,l) = q_rr(k,l) - 0.5 * rhoi(i) * rhoi(j) * mx(i,ii) * mx(j,jj) &
920  * nhb_no(i,ii)*nhb_no(j,jj) *gij_rr(k,l,i,j) *ass_d(i,j,ii,jj)
921  END DO
922  END DO
923  END DO
924  END DO
925  end do
926  end do
927 
928  !Ahb_rr = q_rr - q_Xr_transpose * inv( q_XX ) * q_Xr
929  call matinv( n_dim, m_dim, q_xx, q_xr, determinant ) ! output q_Xr := inv( q_XX ) * q_Xr
930 
931  ahb_rr = matmul( q_xr_transpose, q_xr )
932  ahb_rr = q_rr - ahb_rr
933 
934  do k = 1, ncomp
935  do l = 1, ncomp
936  a_rr(k,l) = a_rr(k,l) + ahb_rr(k,l)
937  end do
938  end do
939 
940  deallocate( q_xx, q_xr, q_xr_transpose, q_rr, ahb_rr )
941 
942  END IF
943 
944  !-----------------------------------------------------------------------------
945  ! polar terms
946  !-----------------------------------------------------------------------------
947  CALL a_polar_drhoi_drhoj ( n_comp, a_polar_rr )
948  a_rr(:,:) = a_rr(:,:) + a_polar_rr(:,:)
949 
950  !-----------------------------------------------------------------------------
951  ! ideal gas term
952  !-----------------------------------------------------------------------------
953  aig_rr(:,:) = 0.0
954  do k = 1, ncomp
955  aig_rr(k,k) = 1.0 / rhoi(k)
956  end do
957 
958  deallocate( gij_r )
959  deallocate( gij_rr )
960  deallocate( m_r )
961  deallocate( ap_r, bp_r )
962  deallocate( i1_r, i2_r )
963  deallocate( ord1_r, ord2_r )
964  deallocate( a_polar_rr )
965 
966 END SUBROUTINE dda_drhoi_drhoj_eos
967 
968 
969 
970 
971 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
972 ! SUBROUTINE P_EOS
973 !
974 ! calculates the pressure in units (Pa).
975 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
976 
977 SUBROUTINE p_eos
978 
979  !-----------------------------------------------------------------------------
980  USE parameters, ONLY: nc, nsite
981  USE eos_variables
982  USE eos_constants
983  USE eos_polar, ONLY: p_polar
984  IMPLICIT NONE
985 
986  !-----------------------------------------------------------------------------
987  INTEGER :: i, j, ii, jj, m
988  INTEGER :: ass_cnt,max_eval
989  LOGICAL :: assoc
990  REAL :: z0, z1, z2, z3
991  REAL :: ome, ome2, ome3, ome4, ome5, z3_m
992  REAL :: m_mean
993  REAL :: zges, zgesdz, zgesd2, zgesd3
994  REAL :: zhs, zhsdz, zhsd2, zhsd3
995  REAL :: zhc, zhcdz, zhcd2, zhcd3
996  REAL, DIMENSION(nc,nc) :: dgijdz, dgijd2, dgijd3, dgijd4
997  REAL :: zdsp, zdspdz, zdspd2, zdspd3
998  REAL :: c1_con, c2_con, c3_con, c4_con, c5_con
999  REAL :: i2, edi1dz, edi2dz, edi1d2, edi2d2
1000  REAL :: edi1d3, edi2d3, edi1d4, edi2d4
1001  REAL :: zhb, zhbdz, zhbd2, zhbd3
1002  REAL, DIMENSION(nc,nc,nsite,nsite) :: delta, dq_dz, dq_d2, dq_d3, dq_d4
1003  REAL, DIMENSION(nc,nsite) :: mx_itr, dmx_dz, ndmxdz, dmx_d2, ndmxd2
1004  REAL, DIMENSION(nc,nsite) :: dmx_d3, ndmxd3, dmx_d4, ndmxd4
1005  REAL :: err_sum, sum0, sum1, sum2, sum3, sum4, attenu, tol
1006  REAL :: sum_d1, sum_d2, sum_d3, sum_d4
1007  REAL :: zdd, zddz, zddz2, zddz3
1008  REAL :: zqq, zqqz, zqqz2, zqqz3
1009  REAL :: zdq, zdqz, zdqz2, zdqz3
1010  !-----------------------------------------------------------------------------
1011 
1012 
1013  !-----------------------------------------------------------------------------
1014  ! abbreviations
1015  !-----------------------------------------------------------------------------
1016  rho = eta/z3t
1017  z0 = z0t*rho
1018  z1 = z1t*rho
1019  z2 = z2t*rho
1020  z3 = z3t*rho
1021 
1022  m_mean = z0t/(pi/6.0)
1023  ome = 1.0 -eta
1024  ome2 = ome * ome
1025  ome3 = ome2 * ome
1026  ome4 = ome2 * ome2
1027  ome5 = ome4 * ome
1028 
1029  ! m_mean2=0.0
1030  ! lij(1,2)= -0.050
1031  ! lij(2,1)=lij(1,2)
1032  ! DO i =1,ncomp
1033  ! DO j =1,ncomp
1034  ! m_mean2=m_mean2+x(i)*x(j) * (mseg(i)+mseg(j))/2.0*(1.0-lij(i,j))
1035  ! ENDDO
1036  ! ENDDO
1037 
1038 
1039  !-----------------------------------------------------------------------------
1040  ! radial distr. function at contact, gij , and derivatives
1041  ! dgijdz=d(gij)/d(eta) and dgijd2 = dd(gij)/d(eta)**2
1042  !-----------------------------------------------------------------------------
1043  DO i = 1, ncomp
1044  DO j=1,ncomp
1045  ! j=i
1046  gij(i,j) = 1.0/ome + 3.0*dij_ab(i,j)*z2/ome2 + 2.0*(dij_ab(i,j)*z2)**2 /ome3
1047  dgijdz(i,j)= 1.0/ome2 + 3.0*dij_ab(i,j)*z2*(1.0+z3)/z3/ome3 &
1048  + (dij_ab(i,j)*z2/ome2)**2 *(4.0+2.0*z3)/z3
1049  dgijd2(i,j) = 2.0/ome3 + 6.0*dij_ab(i,j)*z2/z3/ome4 *(2.0+z3) &
1050  + (2.0*dij_ab(i,j)*z2/z3)**2 /ome5 *(1.0+4.0*z3+z3*z3)
1051  dgijd3(i,j) = 6.0/ome4 + 18.0*dij_ab(i,j)*z2/z3/ome5 *(3.0+z3) &
1052  + 12.0*(dij_ab(i,j)*z2/z3/ome3 )**2 *(3.0+6.0*z3+z3*z3)
1053  dgijd4(i,j) = 24.0/ome5 + 72.0*dij_ab(i,j)*z2/z3/ome**6 *(4.0+z3) &
1054  + 48.0*(dij_ab(i,j)*z2/z3)**2 /ome**7 *(6.0+8.0*z3+z3*z3)
1055  END DO
1056  END DO
1057 
1058 
1059  !-----------------------------------------------------------------------------
1060  ! p : hard sphere contribution
1061  !-----------------------------------------------------------------------------
1062  zhs = m_mean* ( z3/ome + 3.0*z1*z2/z0/ome2 + z2**3 /z0*(3.0-z3)/ome3 )
1063  zhsdz = m_mean*( 1.0/ome2 + 3.0*z1*z2/z0/z3*(1.0+z3)/ome3 &
1064  + 6.0*z2**3 /z0/z3/ome4 )
1065  zhsd2 = m_mean*( 2.0/ome3 + 6.0*z1*z2/z0/z3*(2.0+z3)/ome4 &
1066  + 6.0*z2**3 /z0/z3/z3*(1.0+3.0*z3)/ome5 )
1067  zhsd3 = m_mean*( 6.0/ome4 + 18.0*z1*z2/z0/z3*(3.0+z3)/ome5 &
1068  + 24.0*z2**3 /z0/z3/z3*(2.0+3.0*z3)/ome**6 )
1069 
1070 
1071  !-----------------------------------------------------------------------------
1072  ! p : chain term
1073  !-----------------------------------------------------------------------------
1074  zhc = 0.0
1075  zhcdz = 0.0
1076  zhcd2 = 0.0
1077  zhcd3 = 0.0
1078  DO i= 1, ncomp
1079  zhc = zhc + x(i)*(1.0-mseg(i))*eta/gij(i,i)* dgijdz(i,i)
1080  zhcdz = zhcdz + x(i)*(1.0-mseg(i)) *(-eta*(dgijdz(i,i)/gij(i,i))**2 &
1081  + dgijdz(i,i)/gij(i,i) + eta/gij(i,i)*dgijd2(i,i))
1082  zhcd2 = zhcd2 + x(i)*(1.0-mseg(i)) &
1083  *( 2.0*eta*(dgijdz(i,i)/gij(i,i))**3 &
1084  -2.0*(dgijdz(i,i)/gij(i,i))**2 &
1085  -3.0*eta/gij(i,i)**2 *dgijdz(i,i)*dgijd2(i,i) &
1086  +2.0/gij(i,i)*dgijd2(i,i) +eta/gij(i,i)*dgijd3(i,i) )
1087  zhcd3 = zhcd3 + x(i)*(1.0-mseg(i)) *( 6.0*(dgijdz(i,i)/gij(i,i))**3 &
1088  -6.0*eta*(dgijdz(i,i)/gij(i,i))**4 &
1089  +12.0*eta/gij(i,i)**3 *dgijdz(i,i)**2 *dgijd2(i,i) &
1090  -9.0/gij(i,i)**2 *dgijdz(i,i)*dgijd2(i,i) +3.0/gij(i,i)*dgijd3(i,i) &
1091  -3.0*eta*(dgijd2(i,i)/gij(i,i))**2 &
1092  -4.0*eta/gij(i,i)**2 *dgijdz(i,i)*dgijd3(i,i) &
1093  +eta/gij(i,i)*dgijd4(i,i) )
1094  END DO
1095 
1096  !-----------------------------------------------------------------------------
1097  ! p : PC-SAFT dispersion contribution
1098  ! note: edI1dz is equal to d(eta*I1)/d(eta), analogous for edI2dz
1099  !-----------------------------------------------------------------------------
1100  i2 = 0.0
1101  edi1dz = 0.0
1102  edi2dz = 0.0
1103  edi1d2 = 0.0
1104  edi2d2 = 0.0
1105  edi1d3 = 0.0
1106  edi2d3 = 0.0
1107  edi1d4 = 0.0
1108  edi2d4 = 0.0
1109  DO m = 0, 6
1110  z3_m = z3**m
1111  i2 = i2 + bpar(m) * z3_m
1112  edi1dz= edi1dz + apar(m) * REAL(m+1) * z3_m
1113  edi2dz= edi2dz + bpar(m) * REAL(m+1) * z3_m
1114  edi1d2= edi1d2 + apar(m) * REAL((m+1)*m) * z3**(m-1)
1115  edi2d2= edi2d2 + bpar(m) * REAL((m+1)*m) * z3**(m-1)
1116  edi1d3= edi1d3 + apar(m) * REAL((m+1)*m*(m-1)) * z3**(m-2)
1117  edi2d3= edi2d3 + bpar(m) * REAL((m+1)*m*(m-1)) * z3**(m-2)
1118  edi1d4= edi1d4 + apar(m) * REAL((m+1)*m*(m-1)*(m-2)) * z3**(m-3)
1119  edi2d4= edi2d4 + bpar(m) * REAL((m+1)*m*(m-1)*(m-2)) * z3**(m-3)
1120  END DO
1121 
1122  c1_con= 1.0/ ( 1.0 + m_mean*(8.0*eta-2.0*eta**2 )/ome4 &
1123  + (1.0 - m_mean)*(20.0*eta-27.0*eta**2 &
1124  + 12.0*eta**3 -2.0*eta**4 ) /(ome*(2.0-eta))**2 )
1125  c2_con= - c1_con*c1_con &
1126  *(m_mean*(-4.0*eta**2 +20.0*eta+8.0)/ome5 + (1.0 - m_mean) &
1127  *(2.0*eta**3 +12.0*eta**2 -48.0*eta+40.0) &
1128  /(ome*(2.0-eta))**3 )
1129  c3_con= 2.0 * c2_con*c2_con/c1_con - c1_con*c1_con &
1130  *( m_mean*(-12.0*eta**2 +72.0*eta+60.0)/ome**6 &
1131  + (1.0 - m_mean) &
1132  *(-6.0*eta**4 -48.0*eta**3 +288.0*eta**2 &
1133  -480.0*eta+264.0) /(ome*(2.0-eta))**4 )
1134  c4_con= 6.0*c2_con*c3_con/c1_con -6.0*c2_con**3 /c1_con**2 &
1135  - c1_con*c1_con &
1136  *( m_mean*(-48.0*eta**2 +336.0*eta+432.0)/ome**7 &
1137  + (1.0 - m_mean) &
1138  *(24.0*eta**5 +240.0*eta**4 -1920.0*eta**3 &
1139  +4800.0*eta**2 -5280.0*eta+2208.0) /(ome*(2.0-eta))**5 )
1140  c5_con= 6.0*c3_con**2 /c1_con - 36.0*c2_con**2 /c1_con**2 *c3_con &
1141  + 8.0*c2_con/c1_con*c4_con+24.0*c2_con**4 /c1_con**3 &
1142  - c1_con*c1_con &
1143  *( m_mean*(-240.0*eta**2 +1920.0*eta+3360.0)/ome**8 &
1144  + (1.0 - m_mean) &
1145  *(-120.0*eta**6 -1440.0*eta**5 +14400.0*eta**4 &
1146  -48000.0*eta**3 +79200.0*eta**2 -66240.0*eta+22560.0) &
1147  /(ome*(2.0-eta))**6 )
1148 
1149  zdsp = - 2.0*pi*rho*edi1dz*order1 &
1150  - pi*rho*order2*m_mean*(c2_con*i2*eta + c1_con*edi2dz)
1151  zdspdz= zdsp/eta - 2.0*pi*rho*edi1d2*order1 &
1152  - pi*rho*order2*m_mean*(c3_con*i2*eta + 2.0*c2_con*edi2dz + c1_con*edi2d2)
1153  zdspd2= -2.0*zdsp/eta/eta +2.0*zdspdz/eta &
1154  - 2.0*pi*rho*edi1d3*order1 - pi*rho*order2*m_mean*(c4_con*i2*eta &
1155  + 3.0*c3_con*edi2dz +3.0*c2_con*edi2d2 +c1_con*edi2d3)
1156  zdspd3= 6.0*zdsp/eta**3 -6.0*zdspdz/eta/eta &
1157  + 3.0*zdspd2/eta - 2.0*pi*rho*edi1d4*order1 &
1158  - pi*rho*order2*m_mean*(c5_con*i2*eta &
1159  + 4.0*c4_con*edi2dz +6.0*c3_con*edi2d2 &
1160  + 4.0*c2_con*edi2d3 + c1_con*edi2d4)
1161 
1162 
1163 
1164  !-----------------------------------------------------------------------------
1165  ! p: TPT-1-association accord. to Chapman et al.
1166  !-----------------------------------------------------------------------------
1167  zhb = 0.0
1168  zhbdz = 0.0
1169  zhbd2 = 0.0
1170  zhbd3 = 0.0
1171  assoc = .false.
1172  DO i = 1,ncomp
1173  IF ( nhb_typ(i) /= 0 ) assoc = .true.
1174  END DO
1175  IF (assoc) THEN
1176 
1177  DO i = 1, ncomp
1178  DO ii = 1, nhb_typ(i)
1179  DO j = 1, ncomp
1180  DO jj = 1, nhb_typ(j)
1181  delta(i,j,ii,jj) = gij(i,j) * ass_d(i,j,ii,jj)
1182  dq_dz(i,j,ii,jj) = dgijdz(i,j) * ass_d(i,j,ii,jj)
1183  dq_d2(i,j,ii,jj) = dgijd2(i,j) * ass_d(i,j,ii,jj)
1184  dq_d3(i,j,ii,jj) = dgijd3(i,j) * ass_d(i,j,ii,jj)
1185  dq_d4(i,j,ii,jj) = dgijd4(i,j) * ass_d(i,j,ii,jj)
1186  END DO
1187  END DO
1188  END DO
1189  END DO
1190 
1191  ! --- constants for iteration ---------------------------------------------
1192  attenu = 0.7
1193  tol = 1.e-10
1194  IF ( eta < 0.2 ) tol = 1.e-12
1195  IF ( eta < 0.01 ) tol = 1.e-13
1196  IF ( eta < 1.e-6) tol = 1.e-15
1197  max_eval = 1000
1198 
1199  ! --- initialize mx(i,j) --------------------------------------------------
1200  DO i = 1, ncomp
1201  DO ii = 1, nhb_typ(i)
1202  mx(i,ii) = 1.0
1203  dmx_dz(i,ii) = 0.0
1204  dmx_d2(i,ii) = 0.0
1205  dmx_d3(i,ii) = 0.0
1206  dmx_d4(i,ii) = 0.0
1207  END DO
1208  END DO
1209 
1210  ! --- iterate over all components and all sites ---------------------------
1211  ass_cnt = 0
1212  err_sum = tol + 1.0
1213  DO WHILE ( err_sum > tol .AND. ass_cnt <= max_eval)
1214  ass_cnt = ass_cnt + 1
1215  DO i = 1, ncomp
1216  DO ii = 1, nhb_typ(i)
1217  sum0 = 0.0
1218  sum1 = 0.0
1219  sum2 = 0.0
1220  sum3 = 0.0
1221  sum4 = 0.0
1222  DO j = 1, ncomp
1223  DO jj = 1, nhb_typ(j)
1224  sum0 =sum0 +x(j)*nhb_no(j,jj)* mx(j,jj)* delta(i,j,ii,jj)
1225  sum1 =sum1 +x(j)*nhb_no(j,jj)*( mx(j,jj)* dq_dz(i,j,ii,jj) &
1226  + dmx_dz(j,jj)* delta(i,j,ii,jj))
1227  sum2 =sum2 +x(j)*nhb_no(j,jj)*( mx(j,jj)* dq_d2(i,j,ii,jj) &
1228  + 2.0*dmx_dz(j,jj)* dq_dz(i,j,ii,jj) &
1229  + dmx_d2(j,jj)* delta(i,j,ii,jj))
1230  sum3 =sum3 +x(j)*nhb_no(j,jj)*( mx(j,jj)* dq_d3(i,j,ii,jj) &
1231  + 3.0*dmx_dz(j,jj)* dq_d2(i,j,ii,jj) &
1232  + 3.0*dmx_d2(j,jj)* dq_dz(i,j,ii,jj) &
1233  + dmx_d3(j,jj)* delta(i,j,ii,jj))
1234  sum4 =sum4 + x(j)*nhb_no(j,jj)*( mx(j,jj)* dq_d4(i,j,ii,jj) &
1235  + 4.0*dmx_dz(j,jj)* dq_d3(i,j,ii,jj) &
1236  + 6.0*dmx_d2(j,jj)* dq_d2(i,j,ii,jj) &
1237  + 4.0*dmx_d3(j,jj)* dq_dz(i,j,ii,jj) &
1238  + dmx_d4(j,jj)* delta(i,j,ii,jj))
1239  END DO
1240  END DO
1241  mx_itr(i,ii)= 1.0 / (1.0 + sum0 * rho)
1242  ndmxdz(i,ii)= -(mx_itr(i,ii)*mx_itr(i,ii))* (sum0/z3t +sum1*rho)
1243  ndmxd2(i,ii)= + 2.0/mx_itr(i,ii)*ndmxdz(i,ii)*ndmxdz(i,ii) &
1244  - (mx_itr(i,ii)*mx_itr(i,ii)) * (2.0/z3t*sum1 + rho*sum2)
1245  ndmxd3(i,ii)= - 6.0/mx_itr(i,ii)**2 *ndmxdz(i,ii)**3 &
1246  + 6.0/mx_itr(i,ii)*ndmxdz(i,ii)*ndmxd2(i,ii) - mx_itr(i,ii)*mx_itr(i,ii) &
1247  * (3.0/z3t*sum2 + rho*sum3)
1248  ndmxd4(i,ii)= 24.0/mx_itr(i,ii)**3 *ndmxdz(i,ii)**4 &
1249  -36.0/mx_itr(i,ii)**2 *ndmxdz(i,ii)**2 *ndmxd2(i,ii) &
1250  +6.0/mx_itr(i,ii)*ndmxd2(i,ii)**2 &
1251  +8.0/mx_itr(i,ii)*ndmxdz(i,ii)*ndmxd3(i,ii) - mx_itr(i,ii)**2 &
1252  *(4.0/z3t*sum3 + rho*sum4)
1253  END DO
1254  END DO
1255 
1256  err_sum = 0.0
1257  DO i = 1, ncomp
1258  DO ii = 1, nhb_typ(i)
1259  err_sum = err_sum + abs(mx_itr(i,ii) - mx(i,ii)) &
1260  + abs(ndmxdz(i,ii) - dmx_dz(i,ii)) + abs(ndmxd2(i,ii) - dmx_d2(i,ii))
1261  mx(i,ii) = mx_itr(i,ii)*attenu + mx(i,ii) * (1.0-attenu)
1262  dmx_dz(i,ii) = ndmxdz(i,ii)*attenu + dmx_dz(i,ii) * (1.0-attenu)
1263  dmx_d2(i,ii) = ndmxd2(i,ii)*attenu + dmx_d2(i,ii) * (1.0-attenu)
1264  dmx_d3(i,ii) = ndmxd3(i,ii)*attenu + dmx_d3(i,ii) * (1.0-attenu)
1265  dmx_d4(i,ii) = ndmxd4(i,ii)*attenu + dmx_d4(i,ii) * (1.0-attenu)
1266  END DO
1267  END DO
1268  END DO
1269 
1270  IF ( ass_cnt >= max_eval .AND. err_sum > sqrt(tol) ) THEN
1271  WRITE (*,'(a,2G15.7)') 'P_EOS: Max_eval violated (mx) Err_Sum= ',err_sum,tol
1272  ! stop
1273  END IF
1274 
1275 
1276  ! --- calculate the hydrogen-bonding contribution -------------------------
1277  DO i = 1, ncomp
1278  sum_d1 = 0.0
1279  sum_d2 = 0.0
1280  sum_d3 = 0.0
1281  sum_d4 = 0.0
1282  DO ii = 1, nhb_typ(i)
1283  sum_d1= sum_d1 +nhb_no(i,ii)* dmx_dz(i,ii)*(1.0/mx(i,ii)-0.5)
1284  sum_d2= sum_d2 +nhb_no(i,ii)*(dmx_d2(i,ii)*(1.0/mx(i,ii)-0.5) &
1285  -(dmx_dz(i,ii)/mx(i,ii))**2 )
1286  sum_d3= sum_d3 +nhb_no(i,ii)*(dmx_d3(i,ii)*(1.0/mx(i,ii)-0.5) &
1287  -3.0/mx(i,ii)**2 *dmx_dz(i,ii)*dmx_d2(i,ii) + 2.0*(dmx_dz(i,ii)/mx(i,ii))**3 )
1288  sum_d4= sum_d4 +nhb_no(i,ii)*(dmx_d4(i,ii)*(1.0/mx(i,ii)-0.5) &
1289  -4.0/mx(i,ii)**2 *dmx_dz(i,ii)*dmx_d3(i,ii) &
1290  + 12.0/mx(i,ii)**3 *dmx_dz(i,ii)**2 *dmx_d2(i,ii) &
1291  - 3.0/mx(i,ii)**2 *dmx_d2(i,ii)**2 - 6.0*(dmx_dz(i,ii)/mx(i,ii))**4 )
1292  END DO
1293  zhb = zhb + x(i) * eta * sum_d1
1294  zhbdz = zhbdz + x(i) * eta * sum_d2
1295  zhbd2 = zhbd2 + x(i) * eta * sum_d3
1296  zhbd3 = zhbd3 + x(i) * eta * sum_d4
1297  END DO
1298  zhbdz = zhbdz + zhb/eta
1299  zhbd2 = zhbd2 + 2.0/eta*zhbdz-2.0/eta**2 *zhb
1300  zhbd3 = zhbd3 - 6.0/eta**2 *zhbdz+3.0/eta*zhbd2 + 6.0/eta**3 *zhb
1301  END IF
1302 
1303 
1304  !-----------------------------------------------------------------------------
1305  ! p: polar terms
1306  !-----------------------------------------------------------------------------
1307  CALL p_polar ( zdd, zddz, zddz2, zddz3, zqq, zqqz, zqqz2, zqqz3, zdq, zdqz, zdqz2, zdqz3 )
1308 
1309 
1310  !-----------------------------------------------------------------------------
1311  ! compressibility factor z and total p
1312  ! as well as derivatives d(z)/d(eta) and d(p)/d(eta) with unit [Pa]
1313  !-----------------------------------------------------------------------------
1314  zges = 1.0 + zhs + zhc + zdsp + zhb + zdd + zqq + zdq
1315  zgesdz = zhsdz + zhcdz + zdspdz + zhbdz + zddz + zqqz + zdqz
1316  zgesd2 = zhsd2 + zhcd2 + zdspd2 + zhbd2 + zddz2 +zqqz2+zdqz2
1317  zgesd3 = zhsd3 + zhcd3 + zdspd3 + zhbd3 + zddz3 +zqqz3+zdqz3
1318 
1319  pges = zges *rho *(kbol*t)/1.e-30
1320  pgesdz = ( zgesdz*rho + zges*rho/z3 )*(kbol*t)/1.e-30
1321  pgesd2 = ( zgesd2*rho + 2.0*zgesdz*rho/z3 )*(kbol*t)/1.e-30
1322  pgesd3 = ( zgesd3*rho + 3.0*zgesd2*rho/z3 )*(kbol*t)/1.e-30
1323 
1324 END SUBROUTINE p_eos
1325 
1326 
1327 
1328 
1329 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1330 ! SUBROUTINE F_EOS
1331 !
1332 ! calculates the Helmholtz energy f/kT. The input to the subroutine is
1333 ! (T,eta,x), where eta is the packing fraction.
1334 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1335 
1336 SUBROUTINE f_eos
1337 
1338  USE parameters, ONLY: nc, nsite
1339  USE eos_variables
1340  USE eos_constants
1341  USE eos_polar, ONLY: f_polar
1342  IMPLICIT NONE
1343 
1344  !-----------------------------------------------------------------------------
1345  INTEGER :: i, j, ii, jj, m
1346  REAL :: z0, z1, z2, z3
1347  REAL :: ome, ome2, ome3, m_mean ! ,lij(nc,nc)
1348  REAL :: i1, i2, c1_con
1349  REAL :: fhs, fdsp, fhc
1350 
1351  LOGICAL :: assoc
1352  INTEGER :: ass_cnt,max_eval
1353  REAL :: delta(nc,nc,nsite,nsite)
1354  REAL :: mx_itr(nc,nsite), err_sum, sum, attenu, tol, fhb
1355  REAL :: ass_s1, ass_s2
1356 
1357  REAL :: fdd, fqq, fdq
1358  !-----------------------------------------------------------------------------
1359 
1360 
1361  !-----------------------------------------------------------------------------
1362  ! abbreviations
1363  !-----------------------------------------------------------------------------
1364  rho = eta/z3t
1365  IF ( rho /= rho ) write (*,*) 'F_EOS: error in density',eta, z3t
1366  z0 = z0t*rho
1367  z1 = z1t*rho
1368  z2 = z2t*rho
1369  z3 = z3t*rho
1370 
1371  m_mean = z0t / ( pi / 6.0 )
1372  ome = 1.0 - eta
1373  ome2 = ome * ome
1374  ome3 = ome * ome2
1375 
1376  ! m_mean2 = 0.0
1377  ! lij(1,2) = -0.05
1378  ! lij(2,1) = lij(1,2)
1379  ! DO i = 1, ncomp
1380  ! DO j = 1, ncomp
1381  ! m_mean2 = m_mean2 + x(i)*x(j)*(mseg(i)+mseg(j))/2.0*(1.0-lij(i,j))
1382  ! ENDDO
1383  ! ENDDO
1384 
1385 
1386  !-----------------------------------------------------------------------------
1387  ! radial distr. function at contact, gij
1388  !-----------------------------------------------------------------------------
1389  DO i = 1, ncomp
1390  DO j=1,ncomp
1391  gij(i,j) = 1.0/ome + 3.0*dij_ab(i,j)*z2/ome2 + 2.0*(dij_ab(i,j)*z2)**2 /ome3
1392  END DO
1393  END DO
1394 
1395 
1396  !-----------------------------------------------------------------------------
1397  ! Helmholtz energy : hard sphere contribution
1398  !-----------------------------------------------------------------------------
1399  fhs= m_mean*( 3.0*z1*z2/ome + z2**3 /z3/ome2 + (z2**3 /z3/z3-z0)*log(ome) )/z0
1400 
1401 
1402  !-----------------------------------------------------------------------------
1403  ! Helmholtz energy : chain term
1404  !-----------------------------------------------------------------------------
1405  fhc = 0.0
1406  DO i = 1, ncomp
1407  fhc = fhc + x(i) *(1.0- mseg(i)) *log(gij(i,i))
1408  END DO
1409 
1410 
1411  !-----------------------------------------------------------------------------
1412  ! Helmholtz energy : PC-SAFT dispersion contribution
1413  !-----------------------------------------------------------------------------
1414 
1415  i1 = 0.0
1416  i2 = 0.0
1417  DO m = 0, 6
1418  i1 = i1 + apar(m) * eta**m
1419  i2 = i2 + bpar(m) * eta**m
1420  END DO
1421 
1422  c1_con= 1.0/ ( 1.0 + m_mean*(8.0*eta-2.0*eta**2 )/ome**4 &
1423  + (1.0 - m_mean)*(20.0*eta-27.0*eta**2 + 12.0*eta**3 -2.0*eta**4 ) /(ome*(2.0-eta))**2 )
1424 
1425  fdsp = -2.0*pi * rho * i1 * order1 - pi * rho * c1_con * m_mean * i2 * order2
1426 
1427 
1428 
1429  !-----------------------------------------------------------------------------
1430  ! TPT-1-association according to Chapman et al.
1431  !-----------------------------------------------------------------------------
1432  fhb = 0.0
1433  assoc = .false.
1434  DO i = 1, ncomp
1435  IF (nhb_typ(i) /= 0) assoc = .true.
1436  END DO
1437  IF (assoc) THEN
1438 
1439  DO i = 1, ncomp
1440  DO ii = 1, nhb_typ(i)
1441  IF (mx(i,ii) == 0.0) mx(i,ii) = 1.0 ! Initialize mx(i,j)
1442  DO j = 1, ncomp
1443  DO jj = 1, nhb_typ(j)
1444  delta(i,j,ii,jj) = gij(i,j) * ass_d(i,j,ii,jj)
1445  END DO
1446  END DO
1447  END DO
1448  END DO
1449 
1450 
1451  ! --- constants for iteration ---------------------------------------------
1452  attenu = 0.70
1453  tol = 1.e-10
1454  IF (eta < 0.2) tol = 1.e-12
1455  IF (eta < 0.01) tol = 1.e-13
1456  max_eval = 200
1457  err_sum = 2.0 * tol
1458 
1459  ! --- iterate over all components and all sites ---------------------------
1460  ass_cnt = 0
1461  DO WHILE ( err_sum > tol .AND. ass_cnt <= max_eval )
1462 
1463  ass_cnt = ass_cnt + 1
1464 
1465  DO i = 1, ncomp
1466  DO ii = 1, nhb_typ(i)
1467  sum = 0.0
1468  DO j = 1, ncomp
1469  DO jj = 1, nhb_typ(j)
1470  sum = sum + x(j)* mx(j,jj)*nhb_no(j,jj) *delta(i,j,ii,jj)
1471  ! if (ass_cnt == 1) write (*,*) j,jj,x(j), mx(j,jj)
1472  END DO
1473  END DO
1474  mx_itr(i,ii) = 1.0 / (1.0 + sum * rho)
1475  ! if (ass_cnt == 1) write (*,*) 'B',ass_cnt,sum, rho
1476  END DO
1477  END DO
1478 
1479  err_sum = 0.0
1480  DO i = 1, ncomp
1481  DO ii = 1, nhb_typ(i)
1482  err_sum = err_sum + abs(mx_itr(i,ii) - mx(i,ii)) ! / ABS(mx_itr(i,ii))
1483  mx(i,ii) = mx_itr(i,ii) * attenu + mx(i,ii) * (1.0 - attenu)
1484  END DO
1485  END DO
1486 
1487  END DO
1488 
1489  IF ( err_sum /= err_sum ) write (*,*) 'F_EOS: association "not a number"',ass_cnt, rho, sum
1490  IF ( ass_cnt >= max_eval ) THEN
1491  WRITE(*,'(a,2G15.7)') 'F_EOS: Max_eval violated (mx) Err_Sum = ',err_sum,tol
1492  stop
1493  END IF
1494 
1495 
1496  DO i = 1, ncomp
1497  ass_s1 = 0.0
1498  ass_s2 = 0.0
1499  DO ii = 1, nhb_typ(i)
1500  ass_s1 = ass_s1 + nhb_no(i,ii) * ( 1.0 - mx(i,ii) )
1501  ass_s2 = ass_s2 + nhb_no(i,ii) * log( mx(i,ii) )
1502  END DO
1503  fhb = fhb + x(i) * ( ass_s2 + ass_s1 / 2.0 )
1504  END DO
1505 
1506  END IF
1507 
1508 
1509  !-----------------------------------------------------------------------------
1510  ! polar terms
1511  !-----------------------------------------------------------------------------
1512  CALL f_polar ( fdd, fqq, fdq )
1513 
1514 
1515  !-----------------------------------------------------------------------------
1516  ! resid. Helmholtz energy f/kT
1517  !-----------------------------------------------------------------------------
1518  fres = fhs + fhc + fdsp + fhb + fdd + fqq + fdq
1519 
1520  tfr= fres
1521 
1522 END SUBROUTINE f_eos
1523 
1524 
1525 
1526 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1527 ! SUBROUTINE PERTURBATION_PARAMETER
1528 !
1529 ! calculates density-independent parameters of the equation of state.
1530 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1531 
1532 SUBROUTINE perturbation_parameter
1533 
1534  USE parameters, ONLY: pi, kbol, rgas, nav, tau
1535  USE eos_variables
1536  USE eos_constants
1537  IMPLICIT NONE
1538 
1539  !-----------------------------------------------------------------------------
1540  INTEGER :: i, j, k, l, m, no
1541  LOGICAL :: assoc, qudpole, dipole
1542  REAL :: m_mean
1543  REAL, DIMENSION(nc) :: d00, u
1544  REAL, DIMENSION(nc,nc,nsite,nsite) :: eps_hb
1545  REAL, DIMENSION(nc,nc) :: kap_hb
1546  REAL :: eps_kij, k_kij
1547  !-----------------------------------------------------------------------------
1548 
1549  !-----------------------------------------------------------------------------
1550  ! pure component parameters
1551  !-----------------------------------------------------------------------------
1552  DO i = 1, ncomp
1553  u(i) = parame(i,3)
1554  mseg(i)= parame(i,1)
1555  dhs(i) = parame(i,2) * ( 1.0 - 0.12 * exp( -3.0 * parame(i,3) / t ) )
1556  d00(i) = parame(i,2)
1557  END DO
1558 
1559 
1560  !-----------------------------------------------------------------------------
1561  ! combination rules
1562  !-----------------------------------------------------------------------------
1563  DO i = 1, ncomp
1564  DO j = 1, ncomp
1565  sig_ij(i,j) = 0.5 * ( d00(i) + d00(j) )
1566  uij(i,j) = ( 1.0 - kij(i,j) ) * ( u(i)*u(j) )**0.5
1567  END DO
1568  END DO
1569 
1570 
1571  !-----------------------------------------------------------------------------
1572  ! abbreviations
1573  !-----------------------------------------------------------------------------
1574  z0t = pi / 6.0 * sum( x(1:ncomp) * mseg(1:ncomp) )
1575  z1t = pi / 6.0 * sum( x(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp) )
1576  z2t = pi / 6.0 * sum( x(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**2 )
1577  z3t = pi / 6.0 * sum( x(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
1578 
1579  m_mean = z0t / ( pi / 6.0 )
1580 
1581  DO i = 1, ncomp
1582  DO j = 1, ncomp
1583  dij_ab(i,j) = dhs(i)*dhs(j) / ( dhs(i) + dhs(j) )
1584  END DO
1585  END DO
1586 
1587  !-----------------------------------------------------------------------------
1588  ! dispersion term parameters for chain molecules
1589  !-----------------------------------------------------------------------------
1590  DO m = 0, 6
1591  apar(m) = ap(m,1) + (1.0-1.0/m_mean)*ap(m,2) + (1.0-1.0/m_mean)*(1.0-2.0/m_mean)*ap(m,3)
1592  bpar(m) = bp(m,1) + (1.0-1.0/m_mean)*bp(m,2) + (1.0-1.0/m_mean)*(1.0-2.0/m_mean)*bp(m,3)
1593  END DO
1594 
1595 
1596  !-----------------------------------------------------------------------------
1597  ! van der Waals mixing rules for perturbation terms
1598  !-----------------------------------------------------------------------------
1599  order1 = 0.0
1600  order2 = 0.0
1601  DO i = 1, ncomp
1602  DO j = 1, ncomp
1603  order1 = order1 + x(i)*x(j)* mseg(i)*mseg(j)*sig_ij(i,j)**3 * uij(i,j)/t
1604  order2 = order2 + x(i)*x(j)* mseg(i)*mseg(j)*sig_ij(i,j)**3 * (uij(i,j)/t)**2
1605  END DO
1606  END DO
1607 
1608 
1609  !-----------------------------------------------------------------------------
1610  ! association and polar parameters
1611  !-----------------------------------------------------------------------------
1612  assoc = .false.
1613  qudpole = .false.
1614  dipole = .false.
1615  DO i = 1, ncomp
1616  IF (nint(parame(i,12)) /= 0) assoc = .true.
1617  IF (parame(i,7) /= 0.0) qudpole = .true.
1618  IF (parame(i,6) /= 0.0) dipole = .true.
1619  END DO
1620 
1621  ! --- dipole and quadrupole constants ----------------------------------------
1622  IF (qudpole) CALL qq_const ( qqp2, qqp3, qqp4 )
1623  IF (dipole) CALL dd_const ( ddp2, ddp3, ddp4 )
1624  IF (dipole .AND. qudpole) CALL dq_const ( dqp2, dqp3, dqp4 )
1625 
1626 
1627  ! --- TPT-1-association parameters -------------------------------------------
1628  IF (assoc) THEN
1629 
1630  eps_kij = 0.0
1631  k_kij = 0.0
1632 
1633  DO i = 1, ncomp
1634  IF (nint(parame(i,12)) /= 0) THEN
1635  nhb_typ(i) = nint(parame(i,12))
1636  kap_hb(i,i) = parame(i,13)
1637  no = 0
1638  DO j = 1, nhb_typ(i)
1639  DO k = 1, nhb_typ(i)
1640  eps_hb(i,i,j,k) = parame(i,(14+no))
1641  no=no+1
1642  END DO
1643  END DO
1644  DO j = 1, nhb_typ(i)
1645  nhb_no(i,j) = parame(i,(14+no))
1646  no=no+1
1647  END DO
1648  ELSE
1649  nhb_typ(i) = 0
1650  kap_hb(i,i)= 0.0
1651  DO k = 1, nsite
1652  DO l = 1, nsite
1653  eps_hb(i,i,k,l) = 0.0
1654  END DO
1655  END DO
1656  END IF
1657  END DO
1658 
1659  DO i = 1, ncomp
1660  DO j = 1, ncomp
1661  IF (i /= j .AND. (nhb_typ(i) /= 0 .AND. nhb_typ(j) /= 0)) THEN
1662  kap_hb(i,j)= (kap_hb(i,i)*kap_hb(j,j))**0.5 &
1663  *((parame(i,2)*parame(j,2))**3 )**0.5 &
1664  /(0.5*(parame(i,2)+parame(j,2)))**3
1665  kap_hb(i,j)= kap_hb(i,j)*(1.0-k_kij)
1666  DO k = 1, nhb_typ(i)
1667  DO l = 1, nhb_typ(j)
1668  IF (k /= l) THEN
1669  eps_hb(i,j,k,l) = (eps_hb(i,i,k,l)+eps_hb(j,j,l,k))/2.0
1670  eps_hb(i,j,k,l) = eps_hb(i,j,k,l)*(1.0-eps_kij)
1671  END IF
1672  END DO
1673  END DO
1674  END IF
1675  END DO
1676  END DO
1677  IF (nhb_typ(1) == 3) THEN
1678  ! write(*,*)'eps_hb manuell eingegeben'
1679  eps_hb(1,2,3,1) = 0.5*(eps_hb(1,1,3,1)+eps_hb(2,2,1,2))
1680  eps_hb(2,1,1,3) = eps_hb(1,2,3,1)
1681  END IF
1682  IF (nhb_typ(2) == 3) THEN
1683  eps_hb(2,1,3,1) = 0.5*(eps_hb(2,2,3,1)+eps_hb(1,1,1,2))
1684  eps_hb(1,2,1,3) = eps_hb(2,1,3,1)
1685  END IF
1686 
1687  DO i = 1, ncomp
1688  DO k = 1, nhb_typ(i)
1689  DO j = 1, ncomp
1690  DO l = 1, nhb_typ(j)
1691  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)
1692  END DO
1693  END DO
1694  END DO
1695  END DO
1696 
1697  END IF
1698 
1699 END SUBROUTINE perturbation_parameter
1700 
1701 
1702 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1703 ! SUBROUTINE DENSITY_ITERATION
1704 !
1705 ! iterates the density until the calculated pressure 'pges' is equal to
1706 ! the specified pressure 'p'. A Newton-scheme is used for determining
1707 ! the root to the objective function f(eta) = (pges / p ) - 1.0.
1708 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1709 
1710 SUBROUTINE density_iteration
1711 
1712  USE eos_variables
1713  IMPLICIT NONE
1714 
1715  !-----------------------------------------------------------------------------
1716  INTEGER :: i, start, max_i
1717  REAL :: eta_iteration
1718  REAL :: error, dydx, acc_i, delta_eta
1719  !-----------------------------------------------------------------------------
1720 
1721 
1722  IF ( densav(phas) /= 0.0 .AND. eta_start == denold(phas) ) THEN
1723  denold(phas) = eta_start
1724  eta_start = densav(phas)
1725  ELSE
1726  denold(phas) = eta_start
1727  densav(phas) = eta_start
1728  END IF
1729 
1730 
1731  acc_i = 1.e-9
1732  max_i = 30
1733  density_error(:) = 0.0
1734 
1735  i = 0
1736  eta_iteration = eta_start
1737 
1738  !-----------------------------------------------------------------------------
1739  ! iterate density until p_calc = p
1740  !-----------------------------------------------------------------------------
1741  iterate_density: DO
1742 
1743  i = i + 1
1744  eta = eta_iteration
1745 
1746  CALL p_eos_interface
1747 
1748  error = (pges / p ) - 1.0
1749 
1750  !--------------------------------------------------------------------------
1751  ! correction for instable region
1752  !--------------------------------------------------------------------------
1753  IF ( pgesdz < 0.0 .AND. i < max_i ) THEN
1754  IF ( error > 0.0 .AND. pgesd2 > 0.0 ) THEN ! no liquid density
1755  CALL pressure_spinodal
1756  eta_iteration = eta
1757  error = (pges / p ) - 1.0
1758  IF ( ((pges/p ) -1.0) > 0.0 ) eta_iteration = 0.001 ! no solution possible
1759  IF ( ((pges/p ) -1.0) <=0.0 ) eta_iteration = eta_iteration * 1.1 ! no solution found so far
1760  ELSE IF ( error < 0.0 .AND. pgesd2 < 0.0 ) THEN ! no vapor density
1761  CALL pressure_spinodal
1762  eta_iteration = eta
1763  error = (pges / p ) - 1.0
1764  IF ( ((pges/p ) -1.0) < 0.0 ) eta_iteration = 0.5 ! no solution possible
1765  IF ( ((pges/p ) -1.0) >=0.0 ) eta_iteration = eta_iteration * 0.9 ! no solution found so far
1766  ELSE
1767  eta_iteration = (eta_iteration + eta_start) / 2.0
1768  IF (eta_iteration == eta_start) eta_iteration = eta_iteration + 0.2
1769  END IF
1770  cycle iterate_density
1771  END IF
1772 
1773 
1774  dydx = pgesdz/p
1775  delta_eta = error/ dydx
1776  IF ( delta_eta > 0.05 ) delta_eta = 0.05
1777  IF ( delta_eta < -0.05 ) delta_eta = -0.05
1778 
1779  eta_iteration = eta_iteration - delta_eta
1780 
1781  IF (eta_iteration > 0.9) eta_iteration = 0.6
1782  IF (eta_iteration <= 0.0) eta_iteration = 1.e-16
1783  start = 1
1784 
1785  IF ( abs(error*p/pgesdz) < 1.e-12 ) start = 0
1786  IF ( abs(error) < acc_i ) start = 0
1787  IF ( i > max_i ) THEN
1788  start = 0
1789  density_error(phas) = abs( error )
1790  ! write (*,*) 'density iteration failed'
1791  END IF
1792 
1793  IF (start /= 1) EXIT iterate_density
1794 
1795  END DO iterate_density
1796 
1797  eta = eta_iteration
1798 
1799  IF ((eta > 0.3 .AND. densav(phas) > 0.3) .OR. &
1800  (eta < 0.1 .AND. densav(phas) < 0.1)) densav(phas) = eta
1801 
1802 END SUBROUTINE density_iteration
1803 
1804 
1805 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1806 ! SUBROUTINE PRESSURE_SPINODAL
1807 !
1808 ! iterates the density until the derivative of pressure 'pges' to
1809 ! density is equal to zero. A Newton-scheme is used for determining
1810 ! the root to the objective function.
1811 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1812 
1813 SUBROUTINE pressure_spinodal
1814 
1815  USE eos_variables
1816  IMPLICIT NONE
1817 
1818  !-----------------------------------------------------------------------------
1819  INTEGER :: i, max_i
1820  REAL :: eta_iteration
1821  REAL :: error, acc_i, delta_eta
1822  !-----------------------------------------------------------------------------
1823 
1824  acc_i = 1.e-6
1825  max_i = 30
1826 
1827  i = 0
1828  eta_iteration = eta_start
1829 
1830  !-----------------------------------------------------------------------------
1831  ! iterate density until p_calc = p
1832  !-----------------------------------------------------------------------------
1833 
1834  error = acc_i + 1.0
1835  DO WHILE ( abs(error) > acc_i .AND. i < max_i )
1836 
1837  i = i + 1
1838  eta = eta_iteration
1839 
1840  CALL p_eos_interface
1841 
1842  error = pgesdz
1843 
1844  delta_eta = error/ pgesd2
1845  IF ( delta_eta > 0.02 ) delta_eta = 0.02
1846  IF ( delta_eta < -0.02 ) delta_eta = -0.02
1847 
1848  eta_iteration = eta_iteration - delta_eta
1849  ! write (*,'(a,i3,3G18.10)') 'iter',i, error, eta_iteration, pgesdz
1850 
1851  IF (eta_iteration > 0.9) eta_iteration = 0.5
1852  IF (eta_iteration <= 0.0) eta_iteration = 1.e-16
1853 
1854  END DO
1855 
1856  eta = eta_iteration
1857 
1858 END SUBROUTINE pressure_spinodal
1859 
1860 
1861 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1862 !
1863 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1864 
1865 SUBROUTINE p_dz
1866 
1867  USE parameters, ONLY: nc, nsite
1868  USE eos_variables
1869  IMPLICIT NONE
1870 
1871  !-----------------------------------------------------------------------------
1872  REAL :: eta_0, dist, fact
1873  REAL :: fres1, fres2, fres3, fres4, fres5
1874  REAL :: df_dr, df_dr2, df_dr3, df_dr4
1875  !-----------------------------------------------------------------------------
1876 
1877 
1878  IF (eta > 1.e-1) THEN
1879  fact = 1.0
1880  ELSE IF (eta <= 1.e-1 .AND. eta > 1.e-2) THEN
1881  fact = 10.0
1882  ELSE
1883  fact = 100.0
1884  END IF
1885  dist = eta * 5.e-4 * fact
1886  ! dist = eta * 4.E-3 * fact
1887 
1888 
1889  eta_0 = eta
1890  eta = eta_0 - 2.0*dist
1891  CALL p_eos
1892  fres1 = pges
1893  eta = eta_0 - dist
1894  CALL p_eos
1895  fres2 = pges
1896  eta = eta_0 + dist
1897  CALL p_eos
1898  fres3 = pges
1899  eta = eta_0 + 2.0*dist
1900  CALL p_eos
1901  fres4 = pges
1902  eta = eta_0
1903  CALL p_eos
1904  fres5 = pges
1905 
1906  df_dr = (-fres4+8.0*fres3-8.0*fres2+fres1)/(12.0*dist)
1907  df_dr2 = (-fres4+16.0*fres3-3.d1*fres5+16.0*fres2-fres1) &
1908  /(12.0*(dist**2 ))
1909  df_dr3 = (fres4-2.0*fres3+2.0*fres2-fres1) /(2.0*dist**3 )
1910  df_dr4 = (fres4-4.0*fres3+6.0*fres5-4.0*fres2+fres1) /(1.0*dist**4 )
1911 
1912  WRITE (*,*) 'f` = ',df_dr
1913  WRITE (*,*) 'f`` = ',df_dr2
1914  WRITE (*,*) 'f``` = ',df_dr3
1915  WRITE (*,*) 'f````= ',df_dr4,eta
1916  ! if (eta.gt.0.3) stop
1917 
1918 END SUBROUTINE p_dz
1919 
1920 
1921 
1922 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1923 !
1924 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1925 
1926 SUBROUTINE f_dz
1927 
1928  USE parameters, ONLY: nc, nsite
1929  USE eos_variables
1930  IMPLICIT NONE
1931 
1932  !-----------------------------------------------------------------------------
1933  REAL :: eta_0, dist, fact
1934  REAL :: fres1, fres2, fres3, fres4, fres5
1935  REAL :: df_dr, df_dr2, df_dr3, df_dr4
1936  !-----------------------------------------------------------------------------
1937 
1938 
1939  IF (eta > 1.e-1) THEN
1940  fact = 1.0
1941  ELSE IF (eta <= 1.e-1 .AND. eta > 1.e-2) THEN
1942  fact = 10.0
1943  ELSE
1944  fact = 100.0
1945  END IF
1946  dist = eta*5.e-4 *fact
1947  ! dist = eta*4.E-3 *fact
1948 
1949  eta_0 = eta
1950  eta = eta_0 - 2.0*dist
1951  CALL f_eos
1952  fres1 = tfr
1953  eta = eta_0 - dist
1954  CALL f_eos
1955  fres2 = tfr
1956  eta = eta_0 + dist
1957  CALL f_eos
1958  fres3 = tfr
1959  eta = eta_0 + 2.0*dist
1960  CALL f_eos
1961  fres4 = tfr
1962  eta = eta_0
1963  CALL f_eos
1964  fres5 = tfr
1965 
1966  df_dr = (-fres4+8.0*fres3-8.0*fres2+fres1)/(12.0*dist)
1967  df_dr2 = (-fres4+16.0*fres3-3.d1*fres5+16.0*fres2-fres1) &
1968  /(12.0*(dist**2 ))
1969  df_dr3 = (fres4-2.0*fres3+2.0*fres2-fres1) /(2.0*dist**3 )
1970  df_dr4 = (fres4-4.0*fres3+6.0*fres5-4.0*fres2+fres1) /(1.0*dist**4 )
1971 
1972  WRITE (*,*) 'f` = ',df_dr
1973  WRITE (*,*) 'f`` = ',df_dr2
1974  WRITE (*,*) 'f``` = ',df_dr3
1975  WRITE (*,*) 'f````= ',df_dr4,eta
1976  WRITE (*,*) 'z = ',df_dr*eta
1977  WRITE (*,*) 'z` = ',df_dr2*eta + df_dr
1978  WRITE (*,*) 'z`` = ',df_dr3*eta + 2.0* df_dr2
1979  WRITE (*,*) 'z```= ',df_dr4*eta + 3.0* df_dr3 ,eta
1980 
1981 END SUBROUTINE f_dz
1982 
1983 
1984 
1985 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1986 !
1987 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1988 
1989 SUBROUTINE dfr_eos
1990 
1991  USE parameters, ONLY: nc, nsite
1992  USE eos_variables
1993  IMPLICIT NONE
1994 
1995  !-----------------------------------------------------------------------------
1996  INTEGER :: i, k
1997  REAL :: flsum
1998  REAL, DIMENSION(nc) :: grada, atilde2, atilde, xsav
1999  !-----------------------------------------------------------------------------
2000 
2001  rho = eta/z3t
2002 
2003  CALL f_eos
2004 
2005  flsum = 0.0
2006  DO k = 1, ncomp
2007  flsum = flsum + log(rho*x(k))
2008  END DO
2009  ! flsum = flsum -LOG(p*1.E-30/(rho*kbol*t))
2010  DO k = 1, ncomp
2011  atilde(k) = fres*rho+x(k)*rho*log(rho*x(k))
2012  END DO
2013 
2014 
2015 
2016  DO i=1,ncomp
2017  xsav(i) = x(i)
2018  x(i) = xsav(i) + 1.e-8
2019  CALL f_eos
2020  flsum = 0.0
2021  DO k = 1, ncomp
2022  flsum=flsum+ rho*x(k)*log(x(k))
2023  END DO
2024  ! flsum = flsum -LOG(p*1.E-30/(rho*kbol*t))
2025  DO k = 1, ncomp
2026  atilde2(k) = fres*rho+x(k)*rho*log(rho*x(k))
2027  END DO
2028  x(i) = xsav(i)
2029 
2030  grada(i) = (atilde2(i)-atilde(i) ) /1.e-8
2031  WRITE (*,*) 'grad',i,grada(i)
2032 
2033  END DO
2034 
2035 END SUBROUTINE dfr_eos
determinant
double determinant(double **in_matrix, int n)
Definition: determinant.h:12
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
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29