ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
promat_e_I_f0.f
Go to the documentation of this file.
1  subroutine promat (n,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psim,iter,
2  * kstep,fvac,psi_eav,pspl_av,psibon0,cur_mu)
3 ! diffusion & equilibrium both equations solvinq
4 ! Newton metod with linearisation
5 ! error-correcting
6  implicit real*8 (a-h,o-z),integer*4 (i-n)
7  include 'dimpl1.inc'
8 
9  common /com_fixbon/ key_fixfree
10  common/com_flag/kastr
11 
12  common/selcon/ psi_d(nrp),fi_d(nrp),f_d(nrp),ri_d(nrp),
13  * ps_pnt(nrp),del_psb,psi_bn1
14  common/savt0/ psi0(nrp),fi0(nrp),f0(nrp),ri0(nrp),q0(nrp),
15  * dpsidt(nrp),dfidt(nrp),rm0,ac0n,skcen0
16 ! * ,kstep_prev
17  common/com_but/ sigma(nrp),cbut_b(nrp)
18  dimension dfdpsi(n),dpdpsi(n),f(n),flx_fi(n),psia(n),q(n)
19  dimension wrk(512)
20 ! save kstep_prev
21  integer kstep_prev
22 
23 
24  real*8 y(nrp,2)
25  real*8 a(nrp,4)
26  real*8 b(nrp,4)
27  real*8 c(nrp,4)
28  real*8 h(nrp,2)
29  real*8 zn(nrp,2)
30 
31  real*8 psi(nrp)
32  real*8 fi (nrp)
33  real*8 ri (nrp)
34  real*8 alfa22(nrp)
35  real*8 alfa33(nrp)
36  real*8 alp33k(nrp)
37  real*8 ds (nrp)
38  real*8 dsk (nrp)
39  real*8 dv (nrp)
40  real*8 dvk (nrp)
41  real*8 dpsi(nrp)
42  real*8 dfi (nrp)
43  real*8 dri (nrp)
44  real*8 df (nrp)
45  real*8 fk (nrp)
46  real*8 rik (nrp)
47  real*8 qk (nrp)
48  real*8 delf(nrp)
49  real*8 dfdpsn(nrp)
50 
51 !---------------------------------------------------------------------
52  !return
53  key_boncon=1
54  if(key_boncon.eq.0) then
55  dpdpsi(n)=0.d0
56  endif
57  epsel = 1.0d-10
58  dt = 1.d0
59  cap_psi=0.0d0 !1.d0
60  cap_fi =0.d0
61  !if(iter.lt.10) then
62  !cap_psi=0.d0
63  !cap_fi =0.d0
64  !endif
65 
66  fbnd = fvac
67 
68  !call psib_pla(pspl_av)
69  call cof_bon(cps_bon,bps_bon,dps_bon)
70  pspl_tst=-cps_bon*(psia(n)-psia(n-1))*psim
71  * +bps_bon*dfdpsi(n)+dps_bon*dpdpsi(n)
72 
73  psi_ext=psi_eav
74 
75  psi_bnd=psi_ext+pspl_av
76 
77  !!!psi_bnd=0.d0 !!!*
78  !!!psi_ext=-pspl_av !!!*
79 
80 ! metric coefficients
81  call metcof(alfa22,alfa33,ds,dv,ac0,skcen,ri,
82  * alp33k,dsk,dvk)
83 
84  do i=1,n
85  dfdpsn(i)=dfdpsi(i)
86  psi(i)=psim*psia(i)+psi_bnd
87  fi(i)=flx_fi(i)
88  if(kstep.ne.kstep_prev) then
89  jumpstep=kstep - kstep_prev
90  del_psb=psi_bn1-psibon0
91  if((jumpstep.eq.1 .AnD. kstep.gt.1) .OR. key_fixfree.eq.0)then
92  psibon0=psi_bn1
93  endif
94  psi0(i)=psim*psia(i)+psibon0
95  fi0(i)=fi(i)
96  f0(i) =f(i)
97  ri0(i) =ri(i)
98  q0(i) =q(i)
99  !dpsidt(i)=(psi(i)-psi0(i))/dt
100  dpsidt(i)=0.d0
101  dfidt(i )=0.d0
102  endif
103  !xa=dfloat(i-1)/dfloat(n-1)
104  xa=dsqrt(flx_fi(i)/flx_fi(n))
105  if(kastr.eq.1) then
106  sigma(i) =fun_sig(xa,i)
107  !sigma(i) =fun_sig_r(xa)
108  cbut_b(i) =fun_jb(xa,i)
109  !cbut_b(i) =fun_jb_rus(xa,dvk(i),dsk(i))
110  else
111  sigma(i) =fun_sig_a(xa,i)
112  cbut_b(i) =fun_jb_a(xa,i)
113  endif
114  enddo
115 
116  if(key_boncon.eq.0) then
117  cbut_b(n)=0.d0
118  sigma(n)=-1.d-6
119  endif
120 
121  if(iter.eq.1) then
122  rm0=rm
123  skcen0=skcen
124  ac0n=ac0
125  endif
126 
127 !---------------------------------------------------------------------
128 
129  do i=1,n
130  wrk(i)=dsqrt(dabs(alfa22(i)))
131  !wrk(i)=(alfa22(i)**2-alfa22(i-1)**2)/ dsk(i)
132  enddo
133 
134 
135  !do i=1,n-1
136  !fi(i)=fi(n)*(i-1.d0)**2/(n-1.d0)**2
137  !enddo
138 
139  do i=1,n-1
140  f(i)=alfa33(i)*(fi(i+1)-fi(i))/ds(i)
141  ri(i)=alfa22(i)*(psi(i+1)-psi(i))/ds(i)
142  q(i)=-(fi(i+1)-fi(i))/(psi(i+1)-psi(i))
143  enddo
144 
145  f(n)=fbnd
146 
147  ffprim = rm*( ac0*(psi(1)-psi(2))/skcen - rm*dpdpsi(1) )
148  delf(1) = ffprim - dfdpsi(1)
149  dfdpsi(1)=ffprim
150  do i=2,n-1
151  ffprim = (f(i)**2-f(i-1)**2)/(psi(i+1)-psi(i-1))
152  delf(i) = ffprim - dfdpsi(i)
153  dfdpsi(i) = ffprim
154  !dfdpsn(i) =alp33k(i)*(ri(i)-ri(i-1))/dsk(i)
155  enddo
156 
157  alfa22n=alfa22(n)
158  !rI(n)=dsqrt(ri(n-1)**2+alfa22n/alp33k(n)*(F(N)**2-f(n-1)**2))
159  dfdpsi(n) = (f(n)**2-f(n-1)**2)/(psi(n)-psi(n-1))
160 
161 
162 !---------------------------------------------------------------------
163  do ks = 1,100 ! iteration cycle begins
164 
165  do i=2,n-1
166  dpsi(i)=0.5d0*(psi(i+1)-psi(i-1))
167  dfi(i) =0.5d0*(fi(i+1)-fi(i-1))
168  dri(i) =ri(i)-ri(i-1)
169  df(i) =f(i)-f(i-1)
170  fk(i) =0.5d0*(f(i)+f(i-1))
171  rik(i) =0.5d0*(ri(i)+ri(i-1))
172  !Qk(i) =0.5d0*(q(i)+q(i-1))
173  qk(i) =-(dsk(i)/alp33k(i))*(f(i)+f(i-1))/(psi(i+1)-psi(i-1))
174  enddo
175 
176 ! matrix completing: a(i)*y(i-1)-c(i)*y(i)+b(i)*y(i+1)=-h(i)
177  do i=2,n-1
178 
179  dpsdt=dpsidt(i)
180  if(dpsdt.GE.0.d0) then
181  capfi=-cap_fi
182  else
183  capfi= cap_fi
184  endif
185 
186  dfldt=dfidt(i)
187  if(dfldt.LE.0.d0) then
188  capsi=-cap_psi
189  else
190  capsi= cap_psi
191  endif
192 
193  d2psi=0.5d0*(psi(i+1)-2.d0*psi(i)+psi(i-1))
194  d2fi=0.5d0*(fi(i+1)-2.d0*fi(i)+fi(i-1))
195 
196  a(i,1)= 0.5d0*dfidt(i)-( fk(i)+0.5d0*df(i) )*
197  + ( alfa22(i-1)/ds(i-1) )/sigma(i) !a11
198  + -0.5d0*dfidt(i)*capsi
199 
200  a(i,2)=-0.5d0*dpsidt(i)+( rik(i)+0.5d0*dri(i) )*
201  + ( alfa33(i-1)/ds(i-1) )/sigma(i) !a12
202  + +0.5d0*dpsidt(i)*capfi
203 
204  a(i,3)=-0.5d0*dri(i)/dpsi(i) + alfa22(i-1)/ds(i-1) !a21
205 
206  a(i,4)= qk(i)*alfa33(i-1)/ds(i-1)
207  + +0.5d0*(df(i)/dpsi(i))*(dsk(i)/alp33k(i))*(alfa33(i-1)/ds(i-1)) !a22
208 
209 !---------------------------------------------------------------------
210  b(i,1)=-0.5d0*dfidt(i)+( 0.5d0*df(i)-fk(i) )*
211  + ( alfa22(i)/ds(i) )/sigma(i) !b11
212  + -0.5d0*dfidt(i)*capsi
213 
214  b(i,2)= 0.5d0*dpsidt(i)+( rik(i)-0.5d0*dri(i) )*
215  + ( alfa33(i)/ds(i) )/sigma(i) !b12
216  + +0.5d0*dpsidt(i)*capfi
217 
218  b(i,3)= 0.5d0*dri(i)/dpsi(i)+alfa22(i)/ds(i) !b21
219 
220  b(i,4)= qk(i)*alfa33(i)/ds(i)
221  + -0.5d0*(df(i)/dpsi(i))*(dsk(i)/alp33k(i))*(alfa33(i)/ds(i)) !b22
222 
223 !---------------------------------------------------------------------
224  c(i,1)=-dfi(i)/dt-( fk(i)*(alfa22(i)/ds(i)+alfa22(i-1)/ds(i-1)) !c11
225  + -0.5d0*df(i)*( alfa22(i)/ds(i)-alfa22(i-1)/ds(i-1) ) )/sigma(i) !(-0.5..)
226  + -dfidt(i)*capsi !-capfi*d2fi/dt
227 
228  c(i,2)= dpsi(i)/dt+
229  + ( 0.5d0*dri(i)*(alfa33(i-1)/ds(i-1)-alfa33(i)/ds(i)) + !c12
230  + rik(i)*(alfa33(i)/ds(i)+alfa33(i-1)/ds(i-1)) )/sigma(i)
231  + +dpsidt(i)*capfi !+capsi*d2psi/dt
232 
233  c(i,3)= alfa22(i)/ds(i)+alfa22(i-1)/ds(i-1) !c21
234 
235  c(i,4)= qk(i)*( alfa33(i)/ds(i)+alfa33(i-1)/ds(i-1) ) !c22
236  + +0.5d0*(df(i)/dpsi(i))*(dsk(i)/alp33k(i))
237  + *( alfa33(i-1)/ds(i-1)-alfa33(i)/ds(i) )
238 !---------------------------------------------------------------------
239  fk0 =0.5d0*(f0(i)+f0(i-1))
240  rik0 =0.5d0*(ri0(i)+ri0(i-1))
241  qk0 =0.5d0*(q0(i)+q0(i-1))
242  dri0 =ri0(i)-ri0(i-1)
243  df0 =f0(i)-f0(i-1)
244 
245  h(i,1)= dpsidt(i)*dfi(i)-dfidt(i)*dpsi(i)
246  + +dpsidt(i)*capfi*d2fi-dfidt(i)*capsi*d2psi
247  + +( cbut_b(i)*dvk(i) )/sigma(i)
248 ! + -Fk(i)**2*( ri(i)/f(i)-ri(i-1)/f(i-1) )/sigma(i)
249  + -(f(i)*f(i-1))*( ri(i)/f(i)-ri(i-1)/f(i-1) )/sigma(i)
250 ! + -( Fk(i)*dri(i)-rIk(i)*df(i) )/sigma(i)
251 ! + -( Fk0*dri0-rIk0*df0 )/sigma(i)
252 
253 ! h(i,2)= Qk(i)*df(i)+dri(i)-dpdpsi(i)*(dv(i)+dv(i-1))*0.5d0
254  h(i,2)= dri(i)-dfdpsi(i)*dsk(i)/alp33k(i)-dpdpsi(i)*dvk(i)
255 !---------------------------------------------------------------------
256  enddo
257 
258 ! boundary equations
259  b(1,1)= rm*ac0/(sigma(1)*skcen)
260  b(1,2)= 0.d0
261  b(1,3)= 0.d0
262  b(1,4)= 0.d0
263 
264  c(1,1)= rm*ac0/(sigma(1)*skcen)-1.d0/dt !c11
265  c(1,2)= 0.d0 !c12
266  c(1,3)= 0.d0 !c21
267  c(1,4)= 1.d0 !c22
268 
269  h(1,1)= dpsidt(1)-rm*ac0*(psi(1)-psi(2))/(sigma(1)*skcen)
270  * +cbut_b(1)*rm**2/f(1)/sigma(1)
271 ! * -rm0*ac0n*(psi0(1)-psi0(2))/(sigma(1)*skcen0)
272  h(1,2)= 0.d0
273 
274 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275  if(key_boncon.eq.0) then
276  a(n,1)= cps_bon !a11
277  a(n,2)= 0.d0 !a12
278  a(n,3)= 0.d0 !a21
279  a(n,4)= 1.d0 !a22
280 
281  c(n,1)= cps_bon+1.d0 !c11
282  c(n,2)= 0.d0 !c12
283  c(n,3)= 0.d0 !c21
284  c(n,4)= 1.d0 !c22
285 
286  h(n,1)= psi_ext+cps_bon*psi(n-1)-(cps_bon+1.d0)*psi(n)
287  & +dps_bon*dpdpsi(n)
288  h(n,2)= fbnd*ds(n-1)/alfa33(n-1)-(fi(n)-fi(n-1))
289 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
290  else
291 ch4astra include 'pcof_e_in.f' !!*!
292  include 'pcof_e_in.inc' !!*!
293  endif
294 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
295 ! if(ks.eq.1) then
296 ! do i=1,n
297 ! !h(i,2)=0.d0
298 ! enddo
299 ! !h(1,1)=0.d0
300 ! endif
301 !
302 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
303 ! progonka
304  call mtrx_prog(n,y,a,b,c,h)
305 
306 
307  do i=1,n
308  psi(i)= psi(i) + y(i,1) ! y(i,1) = delta(psi(i))
309  fi(i) = fi(i) + y(i,2) ! y(i,2) = delta(fi (i))
310  dpsidt(i)=(psi(i)-psi0(i))/dt
311  dfidt(i )=(fi(i)-fi0(i))/dt
312  enddo
313 
314  do i=1,n-1
315  f(i)=alfa33(i)*(fi(i+1)-fi(i))/ds(i)
316  ri(i)=alfa22(i)*(psi(i+1)-psi(i))/ds(i)
317  q(i)=-(fi(i+1)-fi(i))/(psi(i+1)-psi(i))
318  enddo
319 
320  ffprim=rm*(ac0*(psi(1)-psi(2))/skcen-rm*dpdpsi(1))
321  delf(1) = ffprim - dfdpsi(1)
322  dfdpsi(1)=ffprim
323  !dfdpsi(1)=0.5d0*(dfdps+dfdpsi(1))
324  do i=2,n-1
325  ffprim = (f(i)**2-f(i-1)**2)/(psi(i+1)-psi(i-1))
326  delf(i) = ffprim - dfdpsi(i)
327  dfdpsi(i) = ffprim
328  !dfdpsn(i) =alp33k(i)*(ri(i)-ri(i-1))/dsk(i)
329  enddo
330  !dfdpsi(n) = dfdpsi(n-1) !<<<<<<<<<<<<<
331  alfa22n=alfa22(n)
332  f(n)=fbnd
333  ri(n)=dsqrt( ri(n-1)**2+alfa22n/alp33k(n)*(f(n)**2-f(n-1)**2)
334  * +alfa22n*dpdpsi(n)*(psi(n)-psi(n-1))*dvk(n)/dsk(n) )
335  dfdpsi(n) = (f(n)**2-f(n-1)**2)/(psi(n)-psi(n-1))
336 
337  deltamax = 0.d0
338  do k=1,2
339  do i=1,n
340  aby = dabs(y(i,k))
341  if (aby.ge.deltamax) deltamax=aby
342  enddo
343  enddo
344 
345  do i=2,n-1
346 
347  dpsdt=dpsidt(i)
348  if(dpsdt.GE.0.d0) then
349  capfi=-cap_fi
350  else
351  capfi= cap_fi
352  endif
353 
354  dfldt=dfidt(i)
355  if(dfldt.LE.0.d0) then
356  capsi=-cap_psi
357  else
358  capsi= cap_psi
359  endif
360 
361  d2psi= (psi(i+1)-2.d0*psi(i)+psi(i-1))
362  d2fi= (fi(i+1)-2.d0*fi(i)+fi(i-1))
363 
364  znv=(dpsidt(i)*(fi(i+1)-fi(i-1))-dfidt(i)*(psi(i+1)-psi(i-1)))
365  *-(f(i)*f(i-1))*( ri(i)/f(i)-ri(i-1)/f(i-1) )*2.0d0
366 ! *-((f(i)+f(i-1))*(ri(i)-ri(i-1))-(f(i)-f(i-1))*(ri(i)+ri(i-1)))
367  * /sigma(i)
368  + +dpsidt(i)*capfi*d2fi-dfidt(i)*capsi*d2psi
369 ! *-( (f0(i)+f0(i-1))*(ri0(i)-ri0(i-1))
370 ! * -(f0(i)-f0(i-1))*(ri0(i)+ri0(i-1)) )
371 ! * /sigma
372  zn(i,1)=0.5d0*znv
373 
374 ! zn(i,2)= (ri(i)-ri(i-1)) + 0.5d0*(q(i)+q(i-1))*(f(i)-f(i-1))
375 ! * - dpdpsi(i)*(dv(i)+dv(i-1))*0.5d0
376  zn(i,2)= (ri(i)-ri(i-1)) - dfdpsi(i)*dsk(i)/alp33k(i)
377  * - dpdpsi(i)*dvk(i)
378  enddo
379 
380  ddivmax = 0.d0
381  do k=1,2
382  do i=1,n
383  abd = dabs(zn(i,k))
384  if (abd.ge.ddivmax) ddivmax=abd
385  enddo
386  enddo
387 
388  cur_mu=ri(n)
389 
390  if (deltamax.le.epsel) exit
391  !if (ddivmax.le.epsel) exit
392 
393  if(ks.eq.199) then
394  write(*,*) 'no newton iterations convergence in promat'
395  write(*,*) 'execution was terminated'
396  write(*,*) 'ks=',ks
397  stop
398  endif
399 
400  enddo ! iteration cycle ends
401 !---------------------------------------------------------------------
402 
403  kstep_prev=kstep
404  psi_bn1=psi(n)
405 
406  10 continue
407  f_wght=0.75d0
408  do i=1,n
409  dfdpsi(i)=f_wght*dfdpsi(i)+(1.d0-f_wght)*dfdpsn(i)
410  enddo
411 
412 ! For selfconsistency checking
413 
414  do i=1,n
415  psi_d(i)=psi(i)-psi(n)
416  fi_d(i)=fi(i)
417  f_d(i)=f(i)
418  ri_d(i)=ri(i)
419  ps_pnt(i)=dpsidt(i)
420  enddo
421 
422 !---------------------------------------------------------------------
423  !write(*,*) 'promat'
424  return
425  end
426 
427 ! +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
428  subroutine promat_i (n,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psim,iter,
429  * kstep,fvac,psi_eav,pspl_av,psibon0,tok)
430 ! diffusion & equilibrium both equations solvinq
431 ! Newton metod with linearisation
432 ! error-correcting
433  implicit real*8 (a-h,o-z),integer*4 (i-n)
434  include 'dimpl1.inc'
435 
436  common/selcon/ psi_d(nrp),fi_d(nrp),f_d(nrp),ri_d(nrp),
437  * ps_pnt(nrp),del_psb,psi_bn1
438  common/savt0/ psi0(nrp),fi0(nrp),f0(nrp),ri0(nrp),q0(nrp),
439  * dpsidt(nrp),dfidt(nrp),rm0,ac0n,skcen0
440 ! * ,kstep_prev
441  common/com_but/ sigma(nrp),cbut_b(nrp)
442  common /com_fixbon/ key_fixfree
443  common/com_flag/kastr
444  dimension dfdpsi(n),dpdpsi(n),f(n),flx_fi(n),psia(n),q(n)
445  dimension wrk(512)
446 ! save kstep_prev
447  integer kstep_prev
448 
449  real*8 y(nrp,2)
450  real*8 a(nrp,4)
451  real*8 b(nrp,4)
452  real*8 c(nrp,4)
453  real*8 h(nrp,2)
454  real*8 zn(nrp,2)
455 
456  real*8 psi(nrp)
457  real*8 fi (nrp)
458  real*8 ri (nrp)
459  real*8 alfa22(nrp)
460  real*8 alfa33(nrp)
461  real*8 alp33k(nrp)
462  real*8 ds (nrp)
463  real*8 dsk (nrp)
464  real*8 dv (nrp)
465  real*8 dvk (nrp)
466  real*8 dpsi(nrp)
467  real*8 dfi (nrp)
468  real*8 dri (nrp)
469  real*8 df (nrp)
470  real*8 fk (nrp)
471  real*8 rik (nrp)
472  real*8 qk (nrp)
473  real*8 delf(nrp)
474  real*8 dfdpsn(nrp)
475 
476 !---------------------------------------------------------------------
477 
478  epsel = 1.0d-10
479  dt = 1.d0
480  cap_psi=0.d0 !1.d0
481  cap_fi =0.d0
482  !if(iter.lt.10) then
483  !cap_psi=0.d0
484  !cap_fi =0.d0
485  !endif
486 
487  fbnd = fvac
488 
489  ! call psib_pla(pspl_av)
490 ! call cof_bon(cps_bon,bps_bon,dps_bon)
491 ! pspl_tst=-cps_bon*(psia(n)-psia(n-1))*psim
492 ! * +bps_bon*dfdpsi(n)+dps_bon*dpdpsi(n)
493 !
494 ! psi_ext=psi_eav
495 !
496 ! psi_bnd=psi_ext+pspl_av
497 
498  if(kstep.eq.1)then
499  psi_bnd=psibon0
500  else
501  psi_bnd=psi_bn1
502  endif
503 
504  psi_eav=psi_bnd-pspl_av
505 
506 ! metric coefficients
507  call metcof(alfa22,alfa33,ds,dv,ac0,skcen,ri,
508  * alp33k,dsk,dvk)
509 
510  do i=1,n
511  dfdpsn(i)=dfdpsi(i)
512  psi(i)=psim*psia(i)+psi_bnd
513  fi(i)=flx_fi(i)
514 
515  if(kstep.ne.kstep_prev) then
516  jumpstep=kstep - kstep_prev
517  del_psb=psi_bn1-psibon0
518  if((jumpstep.eq.1 .AnD. kstep.gt.1) .OR. key_fixfree.eq.0)then
519  psibon0=psi_bn1
520  endif
521  psi0(i)=psim*psia(i)+psibon0
522  fi0(i)=fi(i)
523  f0(i) =f(i)
524  ri0(i) =ri(i)
525  q0(i) =q(i)
526  !dpsidt(i)=(psi(i)-psi0(i))/dt
527  dpsidt(i)=0.d0
528  dfidt(i )=0.d0
529  endif
530 
531 ! xa=dfloat(i-1)/dfloat(n-1)
532 ! sigma(i) =fun_sig(xa,i)
533 ! cbut_b(i) =fun_jb(xa,i)
534 
535  xa=dsqrt(flx_fi(i)/flx_fi(n))
536  if(kastr.eq.1) then
537  sigma(i) =fun_sig(xa,i)
538  !sigma(i) =fun_sig_r(xa)
539  cbut_b(i) =fun_jb(xa,i)
540  !cbut_b(i) =fun_jb_rus(xa,dvk(i),dsk(i))
541  else
542  sigma(i) =fun_sig_a(xa,i)
543  cbut_b(i) =fun_jb_a(xa,i)
544  endif
545 
546  enddo
547 
548 
549  if(iter.eq.1) then
550  rm0=rm
551  skcen0=skcen
552  ac0n=ac0
553  endif
554 
555 !---------------------------------------------------------------------
556 
557  do i=1,n
558  wrk(i)=dsqrt(dabs(alfa22(i)))
559  !wrk(i)=(alfa22(i)**2-alfa22(i-1)**2)/ dsk(i)
560  enddo
561 
562 
563  do i=1,n-1
564  fi(i)=fi(n)*(i-1.d0)**2/(n-1.d0)**2
565  enddo
566 
567  do i=1,n-1
568  f(i)=alfa33(i)*(fi(i+1)-fi(i))/ds(i)
569  ri(i)=alfa22(i)*(psi(i+1)-psi(i))/ds(i)
570  q(i)=-(fi(i+1)-fi(i))/(psi(i+1)-psi(i))
571  enddo
572 
573  f(n)=fbnd
574  ri(n)=tok
575 
576  ffprim=rm*(ac0*(psi(1)-psi(2))/skcen-rm*dpdpsi(1))
577  delf(1) = ffprim - dfdpsi(1)
578  dfdpsi(1)=ffprim
579  do i=2,n-1
580  ffprim = (f(i)**2-f(i-1)**2)/(psi(i+1)-psi(i-1))
581  delf(i) = ffprim - dfdpsi(i)
582  dfdpsi(i) = ffprim
583  dfdpsn(i) =alp33k(i)*(ri(i)-ri(i-1))/dsk(i)
584  enddo
585 
586  alfa22n=alfa22(n)
587  dfdpsi(n) = (f(n)**2-f(n-1)**2)/(psi(n)-psi(n-1))
588 
589 
590 !---------------------------------------------------------------------
591  do ks = 1,100 ! iteration cycle begins
592 
593  do i=2,n-1
594  dpsi(i)=0.5d0*(psi(i+1)-psi(i-1))
595  dfi(i) =0.5d0*(fi(i+1)-fi(i-1))
596  dri(i) =ri(i)-ri(i-1)
597  df(i) =f(i)-f(i-1)
598  fk(i) =0.5d0*(f(i)+f(i-1))
599  rik(i) =0.5d0*(ri(i)+ri(i-1))
600  !Qk(i) =0.5d0*(q(i)+q(i-1))
601  qk(i) =-(dsk(i)/alp33k(i))*(f(i)+f(i-1))/(psi(i+1)-psi(i-1))
602  enddo
603 
604 ! matrix completing: a(i)*y(i-1)-c(i)*y(i)+b(i)*y(i+1)=-h(i)
605  do i=2,n-1
606 
607  dpsdt=dpsidt(i)
608  if(dpsdt.GE.0.d0) then
609  capfi=-cap_fi
610  else
611  capfi= cap_fi
612  endif
613 
614  dfldt=dfidt(i)
615  if(dfldt.LE.0.d0) then
616  capsi=-cap_psi
617  else
618  capsi= cap_psi
619  endif
620 
621  d2psi=0.5d0*(psi(i+1)-2.d0*psi(i)+psi(i-1))
622  d2fi=0.5d0*(fi(i+1)-2.d0*fi(i)+fi(i-1))
623 
624  a(i,1)= 0.5d0*dfidt(i)-( fk(i)+0.5d0*df(i) )*
625  + ( alfa22(i-1)/ds(i-1) )/sigma(i) !a11
626  + -0.5d0*dfidt(i)*capsi
627 
628  a(i,2)=-0.5d0*dpsidt(i)+( rik(i)+0.5d0*dri(i) )*
629  + ( alfa33(i-1)/ds(i-1) )/sigma(i) !a12
630  + +0.5d0*dpsidt(i)*capfi
631 
632  a(i,3)=-0.5d0*dri(i)/dpsi(i) + alfa22(i-1)/ds(i-1) !a21
633 
634  a(i,4)= qk(i)*alfa33(i-1)/ds(i-1)
635  + +0.5d0*(df(i)/dpsi(i))*(dsk(i)/alp33k(i))*(alfa33(i-1)/ds(i-1)) !a22
636 
637 !---------------------------------------------------------------------
638  b(i,1)=-0.5d0*dfidt(i)+( 0.5d0*df(i)-fk(i) )*
639  + ( alfa22(i)/ds(i) )/sigma(i) !b11
640  + -0.5d0*dfidt(i)*capsi
641 
642  b(i,2)= 0.5d0*dpsidt(i)+( rik(i)-0.5d0*dri(i) )*
643  + ( alfa33(i)/ds(i) )/sigma(i) !b12
644  + +0.5d0*dpsidt(i)*capfi
645 
646  b(i,3)= 0.5d0*dri(i)/dpsi(i)+alfa22(i)/ds(i) !b21
647 
648  b(i,4)= qk(i)*alfa33(i)/ds(i)
649  + -0.5d0*(df(i)/dpsi(i))*(dsk(i)/alp33k(i))*(alfa33(i)/ds(i)) !b22
650 
651 !---------------------------------------------------------------------
652  c(i,1)=-dfi(i)/dt-( fk(i)*(alfa22(i)/ds(i)+alfa22(i-1)/ds(i-1)) !c11
653  + -0.5d0*df(i)*( alfa22(i)/ds(i)-alfa22(i-1)/ds(i-1) ) )/sigma(i) !(-0.5..)
654  + -dfidt(i)*capsi !-capfi*d2fi/dt
655 
656  c(i,2)= dpsi(i)/dt+
657  + ( 0.5d0*dri(i)*(alfa33(i-1)/ds(i-1)-alfa33(i)/ds(i)) + !c12
658  + rik(i)*(alfa33(i)/ds(i)+alfa33(i-1)/ds(i-1)) )/sigma(i)
659  + +dpsidt(i)*capfi !+capsi*d2psi/dt
660 
661  c(i,3)= alfa22(i)/ds(i)+alfa22(i-1)/ds(i-1) !c21
662 
663  c(i,4)= qk(i)*( alfa33(i)/ds(i)+alfa33(i-1)/ds(i-1) ) !c22
664  + +0.5d0*(df(i)/dpsi(i))*(dsk(i)/alp33k(i))
665  + *( alfa33(i-1)/ds(i-1)-alfa33(i)/ds(i) )
666 !---------------------------------------------------------------------
667  fk0 =0.5d0*(f0(i)+f0(i-1))
668  rik0 =0.5d0*(ri0(i)+ri0(i-1))
669  qk0 =0.5d0*(q0(i)+q0(i-1))
670  dri0 =ri0(i)-ri0(i-1)
671  df0 =f0(i)-f0(i-1)
672 
673  h(i,1)= dpsidt(i)*dfi(i)-dfidt(i)*dpsi(i)
674  + +dpsidt(i)*capfi*d2fi-dfidt(i)*capsi*d2psi
675  + +( cbut_b(i)*dvk(i) )/sigma(i)
676 ! + -Fk(i)**2*( ri(i)/f(i)-ri(i-1)/f(i-1) )/sigma(i)
677  + -(f(i)*f(i-1))*( ri(i)/f(i)-ri(i-1)/f(i-1) )/sigma(i)
678 ! + -( Fk(i)*dri(i)-rIk(i)*df(i) )/sigma(i)
679 ! + -( Fk0*dri0-rIk0*df0 )/sigma(i)
680 
681 ! h(i,2)= Qk(i)*df(i)+dri(i)-dpdpsi(i)*(dv(i)+dv(i-1))*0.5d0
682  h(i,2)= dri(i)-dfdpsi(i)*dsk(i)/alp33k(i)-dpdpsi(i)*dvk(i)
683 !---------------------------------------------------------------------
684  enddo
685 
686 ! boundary equations
687  b(1,1)= rm*ac0/(sigma(1)*skcen)
688  b(1,2)= 0.d0
689  b(1,3)= 0.d0
690  b(1,4)= 0.d0
691 
692  c(1,1)= rm*ac0/(sigma(1)*skcen)-1.d0/dt !c11
693  c(1,2)= 0.d0 !c12
694  c(1,3)= 0.d0 !c21
695  c(1,4)= 1.d0 !c22
696 
697  h(1,1)= dpsidt(1)-rm*ac0*(psi(1)-psi(2))/(sigma(1)*skcen)
698  * +cbut_b(1)*rm**2/f(1)/sigma(1)
699 ! * -rm0*ac0n*(psi0(1)-psi0(2))/(sigma(1)*skcen0)
700  h(1,2)= 0.d0
701 
702 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
703 ! a(n,1)= cps_bon ! 0.d0 !a11
704 ! a(n,2)= 0.d0 !a12
705 ! a(n,3)= 0.d0 !a21
706 ! a(n,4)=-alfa33(n-1)/ds(n-1) !a22
707 !
708 ! c(n,1)=cps_bon-1.d0 !c11
709 ! c(n,2)= 0.d0 !c12
710 ! c(n,3)= 0.d0 !c21
711 ! c(n,4)=-alfa33(n-1)/ds(n-1) !c22
712 !
713 ! h(n,1)=-psi_ext+cps_bon*psi(n-1)-(cps_bon-1.d0)*psi(n) !0.d0
714 ! h(n,2)=-fbnd+alfa33(n-1)*(fi(n)-fi(n-1))/ds(n-1) !dfi(n-1) ?
715 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
716 
717  !!!include 'cof_in.f'
718 ch4astra include 'pcof_e_I_in.f' !!!
719  include 'pcof_e_I_in.inc' !!!
720 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
721  if(ks.eq.1) then
722  do i=1,n
723  !h(i,2)=0.d0
724  enddo
725  !h(1,1)=0.d0
726  endif
727 
728 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
729 ! progonka
730  call mtrx_prog(n,y,a,b,c,h)
731 
732 
733  do i=1,n
734  psi(i)= psi(i) + y(i,1) ! y(i,1) = delta(psi(i))
735  fi(i) = fi(i) + y(i,2) ! y(i,2) = delta(fi (i))
736  dpsidt(i)=(psi(i)-psi0(i))/dt
737  dfidt(i )=(fi(i)-fi0(i))/dt
738  enddo
739 
740  do i=1,n-1
741  f(i)=alfa33(i)*(fi(i+1)-fi(i))/ds(i)
742  ri(i)=alfa22(i)*(psi(i+1)-psi(i))/ds(i)
743  q(i)=-(fi(i+1)-fi(i))/(psi(i+1)-psi(i))
744  enddo
745 
746  ffprim=rm*(ac0*(psi(1)-psi(2))/skcen-rm*dpdpsi(1))
747  delf(1) = ffprim - dfdpsi(1)
748  dfdpsi(1)=ffprim
749  !dfdpsi(1)=0.5d0*(dfdps+dfdpsi(1))
750  do i=2,n-1
751  ffprim = (f(i)**2-f(i-1)**2)/(psi(i+1)-psi(i-1))
752  delf(i) = ffprim - dfdpsi(i)
753  dfdpsi(i) = ffprim
754  dfdpsn(i) =alp33k(i)*(ri(i)-ri(i-1))/dsk(i)
755  enddo
756  !dfdpsi(n) = dfdpsi(n-1) !<<<<<<<<<<<<<
757  !F(N)=fbnd
758  !rI(n)=tok
759  alfa22n=alfa22(n)
760  dfdpsi(n) = (f(n)**2-f(n-1)**2)/(psi(n)-psi(n-1))
761 
762  deltamax = 0.d0
763  do k=1,2
764  do i=1,n
765  aby = dabs(y(i,k))
766  if (aby.ge.deltamax) deltamax=aby
767  enddo
768  enddo
769 
770  do i=2,n-1
771 
772  dpsdt=dpsidt(i)
773  if(dpsdt.GE.0.d0) then
774  capfi=-cap_fi
775  else
776  capfi= cap_fi
777  endif
778 
779  dfldt=dfidt(i)
780  if(dfldt.LE.0.d0) then
781  capsi=-cap_psi
782  else
783  capsi= cap_psi
784  endif
785 
786  d2psi= (psi(i+1)-2.d0*psi(i)+psi(i-1))
787  d2fi= (fi(i+1)-2.d0*fi(i)+fi(i-1))
788 
789  znv=(dpsidt(i)*(fi(i+1)-fi(i-1))-dfidt(i)*(psi(i+1)-psi(i-1)))
790  *-(f(i)*f(i-1))*( ri(i)/f(i)-ri(i-1)/f(i-1) )*2.0d0
791 ! *-((f(i)+f(i-1))*(ri(i)-ri(i-1))-(f(i)-f(i-1))*(ri(i)+ri(i-1)))
792  * /sigma(i)
793  + +dpsidt(i)*capfi*d2fi-dfidt(i)*capsi*d2psi
794 ! *-( (f0(i)+f0(i-1))*(ri0(i)-ri0(i-1))
795 ! * -(f0(i)-f0(i-1))*(ri0(i)+ri0(i-1)) )
796 ! * /sigma
797  zn(i,1)=0.5d0*znv
798 
799 ! zn(i,2)= (ri(i)-ri(i-1)) + 0.5d0*(q(i)+q(i-1))*(f(i)-f(i-1))
800 ! * - dpdpsi(i)*(dv(i)+dv(i-1))*0.5d0
801  zn(i,2)= (ri(i)-ri(i-1)) - dfdpsi(i)*dsk(i)/alp33k(i)
802  * - dpdpsi(i)*dvk(i)
803  enddo
804 
805  ddivmax = 0.d0
806  do k=1,2
807  do i=1,n
808  abd = dabs(zn(i,k))
809  if (abd.ge.ddivmax) ddivmax=abd
810  enddo
811  enddo
812 
813  if (deltamax.le.epsel) exit
814  !if (ddivmax.le.epsel) exit
815 
816  if(ks.eq.199) then
817  write(*,*) 'no newton iterations convergence in promat'
818  write(*,*) 'execution was terminated'
819  write(*,*) 'ks=',ks
820  stop
821  endif
822 
823  enddo ! iteration cycle ends
824 !---------------------------------------------------------------------
825 
826  kstep_prev=kstep
827  psi_bn1=psi(n)
828 
829  10 continue
830 
831 ! For selfconsistency checking
832 
833  do i=1,n
834  psi_d(i)=psi(i)-psi(n)
835  fi_d(i)=fi(i)
836  f_d(i)=f(i)
837  ri_d(i)=ri(i)
838  ps_pnt(i)=dpsidt(i)
839  enddo
840 
841 !---------------------------------------------------------------------
842  return
843  end
844 
845 ! +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
846 
847  subroutine mtrx_prog (n,y,a,b,c,f)
848 ! matrix progonka for two equations
849  implicit real*8 (a-h,o-z),integer*4 (i-n)
850  include 'dimpl1.inc'
851  dimension y(nrp,2),a(nrp,4),b(nrp,4),c(nrp,4),f(nrp,2)
852  real*8 alfa(nrp,4)
853  real*8 beta(nrp,2)
854 !---------------------------------------------------------------------
855 ! a(i)*y(i-1)-c(i)*y(i)+b(i)*y(i+1)=-f(i), i=2,n-1
856 ! -c(1)*y(1)+b(1)*y(2)=-f(1), a(n)*y(n-1)-c(n)*y(n)=-f(n)
857 !---------------------------------------------------------------------
858 ! a(i), b(i), c(i) -> matrix [2*2]; y(i), f(i) -> vector [2]
859 ! {a(i)11 a(i)12 a(i)21 a(i)22} -> {a(i)1 a(i)2 a(i)3 a(i)4}
860 !---------------------------------------------------------------------
861  det=c(1,1)*c(1,4)-c(1,2)*c(1,3) ! inverse c(1)
862  c1i1= c(1,4)/det
863  c1i2=-c(1,2)/det
864  c1i3=-c(1,3)/det
865  c1i4= c(1,1)/det
866 !---------------------------------------------------------------------
867  alfa(1,1) = c1i1*b(1,1)+c1i2*b(1,3)
868  alfa(1,2) = c1i1*b(1,2)+c1i2*b(1,4)
869  alfa(1,3) = c1i3*b(1,1)+c1i4*b(1,3)
870  alfa(1,4) = c1i3*b(1,2)+c1i4*b(1,4)
871 
872  beta(1,1) = c1i1*f(1,1)+c1i2*f(1,2)
873  beta(1,2) = c1i3*f(1,1)+c1i4*f(1,2)
874 !---------------------------------------------------------------------
875  do j=2,n
876  cmb1 = c(j,1) - a(j,1)*alfa(j-1,1)-a(j,2)*alfa(j-1,3)
877  cmb2 = c(j,2) - a(j,1)*alfa(j-1,2)-a(j,2)*alfa(j-1,4)
878  cmb3 = c(j,3) - a(j,3)*alfa(j-1,1)-a(j,4)*alfa(j-1,3)
879  cmb4 = c(j,4) - a(j,3)*alfa(j-1,2)-a(j,4)*alfa(j-1,4)
880 
881  det=cmb1*cmb4-cmb2*cmb3 ! inverse {c(i)-a(i)*alfa(i-1)}
882  cmbi1= cmb4/det
883  cmbi2=-cmb2/det
884  cmbi3=-cmb3/det
885  cmbi4= cmb1/det
886 
887  alfa(j,1) = cmbi1*b(j,1)+cmbi2*b(j,3)
888  alfa(j,2) = cmbi1*b(j,2)+cmbi2*b(j,4)
889  alfa(j,3) = cmbi3*b(j,1)+cmbi4*b(j,3)
890  alfa(j,4) = cmbi3*b(j,2)+cmbi4*b(j,4)
891 
892  fab1 = f(j,1) + a(j,1)*beta(j-1,1)+a(j,2)*beta(j-1,2)
893  fab2 = f(j,2) + a(j,3)*beta(j-1,1)+a(j,4)*beta(j-1,2)
894 
895  beta(j,1) = cmbi1*fab1+cmbi2*fab2
896  beta(j,2) = cmbi3*fab1+cmbi4*fab2
897  enddo
898 !---------------------------------------------------------------------
899  y(n,1) = beta(n,1)
900  y(n,2) = beta(n,2)
901 
902  do j=n-1,1,-1
903  y(j,1) = alfa(j,1)*y(j+1,1)+alfa(j,2)*y(j+1,2) + beta(j,1)
904  y(j,2) = alfa(j,3)*y(j+1,1)+alfa(j,4)*y(j+1,2) + beta(j,2)
905  enddo
906 !---------------------------------------------------------------------
907 !---------------------------------------------------------------------
908  return
909  end
910 
911 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
912 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
913  real*8 function fun_sig(a,i)
914  implicit real*8(a-h,o-z)
915  include 'dim.inc'
916  common /com_tim/ dtim,ctim
917  common /com_sigcd/ c_sig(nrp),t_el(nrp),c_bts(nrp),c_driv(nrp)
918 
919  fun_sig=-c_sig(i)*amu0/dtim
920  t= t_el(i)
921  fun_sig_tst=-(8.d0*pi/3.d0)*dabs(t)**1.5d0*(1.d-4/dtim)
922 
923  return
924  end
925 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
926 
927  real*8 function fun_sig_(a)
928  implicit real*8(a-h,o-z)
929  pi=3.14159265359d0
930  a0=0.8d0
931  a1=0.97d0
932  sigm0 =-5.5d6 !-1.d3 -1.d3 -5.d5
933  sigm1 =-1.5d3 !-1.d3 -1.d3 -1.d4
934  if(a.le.a0)then
935  sigma= sigm0
936  elseif(a.gt.a0 .AND. a.le.a1)then
937  sigma=sigm1+(sigm0-sigm1)*
938 ! * ( 0.5d0*(1.d0+dcos(pi*(a-a0)/(a1-a0))) ) !**0.5d0
939  * dabs( dcos(0.5d0*pi*(a-a0)/(a1-a0)) )**1.5
940  elseif(a.gt.a1)then
941  sigma= sigm1
942  endif
943  !sigma=sigm1 + sigm0*(1.d0-a**3)
944  fun_sig=sigma
945  return
946  end
947 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
948  real*8 function fun_sig__(a)
949  implicit real*8(a-h,o-z)
950  pi=3.14159265359d0
951  a0=0.95d0
952  a1=0.97d0
953  sigm0 =-5.d6 !-1.d3 -1.d3 -5.d5
954  sigm1 =-1.d2 !-1.d3 -1.d3 -1.d4
955  if(a.le.a0)then
956  sigma= sigm0
957  elseif(a.gt.a0 .AND. a.le.a1)then
958  sigma=sigm1+(sigm0-sigm1)*
959 ! * ( 0.5d0*(1.d0+dcos(pi*(a-a0)/(a1-a0))) ) !**0.5d0
960  * ( dcos(0.5d0*pi*(a-a0)/(a1-a0)) )**3
961  elseif(a.gt.a1)then
962  sigma= sigm1
963  endif
964  !sigma=sigm1 + sigm0*(1.d0-a**3)
965  fun_sig=sigma
966  return
967  end
968 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
969  real*8 function fun___sig(a)
970  implicit real*8(a-h,o-z)
971  pi=3.14159265359d0
972  a0=0.0d0
973  a1=0.5d0 !1.0d0!
974  a2=0.5075d0 !1.0d0!
975  t0 =5.0d3
976  t1 =3.0d3
977  t2 =1.0d3
978  tb =1.d2
979  if(a.ge.a2)then
980  t= tb + (t2-tb)*(1.d0-((a-a2)/(1.d0-a2))**2)
981  elseif(a.le.a0)then
982  t= t0
983  elseif(a.gt.a0 .AND. a.lt.a1)then
984  t=t1 + (t0-t1)*(1.d0-((a-a0)/(a1-a0))**2)
985  elseif(a.ge.a1 .AND. a.lt.a2)then
986  t=t2+(a-a2)*(t1-t2)/(a1-a2)
987  endif
988  fun_sig=-(8.d0*pi/3.d0)*dabs(t)**1.5d0
989  fun_sig=fun_sig/10.d0
990  return
991  end
992 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
993  !real*8 function fun_sig(a,i)
994  real*8 function fun_sig_a(a,i)
995  implicit real*8(a-h,o-z)
996  common /com_tim/ dtim,ctim
997  pi=3.14159265359d0
998  aeps=1.d-8
999  a0=0.0d0
1000  a1=1.0d0 !1.0d0!
1001  a2=1.0d0 !1.0d0!
1002  t0 =3.0d3
1003  t1 =1.0d2
1004  t2 =1.0d1
1005  tb =1.d1
1006 
1007  !if(a.ge.a2)then
1008  ! t= tb + (t2-tb)*(1.d0-((a-a2)/(1.d0-a2+aeps))**2)
1009  !elseif(a.le.a0)then
1010  ! t= t0
1011  !elseif(a.gt.a0 .AND. a.lt.a1)then
1012  ! t=t1 + (t0-t1)*(1.d0-((a-a0)/(a1-a0))**2)
1013  !elseif(a.ge.a1 .AND. a.lt.a2)then
1014  ! t=t2+(a-a2)*(t1-t2)/(a1-a2)
1015  !endif
1016 
1017  if(a.le.a0)then
1018  t= t0
1019  elseif(a.gt.a0 .AND. a.lt.a1)then
1020  t=t1 + (t0-t1)*(1.d0-((a-a0)/(a1-a0))**2)
1021  elseif(a.ge.a1)then
1022  t=t1
1023  endif
1024 
1025  !!t=tb
1026 
1027  fun_sig=-(8.d0*pi/3.d0)*dabs(t)**1.5d0
1028  fun_sig_a=fun_sig*(1.d-4/dtim)
1029 
1030  return
1031  end
1032 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1033  real*8 function fun_sigm(a)
1034  implicit real*8(a-h,o-z)
1035  pi=3.14159265359d0
1036  t0=3500.d0
1037  tp=0.2d0*t0
1038  ts=0.12d0*t0
1039  c=0.02d0
1040  d=0.98d0
1041  te=0.5d0*tp*( dtanh( (d-1.d0)/c ) +1.d0 )
1042  pa=2.d0
1043  pb=2.d0
1044  t=(t0-tp-ts+te)*(1-a**pa)**pb +.5*tp*(dtanh((d-a)/c) + 1)+ts-te
1045 
1046 
1047  fun_sig=-(8.d0*pi/3.d0)*dabs(t)**1.5d0
1048  fun_sig=fun_sig/2.5d0 !
1049  !!fun_sig=-8.d0
1050  return
1051  end
1052 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1053 
1054  real*8 function fun_jb(a,i)
1055  implicit real*8(a-h,o-z)
1056  include 'dim.inc'
1057  common /com_sigcd/ c_sig(nrp),t_el(nrp),c_bts(nrp),c_driv(nrp)
1058 
1059  fun_jb=(c_bts(i)+c_driv(i))*amu0
1060 
1061  return
1062  end
1063 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1064 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1065 
1066  real*8 function fun_jb_a(a,i)
1067  implicit real*8(a-h,o-z)
1068  include 'dimpl1.inc'
1069  common/com_jb/ bj_av(nrp),curfi_av(nrp)
1070 
1071  pi=3.14159265359d0
1072  ampl=0.0d0 !0.d0 !1.5d0 0.5d0
1073  a0=0.50d0
1074  w=0.05d0
1075  fun=bj_av(1)*( 1.d0-dtanh(((a0-a)/w)**2) )
1076  !fun=BJ_av(1)*dsin(pi*a) !**0.5d0
1077  fun_jb_a=ampl*fun
1078 
1079  return
1080  end
1081 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1082 
1083 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1084  subroutine promat_j (n,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,
1085  * psim,iter,kstep,fvac,psi_eav,pspl_av,psibon0)
1086 ! diffusion & equilibrium both equations solvinq
1087 ! Newton metod with linearisation
1088  implicit real*8 (a-h,o-z),integer*4 (i-n)
1089  include 'dimpl1.inc'
1090 
1091  common /com_jb/ bj_av(nrp),curfi_av(nrp)
1092  common/savt0/ psi0(nrp),fi0(nrp),f0(nrp),ri0(nrp),q0(nrp),
1093  * dpsidt(nrp),dfidt(nrp),rm0,ac0n,skcen0
1094 ! * ,kstep_prev
1095  common/com_but/ sigma(nrp),cbut_b(nrp)
1096  dimension dfdpsi(n),dpdpsi(n),f(n),flx_fi(n),psia(n),q(n)
1097  dimension wrk(512)
1098 ! save kstep_prev
1099  integer kstep_prev
1100 
1101  real*8 y(nrp,2)
1102  real*8 a(nrp,4)
1103  real*8 b(nrp,4)
1104  real*8 c(nrp,4)
1105  real*8 h(nrp,2)
1106  real*8 zn(nrp,2)
1107 
1108  real*8 psi(nrp)
1109  real*8 fi (nrp)
1110  real*8 ri (nrp)
1111  real*8 alfa22(nrp)
1112  real*8 alfa33(nrp)
1113  real*8 alp33k(nrp)
1114  real*8 ds (nrp)
1115  real*8 dsk (nrp)
1116  real*8 dv (nrp)
1117  real*8 dvk (nrp)
1118  real*8 dpsi(nrp)
1119  real*8 dfi (nrp)
1120  real*8 dri (nrp)
1121  real*8 df (nrp)
1122  real*8 fk (nrp)
1123  real*8 rik (nrp)
1124  real*8 qk (nrp)
1125  real*8 delf(nrp)
1126  real*8 dfdpsn(nrp)
1127 
1128 !---------------------------------------------------------------------
1129 
1130  epsel = 1.0d-11
1131  !dt = 1.d0
1132  dt = 1.d20
1133  cap_psi=0.d0 !1.d0
1134  cap_fi =0.d0
1135  !if(iter.lt.10) then
1136  !cap_psi=0.d0
1137  !cap_fi =0.d0
1138  !endif
1139 
1140  fbnd = fvac
1141  !call psib_pla(pspl_av)
1142  call cof_bon(cps_bon,bps_bon,dps_bon)
1143  pspl_tst=-cps_bon*(psia(n)-psia(n-1))*psim
1144  * +bps_bon*dfdpsi(n)+dps_bon*dpdpsi(n)
1145 
1146  !!!psi_bnd=psibon
1147  !psi_bnd=psi_ext+pspl_av
1148  psi_bnd=0.d0
1149 
1150 ! metric coefficients
1151  call metcof(alfa22,alfa33,ds,dv,ac0,skcen,ri,
1152  * alp33k,dsk,dvk)
1153 
1154  do i=1,n
1155  dfdpsn(i)=dfdpsi(i)
1156  psi(i)=psim*psia(i)+psi_bnd
1157  fi(i)=flx_fi(i)
1158  if(kstep.ne.kstep_prev) then
1159  psi0(i)=psim*psia(i)+psibon0
1160  fi0(i)=fi(i)
1161  f0(i) =f(i)
1162  ri0(i) =ri(i)
1163  q0(i) =q(i)
1164  !dpsidt(i)=0.d0
1165  dpsidt(i)=(psi(i)-psi0(i))/dt
1166  dfidt(i )=0.d0
1167  endif
1168  xa=dfloat(i-1)/dfloat(n-1)
1169  !sigma(i) =fun_sig(xa,i)
1170  sigma(i) =1.d0
1171  !cbut_b(i) =fun_jb(xa,i)
1172  cbut_b(i) =bj_av(i) !+fun_jb(xa,i)
1173  enddo
1174 
1175  if(iter.eq.1) then
1176  rm0=rm
1177  skcen0=skcen
1178  ac0n=ac0
1179  endif
1180 
1181 !---------------------------------------------------------------------
1182 
1183  do i=1,n
1184  wrk(i)=dsqrt(dabs(alfa22(i)))
1185  !wrk(i)=(alfa22(i)**2-alfa22(i-1)**2)/ dsk(i)
1186  enddo
1187 
1188 
1189  do i=1,n-1
1190  fi(i)=fi(n)*(i-1.d0)**2/(n-1.d0)**2
1191  enddo
1192 
1193  do i=1,n-1
1194  f(i)=alfa33(i)*(fi(i+1)-fi(i))/ds(i)
1195  ri(i)=alfa22(i)*(psi(i+1)-psi(i))/ds(i)
1196  q(i)=-(fi(i+1)-fi(i))/(psi(i+1)-psi(i))
1197  enddo
1198 
1199  f(n)=fbnd
1200 
1201  ffprim=rm*(ac0*(psi(1)-psi(2))/skcen-rm*dpdpsi(1))
1202  delf(1) = ffprim - dfdpsi(1)
1203  dfdpsi(1)=ffprim
1204  do i=2,n-1
1205  ffprim = (f(i)**2-f(i-1)**2)/(psi(i+1)-psi(i-1))
1206  delf(i) = ffprim - dfdpsi(i)
1207  dfdpsi(i) = ffprim
1208  dfdpsn(i) =alp33k(i)*(ri(i)-ri(i-1))/dsk(i)
1209  enddo
1210 
1211  alfa22n=alfa22(n)
1212  !rI(n)=dsqrt(ri(n-1)**2+alfa22n/alp33k(n)*(F(N)**2-f(n-1)**2))
1213  dfdpsi(n) = (f(n)**2-f(n-1)**2)/(psi(n)-psi(n-1))
1214 
1215 
1216 !---------------------------------------------------------------------
1217  do ks = 1,100 ! iteration cycle begins
1218 
1219  do i=2,n-1
1220  dpsi(i)=0.5d0*(psi(i+1)-psi(i-1))
1221  dfi(i) =0.5d0*(fi(i+1)-fi(i-1))
1222  dri(i) =ri(i)-ri(i-1)
1223  df(i) =f(i)-f(i-1)
1224  fk(i) =0.5d0*(f(i)+f(i-1))
1225  rik(i) =0.5d0*(ri(i)+ri(i-1))
1226  !Qk(i) =0.5d0*(q(i)+q(i-1))
1227  qk(i) =-(dsk(i)/alp33k(i))*(f(i)+f(i-1))/(psi(i+1)-psi(i-1))
1228  enddo
1229 
1230 ! matrix completing: a(i)*y(i-1)-c(i)*y(i)+b(i)*y(i+1)=-h(i)
1231  do i=2,n-1
1232 
1233  dpsdt=dpsidt(i)
1234  if(dpsdt.GE.0.d0) then
1235  capfi=-cap_fi
1236  else
1237  capfi= cap_fi
1238  endif
1239 
1240  dfldt=dfidt(i)
1241  if(dfldt.LE.0.d0) then
1242  capsi=-cap_psi
1243  else
1244  capsi= cap_psi
1245  endif
1246 
1247  d2psi=0.5d0*(psi(i+1)-2.d0*psi(i)+psi(i-1))
1248  d2fi=0.5d0*(fi(i+1)-2.d0*fi(i)+fi(i-1))
1249 
1250  a(i,1)= 0.5d0*dfidt(i)-( fk(i)+0.5d0*df(i) )*
1251  + ( alfa22(i-1)/ds(i-1) )/sigma(i) !a11
1252  + -0.5d0*dfidt(i)*capsi
1253 
1254  a(i,2)=-0.5d0*dpsidt(i)+( rik(i)+0.5d0*dri(i) )*
1255  + ( alfa33(i-1)/ds(i-1) )/sigma(i) !a12
1256  + +0.5d0*dpsidt(i)*capfi
1257 
1258  a(i,3)=-0.5d0*dri(i)/dpsi(i) + alfa22(i-1)/ds(i-1) !a21
1259 
1260  a(i,4)= qk(i)*alfa33(i-1)/ds(i-1)
1261  + +0.5d0*(df(i)/dpsi(i))*(dsk(i)/alp33k(i))*(alfa33(i-1)/ds(i-1)) !a22
1262 
1263 !---------------------------------------------------------------------
1264  b(i,1)=-0.5d0*dfidt(i)+( 0.5d0*df(i)-fk(i) )*
1265  + ( alfa22(i)/ds(i) )/sigma(i) !b11
1266  + -0.5d0*dfidt(i)*capsi
1267 
1268  b(i,2)= 0.5d0*dpsidt(i)+( rik(i)-0.5d0*dri(i) )*
1269  + ( alfa33(i)/ds(i) )/sigma(i) !b12
1270  + +0.5d0*dpsidt(i)*capfi
1271 
1272  b(i,3)= 0.5d0*dri(i)/dpsi(i)+alfa22(i)/ds(i) !b21
1273 
1274  b(i,4)= qk(i)*alfa33(i)/ds(i)
1275  + -0.5d0*(df(i)/dpsi(i))*(dsk(i)/alp33k(i))*(alfa33(i)/ds(i)) !b22
1276 
1277 !---------------------------------------------------------------------
1278  c(i,1)=-dfi(i)/dt-( fk(i)*(alfa22(i)/ds(i)+alfa22(i-1)/ds(i-1)) !c11
1279  + -0.5d0*df(i)*( alfa22(i)/ds(i)-alfa22(i-1)/ds(i-1) ) )/sigma(i) !(-0.5..)
1280  + -dfidt(i)*capsi !-capfi*d2fi/dt
1281 
1282  c(i,2)= dpsi(i)/dt+
1283  + ( 0.5d0*dri(i)*(alfa33(i-1)/ds(i-1)-alfa33(i)/ds(i)) + !c12
1284  + rik(i)*(alfa33(i)/ds(i)+alfa33(i-1)/ds(i-1)) )/sigma(i)
1285  + +dpsidt(i)*capfi !+capsi*d2psi/dt
1286 
1287  c(i,3)= alfa22(i)/ds(i)+alfa22(i-1)/ds(i-1) !c21
1288 
1289  c(i,4)= qk(i)*( alfa33(i)/ds(i)+alfa33(i-1)/ds(i-1) ) !c22
1290  + +0.5d0*(df(i)/dpsi(i))*(dsk(i)/alp33k(i))
1291  + *( alfa33(i-1)/ds(i-1)-alfa33(i)/ds(i) )
1292 !---------------------------------------------------------------------
1293  fk0 =0.5d0*(f0(i)+f0(i-1))
1294  rik0 =0.5d0*(ri0(i)+ri0(i-1))
1295  qk0 =0.5d0*(q0(i)+q0(i-1))
1296  dri0 =ri0(i)-ri0(i-1)
1297  df0 =f0(i)-f0(i-1)
1298 
1299  h(i,1)= dpsidt(i)*dfi(i)-dfidt(i)*dpsi(i)
1300  + +dpsidt(i)*capfi*d2fi-dfidt(i)*capsi*d2psi
1301 ! + -Fk(i)**2*( ri(i)/f(i)-ri(i-1)/f(i-1) )/sigma(i)
1302  + -(f(i)*f(i-1))*( ri(i)/f(i)-ri(i-1)/f(i-1) )/sigma(i)
1303  + +( cbut_b(i)*dvk(i) )/sigma(i)
1304 ! + -( Fk(i)*dri(i)-rIk(i)*df(i) )/sigma(i)
1305 ! + -( Fk0*dri0-rIk0*df0 )/sigma(i)
1306 
1307 ! h(i,2)= Qk(i)*df(i)+dri(i)-dpdpsi(i)*(dv(i)+dv(i-1))*0.5d0
1308  h(i,2)= dri(i)-dfdpsi(i)*dsk(i)/alp33k(i)-dpdpsi(i)*dvk(i)
1309 !---------------------------------------------------------------------
1310  enddo
1311 
1312 ! boundary equations
1313  b(1,1)= rm*ac0/(sigma(1)*skcen)
1314  b(1,2)= 0.d0
1315  b(1,3)= 0.d0
1316  b(1,4)= 0.d0
1317 
1318  c(1,1)= rm*ac0/(sigma(1)*skcen)-1.d0/dt !c11
1319  c(1,2)= 0.d0 !c12
1320  c(1,3)= 0.d0 !c21
1321  c(1,4)= 1.d0 !c22
1322 
1323  h(1,1)= dpsidt(1)-rm*ac0*(psi(1)-psi(2))/(sigma(1)*skcen)
1324  * +cbut_b(1)*rm**2/f(1)/sigma(1)
1325 ! * -rm0*ac0n*(psi0(1)-psi0(2))/(sigma(1)*skcen0)
1326  h(1,2)= 0.d0
1327 
1328 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1329 ! a(n,1)= cps_bon ! 0.d0 !a11
1330 ! a(n,2)= 0.d0 !a12
1331 ! a(n,3)= 0.d0 !a21
1332 ! a(n,4)=-alfa33(n-1)/ds(n-1) !a22
1333 !
1334 ! c(n,1)=cps_bon-1.d0 !c11
1335 ! c(n,2)= 0.d0 !c12
1336 ! c(n,3)= 0.d0 !c21
1337 ! c(n,4)=-alfa33(n-1)/ds(n-1) !c22
1338 !
1339 ! h(n,1)=-psi_ext+cps_bon*psi(n-1)-(cps_bon-1.d0)*psi(n) !0.d0
1340 ! h(n,2)=-fbnd+alfa33(n-1)*(fi(n)-fi(n-1))/ds(n-1) !dfi(n-1) ?
1341 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1342 
1343  !!!include 'cof_in.f'
1344 ch4astra include 'pcof_jb_in.f' !!!
1345  include 'pcof_jb_in.inc' !!!
1346 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1347  if(ks.eq.1) then
1348  do i=1,n
1349  !h(i,2)=0.d0
1350  enddo
1351  !h(1,1)=0.d0
1352  endif
1353 
1354 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1355 ! progonka
1356  call mtrx_prog(n,y,a,b,c,h)
1357 
1358 
1359  do i=1,n
1360  psi(i)= psi(i) + y(i,1) ! y(i,1) = delta(psi(i))
1361  fi(i) = fi(i) + y(i,2) ! y(i,2) = delta(fi (i))
1362  dpsidt(i)=(psi(i)-psi0(i))/dt
1363  dfidt(i )=(fi(i)-fi0(i))/dt
1364  enddo
1365 
1366  do i=1,n-1
1367  f(i)=alfa33(i)*(fi(i+1)-fi(i))/ds(i)
1368  ri(i)=alfa22(i)*(psi(i+1)-psi(i))/ds(i)
1369  q(i)=-(fi(i+1)-fi(i))/(psi(i+1)-psi(i))
1370  enddo
1371 
1372  ffprim=rm*(ac0*(psi(1)-psi(2))/skcen-rm*dpdpsi(1))
1373  delf(1) = ffprim - dfdpsi(1)
1374  dfdpsi(1)=ffprim
1375  !dfdpsi(1)=0.5d0*(dfdps+dfdpsi(1))
1376  do i=2,n-1
1377  ffprim = (f(i)**2-f(i-1)**2)/(psi(i+1)-psi(i-1))
1378  delf(i) = ffprim - dfdpsi(i)
1379  dfdpsi(i) = ffprim
1380  dfdpsn(i) =alp33k(i)*(ri(i)-ri(i-1))/dsk(i)
1381  enddo
1382  !dfdpsi(n) = dfdpsi(n-1) !<<<<<<<<<<<<<
1383  alfa22n=alfa22(n)
1384  f(n)=fbnd
1385  !rI(n)=dsqrt(ri(n-1)**2+alfa22n/alp33k(n)*(F(N)**2-f(n-1)**2))
1386  dfdpsi(n) = (f(n)**2-f(n-1)**2)/(psi(n)-psi(n-1))
1387 
1388  deltamax = 0.d0
1389  do k=1,2
1390  do i=1,n
1391  aby = dabs(y(i,k))
1392  if (aby.ge.deltamax) deltamax=aby
1393  enddo
1394  enddo
1395 
1396  do i=2,n-1
1397 
1398  dpsdt=dpsidt(i)
1399  if(dpsdt.GE.0.d0) then
1400  capfi=-cap_fi
1401  else
1402  capfi= cap_fi
1403  endif
1404 
1405  dfldt=dfidt(i)
1406  if(dfldt.LE.0.d0) then
1407  capsi=-cap_psi
1408  else
1409  capsi= cap_psi
1410  endif
1411 
1412  d2psi= (psi(i+1)-2.d0*psi(i)+psi(i-1))
1413  d2fi= (fi(i+1)-2.d0*fi(i)+fi(i-1))
1414 
1415  znv=(dpsidt(i)*(fi(i+1)-fi(i-1))-dfidt(i)*(psi(i+1)-psi(i-1)))
1416  *-(f(i)*f(i-1))*( ri(i)/f(i)-ri(i-1)/f(i-1) )*2.0d0
1417 ! *-((f(i)+f(i-1))*(ri(i)-ri(i-1))-(f(i)-f(i-1))*(ri(i)+ri(i-1)))
1418  * /sigma(i)
1419  + +dpsidt(i)*capfi*d2fi-dfidt(i)*capsi*d2psi
1420 ! *-( (f0(i)+f0(i-1))*(ri0(i)-ri0(i-1))
1421 ! * -(f0(i)-f0(i-1))*(ri0(i)+ri0(i-1)) )
1422 ! * /sigma
1423  zn(i,1)=0.5d0*znv
1424 
1425 ! zn(i,2)= (ri(i)-ri(i-1)) + 0.5d0*(q(i)+q(i-1))*(f(i)-f(i-1))
1426 ! * - dpdpsi(i)*(dv(i)+dv(i-1))*0.5d0
1427  zn(i,2)= (ri(i)-ri(i-1)) - dfdpsi(i)*dsk(i)/alp33k(i)
1428  * - dpdpsi(i)*dvk(i)
1429  enddo
1430 
1431  ddivmax = 0.d0
1432  do k=1,2
1433  do i=1,n
1434  abd = dabs(zn(i,k))
1435  if (abd.ge.ddivmax) ddivmax=abd
1436  enddo
1437  enddo
1438 
1439  if (deltamax.le.epsel) exit
1440  !if (ddivmax.le.epsel) exit
1441 
1442  ri(n)=dsqrt(ri(n-1)**2+alfa22(n)/alp33k(n)*(f(n)**2-f(n-1)**2)
1443  * +alfa22n*dpdpsi(n)*(psi(n)-psi(n-1))*dvk(n)/dsk(n) )
1444 
1445  if(ks.eq.99) then
1446  write(*,*) 'no newton iterations convergence in promat'
1447  write(*,*) 'execution was terminated'
1448  write(*,*) 'ks=',ks
1449  stop
1450  endif
1451 
1452  enddo ! iteration cycle ends
1453 !---------------------------------------------------------------------
1454 
1455  kstep_prev=kstep
1456 
1457  10 continue
1458 !---------------------------------------------------------------------
1459  return
1460  end
1461 
1462 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1463 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1464  subroutine put_key_fix(k)
1465  common /com_fixbon/ key_fixfree
1466  integer k,key_fixfree
1467  key_fixfree=k
1468  return
1469  end
cof_bon
subroutine cof_bon(cps_bon, bps_bon, dps_bon)
Definition: com_sub.f:275
fun_sig
real *8 function fun_sig(a, i)
Definition: exp_promat_e_I_f0.f:913
fun
subroutine fun(X, F)
Definition: Ev2.f:10
fun_sig_a
real *8 function fun_sig_a(a, i)
Definition: exp_promat_e_I_f0.f:994
mtrx_prog
subroutine mtrx_prog(n, y, a, b, c, f)
Definition: exp_promat_e_I_f0.f:847
promat
subroutine promat(n, dfdpsi, dpdpsi, f, flx_fi, psia, q, rm, psim, iter, kstep, fvac, psi_eav, pspl_av, psibon0, cur_mu)
Definition: exp_promat_e_I_f0.f:1
helena21_mod::dpdpsi
real(r8) function dpdpsi(psi_n)
Definition: helena21_mod.f90:122
fun_sig__
real *8 function fun_sig__(a)
Definition: exp_promat_e_I_f0.f:948
put_key_fix
subroutine put_key_fix(k)
Definition: exp_promat_e_I_f0.f:1471
fun_jb
real *8 function fun_jb(a, i)
Definition: exp_promat_e_I_f0.f:1058
fun___sig
real *8 function fun___sig(a)
Definition: exp_promat_e_I_f0.f:969
metcof
subroutine metcof(alp22, alp33, delsc, delv, ac0, skcen, cur_I,
Definition: com_sub.f:413
fun_sig_
real *8 function fun_sig_(a)
Definition: exp_promat_e_I_f0.f:927
fun_jb_a
real *8 function fun_jb_a(a, i)
Definition: exp_promat_e_I_f0.f:1070
promat_i
subroutine promat_i(n, dfdpsi, dpdpsi, f, flx_fi, psia, q, rm, psim, iter, kstep, fvac, psi_eav, pspl_av, psibon0, tok)
Definition: exp_promat_e_I_f0.f:428
promat_j
subroutine promat_j(n, dfdpsi, dpdpsi, f, flx_fi, psia, q, rm, psim, iter, kstep, fvac, psi_eav, pspl_av, psibon0, cur_mu)
Definition: exp_promat_e_I_f0.f:1088
fun_sigm
real *8 function fun_sigm(a)
Definition: exp_promat_e_I_f0.f:1037