ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
_mesh.f
Go to the documentation of this file.
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 !!!! grdef axdef regrid remesh prgrid regrid0 grid
3 !!!! fdefln ada reform cpr bonspl loop95 deriv5 ge
4 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5 
6  subroutine renet
7 
8  include 'double.inc'
9  include 'parevo.inc'
10  parameter(nkp=njlim)
11  include 'dim.inc'
12  include 'compol.inc'
13  include 'compol_add.inc'
14  include 'parrc1.inc'
15  include 'comrec.inc'
16  common
17  * /c_kpr/kpr
18 
19  do i=iplas,1,-1
20  if(psia(i).ge.0.05d0) then
21  i_bc=i
22  exit
23  endif
24  enddo
25 
26  do j=2,nt1
27  ronor(nr,j)=1.57d0
28  ro(nr,j)=ro(i_bc,j)*ronor(nr,j)
29  enddo
30 
31  do i=2,iplas
32  do j=2,nt1
33  ro(i,j)=ro(iplas,j)*ronor(i,j)
34  enddo
35  enddo
36 
37  do i=iplas+1,nr-1
38  do j=2,nt1
39  ropla=ro(iplas,j)
40  robon=ro(nr,j)
41  delro=(robon-ropla)/(nr-iplas)
42  ro(i,j)=ropla+delro*(i-iplas)
43  ronor(i,j)=ro(i,j)/ropla
44  enddo
45  enddo
46 
47  do i=2,nr
48  do j=2,nt1
49  r(i,j)=ro(i,j)*cos(teta(j))+rm
50  z(i,j)=ro(i,j)*sin(teta(j))+zm
51  enddo
52  enddo
53 
54  do i=1,nr
55 
56  ro(i,1)=ro(i,nt1)
57  ro(i,nt)=ro(i,2)
58 
59  ronor(i,1)=ronor(i,nt1)
60  ronor(i,nt)=ronor(i,2)
61 
62  r(i,1)=r(i,nt1)
63  r(i,nt)=r(i,2)
64 
65  z(i,1)=z(i,nt1)
66  z(i,nt)=z(i,2)
67 
68  enddo
69 
70 
71  rb_max=rm
72  rb_min=rm
73  zb_max=zm
74  zb_min=zm
75 
76  do j=2,nt1
77  if(r(nr,j).gt.rb_max) rb_max=r(nr,j)
78  if(r(nr,j).lt.rb_min) rb_min=r(nr,j)
79  if(z(nr,j).gt.zb_max) zb_max=z(nr,j)
80  if(z(nr,j).lt.zb_min) zb_min=z(nr,j)
81  enddo
82 
83  if(rb_max.gt.xmax
84  * .OR.
85  * zb_max.gt.ymax
86  * .OR.
87  * rb_min.lt.xmin
88  * .OR.
89  * zb_min.lt.ymin ) then
90 
91  if(kpr.eq.1) then
92  write(*,*) 'eqq is large then rectangular box'
93  call f_wrd
94  pause 'pause '
95  endif
96  endif
97 
98 !!!!!!!!! angle theta redefinition
99  key_ang=1
100  if(key_ang.eq.1) then
101  tet0=teta(2)
102  do j=2,nt1
103  teta(j)=tet0+2.d0*pi*(j-2)/(nt-2.d0)
104  enddo
105  teta(1)=teta(nt1)-2.d0*pi
106  teta(nt)=teta(2)+2.d0*pi
107 
108  do i=2,nr
109  do j=2,nt1
110  r(i,j)=ro(i,j)*cos(teta(j))+rm
111  z(i,j)=ro(i,j)*sin(teta(j))+zm
112  enddo
113  enddo
114 
115  do i=1,nr
116 
117  ro(i,1)=ro(i,nt1)
118  ro(i,nt)=ro(i,2)
119 
120  ronor(i,1)=ronor(i,nt1)
121  ronor(i,nt)=ronor(i,2)
122 
123  r(i,1)=r(i,nt1)
124  r(i,nt)=r(i,2)
125 
126  z(i,1)=z(i,nt1)
127  z(i,nt)=z(i,2)
128 
129  enddo
130  endif
131 
132 
133  return
134  end
135 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
136 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
137 
138  subroutine f_regrid(erro)
139 
140  include 'double.inc'
141  include 'parevo.inc'
142  parameter(nkp=njlim)
143  include 'dim.inc'
144  parameter(ntp4=ntp+4,ntp6=ntp4*6)
145  !!include 'comlmtr.inc'
146  common /com_sp3/ rrksp(ntp4),ccksp(ntp4)
147  include 'compol.inc'
148  include 'compol_add.inc'
149 
150  dimension ron(ntp),rop(ntp),rn(ntp),zn(ntp),tetn(ntp)
151  dimension rob(ntp),teti(ntp),roi(ntp)
152  dimension roh(nrp),rh(nrp),zh(nrp)
153  dimension roplt(nrp,ntp)
154  dimension psiold(nrp,ntp)
155  dimension dp(5)
156  real*8 cwksp(4)
157  common
158  * /c_kpr/kpr
159 
160 ! xzer(x1,x2,v1,v2) = (x1*v2-x2*v1)/(v2-v1)
161 
162  cos(xx)=dcos(xx)
163  sin(xx)=dsin(xx)
164  atan(xx)=datan(xx)
165  asin(xx)=dasin(xx)
166  sqrt(xx)=dsqrt(xx)
167 
168  npro=nr-iplas+1
169 
170  if(ngav.eq.0) then
171  alfa=0.995d0
172  else
173  alfa=0.5d0
174  endif
175 
176  ! if(itin.gt.1) then
177  ! do i=1,nr
178  ! do j=1,nt
179  ! psi(i,j)=alfa*psi(i,j)+(1.0d0-alfa)*psiold(i,j)
180  ! enddo
181  ! enddo
182  ! endif
183 
184  rmold=rm
185  zmold=zm
186 c write(6,*) 'rm,zm,psim',rm,zm,psim
187 
188  call axdef(rma,zma,psima,dp)
189 c write(6,*) 'rma,zma,psima',rma,zma,psima
190 
191 
192  rm=rma*alfa+rmold*(1.0d0-alfa)
193  zm=zma*alfa+zmold*(1.0d0-alfa)
194  psim=psima
195 
196  erro=dsqrt((rmold-rm)**2+(zmold-zm)**2)
197 
198  if(kpr.eq.1) then
199  write(*,*) 'erro mag.axix',erro
200  endif
201 ! pause ' pau'
202 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
203 
204  rxold=rx0
205  zxold=zx0
206 
207  call f_xpoint(rx1,zx1,ixp1,jxp1,psix1,tetx1,kodex1)
208 
209  call f_xpoint(rx2,zx2,ixp2,jxp2,psix2,tetx2,kodex2)
210 
211  if(kpr.eq.1) then
212  write(*,*) 'kodex1,kodex2',kodex1,kodex2
213  endif
214 ! pause ' pau'
215 
216  if(kodex1.eq.0 .AND. kodex2.eq.0) then
217  if(psix1.ge.psix2) then
218  psix0=psix1
219  rx0=rx1
220  zx0=zx1
221  else
222  psix0=psix2
223  rx0=rx2
224  zx0=zx2
225  endif
226  elseif(kodex1.eq.0) then
227  psix0=psix1
228  rx0=rx1
229  zx0=zx1
230  elseif(kodex2.eq.0) then
231  psix0=psix2
232  rx0=rx2
233  zx0=zx2
234  endif
235  errox=dsqrt((rxold-rx0)**2+(zxold-zx0)**2)
236  if(kpr.eq.1) then
237  write(*,*) 'erro x-point',errox
238  endif
239 c write(*,*) 'rx zx',rx0,zx0,psix0
240 
241  psip=psim-alp*(psim-psix0)
242 
243  if(nctrl.eq.1) then
244  do llim=1,nblm
245  call fdefln(psiblm,rblm(llim),zblm(llim))
246  if(psiblm.gt.psip ) then
247  psip=psiblm
248  numlim=llim
249  endif
250  enddo
251  endif
252 
253 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
254 
255  errpsi=0.d0
256 
257  do i=1,nr
258  do j=1,nt
259 
260  psnn=psin(i,j)
261  psin(i,j)=(psi(i,j)-psip)/(psim-psip)
262  delpsn=dabs(psin(i,j)-psnn)
263 
264  if(delpsn.gt.errpsi) then
265  ierm=i
266  jerm=j
267  errpsi=delpsn
268  endif
269  !errpsi=dmax1(errpsi,delpsn)
270 
271  enddo
272  enddo
273 
274  if(kpr.eq.1) then
275  write(6,*) 'i,j errpsi',ierm,jerm,errpsi,(psim-psip)
276  endif
277 c pause ' pau'
278 
279 cc definition of new angle grid teta(j)
280 
281  do j=1,nt
282  tetn(j)=teta(j)
283  enddo
284 
285  do j=1,nt
286 
287  drx=r(nr,j)-rm
288  dzx=z(nr,j)-zm
289 
290  rob(j)=sqrt(drx**2+dzx**2)
291  teta(j)=atan(dzx/drx)
292 
293  enddo
294 
295  do j=2,nt
296 
297  if(teta(j).lt.teta(j-1)) then
298  teta(j)=teta(j)+pi
299  endif
300 
301  if(teta(j).lt.teta(j-1)) then
302  teta(j)=teta(j)+pi
303  endif
304 
305  enddo
306 
307 
308  teta(1)=teta(nt1)-2.d0*pi
309  teta(nt)=teta(2)+2.d0*pi
310 
311 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
312 
313 ! i=2
314 ! zps =psia(i)
315 !
316 ! ps0=psip+zps*(psim-psip)
317 !
318 ! do j=1,nt
319 !
320 ! ttj=teta(j)
321 ! ro2j=(ps0-psim)/( 0.5d0*dp(3)*cos(ttj)**2
322 ! + + dp(4)*cos(ttj)*sin(ttj)
323 ! + + 0.5d0*dp(5)*sin(ttj)**2 )
324 !
325 ! rop(j)=dsqrt(ro2j)
326 !
327 !
328 ! r(i,j)=rm+rop(j)*cos(teta(j))
329 ! z(i,j)=zm+rop(j)*sin(teta(j))
330 !
331 ! enddo ! j-loop
332 
333 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
334 
335 ccc grid moving along rais
336 
337  do i=2,iplas
338 
339  zps =psia(i)
340 
341  do j=1,nt
342 
343  do is=1,nr1
344 
345  ps0 =psin(is,j)
346  !psmn=psin(is-1,j)
347  pspl=psin(is+1,j)
348 
349  ro0 =ro(is,j)
350  !romn=ro(is-1,j)
351  ropl=ro(is+1,j)
352 
353  if(zps.le.ps0 .and. zps.ge.pspl) then
354  zro=((pspl-zps)*ro0-(ps0-zps)*ropl)/(pspl-ps0)
355  go to 822
356  endif
357  !zro=qvadin(zps, psmn,ps0,pspl, romn,ro0,ropl)
358  enddo
359  822 continue
360 
361 
362  ! ps0 =psin(i,j)
363  ! psmn=psin(i-1,j)
364  ! pspl=psin(i+1,j)
365 
366  ! ro0 =ro(i,j)
367  ! romn=ro(i-1,j)
368  ! ropl=ro(i+1,j)
369 
370  ! gradpl=-(pspl-ps0)/(ropl-ro0)
371  ! gradmn=-(ps0-psmn)/(ro0-romn)
372 
373  ! grad=dmax1(gradpl,gradmn)
374 
375  ! !alfa=1.0d0+0.45d0*(q(i)-q(1))/q(1)
376  ! alfa=1.5d0
377 
378  ! if(ngav.eq.1) grad=grad*alfa
379 
380  ! zro=ro0-(psia(i)-ps0)/grad
381 
382  ron(j)=zro*alfa+ro(i,j)*(1.d0-alfa)
383 
384 
385  enddo
386 
387  do j=1,nt
388 
389  rnj=rmold+ron(j)*cos(tetn(j))
390  znj=zmold+ron(j)*sin(tetn(j))
391 
392  drx=rnj-rm
393  dzx=znj-zm
394 
395  roi(j)=sqrt(drx**2+dzx**2)
396  teti(j)=atan(dzx/drx)
397 
398  enddo
399 
400  do j=2,nt
401 
402  if(teti(j).lt.teti(j-1)) then
403  teti(j)=teti(j)+pi
404  endif
405 
406  if(teti(j).lt.teti(j-1)) then
407  teti(j)=teti(j)+pi
408  endif
409 
410  enddo
411 
412 
413  teti(1)=teti(nt1)-2.d0*pi
414  teti(nt)=teti(2)+2.d0*pi
415 
416  roi(1)=roi(nt1)
417  roi(nt)=roi(2)
418 
419 
420  do j=2,nt1
421 
422  tetv=teta(j)
423 
424  if(tetv.lt.teti(1)) tetv=tetv+2.d0*pi
425  if(tetv.gt.teti(nt)) tetv=tetv-2.d0*pi
426 
427  do jj=1,nt1
428 
429 
430  tt0 =teti(jj)
431  !ttmn=teti(j-1)
432  ttpl=teti(jj+1)
433 
434  ro0 =roi(jj)
435  !romn=roi(j-1)
436  ropl=roi(jj+1)
437 
438  if(tetv.le.ttpl .and. tetv.ge.tt0) then
439  zro=((ttpl-tetv)*ro0+(tetv-tt0)*ropl)/(ttpl-tt0)
440  go to 1047
441  endif
442  enddo
443 
444  1047 continue
445 
446  ! zro=qvadin(teta(j), ttmn,tt0,ttpl, romn,ro0,ropl)
447  rop(j)=zro
448 
449  roerr=dabs(ron(j)-ro(i,j))
450  if(roerr.gt.erro) then
451  erro=roerr
452  ierm=i
453  jerm=j
454  endif
455 
456  roplt(i,j)=ron(j)-ro(i,j)
457 
458 
459  r(i,j)=rop(j)*cos(teta(j))+rm
460  z(i,j)=rop(j)*sin(teta(j))+zm
461 
462  enddo
463 
464  enddo ! i-loop
465 
466  if(kpr.eq.1) then
467  write(*,*) 'erro max.',erro
468  write(*,*) 'i,j erro',ierm,jerm,erro
469  endif
470 c pause 'pau'
471 
472  do i=2,iplas
473  do j=2,nt1
474 
475  ro(i,j)=dsqrt((r(i,j)-rm)**2+(z(i,j)-zm)**2)
476 
477  enddo
478  enddo ! i-loop
479 
480 
481  do j=2,nt1
482 
483  ro(nr,j)=rob(j)
484 
485  ro(1,j)=0.d0
486  r(1,j)=rm
487  z(1,j)=zm
488 
489  enddo
490 
491 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
492 
493  do j=2,nt1
494 
495  spro=ro(nr,j)-ro(iplas-1,j)
496 
497  droj=ro(iplas,j)-ro(iplas-1,j)
498  cpro=cpr(droj,npro,spro)
499  row=ro(iplas,j)
500  dro=droj
501 
502  do 220 i=iplas+1,nr1
503 
504  dro=dro*cpro
505  row=row+dro
506 
507  rh(i)=rm+row*cos(teta(j))
508  zh(i)=zm+row*sin(teta(j))
509 
510  roh(i)=row
511 
512  do is=2,nr1
513 
514  ps0 =psin(is,j)
515  pspl=psin(is+1,j)
516 
517  ro0 =ro(is,j)
518  ropl=ro(is+1,j)
519 
520  if(row.le.ropl .and. row.ge.ro0) then
521  isw=is
522  go to 755
523  endif
524 
525  enddo
526 
527  755 continue
528 
529  ro2=ropl
530  ro1=ro0
531 
532  ps2=pspl
533  ps1=ps0
534 
535  psia(i)=( ps1*(ro2-row)+ps2*(row-ro1) )/(ro2-ro1)
536 
537  220 continue
538 
539  do 250 i=iplas+1,nr1
540 
541  roerr=dabs(roh(i)-ro(i,j))
542  erro=dmax1(erro,roerr)
543 
544  ro(i,j)=roh(i)
545 
546  r(i,j)=roh(i)*cos(teta(j))+rm
547  z(i,j)=roh(i)*sin(teta(j))+zm
548 
549  psin(i,j)=psia(i)
550 
551  250 continue
552 
553  enddo
554 
555  do i=1,nr
556  do j=2,nt1
557 
558  ronor(i,j)=ro(i,j)/ro(iplas,j)
559 
560  enddo
561  enddo
562 
563  do 25 i=1,nr
564 
565  ro(i,1)=ro(i,nt1)
566  ro(i,nt)=ro(i,2)
567 
568  ronor(i,1)=ronor(i,nt1)
569  ronor(i,nt)=ronor(i,2)
570 
571  r(i,1)=r(i,nt1)
572  r(i,nt)=r(i,2)
573 
574  z(i,1)=z(i,nt1)
575  z(i,nt)=z(i,2)
576 
577  roplt(i,1)=roplt(i,nt1)
578  roplt(i,nt)=roplt(i,2)
579 
580  25 continue
581 
582 
583  do 30 i=1,iplas
584  do 30 j=1,nt
585 
586  psin(i,j)=psia(i)
587 
588  30 continue
589 
590  do i=1,nr
591  do j=1,nt
592 
593  psiold(i,j)=psip+psin(i,j)*(psim-psip)
594 
595  enddo
596  enddo
597 
598 
599 
600 ! if(ngav.eq.10) then
601 ! open(1,file='erro.wr')
602 ! write(1,*) iplas,nt
603 ! write(1,*) ((r(i,j),i=1,iplas),j=1,nt)
604 ! write(1,*) ((z(i,j),i=1,iplas),j=1,nt)
605 ! write(1,*) ((roplt(i,j),i=1,iplas),j=1,nt)
606 ! close(1)
607 ! call f_wrd
608 ! pause 'erro.wr '
609 ! endif
610 
611  return
612  end
613 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
614 
615 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
616 
617  subroutine f_remesh(erro)
618 
619  include 'double.inc'
620  include 'parevo.inc'
621  parameter(nkp=njlim)
622  include 'dim.inc'
623  include 'compol.inc'
624  include 'compol_add.inc'
625 
626  psimax=psi(1,2)
627  imax=1
628  jmax=2
629 
630  do i=2,iplas1
631  do j=2,nt1
632 
633  if(psi(i,j).gt.psimax) then
634  psimax=psi(i,j)
635  imax=i
636  jmax=j
637  endif
638 
639  enddo
640  enddo
641 
642 
643 
644  ! if(kstep.ge.5 .AnD. kstep.le.8 .AnD. iter/2*2.eq.iter) then
645  ! call f_prgrid(imax,jmax,erro)
646  ! return
647  ! endif
648 
649 
650 
651  if(imax.ge.2) then
652  call f_prgrid(imax,jmax,erro)
653  else
654  call f_regrid0(erro)
655  endif
656 
657  return
658  end
659 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
660 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
661 
662  subroutine f_prgrid(imax,jmax,erro)
663 
664  include 'double.inc'
665  include 'parevo.inc'
666  parameter(nkp=njlim)
667  include 'dim.inc'
668  parameter(nshp=ntp+1)
669  !!include 'comlmtr.inc'
670  include 'compol.inc'
671  include 'compol_add.inc'
672 
673  dimension ron(ntp),rop(ntp),rn(ntp),zn(ntp),tetn(ntp)
674  dimension rob(ntp),teti(ntp),roi(ntp)
675  dimension xs(nshp),ys(nshp),fun(nshp),dp(5)
676  dimension roh(nrp),rh(nrp),zh(nrp)
677  dimension roplt(nrp,ntp)
678  dimension psiold(nrp,ntp)
679  common
680  * /c_kpr/kpr
681  dimension dp_ex(5)
682 
683  atan(xx)=datan(xx)
684  acos(xx)=dacos(xx)
685  sqrt(xx)=dsqrt(xx)
686 
687 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
688 !!!!!!!!!! new position of magnetic axis !!!!!!!!!!!!!!!!!!!!!
689 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
690  write(17,*) 'prgrid:::'
691  nctrli=nctrl
692 
693  npro=nr-iplas+1
694 
695 
696  rimjm=r(imax,jmax)
697  zimjm=z(imax,jmax)
698 
699  !call axpnt_e(ue_ax,rimjm,zimjm,i_ax,j_ax,dp_ex,kodax)
700 
701 
702  rm0=r(1,1)
703  zm0=z(1,1)
704  psim0=psi(1,1)
705 
706  nsh=1
707 
708  xs(nsh)=r(imax,jmax)
709  ys(nsh)=z(imax,jmax)
710  fun(nsh)=psi(imax,jmax)
711 
712  if(imax.eq.1) then
713  do j=2,nt1
714  nsh=nsh+1
715  xs(nsh)=r(2,j)
716  ys(nsh)=z(2,j)
717  fun(nsh)=psi(2,j)
718  enddo
719  elseif(imax.eq.2) then
720 
721  nsh=nsh+1
722  xs(nsh)=r(1,2)
723  ys(nsh)=z(1,2)
724  fun(nsh)=psi(1,2)
725  do j=2,nt1
726  nsh=nsh+1
727  xs(nsh)=r(3,j)
728  ys(nsh)=z(3,j)
729  fun(nsh)=psi(3,j)
730  enddo
731 
732  else
733 
734  do k=-1,1
735  i= imax+k
736  do l=-1,1
737  j= jmax+l
738  if(i.ne.imax .OR. j.ne.jmax) then
739  nsh=nsh+1
740  xs(nsh)=r(i,j)
741  ys(nsh)=z(i,j)
742  fun(nsh)=psi(i,j)
743  endif
744  enddo
745  enddo
746 
747  endif
748 
749  call deriv5(xs,ys,fun,nsh,5,dp)
750 
751  !do kdp=1,5
752  ! dp(kdp)=dp(kdp) + dp_ex(kdp)
753  !enddo
754 
755  det = dp(3)*dp(5) - dp(4)**2
756  rm = xs(1) + ( dp(2)*dp(4) - dp(1)*dp(5) )/det
757  zm = ys(1) + ( dp(1)*dp(4) - dp(2)*dp(3) )/det
758 
759  erro=dsqrt( (rm-rm0)**2+(zm-zm0)**2 )
760 
761  rmx=r(imax,jmax)
762  zmx=z(imax,jmax)
763 
764  psim=fun(1)+ dp(1)*(rm-rmx) + dp(2)*(zm-zmx)
765  + + 0.5d0*dp(3)*(rm-rmx)*(rm-rmx)
766  + + dp(4)*(rm-rmx)*(zm-zmx)
767  + + 0.5d0*dp(5)*(zm-zmx)*(zm-zmx)
768 
769  if(kpr.eq.1) then
770  write(*,*) 'rm zm psim',rm,zm,psim
771  write(*,*) 'prgrid:errma',erro
772  endif
773 
774 
775 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
776 !!!!!!!!!! new position of x-point !!!!!!!!!!!!!!!!!!!!!
777 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
778 
779  rxold=rx0
780  zxold=zx0
781 
782  call f_xpoint(rx1,zx1,ixp1,jxp1,psix1,tetx1,kodex1)
783 
784  call f_xpoint(rx2,zx2,ixp2,jxp2,psix2,tetx2,kodex2)
785 
786  if(kpr.eq.1) then
787  write(*,*) 'kodex1,kodex2',kodex1,kodex2
788  endif
789 ! pause ' pau'
790  kodxp=kodex1*kodex2
791  if(kodxp.ne.0) then
792 
793  if(kpr.eq.1) then
794  write(*,*) ' all x-points out of box'
795  write(*,*) ' only limiter case '
796  endif
797  nctrli=1
798 
799  !write(6,*) ' program cannot run more '
800  !stop
801 
802  psip=-1.d12
803 
804  else !kodxp=0
805  psipn=psip
806 
807  if(kodex1.eq.0 .AND. kodex2.eq.0) then
808  if(psix1.ge.psix2) then
809  psix0=psix1
810  rx0=rx1
811  zx0=zx1
812  else
813  psix0=psix2
814  rx0=rx2
815  zx0=zx2
816  endif
817  elseif(kodex1.eq.0) then
818  psix0=psix1
819  rx0=rx1
820  zx0=zx1
821  elseif(kodex2.eq.0) then
822  psix0=psix2
823  rx0=rx2
824  zx0=zx2
825  endif
826  errox=dsqrt((rxold-rx0)**2+(zxold-zx0)**2)
827  if(kpr.eq.1) then
828  write(*,*) 'erro x-point',errox
829  endif
830 c write(6,*) 'rx zx',rx0,zx0,psix0
831 
832  psip=psim-alp*(psim-psix0)
833  if(itin.gt.1) then
834  zwgt=1.0d0
835  psipr=zwgt*psip+(1.d0-zwgt)*psipn
836  psip=dmax1(psip,psipr)
837  endif
838  endif
839 
840  if(nctrl.eq.1) then
841  numlim=0
842  do llim=1,nblm
843  call fdefln(psiblm,rblm(llim),zblm(llim))
844  if(psiblm.gt.psip ) then
845  psip=psiblm
846  numlim=llim
847  endif
848  enddo
849  endif
850  if(kpr.eq.1) then
851  write(*,*) 'numlim',numlim
852  endif
853 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
854 !!!!!!!!!!! nomalized poloidal flux definition !!!!!!!!!!!
855 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
856 
857  errpsi=0.d0
858 
859  do i=1,nr
860  do j=1,nt
861 
862  psnn=psin(i,j)
863  psin(i,j)=(psi(i,j)-psip)/(psim-psip)
864  delpsn=dabs(psin(i,j)-psnn)
865 
866  if(delpsn.gt.errpsi) then
867  ierm=i
868  jerm=j
869  errpsi=delpsn
870  endif
871 
872  !errpsi=dmax1(errpsi,delpsn)
873 
874  enddo
875  enddo
876 
877 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
878 !!!!!!!!!! definition of new angle grid teta(j) !!!!!!!!!!
879 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
880  do j=1,nt
881  tetn(j)=teta(j)
882  enddo
883 
884  do j=1,nt
885  drx=r(nr,j)-rm
886  dzx=z(nr,j)-zm
887  robj=sqrt(drx**2+dzx**2)
888  rob(j)=robj
889  tetp=acos(drx/robj)
890  if(dzx.lt.0.d0) then
891  teta(j)=-tetp
892  else
893  teta(j)=tetp
894  endif
895  enddo
896 
897  do j=2,nt
898 
899  if(teta(j).lt.teta(j-1)) then
900  teta(j)=teta(j)+2.d0*pi
901  endif
902 
903  if(teta(j).lt.teta(j-1)) then
904  teta(j)=teta(j)+2.d0*pi
905  endif
906 
907  enddo
908 
909  teta(1)=teta(nt1)-2.d0*pi
910  teta(nt)=teta(2)+2.d0*pi
911 
912 ccc grid moving along rais
913 
914  psinm0=psin(1,2)
915 
916  do i=3,iplas1
917  if(psia(i).lt.psinm0) then
918  !ibeg=i
919  ibeg=i+1
920  go to 10
921  endif
922  enddo
923 
924  10 continue
925 
926  do i=ibeg,iplas
927  zps =psia(i)
928  do j=1,nt
929  do is=1,nr1
930  ps0 =psin(is,j)
931  pspl=psin(is+1,j)
932  ro0 =ro(is,j)
933  ropl=ro(is+1,j)
934  if(zps.le.ps0 .and. zps.ge.pspl) then
935  isc=is
936  go to 822
937  endif
938  enddo
939 
940  822 continue
941 
942  grad=-(pspl-ps0)/(ropl-ro0)
943 
944  zro=ro0-(psia(i)-ps0)/grad
945  ron(j)=zro
946  enddo
947 
948  do j=1,nt
949 
950  rnj=rm0+ron(j)*cos(tetn(j))
951  znj=zm0+ron(j)*sin(tetn(j))
952 
953  drx=rnj-rm
954  dzx=znj-zm
955 
956  roi(j)=sqrt(drx**2+dzx**2)
957  tetp=acos(drx/roi(j))
958  if(dzx.lt.0.d0) then
959  teti(j)=-tetp
960  else
961  teti(j)=tetp
962  endif
963 
964  enddo
965 
966  do j=2,nt
967 
968  if(teti(j).lt.teti(j-1)) then
969  teti(j)=teti(j)+2.d0*pi
970  endif
971 
972  if(teti(j).lt.teti(j-1)) then
973  teti(j)=teti(j)+2.d0*pi
974  endif
975 
976  enddo
977 
978  teti(1)=teti(nt1)-2.d0*pi
979  teti(nt)=teti(2)+2.d0*pi
980 
981  roi(1)=roi(nt1)
982  roi(nt)=roi(2)
983 
984  do j=2,nt1
985 
986  tetv=teta(j)
987 
988  if(tetv.lt.teti(1)) tetv=tetv+2.d0*pi
989  if(tetv.lt.teti(1)) tetv=tetv+2.d0*pi
990  if(tetv.gt.teti(nt)) tetv=tetv-2.d0*pi
991  if(tetv.gt.teti(nt)) tetv=tetv-2.d0*pi
992 
993  do jj=1,nt1
994 
995  tt0 =teti(jj)
996  ttpl=teti(jj+1)
997 
998  ro0 =roi(jj)
999  ropl=roi(jj+1)
1000 
1001  if(tetv.le.ttpl .and. tetv.ge.tt0) then
1002  zro=((ttpl-tetv)*ro0+(tetv-tt0)*ropl)/(ttpl-tt0)
1003  go to 1047
1004  endif
1005  enddo
1006 
1007  1047 continue
1008 
1009  rop(j)=zro
1010 
1011  r(i,j)=rop(j)*cos(teta(j))+rm
1012  z(i,j)=rop(j)*sin(teta(j))+zm
1013 
1014  enddo
1015 
1016  enddo ! i-loop
1017 
1018 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1019 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1020 
1021  i=2
1022  zps =psia(i)
1023  ps0=psip+zps*(psim-psip)
1024 
1025  do j=1,nt
1026 
1027  ttj=teta(j)
1028  ro2j=(ps0-psim)/( 0.5d0*dp(3)*cos(ttj)**2
1029  + + dp(4)*cos(ttj)*sin(ttj)
1030  + + 0.5d0*dp(5)*sin(ttj)**2 )
1031 
1032  rodeb=ro2j
1033 
1034  rop(j)=dsqrt(ro2j)
1035 
1036  r(i,j)=rm+rop(j)*cos(teta(j))
1037  z(i,j)=zm+rop(j)*sin(teta(j))
1038 
1039  enddo ! j-loop
1040 
1041 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1042 
1043  if(ibeg.gt.3) then
1044 
1045  do j=2,nt1
1046 
1047  ros1=(r(2,j)-rm)**2+(z(2,j)-zm)**2
1048  ros2=(r(ibeg,j)-rm)**2+(z(ibeg,j)-zm)**2
1049  psa1=psia(2)
1050  psa2=psia(ibeg)
1051 
1052  do i=3,ibeg-1
1053 
1054  psai=psia(i)
1055  rosi=( (psa2-psai)*ros1+(psai-psa1)*ros2 )/(psa2-psa1)
1056  rosi=dsqrt(rosi)
1057 
1058  r(i,j)=rm+rosi*cos(teta(j))
1059  z(i,j)=zm+rosi*sin(teta(j))
1060 
1061  enddo
1062 
1063  enddo
1064 
1065  endif
1066 
1067  do i=2,iplas
1068  do j=2,nt1
1069 
1070  ro(i,j)=dsqrt((r(i,j)-rm)**2+(z(i,j)-zm)**2)
1071 
1072  enddo
1073  enddo
1074 
1075 
1076 
1077 
1078 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1079 
1080  do j=2,nt1
1081 
1082  spro=rob(j)-ro(iplas-1,j)
1083 
1084  droj=ro(iplas,j)-ro(iplas-1,j)
1085  !cpro=cpr(droj,npro,spro)
1086  row=ro(iplas,j)
1087  !dro=droj
1088  dro=(rob(j)-ro(iplas,j))/(nr-iplas)
1089 
1090  do 220 i=iplas+1,nr1
1091 
1092  !dro=dro*cpro
1093  row=row+dro
1094 
1095  rh(i)=rm+row*cos(teta(j))
1096  zh(i)=zm+row*sin(teta(j))
1097 
1098  roh(i)=row
1099 
1100  do is=2,nr1
1101 
1102  ps0 =psin(is,j)
1103  pspl=psin(is+1,j)
1104 
1105  ro0 =ro(is,j)
1106  ropl=ro(is+1,j)
1107 
1108  if(row.le.ropl .and. row.ge.ro0) then
1109  isw=is
1110  go to 755
1111  endif
1112 
1113  enddo
1114 
1115  755 continue
1116 
1117  ro2=ropl
1118  ro1=ro0
1119 
1120  ps2=pspl
1121  ps1=ps0
1122 
1123  psia(i)=( ps1*(ro2-row)+ps2*(row-ro1) )/(ro2-ro1)
1124 
1125  220 continue
1126 
1127  do 250 i=iplas+1,nr1
1128 
1129  roerr=dabs(roh(i)-ro(i,j))
1130  erro=dmax1(erro,roerr)
1131 
1132  ro(i,j)=roh(i)
1133 
1134  r(i,j)=roh(i)*cos(teta(j))+rm
1135  z(i,j)=roh(i)*sin(teta(j))+zm
1136 
1137  psin(i,j)=psia(i)
1138 
1139  250 continue
1140 
1141  enddo
1142 
1143  do j=2,nt1
1144 
1145  ro(nr,j)=rob(j)
1146  !r(iplas,j)=rob(j)*cos(teta(j))+rm
1147  !z(iplas,j)=rob(j)*sin(teta(j))+zm
1148 
1149  ro(1,j)=0.d0
1150  r(1,j)=rm
1151  z(1,j)=zm
1152 
1153  enddo
1154 
1155  do 25 i=1,nr
1156 
1157  ro(i,1)=ro(i,nt1)
1158  ro(i,nt)=ro(i,2)
1159 
1160  r(i,1)=r(i,nt1)
1161  r(i,nt)=r(i,2)
1162 
1163  z(i,1)=z(i,nt1)
1164  z(i,nt)=z(i,2)
1165 
1166  roplt(i,1)=roplt(i,nt1)
1167  roplt(i,nt)=roplt(i,2)
1168 
1169  25 continue
1170 
1171  do 30 i=1,iplas
1172  do 30 j=1,nt
1173 
1174  psin(i,j)=psia(i)
1175 
1176  30 continue
1177 
1178  do i=1,nr
1179  do j=1,nt
1180 
1181  psiold(i,j)=psip+psin(i,j)*(psim-psip)
1182 
1183  enddo
1184  enddo
1185 
1186  call f_wrd
1187 
1188  return
1189  end
1190 
1191 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1192 
1193 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1194 
1195  subroutine regrid0_(erro)
1196 
1197  include 'double.inc'
1198  include 'parevo.inc'
1199  parameter(nkp=njlim)
1200  include 'dim.inc'
1201  !!include 'comlmtr.inc'
1202  include 'compol.inc'
1203  include 'compol_add.inc'
1204 
1205  dimension ron(ntp),rop(ntp),rn(ntp),zn(ntp),tetn(ntp)
1206  dimension rob(ntp),teti(ntp),roi(ntp)
1207  dimension roplt(nrp,ntp)
1208  dimension dp(5)
1209  dimension roh(nrp),rh(nrp),zh(nrp)
1210  dimension psiold(nrp,ntp)
1211  common
1212  * /c_kpr/kpr
1213 
1214 ! xzer(x1,x2,v1,v2) = (x1*v2-x2*v1)/(v2-v1)
1215 
1216  cos(xx)=dcos(xx)
1217  sin(xx)=dsin(xx)
1218  atan(xx)=datan(xx)
1219  asin(xx)=dasin(xx)
1220  acos(xx)=dacos(xx)
1221  sqrt(xx)=dsqrt(xx)
1222 
1223  ! write(*,*) 'regrid0:::'
1224 
1225  nctrli=nctrl
1226 
1227  npro=nr-iplas+1
1228 
1229  if(ngav.eq.0) then
1230  alfa=1.0d0
1231  else
1232  alfa=1.0d0
1233  endif
1234 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1235 
1236  if(itin.gt.1) then
1237  do i=1,nr1
1238  do j=1,nt
1239  psi(i,j)=alfa*psi(i,j)+(1.d0-alfa)*psiold(i,j)
1240  enddo
1241  enddo
1242  endif
1243 
1244 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1245 
1246  rmold=rm
1247  zmold=zm
1248 
1249  if(kpr.eq.1) then
1250  write(*,*) 'rm,zm,psim',rm,zm,psim
1251  endif
1252  call axdef(rma,zma,psima,dp)
1253 
1254  if(kpr.eq.1) then
1255  write(*,*) 'rma,zma,psima',rma,zma,psima
1256  endif
1257  rm=rma
1258  zm=zma
1259  psim=psima
1260 
1261  erro=dsqrt((rmold-rm)**2+(zmold-zm)**2)
1262 
1263  if(kpr.eq.1) then
1264  write(*,*) 'erro mag.axis',erro
1265  endif
1266  errpsi=0.d0
1267 
1268 
1269 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1270 
1271  rxold=rx0
1272  zxold=zx0
1273 
1274  call f_xpoint(rx1,zx1,ixp1,jxp1,psix1,tetx1,kodex1)
1275 
1276  call f_xpoint(rx2,zx2,ixp2,jxp2,psix2,tetx2,kodex2)
1277 
1278  if(kpr.eq.1) then
1279  write(*,*) 'kodex1,kodex2',kodex1,kodex2
1280  endif
1281 
1282  kodxp=kodex1*kodex2
1283  if(kodxp.ne.0) then
1284 
1285  if(kpr.eq.1) then
1286  write(*,*) ' all x-points out of box'
1287  write(*,*) ' only limiter case '
1288  endif
1289  nctrli=1
1290 
1291  !write(6,*) ' program cannot run more '
1292  !stop
1293  psip=-1.d12
1294 
1295  else !kodxp=0
1296 
1297  if(kodex1.eq.0 .AND. kodex2.eq.0) then
1298  if(psix1.ge.psix2) then
1299  psix0=psix1
1300  rx0=rx1
1301  zx0=zx1
1302  else
1303  psix0=psix2
1304  rx0=rx2
1305  zx0=zx2
1306  endif
1307  elseif(kodex1.eq.0) then
1308  psix0=psix1
1309  rx0=rx1
1310  zx0=zx1
1311  elseif(kodex2.eq.0) then
1312  psix0=psix2
1313  rx0=rx2
1314  zx0=zx2
1315  endif
1316 
1317  errox=dsqrt((rxold-rx0)**2+(zxold-zx0)**2)
1318  if(kpr.eq.1) then
1319  write(*,*) 'erro x-point',errox
1320  write(*,*) 'rx zx',rx0,zx0,psix0
1321  endif
1322 
1323  psip=psim-alp*(psim-psix0)
1324 
1325  endif
1326 
1327  if(nctrl.eq.1) then
1328  numlim=0
1329  do llim=1,nblm
1330 
1331  call fdefln(psiblm,rblm(llim),zblm(llim))
1332 
1333  if(psiblm.gt.psip) then
1334  psip=psiblm
1335  numlim=llim
1336  endif
1337 
1338  enddo
1339  endif
1340  if(kpr.eq.1) then
1341  write(*,*) 'numlim',numlim
1342  endif
1343 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1344 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1345 
1346 
1347  do i=1,nr
1348  do j=1,nt
1349 
1350  psnn=psin(i,j)
1351  psin(i,j)=(psi(i,j)-psip)/(psim-psip)
1352  delpsn=dabs(psin(i,j)-psnn)
1353 
1354  if(delpsn.gt.errpsi) then
1355 
1356  ierm=i
1357  jerm=j
1358  errpsi=delpsn
1359 
1360  endif
1361 
1362  enddo
1363  enddo
1364 
1365  if(kpr.eq.1) then
1366  write(*,*) 'i,j errpsi',ierm,jerm,errpsi,(psim-psip)
1367  endif
1368 cc definition of new angle grid teta(j)
1369 
1370  do j=1,nt
1371  tetn(j)=teta(j)
1372  enddo
1373 
1374  do j=1,nt
1375 
1376  drx=r(nr,j)-rm
1377  dzx=z(nr,j)-zm
1378 
1379  robj=sqrt(drx**2+dzx**2)
1380  rob(j)=robj
1381  tetp=acos(drx/robj)
1382  if(dzx.lt.0.d0) then
1383  teta(j)=-tetp
1384  else
1385  teta(j)=tetp
1386  endif
1387 
1388  enddo
1389 
1390  do j=2,nt
1391 
1392  if(teta(j).lt.teta(j-1)) then
1393  teta(j)=teta(j)+2.d0*pi
1394  endif
1395 
1396  if(teta(j).lt.teta(j-1)) then
1397  teta(j)=teta(j)+2.d0*pi
1398  endif
1399 
1400  enddo
1401 
1402  teta(1)=teta(nt1)-2.d0*pi
1403  teta(nt)=teta(2)+2.d0*pi
1404 
1405 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1406 
1407  i=2
1408 
1409  zps =psia(i)
1410  ps0=psip+zps*(psim-psip)
1411 
1412  do j=1,nt
1413 
1414  ttj=teta(j)
1415 
1416  ro2j=(ps0-psim)/( 0.5d0*dp(3)*cos(ttj)**2
1417  + + dp(4)*cos(ttj)*sin(ttj)
1418  + + 0.5d0*dp(5)*sin(ttj)**2 )
1419 
1420  rop(j)=dsqrt(ro2j)
1421 
1422  r(i,j)=rm+rop(j)*cos(teta(j))
1423  z(i,j)=zm+rop(j)*sin(teta(j))
1424 
1425  enddo ! j-loop
1426 
1427 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1428 
1429 ccc grid moving along rais
1430 
1431  do i=3,iplas
1432 
1433  zps =psia(i)
1434 
1435  do j=1,nt
1436 
1437  do is=1,nr1
1438 
1439  ps0 =psin(is,j)
1440 ! !psmn=psin(is-1,j)
1441  pspl=psin(is+1,j)
1442 
1443  ro0 =ro(is,j)
1444 ! !romn=ro(is-1,j)
1445  ropl=ro(is+1,j)
1446 
1447  if(zps.le.ps0 .and. zps.ge.pspl) then
1448  zro=((pspl-zps)*ro0-(ps0-zps)*ropl)/(pspl-ps0)
1449  go to 822
1450  endif
1451 ! !zro=qvadin(zps, psmn,ps0,pspl, romn,ro0,ropl)
1452  enddo
1453  822 continue
1454 
1455 
1456 
1457  ron(j)=zro
1458 
1459  enddo
1460 
1461  do j=1,nt
1462 
1463  rnj=rmold+ron(j)*cos(tetn(j))
1464  znj=zmold+ron(j)*sin(tetn(j))
1465 
1466  drx=rnj-rm
1467  dzx=znj-zm
1468 
1469  roi(j)=sqrt(drx**2+dzx**2)
1470  tetp=acos(drx/roi(j))
1471  if(dzx.lt.0.d0) then
1472  teti(j)=-tetp
1473  else
1474  teti(j)=tetp
1475  endif
1476 
1477  enddo
1478 
1479  do j=2,nt
1480 
1481  if(teti(j).lt.teti(j-1)) then
1482  teti(j)=teti(j)+2.d0*pi
1483  endif
1484 
1485  if(teti(j).lt.teti(j-1)) then
1486  teti(j)=teti(j)+2.d0*pi
1487  endif
1488 
1489  enddo
1490 
1491  teti(1)=teti(nt1)-2.d0*pi
1492  teti(nt)=teti(2)+2.d0*pi
1493 
1494  roi(1)=roi(nt1)
1495  roi(nt)=roi(2)
1496 
1497  do j=2,nt1
1498 
1499  tetv=teta(j)
1500 
1501  if(tetv.lt.teti(1)) tetv=tetv+2.d0*pi
1502  if(tetv.lt.teti(1)) tetv=tetv+2.d0*pi
1503  if(tetv.gt.teti(nt)) tetv=tetv-2.d0*pi
1504  if(tetv.gt.teti(nt)) tetv=tetv-2.d0*pi
1505 
1506  do jj=1,nt1
1507 
1508  tt0 =teti(jj)
1509  !ttmn=teti(j-1)
1510  ttpl=teti(jj+1)
1511 
1512  ro0 =roi(jj)
1513  !romn=roi(j-1)
1514  ropl=roi(jj+1)
1515 
1516  if(tetv.le.ttpl .and. tetv.ge.tt0) then
1517  zro=((ttpl-tetv)*ro0+(tetv-tt0)*ropl)/(ttpl-tt0)
1518  go to 1047
1519  endif
1520  enddo
1521 
1522  1047 continue
1523 
1524  ! zro=qvadin(teta(j), ttmn,tt0,ttpl, romn,ro0,ropl)
1525 
1526  rop(j)=zro
1527 
1528  roerr=dabs(ron(j)-ro(i,j))
1529 
1530  if(roerr.gt.erro) then
1531  erro=roerr
1532  ierm=i
1533  jerm=j
1534  endif
1535 
1536  roplt(i,j)=ron(j)-ro(i,j)
1537 
1538  r(i,j)=rop(j)*cos(teta(j))+rm
1539  z(i,j)=rop(j)*sin(teta(j))+zm
1540 
1541  enddo
1542 
1543  enddo ! i-loop
1544 
1545  if(kpr.eq.1) then
1546  write(*,*) 'erro max.',erro
1547  write(*,*) 'i,j erro',ierm,jerm,erro
1548  endif
1549 
1550  do i=2,iplas
1551  do j=2,nt1
1552 
1553  ro(i,j)=dsqrt((r(i,j)-rm)**2+(z(i,j)-zm)**2)
1554  ronor(i,j)=ro(i,j)/rob(j)
1555  psin(i,j)=psia(i)
1556 
1557  enddo
1558  enddo
1559 
1560 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1561 
1562  do j=2,nt1
1563 
1564  spro=rob(j)-ro(iplas-1,j)
1565 
1566  droj=ro(iplas,j)-ro(iplas-1,j)
1567  cpro=cpr(droj,npro,spro)
1568  row=ro(iplas,j)
1569  dro=droj
1570 
1571  do 220 i=iplas+1,nr1
1572 
1573  dro=dro*cpro
1574  row=row+dro
1575 
1576  rh(i)=rm+row*cos(teta(j))
1577  zh(i)=zm+row*sin(teta(j))
1578 
1579  roh(i)=row
1580 
1581  do is=2,nr1
1582 
1583  ps0 =psin(is,j)
1584  pspl=psin(is+1,j)
1585 
1586  ro0 =ro(is,j)
1587  ropl=ro(is+1,j)
1588 
1589  if(row.le.ropl .and. row.ge.ro0) then
1590  isw=is
1591  go to 755
1592  endif
1593 
1594  enddo
1595 
1596  755 continue
1597 
1598  ro2=ropl
1599  ro1=ro0
1600 
1601  ps2=pspl
1602  ps1=ps0
1603 
1604  psia(i)=( ps1*(ro2-row)+ps2*(row-ro1) )/(ro2-ro1)
1605 
1606  220 continue
1607 
1608  do 250 i=iplas+1,nr1
1609 
1610  roerr=dabs(roh(i)-ro(i,j))
1611  erro=dmax1(erro,roerr)
1612 
1613  ro(i,j)=roh(i)
1614 
1615  r(i,j)=roh(i)*cos(teta(j))+rm
1616  z(i,j)=roh(i)*sin(teta(j))+zm
1617 
1618  psin(i,j)=psia(i)
1619 
1620  250 continue
1621 
1622  enddo
1623 
1624  do j=2,nt1
1625  ro(nr,j)=rob(j)
1626  !r(nr,j)=rob(j)*cos(teta(j))+rm
1627  !z(nr,j)=rob(j)*sin(teta(j))+zm
1628  ro(1,j)=0.d0
1629  r(1,j)=rm
1630  z(1,j)=zm
1631  enddo
1632 
1633  do i=1,nr
1634  do j=2,nt1
1635 
1636  ronor(i,j)=ro(i,j)/ro(iplas,j)
1637 
1638  enddo
1639  enddo
1640 
1641 
1642  do 25 i=1,nr
1643 
1644  ro(i,1)=ro(i,nt1)
1645  ro(i,nt)=ro(i,2)
1646 
1647  ronor(i,1)=ronor(i,nt1)
1648  ronor(i,nt)=ronor(i,2)
1649 
1650  r(i,1)=r(i,nt1)
1651  r(i,nt)=r(i,2)
1652 
1653  z(i,1)=z(i,nt1)
1654  z(i,nt)=z(i,2)
1655 
1656  roplt(i,1)=roplt(i,nt1)
1657  roplt(i,nt)=roplt(i,2)
1658 
1659  25 continue
1660 
1661  do 30 i=1,iplas
1662  do 30 j=1,nt
1663 
1664  psin(i,j)=psia(i)
1665 
1666  30 continue
1667 
1668  do i=1,nr
1669  do j=1,nt
1670 
1671  psiold(i,j)=psip+psin(i,j)*(psim-psip)
1672 
1673  enddo
1674  enddo
1675 
1676  !f_call wrd
1677 
1678  return
1679  end
1680 
1681 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1682 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1683 
1684  subroutine f_regrid0(erro)
1685 
1686  include 'double.inc'
1687  include 'parevo.inc'
1688  parameter(nkp=njlim)
1689  include 'dim.inc'
1690  !!include 'comlmtr.inc'
1691  include 'compol.inc'
1692  include 'compol_add.inc'
1693 
1694  dimension ron(ntp),rop(ntp),rn(ntp),zn(ntp),tetn(ntp)
1695  dimension rob(ntp),teti(ntp),roi(ntp)
1696  dimension roplt(nrp,ntp)
1697  dimension dp(5)
1698  dimension roh(nrp),rh(nrp),zh(nrp)
1699  dimension psiold(nrp,ntp)
1700  common
1701  * /c_kpr/kpr
1702 
1703 ! save psiold
1704  real*8 psiold
1705 ! xzer(x1,x2,v1,v2) = (x1*v2-x2*v1)/(v2-v1)
1706 
1707  cos(xx)=dcos(xx)
1708  sin(xx)=dsin(xx)
1709  atan(xx)=datan(xx)
1710  asin(xx)=dasin(xx)
1711  acos(xx)=dacos(xx)
1712  sqrt(xx)=dsqrt(xx)
1713 
1714 
1715  nctrli=nctrl
1716 
1717  npro=nr-iplas+1
1718 
1719  if(ngav.eq.0) then
1720  alfa=1.0d0
1721  else
1722  alfa=0.75d0
1723  endif
1724 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1725  if(itin.gt.1) then
1726  wgt=0.5d0
1727  do i=1,iplas
1728  do j=1,nt
1729  !psi(i,j)=wgt*psi(i,j)+(1.d0-wgt)*psiold(i,j)
1730  enddo
1731  enddo
1732  endif
1733 
1734 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1735 
1736  rmold=rm
1737  zmold=zm
1738 
1739  if(kpr.eq.1) then
1740  write(*,*) 'rm,zm,psim',rm,zm,psim
1741  endif
1742 
1743  call axdef(rma,zma,psima,dp)
1744 
1745  if(kpr.eq.1) then
1746  write(*,*) 'rma,zma,psima',rma,zma,psima
1747  endif
1748 
1749  rm=rma
1750  zm=zma
1751  psim=psima
1752 
1753  erro=dsqrt((rmold-rm)**2+(zmold-zm)**2)
1754 
1755  if(kpr.eq.1) then
1756  write(*,*) 'erro mag.axis',erro
1757  endif
1758 
1759  errpsi=0.d0
1760 
1761 
1762 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1763 
1764  rxold=rx0
1765  zxold=zx0
1766 
1767  call f_xpoint(rx1,zx1,ixp1,jxp1,psix1,tetx1,kodex1)
1768 
1769  call f_xpoint(rx2,zx2,ixp2,jxp2,psix2,tetx2,kodex2)
1770 
1771  if(kpr.eq.1) then
1772  write(*,*) 'kodex1,kodex2',kodex1,kodex2
1773  endif
1774 
1775  kodxp=kodex1*kodex2
1776  if(kodxp.ne.0) then
1777 
1778  if(kpr.eq.1) then
1779  write(*,*) ' all x-points out of box'
1780  write(*,*) ' only limiter case '
1781  endif
1782 
1783  nctrli=1
1784 
1785  !write(6,*) ' program cannot run more '
1786  !stop
1787  psip=-1.d12
1788 
1789  else !kodxp=0
1790  psipn=psip
1791 
1792  if(kodex1.eq.0 .AND. kodex2.eq.0) then
1793  if(psix1.ge.psix2) then
1794  psix0=psix1
1795  rx0=rx1
1796  zx0=zx1
1797  else
1798  psix0=psix2
1799  rx0=rx2
1800  zx0=zx2
1801  endif
1802  elseif(kodex1.eq.0) then
1803  psix0=psix1
1804  rx0=rx1
1805  zx0=zx1
1806  elseif(kodex2.eq.0) then
1807  psix0=psix2
1808  rx0=rx2
1809  zx0=zx2
1810  endif
1811 
1812  errox=dsqrt((rxold-rx0)**2+(zxold-zx0)**2)
1813  if(kpr.eq.1) then
1814  write(*,*) 'erro x-point',errox
1815  write(*,*) 'rx zx',rx0,zx0,psix0
1816  endif
1817 
1818  psip=psim-alp*(psim-psix0)
1819 
1820  if(itin.gt.1) then
1821  zwgt=1.0d0
1822  psipr=zwgt*psip+(1.d0-zwgt)*psipn
1823  psip=dmax1(psip,psipr)
1824  alp_pr=(psim-psip)/(psim-psix0)
1825  endif
1826 
1827  endif
1828 
1829  if(nctrl.eq.1) then
1830  numlim=0
1831  do llim=1,nblm
1832 
1833  call fdefln(psiblm,rblm(llim),zblm(llim))
1834 
1835  if(psiblm.gt.psip) then
1836  psip=psiblm
1837  numlim=llim
1838  endif
1839 
1840  enddo
1841  endif
1842  if(kpr.eq.1) then
1843  write(*,*) 'numlim',numlim
1844  endif
1845 
1846 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1847 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1848 
1849 
1850  do i=1,nr
1851  do j=1,nt
1852 
1853  psnn=psin(i,j)
1854  psin(i,j)=(psi(i,j)-psip)/(psim-psip)
1855  delpsn=dabs(psin(i,j)-psnn)
1856 
1857  if(delpsn.gt.errpsi) then
1858 
1859  ierm=i
1860  jerm=j
1861  errpsi=delpsn
1862 
1863  endif
1864 
1865  enddo
1866  enddo
1867 
1868  if(kpr.eq.1) then
1869  write(*,*) 'i,j errpsi',ierm,jerm,errpsi,(psim-psip)
1870  endif
1871 
1872 cc definition of new angle grid teta(j)
1873 
1874  do j=1,nt
1875  tetn(j)=teta(j)
1876  enddo
1877 
1878  do j=1,nt
1879 
1880  drx=r(nr,j)-rm
1881  dzx=z(nr,j)-zm
1882 
1883  robj=sqrt(drx**2+dzx**2)
1884  rob(j)=robj
1885  tetp=acos(drx/robj)
1886  if(dzx.lt.0.d0) then
1887  teta(j)=-tetp
1888  else
1889  teta(j)=tetp
1890  endif
1891 
1892  enddo
1893 
1894  do j=2,nt
1895 
1896  if(teta(j).lt.teta(j-1)) then
1897  teta(j)=teta(j)+2.d0*pi
1898  endif
1899 
1900  if(teta(j).lt.teta(j-1)) then
1901  teta(j)=teta(j)+2.d0*pi
1902  endif
1903 
1904  enddo
1905 
1906  teta(1)=teta(nt1)-2.d0*pi
1907  teta(nt)=teta(2)+2.d0*pi
1908 
1909 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1910 
1911  i=2
1912 
1913  zps =psia(i)
1914  ps0=psip+zps*(psim-psip)
1915 
1916  do j=1,nt
1917 
1918  ttj=teta(j)
1919 
1920  ro2j=(ps0-psim)/( 0.5d0*dp(3)*cos(ttj)**2
1921  + + dp(4)*cos(ttj)*sin(ttj)
1922  + + 0.5d0*dp(5)*sin(ttj)**2 )
1923 
1924  rop(j)=dsqrt(ro2j)
1925 
1926  r(i,j)=rm+rop(j)*cos(teta(j))
1927  z(i,j)=zm+rop(j)*sin(teta(j))
1928 
1929  enddo ! j-loop
1930 
1931 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1932 
1933 ccc grid moving along rais
1934 
1935  do i=3,iplas
1936 
1937  zps =psia(i)
1938 
1939  do j=1,nt
1940 
1941  do is=1,nr1
1942 
1943  ps0 =psin(is,j)
1944 ! !psmn=psin(is-1,j)
1945  pspl=psin(is+1,j)
1946 
1947  ro0 =ro(is,j)
1948 ! !romn=ro(is-1,j)
1949  ropl=ro(is+1,j)
1950 
1951  if(zps.le.ps0 .and. zps.ge.pspl) then
1952  zro=((pspl-zps)*ro0-(ps0-zps)*ropl)/(pspl-ps0)
1953  go to 822
1954  endif
1955 ! !zro=qvadin(zps, psmn,ps0,pspl, romn,ro0,ropl)
1956  enddo
1957  822 continue
1958 
1959 
1960 
1961  ron(j)=zro
1962 
1963  enddo
1964 
1965  do j=1,nt
1966 
1967  rnj=rmold+ron(j)*cos(tetn(j))
1968  znj=zmold+ron(j)*sin(tetn(j))
1969 
1970  drx=rnj-rm
1971  dzx=znj-zm
1972 
1973  roi(j)=sqrt(drx**2+dzx**2)
1974  tetp=acos(drx/roi(j))
1975  if(dzx.lt.0.d0) then
1976  teti(j)=-tetp
1977  else
1978  teti(j)=tetp
1979  endif
1980 
1981  enddo
1982 
1983  do j=2,nt
1984 
1985  if(teti(j).lt.teti(j-1)) then
1986  teti(j)=teti(j)+2.d0*pi
1987  endif
1988 
1989  if(teti(j).lt.teti(j-1)) then
1990  teti(j)=teti(j)+2.d0*pi
1991  endif
1992 
1993  enddo
1994 
1995  teti(1)=teti(nt1)-2.d0*pi
1996  teti(nt)=teti(2)+2.d0*pi
1997 
1998  roi(1)=roi(nt1)
1999  roi(nt)=roi(2)
2000 
2001  do j=2,nt1
2002 
2003  tetv=teta(j)
2004 
2005  if(tetv.lt.teti(1)) tetv=tetv+2.d0*pi
2006  if(tetv.lt.teti(1)) tetv=tetv+2.d0*pi
2007  if(tetv.gt.teti(nt)) tetv=tetv-2.d0*pi
2008  if(tetv.gt.teti(nt)) tetv=tetv-2.d0*pi
2009 
2010  do jj=1,nt1
2011 
2012  tt0 =teti(jj)
2013  !ttmn=teti(j-1)
2014  ttpl=teti(jj+1)
2015 
2016  ro0 =roi(jj)
2017  !romn=roi(j-1)
2018  ropl=roi(jj+1)
2019 
2020  if(tetv.le.ttpl .and. tetv.ge.tt0) then
2021  zro=((ttpl-tetv)*ro0+(tetv-tt0)*ropl)/(ttpl-tt0)
2022  go to 1047
2023  endif
2024  enddo
2025 
2026  1047 continue
2027 
2028  ! zro=qvadin(teta(j), ttmn,tt0,ttpl, romn,ro0,ropl)
2029 
2030  rop(j)=zro
2031 
2032  roerr=dabs(ron(j)-ro(i,j))
2033 
2034  if(roerr.gt.erro) then
2035  erro=roerr
2036  ierm=i
2037  jerm=j
2038  endif
2039 
2040  roplt(i,j)=ron(j)-ro(i,j)
2041 
2042  r(i,j)=rop(j)*cos(teta(j))+rm
2043  z(i,j)=rop(j)*sin(teta(j))+zm
2044 
2045  enddo
2046 
2047  enddo ! i-loop
2048 
2049  if(kpr.eq.1) then
2050  write(6,*) 'erro max.',erro
2051  write(6,*) 'i,j erro',ierm,jerm,erro
2052  endif
2053 
2054  do i=2,iplas
2055  do j=2,nt1
2056 
2057  ro(i,j)=dsqrt((r(i,j)-rm)**2+(z(i,j)-zm)**2)
2058  ronor(i,j)=ro(i,j)/rob(j)
2059  psin(i,j)=psia(i)
2060 
2061  enddo
2062  enddo
2063 
2064 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2065 
2066  do j=2,nt1
2067 
2068  spro=rob(j)-ro(iplas-1,j)
2069 
2070  droj=ro(iplas,j)-ro(iplas-1,j)
2071  !cpro=cpr(droj,npro,spro)
2072  row=ro(iplas,j)
2073  dro=(rob(j)-ro(iplas,j))/(nr-iplas)
2074 
2075  do 220 i=iplas+1,nr1
2076 
2077  !dro=dro*cpro
2078  row=row+dro
2079 
2080  rh(i)=rm+row*cos(teta(j))
2081  zh(i)=zm+row*sin(teta(j))
2082 
2083  roh(i)=row
2084 
2085  do is=2,nr1
2086 
2087  ps0 =psin(is,j)
2088  pspl=psin(is+1,j)
2089 
2090  ro0 =ro(is,j)
2091  ropl=ro(is+1,j)
2092 
2093  if(row.le.ropl .and. row.ge.ro0) then
2094  isw=is
2095  go to 755
2096  endif
2097 
2098  enddo
2099 
2100  755 continue
2101 
2102  ro2=ropl
2103  ro1=ro0
2104 
2105  ps2=pspl
2106  ps1=ps0
2107 
2108  psia(i)=( ps1*(ro2-row)+ps2*(row-ro1) )/(ro2-ro1)
2109 
2110  220 continue
2111 
2112  do 250 i=iplas+1,nr1
2113 
2114  roerr=dabs(roh(i)-ro(i,j))
2115  erro=dmax1(erro,roerr)
2116 
2117  ro(i,j)=roh(i)
2118 
2119  r(i,j)=roh(i)*cos(teta(j))+rm
2120  z(i,j)=roh(i)*sin(teta(j))+zm
2121 
2122  psin(i,j)=psia(i)
2123 
2124  250 continue
2125 
2126  enddo
2127 
2128  do j=2,nt1
2129  ro(nr,j)=rob(j)
2130  !r(nr,j)=rob(j)*cos(teta(j))+rm
2131  !z(nr,j)=rob(j)*sin(teta(j))+zm
2132  ro(1,j)=0.d0
2133  r(1,j)=rm
2134  z(1,j)=zm
2135  enddo
2136 
2137  do i=1,nr
2138  do j=2,nt1
2139 
2140  ronor(i,j)=ro(i,j)/ro(iplas,j)
2141 
2142  enddo
2143  enddo
2144 
2145 
2146  do 25 i=1,nr
2147 
2148  ro(i,1)=ro(i,nt1)
2149  ro(i,nt)=ro(i,2)
2150 
2151  ronor(i,1)=ronor(i,nt1)
2152  ronor(i,nt)=ronor(i,2)
2153 
2154  r(i,1)=r(i,nt1)
2155  r(i,nt)=r(i,2)
2156 
2157  z(i,1)=z(i,nt1)
2158  z(i,nt)=z(i,2)
2159 
2160  roplt(i,1)=roplt(i,nt1)
2161  roplt(i,nt)=roplt(i,2)
2162 
2163  25 continue
2164 
2165  do 30 i=1,iplas
2166  do 30 j=1,nt
2167 
2168  psin(i,j)=psia(i)
2169 
2170  30 continue
2171 
2172  do i=1,nr
2173  do j=1,nt
2174 
2175  psiold(i,j)=psip+psin(i,j)*(psim-psip)
2176 
2177  enddo
2178  enddo
2179 
2180  !f_call wrd
2181 
2182  return
2183  end
2184 
2185 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2186 
2187  subroutine grid_spdr
2188 
2189  include 'double.inc'
2190 
2191  include 'parevo.inc'
2192  parameter(nkp=njlim)
2193  include 'dim.inc'
2194  include 'compol.inc'
2195  include 'parrc1.inc'
2196  include 'comrec.inc'
2197  include 'compol_add.inc'
2198  common
2199  * /c_kpr/kpr
2200 
2201  cos(xx)=dcos(xx)
2202  sin(xx)=dsin(xx)
2203  sqrt(xx)=dsqrt(xx)
2204 
2205 
2206  !call rdfb
2207 
2208  rm=r(1,2)
2209  zm=z(1,2)
2210 
2211  do i=iplas,1,-1
2212  if(psia(i).ge.0.05d0) then
2213  i_bc=i
2214  exit
2215  endif
2216  enddo
2217 
2218  do j=2,nt1
2219  ronor(nr,j)=1.57d0
2220  ro(nr,j)=ro(i_bc,j)*ronor(nr,j)
2221  enddo
2222 
2223  !do i=2,iplas
2224  !do j=2,nt1
2225  ! ro(i,j)=ro(iplas,j)*ronor(i,j)
2226  !enddo
2227  !enddo
2228 
2229  do i=iplas+1,nr-1
2230  do j=2,nt1
2231  ropla=ro(iplas,j)
2232  robon=ro(nr,j)
2233  delro=(robon-ropla)/(nr-iplas)
2234  ro(i,j)=ropla+delro*(i-iplas)
2235  ronor(i,j)=ro(i,j)/ropla
2236  enddo
2237  enddo
2238 
2239  do i=iplas+1,nr
2240  do j=2,nt1
2241  r(i,j)=ro(i,j)*cos(teta(j))+rm
2242  z(i,j)=ro(i,j)*sin(teta(j))+zm
2243  enddo
2244  enddo
2245 
2246  do i=1,nr
2247 
2248  ro(i,1)=ro(i,nt1)
2249  ro(i,nt)=ro(i,2)
2250 
2251  ronor(i,1)=ronor(i,nt1)
2252  ronor(i,nt)=ronor(i,2)
2253 
2254  r(i,1)=r(i,nt1)
2255  r(i,nt)=r(i,2)
2256 
2257  z(i,1)=z(i,nt1)
2258  z(i,nt)=z(i,2)
2259 
2260  enddo
2261 
2262 
2263  rb_max=rm
2264  rb_min=rm
2265  zb_max=zm
2266  zb_min=zm
2267 
2268  do j=2,nt1
2269  if(r(nr,j).gt.rb_max) rb_max=r(nr,j)
2270  if(r(nr,j).lt.rb_min) rb_min=r(nr,j)
2271  if(z(nr,j).gt.zb_max) zb_max=z(nr,j)
2272  if(z(nr,j).lt.zb_min) zb_min=z(nr,j)
2273  enddo
2274 
2275  if(rb_max.gt.xmax
2276  * .OR.
2277  * zb_max.gt.ymax
2278  * .OR.
2279  * rb_min.lt.xmin
2280  * .OR.
2281  * zb_min.lt.ymin ) then
2282 
2283  if(kpr.eq.1) then
2284  write(*,*) 'egg is large then rectangular box'
2285  call f_wrd
2286  pause 'pause '
2287  endif
2288 
2289  endif
2290 
2291  rx1=xx1 ! x-points from rect. grid
2292  zx1=yx1 !
2293  !
2294  rx2=xx2 !
2295  zx2=yx2 !
2296  !
2297  rx0=xx0 !
2298  zx0=yx0 !
2299 
2300  ixp1=0
2301  jxp1=0
2302 
2303  ixp2=0
2304  jxp2=0
2305 
2306  psim=um
2307  psip=up
2308  ! call f_wrd
2309 
2310  return
2311  end
2312 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2313 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2314 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2315 
2316  subroutine f_grid(igdf,nstep)
2317 
2318  include 'double.inc'
2319 
2320  include 'parevo.inc'
2321  parameter(nkp=njlim)
2322  include 'dim.inc'
2323  parameter(nshp=10)
2324  parameter(ntp4=ntp+4,ntp6=ntp4*6)
2325  include 'parrc1.inc'
2326 
2327  include 'compol.inc'
2328  include 'compol_add.inc'
2329  include 'comrec.inc'
2330  common /com_sp3/ rrksp(ntp4),ccksp(ntp4)
2331 
2332  dimension xs(nshp),ys(nshp),fun(nshp),dp(5)
2333  real*8 ron(ntp)
2334  common
2335  * /c_kpr/kpr
2336 
2337  frbon(r0,as,tr,tet,skv)=
2338  * r0+(r0/as)*dcos(tet+tr*dsin(tet)-skv*dsin(2.d0*tet))
2339  fzbon(r0,z0,as,el,tet)=z0+(r0/as)*el*dsin(tet)
2340 
2341  cos(xx)=dcos(xx)
2342  sin(xx)=dsin(xx)
2343  sqrt(xx)=dsqrt(xx)
2344 
2345  if(nstep.eq.0) then
2346  if(igdf.eq.0)then
2347 
2348  do i=1,iplas
2349  psia(i)=(iplas-i)/(iplas-1.d0)
2350  enddo
2351 
2352  do i=1,iplas1
2353  dpsda(i)=-1.d0
2354  enddo
2355 
2356  elseif(igdf.eq.1 .OR. igdf.eq.2)then
2357 
2358  do i=1,iplas
2359  psia(i)=1.d0-((i-1)/(iplas-1.d0))**2
2360  enddo
2361 
2362  do i=1,iplas1
2363  dpsda(i)=(1-2*i)/(iplas-1.d0)
2364  enddo
2365 
2366  endif
2367  endif
2368 
2369  dtet=2.d0*pi/(nt-2)
2370 
2371  teta(1)=-dtet
2372 
2373  do 30 j=2,nt
2374 
2375  teta(j)=teta(j-1)+dtet
2376 
2377  30 continue
2378 
2379  teta(1)=teta(nt1)-2.d0*pi
2380  teta(nt)=teta(2)+2.d0*pi
2381 
2382  rm=xm !
2383  zm=ym ! position of magn. axis
2384  !
2385  rx1=xx1 ! and x-points from rect. grid
2386  zx1=yx1 !
2387  !
2388  rx2=xx2 !
2389  zx2=yx2 !
2390  !
2391  rx0=xx0 !
2392  zx0=yx0 !
2393 
2394  ixp1=0
2395  jxp1=0
2396 
2397  ixp2=0
2398  jxp2=0
2399 
2400 
2401 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2402 
2403 
2404  write(fname,'(a,a)') path(1:kname),'shegg.dat'
2405  open(1,file=fname)
2406  !open(1,file='shegg.dat')
2407 
2408  read(1,*) rc0
2409  read(1,*) zc0
2410 
2411  read(1,*) asp0
2412 
2413  read(1,*) el_up
2414  read(1,*) tr_up
2415 
2416  read(1,*) el_lw
2417  read(1,*) tr_lw
2418 
2419  read(1,*) skv
2420 
2421  close(1)
2422 
2423 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2424 
2425  do j=1,nt
2426  tet=teta(j)
2427  tri=0.5d0*(tr_up+tr_lw+(tr_up-tr_lw)*sin(tet))
2428  ell=0.5d0*(el_up+el_lw+(el_up-el_lw)*sin(tet))
2429  rrr=frbon(rc0,asp0,tri,tet,skv)
2430  r(nr,j)=rrr
2431  zzz=fzbon(rc0,zc0,asp0,ell,tet)
2432  z(nr,j)=zzz
2433  enddo
2434 
2435 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2436 
2437  do j=1,nt
2438  drx=r(nr,j)-rm
2439  dzx=z(nr,j)-zm
2440  teta(j)=atan(dzx/drx)
2441  enddo
2442 
2443  do j=2,nt
2444 
2445  if(teta(j).lt.teta(j-1)) then
2446  teta(j)=teta(j)+pi
2447  endif
2448 
2449  if(teta(j).lt.teta(j-1)) then
2450  teta(j)=teta(j)+pi
2451  endif
2452 
2453  enddo
2454 
2455  teta(1)=teta(nt1)-2.d0*pi
2456  teta(nt)=teta(2)+2.d0*pi
2457 
2458 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2459 
2460 
2461  nsh=1
2462  xs(nsh)=x(imax)
2463  ys(nsh)=y(jmax)
2464  fun(nsh)=u(imax,jmax)
2465 
2466  do k=-1,1
2467 
2468  ii= imax+k
2469 
2470  do l=-1,1
2471 
2472  jj= jmax+l
2473 
2474  if(ii.ne.imax .OR. jj.ne.jmax) then
2475  nsh=nsh+1
2476  xs(nsh)=x(ii)
2477  ys(nsh)=y(jj)
2478  fun(nsh)=u(ii,jj)
2479  endif
2480 
2481  enddo
2482  enddo
2483 
2484  call deriv5(xs,ys,fun,nsh,5,dp)
2485 
2486  !write(6,*) dp
2487 
2488 
2489  stpx=(x(imax+1)-x(imax))
2490  stpy=(y(jmax+1)-y(jmax))
2491 
2492  do i=2,iplas
2493 
2494  u0=psia(i)
2495  ur0=up+u0*(um-up)
2496 
2497  do j=1,nt
2498 
2499  ttj=teta(j)
2500  ro2j=(ur0-um)/( 0.5d0*dp(3)*cos(ttj)**2
2501  + + dp(4)*cos(ttj)*sin(ttj)
2502  + + 0.5d0*dp(5)*sin(ttj)**2 )
2503 
2504  ron(j)=sqrt(ro2j)
2505 
2506 
2507  r(i,j)=rm+ron(j)*cos(teta(j))
2508  z(i,j)=zm+ron(j)*sin(teta(j))
2509 
2510  ro(i,j)=ron(j)
2511 
2512  psin(i,j)=u0
2513 
2514  enddo ! j-loop
2515 
2516  ibeg=i+1
2517 
2518  if(kpr.eq.1) then
2519  write(6,*) 'ro(1)',ro(i,1),ro(i,2),stpx,stpy
2520  endif
2521 
2522  if( ron(2) .GT. dmax1(stpx,stpy) ) go to 2799
2523 
2524  enddo !i-loop
2525 
2526  2799 continue
2527 
2528 c pause 'pausew'
2529 
2530  do 100 i=ibeg,iplas
2531 
2532  u0=psia(i)
2533 
2534  call f_loop95(teta,nt,ron,u0)
2535 
2536  ! write(6,*) 'loop95 i=',i,u0
2537 
2538  do 110 j=1,nt
2539 
2540  r(i,j)=rm+ron(j)*cos(teta(j))
2541  z(i,j)=zm+ron(j)*sin(teta(j))
2542 
2543  ro(i,j)=ron(j)
2544 
2545  psin(i,j)=u0
2546 
2547  110 continue
2548  100 continue
2549 
2550  do 120 j=1,nt
2551 
2552  r(1,j)=rm
2553  z(1,j)=zm
2554  ro(1,j)=0.
2555  psin(1,j)=1.d0
2556 
2557  120 continue
2558 
2559  npro=nr-iplas+1
2560 
2561  do 200 j=1,nt
2562 
2563  droj=sqrt( (r(iplas,j)-r(iplas-1,j))**2+
2564  + (z(iplas,j)-z(iplas-1,j))**2 )
2565  delbn=sqrt( (r(nr,j)-r(iplas,j))**2+
2566  + (z(nr,j)-z(iplas,j))**2 )
2567 
2568  cpro=cpr(droj,npro,delbn+droj)
2569  row=ron(j)
2570  dro=droj
2571 
2572  do 220 i=iplas+1,nr
2573 
2574  dro=dro*cpro
2575  row=row+dro
2576 
2577  r(i,j)=rm+row*cos(teta(j))
2578  z(i,j)=zm+row*sin(teta(j))
2579 
2580  ro(i,j)=row
2581 
2582  psin(i,j)=0.d0
2583 
2584  220 continue
2585  200 continue
2586 
2587 
2588  do 333 i=1,nr
2589 
2590  r(i,1)=r(i,nt1)
2591  z(i,1)=z(i,nt1)
2592  ro(i,1)=ro(i,nt1)
2593 
2594  r(i,nt)=r(i,2)
2595  z(i,nt)=z(i,2)
2596  ro(i,nt)=ro(i,2)
2597 
2598  333 continue
2599 
2600  teta(1)=teta(nt1)-2.d0*pi
2601  teta(nt)=teta(2)+2.d0*pi
2602 
2603  psip=up
2604 
2605  call bonspl
2606  !pause 'bonspl'
2607 
2608  call f_wrd
2609  !stop
2610 c pause 'wrd'
2611 
2612  return
2613  end
2614 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2615 
2616  subroutine ada(erro)
2617 
2618  include 'double.inc'
2619  include 'parevo.inc'
2620  parameter(nkp=njlim)
2621  include 'dim.inc'
2622  parameter(ntp4=ntp+4,ntp6=ntp4*6)
2623  !!include 'comlmtr.inc'
2624  common /com_sp3/ rrksp(ntp4),ccksp(ntp4)
2625  include 'compol.inc'
2626  include 'compol_add.inc'
2627 
2628  dimension ron(nrp),rn(nrp),zn(nrp)
2629  real*8 cwksp(4)
2630  common
2631  * /c_kpr/kpr
2632 
2633  xzer(x1,x2,v1,v2) = (x1*v2-x2*v1)/(v2-v1)
2634 
2635  cos(xx)=dcos(xx)
2636  sin(xx)=dsin(xx)
2637 
2638  nctrli=nctrl
2639  npro=nr-iplas+1
2640  erro=0.d0
2641 
2642  psim=psi(1,1)
2643 
2644  call f_xpoint(rx1,zx1,ixp1,jxp1,psix1,tetx1,kodex1)
2645 
2646  if(kpr.eq.1) then
2647  write(*,*) 'kodex1',kodex1
2648  write(*,*) 'rx1 zx1',rx1,zx1,psix1
2649  endif
2650 
2651  call f_xpoint(rx2,zx2,ixp2,jxp2,psix2,tetx2,kodex2)
2652 
2653  if(kpr.eq.1) then
2654  write(*,*) 'kodex2',kodex2
2655  write(*,*) 'rx2 zx2',rx2,zx2,psix2
2656  endif
2657 
2658  kodxp=kodex1*kodex2
2659  if(kodxp.ne.0) then
2660 
2661  if(kpr.eq.1) then
2662  write(*,*) ' all x-points out of box'
2663  write(*,*) ' only limiter case '
2664  endif
2665 
2666  nctrli=1
2667 
2668  !write(*,*) ' program cannot run more '
2669  !stop
2670  psip=-1.d12
2671 
2672  else !kodxp=0
2673 
2674  if(kodex1.eq.0 .AND. kodex2.eq.0) then
2675  if(psix1.gt.psix2) then
2676  psix0=psix1
2677  rx0=rx1
2678  zx0=zx1
2679  else
2680  psix0=psix2
2681  rx0=rx2
2682  zx0=zx2
2683  endif
2684  elseif(kodex1.eq.0) then
2685  psix0=psix1
2686  rx0=rx1
2687  zx0=zx1
2688  elseif(kodex2.eq.0) then
2689  psix0=psix2
2690  rx0=rx2
2691  zx0=zx2
2692  endif
2693 
2694  if(kpr.eq.1) then
2695  write(*,*) 'kodex1,kodex2',kodex1,kodex2
2696  write(*,*) 'rx zx',rx0,zx0,psix0
2697  endif
2698 
2699 
2700  if(ngav/10*10.eq.ngav) then
2701  psip=psim-alp*(psim-psix0)
2702  else
2703  psip=psi(iplas,jrolim)
2704  endif
2705 
2706  endif
2707 
2708  numlim=0
2709  if(nctrl.eq.1) then
2710  do llim=1,nblm
2711  call fdefln(psiblm,rblm(llim),zblm(llim))
2712  if(psiblm.gt.psip ) then
2713  psip=psiblm
2714  numlim=llim
2715  endif
2716  enddo
2717  endif
2718 
2719  alpnew=(psim-psip)/(psim-psix0)
2720 
2721  do i=1,nr
2722  do j=1,nt
2723 
2724  psin(i,j)=(psi(i,j)-psip)/(psim-psip)
2725 
2726  enddo
2727  enddo
2728 
2729  write(*,*)'ada:'
2730  write(*,*)'ada:nctrli,nctrl',nctrli,nctrl
2731  write(*,*)'ada:alpnew,numlim',alpnew,numlim
2732 
2733  ron(1)=0.d0
2734 
2735  do 20 j=1,nt
2736 
2737  do 200 i=2,iplas
2738 
2739  zps =psia(i)
2740 
2741  do is=1,nr1
2742 
2743  ps0 =psin(is,j)
2744  pspl=psin(is+1,j)
2745 
2746  ro0 =ro(is,j)
2747  ropl=ro(is+1,j)
2748 
2749  if(zps.le.ps0 .and. zps.ge.pspl) then
2750  isc=is
2751  go to 822
2752  endif
2753 
2754  enddo
2755 
2756  822 continue
2757  if(isc.eq.1) then
2758  grad=-(pspl-ps0)/(ropl-ro0)
2759  else
2760  psmn=psin(isc-1,j)
2761  romn=ro(isc-1,j)
2762  gradpl=-(pspl-ps0)/(ropl-ro0)
2763  gradmn=-(ps0-psmn)/(ro0-romn)
2764  grad=dmax1(gradpl,gradmn)
2765  endif
2766 
2767  grad=-(pspl-ps0)/(ropl-ro0)
2768 
2769  zro=ro0-(zps-ps0)/grad
2770  ron(i)=zro
2771 
2772  !ps0 =psin(i,j)
2773  !psmn=psin(i-1,j)
2774  ! pspl=psin(i+1,j)
2775 
2776  ! ro0 =ro(i,j)
2777  ! romn=ro(i-1,j)
2778  ! ropl=ro(i+1,j)
2779 
2780  ! zro=qvadin(zps, psmn,ps0,pspl, romn,ro0,ropl)
2781 
2782  if(ngav.eq.0) then
2783  alfa=0.75d0
2784  else
2785  alfa=0.550d0
2786  endif
2787  ron(i)=alfa*zro+(1.d0-alfa)*ro(i,j)
2788 
2789 
2790  if(ron(i).lt.ron(i-1)) then
2791  ron(i)=ron(i-1)+1.d-8
2792  if(kpr.eq.1) then
2793  write(*,*) 'grid crash i,j',i,j
2794  pause ' '
2795  endif
2796  endif
2797 
2798  200 continue
2799 
2800  if(ngav.gt.0 .AnD. j.eq.jrolim .AnD. kpr.eq.1) then
2801  write(*,*) 'ro(iplas),rolim',ron(iplas),rolim
2802  endif
2803 
2804  spro=ro(nr,j)-ron(iplas-1)
2805 
2806  droj=ron(iplas)-ron(iplas-1)
2807  !cpro=cpr(droj,npro,spro)
2808  row=ron(iplas)
2809  !dro=droj
2810  dro=(ro(nr,j)-row)/(nr-iplas)
2811 
2812  do 220 i=iplas+1,nr1
2813 
2814  !dro=dro*cpro
2815  row=row+dro
2816 
2817  rn(i)=rm+row*cos(teta(j))
2818  zn(i)=zm+row*sin(teta(j))
2819 
2820  ron(i)=row
2821 
2822 
2823 
2824  do is=isc,nr1
2825 
2826  ps0 =psin(is,j)
2827  pspl=psin(is+1,j)
2828 
2829  ro0 =ro(is,j)
2830  ropl=ro(is+1,j)
2831 
2832  if(row.le.ropl .and. row.ge.ro0) then
2833  isw=is
2834  go to 755
2835  endif
2836 
2837  enddo
2838 
2839  755 continue
2840 
2841  ro2=ropl
2842  ro1=ro0
2843 
2844  ps2=pspl
2845  ps1=ps0
2846 
2847 
2848  psia(i)=( ps1*(ro2-row)+ps2*(row-ro1) )/(ro2-ro1)
2849 
2850  220 continue
2851 
2852  do 250 i=2,nr1
2853 
2854  roerr=dabs(ron(i)-ro(i,j))
2855  erro=dmax1(erro,roerr)
2856 
2857  ro(i,j)=ron(i)
2858 
2859  r(i,j)=ron(i)*cos(teta(j))+rm
2860  z(i,j)=ron(i)*sin(teta(j))+zm
2861 
2862  psin(i,j)=psia(i)
2863 
2864  250 continue
2865 
2866  20 continue
2867 
2868 
2869 
2870 
2871  teta(1)=teta(nt1)-2.d0*pi
2872  teta(nt)=teta(2)+2.d0*pi
2873 
2874  do i=1,nr
2875  ro(i,1)=ro(i,nt1)
2876  ro(i,nt)=ro(i,2)
2877  enddo
2878 
2879  do i=1,nr
2880  do j=2,nt1
2881 
2882  ronor(i,j)=ro(i,j)/ro(iplas,j)
2883 
2884  enddo
2885  enddo
2886 
2887 
2888  do 25 i=1,nr
2889 
2890  ronor(i,1)=ronor(i,nt1)
2891  ronor(i,nt)=ronor(i,2)
2892 
2893  r(i,1)=r(i,nt1)
2894  r(i,nt)=r(i,2)
2895 
2896  z(i,1)=z(i,nt1)
2897  z(i,nt)=z(i,2)
2898 
2899  psin(i,1)=psin(i,nt1)
2900  psin(i,nt)=psin(i,2)
2901 
2902  25 continue
2903 
2904 
2905  do 30 i=1,nr
2906  do 30 j=1,nt
2907 
2908  psi(i,j)=psip+psin(i,j)*(psim-psip)
2909 
2910  30 continue
2911 
2912  !call f_wrd
2913  ! stop
2914 
2915  return
2916  end
2917 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2918 
2919  subroutine reform
2920 
2921  include 'double.inc'
2922  include 'parevo.inc'
2923  parameter(nkp=njlim)
2924  include 'dim.inc'
2925  include 'compol.inc'
2926  include 'compol_add.inc'
2927 
2928  dimension rob(ntp)
2929 
2930  sqrt(xx)=dsqrt(xx)
2931  atan(xx)=datan(xx)
2932  cos(xx)=dcos(xx)
2933  sin(xx)=dsin(xx)
2934 
2935  do 100 j=1,nt
2936 
2937  drx=r(nr,j)-rm
2938  dzx=z(nr,j)-zm
2939 
2940  rob(j)=sqrt(drx**2+dzx**2)
2941 
2942 
2943  tetp=dacos(drx/rob(j))
2944  if(dzx.lt.0.d0) then
2945  tet0=2.d0*pi-tetp
2946  else
2947  tet0=tetp
2948  endif
2949 
2950  teta(j)=tet0
2951 
2952  100 continue
2953 
2954  do 110 j=2,nt
2955 
2956  if(teta(j).lt.teta(j-1)) then
2957  teta(j)=teta(j)+pi
2958  endif
2959 
2960  if(teta(j).lt.teta(j-1)) then
2961  teta(j)=teta(j)+pi
2962  endif
2963  110 continue
2964 
2965  teta(1)=teta(nt1)-2.d0*pi
2966  teta(nt)=teta(2)+2.d0*pi
2967 
2968 
2969  do 200 j=1,nt
2970  do 200 i=1,nr
2971 
2972  ro(i,j)=ro(i,j)*rob(j)/ro(nr,j)
2973  r(i,j)=rm+ro(i,j)*cos(teta(j))
2974  z(i,j)=zm+ro(i,j)*sin(teta(j))
2975 
2976  200 continue
2977 
2978  if(ngav/10*10.ne.ngav) then
2979  ro(iplas,jrolim)=rolim
2980  r(iplas,jrolim)=rm+ro(iplas,jrolim)*cos(teta(jrolim))
2981  z(iplas,jrolim)=zm+ro(iplas,jrolim)*sin(teta(jrolim))
2982  endif
2983 
2984  return
2985  end
2986 
2987 
2988 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2989 
2990  subroutine fdefln(psidf,rgv,zgv)
2991 
2992  include 'double.inc'
2993  include 'parevo.inc'
2994  parameter(nkp=njlim)
2995  include 'dim.inc'
2996  include 'compol.inc'
2997  include 'compol_add.inc'
2998 
2999  sqrt(arg)=dsqrt(arg)
3000 
3001  psidf=0.d0
3002 
3003  dr0=rgv-rm
3004  dz0=zgv-zm
3005 
3006  ro0=sqrt(dr0**2+dz0**2)
3007 
3008  tetp=dacos(dr0/ro0)
3009  if(dz0.lt.0.d0) then
3010  tet0=2.d0*pi-tetp
3011  else
3012  tet0=tetp
3013  endif
3014 
3015  if(tet0.lt.teta(1)) tet0=tet0+2.d0*pi
3016  if(tet0.gt.teta(nt)) tet0=tet0-2.d0*pi
3017 
3018  do j=1,nt1
3019 
3020  if(tet0.ge.teta(j) .AND. tet0.lt.teta(j+1)) jc=j
3021 
3022  enddo
3023 
3024  rob1=ro(nr,jc)
3025  rob2=ro(nr,jc+1)
3026 
3027  ropl1=ro(iplas,jc)
3028  ropl2=ro(iplas,jc+1)
3029 
3030  tet1=teta(jc)
3031  tet2=teta(jc+1)
3032 
3033  rob12=(rob1*(tet2-tet0)+rob2*(tet0-tet1))/(tet2-tet1)
3034  ropl12=(ropl1*(tet2-tet0)+ropl2*(tet0-tet1))/(tet2-tet1)
3035 
3036  if(ro0.lt.rob12) then !!!
3037 
3038  do i=nr,2,-1
3039 
3040  ro1=ro(i,jc)
3041  ro2=ro(i,jc+1)
3042 
3043  ro12=(ro1*(tet2-tet0)+ro2*(tet0-tet1))/(tet2-tet1)
3044 
3045  if(ro0.lt.ro12) ic=i-1
3046 
3047  enddo
3048 
3049  i_rlmp=1
3050 
3051  if(ic.gt.iplas) then
3052  do i=ic+1,iplas-1,-1
3053  der_psi=psi(i,jc)-psi(i-1,jc)
3054  if(der_psi.gt.0.d0) then
3055  i_rlmp=0
3056  psidf=-1.0d12
3057  endif
3058  enddo
3059  endif
3060 
3061  if(i_rlmp.eq.1) then
3062 
3063  ro1=ro(ic,jc)
3064  ro2=ro(ic+1,jc)
3065  ro3=ro(ic+1,jc+1)
3066  ro4=ro(ic,jc+1)
3067 
3068  u1=psi(ic,jc)
3069  u2=psi(ic+1,jc)
3070  u3=psi(ic+1,jc+1)
3071  u4=psi(ic,jc+1)
3072 
3073  psidf=blin_tr(tet0,ro0,
3074  * tet1,tet2,ro1,ro2,ro3,ro4,u1,u2,u3,u4)
3075  endif
3076 
3077  else !!!!!!!!!!!!!!!!!!!!!
3078 
3079  psidf=-1.0d12
3080  return
3081 
3082  rr=rgv
3083  zz=zgv
3084 
3085  do 110 j=2,nt1
3086 
3087  r0=r(nr,j)
3088  z0=z(nr,j)
3089 
3090  r1=r(nr,j+1)
3091  z1=z(nr,j+1)
3092 
3093  call bint(rr,zz,r0,z0,r1,z1,fint,1)
3094 
3095  psidf=psidf-fint*(dgdn(j)+dgdn(j+1))*0.5d0
3096 
3097  !write(6,*) 'flux:'
3098  ! write(6,*) 'fint j:',fint
3099  ! write(6,*) 'dgdn:',dgdn(j)
3100  !write(6,*) 'psitok:',psitok(ik)
3101 
3102  110 continue
3103 
3104  endif !!!!!!!!!!!!!!!!!!!!
3105 
3106  return
3107  end
3108 
3109 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3110 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3111  subroutine bonspl
3112 
3113  include 'double.inc'
3114  include 'parevo.inc'
3115  parameter(nkp=njlim)
3116  include 'dim.inc'
3117  parameter(ntp4=ntp+4,ntp6=ntp4*6)
3118  common /com_sp3/ rrksp(ntp4),ccksp(ntp4)
3119  include 'compol.inc'
3120  include 'compol_add.inc'
3121  real*8 wrk(ntp6)
3122  real*8 robb(ntp)
3123 
3124  !pause 'bonspl enter'
3125 
3126  do j=1,nt
3127  robb(j)=ro(nr,j)
3128  enddo
3129 
3130  CALL e01baf(nt,teta,robb,rrksp,ccksp,nt+4,wrk,6*nt+16,ifail)
3131 
3132  return
3133  end
3134 
3135 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3136  subroutine f_loop95(tetpol,ntet,ro0,u0)
3137 
3138  include 'double.inc'
3139  include 'parrc1.inc'
3140  include 'comrec.inc'
3141 
3142  dimension rxb(nbndp2),zxb(nbndp2)
3143 
3144  real*8 ut(nip,njp)
3145 
3146  real*8 roxb(nbndp2),tetxb(nbndp2)
3147 
3148  real*8 rrk(nbndp4),cck(nbndp4),wrk(nbndp6)
3149 
3150  real*8 cwk(4),tetpol(1),ro0(1)
3151 
3152  xzer(x1,x2,v1,v2) = (x1*v2-x2*v1)/(v2-v1)
3153 
3154  ! write(6,*) '***loop95:enter'
3155 
3156  ! write(6,*) 'loop95:imax,jmax',imax,jmax
3157  !write(6,*) 'loop95:rmax,zmax',rm,zm
3158 
3159  pi=3.14159265358d0
3160 
3161  imax1=imax-1
3162  jmax1=jmax-1
3163 
3164  do 10 i=imax1,imax
3165  if(xm.le.x(i+1) .AND. xm.gt.x(i)) icell=i
3166  10 continue
3167 
3168  do 20 j=jmax1,jmax
3169  if(ym.le.y(j+1) .AND. ym.gt.y(j)) jcell=j
3170  20 continue
3171 
3172  ! write(6,*) 'loop95:icell,jcell',icell,jcell
3173 
3174  i=icell
3175  j=jcell+1
3176 
3177 cccc u0=0.5
3178  do 15 ii=1,ni
3179  do 15 jj=1,nj
3180  ut(ii,jj)=un(ii,jj)-u0
3181  15 continue
3182 
3183  ig=1
3184 
3185  888 continue
3186 
3187  do 885 i=imax,ni1
3188 
3189  ! write(6,*) 'un(i,j)',un(i,j),un(i+1,j)
3190 
3191  if(ut(i,j)*ut(i+1,j).le.0.) then
3192 
3193  ic=i
3194  jc=j
3195 
3196  rxb(ig)=xzer(x(i+1),x(i),ut(i+1,j),ut(i,j))
3197  zxb(ig)=y(j)
3198 
3199  go to 886
3200  endif
3201 
3202  885 continue
3203  write(6,*) 'loop95: first point was not find'
3204  stop
3205 
3206  886 continue
3207 
3208  ic1=ic
3209  jc1=jc
3210 
3211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3212 
3213  lin=1
3214 
3215  100 ig=ig+1
3216 
3217  i=ic
3218  j=jc
3219 
3220  if(ut(i+1,j)*ut(i,j).le.0. .AND. lin.ne.1) then
3221 
3222  ic=i
3223  jc=j-1
3224 
3225  lin=3
3226 
3227  rxb(ig)=xzer(x(i+1),x(i),ut(i+1,j),ut(i,j))
3228  zxb(ig)=y(j)
3229 
3230  elseif(ut(i+1,j+1)*ut(i+1,j).le.0..AND. lin.ne.2) then
3231 
3232  ic=i+1
3233  jc=j
3234 
3235  lin=4
3236 
3237  rxb(ig)=x(i+1)
3238  zxb(ig)=xzer(y(j+1),y(j),ut(i+1,j+1),ut(i+1,j))
3239 
3240  elseif(ut(i+1,j+1)*ut(i,j+1).le.0. .AND. lin.ne.3) then
3241 
3242  ic=i
3243  jc=j+1
3244 
3245  lin=1
3246 
3247  rxb(ig)=xzer(x(i+1),x(i),ut(i+1,j+1),ut(i,j+1))
3248  zxb(ig)=y(j+1)
3249 
3250  elseif(ut(i,j+1)*ut(i,j).le.0..AND. lin.ne.4) then
3251 
3252  ic=i-1
3253  jc=j
3254 
3255  lin=2
3256 
3257  rxb(ig)=x(i)
3258  zxb(ig)=xzer(y(j+1),y(j),ut(i,j+1),ut(i,j))
3259 
3260  endif
3261 
3262 c write(6,*) 'loop:ic,jc',ic,jc,ig
3263 
3264  if(jc.eq.jc1 .AND. ic.eq.ic1) then
3265 
3266  ig=ig+1
3267 
3268  rxb(ig)=rxb(2 )
3269  zxb(ig)=zxb(2 )
3270 
3271  go to 1212
3272 
3273  endif
3274 
3275  go to 100
3276 
3277  1212 nxb=ig
3278 
3279  do 200 ig=1,nxb
3280 
3281  drx=rxb(ig)-xm
3282  dzx=zxb(ig)-ym
3283 
3284  tetp=datan(dzx/drx)
3285 
3286  if(ig.ne.1) then
3287  if(drx.lt.0.d0) tetp=tetp+pi
3288  deltet=tetp-tetxb(ig-1)
3289  if(deltet.lt.0.d0) tetp=tetp+2.d0*pi
3290  endif
3291 
3292  tetxb(ig)=tetp
3293 
3294  roxb(ig)=dsqrt(drx**2+dzx**2)
3295 
3296 c write(6,*) 'ig,ro,tet',ig,roxb(ig),tetxb(ig)
3297 
3298  200 continue
3299 
3300 
3301  CALL e01baf(nxb,tetxb,roxb,rrk,cck,nxb+4,wrk,6*nxb+16,ifail)
3302 
3303 c write(6,*) 'ifail=',ifail
3304 
3305  DO 210 j=1,ntet
3306  tetp=tetpol(j)
3307  if(tetp.lt.tetxb(1)) then
3308  tetp=tetp+2.d0*pi
3309  elseif(tetp.gt.tetxb(nxb)) then
3310  tetp=tetp-2.d0*pi
3311  endif
3312  CALL e02bcf(nxb+4,rrk,cck,tetp ,0,cwk,ifail)
3313  ro0(j)=cwk(1)
3314 c write(6,*) 'ifail=',ifail,j
3315  210 CONTINUE
3316 
3317  return
3318  end
3319 
3320 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3321 
3322 
grid_spdr
subroutine grid_spdr
Definition: _mesh.f:2187
e02bcf
subroutine e02bcf(NCAP7, K, C, X, LEFT, S, IFAIL)
Definition: NAG.f:2761
blin_tr
real *8 function blin_tr(tet0, ro0,
Definition: _fluxt.f:443
fdefln
subroutine fdefln(psidf, rgv, zgv)
Definition: _mesh.f:2990
renet
subroutine renet
Definition: _mesh.f:6
fun
subroutine fun(X, F)
Definition: Ev2.f:10
f_wrd
subroutine f_wrd
Definition: _wrd.f:107
f_regrid0
subroutine f_regrid0(erro)
Definition: _mesh.f:1684
reform
subroutine reform
Definition: _mesh.f:2919
f_loop95
subroutine f_loop95(tetpol, ntet, ro0, u0)
Definition: _mesh.f:3136
f_prgrid
subroutine f_prgrid(imax, jmax, erro)
Definition: _mesh.f:662
bint
subroutine bint(X, Y, R0, Z0, r1, z1, F, I)
Definition: scu.f:557
f_remesh
subroutine f_remesh(erro)
Definition: _mesh.f:617
e01baf
subroutine e01baf(M, X, Y, K, C, LCK, WRK, LWRK, IFAIL)
Definition: NAG.f:2511
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
grid
subroutine grid
Definition: EQ1_m.f:926
ada
subroutine ada(erro)
Definition: _mesh.f:2616
f_grid
subroutine f_grid(igdf, nstep)
Definition: _mesh.f:2316
axdef
subroutine axdef(rma, zma, psima, dpm)
Definition: com_sub.f:86
f_regrid
subroutine f_regrid(erro)
Definition: _mesh.f:138
bonspl
subroutine bonspl
Definition: _mesh.f:3111
cpr
REAL *8 function cpr(H, NP, S)
Definition: scu.f:760
regrid0_
subroutine regrid0_(erro)
Definition: _mesh.f:1195
f_xpoint
subroutine f_xpoint(rx, zx, ix, jx, psx, tet0, kodex)
Definition: _extp.f:672