MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
DFT-nMF2.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! SUBROUTINE DFT_nMFT2
3 !
4 ! This subroutine performs a DFT calculation of the PCP-SAFT model using
5 ! the fundamental measure theory (FMT) for the hard-sphere fluid, the
6 ! TPT1-theory for chain formation and a first (or second) order perturba-
7 ! tion theory for the dispersive interactions.
8 !
9 ! Background to the code:
10 ! The functional derivative of d(F)/d(rho) is done in two parts, as
11 ! 'myrho(i) + mu_att(i)' , where 'mu_att' comprises the perturbation
12 ! theory (dispersive attraction) as well as the non-local part of the hard-
13 ! sphere term and chain term. For the hard-sphere term mu_att contains
14 ! contributions of the form mu_FMT = d(F^hs)/d(rho*) - (dF^hs/drho*)_LDA,
15 ! where the index 'hs' is for hard-sphere (analogous for the chain term)
16 ! and 'LDA' is for 'local density approximation', i.e. the values for the
17 ! local density.
18 !
19 ! The term 'myrho' contains the hard-sphere term and the chain term in
20 ! the LDA. Because the PC-SAFT dispersion term is not identical to the
21 ! Barker-Henderson (BH) perturtation theory (first or second order) for
22 ! LJ fluids, the term 'myrho' also contains the difference of the PC-SAFT
23 ! dispersion contribution and the BH perturbation theory used here the
24 ! LDA. This difference is obtained by using the flag {subtract1 = '1PT'}
25 ! the effect of the flag is that the BH-perturbation theory is subtracted
26 ! from the PC-SAFT model for calculating 'myrho'.
27 !
28 ! Caution: the indication for using the second order term in the pertur-
29 ! bation expansion is currently soley given by the flag {subtract1 = '2PT'}
30 ! The second order case needs to be re-tested before usage !
31 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
32 
33 SUBROUTINE dft_nmft2
34 
35  USE parameters, ONLY: pi, rgas, kbol
36  USE basic_variables
37  USE eos_variables, ONLY: fres, eta, rho, x, z3t
38  USE dft_module
39  use utilities
40  use eos_polar, only: f_polar, phi_polar
41  USE eos_numerical_derivatives
42 
43  IMPLICIT NONE
44 
45  !-----------------------------------------------------------------------------
46  INTERFACE
47  REAL FUNCTION dens_invers2 ( rhob, mu_rho )
48  IMPLICIT NONE
49  REAL, INTENT(IN) :: rhob(2)
50  REAL, INTENT(IN OUT) :: mu_rho
51  END FUNCTION dens_invers2
52 
53  SUBROUTINE dens_inv_coeff2 ( rhob )
54  IMPLICIT NONE
55  REAL, INTENT(IN) :: rhob(2)
56  END SUBROUTINE dens_inv_coeff2
57  END INTERFACE
58 
59  !-----------------------------------------------------------------------------
60  INTEGER :: i, j, kk, converg
61  INTEGER :: discret, fa, outpt, irc, maxiter, ih
62  LOGICAL :: diagram
63  REAL, DIMENSION(-NDFT:NDFT) :: zp, rhop, rhop_o
64  REAL, DIMENSION(-NDFT:NDFT) :: df_drho_tot, f_tot
65  REAL, DIMENSION(-NDFT:NDFT) :: df_drho_att, f_att
66  REAL, DIMENSION(2) :: rhob
67  REAL :: f_disp_1pt, f_disp_pcsaft
68  REAL :: mu_disp_1pt, mu_disp_pcsaft
69  REAL :: zges(np), pbulk
70  REAL :: msegm, delta_st, tanhfac
71  REAL :: end_x, steps, my0
72  REAL :: f_int_z, mu_int_z
73  REAL :: f_int_r, mu_int_r
74  REAL :: f_int2_z, mu_int2_z, mu_int3_z
75  REAL :: f_int2_r, mu_int2_r, mu_int3_r
76  REAL :: zz1
77  REAL :: dev, maxdev
78  REAL :: delta_f, free
79  REAL :: surftens(0:200), st_macro(200)
80  REAL :: f01, f02, f03, f04, f05
81  REAL :: c1_con, c2_con
82  REAL :: zms
83  REAL :: tsav, psav, tc, pc, rhoc
84  REAL :: density(np), w(np,nc), lnphi(np,nc)
85  REAL :: damppic
86  REAL :: box_l_no_unit
87  REAL, DIMENSION(-NDFT:NDFT) :: lambda, rhobar, phi_dn0, phi_dn1
88  REAL, DIMENSION(-NDFT:NDFT) :: phi_dn2, phi_dn3, phi_dn4, phi_dn5
89  REAL, DIMENSION(0:5,-NDFT:NDFT) :: ni
90  REAL :: f_ch, df_drho_ch
91  REAL :: dlngijdr, gij
92  REAL :: z0, z1, z2, z3
93  REAL :: f_fmt, df_drho_fmt
94  REAL :: pn_pt(0:ndft), zs
95  REAL :: mu_assoc, f_assoc
96  REAL :: fres_polar, fdd, fqq, fdq
97  REAL :: mu_polar, fdd_rk, fqq_rk, fdq_rk, z3_rk
98  !-----------------------------------------------------------------------------
99  ! note: variable ensemble_flag = 'tp' is set for calculating: mu_res=mu_res(T,p)
100  ! note: variable ensemble_flag = 'tv' is set for calculating: mu_res=mu_res(T,rho)
101  ! note: variable subtract1 = 'no' is set for regular calculations
102  ! note: variable subtract1 = '1PT' is set for calculating properties subtracting a
103  ! first order perturbation term (dispersion)
104 
105  num = 1
106  CALL set_default_eos_numerical
107 
108 
109  wca = .false.
110  shift = .false.
111  diagram = .true.
112 
113  CALL read_input
114  msegm = parame(1,1)
115 
116 
117  z3t_st = pi/6.0* parame(1,1) * ( 1.0 & ! divided by parame(i,2)**3
118  * ( 1.0 - 0.12*exp( -3.0*parame(1,3)/t ) ) )**3
119  d_hs = ( 1.0 - 0.12*exp( -3.0*parame(1,3)/t ) ) ! divided by parame(i,2)
120  dhs_st = d_hs
121 
122  IF ( ncomp /= 1 ) THEN
123  write (*,*)'SPECIFY ONLY ONE COMPONENT IN THE INPUT-FILE:'
124  write (*,*)' ./input_file/INPUT.INP'
125  stop
126  END IF
127  OPEN (68,file = './output_file/DFT_profiles.xlo')
128  OPEN (69,file = './output_file/DFT_sigma.xlo')
129  OPEN (71,file = './output_file/DFT_iteration.xlo')
130  OPEN (72,file = './output_file/DFT_h.xlo')
131 
132 
133  !-----------------------------------------------------------------------------
134  ! The the cut-off distance rc is the max. distance for integrating inter-
135  ! actions. Beyond rc, tail-corrections are added.
136  !-----------------------------------------------------------------------------
137  rc = 9.0
138 
139 
140  !-----------------------------------------------------------------------------
141  ! basic definitions for calculating density profile,
142  !-----------------------------------------------------------------------------
143  ! grid-size dzp = zp(1)-zp(0)
144 
145  box_l_no_unit = 50.0 ! lenth of simulation box (dimensionless, L/sigma)
146  discret= 1000 ! number of spacial discretizations for profile
147  dzp = box_l_no_unit / REAL(discret) ! dimensionless grid-distance dzp*=dzp/sigma
148  fa = nint( 1.0 / dzp ) ! size of the box is discret/fa * sigma (sigma: LJ diameter)
149  WRITE (*,*) fa, 1.0/dzp
150  outpt = 100 ! number output-points = discret/outpt
151 
152 
153  !-----------------------------------------------------------------------------
154  ! definitions for the numerical algorithm
155  !-----------------------------------------------------------------------------
156 
157  maxiter = 180 ! maximum number of iterations per temp.step
158  maxdev = 2.e-7 ! maximum deviation
159  damppic = 0.002 ! damping of Picard-iteration
160 
161 
162  !-----------------------------------------------------------------------------
163  ! get a matrix of values for the pair correlation function g(r,rho)
164  !-----------------------------------------------------------------------------
165 
166  rg = 4.0
167  den_step = 40
168  dzr = dzp / 2.0 ! dimensionless grid-distance for integrating
169  ! the Barker-Henderson attraction term
170  CALL rdf_matrix (msegm)
171 
172 
173  !-----------------------------------------------------------------------------
174  ! prepare for phase equilibrium calculation for given T
175  !-----------------------------------------------------------------------------
176 
177  diagram_for_various_t_loop: DO
178 
179  nphas = 2
180  n_unkw = ncomp ! number of quantities to be iterated
181  it(1) = 'p' ! iteration of pressure
182 
183  val_init(1) = 0.45 ! starting value for liquid density
184  val_init(2) = 1.e-5 ! starting value for vapor density
185  val_init(3) = t ! value of temperature
186  val_init(4) = 10.0 ! default starting value for p in (Pa)
187  IF ( eos >= 4 ) val_init(4) = 1.e-3 ! default starting value for LJ-model
188  val_init(5) = 0.0 ! logar. mole fraction: lnx=0, x=1, phase 1
189  val_init(6) = 0.0 ! logar. mole fraction: lnx=0, x=1, phase 2
190 
191  running = 't' ! T is running variable in PHASE_EQUILIB - here T=const.
192  end_x = t ! end temperature - here T=const.
193  steps = 1.0 ! number of steps of running var. - here 1, since T=const.
194 
195  outp = 0 ! output to terminal
196 
197 
198  !--------------------------------------------------------------------------
199  ! calculate phase equilibrium for given T
200  !--------------------------------------------------------------------------
201 
202  converg = 0
203  DO WHILE ( converg == 0 )
204  ensemble_flag = 'tp' ! this flag is for: 'regular calculation'
205  CALL phase_equilib( end_x, steps, converg )
206  IF ( converg /= 1 ) THEN
207  val_conv(2) = 0.0
208  val_init(3) = t - 10.0 ! value of temperature
209  IF ( val_conv(4) /= 0.0 ) val_init(4) = val_conv(4) / 4.0
210  ! caution: the intermediate points calculated here are not correct, because
211  ! the temperature-dependent props like d_BH, g(r), or z3t are not evaluated
212  ! for every temp-step. The final converged result, however, is correct.
213  steps = 2.0 ! number of steps of running var.
214  WRITE (*,*) 'no VLE found'
215  IF ( val_init(3) <= 10.0 ) RETURN
216  END IF
217  END DO
218 
219  ! rhob: molecular density times sigma**3 (rho_molec*sigma**3)
220  rhob(1) = dense(1)/z3t_st ! coexisting bulk density L
221  rhob(2) = dense(2)/z3t_st ! coexisting bulk density V
222  WRITE (*,*) 'temperature ',t,p/1.e5
223  WRITE (*,*) 'densities ',rhob(1), rhob(2)
224  WRITE (*,*) 'dense ',dense(1),dense(2)
225 
226 
227  !--------------------------------------------------------------------------
228  ! get density in SI-units (kg/m**3)
229  !--------------------------------------------------------------------------
230  CALL si_dens ( density, w )
231 
232 
233  !--------------------------------------------------------------------------
234  ! (re-)calculate residual chemical potential of both phases
235  !--------------------------------------------------------------------------
236  ensemble_flag = 'tv' ! this flag is for: mu_res=mu_res(T,rho)
237  densta(1) = dense(1) ! Index 1 is for liquid density (here: packing fraction eta)
238  densta(2) = dense(2) ! Index 2 is for vapour density (here: packing fraction eta)
239  CALL fugacity (lnphi)
240  my0 = lnphi(1,1) + log( rhob(1)/parame(1,2)**3 ) ! my0 = mu_res(T,rho_bulk_L) + ln(rho_bulk_l)
241  zges(1) = p / ( rgas*t*density(1) ) * mm(1)/1000.0
242  zges(2) = p / ( rgas*t*density(2) ) * mm(1)/1000.0
243 
244  pbulk = zges(1) * rhob(1) ! pressure p/kT (= Z*rho)
245  WRITE (*,*) 'chem. potentials', lnphi(1,1) + log( rhob(1)/parame(1,2)**3 ), &
246  lnphi(2,1) + log( rhob(2)/parame(1,2)**3 )
247  WRITE (*,*) ' '
248  tc = t
249  IF (num == 2) WRITE (*,*) 'enter an estimate for crit. temp.'
250  IF (num == 2) READ (*,*) tc
251  ensemble_flag = 'tp' ! this flag is for: 'regular calculation'
252  tsav = t
253  psav = p
254  IF ( eos < 4 ) CALL critical (tc,pc,rhoc)
255  !IF ( eos >= 4 .AND. eos /= 7 ) CALL lj_critical (tc,pc,rhoc)
256  t = tsav
257  p = psav
258  WRITE (*,'(a,3(f16.4))') 'critical point',tc, pc/1.e5, rhoc
259  WRITE (*,*) ' '
260 
261 
262  !--------------------------------------------------------------------------
263  ! define initial density profile rhop(i)
264  ! and dimensionless space coordinates zp(i)
265  !
266  ! discret : number of grid-points within "the box"
267  ! irc : number of grid-points extending the the box to left and right
268  ! the box is extended in order to allow for numerical integration
269  ! 'irc' is determined by the cut-off distance 'rc'
270  !--------------------------------------------------------------------------
271 
272  irc = nint(rc/dzp) + 1
273  tanhfac = -2.3625*t/tc + 2.4728
274  ! tanhfac = 2.7*(1.0 - t/tc) ! this parameterization was published (Gross, JCP, 2009)
275  DO i = -irc, (discret+irc)
276  zp(i) = REAL(i) * dzp
277  END DO
278  DO i = -irc, (discret+irc)
279  rhop(i) = tanh(-(zp(i)-zp(int(discret/2)))*tanhfac) * (rhob(1)-rhob(2))/2.0 &
280  + (rhob(1)+rhob(2))/2.0
281  ! rhop(i) = rhob(1)
282  END DO
283 
284 
285  !--------------------------------------------------------------------------
286  ! Initialize the DENS_INV subroutine
287  !--------------------------------------------------------------------------
288 
289  nphas = 1
290 
291  ensemble_flag = 'tv' ! this flag is for: mu_res=mu_res(T,rho)
292 
293 
294 
295  !==========================================================================
296  ! Start iterating the density profile
297  !==========================================================================
298 
299  kk = 1
300  ih = 85
301 
302  dft_convergence_loop: DO
303 
304  !-----------------------------------------------------------------------
305  ! Getting auxilliary quantities along the profile
306  !-----------------------------------------------------------------------
307 
308  CALL aux_chain ( discret, fa, dzp, d_hs, zp, rhop, rhobar, lambda )
309 
310  CALL fmt_dens ( discret, fa, dzp, d_hs, zp, rhop, ni, &
311  phi_dn0, phi_dn1, phi_dn2, phi_dn3, phi_dn4, phi_dn5 )
312 
313 
314  !-----------------------------------------------------------------------
315  ! Start loop for density profile
316  !-----------------------------------------------------------------------
317 
318  DO i = 0, discret
319 
320  ! ---------------------------------------------------------------
321  ! some auxiliary quantities used for compensating LDA contributions
322  ! ---------------------------------------------------------------
323 
324  z0 = rhop(i) * pi/6.0 * parame(1,1)
325  z1 = rhop(i) * pi/6.0 * parame(1,1) * d_hs
326  z2 = rhop(i) * pi/6.0 * parame(1,1) * d_hs**2
327  z3 = rhop(i) * pi/6.0 * parame(1,1) * d_hs**3
328  zms = 1.0 - z3
329  gij = 1.0/zms + 3.0*(d_hs/2.0)*z2/zms/zms + 2.0*((d_hs/2.0)*z2)**2 /zms**3
330  dlngijdr = (1.0/zms/zms + 3.0*(d_hs/2.0)*z2*(1.0+z3)/z3/zms**3 &
331  + ((d_hs/2.0)*z2/zms/zms)**2 *(4.0+2.0*z3)/z3 ) /gij*z3t_st
332 
333  ! ---------------------------------------------------------------
334  ! hard sphere contribution
335  ! f_fmt is the Helmholtz energy density F_FMT/(VkT) = F_FMT/(NkT)*rho
336  ! dF_drho_fmt is the functional derivative of F_FMT/(kT) to rho
337  ! ---------------------------------------------------------------
338  zs = 1.0 - ni(3,i)
339  f_fmt = - ni(0,i)*log(zs) + ni(1,i)*ni(2,i)/zs - ni(4,i)*ni(5,i)/zs &
340  + (ni(2,i)**3 -3.0*ni(2,i)*ni(5,i)*ni(5,i)) *(ni(3,i)+zs*zs*log(zs)) &
341  /36.0/pi/zs/zs/ni(3,i)**2
342 
343  CALL df_fmt_dr ( i, fa, dzp, d_hs, zp, phi_dn0, &
344  phi_dn1, phi_dn2, phi_dn3, phi_dn4, phi_dn5, df_drho_fmt )
345 
346 
347  ! ---------------------------------------------------------------
348  ! chain term
349  ! ---------------------------------------------------------------
350  CALL df_chain_dr ( i, fa, dzp, d_hs, zp, rhop, lambda, &
351  rhobar, z3t_st, f_ch, df_drho_ch )
352 
353 
354  ! ---------------------------------------------------------------
355  ! dispersive attraction
356  ! ---------------------------------------------------------------
357  f_int_z = 0.0
358  f_int2_z = 0.0
359  mu_int_z = 0.0
360  mu_int2_z = 0.0
361  mu_int3_z = 0.0
362  f01 = 0.0
363  f02 = 0.0
364  f03 = 0.0
365  f04 = 0.0
366  f05 = 0.0
367  DO j = (i-irc), (i+irc) ! integration in z-coordinate
368  zz1 = abs( zp(j) - zp(i) ) ! distance z12 between 1 and 2
369  IF ( zz1 <= rc-dzp+1.e-9 ) THEN
370  CALL dft_rad_int ( i, j, ih, zz1, rhop, f_int_r, mu_int_r, &
371  f_int2_r, mu_int2_r, mu_int3_r )
372  ELSE
373  f_int_r = 0.0
374  f_int2_r = 0.0
375  mu_int_r = 0.0
376  mu_int2_r= 0.0
377  mu_int3_r= 0.0
378  END IF
379 
380  f_int_z = f_int_z + dzp * (rhop(j)* f_int_r + f01)/2.0
381  mu_int_z = mu_int_z + dzp * (rhop(j)*mu_int_r + f02)/2.0
382  f_int2_z = f_int2_z + dzp * ( f_int2_r + f03)/2.0
383  mu_int2_z= mu_int2_z+ dzp * (mu_int2_r + f04)/2.0
384  mu_int3_z= mu_int3_z+ dzp * (mu_int3_r + f05)/2.0
385  f01 = rhop(j)* f_int_r
386  f02 = rhop(j)* mu_int_r
387  f03 = f_int2_r
388  f04 = mu_int2_r
389  f05 = mu_int3_r
390  END DO
391 
392  ! ---------------------------------------------------------------
393  ! cut-off corrections
394  ! ---------------------------------------------------------------
395 
396  CALL cutoff_corrections ( i, irc, rhop, rhob, f_int_z, mu_int_z, f_int2_z, mu_int2_z )
397 
398 
399  ! ---------------------------------------------------------------
400  ! for second order dispersive attraction ( usually not used !!! )
401  ! ---------------------------------------------------------------
402 
403  z3 = rhop(i) * z3t_st
404  zms = 1.0 - z3
405  c1_con= 1.0/ ( 1.0 + parame(1,1)*(8.0*z3-2.0*z3**2 )/zms**4 &
406  + (1.0 - parame(1,1))*(20.0*z3-27.0*z3**2 &
407  + 12.0*z3**3 -2.0*z3**4 ) /(zms*(2.0-z3))**2 )
408  c2_con= - c1_con*c1_con &
409  * ( parame(1,1)*(-4.0*z3**2 +20.0*z3+8.0)/zms**5 &
410  + (1.0 - parame(1,1)) *(2.0*z3**3 +12.0*z3**2 -48.0*z3+40.0) &
411  / (zms*(2.0-z3))**3 )
412  mu_int2_z = mu_int2_z /4.0 * ( 2.0*rhop(i)*c1_con + rhop(i)*z3*c2_con )
413  mu_int3_z = mu_int3_z /4.0 *rhop(i)*rhop(i)*c1_con
414  f_int2_z = f_int2_z /4.0 *rhop(i)*rhop(i)*c1_con
415 
416  ! ---------------------------------------------------------------
417  ! The integration of the DFT-integral can be done with cubic splines
418  ! ---------------------------------------------------------------
419  ! stepno = discret + irc*2
420  ! CALL SPLINE_PARA (dzp, intgrid, utri, stepno)
421  ! CALL SPLINE_INT (f_int_z, dzp, intgrid, utri, stepno)
422  ! f_int_z = f_int_z + rhob(1)*( 4.0/90.0*rc**-9 -1.0/3.0*rc**-3 )
423  ! f_int_z = f_int_z + rhob(2)*( 4.0/90.0*rc**-9 -1.0/3.0*rc**-3 )
424 
425 
426  ! ---------------------------------------------------------------
427  ! Calculation of 'dF_drho_att', which denotes non-local contributions
428  ! of d(F)/d(rho*). Apart from the dispersive attraction, it contains
429  ! non-local contributions of the hard-sphere term, and the chain term.
430  ! 'dF_drho_att' is defined as: dF_drho_att = d(F)/d(rho*) - (dF/drho*)_LDA
431  ! where LDA is for 'local density approximation', i.e. the values for
432  ! the local density.
433  !
434  ! F is the intrinsic Helmholtz energy
435  ! rho* = rhop(i) denotes the dimensionless density.
436  ! parame(1,3) is epsilon/k (LJ-energy parameter)
437  ! parame(1,2) is sigma (LJ-size parameter)
438  !
439  ! For a non-local second order term uncomment lines starting with '!2'
440  ! ---------------------------------------------------------------
441  ! remember: factor 2*PI because of the cylindrical coordinates
442 
443  df_drho_att(i) = 2.0*pi*parame(1,1)**2 * parame(1,3)/t * mu_int_z
444  f_att(i) = pi*parame(1,1)**2 *parame(1,3)/t*rhop(i)*f_int_z
445 
446 
447  ! --- if a second order perturbation theory is considered -------
448  IF (subtract1 =='2PT') THEN
449  df_drho_att(i)= df_drho_att(i) -2.0*pi*parame(1,1)**3 *(parame(1,3)/t)**2 *(mu_int2_z+mu_int3_z)
450  f_att(i) = f_att(i) -2.0*pi*parame(1,1)**3 *(parame(1,3)/t)**2 *f_int2_z
451  END IF
452 
453  ! ---------------------------------------------------------------
454  ! make the dispersive attraction term consistent with PC-SAFT, by
455  ! adding the difference ( PC-SAFT - 1PT )_dispersion locally (LDA)
456  ! ---------------------------------------------------------------
457 
458  !*****************************************************************************
459  ! under construction
460  !*****************************************************************************
461  dense(1) = rhop(i) * z3t_st
462  densta(1) = dense(1)
463 
464  ensemble_flag = 'tv'
465 
466  call only_one_term_eos_numerical ( 'disp_term', 'PT1 ' )
467  call fugacity ( lnphi )
468  call restore_previous_eos_numerical
469  !call f_PT (..... f_disp_1PT)
470  !call mu_PT ( .... mu_disp_1PT)
471  f_disp_1pt = fres*rhop(i)
472  mu_disp_1pt = lnphi(1,1)
473 
474  call only_one_term_eos_numerical ( 'disp_term', 'PC-SAFT ' )
475  call fugacity ( lnphi )
476  call restore_previous_eos_numerical
477  f_disp_pcsaft = fres*rhop(i)
478  mu_disp_pcsaft = lnphi(1,1)
479 
480  eta = densta(1)
481  rho = eta/z3t
482  z3_rk = z3t ! only for pure substances
483  x(1) = 1.0
484  call f_polar ( fdd, fqq, fdq )
485  call phi_polar ( 1, z3_rk, fdd_rk, fqq_rk, fdq_rk )
486  fres_polar = ( fdd + fqq + fdq ) * rhop(i)
487  mu_polar = fdd_rk + fqq_rk + fdq_rk
488 
489  f_assoc = 0.0
490  mu_assoc = 0.0
491  IF ( sum( nint(parame(1:ncomp,12)) ) > 0) THEN
492  call only_one_term_eos_numerical ( 'hb_term ', 'TPT1_Chap' )
493  call fugacity ( lnphi )
494  call restore_previous_eos_numerical
495  f_assoc = fres*rhop(i)
496  mu_assoc = lnphi(1,1)
497  END IF
498 
499  ensemble_flag = 'tp'
500 
501  df_drho_att(i)= df_drho_att(i) + ( mu_disp_pcsaft - mu_disp_1pt ) + mu_assoc + mu_polar
502  f_att(i) = f_att(i) + ( f_disp_pcsaft - f_disp_1pt ) + f_assoc + fres_polar
503  !*****************************************************************************
504  !*****************************************************************************
505 
506 
507  ! ---------------------------------------------------------------
508  ! apart from the perturbation theory, dF_drho_att and f_att also
509  ! collect the non-local parts of the hs and chain fluid
510  ! ---------------------------------------------------------------
511 
512  df_drho_tot(i) = log( rhop(i)/parame(1,2)**3 ) + df_drho_fmt + df_drho_ch + df_drho_att(i)
513  f_tot(i) = f_fmt + f_ch + f_att(i)
514 
515 
516  ! --- on-the-fly report on local errors in the profile ----------
517  IF ( mod(i, outpt) == 0 )write(*,'(i7,2(G15.6))') i,rhop(i), my0 - df_drho_tot(i)
518 
519  END DO ! end of loop (i = 0, discret) along the profile
520 
521 
522  !-----------------------------------------------------------------------
523  ! update the density profile using either a Picard-iteration scheme
524  ! (i.e. a direct substitution scheme, or a LDA - Inversion Procedure
525  ! similar to R.Evans (Bristol), which is done in function DENS_INVERS2.
526  !-----------------------------------------------------------------------
527  dev = 0.0
528  DO i = 0, discret
529  rhop_o(i) = rhop(i)
530  rhop(i) = rhop(i) * exp( my0 - df_drho_tot(i) )
531  rhop(i) = rhop_o(i) + (rhop(i)-rhop_o(i))*damppic
532  dev = dev + (rhop(i)-rhop_o(i))**2
533  END DO
534 
535 
536  !-----------------------------------------------------------------------
537  ! calculate surface tension
538  !-----------------------------------------------------------------------
539  free = 0.0
540  DO i = 1, discret
541  delta_f = f_tot(i) + rhop(i)*(log(rhop(i)/parame(1,2)**3 )-1.0) - (rhop(i)*my0 - pbulk) ! all quantities .../(kT)
542  pn_pt(i) = delta_f
543  free = free + delta_f*dzp
544  END DO
545  surftens(kk) = kbol * t *1.e20*1000.0 *free / parame(1,2)**2
546 
547 
548  !-----------------------------------------------------------------------
549  ! add an approximate capillary wave contribution to the surface tension
550  !-----------------------------------------------------------------------
551  st_macro(kk) = surftens(kk) / ( 1.0 + 3.0/8.0/pi *t/tc &
552  * (1.0/2.55)**2 / (0.0674*parame(1,1)+0.0045) )
553 
554 
555  delta_st = 1.0
556  IF ( kk > 1 ) delta_st = abs( surftens(kk)-surftens(kk-1) ) / surftens(kk)
557 
558  WRITE (*,*) '-----------------------------------------------------------'
559  WRITE (*,*) ' # error-sum intrinsic ST total ST'
560  WRITE (*,'(i4,3F16.8)') kk, dev, surftens(kk), st_macro(kk)
561  WRITE (*,*) '-----------------------------------------------------------'
562  WRITE (71,'(i4,4G18.8)') kk, dev, surftens(kk), st_macro(kk)
563  kk = kk + 1
564 
565  !-----------------------------------------------------------------------
566  ! add an approximate capillary wave contribution to the surface tension
567  !-----------------------------------------------------------------------
568  ! write (*,*) 'error measure',dev, delta_st
569  IF ( dev < maxdev .OR. delta_st < 1.e-8 ) THEN
570  EXIT dft_convergence_loop
571  ELSE
572  IF ( kk > maxiter ) THEN
573  WRITE(*,*)' no convergence in ',maxiter,' steps'
574  EXIT dft_convergence_loop
575  ENDIF
576  ENDIF
577 
578  ENDDO dft_convergence_loop
579 
580  !--------------------------------------------------------------------------
581  ! write resulting density profile
582  !--------------------------------------------------------------------------
583  WRITE(68,'(i2,6f18.8)') 0, zp(0)-zp(int(discret/2)), rhop(0), t, &
584  val_conv(4)/1.e5, density(1), density(2)
585  DO i = 0, discret
586  IF ( mod(i, outpt) == 0 ) WRITE (68,'(i6,3(f18.12))') i,zp(i)-zp(int(discret/2)), &
587  rhop(i), -( df_drho_tot(i)-my0 )
588  ! IF ( MOD(i, outpt) == 0 ) WRITE(*,'(i6,3(f18.12))') i,zp(i)-zp(INT(discret/2)),rhop(i)
589  ! write (69,*) zp(i)*parame(1,2), pN_pT(i)/parame(1,2)**3 /1E-30*parame(1,3)*KBOL/1.E6 ! in MPa
590  END DO
591  WRITE (68,*) ' '
592 
593 
594  !--------------------------------------------------------------------------
595  ! summarize results
596  !--------------------------------------------------------------------------
597  WRITE (*,*) ' '
598  WRITE (*,*) 'SUMMARY FOR A SINGLE TEMPERATURE'
599  WRITE (*,*) ' '
600  WRITE (*,*) 'Temp. [K], Pressure [bar] ',t,val_conv(4)/1.e5
601  WRITE (*,*) 'Critical point Temp., Press. ',tc,pc/1.e5
602  WRITE (*,*) 'Density [kg/m**3] ',density(1),density(2)
603  WRITE (*,*) 'Dimensionless Density (rho*) ',rhob(1),rhob(2)
604  WRITE (*,*) 'Excess Grand Potential ',free
605  WRITE (*,*) 'Intrinsic Interf. Tension [mN/m] ',surftens(kk-1)
606  WRITE (*,*) 'Macroscop.Interf. Tension [mN/m] ',st_macro(kk-1)
607  WRITE (*,*) '============================================================'
608  WRITE (*,*) ' '
609  WRITE (69,'(9(f18.10))') t, val_conv(4)/1.e5, &
610  rhob(1),rhob(2),surftens(kk-1),st_macro(kk-1),free,dev
611 
612  !--------------------------------------------------------------------------
613  ! when calc. a phase diagram & diagram of surface tens., loop for higher T
614  !--------------------------------------------------------------------------
615  ensemble_flag = 'tp' ! this flag is for 'regular calculations'
616  subtract1 = 'no' ! this flag is for 'regular calculations'
617  IF ( diagram ) THEN
618  IF ( (t+8.0) <= tc ) THEN
619  t = t + 5.0
620  IF ( (t+15.0) <= tc ) t = t + 5.0
621  IF ( (t+25.0) <= tc ) t = t + 10.0
622  IF ( (t+45.0) <= tc ) t = t + 20.0
623  nphas = 2
624  n_unkw = ncomp ! number of quantities to be iterated
625  it(1) = 'p' ! iteration of pressure
626  val_init(3) = t ! value of temperature
627  running = 't' ! T is running variable in PHASE_EQUILIB - here T=const.
628  end_x = t ! end temperature - here T=const.
629  steps = 1.0 ! number of steps of running var. - here 1, since T=const.
630  z3t_st = pi/6.0 * parame(1,1) * ( 1.0 & ! divided by parame(i,2)**3
631  *( 1.0 - 0.12*exp(-3.0*parame(1,3)/t) ) )**3
632  d_hs = ( 1.0 - 0.12*exp(-3.0*parame(1,3)/t) ) ! divided by parame(i,2)
633  dhs_st = d_hs
634  CALL rdf_matrix (msegm)
635  ELSE
636  EXIT diagram_for_various_t_loop
637  END IF
638  WRITE (69,'(7(f18.10))') tc, pc/1.e5, rhoc, rhoc, 0., 0., 0.
639  ELSE
640  EXIT diagram_for_various_t_loop
641  END IF
642 
643  ENDDO diagram_for_various_t_loop
644 
645 END SUBROUTINE dft_nmft2
646 
647 
648 
649 
650 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
651 !
652 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
653 
654 SUBROUTINE cutoff_corrections ( i, irc, rhop, rhob, f_int_z, mu_int_z, f_int2_z, mu_int2_z )
655 
656  USE dft_module, ONLY: rc, ndft
657  IMPLICIT NONE
658 
659  !-----------------------------------------------------------------------------
660  INTEGER, INTENT(IN) :: i
661  INTEGER, INTENT(IN) :: irc
662  REAL, DIMENSION(-NDFT:NDFT), INTENT(IN) :: rhop
663  REAL, DIMENSION(2), INTENT(IN) :: rhob
664  REAL, INTENT(IN OUT) :: f_int_z, f_int2_z
665  REAL, INTENT(IN OUT) :: mu_int_z, mu_int2_z
666 
667  !-----------------------------------------------------------------------------
668  REAL :: cutoffz, cutoffz2, rhotemp
669  !-----------------------------------------------------------------------------
670 
671  cutoffz = 4.0/90.0 * rc**-9 - 1.0/3.0 * rc**-3
672  cutoffz2 = 16.0/22.0/21.0 * rc**-21 - 2.0/15.0 * rc**-15 + 16.0/90.0 * rc**-9
673  IF ( abs( rhop(i+irc)-rhob(2) ) > 1.e-5 ) THEN
674  rhotemp = rhop(i+irc) ! this is a crude approximation !
675  ! rhotemp = rhob(2)
676  f_int_z = f_int_z + rhotemp*cutoffz
677  mu_int_z = mu_int_z + rhotemp*cutoffz
678  f_int2_z = f_int2_z + cutoffz2
679  mu_int2_z= mu_int2_z+ cutoffz2
680  ELSE
681  f_int_z = f_int_z + rhob(2)*cutoffz
682  mu_int_z = mu_int_z + rhob(2)*cutoffz
683  f_int2_z = f_int2_z + cutoffz2
684  mu_int2_z= mu_int2_z+ cutoffz2
685  END IF
686  IF ( abs( rhop(i-irc)-rhob(1) ) > 1.e-3 ) THEN
687  rhotemp = rhop(i-irc)
688  ! rhotemp=rhob(1)
689  f_int_z = f_int_z + rhotemp*cutoffz
690  mu_int_z = mu_int_z + rhotemp*cutoffz
691  f_int2_z = f_int2_z + cutoffz2
692  mu_int2_z= mu_int2_z+ cutoffz2
693  ELSE
694  f_int_z = f_int_z + rhob(1)*cutoffz
695  mu_int_z = mu_int_z + rhob(1)*cutoffz
696  f_int2_z = f_int2_z + cutoffz2
697  mu_int2_z= mu_int2_z+ cutoffz2
698  END IF
699 
700 END SUBROUTINE cutoff_corrections
701 
702 
703 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
704 !
705 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
706 
707 SUBROUTINE aux_chain (discret,fa,dzp,d_hs,zp,rhop,rhobar,lambda)
708 
709  USE basic_variables
710  USE dft_module, ONLY: ndft
711  IMPLICIT NONE
712 
713  !-----------------------------------------------------------------------------
714  INTEGER, INTENT(IN) :: discret
715  INTEGER, INTENT(IN) :: fa
716  REAL, INTENT(IN) :: dzp
717  REAL, INTENT(IN) :: d_hs
718  REAL, INTENT(IN) :: zp(-ndft:ndft)
719  REAL, INTENT(IN) :: rhop(-ndft:ndft)
720  REAL, INTENT(OUT) :: rhobar(-ndft:ndft)
721  REAL, INTENT(OUT) :: lambda(-ndft:ndft)
722 
723  !-----------------------------------------------------------------------------
724  INTEGER :: i, j, nn
725  REAL :: int1, int2, zz1, xl, xh, dz = 0.0
726  REAL, DIMENSION(100) :: y2, rx, ry1, ry2
727  !-----------------------------------------------------------------------------
728 
729  ! DO i=(-fa),(discret+fa)
730  ! lambda(i) = 0.0
731  ! rhobar(i) = 0.0
732  ! f01 = 0.0 ! = 0.75*rhop(i-fa)*( 1.**2 -1.**2 )
733  ! f02 = 0.5 *rhop(i-fa)
734  ! DO j=(i-fa+1),(i+fa) ! fa=NINT(1./dzp)
735  ! zz1 = ABS(zp(j)-zp(i)) ! distance z12 between 1 and 2
736  ! int_r= 0.75*rhop(j)*(1.0-zz1*zz1)
737  ! int_l= 0.5 *rhop(j)
738  ! rhobar(i) = rhobar(i) + dzp * ( int_r+f01 )/2.0
739  ! lambda(i) = lambda(i) + dzp * ( int_l+f02 )/2.0
740  ! f01=int_r
741  ! f02=int_l
742  ! ENDDO
743  ! write(69,*) i,rhobar(i)
744  ! ENDDO
745 
746  DO i = (-fa), (discret+fa)
747  nn = 1
748  DO j = (i-fa), (i+fa) ! fa=NINT(1./dzp)
749  IF ( zp(j+1) > (zp(i)-d_hs).AND.(zp(i)-d_hs) >= zp(j) ) THEN
750  dz = zp(j+1) - (zp(i)-d_hs)
751  rx(1) = 0.0
752  ry1(1) = 0.0 ! = 0.75*rhop(i-fa)*( d.**2 -d.**2 )
753  ry2(1) = 0.5/d_hs*rhop(j) + ((zp(i)-d_hs)-zp(j))/dzp &
754  *( 0.5/d_hs*rhop(j+1) - 0.5/d_hs*rhop(j) )
755  ELSE IF ( zp(j) > (zp(i)-d_hs) .AND. zp(j) <= (zp(i)+d_hs) ) THEN
756  zz1 = abs( zp(j)-zp(i) ) ! distance z12 between 1 and 2
757  nn = nn + 1
758  ry1(nn) = 0.75/d_hs**3 * rhop(j) * (d_hs**2 -zz1*zz1)
759  ry2(nn) = 0.5/d_hs * rhop(j)
760  rx(nn) = rx(nn-1) + dz
761  dz = dzp
762  IF ( (zp(i)+d_hs) < zp(j+1) ) THEN
763  dz = (zp(i)+d_hs) - zp(j)
764  nn = nn + 1
765  rx(nn) = rx(nn-1) + dz
766  ry1(nn) = 0.0
767  ry2(nn) = 0.5/d_hs*rhop(j) + ((zp(i)+d_hs)-zp(j))/dzp &
768  *( 0.5/d_hs*rhop(j+1) - 0.5/d_hs*rhop(j) )
769  END IF
770  END IF
771  END DO
772  xl = rx(1)
773  xh = rx(nn)
774  CALL spline ( rx, ry1, nn, 1.e30, 1.e30, y2 )
775  CALL splint_integral ( rx, ry1, y2, nn, xl, xh, int1 )
776  rhobar(i) = int1
777  CALL spline ( rx, ry2, nn, 1.e30, 1.e30, y2 )
778  CALL splint_integral ( rx, ry2, y2, nn, xl, xh, int2 )
779  lambda(i) = int2
780  END DO
781 
782 END SUBROUTINE aux_chain
783 
784 
785 
786 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
787 !
788 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
789 
790 SUBROUTINE fmt_dens ( discret, fa, dzp, d_hs, zp, rhop, ni, &
791  phi_dn0, phi_dn1, phi_dn2, phi_dn3, phi_dn4, phi_dn5 )
792 
793  USE parameters, ONLY: pi
794  USE basic_variables
795  USE dft_module, ONLY: ndft
796  IMPLICIT NONE
797 
798  !-----------------------------------------------------------------------------
799  INTEGER, INTENT(IN OUT) :: discret
800  INTEGER, INTENT(IN) :: fa
801  REAL, INTENT(IN) :: dzp
802  REAL, INTENT(IN) :: d_hs
803  REAL, INTENT(IN) :: zp(-ndft:ndft)
804  REAL, INTENT(IN) :: rhop(-ndft:ndft)
805  REAL, INTENT(OUT) :: ni(0:5,-ndft:ndft)
806  REAL, INTENT(OUT) :: phi_dn0(-ndft:ndft)
807  REAL, INTENT(OUT) :: phi_dn1(-ndft:ndft)
808  REAL, INTENT(OUT) :: phi_dn2(-ndft:ndft)
809  REAL, INTENT(OUT) :: phi_dn3(-ndft:ndft)
810  REAL, INTENT(OUT) :: phi_dn4(-ndft:ndft)
811  REAL, INTENT(OUT) :: phi_dn5(-ndft:ndft)
812 
813  !-----------------------------------------------------------------------------
814  INTEGER :: i, j, nn
815  REAL :: int2, int3, int5
816  REAL :: zz1, zs, d2, xl, xh, dz = 0.0
817  REAL, DIMENSION(100) :: y2, hx, hy2, hy3, hy5
818  !-----------------------------------------------------------------------------
819 
820  DO i = (-fa), (discret+fa)
821  nn = 1
822  d2 = d_hs / 2.0
823  DO j = (i-fa/2), (i+fa/2)
824  IF ( zp(j+1) > (zp(i)-d2) .AND. (zp(i)-d2) >= zp(j) ) THEN
825  dz = zp(j+1) - (zp(i)-d2)
826  hx(1) = 0.0
827  hy2(1) = rhop(j) *parame(1,1) + ((zp(i)-d2)-zp(j))/dzp &
828  * ( rhop(j+1)*parame(1,1) - rhop(j)*parame(1,1) )
829  hy3(1) = 0.0 ! = rhop(i-fa/2)*parame(1,1)*(0.25-(zp(i-fa/2)-zp(i))**2 )
830  hy5(1) = rhop(j)*parame(1,1)*(zp(j)-zp(i)) +((zp(i)-d2)-zp(j))/dzp &
831  * ( rhop(j+1)*parame(1,1)*(zp(j+1)-zp(i)) &
832  - rhop(j) *parame(1,1)*(zp(j)-zp(i)) )
833  ELSE IF ( zp(j) > (zp(i)-d2) .AND. zp(j) <= (zp(i)+d2) ) THEN
834  zz1 = zp(j)-zp(i) ! distance z12 between 1 and 2
835  nn = nn + 1
836  hy2(nn) = rhop(j)*parame(1,1)
837  hy3(nn) = rhop(j)*parame(1,1) * ( d2**2 - zz1**2 )
838  hy5(nn) = rhop(j)*parame(1,1) * zz1
839  hx(nn) = hx(nn-1) + dz
840  dz = dzp
841  IF ( (zp(i)+d2) < zp(j+1) ) THEN
842  dz = (zp(i)+d2) - zp(j)
843  nn = nn + 1
844  hy2(nn) = rhop(j)*parame(1,1) + ((zp(i)+d2)-zp(j))/dzp &
845  * ( rhop(j+1)*parame(1,1) - rhop(j)*parame(1,1) )
846  hy3(nn) = 0.0
847  hy5(nn) = rhop(j)*parame(1,1)*(zp(j)-zp(i)) +((zp(i)+d2)-zp(j))/dzp &
848  * ( rhop(j+1)*parame(1,1)*(zp(j+1)-zp(i)) &
849  - rhop(j)*parame(1,1)*(zp(j) -zp(i)) )
850  hx(nn) = hx(nn-1) + dz
851  END IF
852  END IF
853  END DO
854  xl = hx(1)
855  xh = hx(nn)
856  CALL spline ( hx, hy2, nn, 1.e30, 1.e30, y2 )
857  CALL splint_integral ( hx, hy2, y2, nn, xl, xh, int2 )
858  CALL spline ( hx, hy3, nn, 1.e30, 1.e30, y2 )
859  CALL splint_integral ( hx, hy3, y2, nn, xl, xh, int3 )
860  CALL spline ( hx, hy5, nn, 1.e30, 1.e30, y2 )
861  CALL splint_integral ( hx, hy5, y2, nn, xl, xh, int5 )
862  ni(2,i) = pi * d_hs *int2 ! corresponds to z2*6
863  ni(3,i) = pi * int3 ! corresponds to z3
864  ni(5,i) = 2.0* pi *int5
865  ni(0,i) = ni(2,i) / (pi*d_hs**2 ) ! corresponds to z0/PI*6.0
866  ni(1,i) = ni(2,i) / (2.0*pi*d_hs) ! corresponds to z1/PI*3.0
867  ni(4,i) = ni(5,i) / (2.0*pi*d_hs)
868  END DO
869 
870  ! derivatives of phi to ni
871  DO i = (-fa), (discret+fa)
872  zs = 1.0 - ni(3,i)
873  phi_dn0(i) = - log(zs) !!!!
874  phi_dn1(i) = ni(2,i)/zs
875  phi_dn2(i) = ni(1,i)/zs + 3.0*(ni(2,i)*ni(2,i) - ni(5,i)*ni(5,i)) &
876  * (ni(3,i)+zs*zs*log(zs))/(36.0*pi*ni(3,i)*ni(3,i)*zs*zs)
877  phi_dn3(i) = ni(0,i)/zs + (ni(1,i)*ni(2,i)-ni(4,i)*ni(5,i))/zs/zs &
878  - (ni(2,i)**3 - 3.0*ni(2,i)*ni(5,i)*ni(5,i)) &
879  *( ni(3,i)*(ni(3,i)*ni(3,i)-5.0*ni(3,i)+2.0) &
880  + 2.0*zs**3 *log(zs) ) / ( 36.0*pi*(ni(3,i)*zs)**3 )
881  phi_dn4(i) = - ni(5,i)/zs
882  phi_dn5(i) = - ni(4,i)/zs - 6.0*ni(2,i)*ni(5,i) &
883  * (ni(3,i)+zs*zs*log(zs))/(36.0*pi*ni(3,i)*ni(3,i)*zs*zs)
884  END DO
885 
886 END SUBROUTINE fmt_dens
887 
888 
889 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
890 !
891 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
892 
893 SUBROUTINE df_fmt_dr ( i, fa, dzp, d_hs, zp, phi_dn0, &
894  phi_dn1, phi_dn2, phi_dn3, phi_dn4, phi_dn5, dF_drho_fmt )
895 
896  USE parameters, ONLY: pi
897  USE basic_variables
898  USE dft_module, ONLY: ndft
899  IMPLICIT NONE
900 
901  !-----------------------------------------------------------------------------
902  INTEGER, INTENT(IN) :: i
903  INTEGER, INTENT(IN) :: fa
904  REAL, INTENT(IN) :: dzp
905  REAL, INTENT(IN) :: d_hs
906  REAL, INTENT(IN) :: zp(-ndft:ndft)
907  REAL, INTENT(IN) :: phi_dn0(-ndft:ndft)
908  REAL, INTENT(IN) :: phi_dn1(-ndft:ndft)
909  REAL, INTENT(IN) :: phi_dn2(-ndft:ndft)
910  REAL, INTENT(IN) :: phi_dn3(-ndft:ndft)
911  REAL, INTENT(IN) :: phi_dn4(-ndft:ndft)
912  REAL, INTENT(IN) :: phi_dn5(-ndft:ndft)
913  REAL, INTENT(OUT) :: df_drho_fmt
914 
915  !-----------------------------------------------------------------------------
916  INTEGER :: j, nn
917  REAL :: int0, int1, int2, int3, int4, int5
918  REAL :: zz1, d2, xl, xh, dz = 0.0
919  REAL, DIMENSION(100) :: y2, hx, hy0, hy1, hy2, hy3, hy4, hy5
920  !-----------------------------------------------------------------------------
921 
922 
923  nn = 1
924  d2 = d_hs / 2.0
925  DO j = (i-fa/2), (i+fa/2)
926  IF ( zp(j+1) > (zp(i)-d2).AND.(zp(i)-d2) >= zp(j) ) THEN
927  hx(1) = 0.0
928  hy0(1) = phi_dn0(j) + ((zp(i)-d2)-zp(j))/dzp *( phi_dn0(j+1) - phi_dn0(j) )
929  hy1(1) = phi_dn1(j) + ((zp(i)-d2)-zp(j))/dzp *( phi_dn1(j+1) - phi_dn1(j) )
930  hy2(1) = phi_dn2(j) + ((zp(i)-d2)-zp(j))/dzp *( phi_dn2(j+1) - phi_dn2(j) )
931  hy3(1) = 0.0 ! = rhop(i-fa/2)*parame(1,1)*(0.25-(zp(i-fa/2)-zp(i))**2)
932  hy4(1) = phi_dn4(j)*(zp(j)-zp(i)) +((zp(i)-d2)-zp(j))/dzp &
933  * ( phi_dn4(j+1)*(zp(j+1)-zp(i)) - phi_dn4(j)*(zp(j)-zp(i)) )
934  hy5(1) = phi_dn5(j)*(zp(j)-zp(i)) +((zp(i)-d2)-zp(j))/dzp &
935  * ( phi_dn5(j+1)*(zp(j+1)-zp(i)) - phi_dn5(j)*(zp(j)-zp(i)) )
936  dz = zp(j+1) - (zp(i)-d2)
937  ELSE IF ( zp(j) > (zp(i)-d2).AND.zp(j) <= (zp(i)+d2) ) THEN
938  zz1 = zp(j) - zp(i) ! distance z12 between 1 and 2
939  nn = nn + 1
940  hy0(nn) = phi_dn0(j)
941  hy1(nn) = phi_dn1(j)
942  hy2(nn) = phi_dn2(j)
943  hy3(nn) = phi_dn3(j) * ( d2**2 -zz1**2 )
944  hy4(nn) = phi_dn4(j) * zz1
945  hy5(nn) = phi_dn5(j) * zz1
946  hx(nn) = hx(nn-1) + dz
947  dz = dzp
948  IF ( (zp(i)+d2) < zp(j+1) ) THEN
949  dz = (zp(i)+d2) - zp(j)
950  nn = nn + 1
951  hx(nn) = hx(nn-1) + dz
952  hy0(nn) = phi_dn0(j) + ( (zp(i)+d2) - zp(j) ) / dzp &
953  * ( phi_dn0(j+1) - phi_dn0(j) )
954  hy1(nn) = phi_dn1(j) + ( (zp(i)+d2) - zp(j) ) / dzp &
955  * ( phi_dn1(j+1) - phi_dn1(j) )
956  hy2(nn) = phi_dn2(j) + ( (zp(i)+d2) - zp(j) ) / dzp &
957  * ( phi_dn2(j+1) - phi_dn2(j) )
958  hy3(nn) = 0.0
959  hy4(nn) = phi_dn4(j) * (zp(j)-zp(i)) + ((zp(i)+d2)-zp(j)) / dzp &
960  * ( phi_dn4(j+1)*(zp(j+1)-zp(i)) -phi_dn4(j)*(zp(j)-zp(i)) )
961  hy5(nn) = phi_dn5(j) * (zp(j)-zp(i)) + ((zp(i)+d2)-zp(j))/dzp &
962  * ( phi_dn5(j+1)*(zp(j+1)-zp(i)) -phi_dn5(j)*(zp(j)-zp(i)) )
963  END IF
964  END IF
965  END DO
966  xl = hx(1)
967  xh = hx(nn)
968  CALL spline ( hx, hy0, nn, 1.e30, 1.e30, y2 )
969  CALL splint_integral( hx, hy0, y2, nn, xl, xh, int0 )
970  CALL spline ( hx, hy1, nn, 1.e30, 1.e30, y2 )
971  CALL splint_integral( hx, hy1, y2, nn, xl, xh, int1 )
972  CALL spline ( hx, hy2, nn, 1.e30, 1.e30, y2 )
973  CALL splint_integral( hx, hy2, y2, nn ,xl, xh, int2 )
974  CALL spline ( hx, hy3, nn, 1.e30, 1.e30, y2 )
975  CALL splint_integral( hx, hy3, y2, nn, xl, xh, int3 )
976  CALL spline ( hx, hy4, nn, 1.e30, 1.e30, y2 )
977  CALL splint_integral( hx, hy4, y2, nn, xl, xh, int4 )
978  CALL spline ( hx, hy5, nn, 1.e30, 1.e30, y2 )
979  CALL splint_integral( hx, hy5, y2, nn, xl, xh, int5 )
980  ! write (*,*) int3,int4,int5
981  ! call paus (' ')
982  ! int0 = int0/d_hs
983  ! int1 = 0.5*parame(1,2)*int1
984  ! int2 = PI*d_hs*parame(1,2)**2 *int2
985  ! int3 = PI*parame(1,2)**3 *int3
986  ! int4 = parame(1,2)*int4
987  ! int5 = 2.0*PI*parame(1,2)**2 *int5
988  int0 = int0 / d_hs
989  int1 = 0.5 * int1
990  int2 = pi * d_hs * int2
991  int3 = pi * int3
992  int4 = int4 / d_hs
993  int5 = 2.0 * pi * int5
994  ! dF_drho_fmt=dF_hs/drho = dF_hs/drho_segment *m_mean
995  df_drho_fmt = (int0+int1+int2+int3+int4+int5)*parame(1,1)
996 
997 END SUBROUTINE df_fmt_dr
998 
999 
1000 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1001 !
1002 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1003 
1004 SUBROUTINE df_chain_dr ( i, fa, dzp, d_hs, zp, rhop, lambda, &
1005  rhobar, z3t_st, f_ch, dF_drho_ch )
1006 
1007  USE parameters, ONLY: pi
1008  USE basic_variables
1009  USE dft_module, ONLY: ndft
1010  IMPLICIT NONE
1011 
1012  !-----------------------------------------------------------------------------
1013  INTEGER, INTENT(IN) :: i
1014  INTEGER :: fa
1015  REAL, INTENT(IN) :: dzp
1016  REAL, INTENT(IN) :: d_hs
1017  REAL, INTENT(IN) :: zp(-ndft:ndft)
1018  REAL, INTENT(IN) :: rhop(-ndft:ndft)
1019  REAL, INTENT(IN) :: lambda(-ndft:ndft)
1020  REAL, INTENT(IN) :: rhobar(-ndft:ndft)
1021  REAL, INTENT(IN OUT) :: z3t_st
1022  REAL, INTENT(OUT) :: f_ch
1023  REAL, INTENT(OUT) :: df_drho_ch
1024 
1025  !-----------------------------------------------------------------------------
1026  INTEGER :: j, nn
1027  REAL :: zz1, xl, xh, dz=0.0
1028  REAL, DIMENSION(100) :: y2, hx, hy0, hy1
1029  REAL :: ycorr, dlnydr, i_ch1, i_ch2
1030  !-----------------------------------------------------------------------------
1031 
1032 
1033  nn = 1
1034  DO j = (i-fa), (i+fa) ! fa=NINT(1./dzp)
1035  IF ( zp(j+1) > (zp(i)-d_hs) .AND. zp(j)-1.e-11 <= (zp(i)-d_hs) ) THEN
1036  dz = zp(j+1) - (zp(i)-d_hs)
1037  hx(1) = 0.0
1038  hy0(1) = 0.0 ! = 0.75*rhop(i-fa)*DLNYDR(d_hs,rhobar(i-fa))*(d.**2 -d.**2 )
1039  hy1(1) = 0.5/d_hs*rhop(j)/lambda(j) + ((zp(i)-d_hs)-zp(j))/dzp &
1040  * ( 0.5/d_hs*rhop(j+1)/lambda(j+1) - 0.5/d_hs*rhop(j)/lambda(j) )
1041  ELSE IF ( zp(j) > (zp(i)-d_hs).AND.zp(j)+1.e-11 < (zp(i)+d_hs) ) THEN
1042  zz1 = zp(j) - zp(i) ! distance z12 between 1 and 2
1043  nn = nn + 1
1044  CALL cavity ( d_hs, z3t_st, rhobar(j), ycorr, dlnydr )
1045  hy0(nn) = 0.75/d_hs**3 *rhop(j)*dlnydr*(d_hs*d_hs-zz1*zz1)
1046  hy1(nn) = 0.5/d_hs*rhop(j)/lambda(j)
1047  hx(nn) = hx(nn-1) + dz
1048  dz = dzp
1049  IF ( zp(j+1)+1.e-11 >= (zp(i)+d_hs) ) THEN
1050  dz = (zp(i)+d_hs) - zp(j)
1051  nn = nn + 1
1052  hy0(nn) = 0.0
1053  hy1(nn) = 0.5/d_hs*rhop(j)/lambda(j) +((zp(i)+d_hs)-zp(j))/dzp &
1054  * ( 0.5/d_hs*rhop(j+1)/lambda(j+1)-0.5/d_hs*rhop(j)/lambda(j) )
1055  hx(nn) = hx(nn-1) + dz
1056  END IF
1057  END IF
1058  END DO
1059  xl = hx(1)
1060  xh = hx(nn)
1061  CALL spline ( hx, hy0, nn, 1.e30, 1.e30, y2 )
1062  CALL splint_integral( hx, hy0, y2, nn, xl, xh, i_ch1 )
1063  CALL spline ( hx, hy1, nn, 1.e30, 1.e30, y2 )
1064  CALL splint_integral( hx, hy1, y2, nn, xl, xh, i_ch2 )
1065  ! write (*,*) '111',I_ch1,I_ch2
1066 
1067  ! DO j=1,nn
1068  ! write (*,*) j,hx(j),hy1(j)
1069  ! ENDDO
1070 
1071  ! I_ch1 = 0.0
1072  ! I_ch2 = 0.0
1073  ! DO j=1,nn-1
1074  ! I_ch1=I_ch1+(hy0(j)+hy0(j+1))/2.0*(hx(j+1)-hx(j))
1075  ! I_ch2=I_ch2+(hy1(j)+hy1(j+1))/2.0*(hx(j+1)-hx(j))
1076  ! ENDDO
1077  ! write (*,*) '222',I_ch1,I_ch2
1078 
1079  CALL cavity ( d_hs, z3t_st, rhobar(i), ycorr, dlnydr )
1080  f_ch = rhop(i) * ( log(ycorr*lambda(i)) -1.0 )
1081  f_ch = - ( parame(1,1) - 1.0 ) * f_ch
1082  f_ch = f_ch + ( parame(1,1) - 1.0 ) * rhop(i) * ( log(rhop(i)) - 1.0 )
1083  df_drho_ch= log( ycorr*lambda(i) ) - 1.0 + i_ch1 + i_ch2
1084  df_drho_ch= - ( parame(1,1) - 1.0 ) * df_drho_ch
1085  df_drho_ch= df_drho_ch + ( parame(1,1) - 1.0 ) * log(rhop(i))
1086 
1087 END SUBROUTINE df_chain_dr
1088 
1089 
1090 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1091 !
1092 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1093 
1094 SUBROUTINE cavity ( d_hs, z3t_st, rho, ycorr, dlnydr )
1095 
1096  USE parameters, ONLY: pi
1097  USE basic_variables
1098  IMPLICIT NONE
1099 
1100  !-----------------------------------------------------------------------------
1101  REAL, INTENT(IN) :: d_hs
1102  REAL, INTENT(IN) :: z3t_st
1103  REAL, INTENT(IN) :: rho
1104  REAL, INTENT(OUT) :: ycorr
1105  REAL, INTENT(OUT) :: dlnydr
1106 
1107  !-----------------------------------------------------------------------------
1108  INTEGER :: i, j
1109  REAL :: z0, z1, z2, z3, zms
1110  REAL, DIMENSION(nc,nc) :: dij_ab, gij, dgijdz
1111  !-----------------------------------------------------------------------------
1112 
1113  ncomp = 1
1114 
1115  !z0 = PI / 6.0 * rho * SUM( x(1:ncomp) * mseg(1:ncomp) )
1116  !z1 = PI / 6.0 * rho * SUM( x(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp) )
1117  !z2 = PI / 6.0 * rho * SUM( x(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**2 )
1118  !z3 = PI / 6.0 * rho * SUM( x(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
1119 
1120  z0 = pi / 6.0 * rho * parame(1,1)
1121  z1 = pi / 6.0 * rho * parame(1,1) * d_hs
1122  z2 = pi / 6.0 * rho * parame(1,1) * d_hs**2
1123  z3 = pi / 6.0 * rho * parame(1,1) * d_hs**3
1124  zms = 1.0 - z3
1125 
1126  ! DO i = 1,ncomp
1127  ! DO j=1,ncomp
1128  ! dij_ab(i,j)=dhs(i)*dhs(j)/(dhs(i)+dhs(j))
1129  ! ENDDO
1130  ! END DO
1131  dij_ab(1,1) = d_hs / 2.0
1132 
1133  DO i = 1, ncomp
1134  DO j = 1, ncomp
1135  gij(i,j) = 1.0/zms + 3.0*dij_ab(i,j)*z2/zms/zms &
1136  + 2.0*(dij_ab(i,j)*z2)**2 /zms**3
1137  dgijdz(i,j)= 1.0/zms/zms + 3.0*dij_ab(i,j)*z2*(1.0+z3)/z3/zms**3 &
1138  + (dij_ab(i,j)*z2/zms/zms)**2 *(4.0+2.0*z3)/z3
1139  END DO
1140  END DO
1141 
1142  ycorr = gij(1,1)
1143  dlnydr = dgijdz(1,1) / gij(1,1) * z3t_st
1144 
1145 END SUBROUTINE cavity
1146 
1147 
1148 
1149 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1150 ! Module DFT_INVERSION
1151 !
1152 ! This module contains variables for a LDA - Inversion Procedure (an algo-
1153 ! rithm for solving for the equilibrium density profile). Procedure is
1154 ! similar to R.Evans (Bristol).
1155 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1156 
1157 Module dft_inversion
1158 
1159  implicit none
1160  save
1161 
1162  INTEGER :: no_step_l, no_step_h
1163  REAL :: mu_eta_div, den_min, den_max, den_div
1164  REAL, DIMENSION(200) :: rho_array1, rho_array2
1165  REAL, DIMENSION(200) :: mu_rho_array1, mu_rho_array2
1166 
1167 End Module dft_inversion
1168 
1169 
1170 
1171 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1172 ! SUBROUTINE DENS_INV_COEFF2
1173 !
1174 ! First step of LDA - Inversion Procedure (algorithm for solving for the
1175 ! equilibrium density profile): Determine coefficients.
1176 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1177 
1178 SUBROUTINE dens_inv_coeff2 ( rhob )
1179 
1180  USE basic_variables, ONLY: np, nc, dense, densta, parame
1181  USE dft_module
1182  USE dft_inversion
1183  USE utilities
1184  IMPLICIT NONE
1185 
1186  !-----------------------------------------------------------------------------
1187  REAL, INTENT(IN) :: rhob(2)
1188 
1189  !-----------------------------------------------------------------------------
1190  INTEGER :: i
1191  REAL :: den_st
1192  REAL :: lnphi(np,nc)
1193  !-----------------------------------------------------------------------------
1194 
1195  den_min = rhob(2)*z3t_st * 0.8
1196  den_max = rhob(1)*z3t_st * 1.1
1197  den_div = 0.04
1198 
1199  no_step_l = 80
1200  no_step_h = 200
1201 
1202  mu_eta_div = 0.0
1203  IF ( den_min < den_div ) THEN
1204  DO i = 1, no_step_l
1205  dense(1) = den_min + (den_div-den_min) * REAL(i-1)/REAL(no_step_l-1)
1206  densta(1) = dense(1)
1207  rho_array1(i) = dense(1)
1208  CALL fugacity (lnphi)
1209  mu_rho_array1(i) = exp( lnphi(1,1)+log(dense(1)/z3t_st / parame(1,2)**3 ) ) ! myrho = mu_res(T,rho) + ln(rho)
1210  END DO
1211  mu_eta_div = mu_rho_array1(no_step_l)
1212  END IF
1213 
1214  den_st = max( den_min, den_div )
1215  DO i = 1, no_step_h
1216  dense(1) = den_st + (den_max-den_st) * REAL(i-1)/REAL(no_step_h-1)
1217  densta(1) = dense(1)
1218  rho_array2(i) = dense(1)
1219  CALL fugacity (lnphi)
1220  mu_rho_array2(i) = lnphi(1,1) + log(dense(1)/z3t_st /parame(1,2)**3 ) ! myrho = mu_res(T,rho) + ln(rho)
1221  IF ( i >= 2 ) THEN
1222  IF ( (mu_rho_array2(i)- mu_rho_array2(i-1)) < 0.0 ) THEN
1223  call paus ('DENS_INV_COEFF2: derivative negative')
1224  END IF
1225  END IF
1226  END DO
1227 
1228 END SUBROUTINE dens_inv_coeff2
1229 
1230 
1231 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1232 ! FUNCTION DENS_INVERS2
1233 !
1234 ! Second step of LDA - Inversion Procedure (algorithm for solving for the
1235 ! equilibrium density profile): Determine new guess of density profile.
1236 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1237 
1238 REAL FUNCTION dens_invers2 ( rhob, mu_rho )
1239 
1240  USE dft_module
1241  USE dft_inversion
1242  IMPLICIT NONE
1243 
1244  !-----------------------------------------------------------------------------
1245  REAL, INTENT(IN) :: rhob(2)
1246  REAL, INTENT(IN OUT) :: mu_rho
1247 
1248  !-----------------------------------------------------------------------------
1249  INTEGER :: i
1250  REAL :: eta_i
1251  LOGICAL :: flag_exit
1252  !-----------------------------------------------------------------------------
1253 
1254  ! den_min = rhob(2)*z3t_st * 0.8
1255  ! den_max = rhob(1)*z3t_st * 1.1
1256  ! den_div = 0.04
1257 
1258  ! no_step_l = 40
1259  ! no_step_h = 100
1260 
1261  IF ( exp(mu_rho) < mu_eta_div ) THEN
1262  i = 2
1263  flag_exit = .false.
1264  DO WHILE (i < no_step_l .AND. .NOT. flag_exit)
1265  i = i + 1
1266  IF ( mu_rho_array1(i) > exp(mu_rho) ) THEN
1267  eta_i = rho_array1(i-1)+(rho_array1(i)-rho_array1(i-1)) &
1268  /(mu_rho_array1(i)-mu_rho_array1(i-1)) &
1269  *(exp(mu_rho)-mu_rho_array1(i-1))
1270  dens_invers2 = eta_i/z3t_st
1271  flag_exit = .true.
1272  END IF
1273  END DO
1274  IF ( .NOT. flag_exit ) WRITE (*,*) 'error 1', exp(mu_rho), mu_eta_div
1275  ! IF ( .NOT. flag_exit ) call paus (' ')
1276  ENDIF
1277 
1278  IF ( exp(mu_rho) >= mu_eta_div ) THEN
1279  i = 2
1280  flag_exit = .false.
1281  DO WHILE (i < no_step_h .AND. .NOT. flag_exit )
1282  i = i + 1
1283  IF ( mu_rho_array2(i) > mu_rho ) THEN
1284  eta_i = rho_array2(i-1)+(rho_array2(i)-rho_array2(i-1)) &
1285  /(mu_rho_array2(i)-mu_rho_array2(i-1)) *(mu_rho-mu_rho_array2(i-1))
1286  dens_invers2 = eta_i/z3t_st
1287  flag_exit = .true.
1288  END IF
1289  END DO
1290  IF ( .NOT. flag_exit ) WRITE (*,*) 'error 2', mu_rho, mu_eta_div
1291  ! IF ( .NOT. flag_exit ) call paus (' ')
1292  END IF
1293 
1294  IF ( dens_invers2 < (rhob(2)*0.8) ) THEN
1295  WRITE (*,*) 'lower bound', dens_invers2, mu_rho
1296  dens_invers2 = rhob(2) * 0.8
1297  END IF
1298  IF ( dens_invers2 > (rhob(1)*1.1) ) THEN
1299  WRITE (*,*) 'upper bound', dens_invers2, mu_rho
1300  dens_invers2 = rhob(1)
1301  END IF
1302 
1303 END FUNCTION dens_invers2
1304 
1305 
1306 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1307 ! SUBROUTINE spline
1308 !
1309 ! Given arrays x(1:n) and y(1:n) containing a tabulated function,
1310 ! i.e., yi = f(xi), with x1 < x2 < .. . < xN, and given values yp1 and
1311 ! ypn for the first derivative of the interpolating function at points 1
1312 ! and n, respectively, this routine returns an array y2(1:n) of length n
1313 ! which contains the second derivatives of the interpolating function at
1314 ! the tabulated points xi. If yp1 and/or ypn are equal to 1 1030 or
1315 ! larger, the routine is signaled to set the corresponding boundary
1316 ! condition for a natural spline, with zero second derivative on that
1317 ! boundary.
1318 ! Parameter: NMAX is the largest anticipated value of n.
1319 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1320 
1321 SUBROUTINE spline ( x, y, n, yp1, ypn, y2 )
1322 
1323  IMPLICIT NONE
1324 
1325  !-----------------------------------------------------------------------------
1326  INTEGER, INTENT(IN) :: n
1327  REAL, INTENT(IN) :: x(n)
1328  REAL, INTENT(IN) :: y(n)
1329  REAL, INTENT(IN) :: yp1
1330  REAL, INTENT(IN) :: ypn
1331  REAL, INTENT(OUT) :: y2(n)
1332 
1333  !-----------------------------------------------------------------------------
1334  INTEGER, PARAMETER :: nmax = 500
1335  INTEGER :: i, k
1336  REAL :: p, qn, sig, un, u(nmax)
1337  !-----------------------------------------------------------------------------
1338 
1339  IF ( yp1 > 0.99e30 ) THEN
1340  y2(1) = 0.0
1341  u(1) = 0.0
1342  ELSE
1343  y2(1) = -0.5
1344  u(1) = ( 3.0/(x(2)-x(1)) ) * ( (y(2)-y(1))/(x(2)-x(1))-yp1 )
1345  END IF
1346  DO i = 2, n-1
1347  IF ( (x(i+1)-x(i)) == 0.0 .OR. (x(i)-x(i-1)) == 0.0 .OR. (x(i+1)-x(i-1)) == 0.0 ) THEN
1348  write (*,*) 'error in spline-interpolation'
1349  stop
1350  END IF
1351  sig = (x(i)-x(i-1)) / (x(i+1)-x(i-1))
1352  p = sig*y2(i-1) + 2.0
1353  y2(i) = (sig-1.0) / p
1354  u(i) = ( 6.0 * ((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1))) / (x(i+1)-x(i-1)) &
1355  - sig * u(i-1) ) / p
1356  END DO
1357  IF ( ypn > 0.99e30 ) THEN
1358  qn = 0.0
1359  un = 0.0
1360  ELSE
1361  qn = 0.5
1362  un = (3.0/(x(n)-x(n-1))) * (ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
1363  END IF
1364  y2(n) = (un-qn*u(n-1)) / (qn*y2(n-1)+1.0)
1365  DO k = n-1, 1, -1
1366  y2(k) = y2(k) * y2(k+1) + u(k)
1367  END DO
1368 
1369 END SUBROUTINE spline
1370 
1371 
1372 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1373 ! SUBROUTINE splint(xa,ya,y2a,n,x,y)
1374 !
1375 ! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a
1376 ! function (with the in order), and given the array y2a(1:n), which is
1377 ! the output from spline above, and given a value of x, this routine
1378 ! returns a cubic-spline interpolated value y.
1379 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1380 
1381 SUBROUTINE splint ( xa, ya, y2a, n, x, y )
1382 
1383  USE utilities
1384  IMPLICIT NONE
1385 
1386  !-----------------------------------------------------------------------------
1387  INTEGER, INTENT(IN) :: n
1388  REAL, INTENT(IN) :: xa(n)
1389  REAL, INTENT(IN) :: ya(n)
1390  REAL, INTENT(IN) :: y2a(n)
1391  REAL, INTENT(IN OUT) :: x
1392  REAL, INTENT(OUT) :: y
1393 
1394  !-----------------------------------------------------------------------------
1395  INTEGER :: k, khi, klo
1396  REAL :: a, b, h
1397  !-----------------------------------------------------------------------------
1398 
1399  klo = 1
1400  khi = n
1401 1 IF ( khi-klo > 1 ) THEN
1402  k = ( khi + klo ) / 2
1403  IF ( xa(k) > x ) THEN
1404  khi = k
1405  ELSE
1406  klo = k
1407  END IF
1408  GO TO 1
1409  END IF
1410  h = xa(khi) - xa(klo)
1411  IF ( h == 0.0 ) call paus ('bad xa input in splint')
1412  a = ( xa(khi) - x ) / h
1413  b = ( x - xa(klo) ) / h
1414  y = a*ya(klo) + b*ya(khi) + ( (a**3-a)*y2a(klo)+(b**3-b)*y2a(khi) )*h**2 / 6.0
1415 
1416 END SUBROUTINE splint
1417 
1418 
1419 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1420 ! SUBROUTINE splint_integral
1421 !
1422 ! Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a
1423 ! function (with the in order), and given the array y2a(1:n), which is
1424 ! the output from spline above, and given a value of x, this routine
1425 ! returns a cubic-spline interpolated value y.
1426 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1427 
1428 SUBROUTINE splint_integral ( xa, ya, y2a, n, xlo, xhi, integral )
1429 
1430  USE utilities
1431  IMPLICIT NONE
1432 
1433  !-----------------------------------------------------------------------------
1434  INTEGER, INTENT(IN) :: n
1435  REAL, INTENT(IN) :: xa(n)
1436  REAL, INTENT(IN) :: ya(n)
1437  REAL, INTENT(IN) :: y2a(n)
1438  REAL, INTENT(IN) :: xlo
1439  REAL, INTENT(IN) :: xhi
1440  REAL, INTENT(OUT) :: integral
1441 
1442  !-----------------------------------------------------------------------------
1443  INTEGER :: k, khi_l, klo_l, khi_h, klo_h
1444  REAL :: xl, xh = 0.0
1445  REAL :: h, int, x0, x1, y0, y1, y20, y21
1446  !-----------------------------------------------------------------------------
1447 
1448  integral = 0.0
1449  klo_l = 1
1450  khi_l = n
1451  do while ( khi_l-klo_l > 1 )
1452  k = ( khi_l + klo_l ) / 2
1453  IF ( xa(k) > xlo ) THEN
1454  khi_l = k
1455  ELSE
1456  klo_l = k
1457  END IF
1458 
1459  end do
1460 
1461  klo_h = 1
1462  khi_h = n
1463  do while ( khi_h-klo_h > 1 )
1464  k = ( khi_h + klo_h ) / 2
1465  IF ( xa(k) > xhi ) THEN
1466  khi_h = k
1467  ELSE
1468  klo_h = k
1469  END IF
1470  end do
1471 
1472  ! integration in spline pieces, the lower interval, bracketed
1473  ! by xa(klo_L) and xa(khi_L) is in steps shifted upward.
1474 
1475  ! first: determine upper integration bound
1476  xl = xlo
1477 3 CONTINUE
1478  IF ( khi_h > khi_l ) THEN
1479  xh = xa(khi_l)
1480  ELSE IF ( khi_h == khi_l ) THEN
1481  xh = xhi
1482  ELSE
1483  call paus ('error in spline-integration')
1484  END IF
1485 
1486  h = xa(khi_l) - xa(klo_l)
1487  IF ( h == 0.0 ) call paus ('bad xa input in splint')
1488  x0 = xa(klo_l)
1489  x1 = xa(khi_l)
1490  y0 = ya(klo_l)
1491  y1 = ya(khi_l)
1492  y20= y2a(klo_l)
1493  y21= y2a(khi_l)
1494  ! int = -xL/h * ( (x1-.5*xL)*y0 + (0.5*xL-x0)*y1 &
1495  ! +y20/6.*(x1**3-1.5*xL*x1*x1+xL*xL*x1-.25*xL**3) &
1496  ! -y20/6.*h*h*(x1-.5*xL) &
1497  ! +y21/6.*(.25*xL**3-xL*xL*x0+1.5*xL*x0*x0-x0**3) &
1498  ! -y21/6.*h*h*(.5*xL-x0) )
1499  ! int = int + xH/h * ( (x1-.5*xH)*y0 + (0.5*xH-x0)*y1 &
1500  ! +y20/6.*(x1**3-1.5*xH*x1*x1+xH*xH*x1-.25*xH**3) &
1501  ! -y20/6.*h*h*(x1-.5*xH) &
1502  ! +y21/6.*(.25*xH**3-xH*xH*x0+1.5*xH*x0*x0-x0**3) &
1503  ! -y21/6.*h*h*(.5*xH-x0) )
1504  int = -1.0/h * ( xl*((x1-.5*xl)*y0 + (0.5*xl-x0)*y1) &
1505  -y20/24.*(x1-xl)**4 + y20/6.*(0.5*xl*xl-x1*xl)*h*h &
1506  +y21/24.*(xl-x0)**4 - y21/6.*(0.5*xl*xl-x0*xl)*h*h )
1507  int = int + 1.0/h * ( xh*((x1-.5*xh)*y0 + (0.5*xh-x0)*y1) &
1508  -y20/24.*(x1-xh)**4 + y20/6.*(0.5*xh*xh-x1*xh)*h*h &
1509  +y21/24.*(xh-x0)**4 - y21/6.*(0.5*xh*xh-x0*xh)*h*h )
1510 
1511  integral = integral + int
1512  ! write (*,*) integral,x1,xH
1513  klo_l = klo_l + 1
1514  khi_l = khi_l + 1
1515  xl = x1
1516  IF (khi_h /= (khi_l-1)) GO TO 3 ! the -1 in (khi_L-1) because khi_L was already counted up
1517 
1518 END SUBROUTINE splint_integral
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
lambda
double lambda
Latent heat of blowing agent, J/kg.
Definition: experimentalInputs.h:118
sig
double sig
correlated to variance of initial distribution
Definition: experimentalInputs.h:138
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
eos_numerical_derivatives
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:220