ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
com_sub.f
Go to the documentation of this file.
1 !!!!Spider subroutines common for fix- and free-boundary
2 
3 
4 
5 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
6 
7  subroutine grdef(igdf)
8 
9  include 'double.inc'
10  include 'dim.inc'
11  include 'compol.inc'
12 
13  dimension h_fi(nrp),fia(nrp)
14 
15  if(igdf.eq.0)then
16 
17  do i=1,iplas
18  psia(i)=(iplas-i)/(iplas-1.d0)
19  enddo
20 
21  do i=1,iplas1
22  dpsda(i)=-1.d0
23  enddo
24 
25  elseif(igdf.eq.1)then
26 
27  do i=1,iplas
28  psia(i)=1.d0-((i-1)/(iplas-1.d0))**2
29  enddo
30 
31  do i=1,iplas1
32  dpsda(i)=(1-2*i)/(iplas-1.d0)
33  enddo
34 
35  elseif(igdf.eq.2)then
36 
37  psia(1)=0.d0
38 
39  do i=2,iplas
40  znam=q(i-1)*(iplas-1)**2
41  psia(i)=psia(i-1)+(2*i-3)/znam
42  enddo
43 
44  ! write(*,*) 'psia(iplas)',psia(iplas)
45  !pause 'pause'
46 
47  psnor=psia(iplas)
48 
49  do i=1,iplas
50  psia(i)=1.d0-psia(i)/psnor
51  enddo
52 
53  do i=1,iplas1
54  dpsda(i)=-(2*i-1)/( (iplas-1.d0)*q(i)*psnor )
55  enddo
56 
57  elseif(igdf.eq.3)then
58 
59  zn_gp=0.95d0
60  fia(1)=0.d0
61  psia(1)=0.d0
62  do i=1,iplas-1
63  if( i.lt.iplas/2 ) then
64  h_fi(i)=1.d0
65  elseif( i.ge.iplas/2 ) then
66  h_fi(i)=h_fi(i-1)*zn_gp
67  endif
68  fia(i+1)=fia(i)+h_fi(i)
69  psia(i+1)=psia(i)+(fia(i+1)**2-fia(i)**2)/q(i)
70  enddo
71 
72  psnor=psia(iplas)
73 
74  do i=1,iplas
75  psia(i)=1.d0-psia(i)/psnor
76  enddo
77 
78  endif
79 
80  return
81  end
82 
83 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
84 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
85 
86  subroutine axdef(rma,zma,psima,dpm)
87 
88  include 'double.inc'
89  include 'parevo.inc'
90  parameter(nkp=njlim)
91  include 'dim.inc'
92  parameter(nshp=ntp+1)
93  include 'compol.inc'
94 
95  dimension xs(nshp),ys(nshp),fun(nshp),dpm(*)
96 
97  rm0=rm
98  zm0=zm
99  psim0=psi(1,1)
100 
101  nsh=1
102 
103  xs(nsh)=rm0
104  ys(nsh)=zm0
105  fun(nsh)=psim0
106 
107  do j=2,nt1
108 
109  nsh=nsh+1
110  xs(nsh)=r(2,j)
111  ys(nsh)=z(2,j)
112  fun(nsh)=psi(2,j)
113 
114  enddo
115 
116  call deriv5(xs,ys,fun,nsh,5,dpm)
117 
118  det = dpm(3)*dpm(5) - dpm(4)**2
119 
120  rma = xs(1) + ( dpm(2)*dpm(4) - dpm(1)*dpm(5) )/det
121  zma = ys(1) + ( dpm(1)*dpm(4) - dpm(2)*dpm(3) )/det
122 
123  erroma=dsqrt((rma-rm0)**2+(zma-zm0)**2)
124  erroma=erroma/(dabs(r(2,2)-rm0))
125 
126  psima=fun(1)+ dpm(1)*(rma-rm0) + dpm(2)*(zma-zm0)
127  + + 0.5d0*dpm(3)*(rma-rm0)*(rma-rm0)
128  + + dpm(4)*(rma-rm0)*(zma-zm0)
129  + + 0.5d0*dpm(5)*(zma-zm0)*(zma-zm0)
130 
131  return
132  end
133 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
134 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
135 
136  subroutine avr2_c(arr2,nro,nteta,arr1)
137  implicit real*8(a-h,o-z)
138  include 'dim.inc'
139  include 'compol.inc'
140  real*8 arr2(nro,nteta)
141  real*8 arr1(nro)
142  real(8), allocatable :: arrw(:)
143  allocate( arrw(nteta) )
144 
145  do i=1,iplas-1
146  do j=2,nt1
147  arrw(j)=arr2(i,j)
148  enddo
149  arr1(i)=avr1_c(arrw,i)
150  enddo
151 
152  !deallocate( arrw )
153  return
154  end
155 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
156 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
157  subroutine avr2_k(arr2,nro,nteta,arr1)
158  implicit real*8(a-h,o-z)
159  include 'dim.inc'
160  include 'compol.inc'
161  real*8 arr2(nro,nteta)
162  real*8 arr1(nro)
163  real(8), allocatable :: arrw(:)
164  allocate( arrw(nteta) )
165  do i=2,iplas
166  do j=2,nt1
167  arrw(j)=arr2(i,j)
168  enddo
169  arr1(i)=avr1_k(arrw,i)
170  enddo
171 
172  !deallocate( arrw )
173  return
174  end
175 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
176 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
177  real*8 function avr1_c(arr1_c,i)
178  implicit real*8(a-h,o-z)
179  include 'dim.inc'
180  include 'compol.inc'
181  real*8 arr1_c(*)
182 
183  avrg=0.d0
184  sum_vol=0.d0
185  do j=2,nt1
186  avrg=avrg+arr1_c(j)*vol(i,j)
187  sum_vol=sum_vol+vol(i,j)
188  !avrg=avrg+0.5d0*(arr(j+1)+arr(j))*(teta(j+1)-teta(j))
189  enddo
190  avr1_c=avrg/sum_vol
191 
192  return
193  end
194 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
195  real*8 function avr1_k(arr1_k,i)
196  implicit real*8(a-h,o-z)
197  include 'dim.inc'
198  include 'compol.inc'
199  real*8 arr1_k(*)
200 
201  avrg=0.d0
202  sum_vol=0.d0
203  do j=2,nt1
204 
205  if(i.eq.1) then
206  volk=vol1(i,j)+vol4(i,j-1)
207  elseif(i.eq.iplas) then
208  volk=vol2(i-1,j)+vol3(i-1,j-1)
209  else
210  volk=vol1(i,j)+vol2(i-1,j)+vol3(i-1,j-1)+vol4(i,j-1)
211  endif
212  avrg=avrg+arr1_k(j)*volk
213  sum_vol=sum_vol+volk
214  enddo
215  avr1_k=avrg/sum_vol
216 ! write(*,*) 'i=',i
217  return
218  end
219 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
220 
221  real*8 function avr_bnd(arr)
222  implicit real*8(a-h,o-z)
223  include 'parevo.inc'
224  parameter(nkp=njlim)
225  include 'dim.inc'
226  include 'compol.inc'
227  real*8 arr(*)
228  avrg=0.d0
229  sum_len=0.d0
230  do j=2,nt1
231  avrg=avrg+0.5d0*(arr(j+1)+arr(j))*dlt(iplas,j)
232  sum_len=sum_len+dlt(iplas,j)
233  !avrg=avrg+0.5d0*(arr(j+1)+arr(j))*(teta(j+1)-teta(j))
234  enddo
235  avr_bnd=avrg/sum_len
236  !avr_bnd=0.5d0*avrg/pi
237  return
238  end
239 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
240 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
241 
242  subroutine bongri
243 cc int(g(r,r')dr')
244  include 'double.inc'
245  include 'dim.inc'
246  include 'compol.inc'
247 
248  common/com_bgr/ bin_adg(ntp,ntp),dg_dn(ntp)
249 
250 
251  do j=1,nt1
252 
253  rr=r(iplas,j)
254  zz=z(iplas,j)
255 
256  do jb=1,nt1
257 
258  r0=r(iplas,jb)
259  z0=z(iplas,jb)
260 
261  r1=r(iplas,jb+1)
262  z1=z(iplas,jb+1)
263 
264  call bint(rr,zz, r0,z0,r1,z1, fint,1)
265 
266  bin_adg(jb,j)=fint
267 
268  enddo
269  enddo
270 
271  return
272  end
273 
274 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
275  subroutine cof_bon(cps_bon,bps_bon,dps_bon)
276 
277  include 'double.inc'
278  include 'dim.inc'
279  include 'compol.inc'
280 
281  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
282  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
283 
284  common/com_bgr/ bin_adg(ntp,ntp),dg_dn(ntp)
285  dimension pspl_b(ntp),aj(ntp),bj(ntp),dj(ntp)
286  dimension aj_g(ntp),bj_g(ntp),dj_g(ntp)
287 
288  i=iplas !!!!!!!!!!!!!!!!!!!!!!
289 
290  do 10 j=2,nt1
291 
292  a1=a13(i-1,j-1)
293  a2=a34(i-1,j-1)+a12(i-1,j)
294  a3=a24(i-1,j)
295 
296  g1=psi(i-1,j-1)-psi(i,j-1)
297  g2=psi(i-1,j)-psi(i,j)
298  g3=psi(i-1,j+1)-psi(i,j+1)
299 
300  sqk=sq2(i-1,j)+sq3(i-1,j-1)
301  dltk=(dlt(i,j-1)+dlt(i,j))*0.5d0
302  dgdnl=-(a1*g1+a2*g2+a3*g3)+ cur(i,j)*sqk !
303  dg_dn(j)=dgdnl/dltk
304  aj(j)=-(a1+a2+a3)/dltk
305  r0=(3.d0*r(i,j)+r(i-1,j))*0.25d0
306  bj(j)=sqk/dltk/r0
307  dj(j)=sqk*r0/dltk
308 
309  10 continue
310 
311  dg_dn(1)=dg_dn(nt1)
312  dg_dn(nt)=dg_dn(2)
313 
314  aj(1)=aj(nt1)
315  aj(nt)=aj(2)
316  bj(1)=bj(nt1)
317  bj(nt)=bj(2)
318  dj(1)=dj(nt1)
319  dj(nt)=dj(2)
320 
321  do l=2,nt1
322 
323  psb=0.d0
324  ajg=0.d0
325  bjg=0.d0
326  djg=0.d0
327 
328  do jb=2,nt1
329 
330  psb=psb+bin_adg(jb,l)*(dg_dn(jb)+dg_dn(jb+1))*0.5d0
331  ajg=ajg+bin_adg(jb,l)*(aj(jb)+aj(jb+1))*0.5d0
332  bjg=bjg+bin_adg(jb,l)*(bj(jb)+bj(jb+1))*0.5d0
333  djg=djg+bin_adg(jb,l)*(dj(jb)+dj(jb+1))*0.5d0
334 
335  enddo
336 
337  aj_g(l)=ajg
338  bj_g(l)=bjg
339  dj_g(l)=djg
340 
341  enddo
342 
343  aj_g(1)=aj_g(nt1)
344  aj_g(nt)=aj_g(2)
345  bj_g(1)=bj_g(nt1)
346  bj_g(nt)=bj_g(2)
347  dj_g(1)=dj_g(nt1)
348  dj_g(nt)=dj_g(2)
349 
350  sum_len=0.d0
351  do j=2,nt1
352  sum_len=sum_len+dlt(iplas,j)
353  enddo
354 
355  avrg=0.d0
356  do j=2,nt1
357  !avrg=avrg+0.5d0*(aj_G(j+1)+aj_G(j))*(teta(j+1)-teta(j))
358  avrg=avrg+0.5d0*(aj_g(j+1)+aj_g(j))*dlt(iplas,j)
359  enddo
360  !cps_bon=0.5d0*avrg/pi
361  cps_bon=avrg/sum_len
362 
363  avrg=0.d0
364  do j=2,nt1
365  !avrg=avrg+0.5d0*(bj_G(j+1)+bj_G(j))*(teta(j+1)-teta(j))
366  avrg=avrg+0.5d0*(bj_g(j+1)+bj_g(j))*dlt(iplas,j)
367  enddo
368  !bps_bon=0.5d0*avrg/pi
369  bps_bon=avrg/sum_len
370 
371  avrg=0.d0
372  do j=2,nt1
373  !avrg=avrg+0.5d0*(dj_G(j+1)+dj_G(j))*(teta(j+1)-teta(j))
374  avrg=avrg+0.5d0*(dj_g(j+1)+dj_g(j))*dlt(iplas,j)
375  enddo
376  !dps_bon=0.5d0*avrg/pi
377  dps_bon=avrg/sum_len
378 
379  return
380  end
381 
382 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
383 
384  subroutine get_psiext(psi_bon_ext)
385 
386  include 'double.inc'
387  include 'dim.inc'
388  include 'compol.inc'
389 
390  psi_bon_ext=psi_eav
391 
392  return
393  end
394  subroutine get_flfi(flfi_m)
395 
396  include 'double.inc'
397  include 'dim.inc'
398  include 'compol.inc'
399 
400  flfi_m=flx_fi(iplas)
401  !flfi_m=flucfm
402 
403  s_flux_fi=0.d0
404 
405  do i=1,iplas-1
406  s_flux_fi=s_flux_fi-q(i)*(psia(i+1)-psia(i))*psim
407  enddo
408  flfi_m=s_flux_fi
409 
410  return
411  end
412 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
413  subroutine metcof(alp22,alp33,delsc,delv,ac0,skcen,cur_I,
414  * alp33k,delsk,delvk)
415 
416  include 'double.inc'
417  include 'dim.inc'
418  include 'compol.inc'
419 
420  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
421  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
422 
423  common/compsf/ psf(nrp),sqtor(nrp)
424 
425  dimension alp22(nrp),alp33(nrp),delsc(nrp),delv(nrp)
426  dimension alp33k(nrp),delvk(nrp),delsk(nrp)
427  dimension cur_i(nrp),curj_i(nrp),deviat(nrp)
428 
429  sqrt(xx)=dsqrt(xx)
430 
431  do 30 i=1,iplas1
432 
433  zdelsc=0.d0
434  zdelv=0.d0
435  zavrc=0.d0
436 
437  do 35 j=2,nt1
438 
439  r1=r(i,j)
440  r2=r(i+1,j)
441  r3=r(i+1,j+1)
442  r4=r(i,j+1)
443  r0=(r1+r2+r3+r4)*0.25d0
444 
445  zdelv=zdelv+vol(i,j)
446  zdelsc=zdelsc+s(i,j)
447  zavrc=zavrc+s(i,j)/r0
448 
449  35 continue
450 
451  alp33(i)=zdelsc/zavrc
452  delsc(i)=zdelsc
453  delv(i)=zdelv
454 
455  30 continue
456 
457  sqtor(1)=0.d0
458 
459  do 33 i=2,iplas
460  sqtor(i)=sqtor(i-1)+delsc(i-1)
461  33 continue
462 
463  do 40 i=2,iplas
464 
465  zdelv=0.d0
466  zdelsk=0.d0
467  zavrk=0.d0
468 
469  do 45 j=2,nt1
470 
471  if(i.ne.iplas) then
472  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
473  r0=r(i,j)
474  else
475  sqk=sq2(i-1,j)+sq3(i-1,j-1)
476  r0=(3.d0*r(i,j)+r(i-1,j))*0.25d0
477  endif
478 
479  zdelsk=zdelsk+sqk
480  zdelv=zdelv +sqk*r0
481  zavrk=zavrk+sqk/r0
482 
483  45 continue
484 
485  alp33k(i)=zdelsk/zavrk
486  delsk(i)=zdelsk
487  delvk(i)=zdelv
488 
489  40 continue
490 
491 
492  do 10 i=2,iplas1
493 
494  samn=0.d0
495  sa0 =0.d0
496  sapl=0.d0
497 
498  do 110 j=2,nt1
499 
500  a1=a13(i-1,j-1)
501  a2=a34(i-1,j-1)+a12(i-1,j)
502  a3=a24(i-1,j)
503  a4=a23(i-1,j-1)+a14(i,j-1)
504  a6=a14(i,j)+a23(i-1,j)
505  a7=a24(i,j-1)
506  a8=a34(i,j-1)+a12(i,j)
507  a9=a13(i,j)
508  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
509 
510  zamn=a1+a2+a3
511  za0 =a4+a5+a6
512  zapl=a7+a8+a9
513 
514  samn=samn+zamn
515  sa0 =sa0 +za0
516  sapl=sapl+zapl
517 
518  110 continue
519 
520  alp22(i-1)=samn*delsc(i-1)
521  if(i.eq.iplas1) alp22(i)=sapl*delsc(i)
522  10 continue
523 
524  do i=1,iplas1
525  dpsids=psipla*(psia(i+1)-psia(i))/delsc(i)
526  cur_i(i)=dpsids*alp22(i)
527  enddo
528 
529 !!!!!!!!!!///////////////////////////////////////////
530 
531 !!! equation in central node
532 
533  ac0=0.d0
534  skcen=0.d0
535 
536  do 20 j=2,nt1
537 
538  acj=a12(1,j)+a24(1,j)+a34(1,j-1)+a13(1,j-1)
539  ac0=ac0-acj
540 
541  skcen=skcen+sq1(1,j)+sq4(1,j-1)
542 
543  20 continue
544 
545 ! dfdpsi(1)=
546 ! * ( (ac0*(psiax-psf(2)))/skcen-rm*dpdpsi(1) )*rm
547 
548 !!!!!!!!!!!!!!!! deviation !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
549 
550  curj_i(1)=cur(1,2)*skcen
551  deviat(1)=cur_i(1)-curj_i(1)
552 
553  do i=2,iplas1
554  dcurj_i=0.d0
555  do j=2,nt1
556  if(i.ne.iplas) then
557  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
558  else
559  sqk=sq2(i-1,j)+sq3(i-1,j-1)
560  endif
561  dcurj_i=dcurj_i+cur(i,j)*sqk
562  enddo
563  curj_i(i)=curj_i(i-1)+dcurj_i
564  dcur_i=cur_i(i)-cur_i(i-1)
565  deviat(i)=(dcur_i-dcurj_i)/delsk(i)
566  enddo
567 
568  diff=(alp22(iplas-1)-alp22(iplas-2))/delsk(iplas-1)
569  alp22(iplas)=alp22(iplas-1)+0.5d0*diff*delsk(iplas)
570 
571  return
572  end
573 
574 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
575 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
576  subroutine qst_b
577 
578  include 'double.inc'
579  include 'dim.inc'
580  parameter(nshp=ntp+1)
581  include 'compol.inc'
582  !common /combsh/ rm0,zm0,rc0,zc0,asp0,el_up,el_lw,tr_up,tr_lw,nbsh
583  common/efites/ fcefit,rc_efit,iefit
584  dimension xs(nshp),ys(nshp),fun(nshp),dp(5)
585 
586  sqrt(xx)=dsqrt(xx)
587 
588  nsh=1
589 
590  xs(nsh)=r(1,1)
591  ys(nsh)=z(1,1)
592  fun(nsh)=psin(1,1)
593 
594  do 100 j=2,nt1
595 
596  nsh=nsh+1
597  xs(nsh)=r(2,j)
598  ys(nsh)=z(2,j)
599  fun(nsh)=psin(2,j)
600 
601  100 continue
602 
603  call deriv5(xs,ys,fun,nsh,5,dp)
604 
605  drr=dp(3)
606  drz=dp(4)
607  dzz=dp(5)
608 
609  tg2a=2.d0*drz/(drr-dzz)
610  cos2a=1.d0/sqrt(1.d0+tg2a**2)
611  sin2a=cos2a*tg2a
612 
613  dxx=0.5d0*(drr+dzz)+0.5d0*cos2a*(drr-dzz)+sin2a*drz
614  dyy=0.5d0*(drr+dzz)+0.5d0*cos2a*(dzz-drr)-sin2a*drz
615 c write(6,*) 'drr,dzz,drz',drr,dzz,drz
616 c write(6,*) 'dxx,dyy',dxx,dyy
617  ! bfcen=qcen*dsqrt(dxx*dyy)*(psim-psip)
618  ! fcen=bfcen*rm
619  ! fcen=fcefit
620 c write(6,*) 'qst:b0,r0',b0ax,rc0
621  if(iefit.eq.1) r0ax=rc_efit
622 
623  fvac=b0ax*r0ax
624  !!fvac=fcefit
625 
626 c write(6,*) 'qst:fvac',fvac
627 
628  dpsi=(psia(iplas)-psia(iplas-1))
629  ps14=-0.25d0*dpsi
630  !ffp=tabf(ps14)
631  ffp=dfdpsi(iplas)
632  f(iplas-1)=sqrt(fvac**2-ffp*(psim-psip)*dpsi)
633  !!!!!!!!!!f(iplas-1)=fvac !!!!!!!!!!!
634 
635  do 200 i=iplas-2,1,-1
636  dpsi=0.5d0*(psia(i+2)-psia(i))*(psim-psip)
637  ffp=dfdpsi(i+1)
638  f(i)=sqrt(f(i+1)**2-2.d0*ffp*dpsi)
639  200 continue
640 
641  f(iplas)=fvac
642 
643  !open(1,file='q.pr')
644 
645  flucfm=0.d0
646  flx_fi(1)=0.d0
647 
648  do 300 i=1,iplas
649 
650  dflufi=0.d0
651 
652  do 310 j=2,nt1
653 
654  if(i.ne.iplas) then
655 
656  r1=r(i,j)
657  r2=r(i+1,j)
658  r3=r(i+1,j+1)
659  r4=r(i,j+1)
660  r0=(r1+r2+r3+r4)*0.25d0
661  dflufi=dflufi+s(i,j)*f(i)/r0
662 
663  else
664 
665  r0=r(i,j)
666  dflufi=dflufi+(sq3(i-1,j-1)+sq2(i-1,j))*f(i)/r0
667 
668  endif
669 
670  310 continue
671 
672  if(i.ne.iplas) then
673 
674  q2pi=-dflufi/((psia(i+1)-psia(i))*(psim-psip))
675  flucfm=flucfm+dflufi
676  flx_fi(i+1)=flx_fi(i)+dflufi
677 
678  else
679 
680  q2pi=-2.d0*dflufi/((psia(i)-psia(i-1))*(psim-psip))
681 
682  endif
683 
684  !q(i)=0.5d0*q2pi/pi
685  q(i)=q2pi
686 
687  !write(1,*) 'q(i),i',q(i),f(i),i
688 
689  300 continue
690  q(iplas)=1.5d0*q(iplas-1)-0.5d0*q(iplas-2)
691 
692  ! close(1)
693 
694  return
695  end
696 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
697 
698 
699  SUBROUTINE extrpc(X,F, X0,X1,X2,X3, F1,F2,F3)
700 
701  IMPLICIT REAL*8(a-h,o-z)
702 
703 
704  x1m = 0.5*(x0+x1)
705  x2m = 0.5*(x1+x2)
706  x3m = 0.5*(x2+x3)
707 
708  d1f2m =(f2-f1)/(x2m-x1m)
709  d1f3m =(f3-f2)/(x3m-x2m)
710 
711  d2f23 = (d1f3m-d1f2m)/(x3m-x1m)
712 c d2f23 = 0.
713  f = f1+(x-x1m)*(d1f2m + (x-x2m)*d2f23)
714 
715  RETURN
716  END
717 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
718 
719  SUBROUTINE extrp2(X,F, X1,X2,X3, F1,F2,F3)
720 
721  IMPLICIT REAL*8(a-h,o-z)
722 
723  dfdx=(f3-f1)/(x3-x1)
724 
725  d2fdx=( (f3-f2)/(x3-x2) - (f2-f1)/(x2-x1) )/(x3-x1)
726 
727  f= f2+(x-x2)*( dfdx + d2fdx*(x-x2) )
728 
729  RETURN
730  END
731 
732 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
733  real*8 function qvadin(x,
734  * xm,x0,xp,
735  * fm,f0,fp )
736  include 'double.inc'
737 
738  f=f0+( fp-fm + ( (fp-f0)/(xp-x0)-(f0-fm)/(x0-xm) )*(x-x0) )
739  * *(x-x0)/(xp-xm)
740 
741  qvadin=f
742 
743  return
744  end
745 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
746 
747  subroutine bt_pol(betpol)
748 
749  include 'double.inc'
750  include 'dim.inc'
751  include 'compol.inc'
752  !!amu0=0.4d0*pi
753  sqcen=0.d0
754 
755  do j=2,nt1
756  sqcen = sqcen + sq1(1,j) + sq4(1,j)
757  enddo
758  psn=psin(1,2)
759 
760  zpres = funppp(psn)
761  pintg = zpres*sqcen
762 
763  do i=2,iplas1
764  do j=2,nt1
765 
766  psn=psin(i,j)
767  zpres=funppp(psn)
768  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
769  pintg=pintg+zpres*sqk
770 
771  enddo
772  enddo
773 
774  betpol=8.d0*pi*pintg/(tokp*tokp)
775  betpol=betpol*(psim-psip)*cnor/amu0**2
776 
777  return
778  end
779 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
780 
781  subroutine skbetp(betplx,betpol)
782 
783  include 'double.inc'
784  include 'dim.inc'
785  include 'compol.inc'
786 
787  call bt_pol(betpol)
788 
789  if(ngav.eq.0) then
790 
791  zcoin = betplx / betpol
792  !write(6,*) 'betplx/betpol',zcoin
793  zcoin = (1.d0/zcoin-1.d0)*tok/(cnor*tokff)+1.d0
794  coin = 1.d0/zcoin
795  coin = (coin-1.d0)*0.33d0+1.d0
796 
797  else
798 
799  coin=betplx/betpol
800  do i=1,iplas
801  dpdpsi(i)=coin*dpdpsi(i)
802  enddo
803  endif
804 
805  call taburs(1,coin,nurs)
806 
807  return
808  end
809 
810 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
811 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
812 
813  subroutine cur_avg
814 
815  include 'double.inc'
816  include 'dim.inc'
817  include 'compol.inc'
818  common/savt0/ psi0(nrp),fi0(nrp),f0(nrp),ri0(nrp),q0(nrp),
819  * dpsidt(nrp),dfidt(nrp),rm0,ac0n,skcen0
820  common/com_flag/kastr
821  common /com_jb/ bj_av(nrp),curfi_av(nrp)
822  common /com_b2/ b2_av(nrp)
823  common /com_eb/ eb(nrp),eb_c(nrp),epar_c(nrp)
824  common/com_heat_dj/ wdj(nrp)
825  dimension bfjf(nrp)
826  dimension alfa22(nrp),alfa33(nrp),ds(nrp),dv(nrp),cj(nrp),
827  * alp33k(nrp),dsk(nrp),dvk(nrp)
828  dimension sigma(nrp)
829  dimension z_nvzk(nrp)
830 
831 
832  bfjf(1)=cur(1,2)*f(1)/rm
833 
834  do i=2,iplas1
835  volk_i=0.d0
836  bj_i=0.d0
837  do j=2,nt1
838 
839  volk=vol1(i,j)+vol2(i-1,j)+vol3(i-1,j-1)+vol4(i,j-1)
840  volk_i=volk_i+volk
841  bj_i=bj_i+cur(i,j)*(f(i)+f(i-1))*0.5d0*volk/r(i,j)
842 
843  enddo
844  bfjf(i)=bj_i/volk_i
845  enddo
846 
847  call metcof(alfa22,alfa33,ds,dv,ac0,skcen,cj,
848  * alp33k,dsk,dvk)
849 
850  curfi_av(1)=cur(1,2)
851  bj_av(1)=bfjf(1)
852  b2_av(1)=(f(1)/rm)**2
853  cj(iplas)=tokp*amu0
854  do i=2,iplas
855  bj_av(i)=(cj(i)*f(i-1)-cj(i-1)*f(i))/dvk(i)
856  curfi_av(i)=(cj(i)-cj(i-1))/dsk(i)
857  enddo
858  do i=1,iplas1
859  b2_av(i)=( f(i)*(flx_fi(i+1)-flx_fi(i))-
860  & cj(i)*psim*(psia(i+1)-psia(i)) )/dv(i)
861  enddo
862 
863  b2_av(1)=(f(1)/rm)**2 !test
864 
865  i=iplas
866 
867 ! B2_av(i)=( F(i)*(flx_fi(i)-flx_fi(i-1))-
868 ! & cJ(i)*psim*(psia(i)-psia(i-1)) )/dvk(i)/2.d0
869 
870  b2_av(i)=(3.d0*b2_av(i-1) - b2_av(i-2))/2.0d0
871 
872  b2_av_m=(3.d0*b2_av(1) - b2_av(2))/2.0d0
873 
874 
875 
876 
877 
878  do i=1,iplas
879  xa=dfloat(i-1)/dfloat(iplas-1)
880  if(kastr.eq.1) then
881  sigma(i) =dabs( fun_sig(xa,i)*dtim /amu0)
882  else
883  sigma(i) =dabs( fun_sig_a(xa,i)*dtim/amu0 )
884  endif
885  enddo
886 
887  do i=2,iplas-1
888  eb(i)= dpsidt(i)*(flx_fi(i+1)-flx_fi(i-1))
889  & -dfidt(i)*(psia(i+1)-psia(i-1))*psim
890  eb(i)=eb(i)/(dv(i)+dv(i-1))/dtim
891  sebeb=sigma(i)*eb(i)**2
892  wdj(i)=sebeb/b2_av(i)
893  enddo
894 
895  i=1
896 
897  fm=(3.d0*f(1)-f(2))/2.0d0
898  ebm=dpsidt(1)*fm/rm**2/dtim
899  eb(1)= ebm
900  em=dpsidt(1)/rm/dtim
901  !SEBEB=sigma(1)*EB**2
902  see=sigma(1)*em**2
903  !Wdj(1)=SEBEB/B2_av_m
904  wdj(1)=see
905 
906  i=iplas
907  eb(i)= dpsidt(i)*(flx_fi(i)-flx_fi(i-1))
908  & -dfidt(i)*(psia(i)-psia(i-1))*psim
909  eb(i)=eb(i)/(dv(i)+dv(i-1))/dtim
910  sebeb=sigma(i)*eb(i)**2
911  wdj(i)=sebeb/b2_av(i)
912 
913 !!!!!!!!! (E,B) in cells
914 
915  do i=1,iplas-1
916  dpsdt05=0.5d0*(dpsidt(i)+dpsidt(i+1))/dtim
917  dfidt05=0.5d0*(dfidt(i)+dfidt(i+1))/dtim
918 
919  eb_c(i)= dpsdt05*(flx_fi(i+1)-flx_fi(i))
920  & -dfidt05*(psia(i+1)-psia(i))*psim
921  !EB_c(i)=EB_c(i)/dv(i)
922  eb_c(i)=-0.5d0*(eb(i)+eb(i+1))
923  epar_c(i)=eb_c(i)/dsqrt(b2_av(i))
924  enddo
925 
926  do i=1,iplas
927  z_nvzk(i)=amu0*sigma(i)*eb(i)+bj_av(i)
928  enddo
929 
930 
931  return
932  end
933 
934 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
935 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
936 
937  subroutine cur_avg_
938 
939  include 'double.inc'
940  include 'dim.inc'
941  include 'compol.inc'
942  common/savt0/ psi0(nrp),fi0(nrp),f0(nrp),ri0(nrp),q0(nrp),
943  * dpsidt(nrp),dfidt(nrp),rm0,ac0n,skcen0
944  common/com_flag/kastr
945  common /com_jb/ bj_av(nrp),curfi_av(nrp)
946  common /com_b2/ b2_av(nrp)
947  common/com_heat_dj/ wdj(nrp)
948  dimension bfjf(nrp)
949  dimension alfa22(nrp),alfa33(nrp),ds(nrp),dv(nrp),cj(nrp),
950  * alp33k(nrp),dsk(nrp),dvk(nrp)
951  dimension sigma(nrp)
952 
953 
954  bfjf(1)=cur(1,2)*f(1)/rm
955 
956  do i=2,iplas1
957  volk_i=0.d0
958  bj_i=0.d0
959  do j=2,nt1
960 
961  volk=vol1(i,j)+vol2(i-1,j)+vol3(i-1,j-1)+vol4(i,j-1)
962  volk_i=volk_i+volk
963  bj_i=bj_i+cur(i,j)*(f(i)+f(i-1))*0.5d0*volk/r(i,j)
964 
965  enddo
966  bfjf(i)=bj_i/volk_i
967  enddo
968 
969  call metcof(alfa22,alfa33,ds,dv,ac0,skcen,cj,
970  * alp33k,dsk,dvk)
971 
972  curfi_av(1)=cur(1,2)
973  bj_av(1)=bfjf(1)
974  b2_av(1)=(f(1)/rm)**2
975  cj(iplas)=tokp*amu0
976  do i=2,iplas
977  bj_av(i)=(cj(i)*f(i-1)-cj(i-1)*f(i))/dvk(i)
978  curfi_av(i)=(cj(i)-cj(i-1))/dsk(i)
979  enddo
980  do i=1,iplas1
981  b2_av(i)=( f(i)*(flx_fi(i+1)-flx_fi(i))-
982  & cj(i)*psim*(psia(i+1)-psia(i)) )/dv(i)
983  enddo
984 
985  b2_av(1)=(f(1)/rm)**2 !test
986 
987  i=iplas
988 
989 ! B2_av(i)=( F(i)*(flx_fi(i)-flx_fi(i-1))-
990 ! & cJ(i)*psim*(psia(i)-psia(i-1)) )/dvk(i)/2.d0
991 
992  b2_av(i)=(3.d0*b2_av(i-1) - b2_av(i-2))/2.0d0
993 
994  b2_av_m=(3.d0*b2_av(1) - b2_av(2))/2.0d0
995 
996 
997 
998 
999 
1000  do i=1,iplas
1001  xa=dfloat(i-1)/dfloat(iplas-1)
1002  if(kastr.eq.1) then
1003  sigma(i) =dabs( fun_sig(xa,i)*dtim /amu0)
1004  else
1005  sigma(i) =dabs( fun_sig_a(xa,i)*dtim/amu0 )
1006  endif
1007  enddo
1008 
1009  do i=2,iplas-1
1010  eb= dpsidt(i)*(flx_fi(i+1)-flx_fi(i-1))
1011  & -dfidt(i)*(psia(i+1)-psia(i-1))*psim
1012  eb=eb/(dv(i)+dv(i-1))/dtim
1013  sebeb=sigma(i)*eb**2
1014  wdj(i)=sebeb/b2_av(i)
1015  enddo
1016 
1017  i=1
1018 
1019  fm=(3.d0*f(1)-f(2))/2.0d0
1020  !EBm=dpsidt(1)*Fm/rm**2/dtim
1021  em=dpsidt(1)/rm/dtim
1022  !SEBEB=sigma(1)*EB**2
1023  see=sigma(1)*em**2
1024  !Wdj(1)=SEBEB/B2_av_m
1025  wdj(1)=see
1026 
1027  i=iplas
1028  eb= dpsidt(i)*(flx_fi(i)-flx_fi(i-1))
1029  & -dfidt(i)*(psia(i)-psia(i-1))*psim
1030  eb=eb/(dv(i)+dv(i-1))/dtim
1031  sebeb=sigma(i)*eb**2
1032  wdj(i)=sebeb/b2_av(i)
1033 
1034  return
1035  end
1036 
1037 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1038 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1039 
1040  function numlin(i,j,nr,nt)
1041  common /comctr/ ngav
1042 
1043 !!! definition of element number using it indexes(i,j) in 2-dim. array
1044  if(i.ne.1) then
1045  if(j.eq.1) then
1046  jj=nt-1
1047  elseif(j.eq.nt) then
1048  jj=2
1049  else
1050  jj=j
1051  endif
1052  numlin=(i-2)*(nt-2)+jj
1053  else
1054  numlin=1
1055  endif
1056  return
1057  end
1058 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1059 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1060 
1061  subroutine matpla
1062 
1063 !!!!!! IL -- number of equation (matrix line)
1064 !!!!!! IM -- element number in one-dimensional array A (a(im))
1065 !!!!!! IA(IL) -- number of first nonzero element in line IL
1066 !!!!!! JA(IM) -- number of matrix column for element a(im)
1067 
1068  include 'double.inc'
1069  include 'dim.inc'
1070  include 'compol.inc'
1071 
1072  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
1073  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
1074 
1075  il=0
1076  im=0
1077 
1078 !!!!!!!!!!! equation for central point
1079 
1080  il=1
1081  im=1
1082 
1083  ia(il)=1
1084  ja(im)=1
1085  a(im)=0.d0
1086 
1087  do 20 j=2,nt1
1088 
1089  im=im+1
1090  ja(im)=numlin(2,j,nr,nt)
1091  a(im)=a12(1,j)+a24(1,j)+a34(1,j-1)+a13(1,j-1)
1092  a(1)=a(1)-a(im)
1093 
1094  20 continue
1095 
1096 !!!!!!!!!!! the main loop
1097 
1098  do 100 i=2,iplas-1
1099 
1100  !!!!!
1101  !
1102  j=2 !
1103  !
1104  !!!!!
1105 
1106  a1=a13(i-1,j-1)
1107  a2=a34(i-1,j-1)+a12(i-1,j)
1108  a3=a24(i-1,j)
1109  a4=a23(i-1,j-1)+a14(i,j-1)
1110  a6=a14(i,j)+a23(i-1,j)
1111  a7=a24(i,j-1)
1112  a8=a34(i,j-1)+a12(i,j)
1113  a9=a13(i,j)
1114 
1115  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
1116 
1117  il=il+1
1118  ia(il)=im+1
1119 
1120  if(i.ne.2) then
1121 
1122 !2!cof. to (i-1,j)
1123 
1124  im=im+1
1125  a(im)=a2
1126  ja(im)=numlin(i-1,j,nr,nt)
1127 
1128 !3!cof. to (i-1,j+1)
1129 
1130  im=im+1
1131  a(im)=a3
1132  ja(im)=numlin(i-1,j+1,nr,nt)
1133 
1134 !1!cof. to (i-1,j-1)
1135 
1136  im=im+1
1137  a(im)=a1
1138  ja(im)=numlin(i-1,j-1,nr,nt)
1139 
1140 !5!cof. to (i,j)
1141 
1142  im=im+1
1143  a(im)=a5
1144  ja(im)=numlin(i,j,nr,nt)
1145 
1146 !6!cof. to (i,j+1)
1147 
1148  im=im+1
1149  a(im)=a6
1150  ja(im)=numlin(i,j+1,nr,nt)
1151 
1152 !4!cof. to (i,j-1)
1153 
1154  im=im+1
1155  a(im)=a4
1156  ja(im)=numlin(i,j-1,nr,nt)
1157 
1158  if(i.eq.iplas-1) go to 71
1159 
1160 !8!cof. to (i+1,j)
1161 
1162  im=im+1
1163  a(im)=a8
1164  ja(im)=numlin(i+1,j,nr,nt)
1165 
1166 !9!cof. to (i+1,j+1)
1167 
1168  im=im+1
1169  a(im)=a9
1170  ja(im)=numlin(i+1,j+1,nr,nt)
1171 
1172 !7!cof. to (i+1,j-1)
1173 
1174  im=im+1
1175  a(im)=a7
1176  ja(im)=numlin(i+1,j-1,nr,nt)
1177 
1178  71 continue
1179 
1180  else
1181 
1182 !1!cof. to central point
1183 
1184  im=im+1
1185  a(im)=a1+a2+a3
1186  ja(im)=numlin(1,j,nr,nt)
1187 
1188 !5!cof. to (2,j)
1189 
1190  im=im+1
1191  a(im)=a5
1192  ja(im)=numlin(2,j,nr,nt)
1193 
1194 !6!cof. to (2,j+1)
1195 
1196  im=im+1
1197  a(im)=a6
1198  ja(im)=numlin(2,j+1,nr,nt)
1199 
1200 !4!cof. to (2,j-1)
1201 
1202  im=im+1
1203  a(im)=a4
1204  ja(im)=numlin(2,j-1,nr,nt)
1205 
1206 !8!cof. to (3,j)
1207 
1208  im=im+1
1209  a(im)=a8
1210  ja(im)=numlin(3,j,nr,nt)
1211 
1212 !9!cof. to (3,j+1)
1213 
1214  im=im+1
1215  a(im)=a9
1216  ja(im)=numlin(3,j+1,nr,nt)
1217 
1218 !7!cof. to (3,j-1)
1219 
1220  im=im+1
1221  a(im)=a7
1222  ja(im)=numlin(3,j-1,nr,nt)
1223 
1224  endif
1225 
1226  do 110 j=3,nt2
1227 
1228  a1=a13(i-1,j-1)
1229  a2=a34(i-1,j-1)+a12(i-1,j)
1230  a3=a24(i-1,j)
1231  a4=a23(i-1,j-1)+a14(i,j-1)
1232  a6=a14(i,j)+a23(i-1,j)
1233  a7=a24(i,j-1)
1234  a8=a34(i,j-1)+a12(i,j)
1235  a9=a13(i,j)
1236 
1237  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
1238 
1239  il=il+1
1240  ia(il)=im+1
1241 
1242  if(i.ne.2) then
1243 
1244 !1!cof. to (i-1,j-1)
1245 
1246  im=im+1
1247  a(im)=a1
1248  ja(im)=numlin(i-1,j-1,nr,nt)
1249 
1250 !2!cof. to (i-1,j)
1251 
1252  im=im+1
1253  a(im)=a2
1254  ja(im)=numlin(i-1,j,nr,nt)
1255 
1256 !3!cof. to (i-1,j+1)
1257 
1258  im=im+1
1259  a(im)=a3
1260  ja(im)=numlin(i-1,j+1,nr,nt)
1261 
1262 !3!cof. to (i,j-1)
1263 
1264  im=im+1
1265  a(im)=a4
1266  ja(im)=numlin(i,j-1,nr,nt)
1267 
1268 !5!cof. to (i,j)
1269 
1270  im=im+1
1271  a(im)=a5
1272  ja(im)=numlin(i,j,nr,nt)
1273 
1274 !6!cof. to (i,j+1)
1275 
1276  im=im+1
1277  a(im)=a6
1278  ja(im)=numlin(i,j+1,nr,nt)
1279 
1280  if(i.eq.iplas-1) go to 7
1281 
1282 !7!cof. to (i+1,j-1)
1283 
1284  im=im+1
1285  a(im)=a7
1286  ja(im)=numlin(i+1,j-1,nr,nt)
1287 
1288 !8!cof. to (i+1,j)
1289 
1290  im=im+1
1291  a(im)=a8
1292  ja(im)=numlin(i+1,j,nr,nt)
1293 
1294 !9!cof. to (i+1,j+1)
1295 
1296  im=im+1
1297  a(im)=a9
1298  ja(im)=numlin(i+1,j+1,nr,nt)
1299 
1300  7 continue
1301 
1302  else
1303 
1304 !1!cof. to central point
1305 
1306  im=im+1
1307  a(im)=a1+a2+a3
1308  ja(im)=numlin(1,j,nr,nt)
1309 
1310 !2!cof. to (2,j-1)
1311 
1312  im=im+1
1313  a(im)=a4
1314  ja(im)=numlin(2,j-1,nr,nt)
1315 
1316 !3!cof. to (2,j)
1317 
1318  im=im+1
1319  a(im)=a5
1320  ja(im)=numlin(2,j,nr,nt)
1321 
1322 !4!cof. to (2,j+1)
1323 
1324  im=im+1
1325  a(im)=a6
1326  ja(im)=numlin(2,j+1,nr,nt)
1327 
1328 !5!cof. to (3,j-1)
1329 
1330  im=im+1
1331  a(im)=a7
1332  ja(im)=numlin(3,j-1,nr,nt)
1333 
1334 !6!cof. to (3,j)
1335 
1336  im=im+1
1337  a(im)=a8
1338  ja(im)=numlin(3,j,nr,nt)
1339 
1340 !7!cof. to (3,j+1)
1341 
1342  im=im+1
1343  a(im)=a9
1344  ja(im)=numlin(3,j+1,nr,nt)
1345 
1346  endif
1347 
1348  110 continue
1349 
1350  !!!!!!!
1351  !
1352  j=nt1 !
1353  !
1354  !!!!!!!
1355 
1356  a1=a13(i-1,j-1)
1357  a2=a34(i-1,j-1)+a12(i-1,j)
1358  a3=a24(i-1,j)
1359  a4=a23(i-1,j-1)+a14(i,j-1)
1360  a6=a14(i,j)+a23(i-1,j)
1361  a7=a24(i,j-1)
1362  a8=a34(i,j-1)+a12(i,j)
1363  a9=a13(i,j)
1364 
1365  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
1366 
1367  il=il+1
1368  ia(il)=im+1
1369 
1370  if(i.ne.2) then
1371 
1372 !3!cof. to (i-1,j+1)
1373 
1374  im=im+1
1375  a(im)=a3
1376  ja(im)=numlin(i-1,j+1,nr,nt)
1377 
1378 !1!cof. to (i-1,j-1)
1379 
1380  im=im+1
1381  a(im)=a1
1382  ja(im)=numlin(i-1,j-1,nr,nt)
1383 
1384 !2!cof. to (i-1,j)
1385 
1386  im=im+1
1387  a(im)=a2
1388  ja(im)=numlin(i-1,j,nr,nt)
1389 
1390 !6!cof. to (i,j+1)
1391 
1392  im=im+1
1393  a(im)=a6
1394  ja(im)=numlin(i,j+1,nr,nt)
1395 
1396 !4!cof. to (i,j-1)
1397 
1398  im=im+1
1399  a(im)=a4
1400  ja(im)=numlin(i,j-1,nr,nt)
1401 
1402 !5!cof. to (i,j)
1403 
1404  im=im+1
1405  a(im)=a5
1406  ja(im)=numlin(i,j,nr,nt)
1407 
1408  if(i.eq.iplas-1) go to 72
1409 
1410 !9!cof. to (i+1,j+1)
1411 
1412  im=im+1
1413  a(im)=a9
1414  ja(im)=numlin(i+1,j+1,nr,nt)
1415 
1416 !7!cof. to (i+1,j-1)
1417 
1418  im=im+1
1419  a(im)=a7
1420  ja(im)=numlin(i+1,j-1,nr,nt)
1421 
1422 !8!cof. to (i+1,j)
1423 
1424  im=im+1
1425  a(im)=a8
1426  ja(im)=numlin(i+1,j,nr,nt)
1427 
1428  72 continue
1429 
1430  else
1431 
1432 !1!cof. to central point
1433 
1434  im=im+1
1435  a(im)=a1+a2+a3
1436  ja(im)=numlin(1,j,nr,nt)
1437 
1438 !6!cof. to (2,j+1)
1439 
1440  im=im+1
1441  a(im)=a6
1442  ja(im)=numlin(2,j+1,nr,nt)
1443 
1444 !4!cof. to (2,j-1)
1445 
1446  im=im+1
1447  a(im)=a4
1448  ja(im)=numlin(2,j-1,nr,nt)
1449 
1450 !5!cof. to (2,j)
1451 
1452  im=im+1
1453  a(im)=a5
1454  ja(im)=numlin(2,j,nr,nt)
1455 
1456 !9!cof. to (3,j+1)
1457 
1458  im=im+1
1459  a(im)=a9
1460  ja(im)=numlin(3,j+1,nr,nt)
1461 
1462 !7!cof. to (3,j-1)
1463 
1464  im=im+1
1465  a(im)=a7
1466  ja(im)=numlin(3,j-1,nr,nt)
1467 
1468 !8!cof. to (3,j)
1469 
1470  im=im+1
1471  a(im)=a8
1472  ja(im)=numlin(3,j,nr,nt)
1473 
1474  endif
1475 
1476  100 continue
1477  !write(6,*) 'matrix:il=',il
1478 
1479  neqpla=il
1480 
1481  !write(6,*) 'number of equations:',neq,neqp
1482  !write(6,*) 'im=',im
1483  !write(6,*) 'lp=',lp
1484 
1485  il=il+1
1486  ia(il)=im+1
1487 
1488  if( (itin/nitdel*nitdel+nitbeg) .EQ. itin ) then
1489 
1490  do km=1,im
1491 
1492  app0(km)=a(km)
1493 
1494  enddo
1495 
1496  !write(6,*) 'matpla:aop0:itin,nitdel',itin,nitdel
1497 
1498  elseif(itin.gt.nitbeg) then
1499 
1500  do km=1,im
1501 
1502  dapp(km)=a(km)-app0(km)
1503 
1504  enddo
1505 
1506  !write(6,*) 'matpla:dapp:itin',itin
1507 
1508  endif
1509 
1510  return
1511  end
1512 
1513 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1514 
1515 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1516  subroutine inter_q(ps,nw,q_ps)
1517 
1518  include 'double.inc'
1519  include 'dim.inc'
1520  parameter(nrpl=nrp+1)
1521  parameter(nrpl4=nrpl+4,nrpl6=nrpl4*6)
1522  include 'compol.inc'
1523  dimension ps(*),q_ps(*)
1524  dimension q_w(nrp+1),ps_w(nrp+1)
1525  real*8 rrk(nrpl4),cck(nrpl4),wrk(nrpl6)
1526  real*8 cwk(4)
1527 
1528 !extrapolation q to axis
1529 
1530  q_m=1.5d0*q(1)-0.5d0*q(2)
1531 
1532 !extrapolation q to boundary
1533 
1534  q_bon=1.5d0*q(iplas-1)-0.5d0*q(iplas-2)
1535 
1536 !arrays q_w,ps_w initialization
1537 
1538  ps_w(1)=0.d0
1539  q_w(1)=q_m
1540 
1541  do i=1,iplas-1
1542 
1543  ps_w(i+1)=1.d0-0.5d0*(psia(i)+psia(i+1))
1544  q_w(i+1)=q(i)
1545 
1546  enddo
1547 
1548  ps_w(iplas+1)=1.d0
1549  q_w(iplas+1)=q_bon
1550 
1551  n3spl=iplas+1
1552 
1553  CALL e01baf(n3spl,ps_w,q_w,rrk,cck,
1554  * n3spl+4,wrk,6*n3spl+16,ifail)
1555 
1556  do i=2,nw-1
1557  zspl=ps(i)
1558  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1559  q_ps(i)=cwk(1)
1560  enddo
1561 
1562  q_ps(1)=q_w(1)
1563  q_ps(nw)=q_w(iplas+1)
1564 
1565  return
1566  end
1567 
1568 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
extrpc
subroutine extrpc(X, F, X0, X1, X2, X3, F1, F2, F3)
Definition: com_sub.f:699
e02bcf
subroutine e02bcf(NCAP7, K, C, X, LEFT, S, IFAIL)
Definition: NAG.f:2761
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
avr1_k
real *8 function avr1_k(arr1_k, i)
Definition: com_sub.f:195
helena21_mod::dpdpsi
real(r8) function dpdpsi(psi_n)
Definition: helena21_mod.f90:122
get_flfi
subroutine get_flfi(flfi_m)
Definition: com_sub.f:394
cur_avg_
subroutine cur_avg_
Definition: com_sub.f:937
numlin
function numlin(i, j, nr, nt)
Definition: com_sub.f:1040
grdef
subroutine grdef(igdf)
Definition: com_sub.f:7
bint
subroutine bint(X, Y, R0, Z0, r1, z1, F, I)
Definition: scu.f:557
e01baf
subroutine e01baf(M, X, Y, K, C, LCK, WRK, LWRK, IFAIL)
Definition: NAG.f:2511
get_psiext
subroutine get_psiext(psi_bon_ext)
Definition: com_sub.f:384
extrp2
subroutine extrp2(X, F, X1, X2, X3, F1, F2, F3)
Definition: com_sub.f:719
skbetp
subroutine skbetp(betplx, betpol)
Definition: com_sub.f:781
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
matpla
subroutine matpla
Definition: com_sub.f:1061
cur_avg
subroutine cur_avg
Definition: com_sub.f:813
avr_bnd
real *8 function avr_bnd(arr)
Definition: com_sub.f:221
bt_pol
subroutine bt_pol(betpol)
Definition: com_sub.f:747
bongri
subroutine bongri
Definition: com_sub.f:242
qst_b
subroutine qst_b
Definition: com_sub.f:576
axdef
subroutine axdef(rma, zma, psima, dpm)
Definition: com_sub.f:86
avr2_c
subroutine avr2_c(arr2, nro, nteta, arr1)
Definition: com_sub.f:136
avr2_k
subroutine avr2_k(arr2, nro, nteta, arr1)
Definition: com_sub.f:157
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
metcof
subroutine metcof(alp22, alp33, delsc, delv, ac0, skcen, cur_I,
Definition: com_sub.f:413
qvadin
real *8 function qvadin(x,
Definition: com_sub.f:733
avr1_c
real *8 function avr1_c(arr1_c, i)
Definition: com_sub.f:177
inter_q
subroutine inter_q(ps, nw, q_ps)
Definition: com_sub.f:1516