ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
_metric.f
Go to the documentation of this file.
1 !!!this file contains the following routines:
2 !!!
3 !!!!!!!! metric
4 !!!!!!!! geom
5 !!!!!!!! funsq
6 !!!!!!!! numlin
7 !!!!!!!! matcof
8 !!!!!!!! matpla
9 !!!!!!!! matrix
10 !!!!!!!! numlin
11 !!!!!!!! prog1d
12 !!!!!!!! procof
13 !!!!!!!! extrpc
14 !!!!!!!!!!!!!!!!
15 
16  subroutine f_metric
17  !-----------------
18 
19  include 'double.inc'
20  include 'parevo.inc'
21  parameter(nkp=njlim)
22  include 'dim.inc'
23  include 'compol.inc'
24  include 'compol_add.inc'
25 
26  sqrt(xx)=dsqrt(xx)
27  !write(*,*) 'metric:enter'
28 
29  do 10 i=1,nr1
30  do 10 j=1,nt
31 
32  dlr(i,j)=sqrt( (r(i+1,j)-r(i,j))**2+(z(i+1,j)-z(i,j))**2 )
33  st(i,j)=dlr(i,j)*(r(i+1,j)+r(i,j))*0.5d0
34 
35  10 continue
36  !write(*,*) 'metric:10 '
37 
38  do 20 i=2,nr
39  do 20 j=1,nt1
40 
41  dlt(i,j)=sqrt( (r(i,j+1)-r(i,j))**2+(z(i,j+1)-z(i,j))**2 )
42  sr(i,j)=dlt(i,j)*(r(i,j+1)+r(i,j))*0.5d0
43 
44  20 continue
45  !write(*,*) 'metric:20 '
46 
47  do 30 i=1,nr1
48  do 35 j=1,nt1
49 
50  r1=r(i,j)
51  r2=r(i+1,j)
52  r3=r(i+1,j+1)
53  r4=r(i,j+1)
54 
55  r0=(r1+r2+r3+r4)*0.25d0
56 
57  r12=(r1+r2)*0.5d0
58  r23=(r3+r2)*0.5d0
59  r34=(r3+r4)*0.5d0
60  r14=(r1+r4)*0.5d0
61 
62  z1=z(i,j)
63  z2=z(i+1,j)
64  z3=z(i+1,j+1)
65  z4=z(i,j+1)
66 
67  z0=(z1+z2+z3+z4)*0.25d0
68 
69  z12=(z1+z2)*0.5d0
70  z23=(z3+z2)*0.5d0
71  z34=(z3+z4)*0.5d0
72  z14=(z1+z4)*0.5d0
73 
74  sq1(i,j)=funsq(r1,r12,r0,r14,z1,z12,z0,z14)
75  sq2(i,j)=funsq(r2,r23,r0,r12,z2,z23,z0,z12)
76  sq3(i,j)=funsq(r3,r34,r0,r23,z3,z34,z0,z23)
77  sq4(i,j)=funsq(r4,r14,r0,r34,z4,z14,z0,z34)
78 
79  s(i,j)=funsq(r1,r2,r3,r4,z1,z2,z3,z4)
80 
81  vol1(i,j)=funsq(r1,r2,r4,r4,z1,z2,z4,z4)*(r1+r2+r4)/6.d0
82  vol2(i,j)=funsq(r1,r2,r3,r3,z1,z2,z3,z3)*(r1+r2+r3)/6.d0
83  vol3(i,j)=funsq(r2,r3,r4,r4,z2,z3,z4,z4)*(r2+r3+r4)/6.d0
84  vol4(i,j)=funsq(r1,r3,r4,r4,z1,z3,z4,z4)*(r1+r3+r4)/6.d0
85 
86  !vol1(i,j)=funsq(r1,r2,r4,r4,z1,z2,z4,z4)*r0*0.5d0
87  !vol2(i,j)=funsq(r1,r2,r3,r3,z1,z2,z3,z3)*r0*0.5d0
88  !vol3(i,j)=funsq(r2,r3,r4,r4,z2,z3,z4,z4)*r0*0.5d0
89  !vol4(i,j)=funsq(r1,r3,r4,r4,z1,z3,z4,z4)*r0*0.5d0
90 
91 
92  vol(i,j)=(vol1(i,j)+vol2(i,j)+vol3(i,j)+vol4(i,j))
93 
94  dr12=r2-r1
95  dz12=z2-z1
96  dl12= (dr12**2+dz12**2)
97 
98  dr14=r4-r1
99  dz14=z4-z1
100  dl14= (dr14**2+dz14**2)
101 
102  dr23=r3-r2
103  dz23=z3-z2
104  dl23= (dr23**2+dz23**2)
105 
106  dr34=r3-r4
107  dz34=z3-z4
108  dl34= (dr34**2+dz34**2)
109 
110  cos2(i,j)=(dr12*dr23+dz12*dz23)/sqrt(dl12*dl23)
111  sin2(i,j)=1.d0-cos2(i,j)**2
112 
113  cos3(i,j)=(dr34*dr23+dz34*dz23)/sqrt(dl34*dl23)
114  sin3(i,j)=1.d0-cos3(i,j)**2
115 
116  if(i.eq.1) go to 35
117 
118  cos4(i,j)=(dr34*dr14+dz34*dz14)/sqrt(dl34*dl14)
119  sin4(i,j)=1.d0-cos4(i,j)**2
120 
121  cos1(i,j)=(dr12*dr14+dz12*dz14)/sqrt(dl12*dl14)
122  sin1(i,j)=1.d0-cos1(i,j)**2
123 
124 
125  !cos1(i,j)=0.d0
126  !sin1(i,j)=1.d0
127  !
128  !cos2(i,j)=0.d0
129  !sin2(i,j)=1.d0
130  !
131  !cos3(i,j)=0.d0
132  !sin3(i,j)=1.d0
133  !
134  !cos4(i,j)=0.d0
135  !sin4(i,j)=1.d0
136 
137  35 continue
138 
139  30 continue
140 
141 
142  return
143  end
144 
145 
146 
147  subroutine f_geom
148 
149 !!!!cylindr geometry
150 
151  include'double.inc'
152  include 'parevo.inc'
153  parameter(nkp=njlim)
154  include'dim.inc'
155  include'compol.inc'
156  include 'compol_add.inc'
157 
158  sqrt(xx)=dsqrt(xx)
159  !write(6,*) 'metric:enter'
160 
161  do 10 i=1,nr1
162  do 10 j=1,nt
163 
164  dlr(i,j)=sqrt( (r(i+1,j)-r(i,j))**2+(z(i+1,j)-z(i,j))**2 )
165  st(i,j)=dlr(i,j)*rm
166 
167  10 continue
168  !write(6,*) 'metric:10 '
169 
170  do 20 i=2,nr
171  do 20 j=1,nt1
172 
173  dlt(i,j)=sqrt( (r(i,j+1)-r(i,j))**2+(z(i,j+1)-z(i,j))**2 )
174  sr(i,j)=dlt(i,j)*rm
175 
176  20 continue
177  !write(6,*) 'metric:20 '
178 
179  do 30 i=1,nr1
180  do 35 j=1,nt1
181 
182  r1=r(i,j)
183  r2=r(i+1,j)
184  r3=r(i+1,j+1)
185  r4=r(i,j+1)
186 
187  r0=(r1+r2+r3+r4)*0.25d0
188 
189  r12=(r1+r2)*0.5d0
190  r23=(r3+r2)*0.5d0
191  r34=(r3+r4)*0.5d0
192  r14=(r1+r4)*0.5d0
193 
194  z1=z(i,j)
195  z2=z(i+1,j)
196  z3=z(i+1,j+1)
197  z4=z(i,j+1)
198 
199  z0=(z1+z2+z3+z4)*0.25d0
200 
201  z12=(z1+z2)*0.5d0
202  z23=(z3+z2)*0.5d0
203  z34=(z3+z4)*0.5d0
204  z14=(z1+z4)*0.5d0
205 
206  sq1(i,j)=funsq(r1,r12,r0,r14,z1,z12,z0,z14)
207  sq2(i,j)=funsq(r2,r23,r0,r12,z2,z23,z0,z12)
208  sq3(i,j)=funsq(r3,r34,r0,r23,z3,z34,z0,z23)
209  sq4(i,j)=funsq(r4,r14,r0,r34,z4,z14,z0,z34)
210 
211  s(i,j)=funsq(r1,r2,r3,r4,z1,z2,z3,z4)
212 
213  vol1(i,j)=funsq(r1,r2,r4,r4,z1,z2,z4,z4)*(rm+rm+rm)/6.
214  vol2(i,j)=funsq(r1,r2,r3,r3,z1,z2,z3,z3)*(rm+rm+rm)/6.
215  vol3(i,j)=funsq(r2,r3,r4,r4,z2,z3,z4,z4)*(rm+rm+rm)/6.
216  vol4(i,j)=funsq(r1,r3,r4,r4,z1,z3,z4,z4)*(rm+rm+rm)/6.
217 
218  vol(i,j)=(vol1(i,j)+vol2(i,j)+vol3(i,j)+vol4(i,j))
219 
220  dr12=r2-r1
221  dz12=z2-z1
222  dl12= (dr12**2+dz12**2)
223 
224  dr14=r4-r1
225  dz14=z4-z1
226  dl14= (dr14**2+dz14**2)
227 
228  dr23=r3-r2
229  dz23=z3-z2
230  dl23= (dr23**2+dz23**2)
231 
232  dr34=r3-r4
233  dz34=z3-z4
234  dl34= (dr34**2+dz34**2)
235 
236  cos2(i,j)=(dr12*dr23+dz12*dz23)/sqrt(dl12*dl23)
237  sin2(i,j)=1.-cos2(i,j)**2
238 
239  cos3(i,j)=(dr34*dr23+dz34*dz23)/sqrt(dl34*dl23)
240  sin3(i,j)=1.-cos3(i,j)**2
241 
242  if(i.eq.1) go to 35
243 
244  cos4(i,j)=(dr34*dr14+dz34*dz14)/sqrt(dl34*dl14)
245  sin4(i,j)=1.-cos4(i,j)**2
246 
247  cos1(i,j)=(dr12*dr14+dz12*dz14)/sqrt(dl12*dl14)
248  sin1(i,j)=1.-cos1(i,j)**2
249 
250 
251  cos1(i,j)=0.d0
252  sin1(i,j)=1.d0
253 
254  cos2(i,j)=0.d0
255  sin2(i,j)=1.d0
256 
257  cos3(i,j)=0.d0
258  sin3(i,j)=1.d0
259 
260  cos4(i,j)=0.d0
261  sin4(i,j)=1.d0
262 
263  35 continue
264 
265  30 continue
266 
267 
268  return
269  end
270 
271 
272 
273 
274 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
276 
277 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
278 
279  subroutine f_matcof
280 
281  include 'double.inc'
282  include 'parevo.inc'
283  parameter(nkp=njlim)
284  include 'dim.inc'
285  include 'compol.inc'
286  include 'compol_add.inc'
287 
288  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
289  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
290 
291 
292  do 10 i=1,nr1
293  do 10 j=1,nt1
294 
295  s12=st(i,j)
296  s14=sr(i,j)
297  s34=st(i,j+1)
298  s23=sr(i+1,j)
299 
300  if(i.ne.1) then
301 
302  a12(i,j)=( (-1.d0/s12+cos1(i,j)/s14)*vol1(i,j)/sin1(i,j) +
303  + (-1.d0/s12-cos2(i,j)/s23)*vol2(i,j)/sin2(i,j) )/s12
304 
305  a23(i,j)=( (-1.d0/s23-cos2(i,j)/s12)*vol2(i,j)/sin2(i,j) +
306  + (-1.d0/s23+cos3(i,j)/s34)*vol3(i,j)/sin3(i,j) )/s23
307 
308  a34(i,j)=( (-1.d0/s34+cos3(i,j)/s23)*vol3(i,j)/sin3(i,j) +
309  + (-1.d0/s34-cos4(i,j)/s14)*vol4(i,j)/sin4(i,j) )/s34
310 
311  a14(i,j)=( (-1.d0/s14+cos1(i,j)/s12)*vol1(i,j)/sin1(i,j) +
312  + (-1.d0/s14-cos4(i,j)/s34)*vol4(i,j)/sin4(i,j) )/s14
313 
314  a13(i,j)= (cos2(i,j)/(s12*s23))*vol2(i,j)/sin2(i,j) +
315  + (cos4(i,j)/(s34*s14))*vol4(i,j)/sin4(i,j)
316 
317  a24(i,j)= (-cos1(i,j)/(s12*s14))*vol1(i,j)/sin1(i,j) +
318  + (-cos3(i,j)/(s34*s23))*vol3(i,j)/sin3(i,j)
319 
320  else
321 
322  a12(i,j)=( (-1.d0/s12-cos2(i,j)/s23)*vol2(i,j)/sin2(i,j) )/s12
323 
324  a23(i,j)=( (-1.d0/s23-cos2(i,j)/s12)*vol2(i,j)/sin2(i,j) +
325  + (-1.d0/s23+cos3(i,j)/s34)*vol3(i,j)/sin3(i,j) )/s23
326 
327  a34(i,j)=( (-1.d0/s34+cos3(i,j)/s23)*vol3(i,j)/sin3(i,j) )/s34
328 
329 
330  a13(i,j)= (cos2(i,j)/(s12*s23))*vol2(i,j)/sin2(i,j)
331 
332  a24(i,j)= (-cos3(i,j)/(s34*s23))*vol3(i,j)/sin3(i,j)
333 
334  endif
335 
336  10 continue
337 
338  return
339  end
340 
341 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
342 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
343 
344  subroutine f_matrix
345 
346 !!!!!! IL -- number of equation (matrix line)
347 !!!!!! IM -- element number in one-dimensional array A (a(im))
348 !!!!!! IA(IL) -- number of first nonzero element in line IL
349 !!!!!! JA(IM) -- number of matrix column for element a(im)
350 
351  include 'double.inc'
352  include 'parevo.inc'
353  parameter(nkp=njlim)
354  include 'dim.inc'
355  include 'compol.inc'
356  include 'compol_add.inc'
357 
358  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
359  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
360 
361 
362 
363 
364 !!!!!!!!!!! equation for central point
365 
366  il=1
367  im=1
368 
369  ia(1)=1
370  ja(1)=1
371  a(1)=0.d0
372 
373  do 20 j=2,nt1
374 
375  im=im+1
376  ja(im)=numlin(2,j,nr,nt)
377  a(im)=a12(1,j)+a24(1,j)+a34(1,j-1)+a13(1,j-1)
378  a(1)=a(1)-a(im)
379 
380  20 continue
381 
382 !!!!!!!!!!! the main loop
383 
384  do 100 i=2,nr1
385 
386  !!!!!
387  !
388  j=2 !
389  !
390  !!!!!
391 
392  a1=a13(i-1,j-1)
393  a2=a34(i-1,j-1)+a12(i-1,j)
394  a3=a24(i-1,j)
395  a4=a23(i-1,j-1)+a14(i,j-1)
396  a6=a14(i,j)+a23(i-1,j)
397  a7=a24(i,j-1)
398  a8=a34(i,j-1)+a12(i,j)
399  a9=a13(i,j)
400 
401  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
402 
403  il=il+1
404  ia(il)=im+1
405 
406  if(i.ne.2) then
407 
408 
409 !2!cof. to (i-1,j)
410 
411  im=im+1
412  a(im)=a2
413  ja(im)=numlin(i-1,j,nr,nt)
414 
415 !3!cof. to (i-1,j+1)
416 
417  im=im+1
418  a(im)=a3
419  ja(im)=numlin(i-1,j+1,nr,nt)
420 
421 !1!cof. to (i-1,j-1)
422 
423  im=im+1
424  a(im)=a1
425  ja(im)=numlin(i-1,j-1,nr,nt)
426 
427 !5!cof. to (i,j)
428 
429  im=im+1
430  a(im)=a5
431  ja(im)=numlin(i,j,nr,nt)
432 
433 !6!cof. to (i,j+1)
434 
435  im=im+1
436  a(im)=a6
437  ja(im)=numlin(i,j+1,nr,nt)
438 
439 !4!cof. to (i,j-1)
440 
441  im=im+1
442  a(im)=a4
443  ja(im)=numlin(i,j-1,nr,nt)
444 
445  if(i.eq.nr1) go to 71
446 
447 !8!cof. to (i+1,j)
448 
449  im=im+1
450  a(im)=a8
451  ja(im)=numlin(i+1,j,nr,nt)
452 
453 !9!cof. to (i+1,j+1)
454 
455  im=im+1
456  a(im)=a9
457  ja(im)=numlin(i+1,j+1,nr,nt)
458 
459 !7!cof. to (i+1,j-1)
460 
461  im=im+1
462  a(im)=a7
463  ja(im)=numlin(i+1,j-1,nr,nt)
464 
465  71 continue
466 
467 
468  else
469 
470 
471 !1!cof. to central point
472 
473  im=im+1
474  a(im)=a1+a2+a3
475  ja(im)=numlin(1,j,nr,nt)
476 
477 !5!cof. to (2,j)
478 
479  im=im+1
480  a(im)=a5
481  ja(im)=numlin(2,j,nr,nt)
482 
483 !6!cof. to (2,j+1)
484 
485  im=im+1
486  a(im)=a6
487  ja(im)=numlin(2,j+1,nr,nt)
488 
489 !4!cof. to (2,j-1)
490 
491  im=im+1
492  a(im)=a4
493  ja(im)=numlin(2,j-1,nr,nt)
494 
495 !8!cof. to (3,j)
496 
497  im=im+1
498  a(im)=a8
499  ja(im)=numlin(3,j,nr,nt)
500 
501 !9!cof. to (3,j+1)
502 
503  im=im+1
504  a(im)=a9
505  ja(im)=numlin(3,j+1,nr,nt)
506 
507 !7!cof. to (3,j-1)
508 
509  im=im+1
510  a(im)=a7
511  ja(im)=numlin(3,j-1,nr,nt)
512 
513 
514  endif
515 
516 
517  do 110 j=3,nt2
518 
519  a1=a13(i-1,j-1)
520  a2=a34(i-1,j-1)+a12(i-1,j)
521  a3=a24(i-1,j)
522  a4=a23(i-1,j-1)+a14(i,j-1)
523  a6=a14(i,j)+a23(i-1,j)
524  a7=a24(i,j-1)
525  a8=a34(i,j-1)+a12(i,j)
526  a9=a13(i,j)
527 
528  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
529 
530  il=il+1
531  ia(il)=im+1
532 
533  if(i.ne.2) then
534 
535 
536 !1!cof. to (i-1,j-1)
537 
538  im=im+1
539  a(im)=a1
540  ja(im)=numlin(i-1,j-1,nr,nt)
541 
542 !2!cof. to (i-1,j)
543 
544  im=im+1
545  a(im)=a2
546  ja(im)=numlin(i-1,j,nr,nt)
547 
548 !3!cof. to (i-1,j+1)
549 
550  im=im+1
551  a(im)=a3
552  ja(im)=numlin(i-1,j+1,nr,nt)
553 
554 !3!cof. to (i,j-1)
555 
556  im=im+1
557  a(im)=a4
558  ja(im)=numlin(i,j-1,nr,nt)
559 
560 !5!cof. to (i,j)
561 
562  im=im+1
563  a(im)=a5
564  ja(im)=numlin(i,j,nr,nt)
565 
566 !6!cof. to (i,j+1)
567 
568  im=im+1
569  a(im)=a6
570  ja(im)=numlin(i,j+1,nr,nt)
571 
572  if(i.eq.nr1) go to 7
573 
574 !7!cof. to (i+1,j-1)
575 
576  im=im+1
577  a(im)=a7
578  ja(im)=numlin(i+1,j-1,nr,nt)
579 
580 !8!cof. to (i+1,j)
581 
582  im=im+1
583  a(im)=a8
584  ja(im)=numlin(i+1,j,nr,nt)
585 
586 !9!cof. to (i+1,j+1)
587 
588  im=im+1
589  a(im)=a9
590  ja(im)=numlin(i+1,j+1,nr,nt)
591 
592  7 continue
593 
594 
595  else
596 
597 
598 !1!cof. to central point
599 
600  im=im+1
601  a(im)=a1+a2+a3
602  ja(im)=numlin(1,j,nr,nt)
603 
604 !2!cof. to (2,j-1)
605 
606  im=im+1
607  a(im)=a4
608  ja(im)=numlin(2,j-1,nr,nt)
609 
610 !3!cof. to (2,j)
611 
612  im=im+1
613  a(im)=a5
614  ja(im)=numlin(2,j,nr,nt)
615 
616 !4!cof. to (2,j+1)
617 
618  im=im+1
619  a(im)=a6
620  ja(im)=numlin(2,j+1,nr,nt)
621 
622 !5!cof. to (3,j-1)
623 
624  im=im+1
625  a(im)=a7
626  ja(im)=numlin(3,j-1,nr,nt)
627 
628 !6!cof. to (3,j)
629 
630  im=im+1
631  a(im)=a8
632  ja(im)=numlin(3,j,nr,nt)
633 
634 !7!cof. to (3,j+1)
635 
636  im=im+1
637  a(im)=a9
638  ja(im)=numlin(3,j+1,nr,nt)
639 
640 
641  endif
642 
643  110 continue
644 
645  !!!!!!!
646  !
647  j=nt1 !
648  !
649  !!!!!!!
650 
651  a1=a13(i-1,j-1)
652  a2=a34(i-1,j-1)+a12(i-1,j)
653  a3=a24(i-1,j)
654  a4=a23(i-1,j-1)+a14(i,j-1)
655  a6=a14(i,j)+a23(i-1,j)
656  a7=a24(i,j-1)
657  a8=a34(i,j-1)+a12(i,j)
658  a9=a13(i,j)
659 
660  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
661 
662  il=il+1
663  ia(il)=im+1
664 
665  if(i.ne.2) then
666 
667 
668 !3!cof. to (i-1,j+1)
669 
670  im=im+1
671  a(im)=a3
672  ja(im)=numlin(i-1,j+1,nr,nt)
673 
674 !1!cof. to (i-1,j-1)
675 
676  im=im+1
677  a(im)=a1
678  ja(im)=numlin(i-1,j-1,nr,nt)
679 
680 !2!cof. to (i-1,j)
681 
682  im=im+1
683  a(im)=a2
684  ja(im)=numlin(i-1,j,nr,nt)
685 
686 !6!cof. to (i,j+1)
687 
688  im=im+1
689  a(im)=a6
690  ja(im)=numlin(i,j+1,nr,nt)
691 
692 !4!cof. to (i,j-1)
693 
694  im=im+1
695  a(im)=a4
696  ja(im)=numlin(i,j-1,nr,nt)
697 
698 !5!cof. to (i,j)
699 
700  im=im+1
701  a(im)=a5
702  ja(im)=numlin(i,j,nr,nt)
703 
704  if(i.eq.nr1) go to 72
705 
706 !9!cof. to (i+1,j+1)
707 
708  im=im+1
709  a(im)=a9
710  ja(im)=numlin(i+1,j+1,nr,nt)
711 
712 !7!cof. to (i+1,j-1)
713 
714  im=im+1
715  a(im)=a7
716  ja(im)=numlin(i+1,j-1,nr,nt)
717 
718 !8!cof. to (i+1,j)
719 
720  im=im+1
721  a(im)=a8
722  ja(im)=numlin(i+1,j,nr,nt)
723 
724  72 continue
725 
726 
727  else
728 
729 
730 !1!cof. to central point
731 
732  im=im+1
733  a(im)=a1+a2+a3
734  ja(im)=numlin(1,j,nr,nt)
735 
736 !6!cof. to (2,j+1)
737 
738  im=im+1
739  a(im)=a6
740  ja(im)=numlin(2,j+1,nr,nt)
741 
742 !4!cof. to (2,j-1)
743 
744  im=im+1
745  a(im)=a4
746  ja(im)=numlin(2,j-1,nr,nt)
747 
748 !5!cof. to (2,j)
749 
750  im=im+1
751  a(im)=a5
752  ja(im)=numlin(2,j,nr,nt)
753 
754 !9!cof. to (3,j+1)
755 
756  im=im+1
757  a(im)=a9
758  ja(im)=numlin(3,j+1,nr,nt)
759 
760 !7!cof. to (3,j-1)
761 
762  im=im+1
763  a(im)=a7
764  ja(im)=numlin(3,j-1,nr,nt)
765 
766 !8!cof. to (3,j)
767 
768  im=im+1
769  a(im)=a8
770  ja(im)=numlin(3,j,nr,nt)
771 
772 
773  endif
774 
775  100 continue
776  !write(6,*) 'matrix:il=',il
777 
778  neq=il
779 
780  !write(6,*) 'number of equations:',neq,neqp
781  !write(6,*) 'im=',im
782  !write(6,*) 'lp=',lp
783 
784  il=il+1
785  ia(il)=im+1
786 
787 
788 
789 
790 
791 c open(1,file='matr')
792 c
793 c do 111 kl=1,neq
794 c
795 c ic1=ia(kl)
796 c ic2=ia(kl+1)-1
797 c
798 c write(1,*) 'il=',kl
799 c write(1,*) (ja(km),km=ic1,ic2)
800 c write(1,*) (a(km),km=ic1,ic2)
801 c
802 c111 continue
803 
804 
805  if( (itin/nitdel*nitdel+nitbeg) .EQ. itin ) then
806 
807  do km=1,im
808 
809  aop0(km)=a(km)
810 
811  enddo
812 
813  !write(6,*) 'matrix:aop0:itin,nitdel',itin,nitdel
814 
815  elseif(itin.gt.nitbeg) then
816 
817  do km=1,im
818 
819  daop(km)=a(km)-aop0(km)
820 
821  enddo
822  !write(6,*) 'matrix:daop:itin',itin
823 
824  endif
825 
826 
827 
828 
829  return
830  end
831 
832 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
833 
834 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
835 
836  subroutine f_procof(icq,cur_mu)
837 
838  include 'double.inc'
839  include 'parevo.inc'
840  parameter(nkp=njlim)
841  include 'dim.inc'
842  include 'compol.inc'
843  include 'compol_add.inc'
844 
845  tok_mu=tok*amu0
846 
847  call bongri
848  call f_psib_ext(psex_av)
849  call f_psib_pla(pspl_av)
850  call cof_bon(cps_bon,bps_bon,dps_bon)
851  psi_bon=psi_eav+pspl_av
852  psibon=psi_bon
853  psi_eav=psex_av
854 
855 
856 
857 
858  if(ngav.eq.1) then
859  call f_procof_i(icq)
860  elseif(ngav.eq.2) then
861  call f_procof_fl(icq)
862  elseif(ngav.eq.-10 .AND. erru.lt.5.d-3) then
863  call promat_j(iplas,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psipla,itin,
864  * kstep,fvac,psi_eav,pspl_av,psibon0,cur_mu )
865  elseif(ngav.eq.-1 .AND. erru.lt.5.d-3) then
866  call promat(iplas,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psipla,itin,
867  * kstep,fvac,psi_eav,pspl_av,psibon0,cur_mu )
868  elseif(ngav.eq.-2 .AND. erru.lt.5.d-2) then
869  call promat(iplas,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psipla,itin,
870  * kstep,fvac,psi_eav,pspl_av,psibon0,cur_mu )
871  elseif(ngav.eq.-3 .AND. erru.lt.5.d-3) then
872  call promat_i(iplas,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psipla,itin,
873  * kstep,fvac,psi_eav,pspl_av,psibon0,tok_mu )
874  endif
875 
876  return
877  end
878 
879 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
880 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
881 
882  subroutine f_procof_fl(icq)
883 
884  include 'double.inc'
885  include 'parevo.inc'
886  parameter(nkp=njlim)
887  include 'dim.inc'
888  include 'compol.inc'
889  include 'compol_add.inc'
890 
891  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
892  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
893 
894  common/compsf/ psf,sqtor
895 
896  dimension amn(nrp),a0(nrp),apl(nrp),capp(nrp)
897  dimension avrc(nrp),delsc(nrp),avrk(nrp),delsk(nrp),delv(nrp)
898  dimension bmn(nrp),b0(nrp),bpl(nrp),wrk1(nrp),wrk2(nrp)
899  dimension psf(nrp),rhs(nrp),sqtor(nrp)
900 
901  sqrt(xx)=dsqrt(xx)
902 
903  psimx=psip+psipla
904 
905  do 30 i=1,iplas1
906  zdelsc=0.d0
907  zavrc=0.d0
908  do 35 j=2,nt1
909 
910  r1=r(i,j)
911  r2=r(i+1,j)
912  r3=r(i+1,j+1)
913  r4=r(i,j+1)
914 
915  r0=(r1+r2+r3+r4)*0.25d0
916 
917  zdelsc=zdelsc+s(i,j)
918  zavrc=zavrc+s(i,j)/r0
919 
920  35 continue
921 
922  avrc(i)=zavrc/zdelsc
923  delsc(i)=zdelsc
924 
925  30 continue
926 
927  sqtor(1)=0.d0
928  do 33 i=2,iplas
929  sqtor(i)=sqtor(i-1)+delsc(i-1)
930  33 continue
931 
932  do 40 i=2,iplas
933 
934  zdelv=0.d0
935  zdelsk=0.d0
936  zavrk=0.d0
937  zapro=0.d0
938 
939  do 45 j=2,nt1
940 
941  if(i.ne.iplas) then
942  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
943  r0=r(i,j)
944  else
945  sqk=sq2(i-1,j)+sq3(i-1,j-1)
946  r0=(3.d0*r(i,j)+r(i-1,j))*0.25d0
947  endif
948 
949  zdelsk=zdelsk+sqk
950  zdelv=zdelv +sqk*r0
951  zavrk=zavrk+sqk/r0
952  zapro=zapro+aex(i,j)
953  45 continue
954 
955  avrk(i)=zavrk/zdelsk
956  delsk(i)=zdelsk
957  delv(i)=zdelv
958  !rhs(i-1)=-zapro
959  rhs(i-1)=0.d0
960 
961  40 continue
962 
963 
964  do 10 i=2,iplas1
965 
966  samn=0.d0
967  sa0 =0.d0
968  sapl=0.d0
969 
970  do 110 j=2,nt1
971 
972  a1=a13(i-1,j-1)
973  a2=a34(i-1,j-1)+a12(i-1,j)
974  a3=a24(i-1,j)
975  a4=a23(i-1,j-1)+a14(i,j-1)
976  a6=a14(i,j)+a23(i-1,j)
977  a7=a24(i,j-1)
978  a8=a34(i,j-1)+a12(i,j)
979  a9=a13(i,j)
980  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
981 
982  zamn=a1+a2+a3
983  za0 =a4+a5+a6
984  zapl=a7+a8+a9
985 
986  samn=samn+zamn
987  sa0 =sa0 +za0
988  sapl=sapl+zapl
989 
990  110 continue
991 
992  qmn=-q(i-1)/(avrc(i-1)*delsc(i-1))
993  qpl=-q(i) / (avrc(i)*delsc(i))
994  qk=0.5d0*(qmn+qpl)
995 
996  afmn=qk*qmn*avrk(i)*delsk(i)
997  afpl=qk*qpl*avrk(i)*delsk(i)
998 
999  af0=-(afmn+afpl)
1000 
1001  amn(i)=samn - afmn
1002  a0(i) =sa0 - af0
1003  apl(i)=sapl - afpl
1004 
1005  capp(i)=sapl*delsc(i)
1006 
1007  if(i.eq.iplas-1) then
1008  cappl=sapl*delsc(i)
1009  capmn=samn*delsc(i-1)
1010  endif
1011 
1012  10 continue
1013 
1014 
1015  npro=iplas-2
1016 
1017  do 100 i=1,npro
1018 
1019  if(i.ne.1) then
1020  bmn(i)=amn(i+1)
1021  else
1022  bmn(i)=0.d0
1023  endif
1024 
1025  b0(i)=-a0(i+1)
1026 
1027  if(i.ne.npro) then
1028  bpl(i)=apl(i+1)
1029  else
1030  bpl(i)=0.d0
1031  endif
1032 
1033  rhs(i)=-dpdpsi(i+1)*delv(i+1)+rhs(i)
1034 
1035  100 continue
1036 
1037  rhs(1)=rhs(1)+amn(2)*psimx
1038  rhs(npro)=rhs(npro)+apl(iplas-1)*psip
1039 
1040  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1041 
1042  call prog1d(npro,psf(2),bpl,bmn,b0,rhs,wrk1,wrk2)
1043 
1044  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1045 
1046  psf(1)=psimx
1047  psf(iplas)=psip
1048 
1049  do 300 i=1,iplas1
1050 
1051  f(i)=-q(i)*(psf(i+1)-psf(i))/(delsc(i)*avrc(i))
1052 
1053  300 continue
1054 
1055 
1056 
1057  do 400 i=2,iplas1
1058 
1059  qmn=-q(i-1)/(avrc(i-1)*delsc(i-1))
1060  qpl=-q(i)/(avrc(i)*delsc(i))
1061  qk=0.5d0*(qmn+qpl)
1062 
1063  dfdpsi(i)=qk*(f(i)-f(i-1))
1064 
1065  400 continue
1066 
1067 !!! equation in central node
1068 
1069  ac0=0.d0
1070  skcen=0.d0
1071 
1072  do 20 j=2,nt1
1073 
1074  acj=a12(1,j)+a24(1,j)+a34(1,j-1)+a13(1,j-1)
1075  ac0=ac0-acj
1076 
1077  skcen=skcen+sq1(1,j)+sq4(1,j-1)
1078 
1079  20 continue
1080 
1081  aexcen=aex(1,2)
1082 
1083  !dfdpsi(1)=((ac0*(psimx-psf(2))-aexcen)/skcen-rm*dpdpsi(1))*rm
1084  dfdpsi(1)=((ac0*(psimx-psf(2)) )/skcen-rm*dpdpsi(1))*rm
1085 
1086 !!! equation on bound node!
1087 !!!!!!!!!!!!!!!!!!!!!!!!!!!
1088  i=iplas
1089 
1090  samn=0.d0
1091  sa0 =0.d0
1092 
1093  do 50 j=2,nt1
1094 
1095  a1=a13(i-1,j-1)
1096  a2=a34(i-1,j-1)+a12(i-1,j)
1097  a3=a24(i-1,j)
1098 
1099  zamn=a1+a2+a3
1100  za0 =-zamn
1101 
1102  samn=samn+zamn
1103  sa0 =sa0 +za0
1104 
1105  50 continue
1106 
1107  !cappa=cappl+delsk(i)*(cappl-capmn)/delsk(i-1)
1108 
1109  call extrpc( sqtor(i),capp(i),
1110  * sqtor(i),sqtor(i-1),sqtor(i-2),sqtor(i-3),
1111  * capp(i-1) ,capp(i-2) ,capp(i-3) )
1112 
1113 ! call extrpc( psia(i),capp(i),
1114 ! * psia(i),psia(i-1),psia(i-2),psia(i-3),
1115 ! * capp(i-1) ,capp(i-2) ,capp(i-3) )
1116 
1117  cappa=capp(i)
1118 
1119 
1120  !!!!!!!!!!!!!!!!!!!!!!!!
1121 
1122  if(icq.eq.0) then
1123 
1124  dpsids=(
1125  * samn*(psip-psf(i-1))+avrk(i)*delsk(i)*dfdpsi(i)
1126  * +delv(i)*dpdpsi(i)
1127  * ) /cappa
1128 
1129  f(i)=sqrt( f(i-1)**2+dfdpsi(i)*(psip-psf(i-1)) )
1130 
1131  q(i)=-f(i)*avrk(i)/dpsids
1132  write(*,*) 'q(iplas)=',q(i)
1133 ! pause ' '
1134  endif
1135  !q(i)=q(i-1)+delsk(i)*(q(i-1)-q(i-2))/delsk(i-1)
1136 
1137  !!!!!!!!!!!!!!!!!!!!!!!!
1138 
1139  qn=-q(i)/avrk(i)
1140  qmn=-q(i-1)/avrc(i-1)
1141 
1142  !Q14= 0.5d0*(Qn+Qmn)
1143  q14=-0.5d0*(q(i)+q(i-1))
1144  avr14=avrk(i)
1145 
1146 ! f(i)=(samn*(psip-psf(i-1))+delv(i)*dpdpsi(i)-Q14*avr14*f(i-1))
1147 ! * /( cappa/Qn-Q14*avr14 )
1148 
1149  f(i)=(samn*(psip-psf(i-1))+delv(i)*dpdpsi(i)-q14*f(i-1))
1150  * /( cappa/qn-q14 )
1151 
1152 ! dfdpsi(i)=Q14*(f(i)-f(i-1))/delsk(i)
1153  dfdpsi(i)=(q14/avr14)*(f(i)-f(i-1))/delsk(i)
1154 
1155 
1156  psiai=psia(i)-0.25d0*(psia(i)- psia(i-1))
1157 
1158  call extrp2( psiai,dfdpsi(i),
1159  * psia(i-3),psia(i-2),psia(i-1),
1160  * dfdpsi(i-3) ,dfdpsi(i-2) ,dfdpsi(i-1) )
1161 
1162  f(i)=dsqrt(f(i-1)**2+dfdpsi(i)*(psip- psf(i-1)))
1163 
1164 
1165 
1166 
1167  ! open(1,file='ddps1.pr')
1168  ! do i=1,iplas
1169  ! write(1,*) dfdpsi(i),f(i),i
1170  ! enddo
1171  ! write(1,*) 'q(iplas)=',q(iplas)
1172  ! close(1)
1173 
1174  ! open(1,file='ddp.wr')
1175  ! write(1,*) iplas
1176  ! write(1,*) (q(i),i=1,iplas)
1177  ! write(1,*) (f(i),i=1,iplas)
1178  ! write(1,*) (dfdpsi(i),i=1,iplas)
1179  ! write(1,*) (dpdpsi(i),i=1,iplas)
1180  !close(1)
1181  return
1182  end
1183 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1184  subroutine f_procof_i(icq)
1185 
1186  include 'double.inc'
1187  include 'parevo.inc'
1188  parameter(nkp=njlim)
1189  include 'dim.inc'
1190  include 'compol.inc'
1191  include 'compol_add.inc'
1192 
1193  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
1194  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
1195 
1196  common/compsf/ psf(nrp),sqtor(nrp)
1197 
1198  dimension amn(nrp),a0(nrp),apl(nrp),qcapp(nrp)
1199  dimension avrc(nrp),delsc(nrp),avrk(nrp),delsk(nrp),delv(nrp)
1200  dimension bmn(nrp),b0(nrp),bpl(nrp),wrk1(nrp),wrk2(nrp)
1201  dimension rhs(nrp)
1202  dimension dfdpsn(nrp),vw(nrp)
1203 
1204  sqrt(xx)=dsqrt(xx)
1205 
1206  do i=1,iplas
1207  dfdpsn(i)=dfdpsi(i)
1208  psf(i)=psia(i)*(psim-psip)+psip
1209  enddo
1210 
1211 
1212 c!!!!!!! dS(i+1/2) and <1/R>(i+1/2) !!!!!!!!!!!
1213 
1214  do 30 i=1,iplas1
1215  zdelsc=0.d0
1216  zavrc=0.d0
1217  do 35 j=2,nt1
1218 
1219  r1=r(i,j)
1220  r2=r(i+1,j)
1221  r3=r(i+1,j+1)
1222  r4=r(i,j+1)
1223 
1224  r0=(r1+r2+r3+r4)*0.25d0
1225 
1226  zdelsc=zdelsc+s(i,j)
1227  zavrc=zavrc+s(i,j)/r0
1228 
1229  35 continue
1230 
1231  avrc(i)=zavrc/zdelsc
1232  delsc(i)=zdelsc
1233 
1234  30 continue
1235 
1236 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1237 
1238  sqtor(1)=0.d0
1239  do 33 i=2,iplas
1240  sqtor(i)=sqtor(i-1)+delsc(i-1)
1241  33 continue
1242 
1243 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1244 c!!!!!!! dS(i), dV(i) and <1/R>(i) !!!!!!!!!!!
1245 
1246  do 40 i=2,iplas
1247 
1248  zdelv=0.d0
1249  zdelsk=0.d0
1250  zavrk=0.d0
1251 
1252  do 45 j=2,nt1
1253 
1254  if(i.ne.iplas) then
1255  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
1256  r0=r(i,j)
1257  else
1258  sqk=sq2(i-1,j)+sq3(i-1,j-1)
1259  r0=(3.d0*r(i,j)+r(i-1,j))*0.25d0
1260  endif
1261 
1262  zdelsk=zdelsk+sqk
1263  zdelv=zdelv +sqk*r0
1264  zavrk=zavrk+sqk/r0
1265 
1266  45 continue
1267 
1268  avrk(i)=zavrk/zdelsk
1269  delsk(i)=zdelsk
1270  delv(i)=zdelv
1271  rhs(i-1)=0.d0 !right hand side for averaged equation
1272 
1273  40 continue
1274 
1275 
1276  do 10 i=2,iplas1
1277 
1278  samn=0.d0
1279  sa0 =0.d0
1280  sapl=0.d0
1281 
1282  do 110 j=2,nt1
1283 
1284  a1=a13(i-1,j-1)
1285  a2=a34(i-1,j-1)+a12(i-1,j)
1286  a3=a24(i-1,j)
1287  a4=a23(i-1,j-1)+a14(i,j-1)
1288  a6=a14(i,j)+a23(i-1,j)
1289  a7=a24(i,j-1)
1290  a8=a34(i,j-1)+a12(i,j)
1291  a9=a13(i,j)
1292  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
1293 
1294  zamn=a1+a2+a3
1295  za0 =a4+a5+a6
1296  zapl=a7+a8+a9
1297 
1298  samn=samn+zamn
1299  sa0 =sa0 +za0
1300  sapl=sapl+zapl
1301 
1302  110 continue
1303 
1304  qmn=-q(i-1)/(avrc(i-1)*delsc(i-1))
1305  qpl=-q(i) / (avrc(i)*delsc(i))
1306  qk=(qmn*dpsda(i-1)+qpl*dpsda(i))/(dpsda(i-1)+dpsda(i))
1307 
1308  afmn=qk*qmn*avrk(i)*delsk(i)
1309  afpl=qk*qpl*avrk(i)*delsk(i)
1310 
1311  af0=-(afmn+afpl)
1312 
1313  amn(i)=samn !- afmn
1314  a0(i) =sa0 !- af0
1315  apl(i)=sapl !- afpl
1316 
1317  qcapp(i)=qpl/sapl
1318 
1319  if(i.eq.iplas-1) then
1320  qcappl=qpl/sapl
1321  qcapmn=qmn/samn
1322  endif
1323 
1324  10 continue
1325 
1326  736 continue
1327 
1328  if(itin.eq.1) then
1329 
1330  errdf=0.d0
1331 
1332  do i=2,iplas1
1333 
1334  !plcu=apl(i)*(psf(i+1)-psf(i))
1335  !write(6,*) 'tok',plcu,i
1336  !pause ' pause'
1337 
1338  zff=amn(i)*psf(i-1)+a0(i)*psf(i)+apl(i)*psf(i+1)
1339 
1340  zdfdps=( zff-dpdpsi(i)*delv(i) )/(avrk(i)*delsk(i))
1341  !znvz=zff-dpdpsi(i)*delv(i)-dfdpsn(i)*avrk(i)*delsk(i)
1342  dfdpsi(i)=zdfdps
1343  deldf=dabs(zdfdps-dfdpsn(i))
1344  !write(6,*) 'deldf',deldf,dfdpsn(i),i
1345  !write(6,*) 'nev.L(psi)',znvz,zff
1346  !pause ' pause'
1347  errdf=dmax1(deldf,errdf)
1348  dfdpsn(i)=dfdpsi(i)
1349 
1350  enddo
1351 
1352  !write(6,*) 'errdf',errdf
1353  !pause ' pause'
1354 
1355  endif
1356  !! return
1357 
1358  !!!!!!!!!!
1359 ! i=iplas
1360 ! !!!!!!!!!!
1361 !
1362 ! zff=tok-apl(i-1)*(psf(i)-psf(i-1))
1363 !
1364 ! dfdpsi(i)=(zff-dpdpsi(i)*delv(i))/(avrk(i)*delsk(i))
1365 ! deldf=dabs(dfdpsi(i)-dfdpsn(i))
1366 ! !write(6,*) 'deldf*',deldf,i
1367 ! !write(6,*) '*',dfdpsi(i),dfdpsn(i)
1368 ! !pause ' pause'
1369 !
1370 !
1371 !
1372 ! Qcappa=Qcappl+delsk(i)*(Qcappl-Qcapmn)/delsk(i-1)
1373 !
1374 ! call extrpc( sqtor(i),Qcapp(i),
1375 ! * sqtor(i),sqtor(i-1), sqtor(i-2), sqtor(i-3),
1376 ! * Qcapp(i-1) ,Qcapp(i-2) ,Qcapp(i-3) )
1377 !
1378 ! call extrpc( psia(i),Qcapp(i),
1379 ! * psia(i),psia(i-1), psia(i-2), psia(i-3),
1380 ! * Qcapp(i-1) ,Qcapp(i-2) ,Qcapp(i-3) )
1381 !
1382 ! !ff2n=f(i-1)**2+dfdpsn(i)*(psf(i)-psf(i-1))
1383 !
1384 ! !f2n=(tok*Qcappa)**2
1385 ! f2n=(tok*Qcapp(i))**2
1386 !
1387 ! !write(6,*) 'fnold fn',ff2n,f2n
1388 !
1389 ! f(i)=dsqrt(f2n)
1390 ! f(i-1)=dsqrt( f2n-dfdpsi(i)*(psf(i)-psf(i-1)) )
1391 !
1392  errdf=0.d0
1393 
1394  f2n1=(tok*qcapp(iplas-1))**2
1395  f(iplas-1)=dsqrt(f2n1)
1396 
1397  do i=iplas-1,2,-1
1398  f(i-1)=dsqrt( f(i)**2-dfdpsi(i)*(psf(i+1)-psf(i-1)) )
1399  enddo
1400 
1401  psf(iplas)=0.d0
1402 
1403  do i=iplas1,1,-1
1404  dftor=f(i)*delsc(i)*avrc(i)
1405  delpsi=dftor/q(i)
1406  psf(i)=psf(i+1)+delpsi
1407  enddo
1408 
1409  do i=2,iplas1
1410 
1411  !plcu=apl(i)*(psf(i+1)-psf(i))
1412  !write(6,*) 'tok',plcu,i
1413  !pause ' pause'
1414 
1415  zff=amn(i)*psf(i-1)+a0(i)*psf(i)+apl(i)*psf(i+1)
1416 
1417  dfdpsi(i)=( zff-dpdpsi(i)*delv(i) )/(avrk(i)*delsk(i))
1418  deldf=dabs(dfdpsi(i)-dfdpsn(i))
1419  !write(6,*) 'deldf',deldf,errdf,i
1420  !pause ' pause'
1421  errdf=dmax1(deldf,errdf)
1422  dfdpsn(i)=dfdpsi(i)
1423 
1424  enddo
1425  !write(6,*) 'errdf',errdf
1426  !pause ' pause'
1427 
1428  if(errdf.gt.1.d-7) go to 736 !!!------***
1429 
1430  i=iplas !!!<<<<<<<<<<<<<
1431 
1432  call extrp2( psia(i),dfdpsi(i),
1433  * psia(i-3),psia(i-2),psia(i-1),
1434  * dfdpsi(i-3) ,dfdpsi(i-2) ,dfdpsi(i-1) )
1435 
1436 
1437 
1438 
1439 !!!!!!!!!!///////////////////////////////////////////
1440 
1441 !!! equation in central node
1442 
1443  ac0=0.d0
1444  skcen=0.d0
1445 
1446  do 20 j=2,nt1
1447 
1448  acj=a12(1,j)+a24(1,j)+a34(1,j-1)+a13(1,j-1)
1449  ac0=ac0-acj
1450 
1451  skcen=skcen+sq1(1,j)+sq4(1,j-1)
1452 
1453  20 continue
1454 
1455 
1456  dfdpsi(1)=( ac0*(psf(1)-psf(2))/skcen-rm*dpdpsi(1))*rm
1457 
1458 
1459  return
1460  end
1461 
1462 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1463 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
extrpc
subroutine extrpc(X, F, X0, X1, X2, X3, F1, F2, F3)
Definition: com_sub.f:699
funsq
REAL *8 function funsq(R1, R2, R3, R4, Z1, Z2, Z3, Z4)
Definition: scu.f:803
cof_bon
subroutine cof_bon(cps_bon, bps_bon, dps_bon)
Definition: com_sub.f:275
f_procof_i
subroutine f_procof_i(icq)
Definition: _metric.f:1184
prog1d
subroutine prog1d(M, U, A, B, C, F, ALF, BET)
Definition: scu.f:822
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
numlin
function numlin(i, j, nr, nt)
Definition: com_sub.f:1040
extrp2
subroutine extrp2(X, F, X1, X2, X3, F1, F2, F3)
Definition: com_sub.f:719
f_psib_ext
subroutine f_psib_ext(psex_av)
Definition: _rig.f:676
bongri
subroutine bongri
Definition: com_sub.f:242
f_metric
subroutine f_metric
Definition: _metric.f:16
f_matcof
subroutine f_matcof
Definition: _metric.f:279
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
f_geom
subroutine f_geom
Definition: _metric.f:147
f_matrix
subroutine f_matrix
Definition: _metric.f:344
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
f_procof_fl
subroutine f_procof_fl(icq)
Definition: _metric.f:882
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
f_psib_pla
subroutine f_psib_pla(pspl_av)
Definition: _rig.f:708
f_procof
subroutine f_procof(icq, cur_mu)
Definition: _metric.f:836