MoDeNa  1.0
Software framework facilitating sequential multi-scale modelling
eos_critical_renorm_mix.f90
1 Module rgt_phase_space_cell
2 
3  use parameters, only: nc
4  implicit none
5  save
6 
7  real, dimension(nc) :: dfdx
8  REAL, DIMENSION(300) :: x1a
9  REAL, DIMENSION(300) :: x2a
10  REAL, DIMENSION(300,300) :: ya, y1a, y2a, y12a
11  REAL, DIMENSION(300,300,4,4) :: c_bicub
12 
13 End Module rgt_phase_space_cell
14 
15 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
16 ! SUBROUTINE F_EOS_RN_PHASE_SPACE_CELL
17 !
18 ! White's recursive procedure is:
19 !
20 ! f = fid + f0 + delta_f1 + delta_f2 + delta_f3 + ...
21 !
22 ! where f is the Helmholtz energy per unit volume, and where f0 is the
23 ! classical residual repulsive part, f0=f_res,rep.
24 !
25 ! PC-SAFT can be considered for the classical residual part, with
26 !
27 ! f0 = f_PCSAFT
28 !
29 ! It is
30 !
31 ! INT_0^rho( exp(-Kn^-1*(fns(rho+x)-2fns(rho)+fns(rho-x))) )
32 ! delta_fn= -Kn*ln----------------------------------------------------------
33 ! INT_0^rho( exp(-Kn^-1*(fnl(rho+x)-2fnl(rho)+fnl(rho-x))) )
34 !
35 ! where x in the integrals INT runs from x=0 to x=rho_x=rho. The counter
36 ! kk in belows code takes care of the integrals INT from 0 to rho.
37 ! An outer loop is established (with index k) in order to calculate f for
38 ! the whole density range. First, the f0 is calculated and approximated
39 ! by a cubic spline. Then delta_fn is calculated for 0<=rho<=rhomax. New
40 ! spline parameters are determined after evey iteration n.
41 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
42 
43 SUBROUTINE f_eos_rn_phase_space_cell
44 
45  USE parameters
46  USE eos_variables
47  USE eos_constants
48  USE eos_numerical, only: f_numerical
49  USE fitting_rgt_parameters
50  USE rgt_phase_space_cell
51  IMPLICIT NONE
52 
53  !-----------------------------------------------------------------------------
54  INTEGER :: n, i, j, m
55  INTEGER :: k1, k2, kk1, kk2
56  INTEGER :: k05max
57  INTEGER :: stepno, niter, pos1, pos2, step1, step2
58  REAL :: rho0, rhomax1, rhomax2, kn(10,0:300,0:300), ll
59  REAL :: xsav1, xsav2
60  REAL :: int1_l, int1_s, int3_l, int3_s
61  REAL :: int4_l, int4_s
62  REAL :: integral_l, integral_s
63  REAL :: del_f(8,0:300,0:300)
64  REAL :: phi_crit, rhoi(nc)
65  REAL :: alph(0:8,0:300,0:300), w_ratio
66  REAL :: rhovec1(0:300), rhovec2(0:300), dzp1, dzp2
67  REAL :: fres_v, sig_m2, m_mean
68  REAL :: chapm, order_chap, fid
69  REAL :: mfrac(2)
70  REAL :: rhot, alphn(0:300,0:300), combi(0:300,0:300)
71  REAL :: combi2(10,0:300,0:300), fact, i1
72  REAL :: rho1, rho2, fdr1, fdr2, fdr11, fdr12, fdr22
73  REAL :: fdr111, fdr112, fdr122, fdr222, f_drho, f_drho2
74  REAL :: f_drho3, f_drho4, z3tsav
75  REAL :: int3_lx, int3_sx, int_o_l, int_o_s
76  REAL :: densav1
77 
78  INTEGER, SAVE :: scan = 0
79  REAL, SAVE :: tempsav = 0.0
80  !-----------------------------------------------------------------------------
81 
82  CALL perturbation_parameter
83 
84  ll = 0.0
85  phi_crit = 0.0
86  sig_m2 = 0.0
87  m_mean = 0.0
88  chapm = 0.0
89  order_chap = 0.0
90  DO i = 1,ncomp
91  IF (lli(i) == 0.0 .OR. phi_criti(i) == 0.0 .OR. chap(i) == 0.0) stop
92  ll = ll + x(i)*lli(i)**3
93  phi_crit = phi_crit + x(i)*phi_criti(i)
94  sig_m2 = sig_m2 + x(i)*mseg(i)*sig_ij(i,i)**2
95  m_mean = m_mean + x(i)*mseg(i)
96  chapm = chapm + x(i)*chap(i)
97  DO j = 1,ncomp
98  order_chap = order_chap + x(i)*x(j)* mseg(i)*mseg(j) &
99  *sig_ij(i,j)**3 * uij(i,j)/t *(chap(i)*chap(j))**0.5
100  END DO
101  END DO
102  ll = ll**(1.0/3.0)
103  sig_m2 = sig_m2 / m_mean
104 
105  stepno = 200 ! for stepno > 200, expand the array dimension of DFT_MODULE
106  niter = 4
107  w_ratio = 2.0
108 
109  xsav1 = x(1)
110  xsav2 = x(2)
111  z3tsav = z3t
112 
113  rho0 = eta/z3t
114  rho1 = x(1)*rho0
115  rho2 = x(2)*rho0
116  ! rhomax = 0.0
117  ! DO i = 1,ncomp
118  ! !rhomax = rhomax + x(i)*mseg(i)*sig_ij(i,i)**3
119  ! rhomax = rhomax + x(i)*mseg(i)*dhs(i)**3 !jg
120  ! ENDDO
121  ! rhomax = SQRT(2.0)/rhomax
122  rhomax1 = mseg(1)*dhs(1)**3 !jg
123  rhomax2 = mseg(2)*dhs(2)**3 !jg
124  rhomax1 = sqrt(2.0)/rhomax1
125  rhomax2 = sqrt(2.0)/rhomax2
126 
127 
128  densav1 = eta
129 
130  IF (tempsav == t .AND. scan == 1) GO TO 5
131 
132  tempsav = t
133  !-----------------------------------------------------------------------------
134 
135 
136  dzp1 = rhomax1/REAL(stepno)
137  dzp2 = rhomax2/REAL(stepno)
138  DO k2 = 0,stepno
139  DO k1 = 0,stepno
140  rhovec1(k1) = REAL(k1)*dzp1
141  rhovec2(k2) = REAL(k2)*dzp2
142  rhot = rhovec1(k1)+rhovec2(k2)
143  x(1) = rhovec1(k1)/rhot
144  x(2) = rhovec2(k2)/rhot
145  alph(0,k1,k2) = 0.0
146  eta = rhot*pi/6.0*(x(1)*mseg(1)*dhs(1)**3 &
147  +x(2)*mseg(2)*dhs(2)**3) !*z3t
148  IF (eta > 0.7) eta = 0.8-0.1*exp((0.7-eta)*10.0)
149  rhoi(1) = rhovec1(k1)
150  rhoi(2) = rhovec2(k2)
151 
152  m_mean = x(1)*mseg(1) + x(2)*mseg(2)
153  ! v_mean = x(1)*mseg(1)*sig_ij(1,1)**3 +x(2)*mseg(2)*sig_ij(2,2)**3
154  ! vfrac(1) = mseg(1)*sig_ij(1,1)**3 / v_mean
155  ! vfrac(2) = mseg(2)*sig_ij(2,2)**3 / v_mean
156  mfrac(1) = mseg(1) / m_mean
157  mfrac(2) = mseg(2) / m_mean
158 
159  combi(k1,k2) = 0.0
160  DO i = 1,ncomp
161  DO j = 1,ncomp
162  combi(k1,k2) = combi(k1,k2)+16.0/9.0*pi*( &
163  rhoi(i)*rhoi(j)*mseg(i)*mseg(j)*sig_ij(i,j)**3 &
164  * uij(i,j)/t *(chap(i)+chap(j))*.5 )
165  !combi(k1,k2) = combi(k1,k2)+16.0/9.0*PI*( &
166  ! rhoi(i)*rhoi(j)*mseg(i)*mseg(j)*sig_ij(i,j)**3 &
167  ! * uij(i,j)/t *(chap(i)*mseg(i)+chap(j)*mseg(j))*0.5 ) &
168  ! / m_mean
169  END DO
170  END DO
171 
172  i1 = 0.0
173  DO m = 0,6
174  i1 = i1 + apar(m)*eta**REAL(m)
175  END DO
176  order1 = 0.0
177  DO i = 1,ncomp
178  DO j = 1,ncomp
179  order1 = order1 + x(i)*x(j)* mseg(i)*mseg(j) &
180  *sig_ij(i,j)**3 * uij(i,j)/t *(chap(i)+chap(j))*.5
181  END DO
182  END DO
183  ! combi(k1,k2) = + 2.0*PI*rhot*rhot*I1*order1
184  ! write (*,*) combi(k1,k2),+ 2.0*PI*rhot*rhot*I1*order1
185 
186  DO n = 1,niter
187  combi2(n,k1,k2) = 0.0
188  DO i = 1,ncomp
189  DO j = 1,ncomp
190  combi2(n,k1,k2) = combi2(n,k1,k2)+16.0/9.0*pi*9.0/7.0*( &
191  rhoi(i)*rhoi(j)*mseg(i)*mseg(j)*sig_ij(i,j)**2 &
192  *(phi_criti(i)*mfrac(i) &
193  +phi_criti(j)*mfrac(j))*0.5 *sig_ij(i,j)**3 &
194  * uij(i,j)/t *(chap(i)+chap(j))*0.5 &
195  /(mfrac(i)**0.33333333*lli(i)) /(mfrac(j)**0.33333333*lli(j)) &
196  )/ ( w_ratio**REAL(2*n) * 2.0 )
197  ! combi2(n,k1,k2) = combi2(n,k1,k2)+16.0/9.0*PI*9.0/7.0*(
198  ! & rhoi(i)*rhoi(j)*mseg(i)*mseg(j)*sig_ij(i,j)**2
199  ! & *(phi_criti(i)/mseg(i) + phi_criti(j)/mseg(j))*0.5
200  ! & *sig_ij(i,j)**3 * uij(i,j)/t *(chap(i)*mseg(i)+chap(j)*mseg(j))*.5
201  ! & /(mfrac(i)**0.33333333*LLi(i))
202  ! & /(mfrac(j)**0.33333333*LLi(j))
203  ! & )/ ( w_ratio**REAL(2*n) * 2.0 )
204  END DO
205  END DO
206  ll = ( x(1)*mfrac(1)*lli(1)**3 &
207  + x(2)*mfrac(2)*lli(2)**3 )**(1.0/3.0)
208  ! LL = 0.0
209  ! DO i = 1,ncomp
210  ! DO j = 1,ncomp
211  ! LL = LL +x(i)*x(j)*mfrac(i)*mfrac(j)*LLi(i)**1.5*LLi(j)**1.5
212  ! ENDDO
213  ! ENDDO
214  ! LL = LL**(1.0/3.0)
215  fact = x(1)*sig_ij(1,1)**2 *phi_criti(1) &
216  +x(2)*sig_ij(2,2)**2 *phi_criti(2)
217  ! combi2(n,k1,k2) = combi(k1,k2)*9.0/7.0*fact /LL/LL / ( 2.0**REAL(2*n) * 2.0 )
218  kn(n,k1,k2) = 1.0/ll**3 / w_ratio**REAL(3*n)
219  END DO
220  i1 = 0.0
221  DO m = 0,6
222  i1 = i1 + apar(m)*eta**REAL(m)
223  END DO
224  order1 = 0.0
225  DO i = 1,ncomp
226  DO j = 1,ncomp
227  order1 = order1 + x(i)*x(j)* mseg(i)*mseg(j) &
228  *sig_ij(i,j)**5 * uij(i,j)/t *(chap(i)+chap(j))*.5 &
229  *(phi_criti(i)*phi_criti(j))**0.5 /lli(i)/lli(j)
230  END DO
231  END DO
232  DO n = 1,niter
233  ! combi2(n,k1,k2) = + 2.0*PI*rhot*rhot*I1*order1 *9.0/7.0 / ( w_ratio**REAL(2*n) * 2.0 )
234  END DO
235 
236 
237  fid = 0.0
238  IF (rhot /= 0.0) THEN
239  DO i = 1,ncomp
240  ! debroglie(i) = 6.62606896d-34 *1d10 & ! in units Angstrom
241  ! * SQRT( 1.0 / (2.0*PI *1.0 /6.022045d23/1000.0*KBOL*T) )
242  IF (rhoi(i) > 0.0) fid = fid + rhoi(i)*(log(rhoi(i))-1.0)
243  END DO
244  CALL f_numerical
245  alph(0,k1,k2) = fres*rhot + fid
246  alphn(k1,k2) = alph(0,k1,k2)
247  ! ya(k1+1,k2+1) = alphn(k1,k2) ! this line can be deleted after debugging
248  END IF
249  ! x1a(k1+1) = rhovec1(k1) ! this line can be deleted after debugging
250  END DO
251  ! x2a(k2+1) = rhovec2(k2) ! this line can be deleted after debugging
252  END DO
253  ! CALL SPLINE_PARA (dzp,alphn,utri,stepno)
254  ! CALL SPLINE_COEFF(beta,gamma,delta,dzp,alphn,utri,stepno)
255  scan = 1
256 
257  x(1) = xsav1
258  x(2) = xsav2
259  z3t = z3tsav
260 
261  !-----------------------------------------------------------------------------
262  ! alpha_l = 16.0/9.0*PI*order_chap !* chapm
263  !c alpha_s = alpha_l*phi_crit* 9.0/7.0*sig_m2 /LL/LL
264  ! alpha_s = phi_crit* 9.0/7.0*sig_m2 /LL/LL
265 
266  ! n = 1 ! for debugging only
267  ! GOTO 33
268  DO n = 1,niter
269 
270  ! Kn(n) = 1.0/LL**3 / 2.0**REAL(3*n)
271  k05max = stepno/2
272 
273 
274  ! calculate for rhovec2 = 0
275  DO k1 = 1,stepno*3/4-1 !!!!!
276  int_o_l = 1.0
277  int_o_s = 1.0
278  integral_l = 0.0
279  integral_s = 0.0
280  step1 = k1 !step1: step# for integrat.(index kk1)
281  IF(k1 > k05max) step1 = stepno-k1
282  DO kk1 = 0,step1-1
283  CALL integr(n,k1,0,kk1,0,kn,combi,combi2, alph,int1_l,int1_s)
284  ! int1_lv(kk1) = int1_l
285  ! int1_sv(kk1) = int1_s
286  ! integral_l = integral_l + dzp1 *(int1_l+int_o_l)/2.0
287  ! integral_s = integral_s + dzp1 *(int1_s+int_o_s)/2.0
288  integral_l = integral_l + dzp1 *int1_l
289  integral_s = integral_s + dzp1 *int1_s
290  int_o_l = int1_l
291  int_o_s = int1_s
292  END DO
293  del_f(n,k1,0) = 0.0
294  IF (integral_s /= 0.0.AND. integral_l /= 0.0) THEN
295  del_f(n,k1,0) = -kn(n,k1,0)*log(integral_s/integral_l)
296  END IF
297  END DO ! enddo k1
298 
299  ! calculate for rhovec1 = 0
300  DO k2 = 1,stepno*3/4-1 !!!!!
301  int_o_l = 1.0
302  int_o_s = 1.0
303  integral_l = 0.0
304  integral_s = 0.0
305  step2 = k2 !step2: step# for integrat.(index kk2)
306  IF(k2 > k05max) step2 = stepno-k2
307  DO kk2 = 0,step2-1
308  CALL integr(n,0,k2,0,kk2,kn,combi,combi2, alph,int3_l,int3_s)
309  ! int3_lv(kk2) = int3_l
310  ! int3_sv(kk2) = int3_s
311  ! integral_l = integral_l + dzp2 *(int3_l+int_o_l)/2.0
312  ! integral_s = integral_s + dzp2 *(int3_s+int_o_s)/2.0
313  integral_l = integral_l + dzp2 *int3_l
314  integral_s = integral_s + dzp2 *int3_s
315  int_o_l = int3_l
316  int_o_s = int3_s
317  END DO
318  del_f(n,0,k2) = 0.0
319  IF (integral_s /= 0.0.AND. integral_l /= 0.0) THEN
320  del_f(n,0,k2) = -kn(n,0,k2)*log(integral_s/integral_l)
321  END IF
322  END DO ! enddo k2
323 
324 
325 
326  ! DO k = 1,stepno-1 !!!!!
327  DO k2 = 0,stepno*3/4-1 !!!!!
328  WRITE (*,*) k2,n
329  DO k1 = 0,stepno*3/4-1 !!!!!
330  int3_lx = 1.0
331  int3_sx = 1.0
332  ! hx(1) = 0.0
333  ! hyl(1) = int3_l
334  ! hys(1) = int3_s
335  ! nn = 1
336  integral_l = 0.0
337  integral_s = 0.0
338  step1 = k1 !step1: step# for integrat.(index kk1)
339  step2 = k2 !step2: step# for integrat.(index kk2)
340  IF(k1 > k05max) step1 = stepno-k1
341  IF(k2 > k05max) step2 = stepno-k2
342 
343  DO kk2 = 0,step2-1
344  DO kk1 = 0,step1-1
345 
346  ! write (*,*)rhovec1(k1)+rhovec2(k2),k1,k2
347  ! CALL INTEGR(n,k1,k2,kk1-1,kk2-1,Kn,combi,combi2,alph,int1_l,int1_s)
348  ! CALL INTEGR(n,k1,k2,kk1,kk2-1,Kn,combi,combi2,alph,int2_l,int2_s)
349  ! CALL INTEGR(n,k1,k2,kk1-1,kk2,Kn,combi,combi2,alph,int3_l,int3_s)
350  CALL integr(n,k1,k2,kk1,kk2,kn,combi,combi2, alph,int4_l,int4_s)
351  ! if (k2.EQ.5.AND.kk1.eq.1)write (*,*) int3_l,int4_l,k1,k2
352  ! if (k2.EQ.5.AND.kk1.eq.1)write (*,*) int1_l,int2_l,kk1,kk2
353  ! write (*,*) 'pp',int3_l,int3_lx,k1,k2,kk1,kk2
354  ! if (k2.EQ.5.AND.kk1.eq.1) read (*,*)
355 
356  ! integral_l = integral_l+dzp1*dzp2*(int4_l+int3_l+int2_l+int1_l)/4.0
357  ! integral_s = integral_s+dzp1*dzp2*(int4_s+int3_s+int2_l+int1_l)/4.0
358  integral_l = integral_l+dzp1*dzp2*int4_l
359  integral_s = integral_s+dzp1*dzp2*int4_s
360  int3_lx = int4_l
361  int3_sx = int4_s
362 
363  IF (int4_l < 1.e-9.AND.int4_s < 1.e-9) THEN
364  ! write (*,*) kk1,kk2,int4_l,int4_s
365  GO TO 15 ! end loop, because no significant contribution is expected
366  END IF
367 
368  END DO ! enddo kk1 (integration over density)
369 15 CONTINUE
370  END DO ! enddo kk2 (integration over density)
371 
372 
373  IF (k1 /= 0.AND.k2 /= 0) del_f(n,k1,k2) = 0.0
374  IF (integral_s /= 0.0.AND. integral_l /= 0.0) THEN
375  del_f(n,k1,k2) = -kn(n,k1,k2)*log(integral_s/integral_l)
376  END IF
377 
378  END DO ! enddo k1
379  END DO ! enddo k2
380 
381  ! 33 CONTINUE ! for debugging only
382 
383  del_f(n,0,0) = 0.0
384  DO k2 = 0,stepno
385  DO k1 = 0,stepno
386  alph(n,k1,k2) = alph(n-1,k1,k2)+del_f(n,k1,k2)
387  alphn(k1,k2) = alph(n,k1,k2)
388  ! subract the ideal gas part in order to have less non-linearity in the low
389  ! density region. The ideal gas part is then added to the final value below.
390  IF(n == niter) THEN
391  fid = 0.0
392  rhoi(1) = rhovec1(k1)
393  rhoi(2) = rhovec2(k2)
394  DO i = 1,ncomp
395  IF (rhoi(i) > 0.0) fid = fid +rhoi(i)*(log(rhoi(i))-1.0)
396  END DO
397  alphn(k1,k2) = alphn(k1,k2)-fid
398  ya(k1+1,k2+1) = alphn(k1,k2)
399  END IF
400  x1a(k1+1) = rhovec1(k1)
401  END DO
402  x2a(k2+1) = rhovec2(k2)
403  END DO
404  ! GOTO 34
405 
406  ! DO k2 = 1,stepno+1
407  ! DO k1 = 1,stepno+1
408  ! ya(k1,k2)
409  ! y1a(k1,k2)
410  ! y2a(k1,k2)
411  ! y12a(k1,k2)
412  ! DO m = 1,4
413  ! DO n = 1,4
414  ! c_bicub(k1,k2,m,n)
415  ! ENDDO
416  ! ENDDO
417  ! x1a(k1)
418  ! ENDDO
419  ! x2a(k2)
420  ! ENDDO
421 
422 
423  END DO ! loop of n cycles
424  ! write (71,*) ' '
425 
426 
427  ! k2 = 0
428  ! DO k1 = 0,stepno/2,1
429  ! write (*,*) k1,k2,alph(0,k1,k2),alph(niter,k1,k2) !,del_f(1,k1,k2)
430  ! ENDDO
431  ! c ENDDO
432  ! read (*,*)
433 
434  ! 34 CONTINUE
435  CALL bicub_derivative2( stepno+1, stepno+1 )
436  CALL bicub_c2 ( stepno+1, stepno+1 )
437 
438 
439 5 CONTINUE
440 
441 
442  pos1 = int(rho1/rhomax1*stepno) + 1 ! plus 1 because the (0:stepno)array is mapped onto a (1:stepno+1)array
443  pos2 = int(rho2/rhomax2*stepno) + 1 ! plus 1 because the (0:stepno)array is mapped onto a (1:stepno+1)array
444  CALL bi_cub_spline_crt (rho1,rho2,fres_v,fdr1,fdr2,fdr11,fdr12,fdr22, &
445  fdr111,fdr112,fdr122,fdr222,stepno+1,stepno+1,pos1,pos2)
446  fres = fres_v/rho0
447  f_drho = x(1)*fdr1 + x(2)*fdr2 ! derivative of Helholtz energy density to rho, d(f'*rho)/d_rho
448 
449  f_drho2 = x(1)*x(1)*fdr11 + 2.0*x(1)*x(2)*fdr12 +x(2)*x(2)*fdr22
450  f_drho3 = x(1)**3 *fdr111 + 3.0*x(1)*x(1)*x(2)*fdr112 &
451  + 3.0*x(1)*x(2)*x(2)*fdr122 + x(2)**3 *fdr222
452  f_drho4 = 0.0 ! not calculated !!!
453 
454  pges = ( f_drho*rho0 - fres_v ) *(kbol*t)/1.e-30 &
455  + rho0 * (kbol*t) / 1.e-30 ! ideal gas contribution
456  pgesdz = ( f_drho2*rho0 ) *(kbol*t)/1.e-30 /z3t &
457  + 1.0/z3t*(kbol*t)/1.e-30 ! ideal gas contribution
458  pgesd2 = ( f_drho3*rho0 + f_drho2 ) *(kbol*t)/1.e-30 /z3t /z3t
459  pgesd3 = ( f_drho4*rho0 + 2.0*f_drho3 ) *(kbol*t)/1.e-30 /z3t /z3t /z3t
460 
461  dfdx(1) = fdr1
462  dfdx(2) = fdr2
463  ! write (*,*) 'chem.p.',fdr1,fdr2,rho1,rho2
464 
465  eta = densav1
466  ! write (*,*) 'fp',fres,pges,pgesdz,pgesd2,pgesd3,eta,rho0*z3t
467  ! read (*,*)
468 
469 END SUBROUTINE f_eos_rn_phase_space_cell
470 
471 
472 
473 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
474 !
475 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
476 
477 SUBROUTINE integr(n,k1,k2,kk1,kk2,kn,combi,combi2, alph,int_l,int_s)
478 
479  IMPLICIT NONE
480 
481  !-----------------------------------------------------------------------------
482  INTEGER, INTENT(IN) :: n
483  INTEGER, INTENT(IN) :: k1
484  INTEGER, INTENT(IN) :: k2
485  INTEGER, INTENT(IN) :: kk1
486  INTEGER, INTENT(IN) :: kk2
487  REAL, INTENT(IN) :: kn(10,0:300,0:300)
488  REAL, INTENT(IN) :: combi(0:300,0:300)
489  REAL, INTENT(IN) :: combi2(10,0:300,0:300)
490  REAL, INTENT(IN) :: alph(0:8,0:300,0:300)
491  REAL, INTENT(OUT) :: int_l
492  REAL, INTENT(OUT) :: int_s
493 
494 
495  REAL :: alp_c,alph_lc,alph_sc,alp_l,alph_ll,alph_sl, &
496  alp_r,alph_lr,alph_sr,gn_l,gn_s
497  !-----------------------------------------------------------------------------
498 
499 
500  alp_c = alph(n-1,k1,k2)
501  alph_lc =alp_c + combi(k1,k2)
502  alph_sc =alp_c + combi2(n,k1,k2)
503  ! write (*,*) 'c',alph(n-1,k1,k2),alph(n-1,k1,k2-kk2) ! ,combi2(n,k1,k2)
504 
505  alp_l = alph(n-1,k1-kk1,k2-kk2)
506  alph_ll = alp_l + combi(k1-kk1,k2-kk2)
507  alph_sl = alp_l + combi2(n,k1-kk1,k2-kk2)
508 
509  alp_r = alph(n-1,k1+kk1,k2+kk2)
510  alph_lr = alp_r + combi(k1+kk1,k2+kk2)
511  alph_sr = alp_r + combi2(n,k1+kk1,k2+kk2)
512 
513  gn_l = (alph_lr+alph_ll)/2.0 - alph_lc
514  gn_s = (alph_sr+alph_sl)/2.0 - alph_sc
515  IF (gn_l < 0.0) gn_l = 0.0
516  IF (gn_s < 0.0) gn_s = 0.0
517  int_l = 0.0
518  int_s = 0.0
519  IF( -gn_l/kn(n,k1,k2) > -300.0.AND.-gn_l/kn(n,k1,k2) < 300.0) int_l = exp( -gn_l/kn(n,k1,k2) )
520  IF( -gn_s/kn(n,k1,k2) > -300.0.AND.-gn_s/kn(n,k1,k2) < 300.0) int_s = exp( -gn_s/kn(n,k1,k2) )
521 
522 END SUBROUTINE integr
523 
524 
525 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
526 !
527 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
528 
529 SUBROUTINE bicub_derivative2 ( i_max, k_max )
530 
531  USE rgt_phase_space_cell
532  IMPLICIT NONE
533 
534  !-----------------------------------------------------------------------------
535  INTEGER, INTENT(IN) :: i_max
536  INTEGER, INTENT(IN) :: k_max
537  !-----------------------------------------------------------------------------
538  INTEGER :: i,k, steps,lm
539  REAL :: d1,d2,gam1,gam2, &
540  d2f_dx1_2(301,301),d2f_dx2_2(301,301)
541  REAL :: yspl(301),xspl(301),y1(301),y2(301),yp1,ypn
542  !-----------------------------------------------------------------------------
543 
544  d1 = x1a(2)-x1a(1)
545  d2 = x2a(2)-x2a(1)
546 
547  gam1 = sqrt(1.0+d2*d2/d1/d1)/(1.0+d2*d2/d1/d1)
548  gam2 = sqrt(1.0+d1*d1/d2/d2)/(1.0+d1*d1/d2/d2)
549 
550  DO i = 2,i_max-1
551  DO k = 2,k_max-1
552  y1a(i,k) = (ya(i+1,k)-ya(i-1,k))/(x1a(i+1)-x1a(i-1))
553  y2a(i,k) = (ya(i,k+1)-ya(i,k-1))/(x2a(k+1)-x2a(k-1))
554  y12a(i,k) = (ya(i+1,k+1)-ya(i+1,k-1)-ya(i-1,k+1)+ya(i-1,k-1)) &
555  /((x1a(i+1)-x1a(i-1))*(x2a(k+1)-x2a(k-1)))
556  END DO
557  END DO
558 
559  i = 1
560  DO k = 2,k_max-1
561  y1a(i,k) = (-ya(i+2,k)+4.0*ya(i+1,k)-3.0*ya(i,k)) /(x1a(i+2)-x1a(i))
562  y2a(i,k) = (ya(i,k+1)-ya(i,k-1))/(x2a(k+1)-x2a(k-1))
563  y12a(i,k) = (ya(i+1,k+1)-ya(i+1,k-1)-ya(i,k+1)+ya(i,k-1)) &
564  /((x1a(i+1)-x1a(i))*(x2a(k+1)-x2a(k-1)))
565  END DO
566 
567  k = 1
568  DO i = 2,i_max-1
569  y1a(i,k) = (ya(i+1,k)-ya(i-1,k))/(x1a(i+1)-x1a(i-1))
570  ! write (*,*) y1a(i,k),ya(i,k),(x1a(i+2)-x1a(i))
571  ! read (*,*)
572  y2a(i,k) = (-ya(i,k+2)+4.0*ya(i,k+1)-3.0*ya(i,k)) /(x2a(k+2)-x2a(k))
573  y12a(i,k) = (ya(i+1,k+1)-ya(i+1,k)-ya(i,k+1)+ya(i,k)) &
574  /((x1a(i+1)-x1a(i))*(x2a(k+1)-x2a(k)))
575  END DO
576 
577 
578  i = i_max
579  DO k = 2,k_max-1
580  y1a(i,k) = (ya(i,k)-ya(i-1,k))/(x1a(i)-x1a(i-1))
581  y2a(i,k) = (ya(i,k+1)-ya(i,k-1))/(x2a(k+1)-x2a(k-1))
582  y12a(i,k) = (ya(i,k+1)-ya(i,k-1)-ya(i-1,k+1)+ya(i-1,k-1)) &
583  /((x1a(i)-x1a(i-1))*(x2a(k+1)-x2a(k-1)))
584  END DO
585 
586 
587  k = k_max
588  DO i = 2,i_max-1
589  y1a(i,k) = (ya(i+1,k)-ya(i-1,k))/(x1a(i+1)-x1a(i-1))
590  y2a(i,k) = (ya(i,k)-ya(i,k-1))/(x2a(k)-x2a(k-1))
591  y12a(i,k) = (ya(i+1,k)-ya(i+1,k-1)-ya(i-1,k)+ya(i-1,k-1)) &
592  /((x1a(i+1)-x1a(i-1))*(x2a(k)-x2a(k-1)))
593  END DO
594 
595  k = k_max
596  i = i_max
597  y1a(i,k) = (ya(i,k)-ya(i-1,k))/(x1a(i)-x1a(i-1))
598  y2a(i,k) = (ya(i,k)-ya(i,k-1))/(x2a(k)-x2a(k-1))
599  y12a(i,k) = (ya(i,k)-ya(i,k-1)-ya(i-1,k)+ya(i-1,k-1)) &
600  /((x1a(i)-x1a(i-1))*(x2a(k)-x2a(k-1)))
601 
602  k = 1
603  i = 1
604  y1a(i,k) = (-ya(i+2,k)+4.0*ya(i+1,k)-3.0*ya(i,k)) /(x1a(i+2)-x1a(i))
605  y2a(i,k) = (-ya(i,k+2)+4.0*ya(i,k+1)-3.0*ya(i,k)) /(x2a(k+2)-x2a(k))
606  y12a(i,k) = (ya(i+1,k+1)-ya(i+1,k)-ya(i,k+1)+ya(i,k)) &
607  /((x1a(i+1)-x1a(i))*(x2a(k+1)-x2a(k)))
608 
609 
610  DO i = 1,i_max
611  DO k = 1,k_max
612  xspl(k) = x2a(k)
613  yspl(k) = ya(i,k)
614  END DO
615  yp1 = (-ya(i,1+2)+4.0*ya(i,1+1)-3.0*ya(i,1))/(x2a(1+2)-x2a(1))
616  ypn = (-ya(i,k_max-2)+4.0*ya(i,k_max-1)-3.0*ya(i,k_max)) &
617  /(x2a(k_max-2)-x2a(k_max))
618  CALL spline2(xspl,yspl,k_max,yp1,ypn,y1,y2)
619  DO k = 1,k_max
620  y2a(i,k) = y1(k)
621  d2f_dx2_2(i,k) = y2(k)
622  END DO
623  END DO
624 
625  DO k = 1,k_max
626  DO i = 1,i_max
627  xspl(i) = x1a(i)
628  yspl(i) = ya(i,k)
629  END DO
630  yp1 = (-ya(1+2,k)+4.0*ya(1+1,k)-3.0*ya(1,k))/(x1a(1+2)-x1a(1))
631  ypn = (-ya(i_max-2,k)+4.0*ya(i_max-1,k)-3.0*ya(i_max,k)) &
632  /(x1a(i_max-2)-x1a(i_max))
633  CALL spline2(xspl,yspl,i_max,yp1,ypn,y1,y2)
634  DO i = 1,i_max
635  y1a(i,k) = y1(i)
636  d2f_dx1_2(i,k) = y2(i)
637  END DO
638  END DO
639 
640  DO lm = 1,k_max-2
641  DO i = lm,i_max
642  k = i - (lm-1)
643  ! steps = (i_max-1) - (lm-1)
644  ! dz = (x1a(2)-x1a(1))**2 + (x2a(2)-x2a(1))**2
645  ! dz = dz**0.5
646  ! fct(k-1) = ya(i,k)
647  ! write (*,*) i-1,k-1,steps
648  xspl(k) = ( x1a(k)**2 + x2a(k)**2 )**0.5
649  yspl(k) = ya(i,k)
650  END DO
651  ! read (*,*)
652  ! CALL SPLINE_PARA (dz,fct,utri,steps)
653  ! CALL SPLINE_COEFF(beta,gamma,delta,dz,fct,utri,steps)
654  steps = i_max - (lm-1)
655  yp1 = y1a(lm,1)*gam1 + y2a(lm,1)*gam2
656  ypn = y1a(i_max,steps)*gam1+y2a(i_max,steps)*gam2
657  CALL spline2(xspl,yspl,steps,yp1,ypn,y1,y2)
658  DO i = lm,i_max
659  k = i - (lm-1)
660  ! secderiv = 2.0*gamma(k-1)
661  ! c write (*,*) secderiv,y2(k)
662  ! secderiv = y2(k)
663  ! write (*,*) y12a(i,k),i,k
664  y12a(i,k) = 0.5*y2(k)/gam1/gam2 -0.5*(d2f_dx1_2(i,k)*gam1/gam2 &
665  +d2f_dx2_2(i,k)*gam2/gam1)
666  ! write (*,*) d2f_dx1_2(i,k),d2f_dx2_2(i,k)
667  ! write (*,*) y12a(i,k),0.5*secderiv, -0.5*(d2f_dx1_2(i,k)+d2f_dx2_2(i,k))
668  ! read (*,*)
669  END DO
670  END DO
671 
672  DO lm = 1,k_max-2
673  DO i = 1,i_max-(lm-1)
674  k = i + (lm-1)
675  ! steps = (i_max-1) - (lm-1)
676  ! dz = (x1a(2)-x1a(1))**2 + (x2a(2)-x2a(1))**2
677  ! dz = dz**0.5
678  ! fct(i-1) = ya(i,k)
679  ! c write (*,*) i,k,steps
680  xspl(i) = ( x1a(i)**2 + x2a(i)**2 )**0.5
681  yspl(i) = ya(i,k)
682  END DO
683  ! read (*,*)
684  ! CALL SPLINE_PARA (dz,fct,utri,steps)
685  ! CALL SPLINE_COEFF(beta,gamma,delta,dz,fct,utri,steps)
686  steps = i_max - (lm-1)
687  yp1 = y1a(1,lm)*gam1 +y2a(1,lm)*gam2
688  ypn = y1a(steps,k_max)*gam1+y2a(steps,k_max)*gam2
689  CALL spline2(xspl,yspl,steps,yp1,ypn,y1,y2)
690  DO i = 1,i_max-(lm-1)
691  k = i + (lm-1)
692  ! secderiv = 2.0*gamma(i-1)
693  ! write (*,*) secderiv,y2(i)
694  ! write (*,*) y12a(i,k),i,k
695  ! y12a(i,k) = -0.5*y2(i) +1.0*(d2f_dx1_2(i,k)+d2f_dx2_2(i,k))
696  y12a(i,k) = 0.5*y2(i)/gam1/gam2 - 0.5*(d2f_dx1_2(i,k)*gam1/gam2 + d2f_dx2_2(i,k)*gam2/gam1)
697  ! write (*,*) y12a(i,k),0.5*y2(i), -0.5*(d2f_dx1_2(i,k)+d2f_dx2_2(i,k))
698  ! read (*,*)
699  END DO
700  END DO
701 
702 END SUBROUTINE bicub_derivative2
703 
704 
705 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
706 ! SUBROUTINE spline2
707 !
708 ! Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., yi = f(xi), with
709 ! x1 < x2 < .. . < xN, and given values yp1 and ypn for the first derivative of the interpolating
710 ! function at points 1 and n, respectively, this routine returns an array y2(1:n) of
711 ! length n which contains the second derivatives of the interpolating function at the tabulated
712 ! points xi. If yp1 and/or ypn are equal to 1 � 1030 or larger, the routine is signaled to set
713 ! the corresponding boundary condition for a natural spline, with zero second derivative on
714 ! that boundary.
715 ! Parameter: NMAX is the largest anticipated value of n.
716 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
717 
718 SUBROUTINE spline2(x,y,n,yp1,ypn,y1,y2)
719 
720  IMPLICIT NONE
721 
722  !-----------------------------------------------------------------------------
723  REAL, INTENT(IN) :: x(n)
724  REAL, INTENT(IN) :: y(n)
725  INTEGER, INTENT(IN) :: n
726  REAL, INTENT(IN OUT) :: yp1
727  REAL, INTENT(IN OUT) :: ypn
728  REAL, INTENT(OUT) :: y1(n)
729  REAL, INTENT(OUT) :: y2(n)
730  !-----------------------------------------------------------------------------
731  REAL :: dx
732  INTEGER, PARAMETER :: nmax = 500
733  INTEGER :: i,k
734  REAL :: p,qn,sig,un,u(nmax)
735  !-----------------------------------------------------------------------------
736 
737 
738  IF (yp1 > .99e30) THEN
739  y2(1) = 0.0
740  u(1) = 0.0
741  ELSE
742  y2(1) = -0.5
743  u(1) = (3.0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
744  END IF
745  DO i = 2,n-1
746  sig = (x(i)-x(i-1))/(x(i+1)-x(i-1))
747  p = sig*y2(i-1)+2.0
748  y2(i) = (sig-1.)/p
749  u(i) = (6.0*((y(i+1)-y(i))/(x(i+1)-x(i)) &
750  -(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
751  END DO
752  IF (ypn > .99d30) THEN
753  qn = 0.0
754  un = 0.0
755  ELSE
756  qn = 0.5
757  un = (3.0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
758  END IF
759  y2(n) = (un-qn*u(n-1))/(qn*y2(n-1)+1.)
760  DO k = n-1,1,-1
761  y2(k) = y2(k)*y2(k+1)+u(k)
762  END DO
763  DO i = 1,n-1
764  dx = x(i+1)-x(i)
765  y1(i) = (y(i+1)-y(i))/dx - 1.0/3.0*dx*y2(i)-y2(i+1)/6.0*dx
766  END DO
767  i = n-1
768  dx = x(i+1)-x(i)
769  y1(n) = (y(i+1)-y(i))/dx + 1.0/6.0*dx*y2(i)+1.0/3.0*y2(i+1)*dx
770  ! y1(1) = yp1
771  ! y1(n) = ypn
772 
773 END SUBROUTINE spline2
774 
775 
776 
777 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
778 !
779 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
780 
781 SUBROUTINE bicub_c2 ( i_max, k_max )
782 
783  USE rgt_phase_space_cell
784  IMPLICIT NONE
785 
786  !-----------------------------------------------------------------------------
787  INTEGER, INTENT(IN) :: i_max
788  INTEGER, INTENT(IN) :: k_max
789  !-----------------------------------------------------------------------------
790  INTEGER :: i,k,m,n
791  REAL :: y(4),y1(4),y2(4),y12(4),x1l,x1u,x2l,x2u
792  REAL :: c(4,4)
793  !-----------------------------------------------------------------------------
794 
795  DO i = 1,i_max-1
796  DO k = 1,k_max-1
797  y(1) = ya(i,k)
798  y(2) = ya(i+1,k)
799  y(3) = ya(i+1,k+1)
800  y(4) = ya(i,k+1)
801 
802  y1(1) = y1a(i,k)
803  y1(2) = y1a(i+1,k)
804  y1(3) = y1a(i+1,k+1)
805  y1(4) = y1a(i,k+1)
806 
807  y2(1) = y2a(i,k)
808  y2(2) = y2a(i+1,k)
809  y2(3) = y2a(i+1,k+1)
810  y2(4) = y2a(i,k+1)
811 
812  y12(1) = y12a(i,k)
813  y12(2) = y12a(i+1,k)
814  y12(3) = y12a(i+1,k+1)
815  y12(4) = y12a(i,k+1)
816 
817  x1l = x1a(i)
818  x1u = x1a(i+1)
819  x2l = x2a(k)
820  x2u = x2a(k+1)
821 
822  CALL bcucof(y,y1,y2,y12,x1u-x1l,x2u-x2l,c)
823  DO m = 1,4
824  DO n = 1,4
825  c_bicub(i,k,m,n) = c(m,n)
826  END DO
827  END DO
828 
829  END DO
830  END DO
831 
832 
833 END SUBROUTINE bicub_c2
834 
835 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
836 ! SUBROUTINE bcuint2
837 !
838 ! this is a variant of bcuint, where higher derivatives are calculated.
839 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
840 
841 SUBROUTINE bcuint2 (y,y1,y2,y12,x1l,x1u,x2l,x2u,x1,x2,c,ansy,ansy1, &
842  ansy2,ansy11,ansy12,ansy22,ansy111,ansy112,ansy122,ansy222)
843 
844  use utilities
845  IMPLICIT NONE
846 
847  !-----------------------------------------------------------------------------
848  REAL, INTENT(IN OUT) :: y(4)
849  REAL, INTENT(IN OUT) :: y1(4)
850  REAL, INTENT(IN OUT) :: y2(4)
851  REAL, INTENT(IN OUT) :: y12(4)
852  REAL, INTENT(IN OUT) :: x1l
853  REAL, INTENT(IN OUT) :: x1u
854  REAL, INTENT(IN OUT) :: x2l
855  REAL, INTENT(IN OUT) :: x2u
856  REAL, INTENT(IN OUT) :: x1
857  REAL, INTENT(IN OUT) :: x2
858  REAL, INTENT(IN) :: c(4,4)
859  REAL, INTENT(OUT) :: ansy
860  REAL, INTENT(OUT) :: ansy1
861  REAL, INTENT(OUT) :: ansy2
862  REAL, INTENT(OUT) :: ansy11
863  REAL, INTENT(OUT) :: ansy12
864  REAL, INTENT(OUT) :: ansy22
865  REAL, INTENT(OUT) :: ansy111
866  REAL, INTENT(OUT) :: ansy112
867  REAL, INTENT(OUT) :: ansy122
868  REAL, INTENT(OUT) :: ansy222
869  !-----------------------------------------------------------------------------
870  REAL :: d1,d2
871  !U USES bcucof
872  INTEGER :: i,j
873  REAL :: t,u
874  !-----------------------------------------------------------------------------
875 
876  ! call bcucof(y,y1,y2,y12,x1u-x1l,x2u-x2l,c)
877  IF( x1u == x1l .OR. x2u == x2l ) call paus ('bad input in bcuint')
878  t = (x1-x1l) / (x1u-x1l)
879  u = (x2-x2l) / (x2u-x2l)
880  IF (u == 0.0) u = 1.e-15
881  IF (t == 0.0) t = 1.e-15
882  ansy = 0.0
883  ansy1 = 0.0
884  ansy2 = 0.0
885  ansy11 = 0.0
886  ansy12 = 0.0
887  ansy22 = 0.0
888  ansy111 = 0.0
889  ansy112 = 0.0
890  ansy122 = 0.0
891  ansy222 = 0.0
892  DO i = 1,4
893  DO j = 1,4
894  ansy = ansy + c(i,j)*u**REAL(j-1) *t**REAL(i-1)
895  ! write (*,*) ansy,c(i,j),u**REAL(j-1),t**REAL(i-1)
896  ansy1 = ansy1+ c(i,j)*REAL(i-1)*u**REAL(j-1) *t**REAL(i-2)
897  ansy2 = ansy2+ c(i,j)*REAL(j-1)*u**REAL(j-2) *t**REAL(i-1)
898 
899  ansy11 = ansy11+c(i,j)*REAL(i-1)*REAL(i-2)*u**REAL(j-1) *t**REAL(i-3)
900  ansy22 = ansy22+c(i,j)*REAL(j-1)*REAL(j-2)*u**REAL(j-3) *t**REAL(i-1)
901  ansy12 = ansy12+c(i,j)*REAL(i-1)*REAL(j-1)*u**REAL(j-2) *t**REAL(i-2)
902 
903  ansy111 = ansy111+c(i,j)*REAL(i-1)*REAL(i-2)*REAL(i-3) *u**REAL(j-1) *t**REAL(i-4)
904  ansy222 = ansy222+c(i,j)*REAL(j-1)*REAL(j-2)*REAL(j-3) *u**REAL(j-4) *t**REAL(i-1)
905  ansy112 = ansy112+c(i,j)*REAL(i-1)*REAL(i-2)*REAL(j-1) *u**REAL(j-2) *t**REAL(i-3)
906  ansy122 = ansy122+c(i,j)*REAL(j-1)*REAL(j-2)*REAL(i-1) *u**REAL(j-3) *t**REAL(i-2)
907  END DO
908  END DO
909  ! stop
910  d1 = (x1u-x1l)
911  d2 = (x2u-x2l)
912  ansy1 = ansy1/d1
913  ansy2 = ansy2/d2
914 
915  ansy11 = ansy11/d1/d1
916  ansy22 = ansy22/d2/d2
917  ansy12 = ansy12/d1/d2
918 
919  ansy111 = ansy111/d1**3
920  ansy222 = ansy222/d2**3
921  ansy112 = ansy112/d1/d1/d2
922  ansy122 = ansy122/d1/d2/d2
923 
924 END SUBROUTINE bcuint2
925 
926 
927 
928 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
929 ! SUBROUTINE bi_cub_spline_crt
930 !
931 ! ............
932 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
933 
934 SUBROUTINE bi_cub_spline_crt (rho1,rho2, fr,fdr1,fdr2,fdr11,fdr12,fdr22, &
935  fdr111,fdr112,fdr122,fdr222,i_max,k_max,ih,k)
936 
937  USE rgt_phase_space_cell
938  IMPLICIT NONE
939 
940  !-----------------------------------------------------------------------------
941  REAL, INTENT(IN OUT) :: rho1
942  REAL, INTENT(IN OUT) :: rho2
943  REAL, INTENT(OUT) :: fr
944  REAL, INTENT(IN OUT) :: fdr1
945  REAL, INTENT(IN OUT) :: fdr2
946  REAL, INTENT(IN OUT) :: fdr11
947  REAL, INTENT(IN OUT) :: fdr12
948  REAL, INTENT(IN OUT) :: fdr22
949  REAL, INTENT(IN OUT) :: fdr111
950  REAL, INTENT(IN OUT) :: fdr112
951  REAL, INTENT(IN OUT) :: fdr122
952  REAL, INTENT(IN OUT) :: fdr222
953  INTEGER, INTENT(IN) :: i_max
954  INTEGER, INTENT(IN) :: k_max
955  INTEGER, INTENT(OUT) :: ih
956  INTEGER, INTENT(IN) :: k
957  !-----------------------------------------------------------------------------
958  INTEGER :: m,n
959  REAL :: y(4),y1(4),y2(4),y12(4),x1l,x1u,x2l,x2u
960  REAL :: c(4,4)
961  !-----------------------------------------------------------------------------
962 
963  IF (rho1 < x1a(1) .AND. rho2 < x2a(1)) THEN
964  fr = 1.0
965  WRITE (*,*) 'error? BI_CUB_SPLINE_CRT'
966  stop
967  RETURN
968  END IF
969  ! write (*,*) ih,k
970  IF (x1a(ih) <= rho1.AND.rho1 < x1a(ih+1).AND. &
971  x2a(k) <= rho2.AND.rho2 < x2a(k+1)) THEN
972  GO TO 10
973  ELSE
974  WRITE (*,*) 'error in BI_CUB_SPLINE_CRT',ih,k
975  WRITE (*,*) rho1,x1a(ih),x1a(ih+1)
976  WRITE (*,*) rho2,x2a(k),x2a(k+1)
977  stop
978  END IF
979  IF (ih > 2) THEN
980  IF (x1a(ih-1) <= rho1.AND.rho1 < x1a(ih)) THEN
981  ih = ih-1
982  GO TO 10
983  END IF
984  END IF
985  WRITE (*,*) 'error in BI_CUB_SPLINE_CRT'
986  CALL hunt(x1a,i_max,rho1,ih)
987 10 CONTINUE
988 
989  y(1) = ya(ih,k)
990  y(2) = ya(ih+1,k)
991  y(3) = ya(ih+1,k+1)
992  y(4) = ya(ih,k+1)
993 
994  y1(1) = y1a(ih,k)
995  y1(2) = y1a(ih+1,k)
996  y1(3) = y1a(ih+1,k+1)
997  y1(4) = y1a(ih,k+1)
998 
999  y2(1) = y2a(ih,k)
1000  y2(2) = y2a(ih+1,k)
1001  y2(3) = y2a(ih+1,k+1)
1002  y2(4) = y2a(ih,k+1)
1003 
1004  y12(1) = y12a(ih,k)
1005  y12(2) = y12a(ih+1,k)
1006  y12(3) = y12a(ih+1,k+1)
1007  y12(4) = y12a(ih,k+1)
1008 
1009  x1l = x1a(ih)
1010  x1u = x1a(ih+1)
1011  x2l = x2a(k)
1012  x2u = x2a(k+1)
1013 
1014  DO m = 1,4
1015  DO n = 1,4
1016  c(m,n) = c_bicub(ih,k,m,n)
1017  END DO
1018  END DO
1019  CALL bcuint2(y,y1,y2,y12,x1l,x1u,x2l,x2u,rho1,rho2,c, &
1020  fr,fdr1,fdr2,fdr11,fdr12,fdr22,fdr111,fdr112,fdr122,fdr222)
1021 
1022 END SUBROUTINE bi_cub_spline_crt
1023 
1024 
1025 
1026 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1027 ! SUBROUTINE PHI_CRITICAL_RENORM_PHASE_SPACE_CELL
1028 !
1029 ! ............
1030 !WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW
1031 
1032 SUBROUTINE phi_critical_renorm_phase_space_cell
1033 
1034  USE parameters, ONLY: kbol
1035  USE eos_constants
1036  USE eos_variables
1037  USE rgt_phase_space_cell
1038  IMPLICIT NONE
1039 
1040  !-----------------------------------------------------------------------------
1041  INTEGER :: k
1042  REAL :: zres, zges
1043  !-----------------------------------------------------------------------------
1044 
1045  CALL density_iteration
1046 
1047  !-----------------------------------------------------------------------------
1048  zges = (p * 1.e-30) / (kbol*t*eta/z3t)
1049  IF ( ensemble_flag == 'tv' ) zges = (pges * 1.e-30) / (kbol*t*eta/z3t)
1050  zres = zges - 1.0
1051 
1052  DO k = 1, ncomp
1053  IF (ensemble_flag == 'tp') lnphi(k) = dfdx(k) - log(zges)
1054  IF (ensemble_flag == 'tv' .AND. eta >= 0.0) lnphi(k) = dfdx(k) ! +LOG(rho)
1055  END DO
1056 
1057 END SUBROUTINE phi_critical_renorm_phase_space_cell
parameters
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains constant...
Definition: modules.f90:6
sig
double sig
correlated to variance of initial distribution
Definition: experimentalInputs.h:138
eos_constants
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:200
eos_variables
WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWW This module contains paramete...
Definition: modules.f90:120