ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
Eq2_m.f
Go to the documentation of this file.
1  SUBROUTINE output(NGRID,betpol,zli3)
2 c
3  include 'double.inc'
4  parameter(nshp=10)
5  include 'param.inc'
6  include 'comblc.inc'
7 c
8  common/equili/ bettot, rmx, zzrmx, rmn, zzrmn,
9  * zmx, rrzmx, zmn, rrzmn,
10  * r0cen, z0cen, radm, aspect,
11  * eupper, elower, delup, dellw, bfvakc,
12  * z_li3
13 
14  common/comlop/ rxb(nbndp2),zxb(nbndp2),nxb
15  common/combrz/ bmr(nip,njp),bmz(nip,njp)
16  common /volpla/ vol_pl
17 c
18  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
19 c----------------------------------------------------------------------
20  if(ngrid.eq.0) go to 999
21 
22  rmx=0.
23  rmn=rmax
24 
25  zmx=zmin
26  zmn=zmax
27 
28 c kavin-----------------------------------------------------------------------
29  slen=0.
30 c------------------------------------------------------------------------------
31  do 100 i=1,nxb
32 
33 c kavin-----------------------------------------------------------------------
34  if(i .eq. nxb) then
35  rxb1=rxb(1)
36  zxb1=zxb(1)
37  else
38  rxb1=rxb(i+1)
39  zxb1=zxb(i+1)
40  endif
41  slen=slen+dsqrt((rxb(i)-rxb1)**2+(zxb(i)-zxb1)**2)
42 c------------------------------------------------------------------------------
43 
44  if(rxb(i).gt.rmx) then
45  rmx=rxb(i)
46  zzrmx=zxb(i)
47  endif
48 
49  if(rxb(i).lt.rmn) then
50  rmn=rxb(i)
51  zzrmn=zxb(i)
52  endif
53 
54  if(zxb(i).gt.zmx) then
55  zmx=zxb(i)
56  rrzmx=rxb(i)
57  endif
58 
59  if(zxb(i).lt.zmn) then
60  zmn=zxb(i)
61  rrzmn=rxb(i)
62  endif
63 
64  100 continue
65 
66  r0cen=0.5*(rmn+rmx)
67  z0cen=0.5*(zmn+zmx)
68  radm=0.5*(rmx-rmn)
69  aspect=r0cen/radm
70  eupper=dabs(zmx-zm)/radm
71  elower=dabs(zmn-zm)/radm
72  delup=(r0cen-rrzmx)/radm
73  dellw=(r0cen-rrzmn)/radm
74 
75 c write(3,*) '___________________________________'
76 c write(3,*) ' max. R on plasma boundary:'
77 c write(3,*) ' rmx = ',rmx
78 c write(3,*) ' Z(rmx)=',zzrmx
79 c write(3,*) ' min. R on plasma boundary:'
80 c write(3,*) ' rmn = ',rmn
81 c write(3,*) ' Z(rmn) ',zzrmn
82 c write(3,*) ' max. Z on plasma boundary:'
83 c write(3,*) ' zmx ',zmx
84 c write(3,*) ' R(zmx) ',rrzmx
85 c write(3,*) ' min. Z on plasma boundary:'
86 c write(3,*) ' zmn ',zmn
87 c write(3,*) ' R(zmn) ',rrzmn
88 c write(3,*) ' major plasma radius:'
89 c write(3,*) ' r0cen=',r0cen
90 c write(3,*) ' z0cen=',z0cen
91 c write(3,*) ' minor plasma radius:'
92 c write(3,*) ' radm= ',radm
93 c write(3,*) ' aspect ratio:'
94 c write(3,*) ' aspect=',aspect
95 c write(3,*) ' upper plasma elongation:'
96 c write(3,*) ' Eupper=',eupper
97 c write(3,*) ' lower plasma elongation:'
98 c write(3,*) ' Elower=',elower
99 c write(3,*) ' upper plasma triang.:'
100 c write(3,*) ' DELup =',delup
101 c write(3,*) ' lower plasma triang.:'
102 c write(3,*) ' DELlw =',dellw
103 
104  call magpol
105  bpolin=0.
106  psres=0.d0
107  volpl=0.
108  sqpl=0.
109  zc=0.
110 
111  do 200 i=2,ni1
112  do 210 j=2,nj1
113  iprr=ipr(i,j)
114  if(iprr.ne.1) go to 210
115 c... bp2ij=(bmz(i-1,j)**2+bmz(i,j)**2+
116 c... + bmr(i,j-1)**2+bmr(i,j)**2)*0.5
117 c... bpolin=bpolin+bp2ij*dri(i)*dzj(j)*r(i)
118 c kavin-----------------------------------------------------
119  volpl = volpl + r(i)*dri(i)*dzj(j)
120  sqpl = sqpl + dri(i)*dzj(j)
121  psres = psres + curf(i,j)*u(i,j) *dri(i)*dzj(j)
122  zc = zc + curf(i,j)*z(j) *dri(i)*dzj(j)/tok
123 c-----------------------------------------------------------
124  210 continue
125  200 continue
126 
127 c... zmx = zc
128  z0cen = zc
129 c... zli3 = 4.*pi*bpolin/(r0cen*tok*tok)
130 
131 c kavin-----------------------------------------------------
132 c psres=(psres/tok-up)/(0.6283185*r0cen*tok)
133  psres=(psres/tok-up)/(0.6283185*tok)*sqpl/volpl
134 c write(3,*) ' li(3) :'
135 c write(3,*) 'ZLI3=',zli3
136 
137  zli3=psres*7.8957
138  z_li3=psres*7.8957
139 c------------------------------------------------------------
140  pintg=0.
141  pintv=0.
142  volpl=0.
143 
144  do 300 i=2,ni1
145  do 310 j=2,nj1
146  iprr=ipr(i,j)
147  if(iprr.ne.1) go to 310
148  psi=un(i,j)
149  zpres=funppp(psi)
150  pintg=pintg+zpres*dri(i)*dzj(j)
151  pintv=pintv+zpres*r(i)*dri(i)*dzj(j)
152  volpl=volpl+ r(i)*dri(i)*dzj(j)
153  310 continue
154  300 continue
155 
156  vol_pl = volpl
157  paverg = pintv/volpl
158  betpol = 8.*pi*pintg/(tok*tok)
159  betpol = betpol*(um-up)*cnor
160 
161 c write(3,*) (paverg*(um-up)*cnor),' li3=', zli3, 'LENGTH :',slen
162 c ---- write(3,*) 'BETPOL :'
163 c ---- write(3,*) 'BETpol=',betpol
164 c ---- write(3,*) 'cnor,um,up',cnor,um,up
165 
166 c write(3,'(3(a,1pe12.4))')
167 c > ' Psmax =',(7.8957*um),
168 c > ' Psbound =',(7.8957*up),
169 c > ' Psxpt =',(7.8957*ux0)
170 
171  nsh = 1
172  xs(nsh) = r(imax)
173  ys(nsh) = z(jmax)
174  fun(nsh) = u(imax,jmax)
175 
176  do 800 k=-1,1
177 
178  i= imax+k
179 
180  do 810 l=-1,1
181 
182  j= jmax+l
183 
184  if(i.eq.imax .AND. j.eq.jmax) go to 810
185  nsh=nsh+1
186  xs(nsh)=r(i)
187  ys(nsh)=z(j)
188  fun(nsh)=u(i,j)
189 
190  810 continue
191  800 continue
192 
193  call deriv5(xs,ys,fun,nsh,5,dp)
194 
195  drr=dp(3)
196  drz=dp(4)
197  dzz=dp(5)
198 
199  tg2a=2.d0*drz/(drr-dzz)
200  cos2a=1.d0/dsqrt(1.d0+tg2a**2)
201  sin2a=cos2a*tg2a
202 
203  dxx=0.5d0*(drr+dzz)+0.5d0*cos2a*(drr-dzz)+sin2a*drz
204  dyy=0.5d0*(drr+dzz)+0.5d0*cos2a*(dzz-drr)-sin2a*drz
205 
206  bfcen=qcen*dsqrt(dxx*dyy)
207 
208  bettot=2.d0*paverg/(bfcen**2)
209 c bettot=bettot*(um-up)*cnor
210  bettot=paverg*(um-up)*cnor
211 
212  bettot=betpol*slen**2/sqpl/4d0/pi
213 
214 c write(6,*) 'bettot ',bettot
215 
216 c _______________________________________________
217 
218  fcen=bfcen*rm
219  fintg=funfff(1.d0)*(um-up)*cnor
220  fxcen=dsqrt(fcen**2-2.*fintg)
221 
222  bfvakc=fxcen
223 
224  999 continue
225 
226  return
227  end
228 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
230 c
231  subroutine loop
232 
233  include 'double.inc'
234  include 'param.inc'
235  include 'comblc.inc'
236 
237  common/comlop/ rxb(nbndp2), zxb(nbndp2), nxb
238  common/comus1/ rus1(nbndp2), zus1(nbndp2), nus1
239  common/comus2/ rus2(nbndp2), zus2(nbndp2), nus2
240  dimension rwork(nbndp2), zwork(nbndp2)
241  dimension us(nip,njp)
242 
243  xzer(x1,x2,v1,v2) = (x1*v2-x2*v1)/(v2-v1+1.d-8)
244 
245 c --- **********************************************
246 c write(5,*) '*** entrance of subr. loop ***'
247 c write(6,*) ' '
248 c write(6,*) '*** entrance of subr. loop ***'
249 c --- **********************************************
250 c --- un(i,j)=0 - plasma boundary
251 
252  !delunb=0.00200d0
253  delunb=0.00001d0
254 
255  do 2212 i=1,ni
256  do 2212 j=1,nj
257  us(i,j) = alpnew * (un(i,j) - 1.d0) + 1.d0
258  un(i,j) = un(i,j) - delunb
259  2212 continue
260 c----------------------------------------------------------
261 c------*****----------------------------------------------------
262  ig=1
263  jc=jmax
264  j=jmax
265 
266  888 continue
267 
268  do 885 i=imax,ni1
269 
270  if(un(i,j)*un(i+1,j).le.0.) then
271 
272  ic=i
273 
274  rxb(ig)=xzer(r(i+1),r(i),un(i+1,j),un(i,j))
275  zxb(ig)=z(j)
276 
277  go to 886
278  endif
279 
280  885 continue
281  !write(6,*) 'loop: first boundary point was not find'
282  stop
283 
284  886 continue
285 
286  ic1=ic
287  jc1=jc
288 
289 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
290 
291  lin=1
292 c---------------------------------------------------------------
293 c
294  100 ig=ig+1
295 c ----
296  i=ic
297  j=jc
298 c ----
299  if(un(i+1,j)*un(i,j).le.0. .AND. lin.ne.1) then
300 
301  ic = i
302  jc = j-1
303  lin = 3
304 
305  rxb(ig) = xzer(r(i+1),r(i),un(i+1,j),un(i,j))
306  zxb(ig) = z(j)
307 
308  elseif(un(i+1,j+1)*un(i+1,j).le.0..AND. lin.ne.2) then
309 
310  ic = i+1
311  jc = j
312  lin = 4
313 
314  rxb(ig) = r(i+1)
315  zxb(ig) = xzer(z(j+1),z(j),un(i+1,j+1),un(i+1,j))
316 
317  elseif(un(i+1,j+1)*un(i,j+1).le.0. .AND. lin.ne.3) then
318 
319  ic = i
320  jc = j+1
321  lin = 1
322 
323  rxb(ig) = xzer(r(i+1),r(i),un(i+1,j+1),un(i,j+1))
324  zxb(ig) = z(j+1)
325 
326  elseif(un(i,j+1)*un(i,j).le.0..AND. lin.ne.4) then
327 
328  ic = i-1
329  jc = j
330  lin = 2
331 
332  rxb(ig) = r(i)
333  zxb(ig) = xzer(z(j+1),z(j),un(i,j+1),un(i,j))
334 
335  endif
336 c------------------------
337 c --- control of finish
338 c
339  if( (jc.eq.jc1) .and. (ic.eq.ic1) ) then
340 
341  !ig=ig+1
342 
343  !rxb(ig)=rxb(1 )
344  !zxb(ig)=zxb(1 )
345 
346  go to 1212
347 
348  endif
349 c------------------------
350 c
351  go to 100
352 c------------------------
353 c
354  1212 nxb=ig
355 c---------------------------------------------------------------
356 c End of plasma boundary treatment
357 c---------------------------------------------------------------
358 c***************************************************************
359 
360  do i=1,ni
361  do j=1,nj
362  un(i,j) = un(i,j) + delunb
363  enddo
364  enddo
365  return
366  end
367 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
368 c**********************************************************
369  subroutine rpoint(ut,rt,zt,dp)
370 
371  include 'double.inc'
372  parameter(nshp=10)
373  include 'param.inc'
374  include 'comblc.inc'
375 
376  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
377 
378  fps(r_,r0,z_,z0,d_r,d_z,d_rr,d_zz,d_rz,ps_0,ps_x)=
379  * (ps_0-ps_x)+d_r*(r_-r0)+d_z*(z_-z0)+
380  * 0.5d0*d_rr*(r_-r0)**2+0.5d0*d_zz*(z_-z0)**2+d_rz*(z_-z0)*(r_-r0)
381 
382  fprim_r(r_,r0,z_,z0,d_r,d_rr,d_rz)=d_r+d_rr*(r_-r0)+d_rz*(z_-z0)
383 
384  niter=0
385  444 continue
386 
387  niter=niter+1
388 c---definition of cell, containig s-point
389 
390  do 300 i=1,ni1
391  itx=i
392  if( (rt.lt.r(i+1)) .AND. (rt.ge.r(i)) ) go to 301
393  300 continue
394  301 continue
395 
396  do 310 j=1,nj1
397  jtx=j
398  if( (zt.lt.z(j+1)) .AND. (zt.ge.z(j)) ) go to 311
399  310 continue
400  311 continue
401 
402 c write(*,*) 'icelt,jcelt',itx,jtx,rt,zt
403 
404 c---define nearest node
405 
406  sdmin=1.d6
407 
408  do 320 k=0,1
409  rr=r(itx+k)
410  do 325 l=0,1
411  zz=z(jtx+l)
412  dlx=dsqrt( (rr-rt)**2+(zz-zt)**2 )
413  if(dlx.lt.sdmin) then
414  sdmin=dlx
415  it=itx+k
416  jt=jtx+l
417  endif
418  325 continue
419  320 continue
420 
421  numit=0
422  555 continue
423  numit=numit+1
424 
425 c write(*,*) 'it,jt',it,jt
426 c
427 c write(*,*) r(it),z(jt)
428 
429  nsh=1
430  xs(nsh)=r(it)
431  ys(nsh)=z(jt)
432 
433 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
434  !! !!
435  fun(nsh)=u(it,jt) !!
436  !! !!
437 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
438 
439 
440  do 400 k=-1,1
441 
442  i= it+k
443 
444  do 410 l=-1,1
445 
446  j= jt+l
447 
448  if(k.eq.0 .AND. l.eq.0 ) go to 410
449  nsh=nsh+1
450  xs(nsh)=r(i)
451  ys(nsh)=z(j)
452 
453 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
454  !! !!
455  fun(nsh)=u(i,j) !!
456  !! !!
457 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
458 
459  410 continue
460  400 continue
461 
462  call deriv5(xs,ys,fun,nsh,5,dp)
463 
464  r_0=xs(1)
465  z_0=ys(1)
466  rrt=rt
467  zzt=zt
468 
469  d_r=dp(1)
470  d_z=dp(2)
471  d_rr=dp(3)
472  d_zz=dp(4)
473  d_rz=dp(5)
474  ps_0=fun(1)
475  ps_x=ut
476 
477 
478 
479  do inw=1,10
480 
481  rt = rrt - fps(rrt,r_0,zt,z_0,d_r,d_z,d_rr,d_zz,d_rz,ps_0,ps_x)/
482  * fprim_r(rrt,r_0,zt,z_0,d_r,d_rr,d_rz)
483  accr=dabs(rt-rrt)
484  if(accr.lt.1.d-6) exit
485  rrt=rt
486 
487  enddo
488 
489 
490  if(niter.lt.4) go to 444
491 
492  return
493  end
494 
495 
496  subroutine zpoint(ut,rt,zt,dp)
497 
498  include 'double.inc'
499  parameter(nshp=10)
500  include 'param.inc'
501  include 'comblc.inc'
502 
503  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
504 
505  fps(r_,r0,z_,z0,d_r,d_z,d_rr,d_zz,d_rz,ps_0,ps_x)=
506  * (ps_0-ps_x)+d_r*(r_-r0)+d_z*(z_-z0)+
507  * 0.5d0*d_rr*(r_-r0)**2+0.5d0*d_zz*(z_-z0)**2+d_rz*(z_-z0)*(r_-r0)
508 
509  fprim_z(r_,r0,z_,z0,d_z,d_zz,d_rz)=d_z+d_zz*(z_-z0)+d_rz*(r_-r0)
510 
511  niter=0
512  444 continue
513 
514  niter=niter+1
515 c---definition of cell, containig s-point
516 
517  do 300 i=1,ni1
518  itx=i
519  if( (rt.lt.r(i+1)) .AND. (rt.ge.r(i)) ) go to 301
520  300 continue
521  301 continue
522 
523  do 310 j=1,nj1
524  jtx=j
525  if( (zt.lt.z(j+1)) .AND. (zt.ge.z(j)) ) go to 311
526  310 continue
527  311 continue
528 
529 c write(*,*) 'icelt,jcelt',itx,jtx,rt,zt
530 
531 c---define nearest node
532 
533  sdmin=rmax-rmin+zmax-zmin
534 
535  do 320 k=0,1
536  rr=r(itx+k)
537  do 325 l=0,1
538  zz=z(jtx+l)
539  dlx=dsqrt( (rr-rt)**2+(zz-zt)**2 )
540  if(dlx.lt.sdmin) then
541  sdmin=dlx
542  it=itx+k
543  jt=jtx+l
544  endif
545  325 continue
546  320 continue
547 
548  numit=0
549  555 continue
550  numit=numit+1
551 
552 c write(*,*) 'it,jt',it,jt
553 c
554 c write(*,*) r(it),z(jt)
555 
556  nsh=1
557  xs(nsh)=r(it)
558  ys(nsh)=z(jt)
559 
560 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
561  !! !!
562  fun(nsh)=u(it,jt) !!
563  !! !!
564 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
565 
566 
567  do 400 k=-1,1
568 
569  i= it+k
570 
571  do 410 l=-1,1
572 
573  j= jt+l
574 
575  if(k.eq.0 .AND. l.eq.0 ) go to 410
576  nsh=nsh+1
577  xs(nsh)=r(i)
578  ys(nsh)=z(j)
579 
580 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
581  !! !!
582  fun(nsh)=u(i,j) !!
583  !! !!
584 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
585 
586  410 continue
587  400 continue
588 
589  call deriv5(xs,ys,fun,nsh,5,dp)
590 
591  r_0=xs(1)
592  z_0=ys(1)
593 
594  rrt=rt
595  zzt=zt
596 
597  d_r=dp(1)
598  d_z=dp(2)
599  d_rr=dp(3)
600  d_zz=dp(4)
601  d_rz=dp(5)
602  ps_0=fun(1)
603  ps_x=ut
604 
605 
606 
607  do inw=1,10
608 
609  zt = zzt - fps(rt,r_0,zzt,z_0,d_r,d_z,d_rr,d_zz,d_rz,ps_0,ps_x)/
610  * fprim_z(rt,r_0,zzt,z_0,d_z,d_zz,d_rz)
611  accr=dabs(zt-zzt)
612  if(accr.lt.1.d-6) exit
613  zzt=zt
614 
615  enddo
616 
617  if(niter.lt.4) go to 444
618 
619  return
620  end
621 
622 
623 c**********************************************************
624  subroutine gap(ut,rt,zt,rg,zg)
625 
626  include 'double.inc'
627  parameter(nshp=10)
628  include 'param.inc'
629  include 'comblc.inc'
630 
631  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
632 
633  fps(r_,r0,z_,z0,d_r,d_z,d_rr,d_zz,d_rz,ps_0,ps_x)=
634  * (ps_0-ps_x)+d_r*(r_-r0)+d_z*(z_-z0)+
635  * 0.5d0*d_rr*(r_-r0)**2+0.5d0*d_zz*(z_-z0)**2+d_rz*(z_-z0)*(r_-r0)
636 
637  fprim_z(r_,r0,z_,z0,d_z,d_zz,d_rz)=d_z+d_zz*(z_-z0)+d_rz*(r_-r0)
638 
639  fprim_r(r_,r0,z_,z0,d_r,d_rr,d_rz)=d_r+d_rr*(r_-r0)+d_rz*(z_-z0)
640 
641  k_iter=0
642 
643 
644  333 continue
645 
646  k_iter=k_iter+1
647 
648  n_iter=0
649  itx_old=0
650  jtx_old=0
651 
652 444 n_iter=n_iter+1
653 c---definition of cell, containig s-point
654 
655  do 300 i=1,ni1
656  itx=i
657  if( (rt.lt.r(i+1)) .AND. (rt.ge.r(i)) ) go to 301
658  300 continue
659  301 continue
660 
661  do 310 j=1,nj1
662  jtx=j
663  if( (zt.lt.z(j+1)) .AND. (zt.ge.z(j)) ) go to 311
664  310 continue
665  311 continue
666 
667 c write(*,*) 'icelt,jcelt',itx,jtx,rt,zt
668 
669 c---define nearest node
670 
671  sdmin=1.d6
672 
673  do 320 k=0,1
674  rr=r(itx+k)
675  do 325 l=0,1
676  zz=z(jtx+l)
677  dlx=dsqrt( (rr-rt)**2+(zz-zt)**2 )
678  if(dlx.lt.sdmin) then
679  sdmin=dlx
680  it=itx+k
681  jt=jtx+l
682  endif
683  325 continue
684  320 continue
685 
686 c write(*,*) 'it,jt',it,jt
687 c
688 c write(*,*) r(it),z(jt)
689 
690  nsh=1
691  xs(nsh)=r(it)
692  ys(nsh)=z(jt)
693 
694 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
695  !! !!
696  fun(nsh)=u(it,jt) !!
697  !! !!
698 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
699 
700 
701  do 400 k=-1,1
702 
703  i= it+k
704 
705  do 410 l=-1,1
706 
707  j= jt+l
708 
709  if(k.eq.0 .AND. l.eq.0 ) go to 410
710  nsh=nsh+1
711  xs(nsh)=r(i)
712  ys(nsh)=z(j)
713 
714 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
715  !! !!
716  fun(nsh)=u(i,j) !!
717  !! !!
718 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
719 
720  410 continue
721  400 continue
722 
723  call deriv5(xs,ys,fun,nsh,5,dp)
724 
725  r_0=xs(1)
726  z_0=ys(1)
727 
728  rrt=rt
729  zzt=zt
730 
731  d_r=dp(1)
732  d_z=dp(2)
733  d_rr=dp(3)
734  d_zz=dp(4)
735  d_rz=dp(5)
736  ps_0=fun(1)
737  ps_x=ut
738 
739  if(k_iter*n_iter.eq.1) then
740  psi_cp=ps_0+d_r*(rg-r_0)+d_z*(zg-z_0)+0.5d0*d_rr*(rg-r_0)**2
741  & +0.5d0*d_zz*(zg-z_0)**2+d_rz*(rg-r_0)*(zg-z_0)
742  endif
743 
744  do inw=1,10
745 
746  dpsidr=fprim_r(rt,r_0,zt,z_0,d_r,d_rr,d_rz)
747  dpsidz=fprim_z(rt,r_0,zt,z_0,d_z,d_zz,d_rz)
748  grad=dsqrt(dpsidr**2+dpsidz**2)
749 
750  dpsi=-fps(rt,r_0,zt,z_0,d_r,d_z,d_rr,d_zz,d_rz,ps_0,ps_x)
751 
752  d_ro=dpsi/grad
753 
754  rtp = rt + d_ro*dpsidr/grad
755  ztp = zt + d_ro*dpsidz/grad
756 
757  accr=dabs(zt-ztp)+dabs(rt-rtp)
758  if(accr.lt.1.d-6) exit
759 
760  rt = rtp
761  zt = ztp
762 
763  enddo
764 
765  if(itx.eq.itx_old .AnD. jtx.eq.jtx_old) go to 555
766  if(n_iter.eq.5) go to 555
767 
768  itx_old=itx
769  jtx_old=jtx
770 
771  go to 444
772  555 continue
773 
774  taur=-dpsidz/dsqrt(dpsidr**2+dpsidz**2)
775  tauz= dpsidr/dsqrt(dpsidr**2+dpsidz**2)
776 
777  dgt_r=rg-rt
778  dgt_z=zg-zt
779  scal_pro=dgt_r*taur+dgt_z*tauz
780 
781  wght=0.8d0
782 
783  rt=rt+scal_pro*taur*wght
784  zt=zt+scal_pro*tauz*wght
785 
786  if(dabs(scal_pro).lt.1.d-6) go to 777
787  if(k_iter.eq.5) go to 777
788  go to 333
789 
790  777 continue
791 
792  ut=psi_cp
793 
794  return
795  end
796 
797 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
798  subroutine f_gap(ut,rt,zt,rg,zg)
799 
800  include 'double.inc'
801  parameter(nshp=20)
802  include 'parevo.inc'
803  parameter(nkp=njlim)
804  include 'dim.inc'
805  include 'compol.inc'
806  include 'compol_add.inc'
807 
808  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
809 
810  fps(r_,r0,z_,z0,d_r,d_z,d_rr,d_zz,d_rz,ps_0,ps_x)=
811  * (ps_0-ps_x)+d_r*(r_-r0)+d_z*(z_-z0)+
812  * 0.5d0*d_rr*(r_-r0)**2+0.5d0*d_zz*(z_-z0)**2+d_rz*(z_-z0)*(r_-r0)
813 
814  fprim_z(r_,r0,z_,z0,d_z,d_zz,d_rz)=d_z+d_zz*(z_-z0)+d_rz*(r_-r0)
815 
816  fprim_r(r_,r0,z_,z0,d_r,d_rr,d_rz)=d_r+d_rr*(r_-r0)+d_rz*(z_-z0)
817 
818  k_iter=0
819 
820 
821  333 continue
822 
823  k_iter=k_iter+1
824 
825  n_iter=0
826  itx_old=0
827  jtx_old=0
828 
829 444 n_iter=n_iter+1
830 c---definition of cell, containig s-point
831 
832  call numcel(rt,zt,itx,jtx)
833 
834 c write(*,*) 'icelt,jcelt',itx,jtx,rt,zt
835  if(itx.eq.nr) then
836  !pause 'f_gap:gap point in bound cell'
837  return
838  endif
839 c---define nearest node
840 
841  sdmin=1.d6
842 
843  do 320 k=0,1
844  do 325 l=0,1
845  rr=r(itx+k,jtx+l)
846  zz=z(itx+k,jtx+l)
847  dlx=dsqrt( (rr-rt)**2+(zz-zt)**2 )
848  if(dlx.lt.sdmin) then
849  sdmin=dlx
850  it=itx+k
851  jt=jtx+l
852  endif
853  325 continue
854  320 continue
855 
856  if(jt.eq.nt) jt=2
857  if(jt.eq.1) jt=nt-1
858  if(it.eq.nr) it=nr-1
859 
860  nsh=1
861  xs(nsh)=r(it,jt)
862  ys(nsh)=z(it,jt)
863 
864 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
865  !! !!
866  fun(nsh)=psi(it,jt) !!
867  !! !!
868 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
869 
870 
871  do 400 k=-1,1
872 
873  i= it+k
874 
875  do 410 l=-1,1
876 
877  j= jt+l
878 
879  if(k.eq.0 .AND. l.eq.0 ) go to 410
880  nsh=nsh+1
881  xs(nsh)=r(i,j)
882  ys(nsh)=z(i,j)
883 
884 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
885  !! !!
886  fun(nsh)=psi(i,j) !!
887  !! !!
888 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
889 
890  410 continue
891  400 continue
892 
893  call deriv5(xs,ys,fun,nsh,5,dp)
894 
895  r_0=xs(1)
896  z_0=ys(1)
897 
898  rrt=rt
899  zzt=zt
900 
901  d_r=dp(1)
902  d_z=dp(2)
903  d_rr=dp(3)
904  d_zz=dp(4)
905  d_rz=dp(5)
906  ps_0=fun(1)
907  ps_x=ut
908 
909  if(k_iter*n_iter.eq.1) then
910  psi_cp=ps_0+d_r*(rg-r_0)+d_z*(zg-z_0)+0.5d0*d_rr*(rg-r_0)**2
911  & +0.5d0*d_zz*(zg-z_0)**2+d_rz*(rg-r_0)*(zg-z_0)
912  endif
913 
914  do inw=1,10
915 
916  dpsidr=fprim_r(rt,r_0,zt,z_0,d_r,d_rr,d_rz)
917  dpsidz=fprim_z(rt,r_0,zt,z_0,d_z,d_zz,d_rz)
918  grad=dsqrt(dpsidr**2+dpsidz**2)
919 
920  dpsi=-fps(rt,r_0,zt,z_0,d_r,d_z,d_rr,d_zz,d_rz,ps_0,ps_x)
921 
922  d_ro=dpsi/grad
923 
924  rtp = rt + d_ro*dpsidr/grad
925  ztp = zt + d_ro*dpsidz/grad
926 
927  accr=dabs(zt-ztp)+dabs(rt-rtp)
928  if(accr.lt.1.d-6) exit
929 
930  rt = rtp
931  zt = ztp
932 
933  enddo
934 
935  if(itx.eq.itx_old .AnD. jtx.eq.jtx_old) go to 555
936  if(n_iter.eq.5) go to 555
937 
938  itx_old=itx
939  jtx_old=jtx
940 
941  go to 444
942  555 continue
943 
944  taur=-dpsidz/dsqrt(dpsidr**2+dpsidz**2)
945  tauz= dpsidr/dsqrt(dpsidr**2+dpsidz**2)
946 
947  dgt_r=rg-rt
948  dgt_z=zg-zt
949  scal_pro=dgt_r*taur+dgt_z*tauz
950 
951  wght=0.8d0
952 
953  rt=rt+scal_pro*taur*wght
954  zt=zt+scal_pro*tauz*wght
955 
956  if(dabs(scal_pro).lt.1.d-6) go to 777
957  if(k_iter.eq.5) go to 777
958  go to 333
959 
960  777 continue
961 
962  ut=psi_cp
963 
964  return
965  end
966 
967 
968 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
969  subroutine adapol
970 c
971  include 'double.inc'
972  parameter(ntet=66,npsi=33)
973  include 'param.inc'
974  include 'comblc.inc'
975  include 'urs.inc'
976 c
977  common/comlop/ rxb(nbndp2),zxb(nbndp2),nxb
978  real*8 teta0(nbndp2),tet(ntet),rpol(ntet),zpol(ntet)
979  real*8 psipol(npsi),pppol(npsi),fppol(npsi)
980 c
981  rbmax=0.
982  rbmin=rmax
983 
984  zbmax=zmin
985  zbmin=zmax
986 
987  do 77 i=1,nxb-1
988 
989  rbmax=dmax1(rxb(i),rbmax)
990  rbmin=dmin1(rxb(i),rbmin)
991 
992  zbmax=dmax1(zxb(i),zbmax)
993  zbmin=dmin1(zxb(i),zbmin)
994 
995  77 continue
996 
997  rcen=0.5*(rbmax+rbmin)
998  zcen=0.5*(zbmax+zbmin)
999 
1000  do 10 i=1,nxb
1001 
1002  drx=rxb(i)-rcen
1003  dzx=zxb(i)-zcen
1004 
1005  tetp=datan(dzx/drx)
1006  if(drx.lt.0) tetp=tetp+pi
1007  if(tetp.lt.0) tetp=tetp+2.*pi
1008  teta0(i)=tetp
1009 
1010 ccc write(6,*) 'teta0(i) i',teta0(i),i
1011 
1012  10 continue
1013 
1014  tet(2)=0.
1015  dtet=2.*pi/(ntet-2.)
1016 
1017  do 100 i=3,ntet-1
1018  tet(i)=tet(i-1)+dtet
1019 
1020  100 continue
1021 
1022  tet(ntet)=2.*pi
1023  tet(1)=2.*pi-dtet
1024 
1025  do 200 j=1,ntet
1026  tetp=tet(j)
1027 
1028 ccc write(6,*) 'tetp j',tetp,j
1029 !!!!! go to 2222 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1030 
1031  do 2000 i=2,nxb
1032  tetx1=teta0(i-1)
1033  tetx2=teta0(i)
1034 
1035  if(tetx1.gt.tetx2) then
1036  if(tetp.gt.pi) then
1037  tetx2=tetx2+2.*pi
1038  elseif(tetp.lt.pi) then
1039  tetx1=tetx1-2.*pi
1040  endif
1041  endif
1042 c write(6,*) 'tetp tetx1 tetx2'
1043 c write(6,*) tetp,tetx1,tetx2
1044 
1045  if(tetp.gt.tetx1 .AND. tetp.le.tetx2) then
1046 
1047 ccc write(6,*) 'tetp tetx1 tetx2'
1048 ccc write(6,*) tetp,tetx1,tetx2
1049 
1050  rpol(j)=( (tetp-tetx1)*rxb(i)+(tetx2-tetp)*rxb(i-1) )/
1051  / (tetx2-tetx1)
1052 
1053  zpol(j)=((tetp-tetx1)*zxb(i)+(tetx2-tetp)*zxb(i-1))/
1054  / (tetx2-tetx1)
1055 
1056  endif
1057 
1058  2000 continue
1059 c2222 !! rrrr=2. !!!!!!!!!!!!
1060  !! rpol(j)=rcen+rrrr*dcos(tetp) !!!!!!!!!!!!
1061  !! zpol(j)=zcen+rrrr*dsin(tetp) !!!!!!!!!!!!
1062 
1063  200 continue
1064 
1065  do 300 i=1,npsi
1066 
1067  psipol(i)=1.-(i-1)/(npsi-1.)
1068 
1069  pppol(i)=cnor*funpp(psipol(i))
1070  fppol(i)=cnor*funfp(psipol(i))
1071 
1072  300 continue
1073 
1074  write(fname,'(a,a)') path(1:kname),'adpol.dat'
1075  open(1,file=fname,form='formatted')
1076  !open(1,file='adpol.dat',form='formatted')
1077 
1078  write(1,*) npsi,ntet,rm,zm,rcen,zcen
1079 
1080  write(1,*) (rpol(i),i=1,ntet)
1081  write(1,*) (zpol(i),i=1,ntet)
1082 
1083  write(1,*) (psipol(i),i=1,npsi)
1084  write(1,*) (pppol(i),i=1,npsi)
1085  write(1,*) (fppol(i),i=1,npsi)
1086 
1087  close(1)
1088 
1089  return
1090  end
1091 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1092 c
1093  subroutine magpol
1094 c
1095  include 'double.inc'
1096  include 'param.inc'
1097  include 'comblc.inc'
1098  common/combrz/ bmr(nip,njp),bmz(nip,njp)
1099 c
1100  do 10 i=1,ni
1101  do 10 j=1,nj1
1102 
1103  bmr(i,j)=(u(i,j+1)-u(i,j))/dz(j)/r(i)
1104 
1105  10 continue
1106 
1107  do 20 i=1,ni1
1108  do 20 j=1,nj
1109 
1110  bmz(i,j)=(u(i+1,j)-u(i,j))/dr(i)/r12(i)
1111 
1112  20 continue
1113 
1114  return
1115  end
1116 c***************************************************************
1117 c
1118  subroutine btpol(betpol)
1119 c
1120  include 'double.inc'
1121  include 'param.inc'
1122  include 'comblc.inc'
1123 c
1124  !!amu0=0.4d0*pi
1125  pintg=0.d0
1126  do 300 i=2,ni1
1127  do 310 j=2,nj1
1128  iprr=ipr(i,j)
1129  if(iprr.ne.1) go to 310
1130  psi=un(i,j)
1131  zpres=funppp(psi)
1132  pintg=pintg+zpres*dri(i)*dzj(j)
1133  310 continue
1134  300 continue
1135  betpol=8.d0*pi*pintg/(tok*tok)
1136  betpol=betpol*(um-up)*cnor/amu0**2
1137 c
1138  return
1139  end
1140 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1141  subroutine bttor(bettor)
1142 c
1143  include 'double.inc'
1144  include 'param.inc'
1145  include 'comblc.inc'
1146  common/equili/ bettot, rmx, zzrmx, rmn, zzrmn,
1147  * zmx, rrzmx, zmn, rrzmn,
1148  * r0cen, z0cen, radm, aspect,
1149  * eupper, elower, delup, dellw, bfvakc,
1150  * z_li3
1151 
1152  !!amu0=0.4d0*pi
1153  pintg=0.d0
1154  volpl=0.d0
1155 
1156  do 300 i=2,ni1
1157  do 310 j=2,nj1
1158  iprr=ipr(i,j)
1159  if(iprr.ne.1) go to 310
1160  psi=un(i,j)
1161  zpres=funppp(psi)
1162  pintg=pintg+zpres*dri(i)*dzj(j)*r(i)
1163  volpl=volpl+ dri(i)*dzj(j)*r(i)
1164  310 continue
1165  300 continue
1166 
1167  paverg = pintg*(um-up)*cnor/volpl
1168 
1169  b0_cen=b0ax*r0ax/r0cen
1170 
1171  bettor=2.d0*paverg/(b0_cen*b0_cen)
1172 
1173  tok_nor=tok/(radm*b0_cen)
1174 
1175  bet_nor=bettor/tok_nor*100.d0
1176 
1177  return
1178  end
1179 c***************************************************************
1180 c
1181  SUBROUTINE arcle ( NBAS, RBAS, ZBAS, ARC, ARCMAX )
1182 c
1183  include 'double.inc'
1184 c
1185  dimension rbas(*), zbas(*), arc(*)
1186 c
1187  sqrt(x) = dsqrt(x)
1188 c
1189  arc(1) = 0.d0
1190  DO 1 i=2,nbas
1191  ds = sqrt( ( rbas(i) - rbas(i-1) )**2 +
1192  * ( zbas(i) - zbas(i-1) )**2 )
1193  1 arc(i) = arc(i-1) + ds
1194  arcmax = arc(nbas)
1195  DO 2 i=1,nbas
1196  2 arc(i) = arc(i) / arcmax
1197 c
1198  RETURN
1199  END
1200 c******************************************************************
1201 c-----------------------------------------------------------------------
1202 c
1203  SUBROUTINE tcontr(KAS,A,JA,IA,NN,Z,B,IPRT,IGRF,BZ,ERR)
1204 c
1205  IMPLICIT REAL*8(a-h,o-z)
1206 c
1207 ccccccccc
1208  dimension a(*),ja(*),ia(*),z(*),b(*),bz(*)
1209 c iprt=0, igrf=0 --- poptpeta het
1210 ccccccccc
1211  abs(x)=dabs(x)
1212 c
1213  nout = 17
1214 c
1215 c ppobepka *drv
1216 c
1217  IF(kas.EQ.1) THEN
1218  CALL sarvec(a,ja,ia,nn,z,b)
1219  END IF
1220 c
1221  IF(kas.EQ.0) THEN
1222  CALL rarvec(a,ja,ia,nn,z,b)
1223  END IF
1224 c
1225  err=0.
1226  DO 201 lx=1,nn
1227  errlx=abs(b(lx)-bz(lx))
1228  IF(err.LT.errlx) err=errlx
1229 201 CONTINUE
1230 c
1231 c print matrix
1232 c
1233  IF(iprt.EQ.1) THEN
1234 c
1235  !WRITE(NOUT,314)
1236  !WRITE(NOUT,299) ERR
1237 c WRITE(nout,404)
1238 c WRITE(nout,102)(z(k),k=1,nn)
1239 c
1240 c WRITE(nout,314)
1241 c WRITE(nout,405)
1242 c WRITE(nout,102)(bz(k),k=1,nn)
1243 c
1244 c WRITE(nout,314)
1245 c WRITE(nout,406)
1246 c WRITE(nout,102)(b(k),k=1,nn)
1247 c WRITE(nout,314)
1248 c
1249 c+++ CALL aprt(nn,a,ia,ja)
1250 c
1251  ! WRITE(NOUT,314)
1252 c
1253  END IF
1254 c
1255 c poptpet matpitsy
1256 c
1257  IF(igrf.EQ.1) THEN
1258 c
1259 c+++ CALL agrf(nn,ia,ja)
1260 c
1261  END IF
1262 c
1263 299 FORMAT(/5x,'PRINT - TCONTR',2x,'ERROR=',e12.5)
1264 404 FORMAT(/5x,'Z - AFTER *DRV')
1265 405 FORMAT(/5x,'B - BEFORE *DRV')
1266 406 FORMAT(/5x,'B - AFTER R(S)ARVEC')
1267 c
1268 102 FORMAT(8e12.5)
1269 314 FORMAT(2x,'* * * * * * * * * * * * * * * * * * * * * * * * * *')
1270 c
1271  RETURN
1272  END
1273 c
1274 c****************************
1275 c
1276 c
1277  SUBROUTINE rarvec(A,JA,IA,N,X,Y)
1278  IMPLICIT REAL*8(a-h,o-z)
1279  dimension a(1),ja(1),ia(1),x(1),y(1)
1280  DO 1 i=1,n
1281  y(i)=0.0
1282  k1=ia(i)
1283  k2=ia(i+1)-1
1284  DO 1 j=k1,k2
1285  y(i)=y(i)+a(j)*x(ja(j))
1286 1 CONTINUE
1287 c
1288  RETURN
1289  END
1290 c****************************
1291 c
1292 c
1293 c****************************
1294  SUBROUTINE sarvec(A,JA,IA,N,X,Y)
1295  IMPLICIT REAL*8(a-h,o-z)
1296  dimension a(1),ja(1),ia(1),x(1),y(1)
1297 c
1298  DO 2 i=1,n
1299 2 y(i)=0.0
1300 c
1301  DO 1 i=1,n
1302 c
1303  k1=ia(i)
1304  k2=ia(i+1)-1
1305 c
1306  IF(ja(k1).EQ.i) THEN
1307  y(i)=y(i)+a(k1)*x(i)
1308  k1=k1+1
1309  END IF
1310 c
1311  IF(k1.LE.k2) THEN
1312  DO 3 j=k1,k2
1313 c
1314  y(i)=y(i)+a(j)*x(ja(j))
1315 c
1316  y(ja(j))=y(ja(j))+a(j)*x(i)
1317 3 CONTINUE
1318  END IF
1319 c
1320 1 CONTINUE
1321 c
1322  RETURN
1323  END
1324 c****************************
1325 c
1326 c
1327  SUBROUTINE aprt(N,A,IA,JA)
1328  IMPLICIT REAL*8(a-h,o-z)
1329  dimension a(1),ia(1),ja(1)
1330 c
1331  nout = 17
1332 c
1333 c DO 1 i=1,n
1334  DO 1 i=1,300
1335  k1=ia(i)
1336  k2=ia(i+1)-1
1337 c
1338  asum=0.
1339  DO 2 k=k1,k2
1340 2 asum=asum+a(k)
1341 c
1342  !WRITE(NOUT,102)
1343  !WRITE(NOUT,100)I,K1
1344  !WRITE(NOUT,104) ASUM
1345  !WRITE(NOUT,101)(JA(K),K=K1,K2)
1346  !WRITE(NOUT,103)(A(K),K=K1,K2)
1347 c
1348 1 CONTINUE
1349 c
1350 100 FORMAT(/5x,'IA(',i7,')=',i10)
1351 101 FORMAT(9i7)
1352 103 FORMAT(9e10.3)
1353 104 FORMAT(/5x,'ASUM=',e12.5)
1354 102 FORMAT(2x,'* * * * * * * * * * * * * * * * * * * * * * * * * *')
1355 c
1356 c
1357  RETURN
1358  END
1359 c*****************************
1360  SUBROUTINE agrf(N,IA,JA)
1361 c
1362  INTEGER n, ia(1), ja(1)
1363  CHARACTER*1 star , blanc , string(132)
1364  DATA star/'*'/, blanc/' '/
1365 c
1366  nout = 17
1367 c
1368  DO 3 i=1,n
1369  DO 1 j=1,132
1370  1 string(j)=blanc
1371  DO 2 j=ia(i),ia(i+1)-1
1372  2 string(ja(j))=star
1373  !WRITE(NOUT,'(132A1)')STRING
1374  3 continue
1375  RETURN
1376  END
1377 c**********************************************************
1378  SUBROUTINE wrd
1379 
1380  include 'double.inc'
1381  include 'param.inc'
1382  include 'comblc.inc'
1383 c
1384  common/comlop/ rxb(nbndp2), zxb(nbndp2), nxb
1385  common/nostep/ kstep1
1386 
1387 c
1388 c---------------------------------------------------------------
1389 c
1390 c++++ keywri = mod( kstep1 , 2 )
1391  keywri = 0
1392 c
1393 c---------------------------------------------------------------
1394 c
1395 c if( keywri .eq. 0 ) then
1396 c
1397  write(fname,'(a,a)') path(1:kname),'out.wr'
1398  open(1,file=fname,form='formatted')
1399  !open(1,file='out.wr')
1400 
1401  write(1,*) ni,nj,ni1,nj1,ni2,nj2,nxb
1402 
1403  write(1,*) (r(i),i=1,ni)
1404  write(1,*) (z(j),j=1,nj)
1405  write(1,*) ((u(i,j),i=1,ni),j=1,nj)
1406  write(1,*) ((curf(i,j),i=1,ni),j=1,nj)
1407  write(1,*) ((ipr(i,j),i=1,ni),j=1,nj)
1408  write(1,*) rm,zm,um,rx0,zx0,ux0,up
1409  write(1,*) (rxb(ig),ig=1,nxb)
1410  write(1,*) (zxb(ig),ig=1,nxb)
1411 
1412  write(1,*) ((ui(i,j),i=1,ni),j=1,nj)
1413  write(1,*) ((ue(i,j),i=1,ni),j=1,nj)
1414  write(1,*) ((un(i,j),i=1,ni),j=1,nj)
1415  close(1)
1416 
1417  write(fname,'(a,a)') path(1:kname),'recbon.wr'
1418  open(1,file=fname,form='formatted')
1419 
1420  write(1,*) nxb
1421  do ig=1,nxb
1422  write(1,*) rxb(ig),zxb(ig)
1423  enddo
1424 
1425  close(1)
1426 
1427 
1428 c else
1429 c
1430 c open(1,file='out1.wr')
1431 c
1432 c write(1,*) ni,nj,ni1,nj1,ni2,nj2,nxb
1433 c
1434 c write(1,*) (r(i),i=1,ni)
1435 c write(1,*) (z(j),j=1,nj)
1436 c write(1,*) ((u(i,j),i=1,ni),j=1,nj)
1437 c write(1,*) ((curf(i,j),i=1,ni),j=1,nj)
1438 c write(1,*) ((ipr(i,j),i=1,ni),j=1,nj)
1439 c write(1,*) rm,zm,um,rx0,zx0,ux0,up
1440 c write(1,*) (rxb(ig),ig=1,nxb)
1441 c write(1,*) (zxb(ig),ig=1,nxb)
1442 c
1443 c write(1,*) ((ui(i,j),i=1,ni),j=1,nj)
1444 c write(1,*) ((ue(i,j),i=1,ni),j=1,nj)
1445 c write(1,*) ((un(i,j),i=1,ni),j=1,nj)
1446 c close(1)
1447 c
1448 c end if
1449 c---------------------------------------------------------------
1450  RETURN
1451  END
1452 
1453  SUBROUTINE wrrec
1454 
1455  include 'double.inc'
1456  include 'param.inc'
1457  include 'comblc.inc'
1458 
1459  write(fname,'(a,a)') path(1:kname),'rect.wr'
1460  open(1,file=fname,form='formatted')
1461  !open(1,file='rect.wr',form='formatted')
1462 
1463  write(1,*) ni,nj,ni1,nj1,ni2,nj2,imax,jmax
1464 
1465  write(1,*) (r(i),i=1,ni)
1466  write(1,*) (z(j),j=1,nj)
1467  write(1,*) ((u(i,j),i=1,ni),j=1,nj)
1468  write(1,*) ((ue(i,j),i=1,ni),j=1,nj)
1469  write(1,*) ((un(i,j),i=1,ni),j=1,nj)
1470  write(1,*) ((ipr(i,j),i=1,ni),j=1,nj)
1471  write(1,*) rm,zm,um,rx0,zx0,ux0,up,qcen,b0ax,r0ax
1472  write(1,*) rx1,zx1,rx2,zx2
1473  write(1,*) rmax,zmax,rmin,zmin
1474 
1475  close(1)
1476  RETURN
1477  END
1478 
1479 
1480  SUBROUTINE wrdfmv(numwr,time)
1481 
1482  include 'double.inc'
1483  include 'param.inc'
1484  include 'comblc.inc'
1485 
1486  common/comlop/ rxb(nbndp2), zxb(nbndp2), nxb
1487  common/nostep/ kstep1
1488  character*40 str,dummy
1489 
1490  write(fname,'(a,a)') path(1:kname),'nw.wr'
1491  open(1,file=fname,form='formatted')
1492  !open(1,file='nw.wr',form='formatted')
1493  write(1,*) numwr
1494  close(1)
1495 
1496  if(numwr.lt.10) then
1497  write(str,'(a,a,i1,a)') path(1:kname),'file',numwr,'.wr'
1498  else
1499  if(numwr.lt.100) then
1500  write(str,'(a,a,i2,a)') path(1:kname),'file',numwr,'.wr'
1501  else
1502  write(str,'(a,a,i3,a)') path(1:kname),'file',numwr,'.wr'
1503  endif
1504  endif
1505 
1506  open(1,file=str,form='formatted')
1507 
1508  write(1,*) ni,nj,ni1,nj1,ni2,nj2,nxb
1509  write(1,*) time
1510 
1511  write(1,*) (r(i),i=1,ni)
1512  write(1,*) (z(j),j=1,nj)
1513  write(1,*) ((u(i,j),i=1,ni),j=1,nj)
1514  write(1,*) ((curf(i,j),i=1,ni),j=1,nj)
1515  write(1,*) ((ipr(i,j),i=1,ni),j=1,nj)
1516  write(1,*) rm,zm,um,rx0,zx0,ux0,up
1517  write(1,*) (rxb(ig),ig=1,nxb)
1518  write(1,*) (zxb(ig),ig=1,nxb)
1519 
1520  !write(1,*) ((ui(i,j),i=1,ni),j=1,nj)
1521  !write(1,*) ((ue(i,j),i=1,ni),j=1,nj)
1522  !write(1,*) ((un(i,j),i=1,ni),j=1,nj)
1523  close(1)
1524 
1525  write(fname,'(a,a)') path(1:kname),'flist.wr'
1526  open(1,file=fname,form='formatted')
1527  !open(1,file='flist.wr',form='formatted')
1528  if(numwr.eq.1) then
1529  write(1,*) str
1530  else
1531  do i=1,numwr-1
1532  read(1,*) dummy
1533  enddo
1534  write(1,*) str
1535  endif
1536  close(1)
1537 
1538 c---------------------------------------------------------------
1539  RETURN
1540  END
1541 c***************************************************************
1542 
1543 c***************************************************************
1544  SUBROUTINE rdd
1545 c
1546  include 'double.inc'
1547  include 'param.inc'
1548  include 'comblc.inc'
1549 
1550  common/comlop/ rxb(nbndp2),zxb(nbndp2),nxb
1551 c
1552  write(fname,'(a,a)') path(1:kname),'out.wr'
1553  open(1,file=fname,form='formatted')
1554  !open(1,file='out.wr')
1555 c open(1,file='out.wr',status='old',form='formatted')
1556 c
1557  read(1,*) ni,nj,ni1,nj1,ni2,nj2,nxb
1558 c
1559  read(1,*) (r(i),i=1,ni)
1560  read(1,*) (z(j),j=1,nj)
1561  read(1,*) ((u(i,j),i=1,ni),j=1,nj)
1562  read(1,*) ((curf(i,j),i=1,ni),j=1,nj)
1563  read(1,*) ((ipr(i,j),i=1,ni),j=1,nj)
1564  read(1,*) rm,zm,um,rx0,zx0,ux0,up
1565  read(1,*) (rxb(ig),ig=1,nxb)
1566  read(1,*) (zxb(ig),ig=1,nxb)
1567 
1568  read(1,*) ((ui(i,j),i=1,ni),j=1,nj)
1569  read(1,*) ((ue(i,j),i=1,ni),j=1,nj)
1570  read(1,*) ((un(i,j),i=1,ni),j=1,nj)
1571  close(1)
1572 c
1573  RETURN
1574  END
1575 c**********************************************************
1576  SUBROUTINE wrdbnd
1577 
1578  include 'double.inc'
1579  include 'param.inc'
1580  include 'comblc.inc'
1581 
1582  write(fname,'(a,a)') path(1:kname),'bnd.wr'
1583  open(1,file=fname,form='formatted')
1584  !open(1,file='bnd.wr',form='formatted')
1585 
1586  write(1,*) ((binadg(i,j),i=1,nbnd),j=1,nbnd)
1587 
1588  close(1)
1589 
1590  RETURN
1591  END
1592 c**********************************************************
1593  SUBROUTINE wrdbin(nk)
1594 
1595  include 'double.inc'
1596  include 'param.inc'
1597  include 'comblc.inc'
1598  integer nk
1599 
1600  write(fname,'(a,a)') path(1:kname),'bndin.wr'
1601  open(1,file=fname,form='formatted')
1602  !open(1,file='bndin.wr',form='formatted')
1603 
1604  write(1,*) nkin,nkout
1605  write(1,*) ((pinadg(ik,j),ik=1,nkout),j=1,nbnd)
1606  write(1,*) (itok(ik),ik=1,nk)
1607  write(1,*) (jtok(ik),ik=1,nk)
1608 
1609  close(1)
1610 
1611 c write(6,*)' nkin,nkout',nkin,nkout
1612 c
1613 c write(6,*) (itok(ik),ik=1,nk)
1614 c write(6,*) (jtok(ik),ik=1,nk)
1615 
1616  RETURN
1617  END
1618 c***************************************************************
1619  SUBROUTINE rddbnd
1620 
1621  include 'double.inc'
1622  include 'param.inc'
1623  include 'comblc.inc'
1624 
1625  write(fname,'(a,a)') path(1:kname),'bnd.wr'
1626  open(1,file=fname,form='formatted')
1627  !open(1,file='bnd.wr',form='formatted')
1628 
1629  read(1,*) ((binadg(i,j),i=1,nbnd),j=1,nbnd)
1630 
1631  close(1)
1632 
1633 
1634  RETURN
1635  END
1636 c**********************************************************
1637 c***************************************************************
1638  SUBROUTINE rddbin(nk)
1639 
1640  include 'double.inc'
1641  include 'param.inc'
1642  include 'comblc.inc'
1643  integer nk
1644 
1645  write(fname,'(a,a)') path(1:kname),'bndin.wr'
1646  open(1,file=fname,form='formatted')
1647  !open(1,file='bndin.wr',form='formatted')
1648 
1649  read(1,*) nkin,nkout
1650  read(1,*) ((pinadg(ik,j),ik=1,nkout),j=1,nbnd)
1651  read(1,*) (itok(ik),ik=1,nk)
1652  read(1,*) (jtok(ik),ik=1,nk)
1653 
1654  close(1)
1655 
1656 c write(6,*)' nkin,nkout',nkin,nkout
1657 c
1658 c write(6,*) (itok(ik),ik=1,nk)
1659 c write(6,*) (jtok(ik),ik=1,nk)
1660 
1661  RETURN
1662  END
1663 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1664 
1665  subroutine eqdsk_rebild
1666 
1667  include 'double.inc'
1668  include 'param.inc'
1669  include 'comblc.inc'
1670 
1671  parameter(np=1000,nbp=np*4)
1672 
1673  character*8 case(6)
1674  dimension rbbbs(nbp),zbbbs(nbp)
1675 
1676  common/efites/ fcefit,rcentr,iefit
1677 
1678  common/com_eqd/ fpol(np),pres(np),qpsi(np),
1679  * ffprim(np),pprime(np),
1680  * rlimtr(np),zlimtr(np),case,simag,sibry,
1681  * idum,nw,nh,limitr
1682 
1683  common/comlop/ rxb(nbndp2), zxb(nbndp2), nxb
1684 
1685 
1686  character*40 eqdfn
1687 c--------------------------------------------------------------------
1688 
1689  write(*,*) '************************* '
1690  write(*,*) ' Entry of subr."eqdsk_rebild":'
1691  write(*,*) '------------------------- '
1692 
1693 
1694  if(ni.ne.nw .or. nj.ne.nh) then
1695 
1696  write(*,*) .ne..or..ne.' error:ninw njnh '
1697  write(*,*) 'nw,nh',nw,nh
1698  write(*,*) 'ni,nj',ni,nj
1699 
1700  endif
1701 
1702  rdim = rmax-rmin
1703  zdim = zmax-zmin
1704 
1705  rleft=rmin
1706  zmid = (zmax+zmin)*0.5d0
1707 
1708  zmaxis=zm
1709  rmaxis=rm
1710 
1711  current=tok*1.d6
1712  !simag=um
1713  !sibry=up
1714 
1715  rcentr=r0ax
1716  bcentr=b0ax
1717 
1718  f2_bon=fpol(nw)**2
1719  pres_bon=pres(nw)
1720  correc=(um-up)/(simag-sibry)
1721 
1722  do i=1,nw
1723 
1724  pres(i)=(pres(i)-pres_bon)*correc + pres_bon
1725 
1726  enddo
1727 
1728  do i=1,nw
1729 
1730  fpol(i)= dsqrt( (fpol(i)**2-f2_bon)*correc + f2_bon )
1731 
1732  enddo
1733 
1734  call get q_spline(qpsi,nw)
1735 
1736  open(1,file='eqdsk_rebild.wr')
1737 
1738  write(1,2000) (case(i),i=1,6),idum,nw,nh
1739  write(*,*) idum,nw,nh
1740 
1741  write(1,2020) rdim,zdim,rcentr,rleft,zmid
1742  write(1,2020) rmaxis,zmaxis,um,up,bcentr
1743  write(1,2020) current,um,xdum,rmaxis,xdum
1744  write(1,2020) zmaxis,xdum,up,xdum,xdum
1745  write(1,2020) (fpol(i),i=1,nw)
1746  write(1,2020) (pres(i),i=1,nw)
1747  write(1,2020) (ffprim(i),i=1,nw)
1748  write(1,2020) (pprime(i),i=1,nw)
1749  write(1,2020) ((u(i,j),i=1,ni),j=1,nj)
1750  write(1,2020) (qpsi(i),i=1,nw)
1751  write(1,2022) nxb,limitr
1752  write(1,2020) (rxb(i),zxb(i),i=1,nxb)
1753  write(1,2020) (rlimtr(i),zlimtr(i),i=1,limitr)
1754 
1755  close(1)
1756 
1757  2000 format(6a8,3i4)
1758  2020 format(5e16.9)
1759  2022 format(2i5)
1760 
1761 
1762  return
1763  end
1764 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1765 
1766 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1767 
1768  subroutine eqdsk_build(eqdfn)
1769 
1770  include 'double.inc'
1771  include 'param.inc'
1772  include 'comblc.inc'
1773 
1774  parameter(np=1000,nbp=np*4)
1775 
1776  character*8 case(6)
1777  dimension rbbbs(nbp),zbbbs(nbp)
1778  dimension ps(np)
1779 
1780  common/efites/ fcefit,rcentr,iefit
1781 
1782  common/com_eqd/ fpol(np),pres(np),qpsi(np),
1783  * ffprim(np),pprime(np),
1784  * rlimtr(np),zlimtr(np),case,simag,sibry,
1785  * idum,nw,nh,limitr
1786 
1787  common/comlop/ rxb(nbndp2), zxb(nbndp2), nxb
1788 
1789  character*40 eqdfn
1790 c--------------------------------------------------------------------
1791 
1792  write(*,*) '************************* '
1793  write(*,*) ' Entry of subr."eqdsk_build":'
1794  write(*,*) '------------------------- '
1795 
1796  case(1) = ' KIAM'
1797  case(2) = ' SPIDER'
1798  idum = 3
1799 
1800  nw=ni
1801  nh=nj
1802 
1803  rdim = rmax-rmin
1804  zdim = zmax-zmin
1805 
1806  rleft=rmin
1807  zmid = (zmax+zmin)*0.5d0
1808 
1809  zmaxis=zm
1810  rmaxis=rm
1811 
1812  current=tok*1.d6
1813  !simag=um
1814  !sibry=up
1815 
1816  rcentr=r0ax
1817  bcentr=b0ax
1818 
1819  f_vcm=bcentr*rcentr
1820 
1821  do i=1,nw
1822  ps(i)= dfloat(i-1)/dfloat(nw-1)
1823  pprime(i)=tabp(1.d0-ps(i))*1.d6/amu0
1824  ffprim(i)=tabf(1.d0-ps(i))
1825  !write(1,*) ps(i),pprime(i)*amu0*1.d-6,ffprim(i)
1826  enddo
1827 
1828  fpol(nw)=f_vcm
1829  pres(nw)=0.d0
1830  dpsi=(um-up)/dfloat(nw-1)
1831 
1832  do i=nw-1,1,-1
1833  pprim_c = 0.5d0*(pprime(i+1)+pprime(i))
1834  fprim_c = 0.5d0*(ffprim(i+1)+ffprim(i))
1835  pres(i)=pres(i+1) + pprim_c*dpsi
1836  fpol(i)=dsqrt(fpol(i+1)**2 + 2.d0*fprim_c*dpsi)
1837  enddo
1838 
1839  call inter_q(ps,nw,qpsi)
1840 
1841  write(fname,'(a,a40)') path(1:kname),eqdfn
1842  open(1,file=fname,form='formatted')
1843  !open(1,file='eqdsk_rebild.wr')
1844 
1845  write(1,2000) (case(i),i=1,6),idum,nw,nh
1846  write(*,*) idum,nw,nh
1847 
1848  write(1,2020) rdim,zdim,rcentr,rleft,zmid
1849  write(1,2020) rmaxis,zmaxis,um,up,bcentr
1850  write(1,2020) current,up,xdum,rmaxis,xdum
1851  write(1,2020) zmaxis,xdum,up,xdum,xdum
1852  write(1,2020) (fpol(i),i=1,nw)
1853  write(1,2020) (pres(i),i=1,nw)
1854  write(1,2020) (ffprim(i),i=1,nw)
1855  write(1,2020) (pprime(i),i=1,nw)
1856  write(1,2020) ((u(i,j),i=1,ni),j=1,nj)
1857  write(1,2020) (0.5d0*qpsi(i)/pi,i=1,nw)
1858  write(1,2022) nxb,nblm
1859  write(1,2020) (rxb(i),zxb(i),i=1,nxb)
1860  write(1,2020) (rblm(i),zblm(i),i=1,nblm)
1861 
1862  close(1)
1863 
1864  2000 format(6a8,3i4)
1865  2020 format(5e16.9)
1866  2022 format(2i5)
1867 
1868 
1869  return
1870  end
1871 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1872  subroutine get q_spline(qpsi,nw)
1873 
1874  include 'double.inc'
1875  include 'dim.inc'
1876  parameter(nspz=nrp+1,nsp4=nspz+4,nsp6=nsp4*6)
1877  include 'compol.inc'
1878 
1879  dimension qpsi(*)
1880 
1881  real*8 qw(nspz),psiw(nspz)
1882  real*8 rrk(nsp4),cck(nsp4),wrk(nsp6)
1883  real*8 cwk(4)
1884 
1885 
1886 ! q_axis and q_bon definition
1887 
1888  xa0=1.d0
1889  xa1=0.5d0*(psia(1)+psia(2))
1890  xa2=0.5d0*(psia(2)+psia(3))
1891  xa3=0.5d0*(psia(3)+psia(4))
1892 
1893  qx1=q(1)
1894  qx2=q(2)
1895  qx3=q(3)
1896 
1897  call extrp2(xa0,qx0, xa1,xa2,xa3, qx1,qx2,qx3)
1898 
1899  xan=psia(iplas)
1900  xa1=0.5d0*(psia(iplas-1)+psia(iplas))
1901  xa2=0.5d0*(psia(iplas-2)+psia(iplas-1))
1902  xa3=0.5d0*(psia(iplas-3)+psia(iplas-2))
1903 
1904  qx1=q(iplas-1)
1905  qx2=q(iplas-2)
1906  qx3=q(iplas-3)
1907 
1908  call extrp2(xan,qxn, xa1,xa2,xa3, qx1,qx2,qx3)
1909 
1910  qw(1)=qx0
1911  psiw(1)=0.d0
1912  do i=2,iplas
1913  qw(i)=q(i-1)
1914  psiw(i)=1.d0-0.5d0*(psia(i-1)+psia(i))
1915  enddo
1916  qw(iplas+1)=qxn
1917  psiw(iplas+1)=1.d0
1918 
1919  nspl=iplas+1
1920 
1921  CALL e01baf(nspl,psiw,qw,rrk,cck,
1922  * nspl+4,wrk,6*nspl+16,ifail)
1923 
1924  qpsi(1)=qw(1)/(2.d0*pi)
1925  do i=2,nw-1
1926  zxx=dfloat((i-1))/dfloat((nw-1))
1927  CALL e02bcf(nspl+4,rrk,cck,zxx,0,cwk,ifail)
1928  qpsi(i)=cwk(1)/(2.d0*pi)
1929  enddo
1930  qpsi(nw)=qw(iplas+1)/(2.d0*pi)
1931 
1932  return
1933  end
arcle
subroutine arcle(NBAS, RBAS, ZBAS, ARC, ARCMAX)
Definition: Eq2_m.f:1181
e02bcf
subroutine e02bcf(NCAP7, K, C, X, LEFT, S, IFAIL)
Definition: NAG.f:2761
rddbin
subroutine rddbin(nk)
Definition: Eq2_m.f:1638
wrdbin
subroutine wrdbin(nk)
Definition: Eq2_m.f:1593
drv
subroutine drv(ZIX1, ZIY1, ZIX2, ZIY2)
Definition: ppplib.f:5671
wrdbnd
subroutine wrdbnd
Definition: Eq2_m.f:1576
fun
subroutine fun(X, F)
Definition: Ev2.f:10
output
subroutine output(NGRID, betpol, zli3)
Definition: Eq2_m.f:1
numcel
subroutine numcel(rrk, zzk, icell, jcell)
Definition: _fluxt.f:259
get
subroutine get(qpsi, nw)
Definition: Eq2_m.f:1872
rddbnd
subroutine rddbnd
Definition: Eq2_m.f:1619
e01baf
subroutine e01baf(M, X, Y, K, C, LCK, WRK, LWRK, IFAIL)
Definition: NAG.f:2511
gap
subroutine gap(ut, rt, zt, rg, zg)
Definition: Eq2_m.f:624
extrp2
subroutine extrp2(X, F, X1, X2, X3, F1, F2, F3)
Definition: com_sub.f:719
wrd
subroutine wrd
Definition: Eq2_m.f:1378
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
rpoint
subroutine rpoint(ut, rt, zt, dp)
Definition: Eq2_m.f:369
agrf
subroutine agrf(N, IA, JA)
Definition: Eq2_m.f:1360
adapol
subroutine adapol
Definition: Eq2_m.f:969
rdd
subroutine rdd
Definition: Eq2_m.f:1544
aprt
subroutine aprt(N, A, IA, JA)
Definition: Eq2_m.f:1327
f_gap
subroutine f_gap(ut, rt, zt, rg, zg)
Definition: Eq2_m.f:798
loop
subroutine loop
Definition: Eq2_m.f:231
current
subroutine current(GEOMETRY, PROFILES, TRANSPORT, SOURCES, EVOLUTION, CONTROL, j_boun, ifail, failstring)
CURRENT TRANSPORT EQUATION.
Definition: main_plasma.f90:205
bttor
subroutine bttor(bettor)
Definition: Eq2_m.f:1141
eqdsk_rebild
subroutine eqdsk_rebild
Definition: Eq2_m.f:1665
wrdfmv
subroutine wrdfmv(numwr, time)
Definition: Eq2_m.f:1480
sarvec
subroutine sarvec(A, JA, IA, N, X, Y)
Definition: Eq2_m.f:1294
zpoint
subroutine zpoint(ut, rt, zt, dp)
Definition: Eq2_m.f:496
eqdsk_build
subroutine eqdsk_build(eqdfn)
Definition: Eq2_m.f:1768
rarvec
subroutine rarvec(A, JA, IA, N, X, Y)
Definition: Eq2_m.f:1277
magpol
subroutine magpol
Definition: Eq2_m.f:1093
wrrec
subroutine wrrec
Definition: Eq2_m.f:1453
btpol
subroutine btpol(betpol)
Definition: Eq2_m.f:1118
tcontr
subroutine tcontr(KAS, A, JA, IA, NN, Z, B, IPRT, IGRF, BZ, ERR)
Definition: Eq2_m.f:1203
inter_q
subroutine inter_q(ps, nw, q_ps)
Definition: com_sub.f:1516
matrix
subroutine matrix
Definition: Matrix.f:2