MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
phase_equilib.f90
1 
2 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
3 ! SUBROUTINE phase_equilib
4 !
5 ! This subroutine varies a predefined "running variable" and
6 ! organizes phase equilibrium calculations. For an isotherm
7 ! calculation e.g., the running variable is often the pressure. The
8 ! code is designed to deliver only converged solutions. In order to
9 ! enforce convergence, a step-width adjustment (reduction) is
10 ! implemented.
11 
12 ! VARIABLE LIST:
13 ! running defines the running variable. For example: if you want
14 ! to calculate the vapor pressure curve of a component
15 ! starting from 100C to 200C, then running is 't'. The
16 ! temperature is step-wise increased until the end-
17 ! -temperature of 200C is reached.
18 ! (in this example end_x=200+273.15)
19 ! end_x end point for running variable
20 ! steps No. of calculation steps towards the end point of calc.
21 ! converg 0 if no convergence achieved, 1 if converged solution
22 !
23 ! PREREQUISITES:
24 ! prior to execution of this routine, the follwing variables have to
25 ! be defined: "val_init" an array containing the starting values for
26 ! this iteration, "it(i)" provides the information, which variable is
27 ! determined iteratively, "sum_rel(i)" indicates, which mole fraction
28 ! is determined from the summation relation sum(xi)=1. Furthermore,
29 ! the number of phases and the variables provided by the subroutine
30 ! INPUT are required.
31 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
32 
33 SUBROUTINE phase_equilib (end_x,steps,converg)
34 
35  USE basic_variables
36  IMPLICIT NONE
37 
38  !-----------------------------------------------------------------------------
39  REAL, INTENT(IN) :: end_x
40  REAL, INTENT(IN) :: steps
41  INTEGER, INTENT(OUT) :: converg
42 
43  !-----------------------------------------------------------------------------
44  INTEGER :: k, count1,count2,runindex,maxiter
45  REAL :: delta_x,delta_org,val_org,runvar
46  CHARACTER (LEN=2) :: compon
47  LOGICAL :: continue_cycle
48  !-----------------------------------------------------------------------------
49 
50  !-----------------------------------------------------------------------------
51  ! determine variable that is varied ("running variable")
52  !-----------------------------------------------------------------------------
53 
54  IF (running(1:2) == 'd1') runindex = 1
55  IF (running(1:2) == 'd2') runindex = 2
56  IF (running(1:1) == 't') runindex = 3
57  IF (running(1:1) == 'p') runindex = 4
58  IF (running(1:2) == 'x1') compon = running(3:3)
59  IF (running(1:2) == 'x1') READ (compon,*) k
60  IF (running(1:2) == 'x1') runindex = 4+k
61  IF (running(1:2) == 'x2') compon = running(3:3)
62  IF (running(1:2) == 'x2') READ (compon,*) k
63  IF (running(1:2) == 'x2') runindex = 4+ncomp+k
64  IF (running(1:2) == 'l1') compon = running(3:3)
65  IF (running(1:2) == 'l1') READ (compon,*) k
66  IF (running(1:2) == 'l1') runindex = 4+k
67  IF (running(1:2) == 'l2') compon = running(3:3)
68  IF (running(1:2) == 'l2') READ (compon,*) k
69  IF (running(1:2) == 'l2') runindex = 4+ncomp+k
70 
71  !-----------------------------------------------------------------------------
72  ! define control variables of the running variable
73  !-----------------------------------------------------------------------------
74 
75  maxiter = 200
76  IF ( ncomp >= 3 ) maxiter = 1000
77  count1 = 0
78  count2 = 0
79  delta_x = ( end_x - val_init(runindex) ) / steps
80  delta_org = ( end_x - val_init(runindex) ) / steps
81  val_org = val_init(runindex)
82  IF ( running(1:1) == 'x' ) THEN
83  delta_x = ( end_x - exp(val_init(runindex)) ) / steps
84  delta_org = ( end_x - exp(val_init(runindex)) ) / steps
85  val_org = exp(val_init(runindex))
86  END IF
87 
88  continue_cycle = .true.
89 
90  !-----------------------------------------------------------------------------
91  ! stepwise change of the running variables
92  !-----------------------------------------------------------------------------
93 
94  DO WHILE ( continue_cycle )
95 
96  count1 = count1 + 1
97  count2 = count2 + 1
98  ! val_org = val_init(runindex)
99 
100 
101  CALL objective_ctrl (converg)
102 
103  IF (converg == 1) THEN
104  val_init( 1:(4+ncomp*nphas) ) = val_conv( 1:(4+ncomp*nphas) )
105  IF (outp >= 1 .AND. (abs(delta_x) > 0.1*abs(delta_org) .OR. count2 == 2)) CALL output
106  ! write (17,'(4G18.8)') rhoi_cal(1,1:2)*32000.,rhoi_cal(2,1:2)*32000.
107  ELSE
108  delta_x = delta_x / 2.0
109  IF (num == 2) delta_x = delta_x / 2.0
110  val_init(runindex) = val_org
111  IF (running(1:1) == 'x') val_init(runindex) = log(val_org)
112  continue_cycle = .true.
113  count2 = 0
114  END IF
115  runvar = val_init(runindex)
116  IF (running(1:1) == 'x') runvar = exp(val_init(runindex))
117 
118  IF ( end_x == 0.0 .AND. running(1:1) /= 'x' ) THEN
119  IF ( abs(runvar-end_x) < 1.e-8 ) continue_cycle = .false.
120  ELSE IF ( abs((runvar-end_x)/end_x) < 1.e-8 ) THEN
121  ! IF(delta_org.NE.0.0) WRITE (*,*)' FINISHED ITERATION',count1
122  continue_cycle = .false.
123  ELSE IF ( count1 == maxiter ) THEN
124  WRITE (*,*) ' MAX. NO OF ITERATIONS',count1
125  converg = 0
126  continue_cycle = .false.
127  ELSE IF ( abs(delta_x) < 1.e-5*abs(delta_org) ) THEN
128  ! WRITE (*,*) ' CLOSEST APPROACH REACHED',count1
129  converg = 0
130  continue_cycle = .false.
131  ELSE
132  continue_cycle = .true.
133  val_org = runvar
134  IF (abs(runvar+delta_x-end_x) > abs(runvar-end_x)) delta_x = end_x - runvar ! if end-point passed
135  val_init(runindex) = runvar + delta_x
136  IF (running(1:1) == 'x') val_init(runindex) = log(runvar+delta_x)
137  END IF
138 
139  IF (abs(delta_x) < abs(delta_org) .AND. count2 >= 5) THEN
140  delta_x = delta_x * 2.0
141  count2 = 0
142  END IF
143 
144  END DO ! continue_cycle
145 
146 END SUBROUTINE phase_equilib
147 
148 
149 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
150 ! SUBROUTINE objective_ctrl
151 !
152 ! This subroutine controls the iso-fugacity iteration. It uses
153 ! the variables defined in the array "val_init". If successfull,
154 ! the converged values are written to "val_conv", and the flag
155 ! converg is set to 1.
156 ! See also above desciption for subroutine PHASE_EQUILIB
157 ! This routine calls SUBROUTINE HYBRID, which is a solver (modified
158 ! POWELL HYBRID METHOD). HYBRID is freely available for non-commercial
159 ! applications. HYBRID requires three definitions:
160 ! 1.the number of equations to be solved (=No. of variables to be
161 ! iterated). The appropriate parameter is: "n_unkw"
162 ! 2.the equations to be iterated, they are here gathered in the SUB-
163 ! ROUTINE OBJEC_FCT (see below). Since HYBRID is a root finder,
164 ! these objective functions are iterated to be zero (essentially,
165 ! OBJEC_FCT contains the iso-fugacity relation.
166 ! 3.an array of variables is required, containing the quatities to be
167 ! iterated. This array is termed "y(i)"
168 !
169 ! INPUT VARIABLES:
170 ! val_init(i) array containing (densities,T,P,lnx's) serving as
171 ! starting values for the phase equilibrium calculation
172 ! it(i) contains the information, which variable is deter-
173 ! mined iteratively. For syntax refer e.g.to SUB BINMIX.
174 ! sum_rel(i) indicates, which mole fraction is determined from the
175 ! summation relation sum(xi)=1
176 !
177 ! OUTPUT VARIABLES:
178 ! val_conv(i) array containing the converged system variables
179 ! analogous to "val_init"
180 ! converg 0 if no convergence achieved, 1 if converged solution
181 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
182 
183 SUBROUTINE objective_ctrl (converg)
184 
185  USE basic_variables
186  USE starting_values
187  USE eos_variables, only: pi, kbol, mseg, dhs
188  USE solve_nonlin
189  USE utilities
190  IMPLICIT NONE
191 
192  !-----------------------------------------------------------------------------
193  INTEGER, INTENT(OUT) :: converg
194 
195  !-----------------------------------------------------------------------------
196  INTERFACE
197  SUBROUTINE objec_fct ( iter_no, y, residu, dummy )
198  INTEGER, INTENT(IN) :: iter_no
199  REAL, INTENT(IN) :: y(iter_no)
200  REAL, INTENT(OUT) :: residu(iter_no)
201  INTEGER, INTENT(IN OUT) :: dummy
202  END SUBROUTINE objec_fct
203  END INTERFACE
204 
205  INTEGER :: info, k, posn, i
206  REAL, ALLOCATABLE :: y(:), diag(:), residu(:)
207  REAL :: x_init, x_solut, r_diff1, r_diff2, totres
208  REAL :: r_thrash, x_thrash
209  CHARACTER (LEN=2) :: compon
210 
211  INTEGER :: ph_split, promising_min
212  REAL :: eta_trial, phi_2_start
213  REAL, dimension(nc) :: xif_save, rhoi_feed, rhoi1, rhoi2
214  LOGICAL :: convergence
215  !-----------------------------------------------------------------------------
216 
217  info=1
218 
219  ALLOCATE( y(n_unkw), diag(n_unkw), residu(n_unkw) )
220 
221  IF (num == 0) acc_a = 1.e-7
222  IF (num == 0) step_a = 2.e-8
223  IF (num == 1) acc_a = 1.e-7
224  IF (num == 1) step_a = 2.e-8
225  IF (num == 2) acc_a = 5.e-7
226  IF (num == 2) step_a = 1.e-7
227 
228  !-----------------------------------------------------------------------------
229  ! assign starting values for the iterated variables (to vector y)
230  !-----------------------------------------------------------------------------
231 
232  posn = 0
233  DO i = 1,n_unkw
234  posn = posn + 1
235  IF (it(i) == 't') y(posn) = val_init(3)
236  IF (it(i) == 'p') y(posn) = val_init(4)
237  IF (it(i) == 'lnp') y(posn) = log( val_init(4) )
238  IF (it(i) == 'fls') y(posn) = alpha
239  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '1') compon = it(i)(3:3)
240  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '1') READ (compon,*) k
241  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '1') y(posn) = val_init(4+k)
242  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '2') compon = it(i)(3:3)
243  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '2') READ (compon,*) k
244  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '2') y(posn) = val_init(4+ncomp+k)
245  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '3') compon = it(i)(3:3)
246  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '3') READ (compon,*) k
247  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '3') y(posn) = val_init(4+ncomp+ncomp+k)
248  END DO
249 
250  !-----------------------------------------------------------------------------
251  ! initialize some variables
252  !-----------------------------------------------------------------------------
253 
254  CALL init_vars
255 
256  x_init = 0.0
257  DO i = 1,ncomp
258  IF (lnx(1,i) /= 0.0 .AND. lnx(2,i) /= 0.0) THEN
259  x_init = x_init + abs( 1.0 - lnx(2,i)/lnx(1,i) )
260  ELSE
261  x_init = x_init + abs( 1.0 - exp(lnx(2,i))/exp(lnx(1,i)) )
262  END IF
263  END DO
264 
265  !-----------------------------------------------------------------------------
266  ! solve the phase equilibrium conditions
267  !-----------------------------------------------------------------------------
268 
269  CALL hbrd (objec_fct, n_unkw, y, residu, step_a, acc_a, info, diag)
270  IF ( info <= 1 .AND. sum( abs( residu(1:n_unkw) ) ) > 100.0*acc_a ) &
271  CALL hbrd (objec_fct, n_unkw, y, residu, step_a, acc_a, info, diag)
272 
273  !-----------------------------------------------------------------------------
274  ! determine, whether the solution is acceptable (if yes: convergence = .true.)
275  !-----------------------------------------------------------------------------
276 
277  x_solut = 0.0
278  DO i = 1,ncomp
279  IF ( lnx(1,i) /= 0.0 .AND. lnx(2,i) /= 0.0 ) THEN
280  x_solut = x_solut + abs( 1.0 - lnx(2,i)/lnx(1,i) )
281  ELSE
282  IF (lnx(1,i) < 300.0 .AND. lnx(1,i) > -300.0 ) &
283  x_solut = x_solut + abs( 1.0 - exp(lnx(2,i))/exp(lnx(1,i)) )
284  END IF
285  END DO
286  r_diff1 = abs( 1.0 - dense(1)/dense(2) )
287  IF ( val_conv(2) > 0.0 ) THEN
288  r_diff2 = abs( 1.0 - val_conv(1)/val_conv(2) )
289  ELSE
290  r_diff2 = 0.0
291  END IF
292 
293  totres = sum( abs( residu(1:n_unkw) ) )
294 
295  r_thrash = 0.0005
296  x_thrash = 0.0005
297  if (num > 0 ) r_thrash = r_thrash * 10.0
298  if (num > 0 ) x_thrash = x_thrash * 100.0
299 
300  convergence = .true.
301 
302  IF ( info >= 2 ) convergence = .false.
303  IF ( abs( 1.0- dense(1)/dense(2) ) < r_thrash .AND. x_solut < x_thrash ) THEN
304  IF ( x_init > 0.050 ) convergence = .false.
305  IF ( ( abs( 1.0- dense(1)/dense(2) ) + x_solut ) < 1.e-7 ) convergence = .false.
306  ENDIF
307  IF ( r_diff2 /= 0.0 .AND. r_diff2 > (4.0*r_diff1) .AND. bindiag == 1 ) convergence = .false.
308  IF ( ncomp == 1 .AND. totres > 100.0*acc_a ) convergence = .false.
309  IF ( ncomp == 1 .AND. r_diff1 < 1.e-5 ) convergence = .false.
310  IF ( totres > 1000.0*acc_a .AND. .NOT. generate_starting_val ) convergence = .false.
311  IF ( totres > 10000.0*acc_a .AND. generate_starting_val ) convergence = .false.
312 
313  !-----------------------------------------------------------------------------
314  ! if convergence: test stability, save result to vector (val_conv) and output result
315  !-----------------------------------------------------------------------------
316 
317  IF ( convergence ) THEN
318 
319  !-----------------------------------------------------------------------------
320  ! store result and output result to terminal
321  !-----------------------------------------------------------------------------
322 
323  converg = 1
324  CALL converged
325  CALL enthalpy_etc
326  ! IF (outp == 2 ) CALL output
327 
328  !-----------------------------------------------------------------------------
329  ! test phase stability
330  !-----------------------------------------------------------------------------
331 
332  if ( check_stability_of_phases ) then
333  my_f( 1:ncomp ) = my_cal(1,1:ncomp)
334  rhoi_feed( 1:ncomp ) = rhoi_cal( maxloc(dense(1:2), 1), 1:ncomp )
335  xif_save(1:ncomp) = xif(1:ncomp)
336  xif(1:ncomp) = 0.95 * rhoi_cal( maxloc(dense(1:2), 1), 1:ncomp ) &
337  / sum( rhoi_cal( maxloc(dense(1:2), 1), 1:ncomp ) ) &
338  + 0.05 * rhoi_cal( minloc(dense(1:2), 1), 1:ncomp ) &
339  / sum( rhoi_cal( minloc(dense(1:2), 1), 1:ncomp ) )
340  if ( minval( dense(1:2) ) < 0.18 ) eta_trial = 0.40
341  if ( minval( dense(1:2) ) >=0.18 ) eta_trial = p / ( kbol*1.e30 * t ) &
342  * pi/6.0 * sum( rhoi_feed( 1:ncomp ) / sum(rhoi_feed( 1:ncomp )) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
343  call phase_stability ( rhoi_feed, eta_trial, ph_split, rhoi_trial )
344 
345  if ( ph_split == 1 ) then
346  rhoi1( 1:ncomp ) = rhoi_trial( 1:ncomp )
347  rhoi2( 1:ncomp ) = rhoi_feed( 1:ncomp )
348  call tangent_plane_line_search ( rhoi1, rhoi2, phi_2_start, promising_min )
349 
350  dense( minloc(dense(1:2), 1) ) = pi/6.0 *sum( rhoi_trial(1:ncomp) * mseg(1:ncomp) * dhs(1:ncomp)**3 )
351  call rachford_rice ( converg, rhoi1, rhoi2 )
352  if ( converg == 1 ) then
353  call occupy_val_init ( rhoi1, rhoi2 )
354  if ( ncomp == 2 ) CALL select_sum_rel (1,0,1)
355  if ( ncomp == 2 ) CALL select_sum_rel (2,0,2)
356  val_conv = val_init
357  if ( outp > 0 ) write (*,*) ' '
358  if ( outp > 0 ) write (*,*) '===== the type (LLE, VLE) of phase equilibrium changed ====='
359  if ( outp > 0 ) write (*,*) ' '
360  else
361  call restore_converged
362  converg = 1
363  if ( outp > 0 .AND. ph_split == 1 ) write (*,*) 'calc. phase equilibr. not stable (press enter)'
364  if ( outp > 0 .AND. ph_split == 1 ) read (*,*)
365  end if
366  else
367  call restore_converged
368  converg = 1
369  end if
370  xif(1:ncomp) = xif_save(1:ncomp)
371  end if
372 
373  ELSE
374 
375  converg = 0
376 
377  END IF
378 
379  DEALLOCATE( y, diag, residu )
380 
381 END SUBROUTINE objective_ctrl
382 
383 
384 
385 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
386 ! SUBROUTINE objec_fct
387 !
388 ! This subroutine contains the equations to be solved numerically
389 ! (iso-fugacity: fi'-fi''=0) as well as other dependent equations,
390 ! which can be solved analytically, namely the summation relation
391 ! xi=1-sum(xj) or the condition of equal charge for electrolyte
392 ! solutions.
393 ! This subroutine is required and controlled by the solver HBRD !
394 ! HBRD varies the variables "y(i)" and eveluates the result of
395 ! these changes from this routine.
396 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
397 
398 SUBROUTINE objec_fct ( iter_no, y, residu, dummy )
399 
400  USE basic_variables
401  USE eos_variables, ONLY: density_error
402  use utilities, only: x_summation
403  IMPLICIT NONE
404 
405  !-----------------------------------------------------------------------------
406  INTEGER, INTENT(IN) :: iter_no
407  REAL, INTENT(IN) :: y(iter_no)
408  REAL, INTENT(OUT) :: residu(iter_no)
409  INTEGER, INTENT(IN OUT) :: dummy
410 
411  !-----------------------------------------------------------------------------
412  INTEGER :: i, ph,k,posn, skip,phase
413  REAL :: lnphi(np,nc),isofugacity
414  CHARACTER (LEN=2) :: compon
415  !-----------------------------------------------------------------------------
416 
417  !-----------------------------------------------------------------------------
418  ! assign values for the iterated quantities
419  !-----------------------------------------------------------------------------
420 
421  posn = 0
422  DO i = 1,n_unkw
423  posn = posn + 1
424  IF (it(i) == 't') t = y(posn)
425  IF (it(i) == 'p') p = y(posn)
426  IF (it(i) == 'lnp') p = exp( y(posn) )
427  IF (it(i) == 'fls') alpha = y(posn)
428  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '1') compon = it(i)(3:3)
429  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '1') READ (compon,*) k
430  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '1') lnx(1,k) = y(posn)
431  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '2') compon = it(i)(3:3)
432  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '2') READ (compon,*) k
433  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '2') lnx(2,k) = y(posn)
434  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '3') compon = it(i)(3:3)
435  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '3') READ (compon,*) k
436  IF (it(i)(1:1) == 'x' .AND. it(i)(2:2) == '3') lnx(3,k) = y(posn)
437  END DO
438 
439  DO k = 1,ncomp
440  IF (lnx(1,k) > 0.0) lnx(1,k) = 0.0
441  IF (lnx(2,k) > 0.0) lnx(2,k) = 0.0
442  END DO
443 
444  !-----------------------------------------------------------------------------
445  ! saveguard real values
446  !-----------------------------------------------------------------------------
447 
448  IF (p < 1.e-100) p = 1.e-12
449  IF ( p /= p ) p = 1000.0 ! rebounce for the case of NaN-solver output
450  IF ( t /= t ) t = 300.0 ! rebounce for the case of NaN-solver output
451  IF ( alpha /= alpha ) alpha = 0.5 ! rebounce for the case of NaN-solver output
452 
453  !-----------------------------------------------------------------------------
454  ! setting of mole fractions
455  !-----------------------------------------------------------------------------
456 
457  DO ph = 1, nphas
458  DO i = 1, ncomp
459  IF ( lnx(ph,i) < -300.0 ) THEN
460  xi(ph,i) = 0.0
461  ELSE
462  xi(ph,i) = exp( lnx(ph,i) )
463  END IF
464  END DO
465  END DO
466 
467  !-----------------------------------------------------------------------------
468  ! use summation relation to determine missing mole fractions
469  !-----------------------------------------------------------------------------
470 
471  IF (ncomp > 1) CALL x_summation
472 
473  !-----------------------------------------------------------------------------
474  ! calculate chemical potentials
475  !-----------------------------------------------------------------------------
476 
477  CALL fugacity (lnphi)
478 
479  !-----------------------------------------------------------------------------
480  ! determine the residuum to isofugacity relation
481  !-----------------------------------------------------------------------------
482 
483  phase = 2
484  DO i = 1,n_unkw
485  skip = 0 !for ions/polymers, the isofug-eq. is not always solved
486  IF (n_unkw < (ncomp*(nphas-1))) skip = ncomp*(nphas-1) - n_unkw
487  IF ((i+skip-ncomp*(phase-2)) > ncomp) phase = phase + 1
488  residu(i) = isofugacity((i+skip-ncomp*(phase-2)),phase,lnphi)
489  if ( density_error(phase) /= 0.0 ) residu(i) = residu(i) + sign( density_error(phase), residu(i) ) * 0.001
490  END DO
491 
492 END SUBROUTINE objec_fct
493 
494 
495 
496 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
497 ! REAL FUNCTION isofugacity
498 !
499 ! calculates the deviation from the condition of equal fugacities in
500 ! logarithmic form.
501 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
502 
503 REAL FUNCTION isofugacity (i,phase,lnphi)
504 
505  USE basic_variables
506  IMPLICIT NONE
507 
508  !-----------------------------------------------------------------------------
509  INTEGER, INTENT(IN) :: i
510  INTEGER, INTENT(IN) :: phase
511  REAL, INTENT(IN) :: lnphi(np,nc)
512 
513  !-----------------------------------------------------------------------------
514  INTEGER :: p1, p2
515  !-----------------------------------------------------------------------------
516 
517 
518  ! p1=1
519  p1 = phase-1
520  p2 = phase
521 
522  isofugacity = scaling(i) *( lnphi(p2,i)+lnx(p2,i)-lnx(p1,i)-lnphi(p1,i) )
523  ! write (*,'(a, 4G18.8)') ' t, p ',t,p,dense(1),dense(2)
524  ! write (*,'(a,i3,i3,3G24.14)') ' phi_V',i,p2,lnx(p2,i),lnphi(p2,i),dense(p2)
525  ! write (*,'(a,i3,i3,3G24.14)') ' phi_L',i,p1,lnx(p1,i),lnphi(p1,i),dense(p1)
526  ! write (*,*) ' isofugacity',i,isofugacity, scaling(i)
527  ! write (*,'(a,i3,4G18.8)') ' isofugacity',i,isofugacity, lnphi(p2,i)+lnx(p2,i), -lnx(p1,i)-lnphi(p1,i)
528  ! read (*,*)
529 
530 END FUNCTION isofugacity
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
starting_values
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW Module STARTING_VALUES This m...
Definition: modules.f90:251
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29