ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
B_eqb.f
Go to the documentation of this file.
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 
3  subroutine eqb( alf0,alf1,alf2, bet0,bet1,bet2, alw0,alw1,alw2,
4  * betplx, i_betp,
5  * keyctr,nstep,platok,rax,zax,b0cen,r0cen,psax,igdf,
6  * n_tht, n_psi, epsro, nurs, i_eqdsk,i_bsh,
7  * psi_bnd,psi0_bnd)
8 
9  include 'double.inc'
10  include 'parurs.inc'
11  parameter(nursp4=nursp+4,nursp6=nursp4*6)
12 
13  include 'dim.inc'
14  include 'compol.inc'
15 
16  common /combsh/ rm0,zm0,rc0,zc0,asp0,el_up,el_lw,tr_up,tr_lw,nbsh
17  common /comabw/ alf0p,alf1p,alf2p, bet0f,bet1f,bet2f
18  * ,alw0p,alw1p,alw2p
19 
20 c -----------------------------------------------------------------
21 
22  integer nstep, ngav
23  real*8 alf0,alf1,alf2, bet0,bet1,bet2, alw0,alw1,alw2
24  real time_beg,time_end, time_b,time_e, dtim1,dtim2,dtim3
25  real*8 q_giv(nrp),q_0(nrp)
26 
27  dimension pstab(nursp), qtab(nursp), psian(nrp)
28  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
29  real*8 cwk(4)
30  dimension bfjf(nrp)
31  common
32  * /c_kpr/kpr
33 
34  abs(xx) = dabs(xx)
35  sqrt(xx) = dsqrt(xx)
36 
37 c--------------------------------------------------------------------
38 
39  nimax=250
40  !NiMax=10
41 
42 c--------------------------------------------------------------------
43 
44  !call cpu_time(time_beg)
45 
46  !!amu0=0.4d0*pi
47  ! if(keyctr.ne.0) return
48  if(keyctr.ge.100) return
49 
50  nt = n_tht
51  iplas = n_psi
52  ngav = keyctr
53  tok = platok
54  nbsh=i_bsh
55  ! qcen = qax
56  b0ax = b0cen
57  r0ax = r0cen
58  psiax = psax
59  if(ngav.eq.-1) then
60  psibon=psi_bnd
61  else
62  psi_eav=psi_bnd
63  endif
64  psibon0=psi0_bnd
65 
66  kstep=nstep
67 
68 c...input initial data
69 
70  rm0 = rax
71  zm0 = zax
72 
73  alf0p = alf0
74  alf1p = alf1
75  alf2p = alf2
76 
77  bet0f = bet0
78  bet1f = bet1
79  bet2f = bet2
80 
81  alw0p = alw0
82  alw1p = alw1
83  alw2p = alw2
84 c--------------------------------------------------------------------
85 
86  if(i_eqdsk.eq.1) then
87  nbsh=1
88  endif
89 
90  nr = iplas
91  nr1 = nr-1
92  nt1 = nt-1
93  nr2 = nr-2
94  nt2 = nt-2
95 
96  iplas1 = iplas-1
97 
98  itrmax = 50
99  nitmax = 5
100  nitdel = 10
101  nitbeg = 5 !+ngav*10000
102 c--------------------------------------------------------------------
103 
104  if(nstep.eq.0) then
105  if(nbsh.eq.0) then
106  write(fname,'(a,a)') path(1:kname),'inpol.dat'
107  open(1,file=fname)
108  !open(1,file='inpol.dat')
109  read(1,*) dummy_nbsh
110  read(1,*) rc0
111  read(1,*) zc0
112  read(1,*) asp0
113  read(1,*) el_up
114  read(1,*) el_lw
115  read(1,*) tr_up
116  read(1,*) tr_lw
117  close(1) ! for "inpol.dat" reading
118  endif
119 
120 ! if( nbsh.eq.0 ) then
121 ! write(* ,*) '********************************'
122 ! write(* ,*) 'nbsh =',nbsh
123 ! write(* ,*) '"inpol.dat" file has been read: '
124 ! write(* ,*) '--------------------------------'
125 ! write(* ,*) ' '
126 ! write(17,*) '********************************'
127 ! write(17,*) 'nbsh =',nbsh
128 ! write(17,*) '"inpol.dat" file has been read: '
129 ! write(17,*) '--------------------------------'
130 ! write(17,*) ' '
131 ! write(17,*) 'rc0 =',rc0
132 ! write(17,*) 'zc0 =',zc0
133 ! write(17,*) 'asp0 =',asp0
134 ! write(17,*) 'el_up =',el_up
135 ! write(17,*) 'el_lw =',el_lw
136 ! write(17,*) 'tr_up =',tr_up
137 ! write(17,*) 'tr_lw =',tr_lw
138 ! write(17,*) ' '
139 ! else
140 ! write(* ,*) '************************************'
141 ! write(* ,*) 'nbsh =',nbsh
142 ! write(* ,*) '"inpol.dat" file has not been read. '
143 ! write(* ,*) '------------------------------------'
144 ! write(* ,*) ' '
145 ! write(17,*) '************************************'
146 ! write(17,*) 'nbsh =',nbsh
147 ! write(17,*) '"inpol.dat" file has not been read. '
148 ! write(17,*) '------------------------------------'
149 ! write(17,*) ' '
150 ! end if
151 
152  igrdf = igdf
153 
154  if(igdf.eq.2) igrdf=1
155 
156  !!!++call arc_x_bnd(nt)
157  if(nbsh.eq.0.or.nbsh.eq.1) then
158  call grid_0(igrdf,nstep)
159  call taburs(0,1.d0,nurs)
160  elseif(nbsh.eq.-1) then
161  call grid_b(igrdf,nstep)
162  !call wrb
163  elseif(nbsh.eq.-2) then
164  call grid_p0(igrdf,nstep)
165  call taburs(0,1.d0,nurs)
166  endif
167 
168  do i=1,iplas
169  dpdpsi(i)=tabp(psia(i))
170  dfdpsi(i)=tabf(psia(i))
171  enddo
172 
173  if(ngav.eq.0) call presol(i_betp,betplx,betpol)
174 
175  else !nstep.ne.0
176 
177  !!!++call taburs(0,1.d0,nurs)
178  !!!++call arc_x_bnd(nt)
179  if(nbsh.eq.0.or.nbsh.eq.1) then
180  call grid_1(igrdf,nstep)
181  elseif(nbsh.eq.-1) then
182  call grid_b1(igrdf,nstep)
183  !call wrb
184  elseif(nbsh.eq.-2) then
185  call grid_p1(igrdf,nstep)
186  endif
187 
188 ! do i=1,iplas
189 ! dpdpsi(i)=tabp(psia(i))
190 ! dfdpsi(i)=tabf(psia(i))
191 ! enddo
192  endif !nstep.eq.0
193 
194  if(ngav.lt.0) then
195  call bongri
196  endif
197 !&&&&&&&&&&&&
198  !if(ngav.eq.-2.AND.nbsh.eq.-1) then
199  if(nbsh.eq.-1) then
200  call psib_ext(psi_eav)
201  endif
202 !&&&&&&&&&&&&
203 
204 c--------------------------------------------------------------------
205 
206 ! call taburs(0, 1.d0, nurs)
207 
208  cnor=1.d0
209 
210 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
212  go to 639 !!<<<<<<+**
213 
214  if(ngav.gt.0) then
215 
216 c------------------------------------------------
217  write(fname,'(a,a)') path(1:kname),'tab_q.dat'
218  open(1,file=fname)
219  !open(1,file='tab_q.dat')
220  read(1,*) nutab
221  do i=1,nutab
222  read(1,*) pstab(i), qtab(i)
223  enddo
224  close(1)
225 c------------------------------------------------
226 
227  do i=1,nutab
228  pstab(i)=pstab(i)/pstab(nutab)
229  enddo
230 
231  call e01baf(nutab,pstab,qtab,rrk,cck,
232  * nutab+4,wrk,6*nutab+16,ifail)
233 
234  2257 continue
235 
236  do i=1,iplas1
237 
238  psia05=1.d0-0.5d0*(psia(i)+psia(i+1))
239 
240  call e02bcf(nutab+4,rrk,cck,psia05,0,cwk,ifail)
241 
242  q_giv(i)=cwk(1)*2.d0*pi
243  q_0(i)=q(i)
244 
245  enddo
246 
247  q_giv0=qtab(1)
248  q_giv(iplas)=qtab(nutab)*2.d0*pi
249 
250  do i=1,iplas
251  q(i) = q_giv(i)
252  enddo
253 
254  if(igdf.eq.2) then
255 
256  do i=1,iplas
257  psian(i)=psia(i)
258  enddo
259 
260  call grdef(igdf)
261 
262  errpsa=0.d0
263 
264  do i=1,iplas
265  errpsn=dabs(psian(i)-psia(i))
266  errpsa=dmax1(errpsa,errpsn)
267  enddo
268 
269  if(errpsa.gt.1.d-9) go to 2257
270 
271  endif
272 
273  qax = q_giv0
274 
275  do i=1,iplas
276  dpdpsi(i)=tabp(psia(i))
277  dfdpsi(i)=tabf(psia(i))
278  enddo
279 
280  endif
281 
282 !!!!!!!!!!!!!!!!!!!!
283  639 continue
284 !!!!!!!!!!!!!!!!!!!!
285 
286 
287 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
288 
289  iter = 0
290  itin = 0
291  imov = 1
292 
293  1000 continue
294 
295  iter = iter+1
296  itin = itin+1
297 
298  !write(6,*)'iter=',iter,itin
299 
300  igrdf = igdf
301 
302  !if(igdf.eq.2 .and. itin.lt.4) igrdf=1
303 
304  !call cpu_time(time_b)
305 ch4astra call metric
306  call metrix
307  !call cpu_time(time_e)
308  !write(*,*) '***metric:time=',time_e-time_b
309  call matcof
310  !call cpu_time(time_b)
311  !write(*,*) '***matcof:time=',time_e-time_b
312  call matpla
313  !call cpu_time(time_e)
314  !write(*,*) '***matpla:time=',time_e-time_b
315  call rightg
316  !call cpu_time(time_b)
317  !w!rite(*,*) '***rightg:time=',time_e-time_b
318 
319  if(i_betp.eq.1) then
320  if(iter.gt.4) call skbetp(betplx,betpol)
321  endif
322 
323  !call bt_pol(betpol)
324  !write(6,*)'betpol',betpol,'*********'
325  !call rigbon
326 
327  !if(ngav.le.0 .and. igrdf.ge.2) then
328  if(ngav.le.0 ) then
329  call qst_b
330  endif
331 
332  call grdef(igrdf)
333  call solint(imov)
334 
335  !write(6,*)'solve(psi)'
336  !call wrb
337  !pause 'wrd after solve'
338 
339  call spider_remesh(erro,errpsi,imov)
340  !call cpu_time(time_e)
341  !write(*,*) '***solve:time=',time_e-time_b
342 
343 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
344 
345  ! if(ngav.eq.1) then
346  ! call wrb
347  ! pause 'wrb after regrid '
348  ! endif
349 
350 !!!!!!!!!!!! accurasy test parameters !!!!!!!!!!!!!!!
351 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
352 
353  errod = 0.5d0*erro/(dabs(z(iplas,2)-zm)+dabs(r(iplas,2)-rm))
354  erru=erro
355  !write(6,*) 'errod,epsro',errod,epsro
356  if(errod.lt.epsro) go to 2000
357  if(itin.ge.nimax) then
358  write(*,*) 'SPIDER: MAX NUMBER OF ITERATIONS IS EXCEEDED'
359  write(*,*) 'ERROD=',errod
360  write(*,*) 'ITER=',itin
361  go to 2000
362  endif
363 
364  go to 1000
365 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
366 
367  2000 continue
368 ! write(*,*) 'OK'
369 
370  if(ngav.le.0) then
371 
372  call qst_b
373  ! write(6,*)'qst'
374 ! open(1,file='ddps0.pr')
375 ! do 567 i=1,iplas
376 ! 567 write(1,*) 'dfdpsi(i), f(i), i', dfdpsi(i),f(i),i
377 ! close(1)
378 
379  psiax=psim
380 
381  !call procof(0)
382 
383  endif
384 c--------------------------------------------------------------------
385 
386  platok = tokp
387  rax = rm
388  zax = zm
389  psax = psim
390 
391  if(ngav.le.0) then
392  psax = psim
393  psipla = psim-psip
394  endif
395  if(ngav.eq.-3) then
396  psi_bnd=psi_eav
397  endif
398  !fvac = f(iplas)
399 
400  if(ngav.gt.0) qax = q_giv0
401 
402  if(ngav.lt.0) call retab_l
403  if(ngav.gt.0) call retab_p
404 
405  call bt_pol(betpol)
406  call bt_tot(bettot)
407 
408  if(kpr.eq.1) then
409  write(* ,*) '*******************************************'
410  write(* ,*) 'Iterations have been converged for epsro =',epsro
411  write(* ,*) 'Achieved accuracy ................ errod =',errod
412  write(* ,*) 'Number of iterations ............. iter =',iter
413  write(* ,*) 'Magnetic axis coordinates ........ rm =', rm
414  write(* ,*) ' zm =', zm
415  write(* ,*) 'Plasma current ................... tokp =', tokp
416  write(* ,*) ' cnor =', cnor
417  write(* ,*) 'Magnetic axis PSI value .......... psim =', psim
418  write(* ,*) 'Poloidal beta value ............. betpol =', betpol
419  write(* ,*) 'Total beta value ............. bettot =', bettot
420  write(* ,*) 'Plas.boun.poloidal current value... fvac =', fvac
421  !write(* ,*) 'Magnetic axis q value ............. qax =', qax
422  write(* ,*) '*******************************************'
423  endif
424 ! write(17,*) '*******************************************'
425 ! write(17,*) 'Iterations have been converged for epsro =',epsro
426 ! write(17,*) 'Achieved accuracy ................ errod =',errod
427 ! write(17,*) 'Number of iterations ............. iter =',iter
428 ! write(17,*) 'Magnetic axis coordinates ........ rm =', rm
429 ! write(17,*) ' zm =', zm
430 ! write(17,*) 'Plasma current ................... tokp =', tokp
431 ! write(17,*) ' cnor =', cnor
432 ! write(17,*) 'Magnetic axis PSI value .......... psim =', psim
433 ! write(17,*) 'Poloidal beta value ............. betpol =', betpol
434 ! write(17,*) 'Total beta value ............. bettot =', bettot
435 ! write(17,*) 'Plas.boun.poloidal current value... fvac =', fvac
436 ! write(17,*) 'Magnetic axis q value ............. qax =', qax
437 ! write(17,*) '*******************************************'
438 
439 
440  call wrb
441  !call out_b
442  !call retab
443  !call retab_f
444  !if(ngav.eq.0) then
445  !call retab_p
446  !endif
447 
448  !call cpu_time(time_end)
449  !write(*,*) 'operation time=',time_end-time_beg
450  !write(*,*) 'eqb done'
451 
452  return
453  end
454 
455 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
456 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
457 
458  subroutine retab
459 
460  include 'double.inc'
461  parameter(nursp=4000,nursp4=nursp+4,nursp6=nursp4*6)
462  include 'dim.inc'
463  include 'compol.inc'
464 
465  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
466  common/comurs/psit(nursp),purs(nursp),furs(nursp),wurs(nursp),
467  * nurs
468  dimension pstab(nursp),pptab(nursp),fptab(nursp)
469  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
470  real*8 cwk(4)
471  common
472  * /c_kpr/kpr
473 
474  write(*,*) 'retab:recalculation p and ff'
475 
476  do i=1,iplas
477  pstab(i)=1.d0-psia(i)
478  enddo
479 
480  call e01baf(iplas,pstab,dfdpsi,rrk,cck,
481  * iplas+4,wrk,6*iplas+16,ifail)
482 
483  if(ifail.ne.0.and.kpr.eq.1) write(*,*) 'ifail=',ifail
484 
485  do i=2,nurs-1
486 
487  zpsi=1.d0-psit(i)
488 
489  call e02bcf(iplas+4,rrk,cck,zpsi,0,cwk,ifail)
490 
491  if(ifail.ne.0.and.kpr.eq.1) write(*,*) 'ifail=',ifail
492 
493  furs(i)=cwk(1)
494 
495  enddo
496 
497  furs(1)=dfdpsi(iplas)
498  furs(nurs)=dfdpsi(1)
499 
500  call e01baf(iplas,pstab,dpdpsi,rrk,cck,
501  * iplas+4,wrk,6*iplas+16,ifail)
502 
503  if(ifail.ne.0.and.kpr.eq.1) write(*,*) 'ifail=',ifail
504 
505  do i=2,nurs-1
506 
507  zpsi=1.d0-psit(i)
508 
509  call e02bcf(iplas+4,rrk,cck,zpsi,0,cwk,ifail)
510 
511  if(ifail.ne.0.and.kpr.eq.1) write(*,*) 'ifail=',ifail
512 
513  purs(i)=cwk(1)
514 
515  enddo
516 
517  purs(1)=dpdpsi(iplas)
518  purs(nurs)=dpdpsi(1)
519 
520  ! do i=1,nurs
521  ! write(6,*) 'p/f,ps:',purs(i)/(furs(i)+1.d-15),psit(i),i
522  ! write(6,*) 'p,f:',purs(i),furs(i),i
523  ! pause ' '
524  ! enddo
525 
526  dpsi=1.d0/(nurs-1.d0)
527  ppp(1)=0.d0
528  fff(1)=0.d0
529 
530  do 20 i=2,nurs
531  ppp(i)=ppp(i-1) +(purs(i-1)+purs(i))*dpsi*0.5d0
532  fff(i)=fff(i-1) +(furs(i-1)+furs(i))*dpsi*0.5d0
533  20 continue
534 
535  return
536  end
537 
538 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
539 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
540 
541  subroutine retab_l
542 
543  include 'double.inc'
544  include 'dim.inc'
545  include 'compol.inc'
546  include 'urs.inc'
547 
548  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
549  dimension pstab(nursp),pptab(nursp),fptab(nursp)
550  common
551  * /c_kpr/kpr
552 
553  !write(*,*) 'retab_L:recalculation p and ff'
554 
555  do i=1,iplas
556  pstab(i)=1.d0-psia(i)
557  enddo
558 
559  do it=2,nurs-1
560  zpsi=1.d0-psit(it)
561 
562  do i=1,iplas-1
563  if(zpsi.gt.pstab(i) .AND. zpsi.le.pstab(i+1)) then
564  ic=i
565  exit
566  endif
567  enddo
568 
569  dpsi=pstab(ic+1)-pstab(ic)
570  furs(it)=( dfdpsi(ic)*(pstab(ic+1)-zpsi)+
571  * dfdpsi(ic+1)*(zpsi-pstab(ic)) )/dpsi
572 
573  purs(it)=( dpdpsi(ic)*(pstab(ic+1)-zpsi)+
574  * dpdpsi(ic+1)*(zpsi-pstab(ic)) )/dpsi
575  enddo
576 
577  furs(1)=dfdpsi(iplas)
578  furs(nurs)=dfdpsi(1)
579 
580  purs(1)=dpdpsi(iplas)
581  purs(nurs)=dpdpsi(1)
582 
583  dpsi=1.d0/(nurs-1.d0)
584  ppp(1)=0.d0
585  fff(1)=0.d0
586 
587  do 20 i=2,nurs
588  ppp(i)=ppp(i-1) +(purs(i-1)+purs(i))*dpsi*0.5d0
589  fff(i)=fff(i-1) +(furs(i-1)+furs(i))*dpsi*0.5d0
590  20 continue
591 
592  return
593  end
594 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
595 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
596 
597  subroutine retab_l_ast(pstab,pptab,fptab,ntab)
598 
599  include 'double.inc'
600  !include 'dim.inc'
601  !include 'compol.inc'
602  include 'urs.inc'
603 
604  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
605  dimension pstab(*),pptab(*),fptab(*)
606  common
607  * /c_kpr/kpr
608 
609  write(*,*) 'retab_L:recalculation p and ff'
610 
611 
612  do it=2,nurs-1
613  zpsi=1.d0-psit(it)
614 
615  do i=1,ntab-1
616  if(zpsi.gt.pstab(i) .AND. zpsi.le.pstab(i+1)) then
617  ic=i
618  exit
619  endif
620  enddo
621 
622  dpsi=pstab(ic+1)-pstab(ic)
623  furs(it)=( fptab(ic)*(pstab(ic+1)-zpsi)+
624  * fptab(ic+1)*(zpsi-pstab(ic)) )/dpsi
625 
626  purs(it)=( pptab(ic)*(pstab(ic+1)-zpsi)+
627  * pptab(ic+1)*(zpsi-pstab(ic)) )/dpsi
628  enddo
629 
630  furs(1)=fptab(ntab)
631  furs(nurs)=fptab(1)
632 
633  purs(1)=pptab(ntab)
634  purs(nurs)=pptab(1)
635 
636  dpsi=1.d0/(nurs-1.d0)
637  ppp(1)=0.d0
638  fff(1)=0.d0
639 
640  do 20 i=2,nurs
641  ppp(i)=ppp(i-1) +(purs(i-1)+purs(i))*dpsi*0.5d0
642  fff(i)=fff(i-1) +(furs(i-1)+furs(i))*dpsi*0.5d0
643  20 continue
644 
645  return
646  end
647 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
648 
649  subroutine retab_f
650 
651  include 'double.inc'
652  include 'dim.inc'
653  include 'compol.inc'
654  include 'urs.inc'
655 
656  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
657  dimension pstab(nursp),pptab(nursp),fptab(nursp)
658  common
659  * /c_kpr/kpr
660 
661  write(*,*) 'retab_F:linear recalculation ff table'
662 
663  do i=1,iplas
664  pstab(i)=1.d0-psia(i)
665  enddo
666 
667  do it=2,nurs-1
668  zpsi=1.d0-psit(it)
669 
670  do i=1,iplas-1
671  if(zpsi.gt.pstab(i) .AND. zpsi.le.pstab(i+1)) then
672  ic=i
673  exit
674  endif
675  enddo
676 
677  dpsi=pstab(ic+1)-pstab(ic)
678  furs(it)=( dfdpsi(ic)*(pstab(ic+1)-zpsi)+
679  * dfdpsi(ic+1)*(zpsi-pstab(ic)) )/dpsi
680 
681 ! purs(it)=( dpdpsi(ic)*(pstab(ic+1)-zpsi)+
682 ! * dpdpsi(ic+1)*(zpsi-pstab(ic)) )/dpsi
683  enddo
684 
685  furs(1)=dfdpsi(iplas)
686  furs(nurs)=dfdpsi(1)
687 
688 ! purs(1)=dpdpsi(iplas)
689 ! purs(nurs)=dpdpsi(1)
690 
691  dpsi=1.d0/(nurs-1.d0)
692 ! ppp(1)=0.d0
693  fff(1)=0.d0
694 
695  do 20 i=2,nurs
696 ! ppp(i)=ppp(i-1) +(purs(i-1)+purs(i))*dpsi*0.5d0
697  fff(i)=fff(i-1) +(furs(i-1)+furs(i))*dpsi*0.5d0
698  20 continue
699 
700  return
701  end
702 
703 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
704 
705  subroutine retab_p
706 
707  include 'double.inc'
708  include 'dim.inc'
709  include 'compol.inc'
710  include 'urs.inc'
711 
712  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
713  dimension pstab(nursp),pptab(nursp),fptab(nursp)
714  common
715  * /c_kpr/kpr
716 
717  write(*,*) 'retab_p:linear recalculation p table'
718 
719  do i=1,iplas
720  pstab(i)=1.d0-psia(i)
721  enddo
722 
723  do it=2,nurs-1
724  zpsi=1.d0-psit(it)
725 
726  do i=1,iplas-1
727  if(zpsi.gt.pstab(i) .AND. zpsi.le.pstab(i+1)) then
728  ic=i
729  exit
730  endif
731  enddo
732 
733  dpsi=pstab(ic+1)-pstab(ic)
734 ! furs(it)=( dfdpsi(ic)*(pstab(ic+1)-zpsi)+
735 ! * dfdpsi(ic+1)*(zpsi-pstab(ic)) )/dpsi
736 
737  purs(it)=( dpdpsi(ic)*(pstab(ic+1)-zpsi)+
738  * dpdpsi(ic+1)*(zpsi-pstab(ic)) )/dpsi
739  enddo
740 
741 ! furs(1)=dfdpsi(iplas)
742 ! furs(nurs)=dfdpsi(1)
743 
744  purs(1)=dpdpsi(iplas)
745  purs(nurs)=dpdpsi(1)
746 
747  dpsi=1.d0/(nurs-1.d0)
748  ppp(1)=0.d0
749 ! fff(1)=0.d0
750 
751  do 20 i=2,nurs
752  ppp(i)=ppp(i-1) +(purs(i-1)+purs(i))*dpsi*0.5d0
753 ! fff(i)=fff(i-1) +(furs(i-1)+furs(i))*dpsi*0.5d0
754  20 continue
755 
756  return
757  end
758 
759 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
760 
761  subroutine bt_tot(bettot)
762 
763  include 'double.inc'
764  include 'dim.inc'
765  include 'compol.inc'
766 
767  volcen=0.d0
768 
769  do j=2,nt1
770  volcen=volcen+vol(1,j)*0.5d0
771  enddo
772 
773  volpl=volcen
774 
775  psn=psin(1,2)
776 
777  zpres=funppp(psn)
778  pintg=zpres*volcen
779 
780  do i=2,iplas1
781  do j=2,nt1
782 
783  psn=psin(i,j)
784  zpres=funppp(psn)
785  volk=vol1(i,j)+vol2(i-1,j)+vol3(i-1,j-1)+vol4(i,j-1)
786  pintg=pintg+zpres*volk
787  volpl=volpl+volk
788 
789  enddo
790  enddo
791 
792  paverg=pintg/volpl
793 
794  bettot=2.d0*paverg/(b0ax*b0ax)
795  bettot=bettot*(psim-psip)*cnor
796 
797  return
798  end
799 
800 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
801 
802 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
803 
804  subroutine presol(i_betp,betplx,betpol)
805 
806  include 'double.inc'
807  include 'dim.inc'
808  include 'compol.inc'
809  common
810  * /c_kpr/kpr
811 
812  imov=0
813  iter=0
814  itin=0
815 
816 ch4astra call metric
817  call metrix
818  call matcof
819  call matpla
820 
821 c write(6,*) 'matpla'
822 
823  1000 continue
824 
825  iter = iter+1
826  itin = itin+1
827  if(kpr.eq.1) then
828  write(*,*)'iter=',iter,itin
829  endif
830  call rightg
831 
832 c write(6,*) 'rightg'
833 
834  if(i_betp.eq.1) then
835 
836  if(iter.gt.4) call skbetp(betplx,betpol)
837 
838  endif
839 
840  ! call bt_pol(betpol)
841  ! write(6,*)'betpol',betpol,'*********'
842 
843  call solint(imov)
844 
845  ! write(6,*)'solve(psi)'
846  ! call wrb
847  ! pause 'wrd after solve'
848 
849  call spider_remesh(erro,errpsi,imov)
850 
851  if(kpr.eq.1) then
852  write(*,*)'presol: errpsi',errpsi
853  endif
854 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
855 
856  ! if(ngav.eq.1) then
857  ! call wrb
858  ! pause 'wrb after regrid '
859  ! endif
860 
861 !!!!!!!!!!!! accurasy test parameters !!!!!!!!!!!!!!!
862 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
863 
864  errod = 0.5d0*erro/dabs(z(iplas,2)-zm)
865 
866  if(errpsi.lt.1.d-4 .or. iter.ge.50) go to 2000
867 
868  go to 1000
869 
870  2000 continue
871 
872  psimax=psi(1,2)
873  imax=1
874  jmax=2
875 
876  do i=2,iplas1
877  do j=2,nt1
878 
879  if(psi(i,j).gt.psimax) then
880  psimax=psi(i,j)
881  imax=i
882  jmax=j
883  endif
884 
885  enddo
886  enddo
887 
888  if(kpr.eq.1) then
889  write(*,*) 'presol: '
890  write(*,*) 'number of iterations ',iter
891  write(*,*) 'accuracy ',errpsi
892  endif
893 c--------------------------------------------------------------------
894 
895  call prgrid(imax,jmax,erro,1,errpsi)
896 
897  platok = tokp
898  rax = rm
899  zax = zm
900  psax = psim
901 
902  if(kpr.eq.1) then
903  write(*,*) 'rax,zax',rax,zax
904  write(*,*) 'plas.current',platok,cnor
905  write(*,*) 'psax',psax
906  endif
907 
908  psax = psim
909  psipla = psim-psip
910  fvac = f(iplas)
911 
912  call bt_pol(betpol)
913  call wrb
914  ! call out_b
915 
916  return
917  end
918 
919 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
920 
921 
922 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
923 
924  subroutine p_pert
925 
926  include 'double.inc'
927  include 'dim.inc'
928  include 'compol.inc'
929 
930  ampl=1.0d0 * dpdpsi(1)
931  !ampl=0.0d0
932  xa0=0.9999d0
933  xw=0.02d0
934 
935  do i=1,iplas
936  xa=dfloat(i-1)/dfloat(iplas-1)
937  fun=( 1.d0-dtanh(((xa0-xa)/xw)**2) )
938  dpdpsi(i)=dpdpsi(i)+ampl*fun
939  enddo
940 
941  return
942  end
943 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
944  subroutine put_tim(dt,time)
945 
946  include 'double.inc'
947  include 'dim.inc'
948  include 'compol.inc'
949 
950  dtim=dt
951  ctim=time
952  return
953  end
954 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
955  subroutine get_tim(dt,time)
956 
957  include 'double.inc'
958  include 'dim.inc'
959  include 'compol.inc'
960 
961  dt=dtim
962  time=ctim
963  return
964  end
965 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
966 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
967 
968  subroutine get_www(roomega)
969 
970  include 'double.inc'
971  include 'dim.inc'
972  include 'compol.inc'
973  include 'urs.inc'
974 
975  common/comppp/ ppp(nursp),fff(nursp),www(nursp)
976 
977  roomega=www(nurs)*psim*cnor
978  write(*,*) 'roomega',roomega
979 
980 
981  return
982  end
983 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
984 
985 
986 
retab_p
subroutine retab_p
Definition: B_eqb.f:705
e02bcf
subroutine e02bcf(NCAP7, K, C, X, LEFT, S, IFAIL)
Definition: NAG.f:2761
psib_ext
subroutine psib_ext(psex_av)
Definition: B_rig.f:614
put_tim
subroutine put_tim(dt, time)
Definition: B_eqb.f:944
fun
subroutine fun(X, F)
Definition: Ev2.f:10
spider_remesh
subroutine spider_remesh(erro, errpsi, imov)
Definition: B_grid.f:5
metrix
subroutine metrix
Definition: B_metric_O.f:16
retab
subroutine retab
Definition: B_eqb.f:458
solint
subroutine solint(imov)
Definition: B_sol.f:2
grid_b
subroutine grid_b(igdf, nstep)
Definition: B_grid.f:2467
eqb
subroutine eqb(alf0, alf1, alf2, bet0, bet1, bet2, alw0, alw1, alw2,
Definition: B_eqb.f:3
helena21_mod::dpdpsi
real(r8) function dpdpsi(psi_n)
Definition: helena21_mod.f90:122
grid_p1
subroutine grid_p1(igdf, nstep)
Definition: B_grid.f:1470
retab_f
subroutine retab_f
Definition: B_eqb.f:649
rightg
subroutine rightg
Definition: B_rig.f:10
get_tim
subroutine get_tim(dt, time)
Definition: B_eqb.f:955
grdef
subroutine grdef(igdf)
Definition: com_sub.f:7
e01baf
subroutine e01baf(M, X, Y, K, C, LCK, WRK, LWRK, IFAIL)
Definition: NAG.f:2511
retab_l_ast
subroutine retab_l_ast(pstab, pptab, fptab, ntab)
Definition: B_eqb.f:597
grid_b1
subroutine grid_b1(igdf, nstep)
Definition: B_grid.f:2704
skbetp
subroutine skbetp(betplx, betpol)
Definition: com_sub.f:781
matpla
subroutine matpla
Definition: com_sub.f:1061
matcof
subroutine matcof
Definition: B_metric_O.f:432
bt_pol
subroutine bt_pol(betpol)
Definition: com_sub.f:747
wrb
subroutine wrb
Definition: B_wrd.f:1
bongri
subroutine bongri
Definition: com_sub.f:242
qst_b
subroutine qst_b
Definition: com_sub.f:576
bt_tot
subroutine bt_tot(bettot)
Definition: B_eqb.f:761
prgrid
subroutine prgrid(imax, jmax, erro, imov, errpsi)
Definition: B_grid.f:60
get_www
subroutine get_www(roomega)
Definition: B_eqb.f:968
ppp
subroutine ppp
Definition: ppplib.f:1
retab_l
subroutine retab_l
Definition: B_eqb.f:541
p_pert
subroutine p_pert
Definition: B_eqb.f:924
grid_p0
subroutine grid_p0(igdf, nstep)
Definition: B_grid.f:1457
presol
subroutine presol(i_betp, betplx, betpol)
Definition: B_eqb.f:804
grid_1
subroutine grid_1(igdf, nstep)
Definition: B_grid.f:1112
grid_0
subroutine grid_0(igdf, nstep)
Definition: B_grid.f:1484