MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
eos_critical_renorm.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 !
3 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
4 
5 SUBROUTINE f_eos_rn
6  USE basic_variables, ONLY: ncomp, rgt_variant
7  IMPLICIT NONE
8 
9  IF (rgt_variant == 'phase_cell' .AND. ncomp >= 2) THEN
10  CALL f_eos_rn_phase_space_cell
11  ELSE
12  CALL f_eos_rn_isomorphic_density
13  END IF
14 
15 END SUBROUTINE f_eos_rn
16 
17 
18 
19 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
20 !
21 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
22 
23 SUBROUTINE phi_critical_renorm
24  USE basic_variables, ONLY: ncomp, rgt_variant
25  IMPLICIT NONE
26 
27  IF (rgt_variant == 'phase_cell' .AND. ncomp >= 2) THEN
28  CALL phi_critical_renorm_phase_space_cell
29  ELSE
30  CALL phi_critical_renorm_isomorphic_density
31  END IF
32 
33 END SUBROUTINE phi_critical_renorm
34 
35 
36 
37 
38 
39 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
40 ! SUBROUTINE F_EOS_RN_ISOMORPHIC_DENSITY
41 !
42 ! White's recursive procedure is:
43 !
44 ! f = fid + f0 + delta_f1 + delta_f2 + delta_f3 + ...
45 !
46 ! where f is the Helmholtz energy per unit volume, and where f0 is the
47 ! classical residual repulsive part, f0 = f_res,rep.
48 !
49 ! PC-SAFT can be considered for the classical residual part, with
50 !
51 ! f0 = f_PCSAFT
52 !
53 ! It is
54 !
55 ! INT_0^rho( exp(-Kn^-1*(fns(rho+x)-2fns(rho)+fns(rho-x))) )
56 ! delta_fn= -Kn*ln----------------------------------------------------------
57 ! INT_0^rho( exp(-Kn^-1*(fnl(rho+x)-2fnl(rho)+fnl(rho-x))) )
58 !
59 ! where x in the integrals INT runs from x=0 to x=rho_x=rho. The counter
60 ! kk in belows code takes care of the integrals INT from 0 to rho.
61 ! An outer loop is established (with index k) in order to calculate f for
62 ! the whole density range. First, the f0 is calculated and approximated
63 ! by a cubic spline. Then delta_fn is calculated for 0<=rho<=rhomax. New
64 ! spline parameters are determined after evey iteration n.
65 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
66 
67 SUBROUTINE f_eos_rn_isomorphic_density
68 
69  USE parameters
70  USE eos_variables
71  USE eos_constants
72  USE eos_numerical, only: f_numerical
73  USE fitting_rgt_parameters
74  USE dft_module
75  IMPLICIT NONE
76 
77  !-----------------------------------------------------------------------------
78  INTEGER :: i, j, k, m, n, kk, k05max
79  INTEGER :: stepno, niter, pos, step2
80  REAL :: rho0, rhomax, kn(10), ll
81  REAL :: drho, densav1
82  REAL :: alpha_l, alpha_s
83  REAL :: integrand_l, integrand_s
84  REAL :: integr_old_l, integr_old_s
85  REAL :: integral_l, integral_s
86  REAL :: gn_l, gn_s
87  REAL :: i1, w_ratio, del_f(8,0:10000)
88  REAL :: phi_crit
89  REAL :: alph(0:8,0:10000)
90  REAL :: rhovec(0:10000), xdev, dzp1, fres_v
91  REAL :: sig_m2
92  REAL :: m_mean
93  REAL :: chapm, order_chap
94  REAL :: fid, rho_i(nc)
95  REAL, DIMENSION(5000) :: utri
96  REAL, DIMENSION(0:5000) :: alphn, beta, gamma, delta
97  REAL :: rho_c, rho_l, rho_r
98  REAL :: alp_c, alph_lc, alph_sc
99  REAL :: alp_r, alph_lr, alph_sr
100  REAL :: alp_l, alph_ll, alph_sl
101  REAL :: dzp2, delta_rho
102 
103  INTEGER, SAVE :: scan = 0
104  REAL, SAVE :: tempsav = 0.0
105  REAL, DIMENSION(5), SAVE :: xsav = 0.0
106  REAL, DIMENSION(0:10000), SAVE :: alphnsav, betasav, gammasav, deltasav, rhovecsav
107  !SAVE alphnsav, betasav, gammasav, deltasav, rhovecsav
108 
109  !SAVE tempsav,xsav,scan
110  !DATA scan /0/
111  !DATA tempsav /0.0/
112  !DATA xsav /0.0,0.0,0.0,0.0,0.0/
113  !-----------------------------------------------------------------------------
114 
115  IF ( critfit == 1 ) THEN
116  IF (t < 0.2*parame(1,3) ) t = 0.2*parame(1,3)
117  lli(1) = llfit*sig_ij(1,1)
118  phi_criti(1) = phifit
119  chap(1) = chapfit
120  END IF
121 
122  CALL perturbation_parameter
123 
124  ll = 0.0
125  phi_crit = 0.0
126  sig_m2 = 0.0
127  m_mean = 0.0
128  chapm = 0.0
129  order_chap = 0.0
130  DO i = 1,ncomp
131  IF (lli(i) == 0.0 .OR. phi_criti(i) == 0.0 .OR. chap(i) == 0.0) stop
132  ll = ll + x(i)*lli(i)**3
133  phi_crit = phi_crit + x(i)*phi_criti(i)
134  sig_m2 = sig_m2 + x(i)*mseg(i)*sig_ij(i,i)**2
135  m_mean = m_mean + x(i)*mseg(i)
136  chapm = chapm + x(i)*chap(i)
137  DO j = 1,ncomp
138  order_chap = order_chap + x(i)*x(j)* mseg(i)*mseg(j) &
139  *sig_ij(i,j)**3 * uij(i,j)/t *(chap(i)*chap(j))**0.5
140  END DO
141  END DO
142  ll = ll**(1.0/3.0)
143  sig_m2 = sig_m2/m_mean
144 
145  stepno = 400
146  niter = 5
147  w_ratio = 2.0
148 
149  rho0 = eta/z3t
150  rhomax = 0.0
151  DO i = 1,ncomp
152  ! rhomax = rhomax + x(i)*mseg(i)*sig_ij(i,i)**3
153  rhomax = rhomax + x(i)*mseg(i)*dhs(i)**3 !jg
154  END DO
155  rhomax = sqrt(2.0)/rhomax
156 
157 
158  densav1 = eta
159 
160  xdev = sum( abs( xsav(1:ncomp)-x(1:ncomp) ) )
161  IF (tempsav == t .AND. xdev <= 1.e-10 .AND. scan == 1) GO TO 5
162 
163  tempsav = t
164  xsav(1:ncomp) = x(1:ncomp)
165  !-----------------------------------------------------------------------------
166 
167 
168 
169  dzp1 = rhomax / REAL(stepno)
170  alph(0,0) = 0.0
171  rhovec(0) = 0.0
172  DO k = 1, stepno
173  drho = REAL(k)*dzp1
174  rhovec(k) = drho
175  alph(0,k) = 0.0
176  eta = drho*z3t
177  fid = 0.0
178  rho_i = 0.0
179  DO i = 1,ncomp
180  ! debroglie(i) = 6.62606896E-34 *1.E10 & ! in units Angstrom
181  ! *SQRT( 1.0 / (2.0*PI *1.0 /6.022045E23/1000.0*KBOL*T) )
182  rho_i(i) = x(i)*drho
183  IF (rho_i(i) > 0.0) fid = fid + x(i)*(log(rho_i(i))-1.0)
184  END DO
185  CALL f_numerical
186  alph(0,k) = (fres+fid)*drho
187  alphn(k) = alph(0,k)
188  END DO
189  ! CALL SPLINE_PARA (dzp1,alphn,utri,stepno)
190  ! CALL SPLINE_COEFF(beta,gamma,delta,dzp1,alphn,utri,stepno)
191  scan = 1
192 
193 
194  ! alpha_l = 2.0/3.0*PI *1.5**3 *order1
195  ! alpha_s = alpha_l*phi_crit* 0.2 *1.5**2 *sig_ij(1,1)**2/LL/LL
196 
197 
198  DO n = 1,niter
199 
200  k05max = stepno/2
201  kn(n) = 1.0 / ll**3 / w_ratio**REAL(3*n)
202  alpha_l = 16.0/9.0 * pi*order_chap !* chapm
203  alpha_s = alpha_l*phi_crit* 9.0/7.0 * sig_m2 /ll/ll /w_ratio**REAL(2*n)/2.0
204 
205  ! DO k = 1,stepno-1 !!!!!
206  DO k = 1,stepno*3/4-1 !!!!!
207  integr_old_l = 1.0
208  integr_old_s = 1.0
209  ! hx(1) = 0.0
210  ! hyl(1) = integr_old_l
211  ! hys(1) = integr_old_s
212  ! nn = 1
213  integral_l = 0.0
214  integral_s = 0.0
215  step2 = k !step2: step# for integrat.(index kk)
216  IF (k > k05max) step2 = stepno-k
217  ! IF (k.LE.20)step2 = 2*step2
218  ! step2 = 2*step2
219  ! IF (step2.GT.stepno) step2 = stepno
220  dzp2 = rhovec(k)/REAL(step2)
221  IF (rhovec(k) > rhomax/2.0) dzp2 = (rhomax-rhovec(k)) / REAL(step2)
222  DO kk = 1, step2
223  delta_rho = REAL(kk)*dzp2
224  rho_c = rhovec(k)
225  rho_l = rhovec(k) - delta_rho
226  rho_r = rhovec(k) + delta_rho
227 
228 
229  eta = rho_c*pi/6.0*mseg(1)*dhs(1)**3
230  i1 = 0.0
231  DO m = 0,6
232  i1 = i1 + apar(m)*eta**REAL(m)
233  END DO
234  ! alpha_l = + 2.0*PI*I1*order_chap
235  ! alpha_s = alpha_l*phi_crit* 9.0/7.0*sig_m2 /LL/LL / w_ratio**REAL(2*n)/2.0
236  ! pos = INT(rho_c/rhomax*stepno + 1.E-9)
237  ! alp_c = alphn(pos) + beta(pos) *(rho_c-rhovec(pos))
238  ! + gamma(pos)*(rho_c-rhovec(pos))**2
239  ! + delta(pos)*(rho_c-rhovec(pos))**3
240  alp_c = alph(n-1,k)
241  alph_lc = alp_c +alpha_l*rho_c*rho_c
242  alph_sc = alp_c +alpha_s*rho_c*rho_c
243 
244  eta = rho_l*pi/6.0*mseg(1)*dhs(1)**3
245  i1 = 0.0
246  DO m = 0,6
247  i1 = i1 + apar(m)*eta**REAL(m)
248  END DO
249  ! alpha_l = + 2.0*PI*I1*order_chap
250  ! alpha_s = alpha_l*phi_crit* 9.0/7.0*sig_m2 /LL/LL / w_ratio**REAL(2*n)/2.0
251  ! pos = INT(rho_l/rhomax*stepno + 1.E-9)
252  ! alp_l = alphn(pos) + beta(pos) *(rho_l-rhovec(pos))
253  ! + gamma(pos)*(rho_l-rhovec(pos))**2
254  ! + delta(pos)*(rho_l-rhovec(pos))**3
255  alp_l = alph(n-1,k-kk)
256  alph_ll = alp_l +alpha_l*rho_l*rho_l
257  alph_sl = alp_l +alpha_s*rho_l*rho_l
258 
259  eta = rho_r*pi/6.0*mseg(1)*dhs(1)**3
260  i1 = 0.0
261  DO m = 0,6
262  i1 = i1 + apar(m)*eta**REAL(m)
263  END DO
264  ! alpha_l = + 2.0*PI*I1*order_chap
265  ! alpha_s = alpha_l*phi_crit* 9.0/7.0*sig_m2 /LL/LL / w_ratio**REAL(2*n)/2.0
266  ! pos = INT(rho_r/rhomax*stepno + 1.E-9)
267  ! alp_r = alphn(pos) + beta(pos) *(rho_r-rhovec(pos))
268  ! + gamma(pos)*(rho_r-rhovec(pos))**2
269  ! + delta(pos)*(rho_r-rhovec(pos))**3
270  alp_r = alph(n-1,k+kk)
271  alph_lr = alp_r +alpha_l*rho_r*rho_r
272  alph_sr = alp_r +alpha_s*rho_r*rho_r
273 
274  gn_l = (alph_lr+alph_ll)/2.0 - alph_lc
275  gn_s = (alph_sr+alph_sl)/2.0 - alph_sc
276  IF (gn_l < 0.0) gn_l = 0.0
277  IF (gn_s < 0.0) gn_s = 0.0
278  ! write (*,*) k,kk,Gn_s
279  integrand_l = 0.0
280  integrand_s = 0.0
281  IF ( -gn_l/kn(n) > -300.0 .AND. -gn_l/kn(n) < 300.0) integrand_l = exp( -gn_l/kn(n) )
282  IF ( -gn_s/kn(n) > -300.0 .AND. -gn_s/kn(n) < 300.0) integrand_s = exp( -gn_s/kn(n) )
283  ! nn = nn+1
284  ! hx(nn) = delta_rho
285  ! hyl(nn) = integrand_l
286  ! hys(nn) = integrand_s
287  integral_l = integral_l + dzp2 *(integrand_l+integr_old_l)/2.0
288  integral_s = integral_s + dzp2 *(integrand_s+integr_old_s)/2.0
289  integr_old_l = integrand_l
290  integr_old_s = integrand_s
291 
292  IF (integrand_l < 1.e-9 .AND. integrand_s < 1.e-9) THEN
293  GO TO 15 ! end loop, because no significant contribution is expected
294  END IF
295 
296  END DO ! enddo kk (integration over density)
297  ! CALL spline(hx,hyl,nn,0.0,0.0,y2)
298  ! CALL splint_integral(hx,hyl,y2,nn,0.0,delta_rho,integral_l)
299  ! CALL spline(hx,hys,nn,0.0,0.0,y2)
300  ! CALL splint_integral(hx,hys,y2,nn,0.0,delta_rho,integral_s)
301 15 CONTINUE
302 
303  del_f(n,k) = 0.0
304  IF (integral_s /= 0.0 .AND. integral_l /= 0.0) THEN
305  del_f(n,k) = -kn(n)*log(integral_s/integral_l)
306  END IF
307  END DO ! enddo k
308 
309  del_f(n,0) = 0.0
310  DO k = 0,stepno
311  alph(n,k) = alph(n-1,k) + del_f(n,k)
312  alphn(k) = alph(n,k)
313  ! subract the ideal gas part in order to have less non-linearity in the low
314  ! density region. The ideal gas part is then added to the final value below.
315  IF(n == niter .AND. k /= 0) THEN
316  fid = 0.0
317  DO i = 1,ncomp
318  rho_i(i) = x(i)*rhovec(k)
319  IF(rho_i(i) > 0.0) fid = fid +rho_i(i)*(log(rho_i(i))-1.0)
320  END DO
321  alphn(k) = alphn(k)-fid
322  END IF
323  END DO
324  CALL spline_para (dzp1,alphn,utri,stepno)
325  CALL spline_coeff(beta,gamma,delta,dzp1,alphn,utri,stepno)
326 
327  DO k = 0,stepno
328  alphnsav(k) = alphn(k)
329  betasav(k) = beta(k)
330  gammasav(k) = gamma(k)
331  deltasav(k) = delta(k)
332  rhovecsav(k) = rhovec(k)
333  END DO
334 
335  END DO ! loop of n cycles
336  ! write (71,*) ' '
337 
338 
339 5 CONTINUE
340  ! CALL SPLINE_PARA (dzp1,alphn,utri,stepno)
341  ! CALL SPLINE_COEFF(beta,gamma,delta,dzp1,alphn,utri,stepno)
342 
343  DO k = 0,stepno
344  alphn(k) = alphnsav(k)
345  beta(k) = betasav(k)
346  gamma(k) = gammasav(k)
347  delta(k) = deltasav(k)
348  rhovec(k) = rhovecsav(k)
349  END DO
350  ! write (*,*) 'beta',beta(70),gamma(70),t
351 
352 
353 
354  ! fid = 0.0
355  ! DO i = 1,ncomp
356  ! ! in units Angstrom
357  ! debroglie(i) = 6.62606896d-34 *1d10 *SQRT( 1.0 / (2.0*PI *1.0 /6.022045d23/1000.0*KBOL*T) )
358  ! rho_i(i) = x(i)*rho0
359  ! fid = fid + x(i)*(LOG(rho_i(i))-1.0)
360  ! ENDDO
361  ! !fid = 0.0
362 
363 
364  pos = int(rho0/rhomax*REAL(stepno))
365  fres_v = alphn(pos) + beta(pos) *(rho0-rhovec(pos)) &
366  + gamma(pos)*(rho0-rhovec(pos))**2 + delta(pos)*(rho0-rhovec(pos))**3
367  fres = fres_v/rho0 !-fid
368 
369  pges = ( (beta(pos)+ 2.0*gamma(pos)*(rho0-rhovec(pos)) &
370  + 3.0* delta(pos)*(rho0-rhovec(pos))**2 )*rho0 -fres_v ) &
371  *(kbol*t)/1.e-30 + rho0 * (kbol*t) / 1.e-30 ! ideal gas contribution
372  ! if(pos.GT.1)write (69,*) pos,delta(pos)
373  ! & (beta(pos)+ 2.0*gamma(pos)*(rho0-rhovec(pos))
374  ! & + 3.0* delta(pos)*(rho0-rhovec(pos))**2 )
375  pgesdz = ( 2.0*gamma(pos)*rho0 &
376  + rho0*6.0* delta(pos)*(rho0-rhovec(pos)) ) *(kbol*t)/1.e-30 /z3t &
377  + 1.0/z3t*(kbol*t)/1.e-30 ! ideal gas contribution
378  pgesd2 = ( 2.0*gamma(pos) + 6.0* delta(pos)*(2.0*rho0-rhovec(pos)) ) *(kbol*t)/1.e-30 /z3t /z3t
379  pgesd3 = ( 12.0* delta(pos) ) *(kbol*t)/1.e-30 /z3t /z3t /z3t
380 
381 
382  eta = densav1
383 
384 END SUBROUTINE f_eos_rn_isomorphic_density
385 
386 
387 
388 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
389 !
390 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
391 
392 SUBROUTINE phi_critical_renorm_isomorphic_density
393 
394  USE parameters, ONLY: nc, nsite, pi, kbol
395  USE eos_constants
396  USE eos_variables
397  IMPLICIT NONE
398 
399  !-----------------------------------------------------------------------------
400  INTEGER :: i, k
401  REAL :: myresq
402  REAL :: zres, zges, sumseg
403  REAL :: fres1, fres2 = 0.0, fres3, fres4
404  REAL :: xconc, d_org, d_new
405  REAL :: dicht, gradi
406  REAL :: myres(nc), mpart(nc)
407  REAL :: gres !,my_pt1, p_pt1
408  REAL :: tfr_1, tfr_2, tfr_3, tfr_4
409  !-----------------------------------------------------------------------------
410 
411  !-----------------------------------------------------------------------------
412  ! density iteration or pressure calculation
413  !-----------------------------------------------------------------------------
414  IF (ensemble_flag == 'tp') THEN
415  CALL density_iteration
416  ELSEIF (ensemble_flag == 'tv') THEN
417  eta = eta_start
418  rho = eta/z3t
419  sumseg = z0t/(pi/6.0) ! mean segment number
420  CALL f_eos_rn
421  ELSE
422  write (*,*) 'PHI_CRITICAL_RENORM: define ensemble, ensemble_flag == (pv) or (pt)'
423  stop
424  END IF
425 
426 !!$!---------------density iteration-------------------------------------
427 !!$IF (dfti >= 98) THEN
428 !!$ eta = eta_start
429 !!$ rho = eta/z3t
430 !!$ sumseg = z0t/(PI/6.0) ! mean segment number
431 !!$ CALL F_EOS_RN
432 !!$ELSE
433 !!$ CALL densitr_rn
434 !!$END IF
435 !!$
436  !-----------------------------------------------------------------------------
437  ! compressibility factor z = p/(kT*rho)
438  !-----------------------------------------------------------------------------
439  zges = (p * 1.e-30) / (kbol*t*eta/z3t)
440  IF ( ensemble_flag == 'tv' ) zges = (pges * 1.e-30) / (kbol*t*eta/z3t)
441  zres = zges - 1.0
442  !z_ges = zges
443 
444 
445  IF (ncomp == 1) THEN
446  gres = fres + zres
447  IF (ensemble_flag == 'tp') lnphi(1) = gres - log(zges)
448  IF (ensemble_flag == 'tv' .AND. eta >= 0.0) lnphi(1) = gres !+LOG(rho)
449  END IF
450 
451 !!$!--------------------------
452 !!$IF (dfti == 99) THEN
453 !!$ ! z3 = rho*PI/6.0*mseg(1)*dhs(1)**3
454 !!$ CALL phipt1 ( my_pt1, p_pt1 )
455 !!$ lnphi(1) = lnphi(1)+my_pt1
456 !!$ z_ges = z_ges+(p_pt1 * 1.E-30)/(KBOL*t*rho)
457 !!$END IF
458 !!$!--------------------------
459 
460 
461  IF (ncomp /= 1) THEN
462  !----------------------derivations with big D--------------------------------
463  ! gradi = 1.E-14
464  gradi = 1.e-13
465 
466  DO k = 1,ncomp
467  xconc = x(k)
468  dicht = eta
469  d_org = z3t/pi*6.0
470 
471  IF (xconc > 1.e-13 .AND. xconc < 0.9999) THEN
472  x(k) = xconc-gradi*10.0**(0.5*(15.0+log10(xconc)))
473  ELSE
474  IF (xconc < 0.9999) THEN
475  x(k) = xconc + 2.e-11
476  ELSE
477  x(k) = xconc - 2.e-10
478  END IF
479  END IF
480  d_new = 0.0
481  DO i = 1,ncomp
482  d_new = d_new +x(i)*mseg(i)*dhs(i)**3
483  END DO
484  eta = dicht*d_new/d_org
485  CALL f_eos_rn
486  fres1 = fres
487  tfr_1 = tfr
488 
489  IF (xconc > 1.e-13 .AND. xconc < 0.9999) THEN
490  !****************************************
491  x(k) = xconc+gradi*10.0**(0.5*(15.0+log10(xconc)))
492  d_new = 0.0
493  DO i = 1,ncomp
494  d_new = d_new +x(i)*mseg(i)*dhs(i)**3
495  END DO
496  eta = dicht*d_new/d_org
497  CALL f_eos_rn
498  fres2 = fres
499  tfr_2 = tfr
500  !****************************************
501  END IF
502 
503  IF (xconc > 1.e-13 .AND. xconc < 0.9999) THEN
504  x(k) = xconc-0.5*gradi*10.0**(0.5*(15.0+log10(xconc)))
505  ELSE
506  x(k) = xconc
507  END IF
508  d_new = 0.0
509  DO i = 1,ncomp
510  d_new = d_new +x(i)*mseg(i)*dhs(i)**3
511  END DO
512  eta = dicht*d_new/d_org
513  CALL f_eos_rn
514  fres3 = fres
515  tfr_3 = tfr
516 
517  IF (xconc > 1.e-13 .AND. xconc < 0.9999) THEN
518  x(k) = xconc+0.5*gradi*10.0**(0.5*(15.0+log10(xconc)))
519  ELSE
520  IF (xconc < 0.9999) THEN
521  x(k) = xconc + 1.e-11
522  ELSE
523  x(k) = xconc - 1.e-10
524  END IF
525  END IF
526  d_new = 0.0
527  DO i = 1,ncomp
528  d_new = d_new +x(i)*mseg(i)*dhs(i)**3
529  END DO
530  eta = dicht*d_new/d_org
531  CALL f_eos_rn
532  fres4 = fres
533  tfr_4 = tfr
534 
535  x(k) = xconc
536  eta = dicht
537 
538 
539  IF (xconc > 1.e-13 .AND. xconc < 0.9999) THEN
540  mpart(k) = (fres1-8.0*fres3+8.0*fres4-fres2) &
541  /(6.0*gradi*10.0**(0.5*(15.0+log10(xconc))))
542  ! tfrdx(k) = (tfr_1-8.0*tfr_3+8.0*tfr_4-tfr_2)
543  ! & /(6.0*gradi*10.0**(0.5d0*(15.0+LOG10(xconc))))
544  ! write (*,*) tfrdx,eta
545  ELSE
546  IF (xconc < 0.9999) THEN
547  mpart(k) = (-3.0*fres3+4.0*fres4-1.0*fres1)/2.e-11
548  ! MPART(k) = (fres4-fres3)/1.E-12
549  ELSE
550  mpart(k) = - (-4.0*fres3+6.0*fres4-2.0*fres1)/2.e-10
551  ! MPART(k) = (fres3-fres4)/1.E-5
552  END IF
553  END IF
554 
555  END DO
556 
557  !--------------------------------------------------------------------------
558  ! f res
559  !--------------------------------------------------------------------------
560  CALL f_eos_rn
561 
562 
563  myresq = 0.0
564  DO i = 1, ncomp
565  myresq = myresq - x(i)* mpart(i)
566  END DO
567 
568  DO k = 1, ncomp
569  myres(k) = myresq + mpart(k) + fres + zres
570  IF (ensemble_flag == 'tp') lnphi(k) = myres(k) - log(zges)
571  IF (ensemble_flag == 'tv' .AND. eta >= 0.0) lnphi(1) = myres(k) !+LOG(rho)
572  ! IF (DFT.GE.98) write (*,*) dft
573  ! write (*,*) 'lnphi',k,LNPHI(k),x(k),MYRES(k), -LOG(ZGES)
574  END DO
575  ! write(*,*) t,p,eta,x(1)
576  ! stop
577 
578  END IF ! ncomp.NE.1
579 
580 END SUBROUTINE phi_critical_renorm_isomorphic_density
581 
582 
583 
584 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
585 !
586 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
587 
588 SUBROUTINE p_eos_rn
589 
590  USE eos_variables
591  IMPLICIT NONE
592 
593  !-----------------------------------------------------------------------------
594  REAL :: z0,z1,z2,z3
595  REAL :: dzetdv,dicht,dist,fact
596  REAL :: pgesdz1,pgesdz4,pgesdz5
597  REAL :: pges2d1,pges2d4,pges2d5
598  !-----------------------------------------------------------------------
599 
600 
601  rho = eta/z3t
602  z0 = z0t*rho
603  z1 = z1t*rho
604  z2 = z2t*rho
605  z3 = z3t*rho
606 
607  !-----------------density-------------------------------------------
608  dzetdv = eta*rho
609 
610  IF (eta > 0.1) THEN
611  fact = 1.0
612  ELSE IF (eta <= 0.1 .AND. eta > 0.01) THEN
613  fact = 10.0
614  ELSE
615  fact = 50.0
616  END IF
617  dist = eta*1.e-2 *fact
618 
619  dicht = eta
620  eta = dicht - 2.0*dist
621  CALL f_eos_rn
622  pgesdz1 = pgesdz
623  pges2d1 = pgesd2
624  ! eta = dicht - dist
625  ! CALL F_EOS_RN
626  ! pgesdz2 = pgesdz
627  ! pges2d2 = pgesd2
628  ! eta = dicht + dist
629  ! CALL F_EOS_RN
630  ! pgesdz3 = pgesdz
631  ! pges2d3 = pgesd2
632  eta = dicht + 2.0*dist
633  CALL f_eos_rn
634  pgesdz4 = pgesdz
635  pges2d4 = pgesd2
636  eta = dicht
637  CALL f_eos_rn
638  pgesdz5 = pgesdz
639  pges2d5 = pgesd2
640 
641  ! pgesd3 = (-pgesdz4+16.0*pgesdz3-30.0*pgesdz5+16.0*pgesdz2-pgesdz1)
642  ! & /(12.0*(dist**2 ))
643  ! write(*,*) pgesd3
644  ! pgesd3 = (pges2d4-pges2d1)/(4.0*dist)
645  ! write(*,*) pgesd3
646  pgesd3 = (pgesdz4-2.0*pgesdz5+pgesdz1)/(2.0*dist)**2
647 
648 END SUBROUTINE p_eos_rn
649 
650 
651 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
652 !
653 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
654 
655 SUBROUTINE h_eos_rn
656 
657  USE parameters, ONLY: rgas
658  USE eos_variables
659  IMPLICIT NONE
660 
661  !-----------------------------------------------------------------------------
662  REAL :: dist,fact
663  REAL :: fres1,fres2,fres3,fres4,fres5
664  REAL :: cv_res,t_tmp,zges,df_dt, df_dtdt
665  !-----------------------------------------------------------------------
666 
667 
668  CALL perturbation_parameter
669  rho = eta/z3t
670 
671 
672  fact = 1.0
673  dist = t * 100.e-5 * fact
674 
675  t_tmp = t
676  t = t_tmp - 2.0*dist
677  CALL f_eos_rn
678  fres1 = fres
679  t = t_tmp - dist
680  CALL f_eos_rn
681  fres2 = fres
682  t = t_tmp + dist
683  CALL f_eos_rn
684  fres3 = fres
685  t = t_tmp + 2.0*dist
686  CALL f_eos_rn
687  fres4 = fres
688  t = t_tmp
689  CALL f_eos_rn
690  fres5 = fres
691  ! *(KBOL*T)/1.E-30
692 
693  zges = (p * 1.e-30) / (kbol*t*rho)
694 
695 
696  df_dt = (-fres4+8.0*fres3-8.0*fres2+fres1)/(12.0*dist)
697  df_dtdt = (-fres4+16.0*fres3-30.0*fres5+16.0*fres2-fres1) / (12.0*(dist**2 ))
698 
699 
700  s_res = (- df_dt -fres/t)*rgas*t + rgas *log(zges)
701  h_res = ( - t*df_dt + zges-1.0 ) * rgas*t
702  cv_res = -( t*df_dtdt + 2.0*df_dt ) * rgas*t
703 
704 END SUBROUTINE h_eos_rn
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
eos_constants
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:200
dft_module
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module DFT_MODULE This module...
Definition: modules.f90:272
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29