ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
_extp.f
Go to the documentation of this file.
1  subroutine psiloop(rlop,zlop,psilop,nlop)
2 
3  include 'double.inc'
4  dimension rlop(*),zlop(*),psilop(*)
5 
6  call f_psloop_e(rlop,zlop,psilop,nlop)
7  call psloop_p(rlop,zlop,psilop,nlop)
8 
9  return
10  end
11 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
12 
13  subroutine b_prob(rprob,zprob,hrprob,hzprob,nprob,bprob)
14 
15  include 'double.inc'
16  dimension rprob(*),zprob(*),hrprob(*),hzprob(*),bprob(*)
17 
18  call bprob_e(rprob,zprob,hrprob,hzprob,nprob,bprob)
19  call bprob_p(rprob,zprob,hrprob,hzprob,nprob,bprob)
20 
21  return
22  end
23 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
24 
25  subroutine bprob_p(rprob,zprob,hrprob,hzprob,nprob,bprob)
26 
27  include 'double.inc'
28  parameter(nshp=10)
29  include 'parevo.inc'
30  parameter(nkp=njlim)
31  include 'dim.inc'
32  include 'compol.inc'
33  include 'compol_add.inc'
34 
35  dimension rprob(*),zprob(*),hrprob(*),hzprob(*),bprob(*)
36  dimension rsh(nshp),zsh(nshp),ush(nshp)
37  dimension dp(5)
38 
39  sqrt(arg)=dsqrt(arg)
40 
41  ik_out=0
42 
43  do 100 ik=1,nprob
44 
45  r0=rprob(ik)
46  z0=zprob(ik)
47 
48  call numcel(r0,z0,ic,jc)
49 
50  if(iprprob(ik).eq.1) then !!!
51 
52  if(ic.eq.nr) then
53  write(*,*) 'probe',ik,'is out of box'
54  write(*,*) 'program is terminated'
55  stop
56  endif
57 
58  dismin=rm
59  do ii=0,1
60  is=ic+ii
61  do jj=0,1
62  js=jc+jj
63  dist=sqrt( (r0-r(is,js))**2 + (z0-z(is,js))**2 )
64  if(dist .lt. dismin) then
65  dismin=dist
66  ick=is
67  jck=js
68  endif
69  enddo
70  enddo
71 
72  if(jck.eq.nt) jck=2
73  if(ick.eq.nr) ick=nr-1
74 
75  nsh=1
76  rsh(nsh)=r(ick,jck)
77  zsh(nsh)=z(ick,jck)
78  ush(nsh)=psii(ick,jck)
79  do ii=-1,1
80  is=ick+ii
81  do jj=-1,1
82  js=jck+jj
83  if(ii.ne.0 .or. jj.ne.0) then
84  nsh=nsh+1
85  rsh(nsh)=r(is,js)
86  zsh(nsh)=z(is,js)
87  ush(nsh)=psii(is,js)
88  endif
89  enddo
90  enddo
91 
92  call deriv5(rsh,zsh,ush,nsh,5,dp)
93 
94  dudr=dp(1)+dp(3)*(r0-rsh(1))+dp(4)*(z0-zsh(1))
95  dudz=dp(2)+dp(5)*(z0-zsh(1))+dp(4)*(r0-rsh(1))
96 
97  br=-dudz/r0
98  bz= dudr/r0
99  bpla=br*hrprob(ik)+bz*hzprob(ik)
100  bprob(ik)=bprob(ik)+bpla
101 
102  else !!!
103 
104  if(ic.ne.nr) then
105  write(*,*) 'bprob_p:prob is in the box',ic,nr
106  write(*,*) 'program is terminated'
107  write(*,*) 'you need to consult program provider'
108  stop
109  endif
110 
111  ik_out=ik_out+1
112 
113  dudr=0.d0
114  dudz=0.d0
115 
116  do 110 j=2,nt1
117 
118  dudr = dudr-adginr(ik_out,j)*(dgdn(j)+dgdn(j+1))*0.5d0
119  dudz = dudz-adginz(ik_out,j)*(dgdn(j)+dgdn(j+1))*0.5d0
120 
121  110 continue
122 
123  br=-dudz/r0
124  bz= dudr/r0
125  bpla=br*hrprob(ik)+bz*hzprob(ik)
126  bprob(ik)=bprob(ik)+bpla
127 
128  endif !!!
129 
130  100 continue
131 
132  if(ik_out.ne.nprob_out) then
133  write(*,*) .ne.'bprob_p:ik_outnprob_out',ik_out,nprob_out
134  write(*,*) 'program is terminated'
135  write(*,*) 'you need to consult program provider'
136  stop
137  endif
138 
139  return
140  end
141 
142 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
143 
144  subroutine blic_d(r0,z0, r1,r2,r3,r4, z1,z2,z3,z4,
145  * u1,u2,u3,u4, u0, dudr,dudz)
146  include 'double.inc'
147 
148  !dimension a(4,4),f(4),x(4),iwrk(4)
149 
150  eps=1.d-9
151 !!!!!!!!test!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
152 
153  r1=1.d0
154  r4=1.d0
155  r2=2.d0
156  r3=2.d0
157 
158  z1=1.d0
159  z4=2.d0
160  z2=1.d0
161  z3=2.d0
162 
163  u1=1.d0
164  u4=1.d0
165  u2=0.d0
166  u3=0.d0
167 
168  r0=1.5d0
169  z0=1.5d0
170 
171 !!!!!!!!!test!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
172 
173  ar=r3-r4+r1-r2+eps
174  br=-r1+r2
175  cr=-r1+r4
176  dr= r1
177 
178  az=z3-z4+z1-z2+eps
179  bz=-z1+z2
180  cz=-z1+z4
181  dz= z1
182 
183  au=u3-u4+u1-u2
184  bu=-u1+u2
185  cu=-u1+u4
186  du= u1
187 
188  cx = dr-r0-cr/(-az*cr+ar*cz)*az*r0+cr/(-az*cr+ar*cz)*az*dr+cr/(-az
189  #*cr+ar*cz)*ar*z0-cr/(-az*cr+ar*cz)*ar*dz
190 
191  bx = -ar/(-az*cr+ar*cz)*az*r0+ar/(-az*cr+ar*cz)*az*dr+ar**2/(-az*c
192  #r+ar*cz)*z0-ar**2/(-az*cr+ar*cz)*dz+cr/(-az*cr+ar*cz)*az*br-cr/(-a
193  #z*cr+ar*cz)*ar*bz+br
194 
195  ax = ar/(-az*cr+ar*cz)*az*br-ar**2/(-az*cr+ar*cz)*bz
196 
197  if(abs(ax/bx).lt.1.d-8) then
198  x=-cx/bx
199  go to 10
200  endif
201 
202  det_=bx**2-4.d0*ax*cx
203  if(det_.lt.0.d0) then
204  write(*,*) 'blic_d: det<0'
205  pause 'pause'
206  endif
207  x1 = 0.5d0*(-bx + sqrt(det_))/ax
208  x2 = 0.5d0*(-bx - sqrt(det_))/ax
209  if(x1.ge.0.d0 .And. x1.le.1.d0) then
210  x=x1
211  elseif(x2.ge.0.d0 .And. x2.le.1.d0) then
212  x=x2
213  else
214  write(*,*) 'blic_d: x<0 or x>1'
215  pause 'pause'
216  endif
217 
218 10 continue
219 
220  y = -(az*r0-az*br*x-az*dr-ar*z0+ar*bz*x+ar*dz)/(-az*cr+ar*cz)
221 
222  u0=au*x*y + bu*x + cu*y + du
223 
224  dudx=au*y + bu
225  dudy=au*x + cu
226 
227  drdx=ar*y + br
228  drdy=ar*x + cr
229 
230  dzdx=az*y + bz
231  dzdy=az*x + cz
232 
233  det = drdx*dzdy - drdy*dzdx
234  det_r = dudx*dzdy - dudy*dzdx
235  det_z = drdx*dudy - drdy*dudx
236 
237  dudr=det_r/det
238  dudz=det_z/det
239 
240  return
241  end
242 
243 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
244 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
245 
246  subroutine blic_d_(r0,z0, r1,r2,r3,r4, z1,z2,z3,z4,
247  * u1,u2,u3,u4, u0, dudr,dudz)
248  include 'double.inc'
249 
250  dimension a(4,4),f(4),x(4),iwrk(4)
251 
252 
253 !!!!!!!!test!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
254 !
255 ! r1=1.d0
256 ! r4=1.d0
257 ! r2=2.d0
258 ! r3=2.d0
259 !
260 ! z1=1.d0
261 ! z4=2.d0
262 ! z2=1.d0
263 ! z3=2.d0
264 !
265 ! u1=0.d0
266 ! u4=1.d0
267 ! u2=1.d0
268 ! u3=0.d0
269 !
270 ! r0=1.5d0
271 ! z0=1.5d0
272 !
273 !!!!!!!!!test!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
274 
275 
276 
277 
278 
279 
280  a(1,1)=r1 !a(1,1) a(1,2) a(1,3) a(1,4)! !x(1)! !f(1)!
281  a(2,1)=r2 !a(2,1) a(2,2) a(2,3) a(2,4)! !x(2)! !f(2)!
282  a(3,1)=r3 !a(3,1) a(3,2) a(3,3) a(3,4)!*!x(3)!=!f(3)!
283  a(4,1)=r4 !a(4,1) a(4,2) a(4,3) a(4,4)! !x(4)! !f(4)!
284 
285  a(1,2)=z1
286  a(2,2)=z2
287  a(3,2)=z3
288  a(4,2)=z4
289 
290  a(1,3)=r1*z1
291  a(2,3)=r2*z2
292  a(3,3)=r3*z3
293  a(4,3)=r4*z4
294 
295  a(1,4)=1.d0
296  a(2,4)=1.d0
297  a(3,4)=1.d0
298  a(4,4)=1.d0
299 
300  f(1)=u1
301  f(2)=u2
302  f(3)=u3
303  f(4)=u4
304 
305  call ge(4,4,a,f,x,iwrk)
306 
307  u0=x(1)*r0+x(2)*z0+x(3)*r0*z0+x(4)
308 
309  dudr=x(1)+x(3)*z0
310  dudz=x(2)+x(3)*r0
311 
312  return
313  end
314 
315 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
316 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
317  subroutine f_psloop_e(rlop,zlop,psilop,nlop)
318 
319  include 'double.inc'
320  include 'parrc1.inc'
321  include 'comrec.inc'
322  dimension rlop(*),zlop(*),psilop(*)
323 
324  ddx=x(2)-x(1)
325  ddy=y(2)-y(1)
326 
327  do i=1,nlop
328 
329  r0=rlop(i)
330  z0=zlop(i)
331 
332  ic=(r0-x(1))/ddx+1
333  jc=(z0-y(1))/ddy+1
334 
335  r1=x(ic)
336  r2=x(ic+1)
337 
338  z1=y(jc)
339  z2=y(jc+1)
340 
341  u1=ue(ic,jc)
342  u2=ue(ic+1,jc)
343  u3=ue(ic+1,jc+1)
344  u4=ue(ic,jc+1)
345 
346  psilop(i)=blin_(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4 )
347 
348  enddo
349 
350  return
351  end
352 
353 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
354 
355  subroutine psloop_p(rlop,zlop,psilop,nlop)
356 
357  include 'double.inc'
358  include 'parevo.inc'
359  parameter(nkp=njlim)
360  include 'dim.inc'
361  include 'compol.inc'
362  include 'compol_add.inc'
363 
364  real*8 psilop(*),rlop(*),zlop(*)
365 
366  sqrt(arg)=dsqrt(arg)
367 
368  ik_out=0
369 
370  do 100 ik=1,nlop
371 
372  r0=rlop(ik)
373  z0=zlop(ik)
374 
375  call numcel(r0,z0,ic,jc)
376 
377  if(iprlop(ik).eq.1) then !!!
378 
379  if(ic.eq.nr) then
380  write(*,*) 'psiloop_p: bug'
381  write(*,*) 'loop',ik,'is out of box'
382  write(*,*) 'iprlop(ik)=',iprlop(ik)
383  write(*,*) 'program is terminated'
384  stop
385  endif
386 
387  r1=r(ic,jc)
388  r2=r(ic+1,jc)
389  r3=r(ic+1,jc+1)
390  r4=r(ic,jc+1)
391 
392  z1=z(ic,jc)
393  z2=z(ic+1,jc)
394  z3=z(ic+1,jc+1)
395  z4=z(ic,jc+1)
396 
397  u1=psii(ic,jc)
398  u2=psii(ic+1,jc)
399  u3=psii(ic+1,jc+1)
400  u4=psii(ic,jc+1)
401 
402  call blic_d(r0,z0,r1,r2,r3,r4,z1,z2,z3,z4,
403  * u1,u2,u3,u4,u0,dudr,dudz)
404 
405  psilop(ik)=psilop(ik)+u0
406 
407  else !!!
408 
409  if(ic.ne.nr) then
410  write(*,*) 'psiloop_p:loop',ik,'is into the box',ic,nr
411  write(*,*) 'program is terminated'
412  write(*,*) 'you need to consult program provider'
413  stop
414  endif
415 
416  ik_out=ik_out+1
417 
418  do 110 j=2,nt1
419 
420  psilop(ik)=psilop(ik)
421  * -adginl(ik_out,j)*(dgdn(j)+dgdn(j+1))*0.5d0
422 
423  110 continue
424 
425  endif !!!
426 
427  100 continue
428 
429  if(ik_out.ne.nlop_out) then
430  write(*,*) .ne.'psiloop_p:ik_outnlop_out',ik_out,nlop_out
431  write(*,*) 'program is terminated'
432  write(*,*) 'you need to consult program provider'
433  stop
434  endif
435 
436  return
437  end
438 
439 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
440 
441  subroutine bprob_e(rprob,zprob,hrprob,hzprob,nprob,bprob)
442 
443  include 'double.inc'
444  include 'parrc1.inc'
445  include 'comrec.inc'
446  dimension rprob(*),zprob(*),hrprob(*),hzprob(*),bprob(*)
447 
448  ddx=x(2)-x(1)
449  ddy=y(2)-y(1)
450 
451  do i=1,nprob
452 
453  r0=rprob(i)
454  z0=zprob(i)
455 
456  ic=(r0-x(1))/ddx+1
457  jc=(z0-y(1))/ddy+1
458 
459  r1=x(ic)
460  r2=x(ic+1)
461 
462  z1=y(jc)
463  z2=y(jc+1)
464 
465  u1=ue(ic,jc)
466  u2=ue(ic+1,jc)
467  u3=ue(ic+1,jc+1)
468  u4=ue(ic,jc+1)
469 
470  call blin_d(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4,dudr,dudz)
471 
472  br=-dudz/r0
473  bz= dudr/r0
474  bprob(i)=br*hrprob(i)+bz*hzprob(i)
475  enddo
476 
477  return
478  end
479 
480  subroutine extpol
481 
482  include 'double.inc'
483  include 'parevo.inc'
484  parameter(nkp=njlim)
485  include 'dim.inc'
486  include 'parrc1.inc'
487  include 'compol.inc'
488  include 'compol_add.inc'
489  include 'comrec.inc'
490 
491 
492  ddx=x(2)-x(1)
493  ddy=y(2)-y(1)
494 
495  do i=1,nr
496  do j=1,nt
497 
498  r0=r(i,j)
499  z0=z(i,j)
500 
501  ic=(r0-x(1))/ddx+1
502  jc=(z0-y(1))/ddy+1
503 
504  r1=x(ic)
505  r2=x(ic+1)
506 
507  z1=y(jc)
508  z2=y(jc+1)
509 
510  u1=ue(ic,jc)
511  u2=ue(ic+1,jc)
512  u3=ue(ic+1,jc+1)
513  u4=ue(ic,jc+1)
514 
515  psie(i,j)=blin_(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4 )
516 
517  enddo
518  enddo
519 
520  return
521  end
522 
523 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
524  subroutine psi_inter
525 
526  include 'double.inc'
527  include 'parevo.inc'
528  parameter(nkp=njlim)
529  include 'dim.inc'
530  include 'parrc1.inc'
531  include 'compol.inc'
532  include 'compol_add.inc'
533  include 'comrec.inc'
534 
535  dimension psi_sur(ntp)
536 
537  ddx=x(2)-x(1)
538  ddy=y(2)-y(1)
539 
540  !do i=1,nr
541  i=iplas
542  do j=1,nt
543 
544  r0=r(i,j)
545  z0=z(i,j)
546 
547  ic=(r0-x(1))/ddx+1
548  jc=(z0-y(1))/ddy+1
549 
550  r1=x(ic)
551  r2=x(ic+1)
552 
553  z1=y(jc)
554  z2=y(jc+1)
555 
556  u1=u(ic,jc)
557  u2=u(ic+1,jc)
558  u3=u(ic+1,jc+1)
559  u4=u(ic,jc+1)
560 
561  psi_sur(j)=blin_(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4 )
562 
563  enddo
564  !enddo
565 
566  return
567  end
568 
569 
570 
571 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
572  subroutine f_psiful
573 
574  include 'double.inc'
575  include 'parevo.inc'
576  parameter(nkp=njlim)
577  include 'dim.inc'
578  include 'compol.inc'
579  include 'compol_add.inc'
580 
581 
582  do i=1,nr
583  do j=1,nt
584 
585  psi(i,j)=psii(i,j)+psie(i,j)
586 
587  enddo
588  enddo
589 
590 
591  !open(3,file='boncon.wr')
592  !write(3,*) (psi(iplas,j),j=2,nt1)
593  !close(3)
594  !stop
595 
596  return
597  end
598 
599 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
600 
601 
602 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
603 
604  subroutine blin_d(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4,dudr,dudz)
605  include 'double.inc'
606 
607  s=(r2-r1)*(z2-z1)
608 
609  dudr=( (u2-u1)*(z2-z0)+(u3-u4)*(z0-z1) )/s
610  dudz=( (u4-u1)*(r2-r0)+(u3-u2)*(r0-r1) )/s
611 
612  return
613  end
614 
615 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
616 
617  subroutine artfil
618 
619  include 'double.inc'
620  include 'parevo.inc'
621  parameter(nkp=njlim)
622  include 'dim.inc'
623  parameter(nshp=ntp+1)
624  include 'compol.inc'
625  include 'compol_add.inc'
626 
627  dimension xs(nshp),ys(nshp),fun(nshp),dpm(5)
628  common
629  * /c_kpr/kpr
630 
631  1234 continue
632 
633  nsh=1
634 
635  xs(nsh)=r(1,1)
636  ys(nsh)=z(1,1)
637  fun(nsh)=psi(1,1)
638  if(kpr.eq.1) then
639  write(*,*) 'artfil:psi(1,1)',fun(nsh)
640  endif
641 
642  do 200 j=2,nt1
643 
644  nsh=nsh+1
645  xs(nsh)=r(2,j)
646  ys(nsh)=z(2,j)
647  fun(nsh)=psi(2,j)
648 
649  200 continue
650 
651  call deriv5(xs,ys,fun,nsh,5,dpm)
652 
653  clrn= -dpm(1)*0.5d0/rm
654  clzn= -dpm(2)
655 
656  sigma=1.d0
657 
658  clr=sigma*clrn+(1.d0-sigma)*clr
659  clz=sigma*clzn+(1.d0-sigma)*clz
660 
661  do 250 i=1,nr
662  do 250 j=1,nt
663 
664  psi(i,j)=psi(i,j)+clz*z(i,j)+clr*r(i,j)*r(i,j)
665 
666  250 continue
667 
668  return
669  end
670 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
671 
672  subroutine f_xpoint(rx,zx,ix,jx,psx,tet0,kodex)
673 
674  include 'double.inc'
675  parameter(nshp=20)
676  include 'parevo.inc'
677  parameter(nkp=njlim)
678  include 'dim.inc'
679  include 'compol.inc'
680  include 'compol_add.inc'
681 
682  dimension dp(10),xs(nshp),ys(nshp),fun(nshp)
683 
684  sqrt(arg)=dsqrt(arg)
685 
686  numbit=0
687  kodex=0
688 
689  100 numbit=numbit+1
690 
691  rxn=rx
692  zxn=zx
693 
694  ixn=ix
695  jxn=jx
696 
697  dr0=rx-rm
698  dz0=zx-zm
699 
700  ro0=sqrt(dr0**2+dz0**2)
701 
702  tetp=dacos(dr0/ro0)
703  if(dz0.lt.0.d0) then
704  tet0=2.d0*pi-tetp
705  else
706  tet0=tetp
707  endif
708 
709  if(tet0.lt.teta(1)) tet0=tet0+2.d0*pi
710  if(tet0.gt.teta(nt)) tet0=tet0-2.d0*pi
711 
712  do j=1,nt1
713 
714  if(tet0.ge.teta(j) .AND. tet0.lt.teta(j+1)) jc=j
715 
716  enddo
717 
718  dro=1.d9
719 
720  !do i=iplas-3,nr1
721  do i=3,nr1
722 
723  dromn=sqrt((r(i+1,jc)-rx)**2+(z(i+1,jc)-zx)**2)
724  dropl=sqrt((r(i+1,jc+1)-rx)**2+(z(i+1,jc+1)-zx)**2)
725 
726  if(dropl.le.dro .AND. dropl.le.dromn) then
727 
728  dro=dropl
729  ix=i+1
730  jx=jc+1
731 
732  elseif(dromn.le.dro .AND. dromn.le.dropl) then
733 
734  dro=dromn
735  ix=i+1
736  jx=jc
737 
738  endif
739 
740  enddo
741 
742  if(jx.eq.nt) jx=2
743  if(jx.eq.1) jx=nt1
744 
745  if(ix.eq.nr) then !xpoint out of box
746 
747  kodex=1
748  return
749 
750  endif
751 
752  nsh=1
753  xs(nsh)=r(ix,jx)
754  ys(nsh)=z(ix,jx)
755  fun(nsh)=psi(ix,jx)
756 
757  do 400 k=-1,1
758 
759  i= ix+k
760 
761  do 410 l=-1,1
762 
763  j= jx+l
764 
765  if(k.eq.0 .AND. l.eq.0 ) go to 410
766  if(i.gt.nr ) go to 410
767  nsh=nsh+1
768  xs(nsh)=r(i,j)
769  ys(nsh)=z(i,j)
770  fun(nsh)=psi(i,j)
771 
772  410 continue
773  400 continue
774 
775 
776  do 500 k=-2,2,2
777 
778  i= ix+k
779 
780  !do 510 l=-2,2,2
781 
782  !j= jx+l
783  j= jx
784 
785  !if(k.eq.0 .AND. l.eq.0 ) go to 510
786  if(k.eq.0) go to 500
787  if(i.gt.nr ) go to 500
788  nsh=nsh+1
789  xs(nsh)=r(i,j)
790  ys(nsh)=z(i,j)
791  fun(nsh)=psi(i,j)
792 
793  !510 continue
794  500 continue
795 
796  call deriv5(xs,ys,fun,nsh,5,dp)
797 
798  det = dp(3)*dp(5) - dp(4)**2
799 
800  rx = xs(1) + ( dp(2)*dp(4) - dp(1)*dp(5) )/det
801  zx = ys(1) + ( dp(1)*dp(4) - dp(2)*dp(3) )/det
802 
803  rrx=xs(1)
804  zzx=ys(1)
805 
806  psx=fun(1)+ dp(1)*(rx-rrx) + dp(2)*(zx-zzx)
807  + + 0.5d0*dp(3)*(rx-rrx)*(rx-rrx)
808  + + dp(4)*(rx-rrx)*(zx-zzx)
809  + + 0.5d0*dp(5)*(zx-zzx)*(zx-zzx)
810 
811  ! write(6,*) 'rx zx psix',rx,zx,psx
812 
813  xxyy=dp(3)*dp(5)
814 
815 
816  if( (ix.ne.ixn .OR. jx.ne.jxn).AND.(numbit.lt.10) ) then
817 
818  go to 100
819 
820  else
821 
822  xxyy=dp(3)*dp(5)
823  !if(xxyy.ge.0.d0) kodex=1
824 
825  return
826 
827  endif
828 
829  end
830 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
831 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
832 
833  subroutine f_xpoint2(rx,zx,ix,jx,psx,tet0,kodex)
834 
835  include 'double.inc'
836  parameter(nshp=20)
837  include 'parevo.inc'
838  parameter(nkp=njlim)
839  include 'dim.inc'
840  include 'compol.inc'
841  include 'compol_add.inc'
842 
843  dimension dp(10),xs(nshp),ys(nshp),fun(nshp)
844  dimension dp_e(10)
845 
846  sqrt(arg)=dsqrt(arg)
847 
848  numbit=0
849  kodex=0
850 
851  100 numbit=numbit+1
852 
853  rxn=rx
854  zxn=zx
855 
856  ixn=ix
857  jxn=jx
858 
859 
860 !!!!!!!!nearest node definition(ix,ix)
861 
862  dr0=rx-rm
863  dz0=zx-zm
864 
865  ro0=sqrt(dr0**2+dz0**2)
866 
867  tetp=dacos(dr0/ro0)
868  if(dz0.lt.0.d0) then
869  tet0=2.d0*pi-tetp
870  else
871  tet0=tetp
872  endif
873 
874  if(tet0.lt.teta(1)) tet0=tet0+2.d0*pi
875  if(tet0.gt.teta(nt)) tet0=tet0-2.d0*pi
876 
877  do j=1,nt1
878 
879  if(tet0.ge.teta(j) .AND. tet0.lt.teta(j+1)) jc=j
880 
881  enddo
882 
883  dro=1.d9
884 
885  do i=iplas-3,nr1
886 
887  dromn=sqrt((r(i+1,jc)-rx)**2+(z(i+1,jc)-zx)**2)
888  dropl=sqrt((r(i+1,jc+1)-rx)**2+(z(i+1,jc+1)-zx)**2)
889 
890  if(dropl.le.dro .AND. dropl.le.dromn) then
891 
892  dro=dropl
893  ix=i+1
894  jx=jc+1
895 
896  elseif(dromn.le.dro .AND. dromn.le.dropl) then
897 
898  dro=dromn
899  ix=i+1
900  jx=jc
901 
902  endif
903 
904  enddo
905 
906  if(jx.eq.nt) jx=2
907  if(jx.eq.1) jx=nt1
908 
909  if(ix.eq.nr) then !xpoint out of box
910 
911  kodex=1
912  return
913 
914  endif
915 
916 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
917 
918 
919  nsh=1
920  xs(nsh)=r(ix,jx)
921  ys(nsh)=z(ix,jx)
922  fun(nsh)=psii(ix,jx)
923 
924  do 400 k=-1,1
925 
926  i= ix+k
927 
928  do 410 l=-1,1
929 
930  j= jx+l
931 
932  if(k.eq.0 .AND. l.eq.0 ) go to 410
933  if(i.gt.nr ) go to 410
934  nsh=nsh+1
935  xs(nsh)=r(i,j)
936  ys(nsh)=z(i,j)
937  fun(nsh)=psii(i,j)
938 
939  410 continue
940  400 continue
941 
942 
943  do 500 k=-2,2,2
944 
945  i= ix+k
946 
947  !do 510 l=-2,2,2
948 
949  !j= jx+l
950  j= jx
951 
952  !if(k.eq.0 .AND. l.eq.0 ) go to 510
953  if(k.eq.0) go to 500
954  if(i.gt.nr ) go to 500
955  nsh=nsh+1
956  xs(nsh)=r(i,j)
957  ys(nsh)=z(i,j)
958  fun(nsh)=psi(i,j)
959 
960  !510 continue
961  500 continue
962 
963  call deriv5(xs,ys,fun,nsh,5,dp)
964 
965 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
966 
967 
968  call xpoint_e(psie_x,rx,zx,r_ix,z_jx,dp_e,kodex_e)
969 
970 
971 
972 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
973 
974  dp1=dp(1)+dp_e(1)
975  dp2=dp(2)+dp_e(2)
976  dp3=dp(3)+dp_e(3)
977  dp4=dp(4)+dp_e(4)
978  dp5=dp(5)+dp_e(5)
979 
980 
981 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
982 
983 
984  det = dp(3)*dp(5) - dp(4)**2
985 
986  rx = xs(1) + ( dp(2)*dp(4) - dp(1)*dp(5) )/det
987  zx = ys(1) + ( dp(1)*dp(4) - dp(2)*dp(3) )/det
988 
989  rrx=xs(1)
990  zzx=ys(1)
991 
992  psx=fun(1)+ dp(1)*(rx-rrx) + dp(2)*(zx-zzx)
993  + + 0.5d0*dp(3)*(rx-rrx)*(rx-rrx)
994  + + dp(4)*(rx-rrx)*(zx-zzx)
995  + + 0.5d0*dp(5)*(zx-zzx)*(zx-zzx)
996 
997  ! write(6,*) 'rx zx psix',rx,zx,psx
998 
999  xxyy=dp(3)*dp(5)
1000 
1001 
1002  if( (ix.ne.ixn .OR. jx.ne.jxn).AND.(numbit.lt.10) ) then
1003 
1004  go to 100
1005 
1006  else
1007 
1008  xxyy=dp(3)*dp(5)
1009  !if(xxyy.ge.0.d0) kodex=1
1010 
1011  return
1012 
1013  endif
1014 
1015  end
1016 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1017 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1018 
1019  subroutine xpoint_e(ue_x,rx,zx,r_ix,z_jx,dp,kodex)
1020 
1021  include 'double.inc'
1022  parameter(nshp=10)
1023  include 'param.inc'
1024  include 'comblc.inc'
1025 
1026  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
1027 
1028  kodex=0
1029  444 continue
1030 
1031 
1032 
1033  if(rx.ge.rmax .OR. rx.le.rmin .OR.
1034  * zx.ge.zmax .OR. zx.le.zmin) then
1035 ccc write(6,*)' xpoint',numxp,'out of box'
1036 ccc write(6,*)' rx,zx ',rx,zx
1037 
1038  ue_x=-1.d6
1039  kodex=1
1040  return
1041  endif
1042 
1043 c---definition of cell, containig x-point
1044 
1045  do 300 i=1,ni1
1046  icx=i
1047  if( (rx.lt.r(i+1)) .AND. (rx.ge.r(i)) ) go to 301
1048  300 continue
1049  301 continue
1050 
1051  do 310 j=1,nj1
1052  jcx=j
1053  if( (zx.lt.z(j+1)) .AND. (zx.ge.z(j)) ) go to 311
1054  310 continue
1055  311 continue
1056 
1057  if(icx.eq.1 .OR. icx.eq.ni1 .OR.
1058  * jcx.eq.1 .OR. jcx.eq.nj1) then
1059 ccc write(6,*)' xpoint',numxp,'in bound cell'
1060 ccc write(6,*)' icell,jcell ',icx,jcx
1061 
1062  ue_x=-1.d6
1063  kodex=1
1064 
1065  return
1066  endif
1067 
1068 ccc write(6,*) 'icelx,jcelx',icx,jcx,rx,zx
1069 
1070 c---define nearest knote
1071 
1072  sdmin=rmax
1073 
1074  do 320 k=0,1
1075  rr=r(icx+k)
1076  do 325 l=0,1
1077  zz=z(jcx+l)
1078  dlx=dsqrt( (rr-rx)**2+(zz-zx)**2 )
1079  if(dlx.lt.sdmin) then
1080  sdmin=dlx
1081  ix=icx+k
1082  jx=jcx+l
1083  endif
1084  325 continue
1085  320 continue
1086 
1087  555 continue
1088 
1089 ccc write(6,*) 'ix,jx',ix,jx
1090 
1091  r_ix=r(ix)
1092  z_jx=z(jx)
1093 
1094  nsh=1
1095  xs(nsh)=r(ix)
1096  ys(nsh)=z(jx)
1097  fun(nsh)=ue(ix,jx)
1098 
1099  do 400 k=-1,1
1100 
1101  i= ix+k
1102 
1103  do 410 l=-1,1
1104 
1105  j= jx+l
1106 
1107  if(k.eq.0 .AND. l.eq.0 ) go to 410
1108  nsh=nsh+1
1109  xs(nsh)=r(i)
1110  ys(nsh)=z(j)
1111  fun(nsh)=ue(i,j)
1112 
1113  410 continue
1114  400 continue
1115 
1116  call deriv5(xs,ys,fun,nsh,5,dp)
1117 
1118  rrx=xs(1)
1119  zzx=ys(1)
1120 
1121  ue_x=fun(1)+ dp(1)*(rx-rrx) + dp(2)*(zx-zzx)
1122  + + 0.5d0*dp(3)*(rx-rrx)*(rx-rrx)
1123  + + dp(4)*(rx-rrx)*(zx-zzx)
1124  + + 0.5d0*dp(5)*(zx-zzx)*(zx-zzx)
1125 
1126  return
1127  end
1128 
1129 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1130 
1131 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1132 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1133 
1134  subroutine axpnt_e(ue_ax,r_ax,z_ax,i_ax,j_ax,dp,kodax)
1135 
1136  include 'double.inc'
1137  parameter(nshp=10)
1138  include 'param.inc'
1139  include 'comblc.inc'
1140 
1141  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
1142 
1143  kodax=0
1144  444 continue
1145 
1146 
1147 
1148  if(r_ax.ge.rmax .OR. r_ax.le.rmin .OR.
1149  * z_ax.ge.zmax .OR. z_ax.le.zmin) then
1150  write(*,*)'subroutine axpnt_e:'
1151  write(*,*)' axis is out of box'
1152  write(*,*)' r_ax,z_ax ',r_ax,z_ax
1153 
1154  kodax=1
1155  return
1156  endif
1157 
1158 c---definition of cell, containig x-point
1159 
1160  do 300 i=1,ni1
1161  ic_ax=i
1162  if( (r_ax.lt.r(i+1)) .AND. (r_ax.ge.r(i)) ) go to 301
1163  300 continue
1164  301 continue
1165 
1166  do 310 j=1,nj1
1167  jc_ax=j
1168  if( (z_ax.lt.z(j+1)) .AND. (z_ax.ge.z(j)) ) go to 311
1169  310 continue
1170  311 continue
1171 
1172  if(ic_ax.eq.1 .OR. ic_ax.eq.ni1 .OR.
1173  * jc_ax.eq.1 .OR. jc_ax.eq.nj1) then
1174  write(*,*)' ax_point in bound cell'
1175  write(*,*)' icell,jcell ',ic_ax,jc_ax
1176 
1177  kodax=1
1178 
1179  return
1180  endif
1181 
1182 
1183 c---define nearest knote
1184 
1185  sdmin=rmax
1186 
1187  do 320 k=0,1
1188  rr=r(ic_ax+k)
1189  do 325 l=0,1
1190  zz=z(jc_ax+l)
1191  dl_ax=dsqrt( (rr-r_ax)**2+(zz-z_ax)**2 )
1192  if(dl_ax.lt.sdmin) then
1193  sdmin=dl_ax
1194  i_ax=ic_ax+k
1195  j_ax=jc_ax+l
1196  endif
1197  325 continue
1198  320 continue
1199 
1200  555 continue
1201 
1202 ccc write(6,*) 'ix,jx',ix,jx
1203 
1204  r_iax=r(i_ax)
1205  z_jax=z(j_ax)
1206 
1207  nsh=1
1208  xs(nsh)=r(i_ax)
1209  ys(nsh)=z(j_ax)
1210  fun(nsh)=ue(i_ax,j_ax)
1211 
1212  do 400 k=-1,1
1213 
1214  i= i_ax+k
1215 
1216  do 410 l=-1,1
1217 
1218  j= j_ax+l
1219 
1220  if(k.eq.0 .AND. l.eq.0 ) go to 410
1221  nsh=nsh+1
1222  xs(nsh)=r(i)
1223  ys(nsh)=z(j)
1224  fun(nsh)=ue(i,j)
1225 
1226  410 continue
1227  400 continue
1228 
1229  call deriv5(xs,ys,fun,nsh,5,dp)
1230 
1231  rrax=xs(1)
1232  zzax=ys(1)
1233 
1234  ue_ax=fun(1)+ dp(1)*(r_ax-rrax) + dp(2)*(z_ax-zzax)
1235  + + 0.5d0*dp(3)*(r_ax-rrax)*(r_ax-rrax)
1236  + + dp(4)*(r_ax-rrax)*(z_ax-zzax)
1237  + + 0.5d0*dp(5)*(z_ax-zzax)*(z_ax-zzax)
1238 
1239  return
1240  end
1241 
1242 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1243 
1244 
1245 
bprob_e
subroutine bprob_e(rprob, zprob, hrprob, hzprob, nprob, bprob)
Definition: _extp.f:441
bprob_p
subroutine bprob_p(rprob, zprob, hrprob, hzprob, nprob, bprob)
Definition: _extp.f:25
fun
subroutine fun(X, F)
Definition: Ev2.f:10
axpnt_e
subroutine axpnt_e(ue_ax, r_ax, z_ax, i_ax, j_ax, dp, kodax)
Definition: _extp.f:1134
extpol
subroutine extpol
Definition: _extp.f:480
numcel
subroutine numcel(rrk, zzk, icell, jcell)
Definition: _fluxt.f:259
artfil
subroutine artfil
Definition: _extp.f:617
blic_d
subroutine blic_d(r0, z0, r1, r2, r3, r4, z1, z2, z3, z4,
Definition: _extp.f:144
f_xpoint2
subroutine f_xpoint2(rx, zx, ix, jx, psx, tet0, kodex)
Definition: _extp.f:833
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
psi_inter
subroutine psi_inter
Definition: _extp.f:524
b_prob
subroutine b_prob(rprob, zprob, hrprob, hzprob, nprob, bprob)
Definition: _extp.f:13
f_psloop_e
subroutine f_psloop_e(rlop, zlop, psilop, nlop)
Definition: _extp.f:317
blin_d
subroutine blin_d(r0, z0, r1, r2, z1, z2, u1, u2, u3, u4, dudr, dudz)
Definition: _extp.f:604
blin_
real *8 function blin_(r0, z0, r1, r2, z1, z2, u1, u2, u3, u4)
Definition: scu.f:857
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
psiloop
subroutine psiloop(rlop, zlop, psilop, nlop)
Definition: _extp.f:1
ge
subroutine ge(N, NZ, A, X, Y, IP)
Definition: scu.f:420
xpoint_e
subroutine xpoint_e(ue_x, rx, zx, r_ix, z_jx, dp, kodex)
Definition: _extp.f:1019
psloop_p
subroutine psloop_p(rlop, zlop, psilop, nlop)
Definition: _extp.f:355
blic_d_
subroutine blic_d_(r0, z0, r1, r2, r3, r4, z1, z2, z3, z4,
Definition: _extp.f:246
f_psiful
subroutine f_psiful
Definition: _extp.f:572
f_xpoint
subroutine f_xpoint(rx, zx, ix, jx, psx, tet0, kodex)
Definition: _extp.f:672