ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
Eq_m.f
Go to the documentation of this file.
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2  subroutine auto(rk,zk,tk,nk,nstep,ngrid,
3  * rk1,zk1, rk2,zk2,
4  * rk3,zk3, rk4,zk4,
5  * ntipe, necon, wecon )
6 
7  include 'double.inc'
8  include 'param.inc'
9  include 'comblc.inc'
10  include 'urs.inc'
11 
12  integer ngrid,nk,nstep
13  real*8 tk(*),rk(*),zk(*), wecon(*)
14 
15  real*8 rk1(*),zk1(*),rk2(*),zk2(*)
16  real*8 rk3(*),zk3(*),rk4(*),zk4(*)
17 
18  integer ntipe(*), necon(*)
19 
20 
21 c...input initial data
22 
23  write(fname,'(a,a)') path(1:kname),'data.dat'
24  open(1,file=fname,form='formatted')
25  !open(1,file='data.dat')
26  read(1,*) icont
27  read(1,*) ni
28  read(1,*) nj
29  read(1,*) rmin
30  read(1,*) rmax
31  read(1,*) zmin
32  read(1,*) zmax
33  read(1,*) alp
34 
35  read(1,*) nctrl
36 
37  read(1,*) qcen
38  read(1,*) rx10
39  read(1,*) zx10
40  read(1,*) rx20
41  read(1,*) zx20
42  read(1,*) rm0
43  read(1,*) zm0
44  close(1)
45 
46  write(fname,'(a,a)') path(1:kname),'data_d.wr'
47  open(1,file=fname,form='formatted')
48  !open(1,file='data_d.wr')
49  write(1,*) icont,ni,nj,nctrl
50  write(1,*) rmin,rmax,zmin,zmax,alp,qcen,
51  * rx10,zx10,rx20,zx20,rm0,zm0
52  close(1)
53 
54  n_ctrl=nctrl
55  psi_bon=psi_bnd
56 
57  rx1=rx10
58  zx1=zx10
59  rx2=rx20
60  zx2=zx20
61 
62  write(fname,'(a,a)') path(1:kname),'limpnt.dat'
63  open(1,file=fname,form='formatted')
64  !open(1,file='limpnt.dat')
65  read(1,*) nblm
66  do i=1,nblm
67  read(1,*) rblm(i),zblm(i)
68  enddo
69  close(1)
70 
71  write(fname,'(a,a)') path(1:kname),'limpnt_d.wr'
72  open(1,file=fname,form='formatted')
73  !open(1,file='limpnt_d.wr')
74  write(1,*) nblm
75  do i=1,nblm
76  write(1,*) rblm(i),zblm(i)
77  enddo
78  close(1)
79 
80 c...index
81 
82  call glbind
83 
84  call grid
85 
86  if(icont.eq.0) then
87 
88  call exfmat(rk,zk,tk,nk,
89  * rk1,zk1,rk2,zk2,
90  * rk3,zk3,rk4,zk4,ntipe, necon, wecon)
91  !call bndint(rk,zk,nk)
92  !call wrdbin(nk)
93 
94  call cfr_mat(rk,zk,tk,nk,
95  * necon, wecon )
96 
97 
98  endif
99 
100  return
101  end
102 
103 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
104 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
105  subroutine noauto
106 
107  include 'double.inc'
108  include 'param.inc'
109  include 'comblc.inc'
110  include 'urs.inc'
111 
112 
113 c...input initial data
114 
115 
116  write(fname,'(a,a)') path(1:kname),'data_d.wr'
117  open(1,file=fname,form='formatted')
118  !open(1,file='data_d.wr')
119  read(1,*) icont,ni,nj,nctrl
120  read(1,*) rmin,rmax,zmin,zmax,alp,qcen,
121  * rx10,zx10,rx20,zx20,rm0,zm0
122  close(1)
123 
124  !icont=1
125 
126  n_ctrl=nctrl
127  psi_bon=psi_bnd
128 
129  rx1=rx10
130  zx1=zx10
131  rx2=rx20
132  zx2=zx20
133 
134  write(fname,'(a,a)') path(1:kname),'limpnt_d.wr'
135  open(1,file=fname,form='formatted')
136  !open(1,file='limpnt_d.wr')
137  read(1,*) nblm
138  do i=1,nblm
139  read(1,*) rblm(i),zblm(i)
140  enddo
141  close(1)
142 
143  return
144  end
145 
146 
147 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
148  subroutine eq_0(pcequi,psitok,ncequi,nstep,ngrid,
149  * alf0, alf1, alf2, bet0, bet1, bet2,
150  * betpol, betplx, zli3,
151  * ngav1,
152  * ftok, tokout, psicen, pscout,
153  * nursb,psi_bnd,alp_b,rax,zax,n_ctrl,b_0,r_0 )
154 
155  include 'double.inc'
156  include 'param.inc'
157  include 'comblc.inc'
158  include 'urs.inc'
159 
160  common /comomg/ omega,sigma
161  common/comsav/ alf0n
162  common /comhel/ helinp, helout
163 
164  integer isol,ngrid,nstep,ngav1
165  real*8 pcequi(*),psitok(*)
166 
167  real*8 alf0,alf1,alf2,bet0,bet1,bet2
168  real*8 betpol,zli3,ftok,psicen,tokout,pscout
169 
170 
171  alf0p=alf0
172  alf1p=alf1
173  alf2p=alf2
174  bet0f=bet0
175  bet1f=bet1
176  bet2f=bet2
177  b0ax=b_0
178  r0ax=r_0
179 c-----------------------------------------------------------------
180  nnstpp=nstep
181 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
182 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
183 !!!!!!!!!!! !!!!
184 !!!!!!!!! !!!
185 !!!!!!!! !!!
186 
187 c if( nstep.eq.0 .AnD. ngrid.eq.0 ) then
188 c
189 c call tab_build
190 c
191 c endif
192 
193 !!!!!!!! !!!!!!
194 !!!!!!!!!!! !!!!!!!!
195 !!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!
196 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
197 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
198 
199 ccc+ if(nstep.gt.0) go to 6789
200 
201 c...input initial data
202 
203 
204  n_ctrl=nctrl
205  psi_bon=psi_bnd
206 
207  if(nstep.gt.0) icont=1
208 
209 c -----------------------
210  tok = ftok
211  ucen = psicen
212 c -----------------------
213 
214  write(fname,'(a,a)') path(1:kname),'itpr.dat'
215  open(1,file=fname,form='formatted')
216  !open(1,file='itpr.dat')
217  read(1,*) eps0
218  read(1,*) eps
219  read(1,*) epsin
220  read(1,*) epscrz
221  read(1,*) iterbf
222  read(1,*) nitl
223  read(1,*) nitin
224  read(1,*) nrun
225  read(1,*) nwr
226  read(1,*) omg
227  read(1,*) sigm
228  read(1,*) ceps
229  close(1)
230 
231 c...index
232 
233  call glbind
234 
235  ix1=1
236  jx1=1
237 
238  ix2=1
239  jx2=1
240 
241  call grid
242  call geom
243 
244  if(icont.eq.0) then
245  call bndmat
246  call wrdbnd
247  icont=1
248  else
249  call rddbnd
250  endif
251 
252  7890 continue
253 
254  call extfil(pcequi,ncequi)
255 
256  coin=1.d0
257  if(nstep.eq.0) then
258  !!!!!!!ien=0
259  ien=1
260  call taburs(ien,coin,nursb)
261  endif
262 
263  isol = 0
264  nflag = 0
265 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
266 
267  if(ngrid.eq.0) then
268  ni=(ni+1)/2
269  nj=(nj+1)/2
270 
271  alp=alp-0.005d0
272 
273  call glbind
274  call botlev(nstep)
275  call geom
276  endif
277 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
278 
279  777 call matrix
280 
281  !if(nstep.eq.0 .AND. nflag.eq.0) then
282 
283  !call zero(isol)
284  call zero_ax(isol)
285 
286  clr=0.d0
287  clz=0.d0
288 
289  !endif
290 
291  if(nflag.eq.0) call psiful
292 
293  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
294  tokout = tok !!
295  pscout = um !!
296  r_ax = rm !!
297  z_ax = zm !!
298  alp_b = alp !!
299  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
300  call wrd
301 
302  iter=0
303  itin=0
304 
305  return
306  end
307 
308 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
309 
310  subroutine eq( pcequi, psicon, ncequi, nstep, ngrid,
311  * alf0, alf1, alf2, bet0, bet1, bet2,
312  * betpol, betplx, zli3,
313  * ngav1,
314  * ftok, tokout, psicen, pscout,
315  * nursb,psi_bnd,alp_b,rax,zax )
316 
317  include 'double.inc'
318  include 'param.inc'
319  include 'comblc.inc'
320  include 'urs.inc'
321 
322  common /comomg/ omega,sigma
323  common/comsav/ alf0n
324  common /comhel/ helinp, helout
325 
326  common
327  * /c_kpr/kpr
328 
329  integer isol,ngrid,ncequi,nstep,ngav1
330 
331  real*8 alf0,alf1,alf2,bet0,bet1,bet2
332  real*8 betpol,zli3,ftok,psicen,tokout,pscout
333  real*8 psicon(*),pcequi(*)
334 
335  real*4 a_print(200)
336  character *30 apr
337 
338 
339 
340  if(nstep.eq.0) then
341  rm=rax
342  zm=zax
343  endif
344 
345 
346  write(fname,'(a,a)') path(1:kname),'itpr.dat'
347  open(1,file=fname,form='formatted')
348  !open(1,file='itpr.dat')
349  read(1,*) eps0
350  read(1,*) eps
351  read(1,*) epsin
352  read(1,*) epscrz
353  read(1,*) iterbf
354  read(1,*) nitl
355  read(1,*) nitin
356  read(1,*) nrun
357  read(1,*) nwr
358  read(1,*) omg
359  read(1,*) sigm
360  read(1,*) ceps
361  close(1)
362 
363  eps0l= eps0/3.d0
364  epsl = eps/3.d0
365  cepsl = ceps/3.d0
366  epsinl= epsin/3.d0
367 c------------------------------------------------
368 
369  if(nstep.gt.0) then
370  iterbf=1
371  go to 7890
372  endif
373 c...index
374 
375  call glbind
376 
377  ix1=1
378  jx1=1
379 
380  ix2=1
381  jx2=1
382 
383  call grid
384  call geom
385 
386  call rddbnd
387 
388  7890 continue
389 
390  call extfil(pcequi,ncequi)
391 
392  coin=1.d0
393  if(nstep.eq.0) then
394  !!!ien=0
395  ien=1
396  call taburs(ien,coin,nursb)
397  endif
398 
399  isol = 0
400  nflag = 0
401 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
402 
403  if(ngrid.eq.0) then
404  ni=(ni+1)/2
405  nj=(nj+1)/2
406 
407  alp=alp-0.005d0
408 
409  call glbind
410  call botlev(nstep)
411  call geom
412  endif
413 
414  777 call matrix
415 
416  if(nstep.eq.0 .AND. nflag.eq.0) then
417 
418  !call zero(isol)
419  call zero_ax(isol)
420 
421  clr=0.d0
422  clz=0.d0
423 
424  endif
425 
426  if(nflag.eq.0) call psiful
427 
428  erru=1.d0
429  omega=1.d0
430  sigma=1.d0
431  ich=0
432  iflag=0
433 
434  call wrd
435 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
436 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
437 
438  itl=-1
439 
440  iter=0
441  itin=0
442 
443  1000 iter=iter+1
444  itin=itin+1
445 
446  if(iter.eq.iterbf .or. nstep.gt.0) then
447  rl=rm
448  zl=zm
449  call shab(il,jl,icelm,jcelm)
450 cccccccc call wrd
451  endif
452  !call wrd
453  !pause 'pause: wrd'
454 
455 ccc if(iter.eq.2) call wrd
456 
457  if(itin.le.3) then
458  rx1=rx10
459  rx2=rx20
460  endif
461 
462  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
463  !
464  if(ich.eq.0) then !!!!!!!!!!!!!!!!!!!!!!!!!
465  !write(6,*) '** 0 **itin=',itin !
466  !
467  if(erru.lt.eps .OR. itin.ge.nitin) then !
468 
469  itin=0 !
470  !
471  omega=1.d0 !
472  sigma=1.d0 !
473 
474  !write(*,*) 'itl=',itl !
475  !pause 'pause'
476  !
477  itl=itl+1 !
478 
479  !eps=dmax1(eps/1.5,epsl) !
480  !ceps=dmax1(ceps/1.5,ceps) !
481 
482  ich=ich+1 !
483  iflag = 1 !
484  !
485  rl0=rl !
486  zl0=zl !
487  !
488  clr0=clr !
489  clz0=clz !
490  !
491  ddzl=dz(jmax)*ceps !
492  zl=zl+ddzl !
493  call shab(il,jl,icelm,jcelm) !
494  endif !
495  !
496  elseif(ich.eq.1) then !!!!!!!!!!!!!!!!!!!!!
497  !write(6,*) '** 1 **itin=',itin !
498  !
499  if(erru.lt.epsin .OR. itin.ge.nitin) then !
500 
501  itin=0 !
502  !
503  omega=1.d0 !
504  sigma=1.d0 !
505  !
506  ich=ich+1 !
507  !
508  clr1=clr !
509  clz1=clz !
510  !
511  !! clr=0.d0 !
512  !! clz=0.d0 !
513  !
514  ddrl=dr(imax)*ceps !
515  rl=rl+ddrl !
516  zl=zl0 !
517  call shab(il,jl,icelm,jcelm) !
518  endif !
519  !
520  elseif(ich.eq.2) then !!!!!!!!!!!!!!!!!!!!!
521  !write(6,*) '** 2 **itin=',itin !
522  !
523  if(erru.lt.epsin .OR. itin.ge.nitin) then !
524  !!! call wrd
525  !!! write(6,*) ' 2 wrd'
526 
527  itin=0 !
528  !
529  omega=1.d0 !
530  sigma=1.d0 !
531  !
532  dcrdr=(clr-clr0)/ddrl !
533  dczdr=(clz-clz0)/ddrl !
534  !
535  dcrdz=(clr1-clr0)/ddzl !
536  dczdz=(clz1-clz0)/ddzl !
537  !
538  det=dcrdr*dczdz-dczdr*dcrdz !
539  !
540  delrl= (clz0*dcrdz-clr0*dczdz)/det !
541  delzl= (clr0*dczdr-clz0*dcrdr)/det !
542  !
543  dll=dsqrt(delrl**2 + delzl**2) !
544  !
545  dllim=0.25d0*dr(imax) !
546  !
547  if(dll .gt. dllim) then !
548  !
549  nstp=dll/dllim
550  !
551  ddrr=delrl/nstp !
552  ddzz=delzl/nstp
553  !
554  if(nstp.gt.10) nstp=10 !
555  !
556  do 234 istep=1,nstp !
557  !
558  !write(6,*) 'slow shift',istep,nstp !
559  !
560  rl=rl0+ ddrr*istep !
561  zl=zl0+ ddzz*istep !
562  call shab(il,jl,icelm,jcelm) !
563  !
564  778 continue !
565  !
566  rx1=rx10
567  rx2=rx20
568 
569  call right0(il,jl,icelm,jcelm,ngav1)
570  call solve(isol,g) !
571  call bound !
572  !!call right1 !! spar !
573  call solve(isol,ui) !
574  call psiful !
575  !
576  234 continue !
577  !
578  !
579  else !
580  !
581  rl=rl0+ delrl !
582  zl=zl0+ delzl !
583  call shab(il,jl,icelm,jcelm) !
584  endif !
585  !
586  !epsin=dmax1(epsin/2.,epsinl) !
587 
588  ich=0 !
589  !
590  rx1=rx10
591  rx2=rx20
592 
593  endif !
594  !
595  endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
596 
597  if(itl.ge.nitl) then
598  write(*,*) 'limit number of ext. iterations is exeeded'
599  write(*,*) 'itl = ',itl,' Nitl = ',nitl
600  go to 1111
601  endif
602 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
603  if( ngav1.eq.4 .OR. ngav1.eq.5 ) then
604  if(itin.ge.5) then
605  call gridpl
606  call qst(qcen,cnor,b0ax,r0ax)
607  endif
608  endif
609 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
610 
611  call right0(il,jl,icelm,jcelm,ngav1)
612  !write(6,*) 'erru = ',erru
613 
614  call solve(isol,g)
615  call bound
616  !!call right1 !!spar
617  call solve(isol,ui)
618  call psiful
619 
620  if(ngav1.eq.1 .OR. ngav1.eq.3 .OR. ngav1.eq.5) then
621  if(itin.ge.3) then
622 
623  call btpol( betpol )
624  zcoin = betplx / betpol
625 c write(6,*) 'betplx/betpol',zcoin
626  zcoin = (1.d0/zcoin-1.d0)*tok/(cnor*f_cur)+1.d0
627  coin = 1.d0/zcoin
628  coin = (coin-1.d0)*0.33d0+1.d0
629  call taburs( 1,coin,nursb)
630 
631  endif
632  endif
633 
634  ! call wrd
635 
636 
637  crz=dabs(clr*rm*rm/(um-up))+dabs(clz*(zm-zx0)/(um-up))
638 
639  if(erru.le.epsin .AND. crz.lt.epscrz .AND. itl.ge.0
640  * .AND. ich.eq.0) go to 1111
641 
642  omega = omg
643  sigma = sigm
644 
645  if(erru.le.eps.and.kpr.eq.1)then
646  write(*,'("iter erru rm zm ", i4,6(1pe12.5))'),
647  * iter,erru,rm,zm
648 
649  write(*,'("ich clr clz ", i4,6(1pe12.5))'),
650  * ich,clr,clz
651 
652  end if
653 
654 
655  if(iter.le.nrun) go to 1000
656 
657  1111 continue
658 
659 
660  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
661  tokout = tok !!
662  pscout = um !!
663  r_ax = rm !!
664  z_ax = zm !!
665  alp_b = alp !!
666  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
667 
668  if(ngrid.eq.1) then
669  call loop
670  call wrd
671 c call wrdcur
672 c call adapol(alf0,alf1,alf2,bet0,bet1,bet2)
673  call psi_fil( psicon,ncequi )
674 
675 coment call gridpl
676 coment call qst(qcen,cnor,b0ax,r0ax)
677  endif
678  !call output(ngrid,betpol,zli3)
679  call btpol(betpol)
680 c
681 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
682  if(ngrid.eq.0) then !
683  !
684  alp=alp+0.00499999d0 !
685  !
686  ni=2*ni-1 !
687  nj=2*nj-1 !
688  !
689  imax=2*imax-1 !
690  jmax=2*jmax-1 !
691  !
692  ix1=2*ix1-1 !
693  jx1=2*jx1-1 !
694  !
695  ix2=2*ix2-1 !
696  jx2=2*jx2-1 !
697  !
698  call glbind !
699  !
700  call grid !
701  call geom !
702  !
703  call uplev !
704  !
705  call extfil(pcequi,ncequi)
706  call rddbnd !
707  !
708  call psiful !
709  !
710  call inter2(imax,jmax)!
711  !
712  call inter2(ix1,jx1) !
713  !
714  call inter2(ix2,jx2) !
715  !
716 c... if(nglev.eq.1) return !
717  !
718  iterbf=1 !
719  !
720  !!!!!!!!!!! !
721  nitl = 15 ! !
722  !!!!!!!!!!! !
723  !
724  ngrid=1 !
725  nflag=1 !
726  isol=0 !
727  !
728  eps = 1.0d-7 !
729  eps0 = 2.0d-7 !
730  epsin = 1.0d-7 !
731  epscrz = 5.0d-7 !
732 c !
733  ceps=0.02d0 !
734  omega=1.d0 !
735  sigma=1.d0 !
736  !
737  go to 777 !
738  !
739  endif !
740 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
741 
742  call tabnor(cnor)
743  !!!!!!call tabper !!!<<<---------
744  alf0n=alf0
745  write(*,*) '..................................................'
746  write(*,*) 'Convergence of free bound. equilibrium iterations:'
747  write(*,*) ' iter = ',iter,' nrun = ',nrun
748  write(*,*) ' itl = ',itl ,' Nitl = ',nitl
749  write(*,*) ' erru = ',erru
750  write(*,*) ' crz = ',crz
751 
752  return
753  end
754 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
755 c
756  subroutine eq_ax( pcequi, psicon, ncequi, nstep,ngrid,
757  * alf0, alf1, alf2, bet0, bet1, bet2,
758  * betpol, betplx, zli3,
759  * ngav1,
760  * ftok,tokout,psicen,pscout,
761  * ereve0, erps,
762  * psi_bnd,alp_b,rax,zax,isymm )
763 c
764  include 'double.inc'
765  include 'param.inc'
766  include 'comblc.inc'
767 c
768  common /comomg/ omega,sigma
769  common /comhel/ helinp, helout
770  common/comsav/ alf0n
771 c
772  integer isol,ngrid,nk,nstep,ngav1
773  real*8 psicon(*),pcequi(*)
774  real*8 alf0,alf1,alf2,bet0,bet1,bet2
775  real*8 betpol,zli3,ftok,psicen,tokout,pscout
776  real timeb, timee,t_start, t_finish
777 c------------------------------------------------------------------
778  !return
779  if(nstep.gt.0) icont=1
780  nnstpp=nstep
781 c -----------------------
782  tok = ftok
783  ucen = psicen
784  psi_bon = psi_bnd
785 c -----------------------
786  eps0 = ereve0
787  eps = 1.d-7
788  epsin = 1.d-7
789  epscrz = 1.d-7
790 
791  nitl = 9
792  nitin = 61
793 
794  nrun = 1000
795  nwr = 1000
796 
797  omg = 0.99d0
798  sigm = 0.95d0
799 
800  iterbf = 99999999
801 
802  !!coin=alf0/(alf0n+1.d-12)
803  !!call taburs( 1,coin,nursb)
804 
805 c------------------------------------------------
806  !call cpu_time(t_start)
807 
808  call ext_fil(pcequi,ncequi)
809 
810  !call cpu_time(t_finish)
811  !write(*,*) 'ext_fil:opertime=',(t_finish-t_start)
812 
813  isol = 1
814 
815 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
816 
817  clr=0.d0
818  clz=0.d0
819 
820  call psiful
821 
822  erru=1.d0
823  omega=1.d0
824  sigma=1.d0
825  ich=0
826  iter=iter+1
827  itin=itin+1
828 
829  !call wrd
830  !pause 'wrd'
831 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
832 
833 c if(itin.le.3) then
834 c rx1=rx10
835 c rx2=rx20
836 c endif
837 
838  rl=rm
839  zl=zm
840 
841  call shab(il,jl,icelm,jcelm)
842 c------------------------------------------------------------------
843  if( ngav1.eq.4 .OR. ngav1.eq.5 ) then
844  call gridpl
845  call qst(qcen,cnor,b0ax,r0ax)
846  endif
847 c------------------------------------------------------------------
848 
849  call right0(il,jl,icelm,jcelm,ngav1)
850  !call cpu_time(t_finish)
851  !write(*,*) 'right0:opertime=',(t_finish-t_start)
852 
853 c write(*,*) 'erru',erru
854 
855  call solve(isol,g)
856  !call cpu_time(t_finish)
857  !write(*,*) 'solv0:opertime=',(t_finish-t_start)
858  call bound
859  !!call right1 !! spar
860  call solve(isol,ui)
861  !call cpu_time(t_finish)
862  !write(*,*) 'solv1:opertime=',(t_finish-t_start)
863 
864 !!!!>>>>>>symmetr.case
865  if(isymm.eq.1) then
866  call updown
867  endif
868 !!!!>>>>>>symmetr.case
869 
870  call psiful
871 
872  if( ngav1.eq.1 .OR. ngav1.eq.3 .OR. ngav1.eq.5 ) then
873 
874  call btpol( betpol )
875  zcoin = betplx / betpol
876 c write(*,*) 'betplx/betpol',zcoin
877  zcoin = (1.d0/zcoin-1.d0)*tok/(cnor*f_cur)+1.d0
878  coin = 1.d0/zcoin
879 
880  coin = (coin-1.d0)*0.33d0+1.d0
881 
882  call taburs( 1,coin,nursb)
883 
884  endif
885 
886  call psi_fil( psicon,ncequi )
887  !call cpu_time(t_finish)
888  !write(*,*) '+flux:opertime=',(t_finish-t_start)
889 
890  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
891  tokout = tok !!
892  pscout = um !!
893  rax = rm !!
894  zax = zm !!
895  alp_b = alp !!
896  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
897 
898  if(erru.le.eps0) then
899  call loop
900  !call wrd
901 c call wrdcur
902  !call output(ngrid,betpol,zli3)
903  ! call adapol(alf0,alf1,alf2,bet0,bet1,bet2)
904 
905 coment call gridpl
906 coment call qst(qcen,cnor,b0ax,r0ax)
907  !call btpol( betpol )
908  !call bttor(bettor)
909 
910  endif
911 
912  erps=erru
913 
914  call tabnor(cnor)
915  alf0n=alf0
916  !call wrd
917  return
918  end
919 
920 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
921 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
922 
923  subroutine shab(il,jl,icelm,jcelm)
924 
925  include 'double.inc'
926  include 'param.inc'
927  include 'comblc.inc'
928 
929  do 10 i=1,ni1
930 
931  ic=i
932  if( rl.ge.r(i) .AND. rl.lt.r(i+1) ) go to 11
933 
934  10 continue
935  11 continue
936 
937  do 20 j=1,nj1
938 
939  jc=j
940  if( zl.ge.z(j) .AND. zl.lt.z(j+1) ) go to 21
941 
942  20 continue
943  21 continue
944 
945  icelm=ic
946  jcelm=jc
947 
948  sdmn=rmax
949 
950  do k=0,1
951  ilp=ic+k
952  ri=r(ilp)
953  do l=0,1
954  jlp=jc+l
955  zj=z(jlp)
956 
957  ddl=dsqrt( (rl-ri)**2 + (zl-zj)**2 )
958  if(ddl.lt.sdmn) then
959  sdmn=ddl
960  il=ilp
961  jl=jlp
962  endif
963 
964  enddo
965  enddo
966 
967  return
968  end
969 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
970 c
971  subroutine glbind
972 c
973  include 'double.inc'
974  include 'param.inc'
975  include 'comblc.inc'
976 c
977  ni1=ni-1
978  nj1=nj-1
979 c
980  ni2=ni-2
981  nj2=nj-2
982 
983  nbnd=2*(ni1+nj1)
984 
985  return
986  end
987 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
988 c
989  subroutine botlev(nstep)
990 c
991  include 'double.inc'
992  include 'param.inc'
993  include 'comblc.inc'
994 c
995  do 10 i=1,ni
996  10 r(i)=r(2*i-1)
997 
998  do 20 j=1,nj
999  20 z(j)=z(2*j-1)
1000 
1001  do 100 i=1,ni
1002  do 110 j=1,nj
1003 
1004  ue(i,j)=ue(2*i-1,2*j-1)
1005 
1006  if(nstep.eq.0) go to 110
1007 
1008  ui(i,j)=ui(2*i-1,2*j-1)
1009 
1010  110 continue
1011  100 continue
1012 
1013  do 200 ib=1,nbnd
1014  do 200 ibs=1,nbnd
1015 
1016  binadg(ibs,ib)=binadg(2*ibs-1,2*ib-1)+binadg(2*ibs,2*ib-1)
1017 
1018  200 continue
1019 
1020  return
1021  end
1022 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1023 c
1024  subroutine uplev
1025 c
1026  include 'double.inc'
1027  include 'param.inc'
1028  include 'comblc.inc'
1029 c
1030  nnil=(ni+1)/2
1031  nnjl=(nj+1)/2
1032 c
1033  do 100 i=nnil,1,-1
1034  do 100 j=nnjl,1,-1
1035 
1036  ui(2*i-1,2*j-1)=ui(i,j)
1037 
1038  100 continue
1039 
1040  do 200 i=1,ni,2
1041  do 200 j=2,nj1,2
1042 
1043  ui(i,j)=(ui(i,j-1)*dz(j-1)+ui(i,j+1)*dz(j))/(dz(j-1)+dz(j))
1044 
1045  200 continue
1046 
1047  do 300 i=2,ni1,2
1048  do 300 j=1,nj
1049 
1050  ui(i,j)=(ui(i-1,j)*dr(i-1)+ui(i+1,j)*dr(i))/(dr(i-1)+dr(i))
1051 
1052  300 continue
1053 
1054  return
1055  end
1056 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1057 c
1058  subroutine inter2(ip,jp)
1059 c
1060  include 'double.inc'
1061  parameter(nshp=10)
1062  include 'param.inc'
1063  include 'comblc.inc'
1064 c
1065  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
1066 c
1067  if(ip.le.2 .OR. ip.gt.ni2 .OR.
1068  + jp.le.2 .OR. jp.gt.nj2) return
1069 
1070 
1071  nsh=1
1072  xs(nsh)=r(ip)
1073  ys(nsh)=z(jp)
1074  fun(nsh)=u(ip,jp)
1075 
1076  do 200 k=-2,2,2
1077 
1078  i= ip+k
1079 
1080  do 210 l=-2,2,2
1081 
1082  j= jp+l
1083 
1084  if(i.eq.ip .AND. j.eq.jp) go to 210
1085  nsh=nsh+1
1086  xs(nsh)=r(i)
1087  ys(nsh)=z(j)
1088  fun(nsh)=u(i,j)
1089 
1090  210 continue
1091  200 continue
1092 
1093  call deriv5(xs,ys,fun,nsh,5,dp)
1094 
1095  rpi=r(ip)
1096  zpj=z(jp)
1097  upij=u(ip,jp)
1098 
1099  do 100 k=-2,2
1100 
1101  i= ip+k
1102  ri=r(i)
1103  do 110 l=-2,2
1104 
1105  j= jp+l
1106  zj=z(j)
1107 
1108  u(i,j)=upij + dp(1)*(ri-rpi) + dp(2)*(zj-zpj)
1109  + + 0.5*dp(3)*(ri-rpi)*(ri-rpi)
1110  + + dp(4)*(ri-rpi)*(zj-zpj)
1111  + + 0.5*dp(5)*(zj-zpj)*(zj-zpj)
1112 
1113  110 continue
1114  100 continue
1115 
1116  return
1117  end
1118 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1119 
1120  subroutine psiful
1121 
1122  include 'double.inc'
1123  include 'param.inc'
1124  include 'comblc.inc'
1125 
1126  erru=0.d0
1127 
1128  do 10 i=1,ni
1129  do 10 j=1,nj
1130  uold=u(i,j)
1131 
1132  u(i,j)=ui(i,j)+ue(i,j)+clz*z(j)+clr*r(i)*r(i)
1133 
1134  delu=dabs(u(i,j)-uold)
1135  del_um=dabs(um-up)+1.d-9
1136  erru=dmax1(delu/del_um,erru)
1137 
1138  10 continue
1139 
1140 c write(6,*)'erru',erru
1141 c write(6,*)'-------'
1142 
1143  return
1144  end
1145 
1146 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1147 
1148  subroutine tab_buil(alf0, alf1, alf2, bet0, bet1, bet2)
1149 
1150  include 'double.inc'
1151  parameter(ntabp=500)
1152  real*8 psi(ntabp),q(ntabp),dp(ntabp),df(ntabp)
1153 
1154  write(fname,'(a,a)') path(1:kname),'eoh2.dat'
1155  open(1,file=fname,form='formatted')
1156  !open(1,file='eoh2.DAT')
1157  read(1,*) ntab
1158  do i=1,ntab
1159  read(1,*) psi(i),q(i),dp(i)
1160  enddo
1161  dpnor=dp(1)
1162  do i=1,ntab
1163  dp(i)=alf0*dp(i)/dpnor
1164  zpsi=1.d0-psi(i)/psi(ntab)
1165  df(i)=funfp(zpsi)
1166  !dp(i)=funpp(zpsi,alf0,alf1,alf2,bet0,bet1,bet2)
1167  enddo
1168  close(1)
1169 
1170  write(fname,'(a,a)') path(1:kname),'tabppf.dat'
1171  open(1,file=fname,form='formatted')
1172  !open(1,file='tabppf.dat')
1173  write(1,*) ntab
1174  do i=1,ntab
1175  write(1,*) psi(i),dp(i),df(i)
1176  enddo
1177  close(1)
1178 
1179  write(fname,'(a,a)') path(1:kname),'tab_q.dat'
1180  open(1,file=fname,form='formatted')
1181  !open(1,file='tab_q.dat')
1182  write(1,*) ntab, (0)
1183  do i=1,ntab
1184  write(1,*) psi(i),q(i)
1185  enddo
1186  close(1)
1187 
1188  return
1189  end
1190 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1191  subroutine tab_build
1192  include 'double.inc'
1193  parameter(np=500)
1194  real*8 ps(np),q(np),p(np),f(np)
1195 
1196  pi=3.14159265358d0
1197 
1198  write(fname,'(a,a)') path(1:kname),'tabpfq.dat'
1199  open(1,file=fname,form='formatted')
1200  !open(1,file='tabpfq.dat')
1201  read(1,*) n
1202  do i=1,n
1203  read(1,*) ps(i),q(i),p(i),f(i)
1204  p(i)=p(i)*4.d-7*pi
1205  enddo
1206  close(1)
1207 
1208  write(fname,'(a,a)') path(1:kname),'tabppf.dat'
1209  open(1,file=fname,form='formatted')
1210  !open(1,file='tabppf.dat')
1211  write(1,*) n
1212  do i=1,n
1213  write(1,*) ps(i),p(i),f(i)
1214  enddo
1215  close(1)
1216 
1217  write(fname,'(a,a)') path(1:kname),'tab_q.dat'
1218  open(1,file=fname,form='formatted')
1219  !open(1,file='tab_q.dat')
1220  write(1,*) n, (0)
1221  do i=1,n
1222  write(1,*) ps(i),q(i)
1223  enddo
1224  close(1)
1225 
1226  !write(6,*) 'tab_build: tabppf.dat,tab_q.dat was created'
1227  ! pause 'pause'
1228 
1229  return
1230  end
1231 
1232 
1233 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1234  subroutine q_pert(psizv)
1235 
1236  include'double.inc'
1237 
1238  parameter(np=500)
1239  real*8 ps(np),q(np)
1240 
1241  write(fname,'(a,a)') path(1:kname),'tab_q.dat'
1242  open(1,file=fname,form='formatted')
1243  !open(1,file='tab_q.dat')
1244  read(1,*) n
1245  do i=1,n
1246  read(1,*) ps(i),q(i)
1247  enddo
1248  close(1)
1249 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1250 ! the first variant of q-perturbation ( psizv = psizv )
1251 !
1252 ! do i=2,n
1253 ! if( ps(i) .gt. psizv ) then
1254 ! izv = i-1
1255 ! qzv = q(i-1)
1256 ! go to 10
1257 ! endif
1258 ! enddo
1259 !
1260 ! 10 continue
1261 !
1262 ! do i=1,izv
1263 ! q(i)=qzv
1264 ! enddo
1265 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1266 ! the second variant of q-perturbation ( psizv = qzv )
1267 !
1268  qzv = psizv
1269  psizv = 0.
1270  izv = 0
1271  do i=1,n
1272  if( q(i) .lt. qzv ) then
1273  q(i) = qzv
1274  izv = i
1275  psizv = ps(i)
1276  endif
1277  enddo
1278 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1279 !
1280  write(fname,'(a,a)') path(1:kname),'tab_q.dat'
1281  open(1,file=fname,form='formatted')
1282  !open(1,file='tab_q.dat')
1283  write(1,*) n, (0)
1284  do i=1,n
1285  write(1,*) ps(i),q(i)
1286  enddo
1287  close(1)
1288 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1289 !
1290 
1291  return
1292  end
1293 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1294  subroutine updown
1295 
1296  include 'double.inc'
1297  include 'param.inc'
1298  include 'comblc.inc'
1299 
1300  nj12= (nj+1)/2
1301 
1302  do 10 j=1,nj12
1303  do 10 i=1,ni
1304 
1305  udown=ui(i,j)
1306  uup=ui(i,nj-j+1)
1307  usym=0.5d0*(uup+udown)
1308  ui(i,j)=usym
1309  ui(i,nj-j+1)=usym
1310  10 continue
1311 
1312  return
1313  end
1314 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1315 
1316  subroutine get_par(psi_bnd)
1317 
1318  include 'double.inc'
1319  include 'param.inc'
1320  include 'comblc.inc'
1321  psi_bnd=up
1322  return
1323  end
1324 
1325 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1326 
1327 
1328 
1329 
1330 
1331 
q_pert
subroutine q_pert(psizv)
Definition: Eq_m.f:1234
updown
subroutine updown
Definition: Eq_m.f:1294
exfmat
subroutine exfmat(rk, zk, tk, nk,
Definition: EQ1_m.f:107
solve
subroutine solve(isol, wdm)
Definition: spider_solver.f:2
wrdbnd
subroutine wrdbnd
Definition: Eq2_m.f:1576
get_par
subroutine get_par(psi_bnd)
Definition: Eq_m.f:1316
fun
subroutine fun(X, F)
Definition: Ev2.f:10
gridpl
subroutine gridpl
Definition: B_MESH.f:156
bound
subroutine bound
Definition: RIG.f:1
eq
subroutine eq(pcequi, psicon, ncequi, nstep, ngrid,
Definition: Eq_m.f:310
tab_build
subroutine tab_build
Definition: Eq_m.f:1191
right0
subroutine right0(ill, jll, icelm, jcelm, ngav1)
Definition: RIG.f:119
rddbnd
subroutine rddbnd
Definition: Eq2_m.f:1619
noauto
subroutine noauto
Definition: Eq_m.f:105
bndmat
subroutine bndmat
Definition: EQ1_m.f:569
tab_buil
subroutine tab_buil(alf0, alf1, alf2, bet0, bet1, bet2)
Definition: Eq_m.f:1148
cfr_mat
subroutine cfr_mat(rk, zk, tk, nk,
Definition: EQ1_m.f:225
auto
subroutine auto(rk, zk, tk, nk, nstep, ngrid,
Definition: Eq_m.f:2
wrd
subroutine wrd
Definition: Eq2_m.f:1378
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
adapol
subroutine adapol
Definition: Eq2_m.f:969
grid
subroutine grid
Definition: EQ1_m.f:926
eq_ax
subroutine eq_ax(pcequi, psicon, ncequi, nstep, ngrid,
Definition: Eq_m.f:756
glbind
subroutine glbind
Definition: Eq_m.f:971
loop
subroutine loop
Definition: Eq2_m.f:231
psiful
subroutine psiful
Definition: Eq_m.f:1120
shab
subroutine shab(il, jl, icelm, jcelm)
Definition: Eq_m.f:923
inter2
subroutine inter2(ip, jp)
Definition: Eq_m.f:1058
phys_functions::p
real(r8) function p(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:383
extfil
subroutine extfil(pcequi, ncequi)
Definition: EQ1_m.f:47
ext_fil
subroutine ext_fil(pcequi, ncequi)
Definition: EQ1_m.f:469
uplev
subroutine uplev
Definition: Eq_m.f:1024
qst
subroutine qst(qcen, cnor, b0ax, r0ax)
Definition: B_MESH.f:2
btpol
subroutine btpol(betpol)
Definition: Eq2_m.f:1118
eq_0
subroutine eq_0(pcequi, psitok, ncequi, nstep, ngrid,
Definition: Eq_m.f:148
psi_fil
subroutine psi_fil(psicon, ncequi)
Definition: EQ1_m.f:3
botlev
subroutine botlev(nstep)
Definition: Eq_m.f:989
geom
subroutine geom
Definition: EQ1_m.f:953
zero_ax
subroutine zero_ax(isol)
Definition: EQ1_m.f:1099
matrix
subroutine matrix
Definition: Matrix.f:2