ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
EQ1_m.f
Go to the documentation of this file.
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 
3  subroutine psi_fil( psicon,ncequi )
4 
5  include 'double.inc'
6  include 'prm.inc'
7  include 'param_.inc'
8  include 'comblc.inc'
9  include 'comevl.inc'
10 
11  common /comext/ zaindk(nip,njp,nkp)
12 
13  real *4 zaindk
14 
15  real*8 psicon(*)
16 
17  !ncpfc=nloc(npfc)
18 
19  !NCEQUI = Nk - NCPFC + NEQUI !!!!
20 
21  do iq=1,ncequi
22 
23  zpscon=0.d0
24 
25  do j=2,nj1
26  do i=2,ni1
27 
28  !if(ipr(i,j).ne.0) then
29  vrftfa=zaindk(i,j,iq)
30  zpscon=zpscon+vrftfa*curf(i,j)*dri(i)*dzj(j)
31  !zpscon=zpscon+vrftfa*curf(i,j)*dri(i)*dzj(j)*amu0
32  !endif
33 
34  enddo
35  enddo
36 
37  psicon(iq)=zpscon
38  enddo
39 
40  return
41  end
42 
43 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
44 
45 
46 
47  subroutine extfil(pcequi,ncequi)
48 
49  include 'double.inc'
50  include 'prm.inc'
51  include 'param_.inc'
52  include 'comblc.inc'
53  include 'comevl.inc'
54 
55 
56  real*8 aindk(nip,njp)
57  real*8 pcequi(*)
58  integer ncequi
59  logical*4 exi
60 
61  !!amu0=0.4d0*pi
62 
63  ncpfc=nloc(npfc)
64  !NCEQUI = Nk - NCPFC + NEQUI !!!!
65 
66  do 500 j=1,nj
67  do 500 i=1,ni
68 
69  ue(i,j)=0.d0
70 
71  500 continue
72 
73  write(fname,'(a,a)') path(1:kname),'exf.wr'
74  open(1,file=fname,form='formatted')
75  !open(1,file='exf.wr',form='formatted')
76 
77  read(1,*) nkread,nequi
78  !if(nk.ne.nkread) then
79  ! write(6,*) '***ERROR:'
80  ! write(6,*)' you change total number of currents'
81  ! write(6,*)' nk_read=',nkread,'nk from trecur=',nk
82  ! write(6,*)' put icont=0 and start again'
83  ! stop
84  !endif
85 
86 
87  do iq=1,ncequi
88 
89  read(1,*) ((aindk(i,j),j=1,nj),i=1,ni)
90 
91  do j=1,nj
92  do i=1,ni
93 
94  ue(i,j)=ue(i,j)+aindk(i,j)*pcequi(iq)*amu0
95 
96  enddo
97  enddo
98 
99  enddo
100 
101 
102  close(1)
103 
104  return
105  end
106 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
107  subroutine exfmat(rk,zk,tk,nk,
108  * rk1,zk1,rk2,zk2,
109  * rk3,zk3,rk4,zk4,ntipe,
110  * necon, wecon )
111 
112  include 'double.inc'
113  include 'prm.inc'
114  include 'param_.inc'
115  include 'comblc.inc'
116  include 'comevl.inc'
117 
118 
119  real*8 aindk(nip,njp)
120  real*8 rk(*),zk(*),tk(*), wecon(*)
121  real*8 rk1(*),zk1(*),rk2(*),zk2(*)
122  real*8 rk3(*),zk3(*),rk4(*),zk4(*)
123  integer nk
124  integer ntipe(*), necon(*)
125  logical*4 exi
126 
127  !!amu0=0.4d0*pi
128 
129  ncpfc=nloc(npfc)
130  nves = nk - ncpfc
131  ncequi = nequi + nves
132 
133  write(fname,'(a,a)') path(1:kname),'exf.wr'
134  open(1,file=fname,form='formatted')
135  !open(1,file='exf.wr',form='formatted')
136 
137  write(1,*) ncequi, nequi
138 
139  do 1010 iq=1,nequi
140 
141  do i=1,ni
142  do j=1,nj
143  aindk(i,j)=0.d0
144 
145  enddo
146  enddo
147 
148  do i=1,ni
149  do j=1,nj
150 
151  do ik=1,ncpfc
152  if( necon(ik) .eq. iq ) then
153  ddx=dsqrt( (r(i)-rk(ik))**2 + (z(j)-zk(ik))**2 )
154  zaindk=greeni(r(i),z(j),rk(ik),zk(ik))/pi
155  aindk(i,j)=aindk(i,j)+zaindk*wecon(ik)
156  end if
157  enddo
158 
159  enddo
160  enddo
161 
162  write(1,*) ((aindk(i,j),j=1,nj),i=1,ni)
163 
164  1010 continue
165 
166  ibeg=ncpfc + 1
167 
168  if(ibeg.gt.nk) go to 11
169 
170 
171 c print *,' ibeg nk ncam==',ibeg,nk,nk-ibeg
172 
173 c read (*,*)
174 
175  ncam=0.
176 
177  do 10 ik=ibeg,nk
178 
179  ntip=ntipe(ik)
180 
181  if(ntip.eq.1) then
182  r1=rk1(ik)
183  z1=zk1(ik)
184  r2=rk2(ik)
185  z2=zk2(ik)
186  dlong=dsqrt( (r2-r1)**2 + (z2-z1)**2 )
187  endif
188 
189  do 100 i=1,ni
190  do 100 j=1,nj
191 
192  ddx=dsqrt( (r(i)-rk(ik))**2 + (z(j)-zk(ik))**2 )
193 
194  if(ntip.eq.1) then
195  call bint(r(i),z(j),r1,z1,r2,z2,fint,1)
196  aindk(i,j)=fint/dlong
197  else
198  aindk(i,j)=greeni(r(i),z(j),rk(ik),zk(ik))/pi
199  endif
200  100 continue
201 
202  write(1,*) ((aindk(i,j),j=1,nj),i=1,ni)
203 
204  ncam=ncam+1
205 
206  !write(6,*) 'wrd ik=',ik
207  10 continue
208  11 continue
209 
210  write(1,*) ncam
211 
212 
213 
214  close(1)
215  !write(6,*) 'wrd nk',nk
216 
217 
218 
219  1000 continue
220 
221  return
222  end
223 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
224 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
225  subroutine cfr_mat(rk,zk,tk,nk,
226  * necon, wecon )
227 
228  include 'double.inc'
229  include 'prm.inc'
230  include 'param_.inc'
231  include 'comblc.inc'
232  include 'comevl.inc'
233 
234 
235  real*8 fr_indk(npfc0,npfc0),fz_indk(npfc0,npfc0)
236  real*8 fpr_indk(nip,njp),fpz_indk(nip,njp)
237  real*8 rk(*),zk(*),tk(*), wecon(*)
238  integer nk
239  integer necon(*)
240  logical*4 exi
241 
242 
243 
244 !!!!!!!!!coil coil interaction
245 
246  do i=1,npfc !2
247  do j=1,npfc !2
248 ! coils number i j
249 
250  sumfrij=0.d0
251  sumfzij=0.d0
252 
253  if(i.ne.j) then
254 
255  ib=nloc(i-1)+1
256  ie=nloc(i)
257 
258  jb=nloc(j-1)+1
259  je=nloc(j)
260 
261  do k=ib,ie
262  do m=jb,je
263  rkk=rk(k) !1.d0
264  zkk=zk(k) !1.d-2
265  rkm=rk(m) !1.d0
266  zkm=zk(m) !-1.d-2
267  !weconk=wecon(k)
268  !weconm=wecon(m)
269  !val_pi=pi
270  !call GRENg(rk(k),zk(k), rk(m),zk(m), fgreen,dGdr,dGdz)
271  !call grGREN(rk(k),zk(k), rk(m),zk(m), dGdr,dGdz)
272  call grgren(rkk,zkk, rkm,zkm, dgdr,dgdz)
273  sumfrij=sumfrij+dgdr*wecon(k)*wecon(m)/rk(m)/pi
274  sumfzij=sumfzij-dgdz*wecon(k)*wecon(m)*2.d0
275  enddo
276  enddo
277 
278  fr_indk(i,j)=sumfrij
279  fz_indk(i,j)=sumfzij
280 
281  endif
282 
283 ! Fr(i,j)= -SumFRij*Ii*Ij !force on coil j
284 ! Fz(i,j)= -SumFZij*Ii*Ij !force on coil j
285 
286  enddo
287  enddo
288 
289 
290  write(fname,'(a,a)') path(1:kname),'forcmat.wr'
291  open(1,file=fname,form='formatted')
292 
293  write(1,*) npfc
294  write(1,*) ((fr_indk(i,j),j=1,npfc),i=1,npfc)
295  write(1,*) ((fz_indk(i,j),j=1,npfc),i=1,npfc)
296 
297 !!!!!!!!!plasma coil interaction
298 
299  do k=1,npfc
300 
301  kb=nloc(k-1)+1
302  ke=nloc(k)
303 
304  do i=1,ni
305  do j=1,nj
306  sumfrpijk=0.d0
307  sumfzpijk=0.d0
308 
309  do ik=kb,ke
310  call grgren(r(i),z(j), rk(ik),zk(ik), dgdr,dgdz)
311  sumfrpijk=sumfrpijk+dgdr*wecon(ik)/rk(ik)/pi
312  sumfzpijk=sumfzpijk-dgdz*wecon(ik)*2.d0
313  enddo
314  fpr_indk(i,j)= sumfrpijk
315  fpz_indk(i,j)= sumfzpijk
316 
317  enddo
318  enddo
319 
320  write(1,*) ((fpr_indk(i,j),j=1,nj),i=1,ni)
321  write(1,*) ((fpz_indk(i,j),j=1,nj),i=1,ni)
322 
323  enddo
324 
325  close(1)
326 
327  return
328  end
329 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
330 
331 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
332  subroutine coil_force
333 
334  include 'double.inc'
335  include 'prm.inc'
336  include 'param_.inc'
337  include 'comblc.inc'
338  include 'comevl.inc'
339 
340  real*8 frcr(npfc0),frcz(npfc0)
341 
342  real*8 fr_indk(npfc0,npfc0),fz_indk(npfc0,npfc0)
343  real*8 fpr_indk(nip,njp),fpz_indk(nip,njp)
344  logical*4 exi
345 
346 ! coil forces calculation
347 ! forces frcz are full Fz forces on coil[MN]
348 ! forces frcr are Fr forces on coil lenth unit[MN/m]
349 
350  write(fname,'(a,a)') path(1:kname),'forcmat.wr'
351  open(1,file=fname,form='formatted')
352 
353  read(1,*) npfc
354  read(1,*) ((fr_indk(i,j),j=1,npfc),i=1,npfc)
355  read(1,*) ((fz_indk(i,j),j=1,npfc),i=1,npfc)
356 
357 
358 
359  do k=1,npfc
360  do iq=1,nequi
361  if(nepfc(k) .eq. iq) then
362  pfcur1(k)=pfceqw(iq)
363  exit
364  endif
365  enddo
366  enddo
367 
368 
369 !!!!!!!!!coil coil interaction
370 
371  do m=1,npfc
372  frsum=0.d0
373  fzsum=0.d0
374  do k=1,npfc
375  if(k.ne.m)then
376  frsum=frsum +fr_indk(k,m)*pfcur1(k)*pfcur1(m)
377  fzsum=fzsum -fz_indk(k,m)*pfcur1(k)*pfcur1(m)
378  endif
379  enddo
380  frcr(m)=frsum*amu0
381  frcz(m)=fzsum*amu0
382  enddo
383 
384 !!!!!!!!!plasma coil interaction
385 
386  fzpl=0.d0
387 
388  do k=1,npfc
389 
390  read(1,*) ((fpr_indk(i,j),j=1,nj),i=1,ni)
391  read(1,*) ((fpz_indk(i,j),j=1,nj),i=1,ni)
392 
393  do i=1,ni
394  do j=1,nj
395  iprij=ipr(i,j)
396  if(iprij.eq.1)then
397  dcurij=curf(i,j)*(dri(i)*dzj(j))
398  frcr(k)=frcr(k)+fpr_indk(i,j)*dcurij*pfcur1(k)
399  frcz(k)=frcz(k)-fpz_indk(i,j)*dcurij*pfcur1(k)
400  fzpl=fzpl-fpr_indk(i,j)*dcurij*pfcur1(k)
401  endif
402  enddo
403  enddo
404 
405  enddo
406 
407  close(1)
408 
409  write(fname,'(a,a)') path(1:kname),'coil_forces.wr'
410  open(1,file=fname,form='formatted')
411  write(1,*) npfc
412  write(1,*) (frcr(k),k=1,npfc)
413  write(1,*) (frcz(k),k=1,npfc)
414  close(1)
415 
416 !force limits for ITER
417 
418  !call iter_force_limits(frcz,frcr)
419 
420 !force limits for ITER
421 
422  open(1,file='coil_forces.pr',form='formatted')
423  write(1,'(a6,e13.5)') 'PF1',frcz(1)
424  write(1,'(a6,e13.5)') 'PF2',frcz(7)
425  write(1,'(a6,e13.5)') 'PF3',frcz(3)
426  write(1,'(a6,e13.5)') 'PF5',frcz(5)
427  write(1,'(a6,e13.5)') 'CS ',frcz(9)
428  close(1)
429 
430 
431 
432  return
433  end
434 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
435 
436  subroutine iter_force_limits(frcz,frcr)
437 
438  include 'double.inc'
439  dimension frcz(*),frcr(*)
440 
441  fzcsup1= frcz(5)
442  fzcsup2= fzcsup1+frcz(3)
443  fzcsup3= fzcsup2+frcz(1)
444  fzcsup4= fzcsup3+frcz(2)
445  fzcsup5= fzcsup4+frcz(4)
446  fzcsup6= fzcsup5+frcz(6)
447  fzmxup= dmax1(fzcsup1,fzcsup2,fzcsup3,fzcsup4,fzcsup5,fzcsup6)
448 
449 
450  fzcsdw1= frcz(6)
451  fzcsdw2= fzcsdw1+frcz(4)
452  fzcsdw3= fzcsdw2+frcz(2)
453  fzcsdw4= fzcsdw3+frcz(1)
454  fzcsdw5= fzcsdw4+frcz(3)
455  fzcsdw6= fzcsdw5+frcz(5)
456  fzmxdw= dmin1(fzcsdw1,fzcsdw2,fzcsdw3,fzcsdw4,fzcsdw5,fzcsdw6)
457 
458  fzpf1= frcz(7)
459  fzpf2= frcz(8)
460  fzpf3= frcz(9)
461  fzpf4= frcz(10)
462  fzpf5= frcz(11)
463  fzpf6= frcz(12)
464 
465  return
466  end
467 
468 
469  subroutine ext_fil(pcequi,ncequi)
470 
471  include 'double.inc'
472  include 'prm.inc'
473  include 'param_.inc'
474  include 'comblc.inc'
475  include 'comevl.inc'
476 
477  common /comext/ zaindk(nip,njp,nkp)
478  real *4 zaindk
479 
480  real*8 pcequi(*)
481  integer ncequi
482 
483  !!amu0=0.4d0*pi
484  !ncpfc=nloc(npfc)
485 
486  1000 continue
487 
488  do j=1,nj
489  do i=1,ni
490  ue(i,j)=0.d0
491  enddo
492  enddo
493 
494  do iq=1,ncequi
495  egcurr=pcequi(iq)
496  do i=1,ni
497  do j=1,nj
498  vrftfa=zaindk(i,j,iq)
499  ue(i,j)=ue(i,j)+vrftfa*egcurr*amu0
500  enddo
501  enddo
502  enddo
503 
504  return
505  end
506 c
507 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
508 c
509  subroutine rdexf(ncequi)
510 
511  include 'double.inc'
512  include 'prm.inc'
513  include 'param_.inc'
514  include 'comblc.inc'
515  include 'comevl.inc'
516 
517  common /comext/ zaindk(nip,njp,nkp)
518 
519  real*4 zaindk
520 
521  real*8 aindk(nip,njp)
522 
523  write(fname,'(a,a)') path(1:kname),'exf.wr'
524  open(1,file=fname,form='formatted')
525  !open(1,file='exf.wr',form='formatted')
526 
527  read(1,*) nk,nequi
528 
529  !if(nk.ne.nkread) then
530  ! write(6,*) '***ERROR:'
531  !write(6,*)' you change total number of currents'
532  !write(6,*)' nk_read=',nkread,'nk from trecur=',nk
533  !write(6,*)' put icont=0 and start again'
534  !stop
535  !endif
536 
537  ncpfc=nloc(npfc)
538 
539  do iq=1,ncequi
540  read(1,*) ((aindk(i,j),j=1,nj),i=1,ni)
541  do j=1,nj
542  do i=1,ni
543  zaindk(i,j,iq)=aindk(i,j)
544  enddo
545  enddo
546  enddo
547 
548 ! ibeg=ncpfc + 1
549 ! if(ibeg.gt.nk) go to 22
550 ! iq=nequi
551 ! do ik=ibeg,nk
552 ! iq=iq+1
553 ! read(1) ((aindk(i,j),i=1,ni),j=1,nj)
554 ! do i=1,ni
555 ! do j=1,nj
556 ! zaindk(i,j,iq)=aindk(i,j)
557 ! enddo
558 ! enddo
559 ! enddo
560 !
561 ! 22 continue
562 
563  close(1)
564 
565  return
566  end
567 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
568 
569  subroutine bndmat
570 
571  include 'double.inc'
572  include 'param.inc'
573  include 'comblc.inc'
574 
575  real*8 xk(nbndp),yk(nbndp)
576 
577  ib=1
578 
579  xk(ib)=r(1)
580  yk(ib)=z(1)
581 
582  do 10 i=2,ni1
583 
584  ib=ib+1
585  xk(ib)=r(i)
586  yk(ib)=z(1)
587 
588  10 continue
589 
590  ib=ib+1
591  xk(ib)=r(ni)
592  yk(ib)=z(1)
593 
594  do 20 j=2,nj1
595 
596  ib=ib+1
597  xk(ib)=r(ni)
598  yk(ib)=z(j)
599 
600  20 continue
601 
602  ib=ib+1
603  xk(ib)=r(ni)
604  yk(ib)=z(nj)
605 
606  do 30 i=ni1,2,-1
607 
608  ib=ib+1
609  xk(ib)=r(i)
610  yk(ib)=z(nj)
611 
612  30 continue
613 
614  ib=ib+1
615  xk(ib)=r(1)
616  yk(ib)=z(nj)
617 
618  do 40 j=nj1,2,-1
619 
620  ib=ib+1
621  xk(ib)=r(1)
622  yk(ib)=z(j)
623 
624  40 continue
625 
626  ib=ib+1
627  xk(ib)=r(1)
628  yk(ib)=z(1)
629 
630 
631  do ib=1,nbnd
632  rr=xk(ib)
633  zz=yk(ib)
634  do ibc=1,nbnd
635 
636  r0=xk(ibc)
637  z0=yk(ibc)
638 
639  r1=xk(ibc+1)
640  z1=yk(ibc+1)
641 
642  call bint(rr,zz,r0,z0,r1,z1,fint,1)
643 
644  binadg(ibc,ib)=fint
645 
646  enddo
647  enddo
648 
649 
650 
651  return
652  end
653 
654 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
655  subroutine bndint(rk,zk,nk)
656 
657  include 'double.inc'
658  include 'param.inc'
659  include 'comblc.inc'
660 
661  real*8 xk(nbndp),yk(nbndp)
662  real*8 rk(*),zk(*)
663  integer nk
664 
665  ib=1
666 
667  xk(ib)=r(1)
668  yk(ib)=z(1)
669 
670  do 10 i=2,ni1
671 
672  ib=ib+1
673  xk(ib)=r(i)
674  yk(ib)=z(1)
675 
676  10 continue
677 
678  ib=ib+1
679  xk(ib)=r(ni)
680  yk(ib)=z(1)
681 
682  do 20 j=2,nj1
683 
684  ib=ib+1
685  xk(ib)=r(ni)
686  yk(ib)=z(j)
687 
688  20 continue
689 
690  ib=ib+1
691  xk(ib)=r(ni)
692  yk(ib)=z(nj)
693 
694  do 30 i=ni1,2,-1
695 
696  ib=ib+1
697  xk(ib)=r(i)
698  yk(ib)=z(nj)
699 
700  30 continue
701 
702  ib=ib+1
703  xk(ib)=r(1)
704  yk(ib)=z(nj)
705 
706  do 40 j=nj1,2,-1
707 
708  ib=ib+1
709  xk(ib)=r(1)
710  yk(ib)=z(j)
711 
712  40 continue
713 
714  ib=ib+1
715  xk(ib)=r(1)
716  yk(ib)=z(1)
717 
718 
719  nkin=0
720  nkout=0
721 
722  do 500 ik=1,nk
723  rr=rk(ik)
724  zz=zk(ik)
725 
726  if(rr.gt.rmax .OR. rr.lt.rmin .OR.
727  * zz.gt.zmax .OR. zz.lt.zmin ) then
728 
729  nkout=nkout+1
730 
731  do 110 ib=1,nbnd
732 
733  r0=xk(ib)
734  z0=yk(ib)
735 
736  r1=xk(ib+1)
737  z1=yk(ib+1)
738 
739  call bint(rr,zz,r0,z0,r1,z1,fint,1)
740 
741  pinadg(nkout,ib)=fint
742 
743  110 continue
744  itok(ik)=0
745  jtok(ik)=0
746 
747  else
748 
749  nkin=nkin+1
750 
751  do 200 i=1,ni1
752 
753  if(rr.ge.r(i) .AND. rr.le.r(i+1)) then
754 
755  itok(ik)=i
756 
757  go to 210
758 
759  endif
760 
761  200 continue
762  210 continue
763 
764  do 300 j=1,nj1
765 
766  if(zz.ge.z(j) .AND. zz.le.z(j+1)) then
767 
768  jtok(ik)=j
769 
770  go to 310
771 
772  endif
773 
774  300 continue
775  310 continue
776 
777  endif
778 
779  500 continue
780 
781  return
782  end
783 
784 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
785 
786  subroutine flux(psitok,rk,zk,nk)
787 
788  include 'double.inc'
789  include 'param.inc'
790  include 'comblc.inc'
791 
792  real*8 psitok(*),rk(*),zk(*)
793 
794  ikout=0
795  ikin=0
796 
797  do 100 ik=1,nk
798 
799  psitok(ik)=0.d0
800 
801  if(itok(ik).eq.0) then
802 
803  ikout=ikout+1
804 
805  do 110 ib=1,nbnd
806 
807  psitok(ik)=psitok(ik)+pinadg(ikout,ib)*dgdn(ib)
808 
809  110 continue
810 
811  else
812 
813  ikin=ikin+1
814 
815  i=itok(ik)
816  j=jtok(ik)
817 
818  psitok(ik)=blin(i,j,rk(ik),zk(ik))
819 
820  endif
821 
822  100 continue
823 
824  if(nkin.ne.ikin .OR. nkout.ne.ikout) then
825  !write(6,*) '*** ERROR: nkin.ne.ikin .OR. nkout.ne.ikout'
826  !write(6,*) ' nkin, ikin = ', nkin, ikin
827  !write(6,*) ' nkout, ikout = ', nkout, ikout
828  stop
829  endif
830 
831  return
832  end
833 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
834 c
835  real*8 function blin(i,j,r0,z0)
836 c
837  include 'double.inc'
838  include 'param.inc'
839  include 'comblc.inc'
840 c
841  z1=z(j)
842  z2=z(j+1)
843 
844  r1=r(i)
845  r2=r(i+1)
846 
847  u1=ui(i,j)
848  u2=ui(i+1,j)
849  u3=ui(i+1,j+1)
850  u4=ui(i,j+1)
851 
852  s1=(r2-r0)*(z2-z0)
853  s3=(r0-r1)*(z0-z1)
854  s2=(r0-r1)*(z2-z0)
855  s4=(r2-r0)*(z0-z1)
856 
857  u0=(s1*u1+s2*u2+s3*u3+s4*u4)/(s1+s2+s3+s4)
858 
859  blin=u0
860 
861  return
862  end
863 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
864 c
865  real*8 function bline(i,j,r0,z0)
866 c
867  include 'double.inc'
868  include 'param.inc'
869  include 'comblc.inc'
870 c
871  z1=z(j)
872  z2=z(j+1)
873 
874  r1=r(i)
875  r2=r(i+1)
876 
877  u1=ue(i,j)
878  u2=ue(i+1,j)
879  u3=ue(i+1,j+1)
880  u4=ue(i,j+1)
881 
882  s1=(r2-r0)*(z2-z0)
883  s3=(r0-r1)*(z0-z1)
884  s2=(r0-r1)*(z2-z0)
885  s4=(r2-r0)*(z0-z1)
886 
887  u0=(s1*u1+s2*u2+s3*u3+s4*u4)/(s1+s2+s3+s4)
888 
889  bline=u0
890 
891  return
892  end
893 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
894 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
895 
896  real*8 function blint(i,j,r0,z0)
897 
898  include 'double.inc'
899  include 'param.inc'
900  include 'comblc.inc'
901 
902  z1=z(j)
903  z2=z(j+1)
904 
905  r1=r(i)
906  r2=r(i+1)
907 
908  u1=u(i,j)
909  u2=u(i+1,j)
910  u3=u(i+1,j+1)
911  u4=u(i,j+1)
912 
913  s1=(r2-r0)*(z2-z0)
914  s3=(r0-r1)*(z0-z1)
915  s2=(r0-r1)*(z2-z0)
916  s4=(r2-r0)*(z0-z1)
917 
918  u0=(s1*u1+s2*u2+s3*u3+s4*u4)/(s1+s2+s3+s4)
919 
920  blint=u0
921 
922  return
923  end
924 
925 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
926  subroutine grid
927 
928  include 'double.inc'
929  include 'param.inc'
930  include 'comblc.inc'
931 
932 
933  ddr=(rmax-rmin)/dfloat(ni1)
934  ddz=(zmax-zmin)/dfloat(nj1)
935 
936  r(1)=rmin
937  z(1)=zmin
938 
939 
940  do 10 i=2,ni
941 
942  10 r(i)=r(i-1)+ddr
943 
944  do 20 j=2,nj
945 
946  20 z(j)=z(j-1)+ddz
947 
948  return
949  end
950 
951 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
952 
953  subroutine geom
954 
955  include 'double.inc'
956  include 'param.inc'
957  include 'comblc.inc'
958 
959  do 30 i=1,ni1
960 
961  dr(i)=r(i+1)-r(i)
962  r12(i)=(r(i+1)+r(i))*0.5d0
963  if(i.eq.1) go to 30
964  dri(i)=(r(i+1)-r(i-1))*0.5d0
965 
966  30 continue
967 
968  do 40 j=1,nj1
969 
970  dz(j)=z(j+1)-z(j)
971  if(j.eq.1) go to 40
972  dzj(j)=(z(j+1)-z(j-1))*0.5d0
973 
974  40 continue
975 
976  do 2223 ilm=1,nblm
977 
978  do 400 i=1,ni1
979 
980  if(rblm(ilm).lt.r(i+1) .AND. rblm(ilm).ge.r(i) ) then
981  ic =i
982  go to 401
983  endif
984 
985  400 continue
986  401 iblm(ilm)=ic
987 
988  do 500 j=1,nj1
989 
990  if(zblm(ilm).lt.z(j+1) .AND. zblm(ilm).ge.z(j) ) then
991  jc =j
992  go to 501
993  endif
994 
995  500 continue
996  501 jblm(ilm)=jc
997 
998 c... write(6,*)' iblm,jblm,ilim',iblm(ilm),jblm(ilm),ilm
999 
1000  2223 continue
1001 
1002  return
1003  end
1004 
1005 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1006 
1007  subroutine zero(isol)
1008 
1009  include 'double.inc'
1010  include 'param.inc'
1011  include 'comblc.inc'
1012 
1013  integer isol
1014 
1015  !!amu0=0.4d0*pi
1016 
1017  do 100 i=1,ni1
1018 
1019  if(rm0.lt.r(i+1) .AND. rm0.ge.r(i) ) then
1020  ic =i
1021  go to 101
1022  endif
1023 
1024  100 continue
1025  101 imax=ic
1026 
1027  do 200 j=1,nj1
1028 
1029  if(zm0.lt.z(j+1) .AND. zm0.ge.z(j) ) then
1030  jc =j
1031  go to 201
1032  endif
1033 
1034  200 continue
1035  201 jmax=jc
1036 
1037  z1=z(jc)
1038  z2=z(jc+1)
1039 
1040  r1=r(ic)
1041  r2=r(ic+1)
1042 
1043  s1=(r2-rm0)*(z2-zm0)
1044  s3=(rm0-r1)*(zm0-z1)
1045  s2=(rm0-r1)*(z2-zm0)
1046  s4=(r2-rm0)*(zm0-z1)
1047 
1048  ssum=s1+s2+s3+s4
1049 
1050  tok1=amu0*tok*s1/ssum
1051  tok2=amu0*tok*s2/ssum
1052  tok3=amu0*tok*s3/ssum
1053  tok4=amu0*tok*s4/ssum
1054 
1055  do il=1,neqp
1056  right(il)=0.d0
1057  enddo
1058 
1059  do i=1,ni
1060  do j=1,nj
1061  curf(i,j)=0.d0
1062  g(i,j)=0.d0
1063  enddo
1064  enddo
1065 
1066  il=nlin(ic ,jc )
1067  right(il)=-tok1
1068  curf(ic ,jc )=tok1/(dri(ic)*dzj(jc))
1069 
1070  il=nlin(ic+1,jc )
1071  right(il)=-tok2
1072  curf(ic+1,jc )=tok2/(dri(ic+1)*dzj(jc))
1073 
1074  il=nlin(ic+1,jc+1)
1075  right(il)=-tok3
1076  curf(ic+1,jc+1)=tok3/(dri(ic+1)*dzj(jc+1))
1077 
1078  il=nlin(ic ,jc+1)
1079  right(il)=-tok4
1080  curf(ic,jc+1 )=tok4/(dri(ic)*dzj(jc+1))
1081 
1082  ux0=-10d9
1083 
1084  call solve(isol,g)
1085 
1086  call bound
1087 
1088  !!++ call right1 !!spar
1089 
1090  call solve(isol,ui)
1091 
1092 cccccc call psiful
1093 
1094  return
1095  end
1096 
1097 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1098 
1099  subroutine zero_ax(isol)
1100 
1101  include 'double.inc'
1102  include 'param.inc'
1103  include 'comblc.inc'
1104 
1105  integer isol
1106 
1107  do 100 i=1,ni1
1108 
1109  if(rm0.lt.r(i+1) .AND. rm0.ge.r(i) ) then
1110  ic =i
1111  go to 101
1112  endif
1113 
1114  100 continue
1115  101 imax=ic
1116 
1117  do 200 j=1,nj1
1118 
1119  if(zm0.lt.z(j+1) .AND. zm0.ge.z(j) ) then
1120  jc =j
1121  go to 201
1122  endif
1123 
1124  200 continue
1125  201 jmax=jc
1126 
1127 
1128  do il=1,neqp
1129  right(il)=0.d0
1130  enddo
1131 
1132  do i=1,ni
1133  do j=1,nj
1134  curf(i,j)=0.d0
1135  g(i,j)=0.d0
1136  enddo
1137  enddo
1138 
1139  tokint=0.d0
1140  do i=2,ni1
1141  r0=r(i)
1142  do j=2,nj1
1143  z0=z(j)
1144  call cur_map(curden,r0,z0)
1145  il=nlin(i,j)
1146  right(il)=-curden*dri(i)*dzj(j)
1147  curf(i,j)=curden
1148  tokint=tokint+curden*dri(i)*dzj(j)
1149  enddo
1150  enddo
1151  tokint=tokint/amu0
1152 
1153  ux0=-10d9
1154 
1155  call solve(isol,g)
1156 
1157  call bound
1158 
1159  !!++ call right1 !!spar
1160 
1161  call solve(isol,ui)
1162 
1163  return
1164  end
1165 
1166 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1167  subroutine cur_map(curden,rk,zk)
1168 
1169  include 'double.inc'
1170  include 'dim.inc'
1171  include 'compol.inc'
1172 
1173  sqrt(arg)=dsqrt(arg)
1174 
1175  !!amu0=0.4d0*pi
1176 
1177  curden=0.d0
1178  dr0=rk-rm
1179  dz0=zk-zm
1180 
1181  ro0=sqrt(dr0**2+dz0**2)
1182 
1183  tetp=dacos(dr0/ro0)
1184  if(dz0.lt.0.d0) then
1185  tet0=2.d0*pi-tetp
1186  else
1187  tet0=tetp
1188  endif
1189 
1190  if(tet0.lt.teta(1)) tet0=tet0+2.d0*pi
1191  if(tet0.gt.teta(nt)) tet0=tet0-2.d0*pi
1192 
1193  do j=2,nt1
1194 
1195  if(tet0.ge.teta(j) .AND. tet0.lt.teta(j+1)) jc=j
1196 
1197  enddo
1198 
1199  rob1=ro(nr,jc)
1200  rob2=ro(nr,jc+1)
1201 
1202  tet1=teta(jc)
1203  tet2=teta(jc+1)
1204 
1205  rob12=(rob1*(tet2-tet0)+rob2*(tet0-tet1))/(tet2-tet1)
1206 
1207  if(ro0.lt.rob12) then !!!
1208 
1209  do i=1,iplas-1
1210 
1211  ro1=ro(i,jc)
1212  ro2=ro(i,jc+1)
1213  ro3=ro(i+1,jc)
1214  ro4=ro(i+1,jc+1)
1215 
1216  ro12=(ro1*(tet2-tet0)+ro2*(tet0-tet1))/(tet2-tet1)
1217  ro34=(ro3*(tet2-tet0)+ro4*(tet0-tet1))/(tet2-tet1)
1218 
1219  if(ro0.gt.ro12 .AND. ro0.le.ro34) then
1220  ic=i
1221  go to 10
1222  endif
1223 
1224  enddo
1225 
1226  10 continue
1227 
1228  ro1=ro(ic,jc)
1229  ro2=ro(ic+1,jc)
1230  ro3=ro(ic+1,jc+1)
1231  ro4=ro(ic,jc+1)
1232 
1233  u1=cur(ic,jc)
1234  u2=cur(ic+1,jc)
1235  u3=cur(ic+1,jc+1)
1236  u4=cur(ic,jc+1)
1237 
1238  curden=blintr(tet0,ro0,
1239  * tet1,tet2,ro1,ro2,ro3,ro4,u1,u2,u3,u4)
1240 
1241  curden=curden !/amu0
1242  endif !!!
1243 
1244  return
1245  end
1246 
1247 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1248 
1249  real*8 function blintr(tet0,ro0,
1250  * tet1,tet2,ro1,ro2,ro3,ro4,u1,u2,u3,u4)
1251 
1252  include 'double.inc'
1253 
1254  ro14=(ro1*(tet2-tet0)+ro4*(tet0-tet1))/(tet2-tet1)
1255  ro23=(ro2*(tet2-tet0)+ro3*(tet0-tet1))/(tet2-tet1)
1256 
1257  u14=(u1*(tet2-tet0)+u4*(tet0-tet1))/(tet2-tet1)
1258  u23=(u2*(tet2-tet0)+u3*(tet0-tet1))/(tet2-tet1)
1259 
1260  blintr=(u14*(ro23-ro0)+u23*(ro0-ro14))/(ro23-ro14)
1261 
1262  return
1263  end
1264 
nlin
function nlin(i, j)
Definition: Matrix.f:173
bline
real *8 function bline(i, j, r0, z0)
Definition: EQ1_m.f:865
exfmat
subroutine exfmat(rk, zk, tk, nk,
Definition: EQ1_m.f:107
solve
subroutine solve(isol, wdm)
Definition: spider_solver.f:2
zero
subroutine zero(x1, y1, x2, y2, func, err, x, y, izero, ll)
Definition: zero.f90:1
blint
real *8 function blint(i, j, r0, z0)
Definition: EQ1_m.f:896
greeni
REAL *8 function greeni(R, Z, RP, ZP)
Definition: scu.f:484
bound
subroutine bound
Definition: RIG.f:1
blintr
real *8 function blintr(tet0, ro0,
Definition: EQ1_m.f:1249
bint
subroutine bint(X, Y, R0, Z0, r1, z1, F, I)
Definition: scu.f:557
bndmat
subroutine bndmat
Definition: EQ1_m.f:569
flux
subroutine flux(psitok, rk, zk, nk)
Definition: EQ1_m.f:786
cfr_mat
subroutine cfr_mat(rk, zk, tk, nk,
Definition: EQ1_m.f:225
grid
subroutine grid
Definition: EQ1_m.f:926
coil_force
subroutine coil_force
Definition: EQ1_m.f:332
bndint
subroutine bndint(rk, zk, nk)
Definition: EQ1_m.f:655
psiful
subroutine psiful
Definition: Eq_m.f:1120
iter_force_limits
subroutine iter_force_limits(frcz, frcr)
Definition: EQ1_m.f:436
extfil
subroutine extfil(pcequi, ncequi)
Definition: EQ1_m.f:47
ext_fil
subroutine ext_fil(pcequi, ncequi)
Definition: EQ1_m.f:469
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
cur_map
subroutine cur_map(curden, rk, zk)
Definition: EQ1_m.f:1167
grgren
subroutine grgren(R0, Z0, R, Z, dGdr, dGdz)
Definition: scu.f:519
blin
real *8 function blin(i, j, r0, z0)
Definition: EQ1_m.f:835
psi_fil
subroutine psi_fil(psicon, ncequi)
Definition: EQ1_m.f:3
rdexf
subroutine rdexf(ncequi)
Definition: EQ1_m.f:509
geom
subroutine geom
Definition: EQ1_m.f:953
zero_ax
subroutine zero_ax(isol)
Definition: EQ1_m.f:1099