ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
B_rig.f
Go to the documentation of this file.
1 !! file contains following routines:
2 !!
3 !! funp
4 !! funf
5 !! rightg
6 !! rigbon
7 
8 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
9 
10  subroutine rightg
11 
12 !! riht-hand side for problem L(g)=J
13 
14  include 'double.inc'
15  include 'dim.inc'
16  include 'compol.inc'
17 
18  common /comabw/ alf0p,alf1p,alf2p,bet0f,bet1f,bet2f
19  * ,alw0p,alw1p,alw2p
20  common
21  * /c_kpr/kpr
22 
23  common/com_flag/kastr
24  common /com_0st/ key_0st,key_prs
25  save cur_mu
26 
27  !!amu0=0.4d0*pi
28 
29  tokp=0.d0
30  tokff=0.d0
31  tokpp=0.d0
32  tokww=0.d0
33 
34  do 1 il=1,neqp
35  right(il)=0.d0
36  1 continue
37 
38 
39  if(kstep.eq.0 .OR. (ngav.eq.0 .AnD. kastr.eq.0)) then
40 
41  !!!!!!!!!!!!!!!!!!!!
42  do 2 i=1,iplas
43  do 2 j=1,nt
44  cur(i,j)=0.d0
45  2 continue
46 
47 !! central point
48 
49  psn=psin(1,2)
50  r0=r(1,2)
51 
52  curp=tabp(psn)
53  curf=tabf(psn)
54  curw=tabw(psn)
55 
56  dpdpsi(1)=curp
57  dfdpsi(1)=curf
58  dwdpsi(1)=curw
59 
60  curcen=r0*curp+curf/r0+curw*r0**3
61 
62  sqcen=0.d0
63 
64  do 10 j=2,nt1
65 
66  sqcen=sqcen+sq1(1,j)+sq4(1,j)
67 
68  10 continue
69 
70  tokp=curcen*sqcen
71  tokff=curf*sqcen/r0
72  tokpp=curp*sqcen*r0
73  tokww=curw*sqcen*r0**3
74 
75  il=numlin(1,1,nr,nt)
76 
77  right(il)=curcen*sqcen
78 
79 !!!!!!!!!!!!! regular points
80 
81  do 20 i=2,iplas !1
82  do 25 j=2,nt1
83 
84  if(i.ne.iplas) then
85 
86  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
87  r0=r(i,j)
88  psn=psin(i,j)
89 
90  else
91 
92  sqk=sq2(i-1,j)+sq3(i-1,j-1)
93  r0=r(i,j)
94  psn=psin(i,j)
95  r0=(3.d0*r(i,j)+r(i-1,j))*0.25d0
96  psn=(3.d0*psin(i,j)+psin(i-1,j))*0.25d0
97 
98  endif
99 
100  curp=tabp(psn)
101  curf=tabf(psn)
102  curw=tabw(psn)
103 
104  cur(i,j)=r0*curp+curf/r0+curw*r0**3
105 
106  tokp=tokp+cur(i,j)*sqk
107  tokff=tokff+curf*sqk/r0
108  tokpp=tokpp+curp*sqk*r0
109  tokww=tokww+curw*sqk*r0**3
110 
111  il=numlin(i,j,nr,nt)
112 
113  right(il)=cur(i,j)*sqk
114 
115  25 continue
116 
117  dpdpsi(i)=curp
118  dfdpsi(i)=curf
119  dwdpsi(i)=curw
120 
121  20 continue
122 
123  dpdpsi(iplas)=tabp(0.d0) !
124  dfdpsi(iplas)=tabf(0.d0) !
125  dwdpsi(iplas)=tabw(0.d0) !
126 
127  cnor=amu0*tok/tokp
128 
129 c write(*,*) 'cnor,tok,tokp'
130 c write(*,*) cnor,tok,tokp
131 
132  if(key_0st.eq.1) cnor=1.d0
133 
134  tok_pl=0.d0
135 
136 
137  do 30 j=1,nt
138 
139  cur(1,j)=curcen*cnor
140 
141 
142  30 continue
143 
144  do 40 i=2,iplas
145  do 40 j=2,nt1
146 
147  cur(i,j)=cur(i,j)*cnor
148 
149 
150  40 continue
151 
152  do 45 i=1,iplas
153 
154  cur(i,1)=cur(i,nt1)
155  cur(i,nt)=cur(i,2)
156 
157  45 continue
158 
159  do 50 il=1,neqpla
160 
161  right(il)=right(il)*cnor
162 
163 !test psi=r**2
164 
165 ! right(il)=0.d0
166 
167 !test psi=r**2
168 
169  50 continue
170 
171  do 55 i=1,iplas
172 
173  dpdpsi(i)=dpdpsi(i)*cnor !
174  dfdpsi(i)=dfdpsi(i)*cnor !
175  dwdpsi(i)=dwdpsi(i)*cnor !
176 
177  55 continue
178 
179  cur_mu=tokp*cnor
180  tokp=tokp/amu0*cnor
181 
182  !!!!!!!!!!!!!!!!!!!!!!!!
183  else !!!!!!!!!!!!!!!!
184  !!!!!!!!!!!!!!!!!!!!!!!!
185 
186  if(kastr.eq.1 .AND. key_prs.eq.1 .AND. erru.lt.5.d-3) then
187  call pres_d_psi !
188  endif
189 
190  !if(ngav.ne.0 ) then
191  call procof(1,cur_mu)
192  !endif
193 
194  dfdpsi_sur=-(fvac**2-f(iplas-1)**2)/psia(iplas-1)/psim
195 
196 
197 !!!!!!!!! central point
198 
199  curcen=rm*dpdpsi(1)+dfdpsi(1)/rm+dwdpsi(1)*rm**3
200 
201  sqcen=0.d0
202 
203  do 100 j=2,nt1
204 
205  sqcen=sqcen+sq1(1,j)+sq4(1,j)
206 
207  100 continue
208 
209  tokp=curcen*sqcen
210 
211  il=numlin(1,1,nr,nt)
212 
213  right(il)=curcen*sqcen
214 
215  tok_pl=0.d0
216 
217 
218  do 300 j=1,nt
219 
220  cur(1,j)=curcen
221 
222 
223  300 continue
224 
225 !! regular points
226 
227 
228  do 200 i=2,iplas !-1
229  do 200 j=2,nt1
230 
231  if(i.ne.iplas) then
232 
233  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
234  r0=r(i,j)
235 
236  curp=dpdpsi(i)
237  curf=dfdpsi(i)
238  curw=dwdpsi(i)
239 
240  else
241 
242  sqk=sq2(i-1,j)+sq3(i-1,j-1)
243  ! r0=r(i,j)
244  r0=(3.d0*r(i,j)+r(i-1,j))*0.25d0
245 
246  curp=dpdpsi(i)
247  curf=dfdpsi(i) !+dfdpsi_sur
248  curw=dwdpsi(i)
249 
250  endif
251 
252  cur(i,j)=r0*curp+curf/r0+curw*r0**3
253 
254  tokp=tokp+cur(i,j)*sqk
255 
256  il=numlin(i,j,nr,nt)
257 
258  right(il)=cur(i,j)*sqk
259 
260  200 continue
261 
262 c write(6,*) 'rightp:tokp=',tokp
263 
264  do 450 i=1,iplas
265 
266  cur(i,1)=cur(i,nt1)
267  cur(i,nt)=cur(i,2)
268 
269 
270  450 continue
271  cnor=cur_mu/tokp
272 
273  tokp=tokp/amu0
274 
275  if(ngav.eq.0 ) then
276 
277  !tokp=tok
278  tokp=cur_mu/amu0
279  do il=1,neqpla
280  right(il)=right(il)*cnor
281  enddo
282  do i=1,iplas
283  do j=1,nt
284  cur(i,j)=cur(i,j)*cnor
285  enddo
286  enddo
287 
288  endif
289 
290  if(kpr.eq.1) then
291  write(*,*) 'cnor=',cnor
292  endif
293  endif
294 
295 
296 
297  return
298  end
299 
300 
301 
302 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
303 
304  subroutine rigbon
305 
306  include 'double.inc'
307  include 'dim.inc'
308  include 'compol.inc'
309 
310  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
311  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
312 
313  real*8 psib(ntp)
314 
315  do 20 j=2,nt1
316 
317  !psib(j)=psi(iplas,j)
318  psib(j)=psip
319 
320 !test psi=r*2
321 ! psib(j)=z(iplas,j) !**2
322 !test psi=r*2
323 
324  20 continue
325 
326  !open(3,file='boncon.wr')
327  !read(3,*) (psib(j),j=2,nt1)
328  !close(3)
329 
330  psib(1)=psib(nt1)
331  psib(nt)=psib(2)
332  !
333  !!
334  ! !!!
335  i=iplas-1 !!!!!!!!!!!!!!!!!!!!!
336  ! !!!
337  !!
338  !
339  do 100 j=2,nt1
340 
341  a7=a24(i,j-1)
342  a8=a34(i,j-1)+a12(i,j)
343  a9=a13(i,j)
344 
345  il=numlin(i,j,nr,nt)
346 
347  right(il)=right(il)-(a7*psib(j-1)+a8*psib(j)+a9*psib(j+1))
348 
349  100 continue
350 
351  return
352  end
353 
354 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
355 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
356 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
357 
358  subroutine qst_bq
359 
360  include 'double.inc'
361  include 'dim.inc'
362  parameter(nshp=ntp+1)
363  include 'compol.inc'
364 
365  common/efites/ fcefit,rc_efit,iefit
366  dimension xs(nshp),ys(nshp),fun(nshp),dp(5)
367 
368  sqrt(xx)=dsqrt(xx)
369 
370  nsh=1
371  xs(nsh)=r(1,1)
372  ys(nsh)=z(1,1)
373  fun(nsh)=psin(1,1)
374 
375  do 100 j=2,nt1
376 
377  nsh=nsh+1
378  xs(nsh)=r(2,j)
379  ys(nsh)=z(2,j)
380  fun(nsh)=psin(2,j)
381 
382  100 continue
383 
384  call deriv5(xs,ys,fun,nsh,5,dp)
385 
386  drr=dp(3)
387  drz=dp(4)
388  dzz=dp(5)
389 
390  tg2a=2.d0*drz/(drr-dzz)
391  cos2a=1.d0/sqrt(1.d0+tg2a**2)
392  sin2a=cos2a*tg2a
393 
394  dxx=0.5d0*(drr+dzz)+0.5d0*cos2a*(drr-dzz)+sin2a*drz
395  dyy=0.5d0*(drr+dzz)+0.5d0*cos2a*(dzz-drr)-sin2a*drz
396 
397 c write(6,*) 'drr,dzz,drz',drr,dzz,drz
398 c write(6,*) 'dxx,dyy',dxx,dyy
399 
400  bfcen=qcen*dsqrt(dxx*dyy)*(psim-psip)
401 
402  fcen=bfcen*rm
403  fcen=fcefit
404 
405  dpsi=(psia(1)-psia(2))
406  ps14=1.d0-0.25d0*dpsi
407 
408  ffp=tabf(ps14)
409 
410  f(1)=sqrt(fcen**2-ffp*cnor*(psim-psip)*dpsi)
411 
412  do 200 i=2,iplas
413 
414  if(i.ne.iplas) then
415 
416  pspl=0.5d0*(psia(i)+psia(i+1))*(psim-psip)
417  psmn=0.5d0*(psia(i)+psia(i-1))*(psim-psip)
418 
419  else
420 
421  pspl=0.5d0*psia(i)*(psim-psip)
422  psmn=0.5d0*psia(i-1)*(psim-psip)
423 
424  endif
425 
426  psn=psia(i)
427 
428  !ffp=funf(psn,alf0p,alf1p,alf2p,bet0f,bet1f,bet2f)
429  ffp=dfdpsi(i)
430 
431  f(i)=sqrt(f(i-1)**2+2.d0*ffp*(pspl-psmn))
432 
433  200 continue
434 
435  fvac=f(iplas)
436 
437  !open(1,file='q.pr')
438 
439  flucfm=0.d0
440 
441  do 300 i=1,iplas
442 
443  flucfi=0.d0
444 
445  do 310 j=2,nt1
446 
447  if(i.ne.iplas) then
448 
449  r1=r(i,j)
450  r2=r(i+1,j)
451  r3=r(i+1,j+1)
452  r4=r(i,j+1)
453  r0=(r1+r2+r3+r4)*0.25d0
454  flucfi=flucfi+s(i,j)*f(i)/r0
455 
456  else
457 
458  r0=r(i,j)
459  flucfi=flucfi+(sq3(i-1,j-1)+sq2(i-1,j))*f(i)/r0
460 
461  endif
462 
463  310 continue
464 
465  if(i.ne.iplas) then
466 
467  q2pi=-flucfi/((psia(i+1)-psia(i))*(psim-psip))
468  flucfm=flucfm+flucfi
469 
470  !write(1,*) 0.5d0*(psia(i)+psia(i+1)),0.5d0*q(i)/pi,i
471  else
472 
473  q2pi=-2.d0*flucfi/((psia(i)-psia(i-1))*(psim-psip))
474 
475  !write(1,*) psia(i),0.5d0*q(i)/pi,i
476  endif
477 
478  !q(i)=0.5d0*q2pi/pi
479  q(i)=q2pi
480 
481  300 continue
482 
483  !close(1)
484 
485  return
486  end
487 
488 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
489 
490 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
491 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
492 
493 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
494 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
495 
496 
497 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
498 
499  subroutine psib_pla(pspl_av)
500 
501  include 'double.inc'
502  include 'dim.inc'
503  include 'parcur.inc'
504  include 'compol.inc'
505 
506  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
507  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
508 
509  common/com_bgr/ binadg(ntp,ntp),dgdn(ntp)
510  dimension pspl_b(ntp)
511 
512  i=iplas !!!!!!!!!!!!!!!!!!!!!!1
513 
514  do 10 j=2,nt1
515 
516  a1=a13(i-1,j-1)
517  a2=a34(i-1,j-1)+a12(i-1,j)
518  a3=a24(i-1,j)
519 
520  g1=psi(i-1,j-1)
521  g2=psi(i-1,j)
522  g3=psi(i-1,j+1)
523 
524  sqk=sq2(i-1,j)+sq3(i-1,j-1)
525  dltk=(dlt(i,j-1)+dlt(i,j))*0.5d0
526  dgdnl=-(a1*g1+a2*g2+a3*g3) + cur(i,j)*sqk !
527  dgdn(j)=dgdnl/dltk
528 
529  10 continue
530 
531  dgdn(1)=dgdn(nt1)
532  dgdn(nt)=dgdn(2)
533 
534  do l=2,nt1
535 
536  psb=0.d0
537 
538  do jb=2,nt1
539 
540  psb=psb+binadg(jb,l)*(dgdn(jb)+dgdn(jb+1))*0.5d0
541 
542  enddo
543 
544  pspl_b(l)=psb
545 
546  enddo
547 
548  pspl_b(1)=pspl_b(nt1)
549  pspl_b(nt)=pspl_b(2)
550 
551  pspl_av = avr_bnd(pspl_b)
552 
553  return
554  end
555 
556 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
557 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
558 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
559 
560 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
561 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
562  subroutine psloop_e(rlop,zlop,psilop,nlop)
563 
564  include 'double.inc'
565  include 'parrc1.inc'
566  !include 'comrec.inc'
567  common /comind/ ni,nj,ni1,nj1,ni2,nj2,nbnd,nkin,nkout
568  common /comrz/ x(nip),y(njp),dx(nip),dy(njp),
569  + dxi(nip),dyj(njp),x12(nip)
570 
571  common /compot/ u(nip,njp),ue(nip,njp),un(nip,njp),
572  + ui(nip,njp),g(nip,njp),
573  + ux0,ux1,ux2,up,um,xm,ym,
574  + xx0,yx0,xx1,yx1,xx2,yx2,imax,jmax,
575  + ix1,jx1,ix2,jx2,
576  + xx10,yx10,xx20,yx20,xm0,ym0,
577  + psi_bon
578  dimension rlop(*),zlop(*),psilop(*)
579 
580  ddx=x(2)-x(1)
581  ddy=y(2)-y(1)
582 
583  do i=1,nlop
584 
585  r0=rlop(i)
586  z0=zlop(i)
587 
588  ic=(r0-x(1))/ddx+1
589  jc=(z0-y(1))/ddy+1
590 
591  r1=x(ic)
592  r2=x(ic+1)
593 
594  z1=y(jc)
595  z2=y(jc+1)
596 
597  u1=ue(ic,jc)
598  u2=ue(ic+1,jc)
599  u3=ue(ic+1,jc+1)
600  u4=ue(ic,jc+1)
601 
602  ! psilop(i)=blin(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4 )
603  psilop(i)=blin_(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4 )
604 
605  enddo
606 
607  return
608  end
609 
610 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
611 
612 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
613 
614  subroutine psib_ext(psex_av)
615 
616  include 'double.inc'
617  include 'dim.inc'
618  !include 'parcur.inc'
619  include 'compol.inc'
620 
621  dimension psex_b(ntp),rbon(ntp),zbon(ntp)
622 
623  i=iplas !!!!!!!!!!!!!!!!!!!!!!1
624 
625  do j=1,nt
626  rbon(j)=r(i,j)
627  zbon(j)=z(i,j)
628  enddo
629 
630  call psloop_e(rbon,zbon,psex_b,nt)
631 
632  psex_b(1)=psex_b(nt1)
633  psex_b(nt)=psex_b(2)
634 
635  psex_av = avr_bnd(psex_b)
636 
637  return
638  end
639 
640 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
641 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
642 
643  subroutine field_c
644 
645  include 'double.inc'
646  include 'dim.inc'
647  include 'compol.inc'
648 
649  common /com_mf/ bpol(nrp,ntp),btor(nrp,ntp),
650  & br(nrp,ntp),bz(nrp,ntp),
651  & rc(nrp,ntp),zc(nrp,ntp)
652 
653  common /com_trap/ trap(nrp)
654 
655  dimension btot(nrp,ntp)
656 
657 !magnetic field at cells
658 
659  do i=1,iplas-1
660  do j=2,nt-1
661 
662  bp_ij=(psi(i+1,j)-psi(i,j))/st(i,j)
663  bp_ij1=(psi(i+1,j+1)-psi(i,j+1))/st(i,j+1)
664 
665  if(i.ne.1) then
666  bp_2= (
667  & bp_ij**2*( vol1(i,j)/sin1(i,j)+vol2(i,j)/sin2(i,j) )+
668  & bp_ij1**2*( vol3(i,j)/sin3(i,j)+vol4(i,j)/sin4(i,j) )
669  & )/vol(i,j)
670  else
671  bp_2= (
672  & bp_ij**2*( vol2(i,j)/sin2(i,j) )+
673  & bp_ij1**2*( vol3(i,j)/sin3(i,j) )
674  & )/vol(i,j)
675  endif
676  bpol_ij=dsqrt(bp_2)
677  btor_ij=f(i)*s(i,j)/vol(i,j)
678 
679  bpol(i,j)=bpol_ij
680  btor(i,j)=btor_ij
681  btot_ij=dsqrt(bp_2+btor_ij**2)
682 
683  r1=r(i,j)
684  r2=r(i+1,j)
685  r3=r(i+1,j+1)
686  r4=r(i,j+1)
687  r0=(r1+r2+r3+r4)*0.25d0
688 
689  r12=(r1+r2)*0.5d0
690  r23=(r3+r2)*0.5d0
691  r34=(r3+r4)*0.5d0
692  r14=(r1+r4)*0.5d0
693 
694  z1=z(i,j)
695  z2=z(i+1,j)
696  z3=z(i+1,j+1)
697  z4=z(i,j+1)
698  z0=(z1+z2+z3+z4)*0.25d0
699 
700  z12=(z1+z2)*0.5d0
701  z23=(z3+z2)*0.5d0
702  z34=(z3+z4)*0.5d0
703  z14=(z1+z4)*0.5d0
704 
705  dr=r34-r12
706  dz=z34-z12
707  drz= dsqrt(dr**2+dz**2)
708  dr=dr/drz
709  dz=dz/drz
710 
711  br(i,j)=bpol_ij*dr
712  bz(i,j)=bpol_ij*dz
713 
714  rc(i,j)=r0
715  zc(i,j)=z0
716 
717  enddo
718  enddo
719 
720  do i=1,iplas-1
721  bpol(i,1)=bpol(i,nt-1)
722  btor(i,1)=btor(i,nt-1)
723  br(i,1)=br(i,nt-1)
724  bz(i,1)=bz(i,nt-1)
725  rc(i,1)=rc(i,nt-1)
726  zc(i,1)=zc(i,nt-1)
727  enddo
728 
729  do i=1,iplas-1
730  bmax=0.d0
731  svol=0.d0
732  do j=2,nt-1
733  b_ij=dsqrt(bpol(i,j)**2+btor(i,j)**2)
734  bmax=dmax1(bmax,b_ij)
735  btot(i,j)=b_ij
736  svol=svol+vol(i,j)
737  enddo
738 
739  avr_v=0.d0
740  do j=2,nt-1
741  bnor=btot(i,j)/bmax
742  avr_v=avr_v+
743  & (b0ax/btot(i,j))**2
744  & *( 1.d0-dsqrt(1.d0-bnor)*(1.d0+.5d0*bnor) )*vol(i,j)
745  enddo
746 
747  trap(i)=avr_v/svol
748 
749  enddo
750 
751  write(fname,'(a,a)') path(1:kname),'fields.wr'
752  open(1,file=fname)
753  !open(1,file='fields.wr',form='formatted')
754  write(1,*) iplas-1,nt-1
755  write(1,*) ((rc(i,j),i=1,iplas-1),j=1,nt-1)
756  write(1,*) ((zc(i,j),i=1,iplas-1),j=1,nt-1)
757  write(1,*) ((br(i,j),i=1,iplas-1),j=1,nt-1)
758  write(1,*) ((bz(i,j),i=1,iplas-1),j=1,nt-1)
759  write(1,*) ((btor(i,j),i=1,iplas-1),j=1,nt-1)
760  write(1,*) ((bpol(i,j),i=1,iplas-1),j=1,nt-1)
761  write(1,*) (trap(i),i=1,iplas-1)
762  close(1)
763 
764  return
765  end
766 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
767 
768 
psloop_e
subroutine psloop_e(rlop, zlop, psilop, nlop)
Definition: B_rig.f:562
qst_bq
subroutine qst_bq
Definition: B_rig.f:358
psib_ext
subroutine psib_ext(psex_av)
Definition: B_rig.f:614
fun
subroutine fun(X, F)
Definition: Ev2.f:10
procof
subroutine procof(icq, cur_mu)
Definition: B_metric_O.f:491
helena21_mod::dpdpsi
real(r8) function dpdpsi(psi_n)
Definition: helena21_mod.f90:122
rightg
subroutine rightg
Definition: B_rig.f:10
numlin
function numlin(i, j, nr, nt)
Definition: com_sub.f:1040
field_c
subroutine field_c
Definition: B_rig.f:643
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
avr_bnd
real *8 function avr_bnd(arr)
Definition: com_sub.f:221
psib_pla
subroutine psib_pla(pspl_av)
Definition: B_rig.f:499
pres_d_psi
subroutine pres_d_psi
Definition: interfac.f:2772
blin_
real *8 function blin_(r0, z0, r1, r2, z1, z2, u1, u2, u3, u4)
Definition: scu.f:857
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
rigbon
subroutine rigbon
Definition: B_rig.f:304