MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
pure_par_fit.f90
1 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
2 ! MODULE pure_fit_parameters
3 !
4 ! This module contains parameters and variables needed for the pure
5 ! component parameter regression.
6 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
7 
8 Module pure_fit_parameters
9 
10  implicit none
11  save
12 
13  !-----------------------------------------------------------------------------
14  INTEGER, PARAMETER :: cont = 4 ! number of data types (e.g. p^sat., rho^liq,...)
15  INTEGER, PARAMETER :: npmax = 40
16  INTEGER, PARAMETER :: mmx = 3500
17 
18  INTEGER :: nlv, nliq, n_vir, n_en
19  REAL, DIMENSION(npmax) :: rel
20  REAL, DIMENSION(cont) :: type_weight
21 
22  REAL, DIMENSION(MMX) :: plvcal, v_cal, b_cal, u_cal
23  REAL, DIMENSION(MMX) :: plv, tlv
24  REAL, DIMENSION(MMX) :: pliq, tliq, vliq
25  REAL, DIMENSION(MMX) :: b_vir, t_vir
26  REAL, DIMENSION(MMX) :: u_en, t_en, rho_en
27  REAL :: v_crit
28 
29  INTEGER, DIMENSION(10) :: fit_array
30  CHARACTER (LEN=25), DIMENSION(npmax) :: aa_txt
31  INTEGER, DIMENSION(cont) :: eqa ! flag, indicating, whether a data type is considered
32 
33  CHARACTER (LEN=50) :: pure_fit_file
34 
35 End Module pure_fit_parameters
36 
37 
38 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
39 ! SUBROUTINE PURE_RESIDUAL
40 !
41 ! Routine for pure component parameter fitting. The residual between
42 ! calculated values and experimental data is calculated and stored in
43 ! vector 'fvec'. The vector 'fvec' serves as the objective function for
44 ! the parameter.
45 ! PURE_RESIDUAL is called by FITTING via the Levenberg-Marquardt routine.
46 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
47 
48 SUBROUTINE pure_residual (lm_m_dat, lm_n_par, para, fvec, iflag)
49 
50  USE parameters, ONLY: rgas
51  USE basic_variables
52  USE pure_fit_parameters
53  USE utilities, only: si_dens, dens_calc, p_calc, error_message
54  IMPLICIT NONE
55 
56  !-----------------------------------------------------------------------------
57  ! INTEGER, PARAMETER :: dp = SELECTED_REAL_KIND(15, 307)
58  INTEGER, INTENT(IN) :: lm_m_dat
59  INTEGER, INTENT(IN) :: lm_n_par
60  REAL, INTENT(IN) :: para(:)
61  REAL, INTENT(IN OUT) :: fvec(:)
62  INTEGER, INTENT(IN OUT) :: iflag
63 
64  !-----------------------------------------------------------------------------
65  INTEGER :: converg, crit_dat
66  INTEGER :: i, j, ps
67  INTEGER :: j_dat = 0
68  LOGICAL :: scan,vapor,nocon
69  REAL :: weigh, tc, v_dum, plv_dum
70  REAL :: pges,zges
71  REAL :: density(np), w(np,nc)
72  REAL :: rho_phas(np)
73  CHARACTER (LEN=4) :: char_len
74  ! REAL :: tclit,pclit,pc,rhoc
75  !-----------------------------------------------------------------------------
76 
77  crit_dat = 0
78  tc = 0.0
79 
80  IF ( (nliq >= 1) .OR. (nlv >= 1) .OR. (n_vir >= 1) .OR. (n_en >= 1) ) THEN
81  DO i=1,11
82  ps = 0
83  DO j = 1, lm_n_par
84  IF (fit_array(j) == i) ps=j
85  END DO
86  IF (i == 1 .AND. ps >= 1) parame(1,1) = mm(1)*para(ps)*rel(ps)
87  IF (i == 2 .AND. ps >= 1) parame(1,2) = para(ps)*rel(ps)
88  IF (i == 3 .AND. ps >= 1) parame(1,3) = para(ps)*rel(ps)
89  IF (i == 5 .AND. ps >= 1) parame(1,13)= para(ps)*rel(ps)
90  IF (i == 6 .AND. ps >= 1) parame(1,6) = para(ps)*rel(ps)
91  IF (i == 7 .AND. ps >= 1) parame(1,8) = para(ps)*rel(ps)
92  IF (i == 8 .AND. ps >= 1) parame(1,7) = para(ps)*rel(ps)
93  IF (i == 9 .AND. ps >= 1) parame(1,9) = para(ps)*rel(ps)
94  IF (i == 11.AND. ps >= 1) parame(1,11)= para(ps)*rel(ps)
95  IF (i == 4 .AND. ps >= 1) THEN
96  IF ( nint(parame(1,12)) == 2 ) THEN
97  parame(1,14) = 0.0 ! eps_hb(1,1,1,1)
98  parame(1,15) = para(ps)*rel(ps) ! eps_hb(1,1,1,2)
99  parame(1,16) = para(ps)*rel(ps) ! eps_hb(1,1,2,1)
100  parame(1,17) = 0.0 ! eps_hb(1,1,2,2)
101  ELSE IF ( nint(parame(1,12)) == 1 ) THEN
102  parame(1,14) = para(ps)*rel(ps) ! eps_hb(1,1,1,1)
103  ELSE
104  call error_message('extend pure_par_fit for this association case')
105  END IF
106  END IF
107  END DO
108 
109  ! IF ( NINT(parame(1,12)) == 0) THEN ! no association
110  ! ELSEIF ( NINT(parame(1,12)) == 1) THEN ! No of assoc sites = 1
111  ! parame(1,15)=1.0
112  ! ELSEIF ( NINT(parame(1,12)) == 2) THEN ! No of assoc sites = 2
113  ! parame(1,18)=1.0
114  ! parame(1,19)=1.0
115  ! ELSEIF ( NINT(parame(1,12)) == 3) THEN ! No of assoc sites = 3
116  ! parame(1,18)=2.0
117  ! parame(1,19)=1.0
118  ! ELSEIF ( NINT(parame(1,12)) == 4) THEN ! No of assoc sites = 4
119  ! parame(1,18)=2.0
120  ! parame(1,19)=2.0
121  ! ELSE
122  ! write (*,*) ' define your association case'
123  ! stop
124  ! ENDIF
125 
126  parame(1,4) = 0.0
127  IF (eos == 0) parame(1,4) = 10.0
128 
129  END IF
130 
131 
132  WRITE (char_len,'(I3)') lm_n_par
133  IF (iflag == 1) WRITE (*,'(a,'//char_len//'(G15.8))') ' parameter ',(para(i)*rel(i),i=1,lm_n_par)
134 
135 
136  IF ( minval( parame(1,1:25) ) < 0.0 ) THEN
137  fvec(:) = fvec(:) * 2.0 ! penalty function for negative parameters
138  write (*,*) 'warning negative pure component parameters encountered',parame(1,minloc(parame(1,1:25)))
139  RETURN
140  ENDIF
141 
142 
143  fvec = 0.0
144 
145 
146  !-----------------------------------------------------------------------------
147  ! liquid densit data
148  !-----------------------------------------------------------------------------
149  DO i = 1, nliq
150 
151  IF ( pliq(i) > 0.0 ) THEN
152 
153  !-----------------------------------------------------------------------
154  ! case 1: calculate rho to given (p,T)
155  !-----------------------------------------------------------------------
156 
157  nocon = .false.
158  scan = .false.
159 
160  j_dat = i
161  nphas = 1
162  t = tliq(i)
163  p = pliq(i)
164  xi(1,1) = 1.0
165  IF ( vliq(i) <= v_crit ) THEN
166  vapor = .false.
167  densta(1) = 0.5
168  IF (d_kond(1,i) /= 0.0) densta(1) = d_kond(1,i)
169  ELSE
170  vapor = .true.
171  densta(1) = 0.002
172  IF (d_kond(1,i) /= 0.0) densta(1) = d_kond(1,i)
173  END IF
174 
175  vapor_scan: DO
176  IF (scan) densta(1) = densta(1) + 0.002
177 
178  CALL dens_calc ( rho_phas )
179  CALL si_dens ( density, w )
180  v_cal(i) = 1.0 / density(1)
181 
182  IF ((v_cal(i)-vliq(i))*1.e2/vliq(i) < 5.0 .AND. .NOT. scan) THEN
183  d_kond(1,i) = dense(1)
184  ELSE
185  d_kond(1,i) = 0.0
186  END IF
187 
188 
189  IF ( vapor .AND. (1.0/density(1)) <= (v_crit/1.1) ) THEN
190  IF ( densta(1) <= 0.228 ) THEN
191  scan = .true.
192  cycle vapor_scan
193  ELSE
194  nocon = .true.
195  END IF
196  END IF
197  EXIT vapor_scan
198  END DO vapor_scan
199  IF (dense(1) > 0.55) WRITE (*,*)'density ',i,dense(1)
200 
201  weigh = 0.5
202  IF (nocon) weigh = weigh / 4.0
203 
204  fvec(j_dat) = weigh * sqrt(type_weight(1))*(v_cal(i)-vliq(i))/ &
205  (sqrt(abs(vliq(i)))*sqrt(abs(v_cal(i))))
206 
207 
208 
209  ELSE
210 
211  !-----------------------------------------------------------------------
212  ! case 2, when pliq(i) = 0, indicating rho_L at coexistence
213  !-----------------------------------------------------------------------
214 
215  nphas = 2
216 
217  j_dat = i
218 
219  val_init(1) = 0.5
220  val_init(2) = 1.e-8
221  val_init(3) = tliq(i)
222  val_init(4) = 10.0
223  IF (d_kond(1,j_dat) /= 0.0) val_init(1) = d_kond(1,j_dat)
224  IF (d_kond(2,j_dat) /= 0.0) val_init(2) = d_kond(2,j_dat)
225  IF (d_kond(1,j_dat) /= 0.0 .AND. d_kond(2,j_dat) /= 0.0) val_init(4) = plv_kon(j_dat)
226  val_init(5) = 0.0 ! mole fraction lnx(1,1)
227  val_init(6) = 0.0 ! mole fraction lnx(2,1)
228 
229  CALL pure_equilibrium_fit ( crit_dat, tc, tliq(i), pliq(i), converg, v_cal(i), plv_dum )
230 
231  IF ( converg == 1 ) THEN
232  d_kond(1,j_dat) = dense(1)
233  d_kond(2,j_dat) = dense(2)
234  plv_kon(j_dat) = val_conv(4)
235  ELSE
236  d_kond(1,j_dat) = 0.0
237  d_kond(2,j_dat) = 0.0
238  plv_kon(j_dat) = 10.0
239  plv_dum = 1.e5
240  v_cal(i) = vliq(i) * 0.75
241  IF (crit_dat == 1 .AND. tliq(i) > tc ) v_cal(i) = vliq(i)*1.1
242  END IF
243 
244  weigh = 0.5
245 
246  fvec(j_dat) = weigh * sqrt(type_weight(1))*(v_cal(i)-vliq(i)) &
247  / (sqrt(abs(vliq(i)))*sqrt(abs(v_cal(i))))
248 
249  END IF
250 
251  END DO
252  !-----------------------------------------------------------------------------
253  ! end: liquid density data
254  !-----------------------------------------------------------------------------
255 
256 
257 
258  !-----------------------------------------------------------------------------
259  ! vapor pressure data
260  !-----------------------------------------------------------------------------
261  DO i = 1, nlv
262 
263  nphas = 2
264 
265  j_dat = i + nliq
266 
267  val_init(1) = 0.5
268  val_init(2) = 1.e-8
269  val_init(3) = tlv(i)
270  val_init(4) = plv(i)
271  IF (d_kond(1,j_dat) /= 0.0) val_init(1) = d_kond(1,j_dat)
272  IF (d_kond(2,j_dat) /= 0.0) val_init(2) = d_kond(2,j_dat)
273  IF (d_kond(1,j_dat) /= 0.0 .AND. d_kond(2,j_dat) /= 0.0) val_init(4) = plv_kon(j_dat)
274  val_init(5) = 0.0 ! mole fraction lnx(1,1)
275  val_init(6) = 0.0 ! mole fraction lnx(2,1)
276 
277  CALL pure_equilibrium_fit ( crit_dat, tc, tlv(i), plv(i), converg, v_dum, plvcal(i) )
278 
279  IF ( converg == 1 ) THEN
280  d_kond(1,j_dat) = dense(1)
281  d_kond(2,j_dat) = dense(2)
282  plv_kon(j_dat) = val_conv(4)
283  ELSE
284  d_kond(1,j_dat) = 0.0
285  d_kond(2,j_dat) = 0.0
286  plv_kon(j_dat) = plv(i)
287  IF ( crit_dat == 1 .AND. tc > tlv(i) ) THEN
288  plvcal(i) = plv(i) * 0.75
289  ! IF (i > 2) plvcal(i) = EXP( LOG(plvcal(i-1)) + (LOG(plvcal(i-1))-LOG(plvcal(i-2)))
290  ! /(1.0/tlv(i-1)-1.0/tlv(i-2)) * (1.0/tlv(i)-1.0/tlv(i-1)) )
291  WRITE (*,'(a,i3,2G14.5,G15.8)') ' loose ! ',i,tlv(i),t,plv(i)
292  ELSE
293  plvcal(i) = plv(i)
294  IF (i > 2) THEN
295  IF (plvcal(i-2) > 0.0 .AND. plvcal(i-1) > 0.0 .AND. &
296  tlv(i-1) /= tlv(i-2)) THEN
297  plvcal(i) = exp( log(plvcal(i-1)) + (log(plvcal(i-1))-log(plvcal(i-2))) &
298  /(1.0/tlv(i-1)-1.0/tlv(i-2)) *(1.0/tlv(i)-1.0/tlv(i-1)) )
299  END IF
300  END IF
301  ! WRITE (*,*) ' supercritical at point ', i, tlv(i), tc
302  END IF
303  END IF
304 
305  fvec(j_dat) = 0.2*sqrt(type_weight(2))*abs(plvcal(i)-plv(i)) &
306  / ((sqrt(abs(plv(i)))*sqrt(abs(plvcal(i)))))**0.85
307 
308  IF (tc < tlv(i) .AND. crit_dat == 1) THEN
309  fvec(j_dat)=fvec(j_dat) + 5.0*(tlv(i)-tc)
310  WRITE (*,'(a,i3,3(f15.5))') ' skipped ! ',i,tlv(i),tc
311  END IF
312 
313  END DO
314  !-----------------------------------------------------------------------------
315  ! end: vapor pressure data
316  !-----------------------------------------------------------------------------
317 
318 
319  !-----------------------------------------------------------------------------
320  ! 2nd virial coefficients
321  !-----------------------------------------------------------------------------
322 
323  DO i = 1, n_vir
324 
325  j_dat = i + nliq + nlv
326 
327  nphas=1
328  dense(1)= 1.e-7
329  t = t_vir(i)
330  p = 1.0
331  xi(1,1) = 1.0
332 
333  CALL p_calc (pges,zges)
334  CALL si_dens (density,w)
335 
336  b_cal(i) = (zges-1.0)/density(1)*1000.0*mm(1)
337 
338  weigh = (parame(1,3)/t)**3
339  fvec(j_dat) = type_weight(3)*weigh*(b_cal(i)-b_vir(i)) &
340  / (sqrt(abs(b_vir(i)))*sqrt(abs(b_cal(i))))
341 
342  END DO
343 
344 
345 
346  !-----------------------------------------------------------------------------
347  ! internal energy data
348  !-----------------------------------------------------------------------------
349 
350  DO i = 1, n_en
351 
352  j_dat = i + nliq + nlv + n_vir
353 
354  nphas = 1
355  dense(1) = rho_en(i)
356  t = t_en(i)
357  xi(1,1) = 1.0
358 
359  ! CALL P_CALC (pges,zges)
360  ! CALL F_EOS (...)
361 
362  ! u_res = f_res + T*S
363 
364  weigh = 1.0
365  fvec(j_dat) = type_weight(4)*weigh*(u_cal(i)-u_en(i))/ abs(u_en(i))
366 
367  END DO
368 
369  IF (j_dat /= lm_m_dat) write (*,*) 'array length not matching!'
370 
371 
372  !-----------------------------------------------------------------------------
373  ! critical point data (work in progress)
374  !-----------------------------------------------------------------------------
375 
376  !-----------critical point---------------------------------------------
377  ! naphtalene
378  ! tclit = 748.210
379  ! pclit = 4100000.0
380  ! rhoclit = 246.0
381  ! co2
382  ! tclit = 304.210
383  ! pclit = 7382500.0
384  ! rhoclit = 466.01015160
385  ! propane
386  ! tclit = 370.00
387  ! pclit = 4.265d+06
388  ! rhoclit = 225.00
389  ! butane
390  ! tclit = 425.20
391  ! pclit = 3.796d+06
392  ! rhoclit = 227.00
393  ! PC-SAFT: T=432.1 P=41.95
394  ! hexane
395  ! tclit = 507.40
396  ! pclit = 29.7d5
397  ! rhoclit = 233.50
398  ! heptane
399  ! tclit = 539.70
400  ! pclit = 2.736d+06
401  ! rhoclit = 234.10
402  ! octane
403  ! tclit = 569.40
404  ! pclit = 2.496d+06
405  ! rhoclit = 234.70
406  ! acetone
407  ! tclit = 508.15
408  ! pclit = 4.7621E+06
409  ! rhoclit = 234.70
410 
411  ! tc = tclit
412  ! pc = pclit
413  ! rhoc = rhoclit
414 
415  ! CALL CRITICAL (tc,pc,rhoc)
416 
417  ! write (*,'(t5,a,f7.2,a,f7.3,a,f7.2)') 'tc =',tc,' pc =',pc/1.E5,' rho =',rhoc
418  ! write (*,'(t5,a,f7.2,a,f7.3,a,f7.2)') 'tc =',tclit,' pc =',pclit/1.E5,' rho =',rhoclit
419 
420  !c j_dat = nliq + nlv + n_vir + n_en + 1
421  !c fvec(j_dat)=0.020*SQRT(ABS(rhoc/rhoclit-1.0))
422  ! j_dat = nliq + nlv + n_vir + n_en + 2
423  ! fvec(j_dat) = fvec(j_dat) + 1.0 * SQRT(ABS(tc/tclit-1.0))
424  ! j_dat = nliq + nlv + n_vir + n_en + 3
425  ! fvec(j_dat)=1.0*SQRT(ABS(pc/pclit-1.0))
426 
427  !-----------critical point---------------------------------------------
428 
429  ! DO i = 1, j_dat
430  ! write (*,*) i,fvec(i)
431  ! ENDDO
432  ! call pause
433 
434 END SUBROUTINE pure_residual
435 
436 
437 
438 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
439 ! SUBROUTINE pure_equilibrium_fit
440 !
441 ! input to the subroutine: val_init has to be available
442 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
443 
444 SUBROUTINE pure_equilibrium_fit (crit_dat, tc, t_in, p_in, converg, v_cal, plvcal)
445 
446  USE basic_variables
447  USE eos_variables, ONLY: densav, eta_start, pges
448  use utilities, only: si_dens, critical
449  IMPLICIT NONE
450 
451  !-----------------------------------------------------------------------------
452  INTEGER, INTENT(IN OUT) :: crit_dat
453  REAL, INTENT(IN OUT) :: tc
454  REAL, INTENT(IN) :: t_in
455  REAL, INTENT(IN) :: p_in
456  INTEGER, INTENT(OUT) :: converg
457  REAL, INTENT(OUT) :: v_cal
458  REAL, INTENT(OUT) :: plvcal
459  !-----------------------------------------------------------------------------
460  REAL :: temp, press
461  REAL :: pc, rhoc
462  REAL :: density(np), w(np,nc)
463  REAL :: p_high, p_low
464  !-----------------------------------------------------------------------------
465 
466 
467  n_unkw = ncomp ! number of quantities to be iterated
468  it(1) = 'lnp' ! iteration of pressure
469  val_conv(2) = 0.0
470  densav(:) = 0.0
471 
472  temp = t_in
473  press = p_in
474 
475  plvcal = 0.0
476  v_cal = 0.0
477 
478 
479  !-----------------------------------------------------------------------------
480  ! phase equilibrium calculation
481  !-----------------------------------------------------------------------------
482 
483  CALL objective_ctrl (converg)
484 
485 
486  !-----------------------------------------------------------------------------
487  ! for the case, where no convergence was found initially
488  !-----------------------------------------------------------------------------
489 
490  ! --- calc critical pt. in order to assess, whether a solution is possible
491  IF ( converg /= 1 .AND. crit_dat == 0 ) THEN
492  tc = temp
493  pc = press
494  CALL critical ( tc, pc, rhoc )
495  IF ( tc < 0.0 ) write (*,*) 'error: negative Tc calculated'
496  IF ( tc < 0.0 ) tc = 500.0
497  crit_dat = 1
498  END IF
499 
500  !-----------------------------------------------------------------------------
501  ! 1st step to achieve convergence: select feasible p (betw. spinodal p's)
502  !-----------------------------------------------------------------------------
503  IF ( converg /= 1 .AND. tc > temp ) THEN
504  t = temp
505  CALL perturbation_parameter
506  eta_start = 1.e-8
507  CALL pressure_spinodal ! vapor spinodal: max p for which a vapor exists
508  p_high = pges
509  eta_start = 0.45
510  CALL pressure_spinodal ! liquid spinodal: min p for which a liquid exists
511  p_low = pges
512  val_init(1) = 0.5
513  val_init(2) = 1.e-8
514  val_init(3) = temp
515  val_init(4) = 0.8 * p_high + 0.2 * max( p_low, 0.0 ) ! choose p betw. the 2 spinodals
516  CALL objective_ctrl (converg)
517  ENDIF
518 
519  !-----------------------------------------------------------------------------
520  ! 2nd step to achieve convergence: approach t from lower t in a few steps
521  !-----------------------------------------------------------------------------
522  IF ( converg /= 1 .AND. tc > temp ) THEN
523  val_init(1) = 0.5
524  val_init(2) = 1.e-6
525  val_init(3) = temp - 0.1*temp
526  val_init(4) = press / 10.0
527  IF (val_init(4) == 0.0) val_init(4) = exp( ( 1.0 - tc/val_init(3) )*7.0 ) * pc ! simple correspondence principle
528  outp = 0 ! output to terminal (set u_out_p before turning output on)
529  running = 't' ! Temperature is running var. in PHASE_EQUILIB
530  CALL phase_equilib ( temp, 4.0, converg )
531  IF ( converg == 1 ) write (*,'(a,G14.5,G15.8)') ' win ! ',temp, p
532  END IF
533 
534  IF ( converg == 1 ) THEN
535  plvcal = val_conv(4)
536  CALL si_dens (density,w)
537  v_cal = 1.0 / density(1)
538  END IF
539 
540 END SUBROUTINE pure_equilibrium_fit
541 
542 
543 
544 
545 
546 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
547 ! SUBROUTINE fitting
548 !
549 ! Driver for the fitting of pure component parameters. This subroutine
550 ! uses LMDIF1 (Levenberg-Marquardt MINPACK routine) in order to minimze
551 ! the deviation of calculation results to experimental values. The
552 ! deviations are calculated in PURE_RESIDUAL and stored in 'fvec'.
553 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
554 
555 SUBROUTINE fitting
556 
557  USE basic_variables
558  USE levenberg_marquardt
559  USE pure_fit_parameters, ONLY: pure_fit_file
560  use utilities, only: file_open
561  IMPLICIT NONE
562 
563  !-----------------------------------------------------------------------------
564  INTERFACE
565  SUBROUTINE pure_residual(lm_m_dat, lm_n_par, para, fvec, iflag)
566  IMPLICIT NONE
567  INTEGER, INTENT(IN) :: lm_m_dat, lm_n_par
568  REAL, INTENT(IN) :: para(:)
569  REAL, INTENT(IN OUT) :: fvec(:)
570  INTEGER, INTENT(IN OUT) :: iflag
571  END SUBROUTINE pure_residual
572 
573  SUBROUTINE paini( lm_n_par, para )
574  IMPLICIT NONE
575  INTEGER, INTENT(OUT) :: lm_n_par
576  REAL, ALLOCATABLE, INTENT(OUT) :: para(:)
577  END SUBROUTINE paini
578 
579  SUBROUTINE pure_output( lm_n_par, para )
580  IMPLICIT NONE
581  INTEGER, INTENT(IN) :: lm_n_par
582  REAL, INTENT(IN) :: para(:)
583  END SUBROUTINE pure_output
584  END INTERFACE
585 
586  !-----------------------------------------------------------------------------
587  INTEGER :: info
588  INTEGER :: lm_n_par ! number of adjustable parameters
589  INTEGER :: lm_m_dat ! number of data points
590  INTEGER, ALLOCATABLE :: ipvt(:)
591  REAL, ALLOCATABLE :: para(:)
592  REAL, ALLOCATABLE :: fvec(:)
593  REAL :: tol = 1.0e-8
594  REAL :: epsfcn
595  REAL :: lm_factor
596  !-----------------------------------------------------------------------------
597 
598 
599  ncomp = 1
600 
601  ! --- numerical constants for Levenberg-Marquardt algorithm ------------
602  scaling(1) = 1.e2
603  IF ( num == 0 ) epsfcn = 1.0e-6**2 ! sqrt of relat. step size (finite differences)
604  IF ( num == 1 ) epsfcn = 1.0e-5**2 ! sqrt of relat. step size (finite differences)
605  lm_factor = 0.1 ! maximum initial step size (parameter iteration)
606 
607  IF ( num == 2 ) write (*,*) 'fitting with RGT is not yet supported. Set num <= 1'
608  IF ( num == 2 ) stop
609 
610  !-----------------------------------------------------------------------------
611  ! read input-file (DATIN), read parameters (PAINI), prepare output
612  !-----------------------------------------------------------------------------
613 
614  WRITE (*,*) ' SPECIFY INPUT FILE: ( ./input_file/pure_comp/ )'
615  READ (*,*) pure_fit_file
616  ! pure_fit_file = '3939.dat' !JG
617 
618  CALL file_open('./input_file/pure_comp/'//pure_fit_file,22)
619  OPEN (23,file='./output_file/pure_comp/'//pure_fit_file)
620  pure_fit_file='./input_file/pure_comp/'//pure_fit_file
621 
622  CALL datin( lm_m_dat )
623  CALL paini( lm_n_par, para )
624 
625  ALLOCATE( fvec(lm_m_dat) )
626  ALLOCATE( ipvt(lm_n_par) )
627 
628 
629  !-----------------------------------------------------------------------------
630  ! adjust parameters (Levenberg-Marquardt scheme)
631  !-----------------------------------------------------------------------------
632 
633  CALL lmdif1 (pure_residual,lm_m_dat,lm_n_par,para,fvec,tol,epsfcn,lm_factor,info,ipvt)
634 
635 
636  !-----------------------------------------------------------------------------
637  ! algorithm output: optimization status
638  !-----------------------------------------------------------------------------
639 
640  WRITE (*,*) ' '
641  WRITE (23,*) ' '
642 
643  SELECT CASE (info)
644  CASE (:-1)
645  WRITE (23,*) 'SOLVER STATUS: Users FCN returned INFO = ', -info
646  WRITE (*, *) 'SOLVER STATUS: Users FCN returned INFO = ', -info
647  CASE (0)
648  WRITE (23,*) 'SOLVER STATUS: Improper values for input parameters'
649  WRITE (*, *) 'SOLVER STATUS: Improper values for input parameters'
650  CASE (1:3)
651  WRITE (23,*) 'SOLVER STATUS: Convergence'
652  WRITE (*, *) 'SOLVER STATUS: Convergence'
653  CASE (4)
654  WRITE (23,*) 'SOLVER STATUS: Residuals orthogonal to the Jacobian'
655  WRITE (*, *) 'SOLVER STATUS: Residuals orthogonal to the Jacobian'
656  WRITE (*, *) 'There may be an error in FCN'
657  CASE (5)
658  WRITE (23,*) 'SOLVER STATUS: Too many calls to FCN'
659  WRITE (*, *) 'SOLVER STATUS: Too many calls to FCN'
660  WRITE (*, *) 'Either slow convergence, or an error in FCN'
661  CASE (6:7)
662  WRITE (23,*) 'SOLVER STATUS: TOL was set too small'
663  WRITE (*, *) 'SOLVER STATUS: TOL was set too small'
664  CASE DEFAULT
665  WRITE (23,*) 'SOLVER STATUS: INFO =', info, ' ???'
666  WRITE (*, *) 'SOLVER STATUS: INFO =', info, ' ???'
667  END SELECT
668 
669  WRITE (*,*) '-------------'
670  WRITE (23,*) '-------------'
671 
672  CALL pure_output ( lm_n_par, para )
673 
674  CLOSE (22)
675  CLOSE (23)
676 
677  DEALLOCATE( fvec, ipvt )
678 
679 END SUBROUTINE fitting
680 
681 
682 
683 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
684 ! SUBROUTINE DATIN
685 !
686 ! Read experimental pure component data for the regression of pure
687 ! component parameters. This file also writes the experimental data to
688 ! the output file.
689 ! So far, the experimental data is stored in vectors of fixed length
690 ! (see module 'pure_fit_parameters'). At some point this should be
691 ! changed - allocate the vector length here.
692 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
693 
694 SUBROUTINE datin ( lm_m_dat )
695 
696  USE basic_variables
697  USE pure_fit_parameters
698  IMPLICIT NONE
699 
700  !-----------------------------------------------------------------------------
701  INTEGER, INTENT(OUT) :: lm_m_dat
702 
703  !-----------------------------------------------------------------------------
704  INTEGER :: mak, i, k
705  INTEGER :: datatype_available(cont)
706  INTEGER :: data_no, data_low, data_up, alldata
707  REAL :: mmpliq
708  CHARACTER (LEN=80) :: info_line
709  !-----------------------------------------------------------------------------
710 
711  nlv = 0
712  nliq = 0
713  n_vir = 0
714  n_en = 0
715 
716  eqa = 0
717 
718  !-----------------------------------------------------------------------------
719  ! reading info-lines (input-file) and writing them (output-file)
720  !-----------------------------------------------------------------------------
721 
722  DO i = 1, 15
723  READ(22,'(A80)') info_line
724  WRITE(23,'(a)') info_line
725  END DO
726 
727  ! --- reading molecular mass and an estimate of the crit. density ------------
728  READ(22,'(A80)') info_line
729  WRITE(23,'(A)') info_line
730  READ(22,*) mmpliq, v_crit
731  WRITE(23, '(2(2x,G13.6))') mmpliq, v_crit
732  mm(1) = mmpliq
733 
734  DO i=1,14
735  READ(22,'(A80)') info_line
736  WRITE(23,'(A)') info_line
737  END DO
738 
739  ! --- reading the number of experimental data points -------------------------
740  READ(22,'(A80)') info_line
741  READ(22,*) (datatype_available(i),i=1,cont)
742  WRITE (*,*) ' '
743 
744  eos = 1
745 
746  !-----------------------------------------------------------------------------
747  ! selection of data sets by user
748  !-----------------------------------------------------------------------------
749 
750  WRITE (*,*), ' DATA SELECTION'
751  WRITE (*,*) ' '
752  WRITE (*,*), ' NUMBER OF DATA-SETS:'
753  WRITE (*,'(T2,A40,I3)') 'PVT - data :', datatype_available(1)
754  WRITE (*,'(T2,A40,I3)') 'vapor pressure data :', datatype_available(2)
755  WRITE (*,'(T2,A40,I3)') '2nd virial coefficients :', datatype_available(3)
756  WRITE (*,'(T2,A40,I3)') 'internal res. energy data :', datatype_available(4)
757  WRITE (*,*) ' '
758 
759  ! --- option to select all data sets -----------------------------------------
760  WRITE (*,*), ' SELECT DATA-SET !'
761  WRITE (*,*), ' consider all available data: 1'
762  WRITE (*,*), ' select data-sets individually: 0'
763  READ (*,*) alldata
764  ! alldata = 1 !JG
765  IF ( alldata == 0 ) WRITE (*,*), ' SELECT DATA-TYPES !'
766  WRITE (*,*) ' '
767  IF (datatype_available(1) /= 0) THEN
768  eqa(1)=1
769  IF (alldata == 0) WRITE (*,*), ' PVT-data ? (0/1) '
770  IF (alldata == 0) READ (*,*) eqa(1)
771  END IF
772  IF (datatype_available(2) /= 0) THEN
773  eqa(2)=1
774  IF (alldata == 0) WRITE (*,*), ' vapor pressure data ? (0/1)'
775  IF (alldata == 0) READ (*,*) eqa(2)
776  END IF
777  IF (datatype_available(3) /= 0) THEN
778  eqa(3)=1
779  IF (alldata == 0) WRITE (*,*), ' 2nd virial coefficients ? (0/1)'
780  IF (alldata == 0) READ (*,*) eqa(3)
781  END IF
782  IF (datatype_available(4) /= 0) THEN
783  eqa(4)=1
784  IF (alldata == 0) WRITE (*,*), ' internal res. energy data? (0/1)'
785  IF (alldata == 0) READ (*,*) eqa(4)
786  END IF
787 
788  WRITE (*,*) ' '
789 
790  !-----------------------------------------------------------------------------
791  ! reading data sets from input file
792  !-----------------------------------------------------------------------------
793 
794  DO i = 1, cont ! --- cont: number of data types (p^sat, rho_L, ...)
795 
796  ! --- two headers ---------------------------------------------------------
797  IF (datatype_available(i) /= 0) READ(22,'(A80)') info_line
798  IF (eqa(i) == 1) WRITE(23,'(a)') info_line
799  IF (datatype_available(i) /= 0) READ(22,'(A80)') info_line
800  IF (eqa(i) == 1) WRITE(23,'(a)') info_line
801 
802  ! --- reading data sets if a certain data type i was chosen ---------------
803  WRITE (*,*) ' '
804  IF (eqa(i) == 1) THEN
805  IF (i == 1) WRITE (*,*), ' CHOOSE FROM PVT-DATA:'
806  IF (i == 2) WRITE (*,*), ' CHOOSE FROM VAPOR PRESSURE DATA:'
807  IF (i == 3) WRITE (*,*), ' CHOOSE FROM 2nd VIRIAL COEFFICIENT DATA:'
808  IF (i == 4) WRITE (*,*), ' CHOOSE FROM INTERNAL RESID. ENERGY DATA:'
809 
810  data_low = 1
811  data_up = datatype_available(i)
812 
813  ! --- user selects the considered data sets ----------------------------
814  WRITE (*,'(T2,A,I3)') ' Total number of data-sets: ',datatype_available(i)
815  WRITE (*,*), ' Specify the data-sets to be considered'
816  WRITE (*,*), ' Lower number of data-set (e.g. 1):'
817  IF (alldata == 0) READ (*,*) data_low
818  WRITE (*,'(a,I3,a)') ' Upper number of data-set (e.g.',datatype_available(i),'):'
819  IF (alldata == 0) READ (*,*) data_up
820 
821  ! --- check the input --------------------------------------------------
822  IF ((data_low > data_up) .OR. (data_up > datatype_available(i))) THEN
823  WRITE (*,*) ' Erroneous input! The lower data-set number must be smaller'
824  WRITE (*,*) ' than the upper value.'
825  WRITE (*,*) ' The upper data-set number must not be greater than ',datatype_available(i)
826  WRITE (*,*) ' Lower number of data-set (e.g. 1):'
827  READ (*,*) data_low
828  WRITE (*,'(a,I3,a)') ' Upper number of data-set (e.g.',datatype_available(i),'):'
829  READ (*,*) data_up
830  END IF
831 
832  mak = 0
833 
834  ! --- reading data -----------------------------------------------------
835  do k = 1, data_low-1
836  read (22,*)
837  end do
838  IF (i == 1) THEN
839  nliq = 0
840  DO k = data_low, data_up
841  nliq = nliq + 1
842  READ (22,*) data_no, pliq(nliq), tliq(nliq), vliq(nliq)
843  WRITE (*,'(i3,3(2x,G12.5))') data_no, pliq(nliq), tliq(nliq), vliq(nliq)
844  WRITE(23,'(i3,3(2x,G12.5))') data_no, pliq(nliq), tliq(nliq), vliq(nliq)
845  END DO
846  ELSE IF (i == 2) THEN
847  nlv = 0
848  DO k = data_low, data_up
849  nlv = nlv + 1
850  READ (22,*) data_no, plv(nlv), tlv(nlv)
851  WRITE (*,'(i3,2(2x,G12.5))') data_no, plv(nlv), tlv(nlv)
852  WRITE(23,'(i3,2(2x,G12.5))') data_no, plv(nlv), tlv(nlv)
853  END DO
854  ELSE IF (i == 3) THEN
855  n_vir = 0
856  DO k = data_low, data_up
857  n_vir = n_vir + 1
858  READ (22,*) data_no, b_vir(n_vir), t_vir(n_vir)
859  WRITE (*,'(i3,2(2x,G12.5))') data_no, b_vir(n_vir), t_vir(n_vir)
860  WRITE(23,'(i3,2(2x,G12.5))') data_no, b_vir(n_vir), t_vir(n_vir)
861  END DO
862  ELSE IF (i == 4) THEN
863  n_en = 0
864  DO k = data_low, data_up
865  n_en = n_en + 1
866  READ (22,*) data_no, u_en(n_en), rho_en(n_en), t_en(n_en)
867  WRITE (*,'(i3,3(2x,G12.5))') data_no, u_en(n_en), rho_en(n_en), t_en(n_en)
868  WRITE(23,'(i3,3(2x,G12.5))') data_no, u_en(n_en), rho_en(n_en), t_en(n_en)
869  END DO
870  END IF
871  do k = data_up+1, datatype_available(i)
872  read (22,*)
873  end do
874 
875  ELSE
876 
877  ! --- a certain data type was not selected for the fitting -------------
878  IF (datatype_available(i) /= 0) THEN
879  DO k = 1, datatype_available(i)
880  READ(22,*) data_no
881  END DO
882  END IF
883 
884  END IF
885 
886  END DO
887 
888  !-----------------------------------------------------------------------------
889  ! total number of data sets considered in parameter fitting
890  !-----------------------------------------------------------------------------
891 
892  lm_m_dat = nliq + nlv + n_vir + n_en ! excluding crit. point
893  ! lm_m_dat = lm_m_dat + 3 ! including crit. point
894 
895  rewind(22)
896 
897 END SUBROUTINE datin
898 
899 
900 
901 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
902 ! SUBROUTINE paini
903 !
904 ! Subroutine for reading in starting values of the pure component
905 ! parameters
906 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
907 
908 SUBROUTINE paini ( lm_n_par, para )
909 
910  USE basic_variables
911  USE pure_fit_parameters
912  USE utilities
913  IMPLICIT NONE
914 
915  !-----------------------------------------------------------------------------
916  INTEGER, INTENT(OUT) :: lm_n_par
917  REAL, ALLOCATABLE, INTENT(OUT) :: para(:)
918 
919  !-----------------------------------------------------------------------------
920  INTEGER :: fitopt, position, i_fit
921  INTEGER :: n_sites
922  INTEGER :: i, j, assoc, quadrup, dipole
923  REAL :: dummy( npmax, 3 )
924  INTEGER :: k, no, nhb_typ(nc)
925  REAL :: nhb_no(nc,nsite)
926  REAL :: eps_hb(nc,nc,nsite,nsite)
927  !-----------------------------------------------------------------------------
928 
929  type_weight = 1.0
930 
931  assoc = 0
932  dipole = 0
933  quadrup= 0
934 
935  nhb_typ(:) = 0
936  nhb_no(:,:) = 0.0
937 
938  !-----------------------------------------------------------------------------
939  ! selecting adjustable parameters
940  !-----------------------------------------------------------------------------
941 
942  WRITE (*,*) ' '
943  WRITE (*,*), '*********** PARAMETER SPECIFICATION ***********'
944  WRITE (*,*) ' '
945 
946  WRITE (*,*) ' You may choose from two default fitting options'
947  WRITE (*,*) ' or individually specify all parameters.'
948  WRITE (*,*) ' (Dipolar and quadrupolar moments are set to zero'
949  WRITE (*,*) ' for both default options; for dipolar or'
950  WRITE (*,*) ' quadrupolar compounds choose option 3)'
951  WRITE (*,*) ' '
952  WRITE (*,*) ' Choose one the following options'
953  WRITE (*,*) ' Default 1: Fitting 3 parameters'
954  WRITE (*,*) ' * segment number m [ ]'
955  WRITE (*,*) ' * segment diameter sigma [A]'
956  WRITE (*,*) ' * segment energy parameter eps/k [K]'
957  WRITE (*,*) ' Default 2: Fitting 5 parameters for associating'
958  WRITE (*,*) ' compounds with 2 association sites'
959  WRITE (*,*) ' * segment number m [ ]'
960  WRITE (*,*) ' * segment diameter sigma [A]'
961  WRITE (*,*) ' * segment energy parameter eps/k [K]'
962  WRITE (*,*) ' * assoc.-energy parameter eps_AB/k [K]'
963  WRITE (*,*) ' * assoc.-volume parameter kappa_AB [ ]'
964  WRITE (*,*) ' Option 3: Individual specification of all'
965  WRITE (*,*) ' parameters to be fitted.'
966  WRITE (*,*) ' Choose 1, 2, or 3.'
967  READ (*,*) fitopt
968  ! fitopt = 3 !JG
969 
970  IF (fitopt == 1) THEN
971  lm_n_par = 3
972  !parame(1,12) = 0.0 ! number of assoc. sites (real-type !)
973  nhb_typ(1) = 0 ! 11.03.2015
974  nhb_no(1,1) = 0.0 ! 11.03.2015
975  aa_txt(1) = ' Segm.No.m / Molar Mass :'
976  aa_txt(2) = ' Segm.Diameter sigma :'
977  aa_txt(3) = ' Segm.Energy Param eps/k:'
978  fit_array(1) = 1
979  fit_array(2) = 2
980  fit_array(3) = 3
981  DO i = 4, 10
982  fit_array(i) = 0
983  ENDDO
984  ELSEIF (fitopt == 2) THEN
985  lm_n_par = 5
986  !parame(1,12) = 2.0 ! number of assoc. sites (real-type !)
987  aa_txt(1) = ' Segm.No.m / Molar Mass :'
988  aa_txt(2) = ' Segm.Diameter sigma :'
989  aa_txt(3) = ' Segm.Energy Param eps/k:'
990  aa_txt(4) = ' Assoc.-energy eps_AB/k :'
991  aa_txt(5) = ' Assoc.-volume kappa_AB :'
992  fit_array(1) = 1
993  fit_array(2) = 2
994  fit_array(3) = 3
995  fit_array(4) = 4
996  fit_array(5) = 5
997  DO i = 6, 10
998  fit_array(i)=0
999  ENDDO
1000  nhb_typ(1) = 2
1001  nhb_no(1,1) = 1.0
1002  nhb_no(1,2) = 1.0
1003  ELSE
1004  WRITE (*,*) ' The following pure component parameters can'
1005  WRITE (*,*) ' be fitted'
1006  WRITE (*,*) ' (1) segment number m [ ]'
1007  WRITE (*,*) ' (2) segment diameter sigma [A]'
1008  WRITE (*,*) ' (3) segment energy parameter eps/k [K]'
1009  WRITE (*,*) ' (4) assoc.-energy parameter eps_AB/k [K]'
1010  WRITE (*,*) ' (5) assoc.-volume parameter kappa_AB [ ]'
1011  WRITE (*,*) ' (6) dipolar moment my [D]'
1012  WRITE (*,*) ' (7) fraction of segm.with dipolar moment xp [ ]'
1013  WRITE (*,*) ' (8) quadrupolar moment Q [D]'
1014  WRITE (*,*) ' (9) fraction of segm.with quadrup.moment xp [ ]'
1015  WRITE (*,*) ' (11)polarizability alpha [A^3]'
1016  WRITE (*,*) ' '
1017  WRITE (*,*) ' Note, that the number of association sites has'
1018  WRITE (*,*) ' to be user-specified'
1019  WRITE (*,*) ' '
1020  WRITE (*,*) ' Choose the total number of fitting-parameters'
1021  READ (*,*) lm_n_par
1022  ! lm_n_par = 3 !JG
1023  WRITE (*,*) ' Choose the position-number of fitting-param.'
1024  WRITE (*,*) ' (see above list of parameter 1 to 9)'
1025  DO i = 1, lm_n_par
1026  WRITE (*,'(a,i1,a,i2)') ' position-number for param. ',i,' of ',lm_n_par
1027  READ (*,*) position
1028  IF (position == 1) aa_txt(i) = ' Segm.No.m / Molar Mass :'
1029  IF (position == 2) aa_txt(i) = ' Segm.Diameter sigma :'
1030  IF (position == 3) aa_txt(i) = ' Segm.Energy Param eps/k:'
1031  IF (position == 4 .OR. position == 5) assoc = 1
1032  IF (position == 4) aa_txt(i) = ' Assoc.-energy eps_AB/k :'
1033  IF (position == 5) aa_txt(i) = ' Assoc.-volume kappa_AB :'
1034  IF (position == 6 .OR. position == 7) dipole = 1
1035  IF (position == 6) aa_txt(i) = ' Dipolar moment :'
1036  IF (position == 7) aa_txt(i) = ' No. of segm.with dipole:'
1037  IF (position == 8 .OR. position == 9) quadrup = 1
1038  IF (position == 8) aa_txt(i) = ' Quadrupolar moment :'
1039  IF (position == 9) aa_txt(i) = ' No.of segm.with quadrup:'
1040  IF (position == 11)aa_txt(i) = ' polarizability alpha'
1041  fit_array(i) = position
1042  ENDDO
1043  IF ( assoc == 1 ) THEN
1044  WRITE (*,*) ' Give number of association sites'
1045  READ (*,*) n_sites
1046  IF (n_sites == 1) nhb_typ(1) = 1
1047  IF (n_sites == 1) nhb_no(1,1)= 1.0
1048  IF (n_sites == 2) nhb_typ(1) = 2
1049  IF (n_sites == 2) nhb_no(1,1)= 1.0
1050  IF (n_sites == 2) nhb_no(1,2)= 1.0
1051  IF (n_sites == 3) nhb_typ(1) = 2
1052  IF (n_sites == 3) nhb_no(1,1)= 1.0
1053  IF (n_sites == 3) nhb_no(1,2)= 2.0
1054  IF (n_sites == 4) nhb_typ(1) = 2
1055  IF (n_sites == 4) nhb_no(1,1)= 2.0
1056  IF (n_sites == 4) nhb_no(1,2)= 2.0
1057  ! parame(1,12) = REAL( nhb_typ(1) )
1058 
1059  ENDIF
1060  ENDIF
1061 
1062  ALLOCATE( para(lm_n_par) )
1063 
1064 
1065  !-----------------------------------------------------------------------------
1066  ! specifying starting values of adjustable parameters
1067  !-----------------------------------------------------------------------------
1068 
1069  WRITE (*,*) ' '
1070  WRITE (*,*), ' Input of starting values for fitting-parameters:'
1071  WRITE (*,*) ' '
1072 
1073  WRITE (*,*), ' GIVE STARTING VALUES FOR PARAMETERS '
1074  WRITE (*,*) ' '
1075 
1076  DO i = 1, lm_n_par
1077  WRITE (*,'(2a,i2,a,i2,a)') aa_txt(i),' (param.',i,' of',lm_n_par,')'
1078  READ (*,*) (dummy(i,j), j = 1,1)
1079  dummy(i,2) = dummy(i,1) / 10.0
1080  dummy(i,3) = dummy(i,1) / dummy(i,2)
1081  END DO
1082  WRITE (*,*) ' '
1083 
1084  !-----------------------------------------------------------------------------
1085  ! for case, where param. are selected individually, determine rest
1086  !-----------------------------------------------------------------------------
1087 
1088  n_sites = nint( sum( nhb_no( 1, 1:nhb_typ(1) ) ) )
1089  IF ( fitopt >= 3 ) THEN
1090  WRITE (*,*) ' '
1091  WRITE (*,*) ' The remaining parameters (not-fitted) now'
1092  WRITE (*,*) ' have to be defined:'
1093  WRITE (*,*) ' '
1094  DO i = 1,11
1095  i_fit = 0
1096  DO j = 1,lm_n_par
1097  IF (fit_array(j) == i) i_fit = 1
1098  ENDDO
1099  IF (i_fit == 0) THEN
1100  IF (i == 1) THEN
1101  write (*,*) ' Give segm.No.m / Molar Mass :'
1102  READ (*,*) parame(1,1)
1103  parame(1,1) = parame(1,1)*mm(1)
1104  ELSEIF (i == 2) THEN
1105  write (*,*) ' Give segm.Diameter sigma :'
1106  READ (*,*) parame(1,2)
1107  ELSEIF (i == 3) THEN
1108  write (*,*) ' Give segm.Energy Param eps/k:'
1109  READ (*,*) parame(1,3)
1110  ELSEIF (i == 4) THEN
1111  IF ( n_sites == 0) THEN
1112  write (*,*) ' Give number of association sites,'
1113  write (*,*) ' choose 0 for non-associating compounds'
1114  READ (*,*) n_sites
1115  ENDIF
1116  IF (n_sites /= 0) THEN
1117  write (*,*) ' Give assoc.-energy eps_AB/k :'
1118  IF (n_sites == 1) nhb_typ(1) = 1
1119  IF (n_sites == 1) nhb_no(1,1)= 1.0
1120  IF (n_sites == 1) READ (*,*) eps_hb(1,1,1,1)
1121  IF (n_sites == 2) nhb_typ(1) = 2
1122  IF (n_sites == 2) nhb_no(1,1)= 1.0
1123  IF (n_sites == 2) nhb_no(1,2)= 1.0
1124  IF (n_sites == 2) READ (*,*) eps_hb(1,1,1,2)
1125  IF (n_sites == 2) eps_hb(1,1,2,1) = eps_hb(1,1,1,2)
1126  IF (n_sites == 3) nhb_typ(1) = 2
1127  IF (n_sites == 3) nhb_no(1,1)= 1.0
1128  IF (n_sites == 3) nhb_no(1,2)= 2.0
1129  IF (n_sites == 3) READ (*,*) eps_hb(1,1,1,2)
1130  IF (n_sites == 3) READ (*,*) eps_hb(1,1,2,1)
1131  IF (n_sites == 4) nhb_typ(1) = 2
1132  IF (n_sites == 4) nhb_no(1,1)= 2.0
1133  IF (n_sites == 4) nhb_no(1,2)= 2.0
1134  IF (n_sites == 4) READ (*,*) eps_hb(1,1,1,2)
1135  IF (n_sites == 4) READ (*,*) eps_hb(1,1,2,1)
1136  ! parame(1,14) = 0.0 ! eps_hb(1,1,1,1)
1137  ! READ (*,*) parame(1,15) ! eps_hb(1,1,1,2)
1138  ! parame(1,16) = parame(1,15) ! eps_hb(1,1,2,1)
1139  ! parame(1,17) = 0.0 ! eps_hb(1,1,2,2)
1140  ENDIF
1141  ELSEIF (i == 5 .AND. n_sites /= 0.0) THEN
1142  write (*,*) ' Give assoc.-volume kappa_AB :'
1143  READ (*,*) parame(1,13)
1144  ! parame(1,13) = 0.01 !JG
1145  ELSEIF (i == 6) THEN
1146  write (*,*) ' Give dipolar moment :'
1147  READ (*,*) parame(1,6)
1148  ! parame(1,6) = 4.886577 !JG
1149  ELSEIF (i == 7 .AND. (parame(1,6) /= 0.0 .OR. dipole == 1))THEN
1150  ! write (*,*) ' Give No. of segm.with dipole:'
1151  ! READ (*,*) parame(1,8)
1152  ELSEIF (i == 8) THEN
1153  write (*,*) ' Give quadrupolar moment :'
1154  READ (*,*) parame(1,7)
1155  ! parame(1,7) = 0.0 !JG
1156  ELSEIF(i == 9 .AND. (parame(1,7) /= 0.0 .OR. quadrup == 1))THEN
1157  ! write (*,*) ' Give No. of segm.with quadrp:'
1158  ! READ (*,*) parame(1,9)
1159  ELSEIF (i == 11 .AND. (parame(1,6) /= 0.0 .OR. dipole == 1))THEN
1160  write (*,*) ' Specify polarizability [A**3]'
1161  READ (*,*) parame(1,11)
1162  ! parame(1,11) = 0.0 !JG
1163  ELSEIF(i == 11 .AND. (parame(1,7) /= 0.0 .OR. quadrup == 1))THEN
1164  write (*,*) ' Specify polarizability [A**3]'
1165  READ (*,*) parame(1,11)
1166  ENDIF
1167  ENDIF
1168  ENDDO
1169  ENDIF
1170 
1171  !-----------------------------------------------------------------------------
1172  ! write association parameters
1173  !-----------------------------------------------------------------------------
1174 
1175  IF ( nhb_typ(1) /= 0 ) THEN
1176  parame(1,12) = REAL(nhb_typ(1))
1177  ! parame(1,13) = kap_hb(1,1)
1178  no = 0
1179  DO j = 1,nhb_typ(1)
1180  DO k = 1,nhb_typ(1)
1181  parame(1,(14+no)) = eps_hb(1,1,j,k)
1182  no = no + 1
1183  ENDDO
1184  ENDDO
1185  DO j = 1,nhb_typ(1)
1186  parame(1,(14+no)) = nhb_no(1,j)
1187  no = no + 1
1188  ENDDO
1189  ENDIF
1190 
1191 
1192  !-----------------------------------------------------------------------------
1193  ! initialize the array of adjustable parameters
1194  !-----------------------------------------------------------------------------
1195 
1196  DO i = 1, lm_n_par
1197  rel(i) = dummy(i,2)
1198  para(i) = dummy(i,3)
1199  END DO
1200 
1201  !-----------------------------------------------------------------------------
1202  ! weighting for the different data types (currently all =1.0)
1203  !-----------------------------------------------------------------------------
1204 
1205  WRITE (*,*) ' '
1206  WRITE (*,*), 'Input of weighting factors for the different data-sets'
1207  WRITE (*,*) ' '
1208 
1209  IF ( eqa(1) /= 0 ) THEN
1210  WRITE (*,*), 'weights for density data ?'
1211  type_weight(1) = 1.0
1212  write (*,*) type_weight(1)
1213  ! READ (*,*) type_weight(1)
1214  END IF
1215  IF (eqa(2) /= 0) THEN
1216  WRITE (*,*), 'weights for vapor pressure data ?'
1217  type_weight(2) = 1.0
1218  write (*,*) type_weight(2)
1219  ! READ (*,*) type_weight(2)
1220  END IF
1221  IF (eqa(3) /= 0) THEN
1222  WRITE (*,*), 'weights for 2nd virial coefficients ?'
1223  READ (*,*) type_weight(3)
1224  END IF
1225  IF (eqa(4) /= 0) THEN
1226  WRITE (*,*), 'weights for energy data ?'
1227  READ (*,*) type_weight(4)
1228  END IF
1229 
1230  ! --- summarizing weighting of the different data types ----------------------
1231  WRITE (*,*), 'the following weights were specified:'
1232  WRITE (*,*) ' '
1233  IF (eqa(1) == 1) WRITE (*,*), ' density data : ', type_weight(1)
1234  IF (eqa(2) == 1) WRITE (*,*), ' vapor pressure data : ', type_weight(2)
1235  IF (eqa(3) == 1) WRITE (*,*), ' 2nd virial coefficients : ', type_weight(3)
1236  IF (eqa(4) == 1) WRITE (*,*), ' int. res. energy data : ', type_weight(4)
1237 
1238  !-----------------------------------------------------------------------------
1239  ! summarizing the starting values and the scaling of parameters
1240  !-----------------------------------------------------------------------------
1241 
1242  WRITE (*,*) ' '
1243  WRITE(23,*) ' '
1244  WRITE(23,*) ' '
1245  WRITE(23,'(a)') ' PARAMETER-TYPE PARAM.-START SCALING para'
1246  WRITE(23,'(a)') ' --------------------------------------------------------------'
1247  DO i=1,lm_n_par
1248  WRITE(23,'(a,t26,G13.6,t41,G13.6,t56,G13.6)') aa_txt(i),dummy(i,1),rel(i),para(i)
1249  END DO
1250  WRITE(23,*) ' '
1251  WRITE(23,'(a)') ' WEIGHTS'
1252  WRITE(23,'(a)') ' -------'
1253  IF (eqa(1) == 1) WRITE(23,'(A,T26,F5.2)')' density data :',type_weight(1)
1254  IF (eqa(2) == 1) WRITE(23,'(A,T26,F5.2)')' vapor pressure data :',type_weight(2)
1255  IF (eqa(3) == 1) WRITE(23,'(A,T26,F5.2)')' 2nd virial coeff. :',type_weight(3)
1256  IF (eqa(4) == 1) WRITE(23,'(A,T26,F5.2)')' energy data :',type_weight(4)
1257  WRITE(23,*) ' '
1258 
1259 END SUBROUTINE paini
1260 
1261 
1262 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1263 !
1264 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1265 
1266 SUBROUTINE pure_output ( lm_n_par, para )
1267 
1268  USE basic_variables
1269  USE pure_fit_parameters
1270  IMPLICIT NONE
1271 
1272  !-----------------------------------------------------------------------------
1273  INTEGER, INTENT(IN) :: lm_n_par
1274  REAL, INTENT(IN) :: para(:)
1275 
1276  !-----------------------------------------------------------------------------
1277  INTEGER :: i,j,ps
1278  REAL :: devi, devimax, devisum, deviav, rms
1279  REAL :: p_deviav,p_devimax,r_deviav,r_devimax
1280  CHARACTER (LEN=4) :: char_len
1281  !-----------------------------------------------------------------------------
1282 
1283 
1284  !-----------------------------------------------------------------------------
1285  ! writing resulting parameters to screen and to output file
1286  !-----------------------------------------------------------------------------
1287 
1288  WRITE (*,*) ' '
1289  WRITE (*,*) '-----------------------------------------------'
1290  DO i = 1, lm_n_par
1291  WRITE (*,'(a,2x,f12.6)') aa_txt(i), para(i)*rel(i)
1292  END DO
1293  !WRITE (*,'(a,3(2x,f15.5))') ' STD ',(stdval(i), i=1,3)
1294  WRITE (*,*) '-----------------------------------------------'
1295 
1296  WRITE (23,*) ' '
1297  DO i = 1, lm_n_par
1298  WRITE (23,'(T2,A,I2,T26,G16.9)') 'PARAMETER ', i, para(i)*rel(i)
1299  END DO
1300 
1301 
1302  !-----------------------------------------------------------------------------
1303  ! writing resulting parameters in the format of para_input.f90
1304  !-----------------------------------------------------------------------------
1305 
1306  WRITE(23, *) ' '
1307  WRITE(23, *)'-----------------------------------------------'
1308  WRITE(23, *)' mm(i) =',mm(1)
1309  DO i = 1, 11
1310  ps = 0
1311  DO j = 1, lm_n_par
1312  IF (fit_array(j) == i) ps = j
1313  END DO
1314 
1315  IF ( i == 1 .AND. ps >= 1 ) THEN
1316  WRITE(23, '(a,G16.8)') '! parame(i,1) = mm(i)*',para(ps)*rel(ps)
1317  WRITE(23, *)' parame(i,1) =',mm(1)*para(ps)*rel(ps)
1318  ELSE IF ( i == 1 ) THEN
1319  WRITE(23, '(a,G16.8)')'! parame(i,1) = mm(i)*',parame(1,1)
1320  WRITE(23, *)' parame(i,1) =',mm(1)*parame(1,1)
1321  END IF
1322 
1323  IF (i == 2 .AND. ps >= 1) THEN
1324  WRITE(23, *)' parame(i,2) =',para(ps)*rel(ps)
1325  ELSE IF (i == 2) THEN
1326  WRITE(23, *)' parame(i,2) =',parame(1,2)
1327  END IF
1328 
1329  IF (i == 3 .AND. ps >= 1) THEN
1330  WRITE(23, *)' parame(i,3) =',para(ps)*rel(ps)
1331  ELSE IF (i == 3) THEN
1332  WRITE(23, *)' parame(i,3) =',parame(1,3)
1333  END IF
1334 
1335  IF (parame(1,12) == 2.0) THEN
1336  IF (i == 4 .AND. ps >= 1) THEN
1337  WRITE(23, *)' nhb_typ(i) = ',nint(parame(1,12))
1338  WRITE(23, *)' nhb_no(i,1) = ',nint(parame(1,18))
1339  WRITE(23, *)' nhb_no(i,2) = ',nint(parame(1,19))
1340  WRITE(23, *)' eps_hb(i,i,1,2)= ',para(ps)*rel(ps)
1341  WRITE(23, *)' eps_hb(i,i,2,1)= ',para(ps)*rel(ps)
1342  WRITE(23, *)' eps_hb(i,i,1,1)= 0.0'
1343  WRITE(23, *)' eps_hb(i,i,2,2)= 0.0'
1344  ELSE IF (i == 4) THEN
1345  WRITE(23, *)' nhb_typ(i) = ',nint(parame(1,12))
1346  WRITE(23, *)' nhb_no(i,1) = ',nint(parame(1,18))
1347  WRITE(23, *)' nhb_no(i,2) = ',nint(parame(1,19))
1348  WRITE(23, *)' eps_hb(i,i,1,2)= ',parame(1,15)
1349  WRITE(23, *)' eps_hb(i,i,2,1)= ',parame(1,15)
1350  WRITE(23, *)' eps_hb(i,i,1,1)= 0.0'
1351  WRITE(23, *)' eps_hb(i,i,2,2)= 0.0'
1352  END IF
1353  IF (i == 5 .AND. ps >= 1) THEN
1354  WRITE(23, *)' kap_hb(i,i) = ',para(ps)*rel(ps)
1355  ELSE IF (i == 5) THEN
1356  WRITE(23, *)' kap_hb(i,i) = ',parame(1,13)
1357  END IF
1358  ELSE IF ( nint(parame(1,12)) == 1 ) THEN
1359  IF (i == 4 .AND. ps >= 1) THEN
1360  WRITE(23, *)' nhb_typ(i) = ',nint(parame(1,12))
1361  WRITE(23, *)' eps_hb(i,i,1,1)= ',para(ps)*rel(ps)
1362  ELSE IF (i == 4) THEN
1363  WRITE(23, *)' nhb_typ(i) = ',nint(parame(1,12))
1364  WRITE(23, *)' eps_hb(i,i,1,1)= ',parame(1,14)
1365  END IF
1366  IF (i == 5 .AND. ps >= 1) THEN
1367  WRITE(23, *)' kap_hb(i,i) = ',para(ps)*rel(ps)
1368  ELSE IF (i == 5) THEN
1369  WRITE(23, *)' kap_hb(i,i) = ',parame(1,13)
1370  END IF
1371  END IF
1372 
1373  IF (i == 6 .AND. ps >= 1) THEN
1374  WRITE(23, *)' parame(i,6) = ',para(ps)*rel(ps)
1375  ELSE IF (i == 6 .AND. parame(1,6) /= 0.0) THEN
1376  WRITE(23, *)' parame(i,6) = ',parame(1,6)
1377  END IF
1378  IF (i == 7 .AND. ps >= 1) THEN
1379  WRITE(23, *)' parame(i,8) = ',para(ps)*rel(ps)
1380  ELSE IF (i == 7 .AND. parame(1,8) /= 0.0) THEN
1381  WRITE(23, *)' parame(i,8) = ',parame(1,8)
1382  END IF
1383 
1384  IF (i == 8 .AND. ps >= 1) THEN
1385  WRITE(23, *)' parame(i,7) = ',para(ps)*rel(ps)
1386  ELSE IF (i == 8 .AND. parame(1,7) /= 0.0) THEN
1387  WRITE(23, *)' parame(i,7) = ',parame(1,7)
1388  END IF
1389  IF (i == 9 .AND. ps >= 1) THEN
1390  WRITE(23, *)' parame(i,9) = ',para(ps)*rel(ps)
1391  ELSE IF (i == 9 .AND. parame(1,9) /= 0.0) THEN
1392  WRITE(23, *)' parame(i,9) = ',parame(1,9)
1393  END IF
1394  IF (i == 11 .AND. ps >= 1) THEN
1395  WRITE(23, *)' parame(i,11) = ',para(ps)*rel(ps)
1396  ELSE IF (i == 11 .AND. parame(1,11) /= 0.0) THEN
1397  WRITE(23, *)' parame(i,11) = ',parame(1,11)
1398  END IF
1399  END DO
1400  WRITE(23, *) '-----------------------------------------------'
1401 
1402  !-----------------------------------------------------------------------------
1403  ! analyizing result
1404  !-----------------------------------------------------------------------------
1405 
1406  !WRITE (*,*) ' '
1407  !WRITE (*,'(T2, A, I1, 2X, A, 1PD12.5)') 'Convergence= ', ic, 'SUMERR= ', sse
1408  !WRITE(23, *)
1409  !WRITE(23, '(A, I1, 2X, A, 1PD12.5)') 'Convergence= ', ic, 'SUMERR= ', sse
1410  !WRITE (*,*) ' '
1411  !WRITE (*,*), 'Normalized correlation of parameters:'
1412  !WRITE (*,*) ' '
1413  !WRITE(23,*)
1414  !WRITE(23,*) 'Normalized correlation of parameters:'
1415  !FORM='(T2, 2X, '//mul//'(3X, A, 5X))'
1416  !WRITE (*,FORM) (namp(i), i=1,lm_n_par)
1417  !WRITE(23, FORM) (namp(i), i=1,lm_n_par)
1418  !FORM='(T2, I2, '//mul//'(2X, D10.4)/)'
1419  !DO i=1,lm_n_par
1420  ! WRITE (*,FORM) i, (corval(i,k), k=1,lm_n_par)
1421  ! WRITE(23, FORM) i, (corval(i,k), k=1,lm_n_par)
1422  !END DO
1423  WRITE (*,*) ' '
1424  WRITE(23, *)
1425 
1426 
1427  !-----------------------------------------------------------------------------
1428  ! comparing calculated values to experimental data
1429  !-----------------------------------------------------------------------------
1430 
1431  !-----------------------------------------------------------------------------
1432  ! 1. density data
1433  !-----------------------------------------------------------------------------
1434 
1435  IF ( eqa(1) == 1 ) THEN
1436  WRITE (23,*) '------------- fluid density data --------------'
1437  WRITE(23,'(/A)') 'DATNR AD (%) T/K v_exp/(m3/kg) v_calc/(m3/kg) '
1438  devisum = 0.0
1439  devimax = 0.0
1440  devi = 0.0
1441  rms = 0.0
1442  ! --- write data, calc. RMS-%, AAD-%, and maximum observed deviation -
1443  DO i = 1, nliq
1444  devi = abs(v_cal(i)-vliq(i))/ (sqrt(abs(vliq(i)))*sqrt(abs(v_cal(i))))*100.0
1445  rms = rms + ((v_cal(i)-vliq(i))/vliq(i))**2
1446  IF (devi > devimax) devimax=devi
1447  devisum = devisum + devi
1448  WRITE(23,'(i3,2x,G10.3,3(2x,G12.5))') i, devi, tliq(i), vliq(i), v_cal(i)
1449  END DO
1450  deviav = devisum / REAL(nliq)
1451  rms = sqrt( rms / REAL(nliq) )*100.0
1452  WRITE(23,*) ' '
1453  WRITE(23,'(A)') 'Deviation of calculated to exp. volume data'
1454  WRITE(23,*) ' '
1455  WRITE(23,'(A, t35, F7.3, A)') 'Max. Deviation:',devimax, ' %'
1456  WRITE(23,'(A, t35, F7.3, A)') 'Average Deviation: ',deviav, ' %'
1457  WRITE(23,'(A, t35, F7.3, A)') 'RMS: ',rms, ' %'
1458  WRITE(23,*) ' '
1459  r_devimax = devimax
1460  r_deviav = deviav
1461  END IF
1462 
1463  !-----------------------------------------------------------------------------
1464  ! 2. vapor pressure data
1465  !-----------------------------------------------------------------------------
1466 
1467  IF ( eqa(2) == 1 ) THEN
1468  WRITE (23,*) '------------- vapor pressure data -------------'
1469  WRITE(23,'(/A)') 'DATNR AD (%) T/K PEXP/PA PCAL/PA'
1470  devisum = 0.0
1471  devimax = 0.0
1472  devi = 0.0
1473  rms = 0.0
1474  DO i = 1, nlv
1475  ! --- write data, calc. RMS-%, AAD-%, and maximum observed deviation
1476  devi = abs(plvcal(i)-plv(i))/ (sqrt(abs(plv(i)))*sqrt(abs(plvcal(i))))*100.0
1477  rms = rms + ((plvcal(i)-plv(i))/plv(i))**2
1478  IF (devi > devimax) devimax = devi
1479  devisum = devisum + devi
1480  WRITE (23,'(i3,2x,G10.3,2x,G12.5,2(2x,G14.7))') i, devi, tlv(i), plv(i), plvcal(i)
1481  END DO
1482  deviav = devisum / REAL(nlv)
1483  rms = sqrt( rms / REAL(nlv)) * 100.0
1484  WRITE(23,*) ' '
1485  WRITE(23,'(A)') 'Deviation of calculated to exp. vapor pressure data'
1486  WRITE(23,*) ' '
1487  WRITE(23, '(3(A, t35, F7.3, A/))') 'Max. Deviation:',devimax, ' %', &
1488  'Average Deviation: ',deviav, ' %', 'RMS: ',rms, ' %'
1489  WRITE(23,*) ' '
1490  p_devimax = devimax
1491  p_deviav = deviav
1492  END IF
1493 
1494  !-----------------------------------------------------------------------------
1495  ! 3.: 2nd virial coefficient data
1496  !-----------------------------------------------------------------------------
1497 
1498  IF ( eqa(3) == 1 ) THEN
1499  WRITE (23,*) '---------------2nd virial coeff. --------------'
1500  WRITE (23,'(/A,A)') 'DATNR AD (%) T/K B_exp B_cal'
1501  devisum = 0.0
1502  devimax = 0.0
1503  devi = 0.0
1504  rms = 0.0
1505  ! --- write data, calc. RMS-%, AAD-%, and maximum observed deviation -
1506  DO i = 1, n_vir
1507  devi = abs(b_cal(i)-b_vir(i))/ (sqrt(abs(b_vir(i)))*sqrt(abs(b_cal(i))))*100.0
1508  rms = rms + ((b_cal(i)-b_vir(i))/b_vir(i))**2
1509  IF (devi > devimax) devimax = devi
1510  devisum = devisum + devi
1511  WRITE(23,'(i3,2x,4(2x,G12.5))') i, devi, t_vir(i), b_vir(i), b_cal(i)
1512  END DO
1513  deviav = devisum / REAL(n_vir)
1514  rms = sqrt( rms / REAL(n_vir) ) * 100.0
1515  WRITE(23,*) ' '
1516  WRITE(23,'(A)') 'Deviation of calculated to exp. 2nd virial coeff. data'
1517  WRITE(23,*) ' '
1518  WRITE(23, '(3(A, t35, F6.2, A/))') 'Max. Deviation:',devimax, ' %', &
1519  'Average Deviation: ',deviav, ' %', 'RMS: ',rms, ' %'
1520  WRITE(23,*) ' '
1521  END IF
1522 
1523 
1524  !-----------------------------------------------------------------------------
1525  ! 4.: internal energy data
1526  !-----------------------------------------------------------------------------
1527 
1528  IF ( eqa(4) == 1 ) THEN
1529  WRITE (23,*) '---------------- energy data ------------------'
1530  WRITE (23,'(/A,A)') 'DATNR AD (%) U_EXP U_CAL '
1531  devisum = 0.0
1532  devimax= 0.0
1533  devi = 0.0
1534  rms = 0.0
1535  ! --- write data, calc. RMS-%, AAD-%, and maximum observed deviation -
1536  DO i = 1, n_en
1537  devi = abs(u_cal(i)-u_en(i))/(sqrt(abs(u_en(i)))*sqrt(abs(u_cal(i))))*100.0
1538  rms = rms + ((u_cal(i)-u_en(i))/u_en(i))**2
1539  IF (devi > devimax) devimax=devi
1540  devisum = devisum + devi
1541  WRITE(23,'(i3,2x,4(2x,G12.5))') i, devi, t_en(i), u_en(i), u_cal(i)
1542  END DO
1543  deviav = devisum / REAL(n_en)
1544  rms = sqrt( rms / REAL(n_en) ) * 100.0
1545  WRITE(23,*) ' '
1546  WRITE(23,'(A)') 'Deviation of calculated to exp. internal energy data'
1547  WRITE(23,*) ' '
1548  WRITE(23, '(3(A, t35, F6.2, A/))') 'Max. Deviation:',devimax, ' %', &
1549  'Average Deviation: ',deviav, ' %', 'RMS: ',rms, ' %'
1550  WRITE(23,*) ' '
1551  END IF
1552 
1553  WRITE(23, *)' '
1554  WRITE (char_len,'(I3)') lm_n_par+5
1555  WRITE(23,'(a,'//char_len//'G15.8)') pure_fit_file, mm(1), &
1556  (para(i)*rel(i),i=1,lm_n_par),r_deviav,r_devimax,p_deviav,p_devimax
1557  WRITE(23, *)' '
1558 
1559 
1560 END SUBROUTINE pure_output
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120
basic_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:29