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