ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
RIG.f
Go to the documentation of this file.
1  subroutine bound
2 
3  include 'double.inc'
4  include 'param.inc'
5  include 'comblc.inc'
6 
7  real*8 psib(nbndp)
8 
9 
10 c--- array dgdn initialization(calculation dg/dn)
11 
12  ib=1
13 
14  dgdn(ib)=0.d0
15  g(1,1)=0.d0
16 
17  do 10 i=2,ni1
18 
19  ib=ib+1
20  dgdn(ib) = g(i,2)/dz(1)/r(i)
21  g(i,1)=0.d0
22 
23  10 continue
24 
25  ib=ib+1
26  dgdn(ib)=0.d0
27  g(ni,1)=0.d0
28 
29  do 20 j=2,nj1
30 
31  ib=ib+1
32  dgdn(ib) = g(ni1,j)/( dr(ni1)*r12(ni1))
33  g(ni,j)=0.d0
34 
35  20 continue
36 
37  ib=ib+1
38  dgdn(ib)=0.d0
39  g(ni,nj)=0.d0
40 
41  do 30 i=ni1,2,-1
42 
43  ib=ib+1
44  dgdn(ib) = g(i,nj1)/dz(nj1)/r(i)
45  g(i,nj)=0.d0
46 
47  30 continue
48 
49  ib=ib+1
50  dgdn(ib)=0.d0
51  g(1,nj)=0.d0
52 
53  do 40 j=nj1,2,-1
54 
55  ib=ib+1
56  dgdn(ib) = g(2,j)/(dr(1)*r12(1))
57  g(1,j)=0.d0
58 
59  40 continue
60 
61  ib=ib+1
62  dgdn(ib)=0.d0
63 
64  do 100 ib=1,nbnd
65  dgdn(ib)= (dgdn(ib)+dgdn(ib+1))*0.5
66  100 continue
67 
68  do 110 ib=1,nbnd
69  zpsi=0.d0
70  do 111 ibc=1,nbnd
71 
72  zpsi= zpsi+dgdn(ibc)*binadg(ibc,ib)
73 
74  111 continue
75 
76  psib(ib)=zpsi
77  110 continue
78 
79 c...boundary condition for ui
80 
81  ib=0
82 
83  do 200 i=1,ni
84 
85  ib=ib+1
86  ui(i,1)=psib(ib)
87 
88  200 continue
89 
90 
91  do 300 j=2,nj
92 
93  ib=ib+1
94  ui(ni,j)=psib(ib)
95 
96  300 continue
97 
98  do 210 i=ni1,1,-1
99 
100  ib=ib+1
101  ui(i,nj)=psib(ib)
102 
103  210 continue
104 
105  do 310 j=nj1,2,-1
106 
107  ib=ib+1
108  ui(1,j)=psib(ib)
109 
110  310 continue
111 
112  return
113  end
114 
115 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
116 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
117 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
118 
119  subroutine right0(ill,jll,icelm,jcelm,ngav1)
120 
121  parameter(nshp=10)
122 
123  include 'double.inc'
124  include 'param.inc'
125  include 'comblc.inc'
126 
127  common/comomg/ omg,sigm
128  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5),dp1(5),dp2(5),clsq(6)
129  common/comkjf/ curn,nctrli
130  real*8 curn(nip,njp)
131 
132  common /comhel/ helinp, helout
133 
134  common
135  * /c_kpr/kpr
136 
137 
138  real*4 a_print(200)
139  character*30 apr
140 
141 c---------------------------------------------------------------
142  !!amu0=0.4d0*pi
143 
144  nctrli=nctrl
145 
146 ccc--- right hand side for eq. lg=jf in rectangular box.
147 
148  do il=1,neqp
149  right(il)=0.d0
150  enddo
151 c---definition of umax- maximum poloidal flux function value on
152 c grid in plasma.
153 
154  umax=u(imax,jmax)
155 
156  if(iter.ge.iterbf) then
157  imax=ill
158  jmax=jll
159  uem=bline(icelm,jcelm,rl,zl)
160  go to 2000
161  endif
162 
163  1000 continue
164 
165  ic=imax
166  jc=jmax
167 
168  do 100 k=-1,1
169 
170  imp= imax+k
171 
172  do 110 l=-1,1
173 
174  jmp= jmax+l
175 
176  if( u(imp,jmp).gt.umax) then
177 
178  ic=imp
179  jc=jmp
180  umax=u(imp,jmp)
181 
182  endif
183 
184  110 continue
185  100 continue
186 
187  if(ic.ne.imax .OR. jc.ne.jmax) then
188 
189  imax=ic
190  jmax=jc
191  umax=u(imax,jmax)
192  go to 1000
193 
194  endif
195 
196  2000 continue
197 
198  if(kpr.eq.1)write(*,*) 'imax jmax umax',imax,jmax,umax
199  if(kpr.eq.1)write(*,*) 'iclm jclm ',icelm,jcelm
200 
201 c---definition of um - poloidal flux function value at magn. axis.
202 
203  nsh=1
204  xs(nsh)=r(imax)
205  ys(nsh)=z(jmax)
206  fun(nsh)=u(imax,jmax)
207 
208  do 200 k=-1,1
209 
210  i= imax+k
211 
212  do 210 l=-1,1
213 
214  j= jmax+l
215 
216  if(i.eq.imax .AND. j.eq.jmax) go to 210
217  nsh=nsh+1
218  xs(nsh)=r(i)
219  ys(nsh)=z(j)
220  fun(nsh)=u(i,j)
221 
222  210 continue
223  200 continue
224 
225 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
226  call lsq_sur6(xs,ys,fun,nsh,clsq,rrmm,zzmm,uumm,dp)
227  rmold=rm
228  zmold=zm
229  rm=rrmm
230  zm=zzmm
231  um=uumm
232  errm=dsqrt( (rm-rmold)**2+(zm-zmold)**2 )
233 
234 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
235 ! call deriv5(xs,ys,fun,nsh,5,dp)
236 !
237 ! DET = dp(3)*dp(5) - dp(4)**2
238 ! rmold=rm
239 ! zmold=zm
240 ! Rm = Xs(1) + ( dp(2)*dp(4) - dp(1)*dp(5) )/DET
241 ! Zm = Ys(1) + ( dp(1)*dp(4) - dp(2)*dp(3) )/DET
242 !
243 ! errm=dsqrt( (rm-rmold)**2+(zm-zmold)**2 )
244 !
245 ! rmi=r(imax)
246 ! zmj=z(jmax)
247 !
248 ! um=fun(1)+ dp(1)*(rm-rmi) + dp(2)*(zm-zmj)
249 ! + + 0.5d0*dp(3)*(rm-rmi)*(rm-rmi)
250 ! + + dp(4)*(rm-rmi)*(zm-zmj)
251 ! + + 0.5d0*dp(5)*(zm-zmj)*(zm-zmj)
252 !
253 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
254  !write(6,*) 'rm zm um',rm,zm,um
255 
256  if(iter.lt.iterbf) then
257  clr=0.d0
258  clz=0.d0
259  rl=rm
260  zl=zm
261 
262  call shab(ilma,jlma,icelma,jcelma)
263  uem=bline(icelma,jcelma,rl,zl)
264 
265  go to 1010
266  endif
267 c...stabilization by artifical field
268 
269 
270 ! delcr= -( dp(1) + dp(3)*(rl-rmi) + dp(4)*(zl-zmj) )*0.5d0/rl
271 ! delcz= -( dp(2) + dp(5)*(zl-zmj) + dp(4)*(rl-rmi) )
272  delcr= -( dp(3)*(rl-rm) + dp(4)*(zl-zm) )*0.5d0/rl
273  delcz= -( dp(5)*(zl-zm) + dp(4)*(rl-rm) )
274  clr=clr+delcr
275  clz=clz+delcz
276 
277 
278  do 250 i=1,ni
279  do 250 j=1,nj
280 
281  u(i,j)=u(i,j)+delcz*z(j)+delcr*r(i)*r(i)
282 
283  250 continue
284 
285  delc = dabs(delcr)+dabs(delcz)
286 
287  if(delc.lt.1.d-9) go to 1010
288 
289  go to 2000
290 
291  1010 continue
292 
293  if(kpr.eq.1)write(*,*) '*****rm,zm,rl,zl',rm,zm,rl,zl
294  !write(6,*) '*****clr,delcr ',clr,delcr
295  !write(6,*) '*****clz,delcz ',clz,delcz
296  if(kpr.eq.1)write(*,*) '*****clr,clz ',clr,clz
297 
298 
299 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
300 c---definition ux- poloidal flux function value on separatris.!
301 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
302 
303  uxold=ux0
304 
305  call xpoint(ux1,rx1,zx1,ix1,jx1,dp1,1,kodex1)
306  if(kodex1.gt.0) then
307  rx1=rx10
308  zx1=zx10
309  endif
310  if(kpr.eq.1)write(*,*) 'xpoint1',kodex1
311  call xpoint(ux2,rx2,zx2,ix2,jx2,dp2,2,kodex2)
312  if(kodex2.gt.0) then
313  rx2=rx20
314  zx2=zx20
315  endif
316  if(kpr.eq.1)write(*,*) 'xpoint2',kodex2
317 
318  kodxp=kodex1*kodex2
319  if(kodxp.ne.0) then
320 
321  if(kpr.eq.1)write(*,*) ' all x-points out of box'
322  if(kpr.eq.1)write(*,*) ' only limiter case '
323 
324  nctrli=1
325 
326  !write(6,*) ' program cannot run more '
327  !stop
328 
329  endif
330 
331  if(ux1.gt.ux2) then
332 
333  if(kpr.eq.1)write(*,*) 'xpoint=1'
334  ux0= ux1
335 
336  rx0= rx1
337  zx0= zx1
338 
339  ix = ix1
340  jx = jx1
341 
342  rix = r(ix)
343  zjx = z(jx)
344 
345  uix = u(ix,jx)
346 
347  do 1840 is=1,5
348  1840 dp(is)=dp1(is)
349 
350  else
351 
352  if(kpr.eq.1)write(*,*) 'xpoint=2'
353  ux0= ux2
354 
355  rx0= rx2
356  zx0= zx2
357 
358  ix = ix2
359  jx = jx2
360 
361  rix = r(ix)
362  zjx = z(jx)
363 
364  uix = u(ix,jx)
365 
366  do 1940 is=1,5
367  1940 dp(is)=dp2(is)
368 
369  endif
370 
371  if(kpr.eq.1)write(*,*)'rx0,zx0,ux0',rx0,zx0,ux0
372 
373  if(uxold.ge.ux0) ux0=(1.d0-sigm)*uxold+sigm*ux0
374 
375  zvmax=dmax1(zx1,zx2)
376  zvmin=dmin1(zx1,zx2)
377 
378  zver=dabs(zm-zx0)
379 
380 c---definition up- poloidal flux function value on plasma boundary.
381 
382  if(kodxp.eq.0) then
383  ups=um - alp*(um-ux0)
384  else !!!!<- no x-points
385  ups=-1.d12
386  endif
387 
388  if(nctrli.eq.0 ) then
389 
390  up=um - alp*(um-ux0)
391  alpnew=alp
392  numlim=0
393 !!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
394  elseif(nctrli.lt.0 ) then
395 
396  upp=um - alp*(um-ux0)
397  if(psi_bon.lt.ux0) then
398  up=upp
399  alpnew=alp
400  else
401  up=psi_bon
402  alpnew=(um-up)/(um-ux0)
403  endif
404  numlim=0
405 !!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
406 
407  else !!!if(nctrli.gt.0 ) <-limiter case
408 
409  ulmax1=ups
410  ulmax2=ups
411  nulim1=0
412  nulim2=0
413  neigb=0
414 
415  do llim=1,nblm
416 
417  if(kodex1.eq.0) then
418  rlim1=frlim(dp1,zblm(llim),rx1,zx1,rm,zm)
419  rmlim1=frlim(dp1,zm,rx1,zx1,rm,zm)
420  else
421  rlim1=rmax
422  rmlim1=rmax
423  endif
424 
425  if(kodex2.eq.0) then
426  rlim2=frlim(dp2,zblm(llim),rx2,zx2,rm,zm)
427  rmlim2=frlim(dp2,zm,rx2,zx2,rm,zm)
428  else
429  rlim2=rmax
430  rmlim2=rmax
431  endif
432 
433 
434  if( (rblm(llim)-rlim1)*(rm-rmlim1).le.0.d0) then
435 
436  go to 1246
437 
438  elseif( (rblm(llim)-rlim2)*(rm-rmlim2).le.0.d0) then
439 
440  go to 1246
441 
442  endif
443 
444 
445  ilm=iblm(llim)
446  jlm=jblm(llim)
447  ublm=blint(iblm(llim),jblm(llim),rblm(llim),zblm(llim))
448  iprlm=ipr(ilm,jlm)+ipr(ilm+1,jlm)
449  * +ipr(ilm,jlm+1)+ipr(ilm+1,jlm+1)
450 
451 
452  if(iprlm.gt.0) then
453 
454  if(ublm.gt.ulmax1)then
455  ulmax1=ublm
456  nulim1=llim
457  neigb=1
458  endif
459 
460  else
461 
462  if(ublm.gt.ulmax2)then
463  ulmax2=ublm
464  nulim2=llim
465  endif
466 
467  endif
468 
469  1246 continue
470  enddo
471 
472  if(neigb.eq.1) then
473  ublmax=ulmax1
474  numlim=nulim1
475  elseif(neigb.eq.0) then
476  ublmax=ulmax2
477  numlim=nulim2
478  endif
479 
480 
481  up=dmax1(ups,ublmax)
482 
483  if(kodxp.eq.0) alpnew=(um-up)/(um-ux0)
484 
485  endif
486 
487  if(kpr.eq.1) then
488  write(*,*) 'up,numlim',up,numlim,neigb
489  write(*,*) 'psi_bon',psi_bon
490  write(*,*) 'um-up',um-up
491  write(*,*) 'alpnew',alpnew
492  endif
493 
494 
495 
496 c--- un- normal poloidal flux.
497 
498  erru=0.d0
499 
500  do 600 i=1,ni
501  do 600 j=1,nj
502 
503  unold=un(i,j)
504  un(i,j)=(u(i,j)-up)/(um-up)
505  delun=dabs(un(i,j)-unold)
506  erru=dmax1(delun,erru)
507  600 continue
508  if(kpr.eq.1)write(*,*)'erru',erru
509 
510 c--- dimension ipr initialization.
511 c--- if ipr(i,j)=1,then point(i,j) is in plasma.
512 
513  do 700 i=1,ni
514  do 700 j=1,nj
515 
516  ipr(i,j)=0
517 
518  700 continue
519 
520  ipr(imax-1,jmax-1)=1
521 
522  do 800 j=jmax,nj1
523 cccc+ zlim=zvmax
524 cccc+ if(z(j).gt.zlim) go to 800
525 
526  if(kodex1.eq.0) then
527  rlim1=frlim(dp1,z(j),rx1,zx1,rm,zm)
528  rmlim1=frlim(dp1,zm,rx1,zx1,rm,zm)
529  else
530  rlim1=rmax
531  rmlim1=rmax
532  endif
533 
534  if(kodex2.eq.0) then
535  rlim2=frlim(dp2,z(j),rx2,zx2,rm,zm)
536  rmlim2=frlim(dp2,zm,rx2,zx2,rm,zm)
537  else
538  rlim2=rmax
539  rmlim2=rmax
540  endif
541 
542  do 801 i=imax,ni1
543 
544  if(un(i,j).lt.0.d0) then
545 
546  ipr(i,j)=0
547 
548  elseif( (r(i)-rlim1)*(rm-rmlim1).le.0.d0) then
549 
550  ipr(i,j)=0
551 
552  elseif( (r(i)-rlim2)*(rm-rmlim2).le.0.d0) then
553 
554  ipr(i,j)=0
555 
556  else
557 
558  isum=ipr(i-1,j-1)+ipr(i,j-1)+ipr(i+1,j-1)+ipr(i-1,j)
559  if(isum.gt.0) ipr(i,j)=1
560 
561  endif
562 
563  801 continue
564  800 continue
565 
566  do 810 j=jmax,nj1
567 ccc+ zlim=zvmax
568 ccc+ if(z(j).gt.zlim) go to 810
569 
570  if(kodex1.eq.0) then
571  rlim1=frlim(dp1,z(j),rx1,zx1,rm,zm)
572  rmlim1=frlim(dp1,zm,rx1,zx1,rm,zm)
573  else
574  rlim1=rmin
575  rmlim1=rmin
576  endif
577 
578  if(kodex2.eq.0) then
579  rlim2=frlim(dp2,z(j),rx2,zx2,rm,zm)
580  rmlim2=frlim(dp2,zm,rx2,zx2,rm,zm)
581  else
582  rlim2=rmin
583  rmlim2=rmin
584  endif
585 
586  do 811 i=imax,2,-1
587 
588  if(un(i,j).lt.0.d0) then
589 
590  ipr(i,j)=0
591 
592  elseif( (r(i)-rlim1)*(rm-rmlim1).le.0.d0) then
593 
594  ipr(i,j)=0
595 
596  elseif( (r(i)-rlim2)*(rm-rmlim2).le.0.d0) then
597 
598  ipr(i,j)=0
599 
600  else
601 
602  isum=ipr(i-1,j-1)+ipr(i,j-1)+ipr(i+1,j-1)+ipr(i+1,j)
603  if(isum.gt.0.) ipr(i,j)=1
604 
605  endif
606 
607  811 continue
608  810 continue
609 
610  do 820 j=jmax,2,-1
611 ccc+ zlim=zvmin
612 ccc+ if(z(j).lt.zlim) go to 820
613 
614  if(kodex1.eq.0) then
615  rlim1=frlim(dp1,z(j),rx1,zx1,rm,zm)
616  rmlim1=frlim(dp1,zm,rx1,zx1,rm,zm)
617  else
618  rlim1=rmax
619  rmlim1=rmax
620  endif
621 
622  if(kodex2.eq.0) then
623  rlim2=frlim(dp2,z(j),rx2,zx2,rm,zm)
624  rmlim2=frlim(dp2,zm,rx2,zx2,rm,zm)
625  else
626  rlim2=rmax
627  rmlim2=rmax
628  endif
629 
630  do 821 i=imax,ni1
631 
632  if(un(i,j).lt.0.) then
633 
634  ipr(i,j)=0
635 
636  elseif( (r(i)-rlim1)*(rm-rmlim1).le.0.d0) then
637 
638  ipr(i,j)=0
639 
640  elseif( (r(i)-rlim2)*(rm-rmlim2).le.0.d0) then
641 
642  ipr(i,j)=0
643 
644  else
645 
646  isum=ipr(i-1,j+1)+ipr(i,j+1)+ipr(i+1,j+1)+ipr(i-1,j)
647  if(isum.gt.0.) ipr(i,j)=1
648 
649  endif
650 
651  821 continue
652  820 continue
653 
654  do 830 j=jmax,2,-1
655 ccc+ zlim=zvmin
656 ccc+ if(z(j).lt.zlim) go to 830
657 
658  if(kodex1.eq.0) then
659  rlim1=frlim(dp1,z(j),rx1,zx1,rm,zm)
660  rmlim1=frlim(dp1,zm,rx1,zx1,rm,zm)
661  else
662  rlim1=rmin
663  rmlim1=rmin
664  endif
665 
666  if(kodex2.eq.0) then
667  rlim2=frlim(dp2,z(j),rx2,zx2,rm,zm)
668  rmlim2=frlim(dp2,zm,rx2,zx2,rm,zm)
669  else
670  rlim2=rmin
671  rmlim2=rmin
672  endif
673 
674  do 831 i=imax,2,-1
675 
676  if(un(i,j).lt.0.d0) then
677 
678  ipr(i,j)=0
679 
680  elseif( (r(i)-rlim1)*(rm-rmlim1).le.0.) then
681 
682  ipr(i,j)=0
683 
684  elseif( (r(i)-rlim2)*(rm-rmlim2).le.0.) then
685 
686  ipr(i,j)=0
687 
688  else
689 
690  isum=ipr(i-1,j+1)+ipr(i,j+1)+ipr(i+1,j+1)+ipr(i+1,j)
691  if(isum.gt.0) ipr(i,j)=1
692 
693  endif
694 
695  831 continue
696  830 continue
697 
698  !if(nnstpp.eq.4) then
699  ! call wrd
700  ! pause 'wrd'
701  !endif
702 
703 c---troidal current density in plasma
704 
705  do 999 i=1,ni
706  do 999 j=1,nj
707 
708  curf(i,j)=0.d0
709 
710  999 continue
711 
712 
713  tokp=0.d0
714 
715  do 900 i=2,ni2
716 
717  x1=r(i)
718  x2=r(i+1)
719  x3=r(i+1)
720  x4=r(i)
721 
722  x12=0.5d0*(x1+x2)
723  x23=0.5d0*(x2+x3)
724  x34=0.5d0*(x3+x4)
725  x14=0.5d0*(x1+x4)
726 
727  x0=0.5d0*(x12+x34)
728 
729  do 901 j=2,nj2
730 
731  cur1=0.d0
732  cur2=0.d0
733  cur3=0.d0
734  cur4=0.d0
735 
736  if(kodex1.ne.0) go to 2080
737 
738  if(i.gt.ix1.OR. i.lt.ix1-1) go to 2080
739  if(j.gt.jx1 .OR. j.lt.jx1-1) go to 2080
740 
741  nloc=21
742  call xdetal(i,j,ix1,jx1,rx1 ,zx1 ,dp1,nloc,
743  * cur1,cur2,cur3,cur4)
744 
745  go to 9010
746 
747  2080 continue
748 
749  if(kodex2.ne.0) go to 2090
750 
751  if(i.gt.ix2 .OR. i.lt.ix2-1) go to 2090
752  if(j.gt.jx2 .OR. j.lt.jx2-1) go to 2090
753 
754  nloc=21
755  call xdetal(i,j,ix2,jx2,rx2 ,zx2 ,dp2,nloc,
756  * cur1,cur2,cur3,cur4)
757 
758  go to 9010
759 
760  2090 continue
761 
762  ipr1=ipr(i,j)
763  ipr2=ipr(i+1,j)
764  ipr3=ipr(i+1,j+1)
765  ipr4=ipr(i,j+1)
766 
767  isum=ipr1+ipr2+ipr3+ipr4
768  if(isum.eq.0) go to 901
769 
770  z1=z(j)
771  z2=z(j)
772  z3=z(j+1)
773  z4=z(j+1)
774 
775  z12=0.5d0*(z1+z2)
776  z23=0.5d0*(z2+z3)
777  z34=0.5d0*(z3+z4)
778  z14=0.5d0*(z1+z4)
779 
780  z0=0.5d0*(z12+z34)
781 
782  u1=un(i,j)
783  u2=un(i+1,j)
784  u3=un(i+1,j+1)
785  u4=un(i,j+1)
786 
787  if(ipr1.eq.0 .AND. u1.gt.0.d0) u1=0.d0
788  if(ipr2.eq.0 .AND. u2.gt.0.d0) u2=0.d0
789  if(ipr3.eq.0 .AND. u3.gt.0.d0) u3=0.d0
790  if(ipr4.eq.0 .AND. u4.gt.0.d0) u4=0.d0
791 
792  u12=0.5d0*(u1+u2)
793  u23=0.5d0*(u2+u3)
794  u34=0.5d0*(u3+u4)
795  u14=0.5d0*(u1+u4)
796 
797  u0=0.5d0*(u12+u34)
798 
799  sij4=0.25d0*dr(i)*dz(j)
800 
801  cur1=0.d0
802  cur2=0.d0
803  cur3=0.d0
804  cur4=0.d0
805 
806  if(isum.eq.4) then
807 
808  cur1=sij4*0.25d0*(funcur(x1,u1)+funcur(x12,u12)
809  + +funcur(x0,u0)+funcur(x14,u14) )
810 
811  cur2=sij4*0.25d0*(funcur(x2,u2)+funcur(x23,u23)
812  + +funcur(x0,u0)+funcur(x12,u12) )
813 
814  cur3=sij4*0.25d0*(funcur(x3,u3)+funcur(x34,u34)
815  + +funcur(x0,u0)+funcur(x23,u23) )
816 
817  cur4=sij4*0.25d0*(funcur(x4,u4)+funcur(x14,u14)
818  + +funcur(x0,u0)+funcur(x34,u34) )
819 
820  !cur1= funcur(x1,u1)*sij4
821  !cur2= funcur(x2,u2)*sij4
822  !cur3= funcur(x3,u3)*sij4
823  !cur4= funcur(x4,u4)*sij4
824 
825  go to 9010
826  endif
827 
828 c...quadr..(i,j)...
829 
830  if(ipr1.eq.1) then
831 
832  cur1=bcur(x1,x12,x0,x14, z1,z12,z0,z14, u1,u12,u0,u14)
833 
834  elseif(u12.gt.0.d0) then
835 
836  cur1=bcur(x12,x0,x14,x1, z12,z0,z14,z1, u12,u0,u14,u1)
837 
838  elseif(u0.gt.0.d0) then
839 
840  cur1=bcur(x0,x14,x1,x12, z0,z14,z1,z12, u0,u14,u1,u12)
841 
842  elseif(u14.gt.0.d0) then
843 
844  cur1=bcur(x14,x1,x12,x0, z14,z1,z12,z0, u14,u1,u12,u0)
845 
846  endif
847 
848 c...quadr..(i+1,j)...
849 
850  if(ipr2.eq.1) then
851 
852  cur2=bcur(x2,x23,x0,x12, z2,z23,z0,z12, u2,u23,u0,u12)
853 
854  elseif(u23.gt.0.d0) then
855 
856  cur2=bcur(x23,x0,x12,x2, z23,z0,z12,z2, u23,u0,u12,u2)
857 
858  elseif(u0.gt.0.d0) then
859 
860  cur2=bcur(x0,x12,x2,x23, z0,z12,z2,z23, u0,u12,u2,u23)
861 
862  elseif(u12.gt.0) then
863 
864  cur2=bcur(x12,x2,x23,x0, z12,z2,z23,z0, u12,u2,u23,u0)
865 
866  endif
867 
868 c...quadr..(i+1,j+1)...
869 
870  if(ipr3.eq.1) then
871 
872  cur3=bcur(x3,x34,x0,x23, z3,z34,z0,z23, u3,u34,u0,u23)
873 
874  elseif(u34.gt.0.d0) then
875 
876  cur3=bcur(x34,x0,x23,x3, z34,z0,z23,z3, u34,u0,u23,u3)
877 
878  elseif(u0.gt.0.d0) then
879 
880  cur3=bcur(x0,x23,x3,x34, z0,z23,z3,z34, u0,u23,u3,u34)
881 
882  elseif(u23.gt.0.d0) then
883 
884  cur3=bcur(x23,x3,x34,x0, z23,z3,z34,z0, u23,u3,u34,u0)
885 
886  endif
887 
888 c...quadr..(i,j+1)...
889 
890  if(ipr4.eq.1) then
891 
892  cur4=bcur(x4,x14,x0,x34, z4,z14,z0,z34, u4,u14,u0,u34)
893 
894  elseif(u14.gt.0.d0) then
895 
896  cur4=bcur(x14,x0,x34,x4, z14,z0,z34,z4, u14,u0,u34,u4)
897 
898  elseif(u0.gt.0.d0) then
899 
900  cur4=bcur(x0,x34,x4,x14, z0,z34,z4,z14, u0,u34,u4,u14)
901 
902  elseif(u34.gt.0.d0) then
903 
904  cur4=bcur(x34,x4,x14,x0, z34,z4,z14,z0, u34,u4,u14,u0)
905 
906  endif
907 
908  9010 continue
909 
910  tokp=tokp+cur1+cur2+cur3+cur4
911 
912  curf(i,j) = curf(i,j) + cur1
913  curf(i+1,j) = curf(i+1,j) + cur2
914  curf(i+1,j+1) = curf(i+1,j+1)+ cur3
915  curf(i,j+1) = curf(i,j+1) + cur4
916 
917  901 continue
918  900 continue
919 c----------------------------------------------------------
920 
921  tokn=tok
922 
923  if(ngav1.eq.2 .OR. ngav1.eq.3) then
924  if(iter.ge.iterbf .OR. nnstpp.ge.1) then
925  uartm=clr*rl*rl+clz*zl
926  tok=tok*(ucen-uem-uartm)/(um-uem-uartm)
927  endif
928  endif
929 
930  wght=0.25d0
931 
932  if(ngav1.eq.4 .OR. ngav1.eq.5) then
933  if(iter.ge.iterbf .OR. nnstpp.ge.1) then
934  tok=tok*helinp/helout
935  tok=tok*wght+tokn*(1.d0-wght)
936  endif
937  endif
938 
939  cnor=amu0*tok/tokp
940  !if(ngav1.eq.10) cnor=1.d0 !!<<<<<<<<
941 
942  !write(6,*) 'right:tok,tokp,cnor'
943  !write(6,*) tok,tokp,cnor
944  !write(6,*)'tokn,tok', tokn, tok
945 c write(6,*)'h0, h ', helinp, helout
946 
947 c----------------------------------------------------------
948 
949 !!!!!!!!!!!!!like Rus!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
950 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
951  i_rus=0
952  if(i_rus.eq.1)then
953 
954  tokp=0.d0
955 
956  do i=2,ni1
957  do j=2,nj1
958 
959  curf(i,j) = 0.d0
960 
961  iprij=ipr(i,j)
962  xi=r(i)
963  psi=un(i,j)
964  sij=dri(i)*dzj(j)
965 
966  if(iprij.eq.1) then
967  curf(i,j) = funcur( xi,psi )*sij
968  tokp=tokp+curf(i,j)
969  endif
970 
971  enddo
972  enddo
973  cnor=amu0*tok/tokp
974 
975  end if
976 
977 
978 !
979 !!!!!!!!!!!!!!!!!!like Rus!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
980 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
981 !
982 ! f_cur=0.d0
983 !
984 ! do i=2,ni1
985 ! do j=2,nj1
986 !
987 ! iprij=ipr(i,j)
988 ! xi=r(i)
989 ! psi=un(i,j)
990 ! sij=dri(i)*dzj(j)
991 !
992 ! if(iprij.eq.1) then
993 ! f_cur=f_cur+funcur_f( xi,psi )*sij
994 ! endif
995 !
996 ! enddo
997 ! enddo
998 !
999 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1000 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1001 
1002 c---right hand side definition
1003 
1004  do 950 i=2,ni1
1005  do 951 j=2,nj1
1006  curf(i,j) = cnor*curf(i,j)
1007  il=nlin(i,j)
1008 
1009  if(iter.eq.1) then
1010  right(il)=-curf(i,j)
1011  else
1012  right(il)=-curf(i,j)*omg-(1.d0-omg)*curn(i,j)
1013  endif
1014 
1015  curn(i,j)=curf(i,j)
1016  curf(i,j)= curf(i,j)/(dri(i)*dzj(j))
1017 
1018  ! write(6,*) 'i,j,ipr(i,j),curf(i,j)'
1019  ! write(6,*) i,j,ipr(i,j),curf(i,j)
1020 
1021  951 continue
1022  950 continue
1023 
1024  return
1025  end
1026 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1027 
1028  subroutine xpoint(ux,rx,zx,ix,jx,dp,numxp,kodex)
1029 
1030  include 'double.inc'
1031  parameter(nshp=10)
1032  include 'param.inc'
1033  include 'comblc.inc'
1034 
1035  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
1036 
1037  ix00=ix
1038  jx00=jx
1039 
1040  numitc=0
1041  kodex=0
1042  444 continue
1043 
1044  numitc=numitc+1
1045 
1046  if(numitc.ge.30) then
1047 c write(6,*) 'not correct definition x-point location'
1048 c write(6,*) 'code=444 (cell,containing x-point)'
1049  return
1050  endif
1051 
1052  if(rx.ge.rmax .OR. rx.le.rmin .OR.
1053  * zx.ge.zmax .OR. zx.le.zmin) then
1054 ccc write(6,*)' xpoint',numxp,'out of box'
1055 ccc write(6,*)' rx,zx ',rx,zx
1056 
1057  ux=-1.d6
1058  kodex=1
1059  return
1060  endif
1061 
1062 c---definition of cell, containig x-point
1063 
1064  do 300 i=1,ni1
1065  icx=i
1066  if( (rx.lt.r(i+1)) .AND. (rx.ge.r(i)) ) go to 301
1067  300 continue
1068  301 continue
1069 
1070  do 310 j=1,nj1
1071  jcx=j
1072  if( (zx.lt.z(j+1)) .AND. (zx.ge.z(j)) ) go to 311
1073  310 continue
1074  311 continue
1075 
1076  if(icx.eq.1 .OR. icx.eq.ni1 .OR.
1077  * jcx.eq.1 .OR. jcx.eq.nj1) then
1078 ccc write(6,*)' xpoint',numxp,'in bound cell'
1079 ccc write(6,*)' icell,jcell ',icx,jcx
1080 
1081  ux=-1.d6
1082  kodex=1
1083 
1084  return
1085  endif
1086 
1087 ccc write(6,*) 'icelx,jcelx',icx,jcx,rx,zx
1088 
1089 c---define nearest knote
1090 
1091  sdmin=rmax
1092 
1093  do 320 k=0,1
1094  rr=r(icx+k)
1095  do 325 l=0,1
1096  zz=z(jcx+l)
1097  dlx=dsqrt( (rr-rx)**2+(zz-zx)**2 )
1098  if(dlx.lt.sdmin) then
1099  sdmin=dlx
1100  ix=icx+k
1101  jx=jcx+l
1102  endif
1103  325 continue
1104  320 continue
1105 
1106  numit=0
1107  555 continue
1108  numit=numit+1
1109 
1110 ccc write(6,*) 'ix,jx',ix,jx
1111 ccc write(6,*) 'ix0,jx0',ix00,jx00
1112 
1113  nsh=1
1114  xs(nsh)=r(ix)
1115  ys(nsh)=z(jx)
1116  fun(nsh)=u(ix,jx)
1117 
1118  do 400 k=-1,1
1119 
1120  i= ix+k
1121 
1122  do 410 l=-1,1
1123 
1124  j= jx+l
1125 
1126  if(k.eq.0 .AND. l.eq.0 ) go to 410
1127  nsh=nsh+1
1128  xs(nsh)=r(i)
1129  ys(nsh)=z(j)
1130  fun(nsh)=u(i,j)
1131 
1132  410 continue
1133  400 continue
1134 
1135  call deriv5(xs,ys,fun,nsh,5,dp)
1136 
1137  det = dp(3)*dp(5) - dp(4)**2
1138 
1139  ! if(det.gt.0.d0) then
1140 ccc ! write(*,*)' xpoint',numxp,'det>0'
1141  ! ux=-1.d6
1142  ! kodex=1
1143  ! return
1144  ! endif
1145 
1146  rx = xs(1) + ( dp(2)*dp(4) - dp(1)*dp(5) )/det
1147  zx = ys(1) + ( dp(1)*dp(4) - dp(2)*dp(3) )/det
1148 
1149  rrx=xs(1)
1150  zzx=ys(1)
1151 
1152  ux=fun(1)+ dp(1)*(rx-rrx) + dp(2)*(zx-zzx)
1153  + + 0.5d0*dp(3)*(rx-rrx)*(rx-rrx)
1154  + + dp(4)*(rx-rrx)*(zx-zzx)
1155  + + 0.5d0*dp(5)*(zx-zzx)*(zx-zzx)
1156 
1157 ccc write(6,*) 'rx zx ux',rx,zx,ux
1158 
1159  if( rx.gt.r(ix+1) .OR. rx.lt.r(ix-1) ) go to 444
1160  if( zx.gt.z(jx+1) .OR. zx.lt.z(jx-1) ) go to 444
1161 
1162  sdmin=rmax
1163 
1164  do 500 k=-1,1
1165  rr=r(ix+k)
1166  do 510 l=-1,1
1167  zz=z(jx+l)
1168  dlx=dsqrt( (rr-rx)**2+(zz-zx)**2 )
1169  if(dlx.lt.sdmin) then
1170  sdmin=dlx
1171  ixp=ix+k
1172  jxp=jx+l
1173  endif
1174  510 continue
1175  500 continue
1176 
1177  if(ixp.ne.ix .OR. jxp.ne.jx) then
1178 
1179  if(numit.ge.20) then
1180 c write(6,*) 'not correct definition x-point location'
1181 c write(6,*) 'code=555 (nearest knote to x-point)'
1182  return
1183  endif
1184 
1185  ix=ixp
1186  jx=jxp
1187 
1188  go to 555
1189 
1190  endif
1191 
1192  if(ix.ne.ix00 .OR. jx.ne.jx00) then
1193 
1194  dold=dsqrt( (rx-r(ix00))**2 + (zx-z(jx00))**2 )
1195  dnew=dsqrt( (rx-r(ix))**2 + (zx-z(jx))**2 )
1196  deld01=dabs(dold-dnew)
1197 
1198  if(deld01.lt.0.001d0*dr(ix)) then
1199  ix=ix00
1200  jx=jx00
1201  go to 555
1202  endif
1203 
1204  endif
1205 
1206  return
1207  end
1208 
1209 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1210 
1211  real*8 function bcur( r1,r2,r3,r4,z1,z2,z3,z4,
1212  + psi1,psi2,psi3,psi4 )
1213  implicit real*8(a-h,o-z)
1214 
1215  squtr(x1,x2,x3,y1,y2,y3)=
1216  = 0.5d0*(y3*(x2-x1)+y1*(x3-x2)+y2*(x1-x3))
1217 
1218  xzer(x1,x2,v1,v2) = (x1*v2-x2*v1)/(v2-v1+1.d-8)
1219 
1220  iflag1=0
1221  iflag2=0
1222 
1223  if(psi2.gt.0.d0 .AND. psi3.gt.0.d0 .AND. psi4.gt.0.d0) then
1224 
1225  !r10=r3
1226  !z10=z3
1227 
1228  !r20=r3
1229  !z20=z3
1230 
1231  !psi10=psi3
1232  !psi20=psi3
1233 
1234  !go to 100
1235 
1236  cur1=funcur(r1,psi1)
1237  cur2=funcur(r2,psi2)
1238  cur3=funcur(r3,psi3)
1239  cur4=funcur(r4,psi4)
1240 
1241  sij=squtr(r1,r2,r3,z1,z2,z3)+squtr(r1,r3,r4,z1,z3,z4)
1242  bcur=0.25d0*sij*(cur1+cur2+cur3+cur4)
1243  return
1244 
1245  endif
1246 
1247 c.....definition of first zero
1248 
1249  if(psi2.le.0.d0) then
1250 
1251  r10=xzer(r1,r2,psi1,psi2)
1252  z10=xzer(z1,z2,psi1,psi2)
1253  psi10=0.d0
1254 
1255  elseif(psi3.le.0.d0) then
1256 
1257  r10=xzer(r2,r3,psi2,psi3)
1258  z10=xzer(z2,z3,psi2,psi3)
1259  psi10=0.
1260 
1261  elseif(psi4.le.0.d0) then
1262 
1263  iflag1=1
1264  r10=xzer(r3,r4,psi3,psi4)
1265  z10=xzer(z3,z4,psi3,psi4)
1266  psi10=0.d0
1267 
1268  endif
1269 
1270 c.....definition of second zero
1271 
1272  if(psi4.le.0.d0) then
1273 
1274  r20=xzer(r1,r4,psi1,psi4)
1275  z20=xzer(z1,z4,psi1,psi4)
1276  psi20=0.d0
1277 
1278  elseif(psi3.le.0.d0) then
1279 
1280  r20=xzer(r4,r3,psi4,psi3)
1281  z20=xzer(z4,z3,psi4,psi3)
1282  psi20=0.d0
1283 
1284  elseif(psi2.le.0.d0) then
1285 
1286  iflag2=1
1287  r20=xzer(r3,r2,psi3,psi2)
1288  z20=xzer(z3,z2,psi3,psi2)
1289  psi20=0.d0
1290 
1291  endif
1292 
1293  100 continue
1294 
1295  cur1=funcur(r1,psi1)
1296  cur2=funcur(r2,psi2)
1297  cur3=funcur(r3,psi3)
1298  cur4=funcur(r4,psi4)
1299  cur10=funcur(r10,psi10)
1300  cur20=funcur(r20,psi20)
1301 
1302  zbcur = squtr(r1,r20,r4,z1,z20,z4)*(cur1+cur20+cur4) +
1303  + squtr(r1,r10,r20,z1,z10,z20)*(cur1+cur10+cur20) +
1304  + squtr(r1,r2,r10,z1,z2,z10)*(cur1+cur2+cur10) +
1305  + iflag1*squtr(r2,r3,r10,z2,z3,z10)*(cur2+cur3+cur10) +
1306  + iflag2*squtr(r3,r4,r20,z3,z4,z20)*(cur3+cur4+cur20)
1307 
1308 c if(psi3.gt.0. .AND. psi2.lt.0. .AND. psi4.lt.0.) then
1309 c
1310 c r30=xzer(r3,r2,psi3,psi2)
1311 c z30=xzer(z3,z2,psi3,psi2)
1312 c
1313 c r40=xzer(r4,r3,psi4,psi3)
1314 c z40=xzer(z4,z3,psi4,psi3)
1315 c
1316 c psi30=0.
1317 c psi40=0.
1318 c
1319 c cur30=funcur(r30,psi30)
1320 c cur40=funcur(r40,psi40)
1321 c
1322 c zbcur = zbcur +
1323 c + squtr(r3,r40,r30,z3,z40,z30)*(cur3+cur40+cur30) +
1324 c + squtr(r10,r30,r20,z10,z30,z20)*(cur10+cur30+cur20) +
1325 c + squtr(r30,r40,r20,z30,z40,z20)*(cur30+cur40+cur20)
1326 c
1327 c endif
1328 
1329  bcur=zbcur/3.d0
1330 !!!tst
1331  !bcur=0.d0 !!!tst!!!!
1332 !!!tst
1333 
1334  return
1335  end
1336 
1337 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1338 
1339  subroutine xdetal(i,j,ix,jx,rxpn ,zxpn ,dp,nloc,
1340  * cur1,cur2,cur3,cur4)
1341 
1342  include 'double.inc'
1343  parameter(nlocp=64)
1344  include 'param.inc'
1345  include 'comblc.inc'
1346 
1347  real*8 dp(5)
1348  real*8 rloc(nlocp),zloc(nlocp)
1349  real*8 uloc(nlocp,nlocp),curloc(nlocp,nlocp)
1350 
1351  psilim=0.05d0
1352 
1353  nlc=(nloc+1)/2
1354 
1355  drloc=(r(i+1)-r(i))/(nloc-1.d0)
1356  dzloc=(z(j+1)-z(j))/(nloc-1.d0)
1357 
1358  rloc(1)=r(i)
1359  zloc(1)=z(j)
1360 
1361  do 10 is=2,nloc
1362  rloc(is)=rloc(is-1)+drloc
1363  zloc(is)=zloc(is-1)+dzloc
1364  10 continue
1365 
1366  uix=u(ix,jx)
1367  rix=r(ix)
1368  zjx=z(jx)
1369 
1370  do 100 is=1,nloc
1371  ris=rloc(is)
1372  do 110 js=1,nloc
1373  zjs=zloc(js)
1374  uisjs=uix +dp(1)*(ris-rix) + dp(2)*(zjs-zjx)
1375  + + 0.5*dp(3)*(ris-rix)*(ris-rix)
1376  + + dp(4)*(ris-rix)*(zjs-zjx)
1377  + + 0.5*dp(5)*(zjs-zjx)*(zjs-zjx)
1378 
1379  uloc(is,js)=(uisjs-up)/(um-up)
1380 
1381  110 continue
1382  100 continue
1383 
1384  do 200 is=1,nloc
1385  ris=rloc(is)
1386  do 210 js=1,nloc
1387  zjs=zloc(js)
1388  psiloc=uloc(is,js)
1389 
1390  if(psiloc.le.0.d0) then
1391  curloc(is,js)=0.d0
1392  go to 210
1393  endif
1394 
1395  rlim=frlim(dp,zjs,rxpn,zxpn,rm,zm)
1396  rmlim=frlim(dp,zm,rxpn,zxpn,rm,zm)
1397  if( (ris -rlim)*(rm-rmlim).le.0.) then
1398 ccc+ if(zjs.gt.zvmax .OR. zjs.lt.zvmin) then
1399  curloc(is,js)=0.d0
1400  go to 210
1401  endif
1402 
1403 c if(psiloc.le.psilim) then
1404 c curlim=funcur(ris,psilim)
1405 c curloc(is,js)=(curlim*psiloc)/psilim
1406 c go to 210
1407 c endif
1408  curloc(is,js)=funcur(ris,psiloc)
1409 
1410  210 continue
1411  200 continue
1412 
1413  sloc=drloc*dzloc
1414 
1415 c..quadr(i,j)
1416 
1417  do 1000 is=1,nlc-1
1418  do 1000 js=1,nlc-1
1419 
1420  cur1=cur1+curloc(is,js)+curloc(is+1,js)+
1421  + curloc(is+1,js+1)+curloc(is,js+1)
1422 
1423  1000 continue
1424 
1425  cur1=cur1*0.25d0*sloc
1426 
1427 c..quadr(i+1,j)
1428 
1429  do 2000 is=nlc,nloc-1
1430  do 2000 js=1,nlc-1
1431 
1432  cur2=cur2+curloc(is,js)+curloc(is+1,js)+
1433  + curloc(is+1,js+1)+curloc(is,js+1)
1434 
1435  2000 continue
1436 
1437  cur2=cur2*0.25d0*sloc
1438 
1439 c..quadr(i+1,j+1)
1440 
1441  do 3000 is=nlc,nloc-1
1442  do 3000 js=nlc,nloc-1
1443 
1444  cur3=cur3+curloc(is,js)+curloc(is+1,js)+
1445  + curloc(is+1,js+1)+curloc(is,js+1)
1446 
1447  3000 continue
1448 
1449  cur3=cur3*0.25d0*sloc
1450 
1451 c..quadr(i,j+1)
1452 
1453  do 4000 is=1,nlc-1
1454  do 4000 js=nlc,nloc-1
1455 
1456  cur4=cur4+curloc(is,js)+curloc(is+1,js)+
1457  + curloc(is+1,js+1)+curloc(is,js+1)
1458 
1459  4000 continue
1460 
1461  cur4=cur4*0.25d0*sloc
1462 
1463  return
1464  end
1465 
1466 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1467 
1468  subroutine right1
1469 
1470  include 'double.inc'
1471  include 'param.inc'
1472  include 'comblc.inc'
1473 
1474 ccc--- right hand side for eq. lu=jf in rectangular box.
1475 
1476  do 100 i=2,ni1
1477 
1478  ri=r(i)
1479 
1480  j=2
1481 
1482  dzm=dz(1)
1483  il=nlin(i,j)
1484 
1485  right(il)=right(il)-ui(i,j-1)*dri(i)/(dzm*ri)
1486 
1487  j=nj1
1488 
1489  dzp=dz(nj1)
1490  il=nlin(i,j)
1491 
1492  right(il)=right(il)-ui(i,j+1)*dri(i)/(dzp*ri)
1493 
1494  100 continue
1495 
1496 
1497  do 200 j=2,nj1
1498 
1499 
1500  i=ni1
1501 
1502  drp=dr(ni1)
1503  rpl=r12(ni1)
1504  il=nlin(i,j)
1505 
1506  right(il)=right(il)-ui(i+1,j)*dzj(j)/(rpl*drp)
1507 
1508  i=2
1509 
1510  drm=dr(1)
1511  rmn=r12(1)
1512  il=nlin(i,j)
1513 
1514  right(il)=right(il)-ui(i-1,j)*dzj(j)/(rmn*drm)
1515 
1516  200 continue
1517 
1518  return
1519  end
1520 
1521 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1522 
1523  real*8 function frlim(dp,ylim,rx,zx,rm,zm)
1524 
1525  implicit real*8(a-h,o-z)
1526  real*8 dp(5)
1527 
1528  dxx=dp(3)
1529  dxy=dp(4)
1530  dyy=dp(5)
1531 
1532  disc=(dxy/dyy)**2 - dxx/dyy
1533 
1534  if(disc.lt.0.) then
1535 
1536  ! write(6,*) ' FRLIM: ***ERROR disc lt.0'
1537  ! write(6,*) ' rx zx ',rx,zx
1538  ! stop
1539 
1540  cc=(zm-zx)/(rm-rx)
1541 
1542  cc=-1.d0/cc
1543  go to 100
1544 
1545  endif
1546 
1547  cdpls = -dxy/dyy + dsqrt(disc)
1548  cdmns = -dxy/dyy - dsqrt(disc)
1549 
1550  ang1 = 0.5d0*(datan(cdpls)+datan(cdmns))
1551  ang2 =-0.5d0*(datan(1.d0/cdpls)+datan(1.d0/cdmns))
1552 
1553 
1554 c.........calculation d2u/dl2(direction ang1 )
1555 
1556  c1=dtan(ang1)
1557  dl2x=dyy*c1*c1+2.d0*dxy*c1+dxx
1558 
1559 c.........calculation d2u/dl2(direction ang2 )
1560 
1561  c2=dtan(ang2)
1562  dl2y=dyy*c2*c2+2.*dxy*c2+dxx
1563 
1564  if(dl2x.lt.0.) then
1565 
1566  cc=c1
1567 
1568  elseif(dl2y.lt.0.) then
1569 
1570  cc=c2
1571 
1572  else
1573 
1574  ! write(6,*) ' FRLIM: ***ERROR both deriv. ge.0'
1575  ! write(6,*) ' rx zx ',rx,zx
1576  ! call wrd
1577  ! stop
1578 
1579  cc=(zm-zx)/(rm-rx)
1580 
1581  cc=-1.d0/cc
1582 
1583  endif
1584 
1585  100 continue
1586 
1587  frlim=rx+(ylim-zx)/cc
1588 
1589  return
1590  end
1591 
1592 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1593 
1594 
nlin
function nlin(i, j)
Definition: Matrix.f:173
bline
real *8 function bline(i, j, r0, z0)
Definition: EQ1_m.f:865
drp
subroutine drp(ZIX, ZIY)
Definition: ppplib.f:5738
right1
subroutine right1
Definition: RIG.f:1468
zero
subroutine zero(x1, y1, x2, y2, func, err, x, y, izero, ll)
Definition: zero.f90:1
xdetal
subroutine xdetal(i, j, ix, jx, rxpn,zxpn,dp, nloc,
Definition: RIG.f:1339
blint
real *8 function blint(i, j, r0, z0)
Definition: EQ1_m.f:896
fun
subroutine fun(X, F)
Definition: Ev2.f:10
bcur
real *8 function bcur(r1, r2, r3, r4, z1, z2, z3, z4,
Definition: RIG.f:1211
bound
subroutine bound
Definition: RIG.f:1
eq
subroutine eq(pcequi, psicon, ncequi, nstep, ngrid,
Definition: Eq_m.f:310
right0
subroutine right0(ill, jll, icelm, jcelm, ngav1)
Definition: RIG.f:119
flux
subroutine flux(psitok, rk, zk, nk)
Definition: EQ1_m.f:786
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
grid
subroutine grid
Definition: EQ1_m.f:926
frlim
real *8 function frlim(dp, ylim, rx, zx, rm, zm)
Definition: RIG.f:1523
axis
subroutine axis(n, h, r, f)
Definition: solution2.f90:586
current
subroutine current(GEOMETRY, PROFILES, TRANSPORT, SOURCES, EVOLUTION, CONTROL, j_boun, ifail, failstring)
CURRENT TRANSPORT EQUATION.
Definition: main_plasma.f90:205
shab
subroutine shab(il, jl, icelm, jcelm)
Definition: Eq_m.f:923
lsq_sur6
subroutine lsq_sur6(r, z, u, n, c, rm, zm, um, dp)
Definition: scu.f:3
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
box
subroutine box(IX1, IX2, IY1, IY2)
Definition: ppplib.f:4360
xpoint
subroutine xpoint(ux, rx, zx, ix, jx, dp, numxp, kodex)
Definition: RIG.f:1028