ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
interfac.f
Go to the documentation of this file.
1  subroutine astra2spider(neql,nteta,nbnd,rzbnd,key_dmf,
2  * na1,yeqpf,yeqff,fp,ipl,
3  * rtor,btor,rho,roc,nstep,yreler,mu,
4  * cc,te,cubs,cd,key_ini,eqdfn,pres,cu,
5  * key_0stp,key_pres )
6 
7 
8  use durs_d_modul
9  use ppf_modul
10  use bnd_modul
11 
12  include 'double.inc'
13  include 'dimpl1.inc'
14  parameter(nbtabp=1000)
15  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
16 ! common /com_bas/ rbtab(nbtabp),zbtab(nbtabp),nbtab
17 ! common/com_pas/ pstab(nursp),pptab(nursp),fptab(nursp),nutab
18  common /creler/ relerr
19  common /com_0st/ key_0st,key_prs
20 
21  real*8 rout(nrp,ntp),zout(nrp,ntp)
22 
23  real t_start, t_finish
24 
25  real*8 rzbnd(*),yeqpf(*),yeqff(*),fp(*),ipl,ybetpl,yli3,
26  * rtor,btor,rho(*),roc,yreler,mu(*),
27  * cc(*),te(*),cubs(*),cd(*),pres(*),cu(*)
28 
29  character*40 eqdfn
30 
31  real(8), allocatable :: eqpf(:), eqff(:)
32 
33  allocate( eqpf(na1), eqff(na1) )
34 
35 
36  do i=1,na1
37  eqpf(i)=yeqpf(i)
38  eqff(i)=yeqff(i)
39  enddo
40 
41  pi=3.14159265359d0
42  amu0=0.4d0*pi
43 
44  relerr =yreler
45  !write(*,*) 'b_eqb, rz', (rzbnd(j),j=1,nbnd)
46  !write(*,*) 'nstep', nstep
47  if(nbnd.gt.nbtabp) then
48  write(*,*) 'spider:nbnd exeeds maximum value',nbtabp
49  write(*,*) 'nbnd=',nbnd
50  write(*,*) 'program is interrupted'
51  stop
52  endif
53 
54  if(na1.gt.nursp) then
55  write(*,*) 'spider:na1 exeeds maximum value',nursp
56  write(*,*) 'na1=',na1
57  write(*,*) 'program is interrupted'
58  stop
59  endif
60 
61  if( neql.gt.nrp .or. nteta.gt.ntp) then
62  write(*,*) 'spider: neql or nteta exeeds maximum value'
63  write(*,*) 'neql=',neql,'nrp=',nrp
64  write(*,*) 'nteta=',nteta,'ntp=',ntp
65  write(*,*) 'program is interrupted'
66  stop
67  endif
68 
69  key_0st=key_0stp
70  key_prs=key_pres
71  nurs=-399
72  i_betp=0
73  i_bsh=1
74  igdf=2
75  i_eqdsk=0
76  epsro=1.d-6
77  maxit=250
78  ! nstep=0
79 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
80  n_tht=nteta
81  n_psi=neql
82  platok=ipl
83  tokf=ipl
84 
85  call put_ipl(tokf)
86 
87  if(nstep.eq.0 .AnD. key_ini.eq.0) then
88  call tab_efit(tokf,psax,eqdfn,rax,zax,b0,r0)
89  keyctr=0
90  key_0st=0
91 ! write(fname,'(a,a)') path(1:kname),'durs_d.dat'
92  if(i_betp.eq.0) then
93  betplx = 0.d0
94  endif
95  if(keyctr.ne.2) then
96  psax = 1.d0
97  endif
98  if(nurs.lt.0) then
99  alf0 = 0.d0
100  alf1 = 0.d0
101  alf2 = 0.d0
102  bet0 = 0.d0
103  bet1 = 0.d0
104  bet2 = 0.d0
105  endif
106 ! open(1,file=fname,form='formatted')
107 ! !open(1,file='durs_d.dat')
108 ! write(1,*) n_tht,n_psi,igdf,nurs,keyctr,i_eqdsk,i_betp
109 ! write(1,*) epsro,betplx,tokf,psax,b0,r0,rax,zax
110 ! write(1,*) alf0,alf1,alf2,bet0,bet1,bet2
111 ! close(1)
112  return
113  endif
114 
115  nutab=na1+1
116 
117  if(allocated(pstab))
118  %deallocate( pstab, pptab, fptab )
119  allocate( pstab(na1+1), pptab(na1+1), fptab(na1+1) )
120 
121  nbtab=nbnd
122 
123  if(allocated(rbtab))
124  %deallocate( rbtab, zbtab )
125  allocate( rbtab(nbtab), zbtab(nbtab) )
126 
127  do ib=1,nbtab
128 
129  rbtab(ib)=rzbnd(ib)
130  zbtab(ib)=rzbnd(ib+nbnd)
131 
132  enddo
133 
134 
135 ! write(fname,'(a,a)') path(1:kname),'tab_bnd.dat'
136 ! open(1,file=fname,form='formatted')
137 ! !open(1,file='tab_bnd.dat')
138 ! write(1,*) nbtab
139 ! do ib=1,nbtab
140 ! write(1,*) rbtab(ib),zbtab(ib)
141 ! enddo
142 ! close(1)
143 
144 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
145 !__CH
146 
147  rbmax=rbtab(1)
148  rbmin=rbtab(1)
149  zbmax=zbtab(1)
150  zbmin=zbtab(1)
151 
152  do ib=1,nbtab
153 
154  if(rbtab(ib).ge.rbmax) rbmax=rbtab(ib)
155  if(rbtab(ib).le.rbmin) rbmin=rbtab(ib)
156  if(zbtab(ib).ge.zbmax) zbmax=zbtab(ib)
157  if(zbtab(ib).le.zbmin) zbmin=zbtab(ib)
158 
159  enddo
160 
161  rc0=0.5d0*(rbmax+rbmin)
162  zc0=0.5d0*(zbmax+zbmin)
163 
164  rax=rc0
165  zax=zc0
166 !__CH
167  b0=btor
168  r0=rtor
169 !__CH
170 
171 
172 !! p' and ff' profiles
173 
174  if(nstep.eq.0) then
175 
176  if(key_0st.eq.1) then
177  do i=1,na1
178  eqpf(i)=0.d0
179  eqff(i)=cu(i)
180  enddo
181  tokf=1.d0
182 
183  do i=1,na1
184  eqpf(i)=eqpf(i)*amu0
185  eqff(i)=eqff(i)*amu0
186  enddo
187 
188  !nutab=na1+1
189  pstab(1)=0.d0
190  do i=2,nutab
191  pstab(i)=(rho(i-1)/rho(na1))**2
192  enddo
193 
194  call profro_cu(na1,neql,rtor,btor,cu,rho,roc)
195  call proffi_pres(na1,neql,pres,rho,roc)
196 
197 
198  elseif(key_0st.eq.0) then
199 
200  do i=1,na1
201  eqpf(i)=eqpf(i)*amu0
202  eqff(i)=eqff(i)*amu0
203  enddo
204 !! exstrapolation psi on magn. axis
205 
206  rh0=0.d0
207 
208  rh1=rho(1)
209  rh2=rho(2)
210  rh3=rho(3)
211 
212  ps1=fp(1)
213  ps2=fp(2)
214  ps3=fp(3)
215 
216 
217  !call EXTRP2(rh0,ps0, rh1,rh2,rh3, ps1,ps2,ps3)
218 
219  ps0=(ps1*rh2**2-ps2*rh1**2)/(rh2**2-rh1**2)
220 
221 
222  psb=fp(na1)
223 
224  pstab(1)=0.d0
225 
226 
227  !nutab=na1+1
228 
229  do i=2,nutab
230  pstab(i)=(fp(i-1)-ps0)/(psb-ps0)
231  if(pstab(i).le.pstab(i-1)) then
232  write(*,*) 'psi is nonmonotonic!!! ', pstab(i),i
233  endif
234  enddo
235 
236 
237  endif !(key_0st
238 
239 !! exstrapolation on magn. axis
240 
241  rh0=0.d0
242 
243  rh1=rho(1)
244  rh2=rho(2)
245  rh3=rho(3)
246 
247  ppx1=eqpf(1)
248  ppx2=eqpf(2)
249  ppx3=eqpf(3)
250 
251  ffx1=eqff(1)
252  ffx2=eqff(2)
253  ffx3=eqff(3)
254 
255  call extrp2(rh0,ppx0, rh1,rh2,rh3, ppx1,ppx2,ppx3)
256  call extrp2(rh0,ffx0, rh1,rh2,rh3, ffx1,ffx2,ffx3)
257 
258  !pptab(1)=ppx0
259  !fptab(1)=ffx0
260  pptab(1)=ppx1
261  fptab(1)=ffx1
262 
263  !nutab=na1+1
264 
265  do i=2,nutab
266  pptab(i)=eqpf(i-1)
267  fptab(i)=eqff(i-1)
268  enddo
269 
270 !normalization
271 
272  do i=1,nutab
273  pptab(i)=pptab(i)/rtor
274  fptab(i)=fptab(i)*rtor
275  enddo
276 
277 ! write(fname,'(a,a)') path(1:kname),'tabppf.dat'
278 ! open(1,file=fname,form='formatted')
279 ! !open(1,file='tabppf.dat')
280 
281 ! write(1,*) nutab
282 
283 ! do i=1,nutab
284 ! write(1,*) pstab(i),pptab(i),fptab(i)
285 ! enddo
286 
287 ! close(1)
288 
289  endif ! (nstep.eq.0
290 
291 !__CH
292 
293  if(nstep.ne.0) then
294 
295  !call profro_pres_L(na1,neql,pres,rho,roc)
296  !call profro_pres(na1,neql,pres,rho,roc)
297  !call proffi_pres(na1,neql,pres,rho,roc)
298 
299  if(key_dmf.ne.0) then
300 
301  call profro_cu(na1,neql,rtor,btor,cu,rho,roc)
302 
303  if(key_pres.eq.0) then
304  call profro_p(na1,neql,rtor,eqpf,rho,roc)
305  else
306  call proffi_pres(na1,neql,pres,rho,roc)
307  endif
308 
309  if(key_dmf.eq.-2 .OR. key_dmf.eq.-3) then
310  call profro_sbc(na1,neql,rtor,btor,cc,te,cubs,cd,rho,roc)
311  endif
312 
313  else !(key_dmf=0)
314 
315  !call retab_L_ast(pstab,pptab,fptab,nutab) !
316  call profro(na1,neql,rtor,eqpf,eqff,rho,roc)
317 
318  endif
319 
320  endif
321 !__CH
322  !do i=1,nutab
323  ! pptab(i)=eqpf(i)
324  ! fptab(i)=eqff(i)
325  ! pstab(i)=fp(i)
326  !enddo
327 
328  !open(17,file='out.pr')
329 
330 
331 
332  1149 format(a40)
333 
334  if(nstep.eq.0) then
335  keyctr=0
336 ! write(fname,'(a,a)') path(1:kname),'durs_d.dat'
337  if(i_betp.eq.0) then
338  betplx = 0.d0
339  endif
340  if(keyctr.ne.2) then
341  psax = 1.d0
342  endif
343  if(nurs.lt.0) then
344  alf0 = 0.d0
345  alf1 = 0.d0
346  alf2 = 0.d0
347  bet0 = 0.d0
348  bet1 = 0.d0
349  bet2 = 0.d0
350  endif
351 ! open(1,file=fname,form='formatted')
352 ! !open(1,file='durs_d.dat')
353 ! write(1,*) n_tht,n_psi,igdf,nurs,keyctr,i_eqdsk,i_betp
354 ! write(1,*) epsro,betplx,tokf,psax,b0,r0,rax,zax
355 ! write(1,*) alf0,alf1,alf2,bet0,bet1,bet2
356 ! close(1)
357  endif
358 
359 
360  do i=1,na1
361  eqpf(i)=eqpf(i)/amu0
362  eqff(i)=eqff(i)/amu0
363  enddo
364 
365  return
366  end
367 
368 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
369  subroutine spider2astra(rout,zout,rtor,btor,rho,roc,na1,
370  * g11,g22,g33,vr,vrs,slat,gradro,rocnew,
371  * mu,ipol,bmaxt,bmint,bdb02,b0db2,bdb0,droda,
372  * yreler,yli3,ni_p,nj_p,platok,cu,fp,pres,w_dj,
373  * yfofb)
374 
375 
376  use ppf_modul
377  use bnd_modul
378 
379  include 'double.inc'
380  include 'dimpl1.inc'
381  parameter(nbtabp=1000)
382  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
383 ! common /com_bas/ rbtab(nbtabp),zbtab(nbtabp),nbtab
384 ! common/com_pas/ pstab(nursp),pptab(nursp),fptab(nursp),nutab
385  common /creler/ relerr
386  integer na1,ni_p,nj_p
387  real*8 rout(ni_p,nj_p),zout(ni_p,nj_p)
388 
389  real*8
390  * g11(na1),g22(na1),g33(na1),
391  * vr(na1),vrs(na1),slat(na1),gradro(na1),
392  * rocnew,ybetpl,betpol,platok,
393  * mu(na1),ipol(na1),bmaxt(na1),bmint(na1),
394  * bdb02(na1),bdb0(na1),b0db2(na1),droda(na1),
395  * rtor,btor,rho(na1),roc,cu(na1),fp(na1),pres(na1),
396  * w_dj(na1),yfofb(na1),traps(na1)
397 
398  character*40 eqdfn
399 
400  betpol=0.0 !!DPC!!
401 
402  call get_rz(rout,zout,ni_p,nj_p)
403 
404  do i=1,na1
405  bmint(i) =0.0d0
406  enddo
407 cw write(*,*) 'after eqb'
408 
409  ybetpl=betpol
410 
411  call comet(na1,platok,
412  * rtor,btor,rho,roc,nstep,
413  * g11,g22,g33,vr,vrs,slat,gradro,rocnew,
414  * mu,ipol,bmaxt,bmint,bdb02,b0db2,bdb0,droda,
415  * ybetpl,yli3,cu,fp,pres,w_dj,yfofb,traps)
416 
417 
418 ! deallocate( pstab, pptab, fptab )
419 ! deallocate( rbtab, zbtab )
420 
421  return
422  end
423 
424 
425 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!last change: 27.05.03
426 
427  subroutine comet(na1,platok,
428  * rtor,btor,rho,roc,nstep,
429  * g11,g22,g33,vr,vrs,slat,gradro,rocnew,
430  * mu,ipol,bmaxt,bmint,bdb02,b0db2,bdb0,droda,
431  * ybetpl,yli3,cu,fp,pres,w_dj,yfofb,traps)
432 
433  include 'double.inc'
434  include 'dim.inc'
435  parameter(nrpl=nrp+1)
436  parameter(nrpl4=nrpl+4,nrpl6=nrpl4*6)
437  include 'compol.inc'
438  common /combsh/ rm0,zm0,rc0,zc0,asp0,el_up,el_lw,tr_up,tr_lw,nbsh
439  common /com_jb/ bj_av(nrp),curfi_av(nrp)
440  common/com_heat_dj/ wdj(nrp)
441  common /com_trap/ trap(nrp)
442 
443  real*8 rhos(nrp),vols(nrpl),g11s(nrpl),g22s(nrpl),g33s(nrpl)
444  real*8 gradrs(nrpl),amus(nrpl),fnors(nrpl),drodas(nrpl)
445  real*8 dvdro(nrpl),sa(nrp),rhocs(nrpl),rhocn(nrpl),rhosn(nrpl)
446  real*8 b_maxt(nrpl),b_mint(nrpl),b_db02(nrpl),b_db0(nrpl),
447  * b_0db2(nrpl),d_roda(nrpl)
448  real*8 psi_1(nrp)
449  real*8 prs_1(nrp)
450 
451  real*8 rrk(nrpl4),cck(nrpl4),wrk(nrpl6)
452  real*8 cwk(4)
453  !real*8 rhowr(1000)
454  real(8), allocatable :: rhowr(:), rhowrh(:)
455 
456  real*8 rtor,btor,rho(na1),roc,ybetpl,yli3,
457  * g11(na1),g22(na1),g33(na1),
458  * vr(na1),vrs(na1),slat(na1),gradro(na1),rocnew,droda(na1),
459  * mu(na1),ipol(na1),bmaxt(na1),bmint(na1),
460  * bdb02(na1),b0db2(na1),bdb0(na1),
461  * cu(na1),fp(na1),pres(na1),w_dj(na1)
462 
463  real*8 yfofb(na1)
464 
465  real*8 traps(nrpl)
466  dimension btot(nrp,ntp)
467 
468  !common /fp_sav/ fp_0(1000)
469  !common /fp_dot/ dfpdt(1000),nna1
470 
471  sqrt(xx)=dsqrt(xx)
472 
473  !do i=1,na1
474  !fp_0(i)=fp(i)
475  !enddo
476  !nna1=na1
477 
478  allocate( rhowr(na1), rhowrh(na1) )
479 
480  na=na1-1
481  platok=tokp
482  do i=1,iplas
483  !rhos(i)=sqrt(flucf(i)/(pi*btor))
484  rhos(i)=sqrt(flx_fi(i)/(pi*btor))
485  enddo
486 
487  call get_psibon(psi_bn1)
488 
489  do i=1,iplas
490  psi_1(i)=(psi(i,2)+psi_bn1)*2.d0*pi
491  enddo
492 
493  do i=1,iplas
494  psn=psia(i)
495  prs_1(i)=funppp(psn)*(psim-psip)*cnor/amu0
496  enddo
497 cp
498 cp
499  rhocs(1)=0.d0
500  rhos(1)=0.d0
501 
502  do i=2,iplas
503  rhocs(i)=0.5d0*(rhos(i-1)+rhos(i))
504  enddo
505 
506  rocnew=rhos(iplas)
507  rhocs(iplas+1)=rhos(iplas)
508 
509 c=============================andrey=02.03.2003vvvvvvvvvvvvv
510  arc0=rhocs(iplas+1)
511  arc1=rhocs(iplas)
512  arc2=rhocs(iplas-1)
513  arc3=rhocs(iplas-2)
514  spw1=(arc0-arc2)*(arc0-arc3)
515  & /(arc1-arc2)/(arc1-arc3)
516  spw2=(arc0-arc1)*(arc0-arc3)
517  & /(arc2-arc1)/(arc2-arc3)
518  spw3=(arc0-arc1)*(arc0-arc2)
519  & /(arc3-arc1)/(arc3-arc2)
520  q(iplas)=spw1*q(iplas-1)+spw2*q(iplas-2)+spw3*q(iplas-3)
521 ! trap(iplas)=SPW1*trap(iplas-1)+SPW2*trap(iplas-2)+
522 ! & SPW3*trap(iplas-3)
523 
524 c=============================andrey=02.03.2003^^^^^^^^^^^^^^
525 
526 
527 
528  do i=1,iplas+1
529  rhocn(i)=rhocs(i)/rhos(iplas)
530  enddo
531  rhocn(1)=0.d0
532 
533  do i=1,iplas
534  rhosn(i)=rhos(i)/rhos(iplas)
535  drodas(i)=rhosn(i)
536  enddo
537  rhosn(1)=0.d0
538 cp
539 
540  rhocn(iplas+1)=1.d0
541  rhosn(iplas)=1.d0
542 
543  do i=1,na
544  rhowr(i)=rho(i)/roc
545  enddo
546  rhowr(na1)=1.d0-1.d-9
547  !assuming main astra grid: rho(j)=h*(j-0.5) j=1,2,...,na1-1
548  ! intermediate astra grid: rho_h(j)=h*j
549  do i=1,na
550  rhowrh(i)=(rhowr(i)+rhowr(i+1))*0.5d0
551  enddo
552  rhowrh(na)=rhowrh(na-1)+(rhowr(2)-rhowr(1))
553  rhowrh(na1)=1.d0-1.d-9
554 
555  sa(1)=0.d0
556 cw write(*,*) (rhocn(j),j=1,iplas)
557 cw write(*,*) (rhosn(j),j=1,iplas)
558 
559  do i=1,iplas-1
560 
561  dsqi=0.d0
562  dvoli=0.d0
563  av_r2=0.d0
564  av_gr1=0.d0
565  av_gr2=0.d0
566  av_grr2=0.d0
567 
568  u1=rhos(i)
569  u2=rhos(i+1)
570  u3=rhos(i+1)
571  u4=rhos(i)
572 
573  do j=2,nt1
574 
575  dvoli=dvoli+vol(i,j)
576  dsqi=dsqi+sr(i+1,j)
577 
578  r1=r(i,j)
579  r2=r(i+1,j)
580  r3=r(i+1,j+1)
581  r4=r(i,j+1)
582  r0=(r1+r2+r3+r4)*0.25d0
583 
584  z1=z(i,j)
585  z2=z(i+1,j)
586  z3=z(i+1,j+1)
587  z4=z(i,j+1)
588 
589  sss=s(i,j)
590 
591  dadr= 0.5d0*((u1+u2)*(z2-z1)+(u2+u3)*(z3-z2)+
592  * (u3+u4)*(z4-z3)+(u4+u1)*(z1-z4))/sss
593 
594  dadz=-0.5d0*((u1+u2)*(r2-r1)+(u2+u3)*(r3-r2)+
595  * (u3+u4)*(r4-r3)+(u4+u1)*(r1-r4))/sss
596 
597  gr2=dadr**2+dadz**2
598 
599  av_gr1=av_gr1+vol(i,j)*sqrt(gr2)
600  av_gr2=av_gr2+vol(i,j)*gr2
601  av_grr2=av_grr2+vol(i,j)*gr2/r0**2
602  av_r2=av_r2+vol(i,j)/r0**2
603 
604  enddo
605  sa(i+1)=dsqi*2.d0*pi
606 c vols(i+1)=dvoli
607  gradrs(i+1)=av_gr1/dvoli
608  amus(i+1)=(2.d0*pi)/q(i)
609  fnors(i+1)=f(i)/(rtor*btor)
610 
611  dvdro(i+1)=2.d0*pi*dvoli/(rhos(i+1)-rhos(i))
612 cc g11s(i+1)=av_gr2/dvoli
613 cc g22s(i+1)=av_grr2/(dvoli*4.d0*pi*pi)*rtor
614  g11s(i+1)=av_gr2*dvdro(i+1)/dvoli
615  g22s(i+1)=av_grr2*dvdro(i+1)/(dvoli*4.d0*pi*pi)*rtor/fnors(i+1)
616  g33s(i+1)=av_r2*rtor**2/dvoli
617  enddo
618 
619 cccc magnetic field
620 
621 cw write(*,*) 'before b'
622 
623  bp2av=0.d0
624  bp2=0.d0
625 cc do i=2,iplas-1
626  do i=1,iplas-1
627 
628  dvoli=0.d0
629  bmx=0.d0
630  bmn=1.d9*btor
631  bavr=0.d0
632  bavr2=0.d0
633  bavrd2=0.d0
634 
635  do j=2,nt1
636 
637  dvoli=dvoli+vol(i,j)
638 
639  r1=r(i,j)
640  r2=r(i+1,j)
641  r3=r(i+1,j+1)
642  r4=r(i,j+1)
643  r0=(r1+r2+r3+r4)*0.25d0
644 
645  bm_pol=psim*(psia(i+1)-psia(i))/st(i,j)
646  bp_pol=psim*(psia(i+1)-psia(i))/st(i,j+1)
647 
648  if(i.ne.1) then
649  bpol2v=bm_pol**2*( vol1(i,j)/sin1(i,j)+vol2(i,j)/sin2(i,j))+
650  + bp_pol**2*( vol3(i,j)/sin3(i,j)+vol4(i,j)/sin4(i,j))
651  else
652  bpol2v=bm_pol**2*( vol2(i,j)/sin2(i,j))+
653  + bp_pol**2*( vol3(i,j)/sin3(i,j) )
654  endif
655 
656  bpol2=bpol2v /vol(i,j)
657  b_fi2=(f(i)/r0)**2
658 
659  b_tot2=bpol2+b_fi2
660  b_tot=dsqrt(b_tot2)
661 !! aai 06/02/10{{{{{{{{{{{
662  btot(i,j)=b_tot
663 !! }}}}}}}}}}}
664 
665  bmx=dmax1(b_tot,bmx)
666  bmn=dmin1(b_tot,bmn)
667 
668  bavr2=bavr2+b_tot2*vol(i,j)
669  bavr =bavr +b_tot *vol(i,j)
670  bavrd2=bavrd2+vol(i,j)/(b_tot2+1.d-8)
671  bp2av=bp2av+bpol2v
672  bp2=bp2+bpol2v
673  enddo
674  b_maxt(i+1)=bmx
675  b_mint(i+1)=bmn
676  b_db02(i+1)=bavr2/(dvoli*btor**2)
677  b_0db2(i+1)=(bavrd2*btor**2)/(dvoli)
678  b_db0(i+1)=bavr/(dvoli*btor)
679 
680 
681  enddo
682  b_maxt(1)=.5d0*(b_maxt(2)+b_mint(2))
683  b_mint(1)=b_maxt(1)
684  b_db02(1)=b_db02(2)
685  b_0db2(1)=b_0db2(2)
686  b_db0(1)=b_db0(2)
687 
688 
689 c yli3 =4.d0*pi*bp2av/(rtor*tokp*tokp)
690  yli3 = 4.d0*pi*bp2/(rtor*(.4d0*pi*tokp)**2)
691 
692 
693 
694 !! aai 06/02/10{{{{{{{{{{{
695  do i=1,iplas-1
696  bmax=0.d0
697  svol=0.d0
698  do j=2,nt-1
699  b_ij=btot(i,j)
700  bmax=dmax1(bmax,b_ij)
701  svol=svol+vol(i,j)
702  enddo
703 
704  avr_v=0.d0
705  do j=2,nt-1
706  bnor=btot(i,j)/bmax
707  avr_v=avr_v+
708  & (b0ax/btot(i,j))**2
709  & *( 1.d0-dsqrt(1.d0-bnor)*(1.d0+.5d0*bnor) )*vol(i,j)
710  enddo
711 
712  trap(i)=avr_v/svol
713  traps(i+1)=trap(i)
714 
715  enddo
716 !! }}}}}}}}}}}
717 
718 c*new
719 
720 c goto 9991
721  ip0=1
722  ip1=2
723  ip2=3
724  ip3=4
725 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
726 
727  x3=rhocs(ip3)
728  x2=rhocs(ip2)
729  x1=rhocs(ip1)
730  x0=rhocs(ip0)
731 !!!!!!!!!!!!g11(1)
732 call b_extrp(x0, x1,x2,x3, ip1,ip2,ip3,ip0,g11s,jerr)
733 !!!!!!!!!!!!g22(1)
734 call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,g22s,jerr)
735 
736 !!!!!!!!!!!!g33(1)
737  call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,g33s,jerr)
738 
739 !!!!!!!!!!!!fnors(1)
740  call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,fnors,jerr)
741 
742 !!!!!!!!!!!!gradrs(1)
743  call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,gradrs,jerr)
744 
745 !!!!!!!!!!!!amus(1)
746  call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,amus,jerr)
747 !!!!!!!!!!!!tpaps(1)
748  call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,traps,jerr)
749 
750 
751 ! x0=rhocs(1)
752 ! x1=rhocs(2)
753 ! x2=rhocs(3)
754 ! x3=rhocs(4)
755 
756 !
757  9991 continue
758 
759  ip0=iplas+1
760  ip1=iplas-2
761  ip2=iplas-1
762  ip3=iplas
763 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
764 
765  x1=rhocs(ip3)
766  x2=rhocs(ip2)
767  x3=rhocs(ip1)
768  x0=rhocs(ip0)
769 
770 
771 cw write(*,*) 'passed 1'
772 !!!!!!!!!!!gradrs(iplas+1)
773  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,gradrs,jerr)
774 
775 !!!!!!!!!!!amus(iplas+1)
776  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,amus,jerr)
777 
778 
779 !!!!!!!!!!!dvdro(iplas+1)
780  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,dvdro,jerr)
781 
782 
783 !!!!!!!!!!!!g11(iplas+1)
784  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,g11s,jerr)
785 
786 !!!!!!!!!!!!g22(iplas+1)
787  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,g22s,jerr)
788 
789 !!!!!!!!!!!!g33(iplas+1)
790  call b_extrp(x0,x1,x2,x3, ip3,ip2,ip1,ip0,g33s,jerr)
791 
792 !!!!!!!!!!!!b_maxt(iplas+1)
793  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_maxt,jerr)
794 
795 !!!!!!!!!!!!b_mint(iplas+1)
796  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_mint,jerr)
797 
798 !!!!!!!!!!!!b_db02(iplas+1)
799  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_db02,jerr)
800 
801 
802 !!!!!!!!!!!!b_0db2(iplas+1)
803  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_0db2,jerr)
804 
805 
806 !!!!!!!!!!!!b_db0(iplas+1)
807  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_db0,jerr)
808 
809 !!!!!!!!!!!!tpaps(iplas+1)
810  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,traps,jerr)
811 
812 !!!!!!!!!!!!!!!!!!!!!splining to astra grids
813 
814  rhocn(1)=0.d0
815  rhocn(iplas+1)=1.d0
816  dvdro(1)=0.d0
817  sa(1)=0.d0
818  g22s(1)=0.d0
819  g11s(1)=0.d0
820  fnors(iplas+1)=fvac/(rtor*btor)
821 
822 cp n3spl=iplas
823  n3spl=iplas+1
824 
825 cw write(*,*) (b_mint(i), i=1,iplas)
826 !!!!!!!!!!!!!!!!!!!!! bmint !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
827 
828  CALL e01baf(n3spl,rhocn,b_mint,rrk,cck,
829  * n3spl+4,wrk,6*n3spl+16,ifail)
830 
831  do i=1,na
832  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
833  zspl=rhowrh(i)
834  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
835  bmint(i)=cwk(1)
836  enddo
837  bmint(na1)=b_mint(iplas)
838 
839 cw write(*,*) 'bmint'
840 
841 !!!!!!!!!!!!!!!!!!!!! bmaxt !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
842 
843  CALL e01baf(n3spl,rhocn,b_maxt,rrk,cck,
844  * n3spl+4,wrk,6*n3spl+16,ifail)
845 
846  do i=1,na
847  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
848  zspl=rhowrh(i)
849  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
850  bmaxt(i)=cwk(1)
851  enddo
852  bmaxt(na1)=b_maxt(iplas)
853 
854 cw write(*,*) 'bmaxt'
855 
856 !!!!!!!!!!!!!!!!!!!!! bdb02 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
857 
858  CALL e01baf(n3spl,rhocn,b_db02,rrk,cck,
859  * n3spl+4,wrk,6*n3spl+16,ifail)
860 
861  do i=1,na
862  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
863  zspl=rhowrh(i)
864  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
865  bdb02(i)=cwk(1)
866  enddo
867  bdb02(na1)=b_db02(iplas)
868 
869 cw write(*,*) 'bdb02'
870 
871 !!!!!!!!!!!!!!!!!!!!! b0db2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
872 
873  CALL e01baf(n3spl,rhocn,b_0db2,rrk,cck,
874  * n3spl+4,wrk,6*n3spl+16,ifail)
875 
876  do i=1,na
877  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
878  zspl=rhowrh(i)
879  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
880  b0db2(i)=cwk(1)
881  enddo
882  b0db2(na1)=b_0db2(iplas)
883 
884 cw write(*,*) 'b0db2'
885 
886 !!!!!!!!!!!!!!!!!!!!! bdb0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
887 
888  CALL e01baf(n3spl,rhocn,b_db0,rrk,cck,
889  * n3spl+4,wrk,6*n3spl+16,ifail)
890 
891  do i=1,na
892  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
893  zspl=rhowrh(i)
894  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
895  bdb0(i)=cwk(1)
896  enddo
897  bdb0(na1)=b_db0(iplas)
898 
899 cw write(*,*) 'bdb0'
900 
901 !!!!!!!!!!!!!!!!!!!!!! cu !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
902  n3spl=iplas
903 
904  CALL e01baf(n3spl,rhosn,bj_av,rrk,cck,
905  * n3spl+4,wrk,6*n3spl+16,ifail)
906 
907  do i=1,na
908  zspl=rhowr(i)
909  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
910  cu(i)=cwk(1)/btor/amu0
911  enddo
912  cu(na1)=bj_av(iplas)/btor/amu0
913 
914 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
915 !!!!!!!!!!!!!!!!!!!!!! pres !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
916  n3spl=iplas
917 
918  CALL e01baf(n3spl,rhosn,prs_1,rrk,cck,
919  * n3spl+4,wrk,6*n3spl+16,ifail)
920 
921  do i=1,na
922  zspl=rhowr(i)
923  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
924  pres(i)=cwk(1)
925  enddo
926  pres(na1)=prs_1(iplas)
927 
928 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
929 !!!!!!!!!!!!!!!!!!!!!! fp !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
930  n3spl=iplas
931 
932  CALL e01baf(n3spl,rhosn,psi_1,rrk,cck,
933  * n3spl+4,wrk,6*n3spl+16,ifail)
934 
935  do i=1,na
936  zspl=rhowr(i)
937  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
938  fp(i)=-cwk(1)
939  enddo
940  fp(na1)=-psi_1(iplas)
941 
942 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
943 !!!!!!!!!!!!!!!!!!!!!! W_Dj !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
944  n3spl=iplas
945 
946  CALL e01baf(n3spl,rhosn,wdj,rrk,cck,
947  * n3spl+4,wrk,6*n3spl+16,ifail)
948 
949  do i=1,na
950  zspl=rhowr(i)
951  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
952  w_dj(i)=-cwk(1)
953  enddo
954  w_dj(na1)=wdj(iplas)
955 
956 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
957 
958  !do i=1,na1
959  ! dfpdt(i)=(fp_0(i)-fp(i))/dtim
960  !enddo
961 
962 !!!!!!!!!!!!!!!!!!!!!! slat !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
963  n3spl=iplas
964 
965  CALL e01baf(n3spl,rhosn,sa,rrk,cck,
966  * n3spl+4,wrk,6*n3spl+16,ifail)
967 
968  do i=1,na
969  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
970  zspl=rhowrh(i)
971  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
972  slat(i)=cwk(1)
973  enddo
974  slat(na1)=sa(iplas)
975 
976 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
977  n3spl=iplas+1
978 
979 !!!!!!!!!!!!!!!!!!!!!!! ipol !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
980 
981  CALL e01baf(n3spl,rhocn,fnors,rrk,cck,
982  * n3spl+4,wrk,6*n3spl+16,ifail)
983 
984  do i=1,na1
985  zspl=rhowr(i)
986  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
987  ipol(i)=cwk(1)
988  enddo
989  ipol(na1)=fnors(iplas+1)
990 !!!!!!!!!!!!!!!!!!!!!!! mu !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
991 
992  CALL e01baf(n3spl,rhocn,amus,rrk,cck,
993  * n3spl+4,wrk,6*n3spl+16,ifail)
994 
995  do i=1,na
996  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
997  zspl=rhowrh(i)
998  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
999  mu(i)=cwk(1)
1000  enddo
1001  mu(na1)=amus(iplas+1)
1002 
1003 !!!!!!!!!!!!!!!!!!!!!!! yFOFB !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1004 
1005  CALL e01baf(n3spl,rhocn,traps,rrk,cck,
1006  * n3spl+4,wrk,6*n3spl+16,ifail)
1007 
1008  do i=1,na
1009  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1010  zspl=rhowrh(i)
1011  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1012  yfofb(i)=cwk(1)
1013  enddo
1014  yfofb(na1)=traps(iplas+1)
1015 
1016 
1017 !!!!!!!!!!!!!!!! dvdro !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1018  CALL e01baf(n3spl,rhocn,dvdro,rrk,cck,
1019  * n3spl+4,wrk,6*n3spl+16,ifail)
1020 
1021  do i=1,na1
1022  zspl=rhowr(i)
1023  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1024  vr(i)=cwk(1)
1025  enddo
1026 
1027  do i=1,na
1028  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1029  zspl=rhowrh(i)
1030  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1031  vrs(i)=cwk(1)
1032  enddo
1033  vrs(na1)=dvdro(iplas)
1034 
1035 
1036 !!!!!!!!!!!!!!!!!! g11 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1037 cw write(*,*) 'g11s',(g11s(j),j=1,iplas)
1038 
1039  CALL e01baf(n3spl,rhocn,g11s,rrk,cck,
1040  * n3spl+4,wrk,6*n3spl+16,ifail)
1041 
1042  do i=1,na
1043  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1044  zspl=rhowrh(i)
1045  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1046  g11(i)=cwk(1)
1047  enddo
1048  g11(na1)=g11s(iplas+1)
1049 
1050 
1051 
1052 !!!!!!!!!!!!!!!!!!!! g22 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1053 
1054  CALL e01baf(n3spl,rhocn,g22s,rrk,cck,
1055  * n3spl+4,wrk,6*n3spl+16,ifail)
1056 
1057  do i=1,na
1058  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1059  zspl=rhowrh(i)
1060  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1061  g22(i)=cwk(1)
1062 
1063  enddo
1064 
1065  g22(na1)=g22s(iplas+1)
1066 
1067 !!!!!!!!!!!!!!!!!!!!!! g33 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1068 
1069  CALL e01baf(n3spl,rhocn,g33s,rrk,cck,
1070  * n3spl+4,wrk,6*n3spl+16,ifail)
1071 
1072  do i=1,na1
1073  zspl=rhowr(i)
1074  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1075  g33(i)=cwk(1)
1076  enddo
1077 
1078 
1079 !!!!!!!!!!!!!!!!!!!!!!! gradro !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1080 
1081 
1082  CALL e01baf(n3spl,rhocn,gradrs,rrk,cck,
1083  * n3spl+4,wrk,6*n3spl+16,ifail)
1084 
1085  do i=1,na
1086  !zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1087  zspl=rhowrh(i)
1088  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1089  gradro(i)=cwk(1)
1090  enddo
1091  gradro(na1)=gradrs(iplas+1)
1092 
1093 !!!!!!!!!!!!!!!!!!!!!! droda !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1094 
1095  CALL e01baf(n3spl,rhocn,drodas,rrk,cck,
1096  * n3spl+4,wrk,6*n3spl+16,ifail)
1097 
1098  do i=1,na1
1099  zspl=rhowr(i)
1100  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1101  droda(i)=cwk(1)
1102  enddo
1103 
1104 
1105 cw write(*,*) 'end'
1106 
1107 
1108  return
1109  end
1110 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1111  subroutine comet_(na1,platok,
1112  * rtor,btor,rho,roc,nstep,
1113  * g11,g22,g33,vr,vrs,slat,gradro,rocnew,
1114  * mu,ipol,bmaxt,bmint,bdb02,b0db2,bdb0,droda,
1115  * ybetpl,yli3)
1116 
1117  include 'double.inc'
1118  include 'dim.inc'
1119  parameter(nrpl=nrp+1)
1120  parameter(nrpl4=nrpl+4,nrpl6=nrpl4*6)
1121  include 'compol.inc'
1122  common /combsh/ rm0,zm0,rc0,zc0,asp0,el_up,el_lw,tr_up,tr_lw,nbsh
1123  real*8 rhos(nrp),vols(nrpl),g11s(nrpl),g22s(nrpl),g33s(nrpl)
1124  real*8 gradrs(nrpl),amus(nrpl),fnors(nrpl),drodas(nrpl)
1125  real*8 dvdro(nrpl),sa(nrp),rhocs(nrpl),rhocn(nrpl),rhosn(nrpl)
1126  real*8 b_maxt(nrpl),b_mint(nrpl),b_db02(nrpl),b_db0(nrpl),
1127  * b_0db2(nrpl),d_roda(nrpl)
1128  real*8 rrk(nrpl4),cck(nrpl4),wrk(nrpl6)
1129  real*8 cwk(4)
1130  real*8 rhowr(1000)
1131  real*8 rtor,btor,rho(*),roc,ybetpl,yli3,
1132  * g11(na1),g22(na1),g33(na1),
1133  * vr(na1),vrs(na1),slat(na1),gradro(na1),rocnew,droda(na1),
1134  * mu(na1),ipol(na1),bmaxt(na1),bmint(na1),
1135  * bdb02(na1),b0db2(na1),bdb0(na1)
1136 
1137  sqrt(xx)=dsqrt(xx)
1138 
1139  na=na1-1
1140 
1141  platok=tokp
1142  do i=1,iplas
1143  !rhos(i)=sqrt(flucf(i)/(pi*btor))
1144  rhos(i)=sqrt(flx_fi(i)/(pi*btor))
1145  enddo
1146 cp
1147 cp
1148  rhocs(1)=0.d0
1149  rhos(1)=0.d0
1150 
1151  do i=2,iplas
1152  rhocs(i)=0.5d0*(rhos(i-1)+rhos(i))
1153  enddo
1154 
1155  rocnew=rhos(iplas)
1156  rhocs(iplas+1)=rhos(iplas)
1157 
1158 c=============================andrey=02.03.2003vvvvvvvvvvvvv
1159  arc0=rhocs(iplas+1)
1160  arc1=rhocs(iplas)
1161  arc2=rhocs(iplas-1)
1162  arc3=rhocs(iplas-2)
1163  spw1=(arc0-arc2)*(arc0-arc3)
1164  & /(arc1-arc2)/(arc1-arc3)
1165  spw2=(arc0-arc1)*(arc0-arc3)
1166  & /(arc2-arc1)/(arc2-arc3)
1167  spw3=(arc0-arc1)*(arc0-arc2)
1168  & /(arc3-arc1)/(arc3-arc2)
1169  q(iplas)=spw1*q(iplas-1)+spw2*q(iplas-2)+spw3*q(iplas-3)
1170 
1171 c=============================andrey=02.03.2003^^^^^^^^^^^^^^
1172 
1173 
1174 
1175  do i=1,iplas+1
1176  rhocn(i)=rhocs(i)/rhos(iplas)
1177  enddo
1178  rhocn(1)=0.d0
1179 
1180  do i=1,iplas
1181  rhosn(i)=rhos(i)/rhos(iplas)
1182  drodas(i)=rhosn(i)
1183  enddo
1184  rhosn(1)=0.d0
1185 cp
1186 
1187  rhocn(iplas+1)=1.d0
1188  rhosn(iplas)=1.d0
1189 
1190  do i=1,na
1191  rhowr(i)=rho(i)/roc
1192  enddo
1193  rhowr(na1)=1.d0-1.d-9
1194 
1195  sa(1)=0.d0
1196 cw write(*,*) (rhocn(j),j=1,iplas)
1197 cw write(*,*) (rhosn(j),j=1,iplas)
1198 
1199  do i=1,iplas-1
1200 
1201  dsqi=0.d0
1202  dvoli=0.d0
1203  av_r2=0.d0
1204  av_gr1=0.d0
1205  av_gr2=0.d0
1206  av_grr2=0.d0
1207 
1208  u1=rhos(i)
1209  u2=rhos(i+1)
1210  u3=rhos(i+1)
1211  u4=rhos(i)
1212 
1213  do j=2,nt1
1214 
1215  dvoli=dvoli+vol(i,j)
1216  dsqi=dsqi+sr(i+1,j)
1217 
1218  r1=r(i,j)
1219  r2=r(i+1,j)
1220  r3=r(i+1,j+1)
1221  r4=r(i,j+1)
1222  r0=(r1+r2+r3+r4)*0.25d0
1223 
1224  z1=z(i,j)
1225  z2=z(i+1,j)
1226  z3=z(i+1,j+1)
1227  z4=z(i,j+1)
1228 
1229  sss=s(i,j)
1230 
1231  dadr= 0.5d0*((u1+u2)*(z2-z1)+(u2+u3)*(z3-z2)+
1232  * (u3+u4)*(z4-z3)+(u4+u1)*(z1-z4))/sss
1233 
1234  dadz=-0.5d0*((u1+u2)*(r2-r1)+(u2+u3)*(r3-r2)+
1235  * (u3+u4)*(r4-r3)+(u4+u1)*(r1-r4))/sss
1236 
1237  gr2=dadr**2+dadz**2
1238 
1239  av_gr1=av_gr1+vol(i,j)*sqrt(gr2)
1240  av_gr2=av_gr2+vol(i,j)*gr2
1241  av_grr2=av_grr2+vol(i,j)*gr2/r0**2
1242  av_r2=av_r2+vol(i,j)/r0**2
1243 
1244  enddo
1245  sa(i+1)=dsqi*2.d0*pi
1246 c vols(i+1)=dvoli
1247  gradrs(i+1)=av_gr1/dvoli
1248  amus(i+1)=(2.d0*pi)/q(i)
1249  fnors(i+1)=f(i)/(rtor*btor)
1250 
1251  dvdro(i+1)=2.d0*pi*dvoli/(rhos(i+1)-rhos(i))
1252 cc g11s(i+1)=av_gr2/dvoli
1253 cc g22s(i+1)=av_grr2/(dvoli*4.d0*pi*pi)*rtor
1254  g11s(i+1)=av_gr2*dvdro(i+1)/dvoli
1255  g22s(i+1)=av_grr2*dvdro(i+1)/(dvoli*4.d0*pi*pi)*rtor/fnors(i+1)
1256  g33s(i+1)=av_r2*rtor**2/dvoli
1257  enddo
1258 
1259 cccc magnetic field
1260 
1261 cw write(*,*) 'before b'
1262 
1263  bp2av=0.d0
1264  bp2=0.d0
1265 cc do i=2,iplas-1
1266  do i=1,iplas-1
1267 
1268  dvoli=0.d0
1269  bmx=0.d0
1270  bmn=1.d9*btor
1271  bavr=0.d0
1272  bavr2=0.d0
1273  bavrd2=0.d0
1274 
1275  do j=2,nt1
1276 
1277  dvoli=dvoli+vol(i,j)
1278 
1279  r1=r(i,j)
1280  r2=r(i+1,j)
1281  r3=r(i+1,j+1)
1282  r4=r(i,j+1)
1283  r0=(r1+r2+r3+r4)*0.25d0
1284 
1285  bm_pol=psim*(psia(i+1)-psia(i))/st(i,j)
1286  bp_pol=psim*(psia(i+1)-psia(i))/st(i,j+1)
1287 
1288  bpol2v=bm_pol**2*( vol1(i,j)/sin1(i,j)+vol2(i,j)/sin2(i,j))+
1289  + bp_pol**2*( vol3(i,j)/sin3(i,j)+vol4(i,j)/sin4(i,j))
1290 
1291  bpol2=bpol2v /vol(i,j)
1292  b_fi2=(f(i)/r0)**2
1293 
1294  b_tot2=bpol2+b_fi2
1295  b_tot=dsqrt(b_tot2)
1296 
1297  bmx=dmax1(b_tot,bmx)
1298  bmn=dmin1(b_tot,bmn)
1299 
1300  bavr2=bavr2+b_tot2*vol(i,j)
1301  bavr =bavr +b_tot *vol(i,j)
1302  bavrd2=bavrd2+vol(i,j)/(b_tot2+1.d-8)
1303  bp2av=bp2av+bpol2v
1304  bp2=bp2+bpol2v
1305  enddo
1306  b_maxt(i+1)=bmx
1307  b_mint(i+1)=bmn
1308  b_db02(i+1)=bavr2/(dvoli*btor**2)
1309  b_0db2(i+1)=(bavrd2*btor**2)/(dvoli)
1310  b_db0(i+1)=bavr/(dvoli*btor)
1311 
1312 
1313  enddo
1314  b_maxt(1)=.5d0*(b_maxt(2)+b_mint(2))
1315  b_mint(1)=b_maxt(1)
1316  b_db02(1)=b_db02(2)
1317  b_0db2(1)=b_0db2(2)
1318  b_db0(1)=b_db0(2)
1319 
1320 
1321 c yli3 =4.d0*pi*bp2av/(rtor*tokp*tokp)
1322  yli3 = 4.d0*pi*bp2/(rtor*(.4d0*pi*tokp)**2)
1323 c*new
1324 
1325 c goto 9991
1326  ip0=1
1327  ip1=2
1328  ip2=3
1329  ip3=4
1330 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1331 
1332  x3=rhocs(ip3)
1333  x2=rhocs(ip2)
1334  x1=rhocs(ip1)
1335  x0=rhocs(ip0)
1336 !!!!!!!!!!!!g11(1)
1337 call b_extrp(x0, x1,x2,x3, ip1,ip2,ip3,ip0,g11s,jerr)
1338 !!!!!!!!!!!!g22(1)
1339 call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,g22s,jerr)
1340 
1341 !!!!!!!!!!!!g33(1)
1342  call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,g33s,jerr)
1343 
1344 !!!!!!!!!!!!fnors(1)
1345  call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,fnors,jerr)
1346 
1347 !!!!!!!!!!!!gradrs(1)
1348  call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,gradrs,jerr)
1349 
1350 !!!!!!!!!!!!amus(1)
1351  call b_extrp(x0,x1,x2,x3, ip1,ip2,ip3,ip0,amus,jerr)
1352 
1353 
1354 ! x0=rhocs(1)
1355 ! x1=rhocs(2)
1356 ! x2=rhocs(3)
1357 ! x3=rhocs(4)
1358 
1359 !
1360  9991 continue
1361 
1362  ip0=iplas+1
1363  ip1=iplas-2
1364  ip2=iplas-1
1365  ip3=iplas
1366 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1367 
1368  x1=rhocs(ip3)
1369  x2=rhocs(ip2)
1370  x3=rhocs(ip1)
1371  x0=rhocs(ip0)
1372 
1373 
1374 cw write(*,*) 'passed 1'
1375 !!!!!!!!!!!gradrs(iplas+1)
1376  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,gradrs,jerr)
1377 
1378 !!!!!!!!!!!amus(iplas+1)
1379  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,amus,jerr)
1380 
1381 
1382 !!!!!!!!!!!dvdro(iplas+1)
1383  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,dvdro,jerr)
1384 
1385 
1386 !!!!!!!!!!!!g11(iplas+1)
1387  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,g11s,jerr)
1388 
1389 !!!!!!!!!!!!g22(iplas+1)
1390  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,g22s,jerr)
1391 
1392 !!!!!!!!!!!!g33(iplas+1)
1393  call b_extrp(x0,x1,x2,x3, ip3,ip2,ip1,ip0,g33s,jerr)
1394 
1395 !!!!!!!!!!!!b_maxt(iplas+1)
1396  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_maxt,jerr)
1397 
1398 !!!!!!!!!!!!b_mint(iplas+1)
1399  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_mint,jerr)
1400 
1401 !!!!!!!!!!!!b_db02(iplas+1)
1402  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_db02,jerr)
1403 
1404 
1405 !!!!!!!!!!!!b_0db2(iplas+1)
1406  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_0db2,jerr)
1407 
1408 
1409 !!!!!!!!!!!!b_db0(iplas+1)
1410  call b_extrp(x0, x1,x2,x3, ip3,ip2,ip1,ip0,b_db0,jerr)
1411 
1412 
1413 !!!!!!!!!!!!!!!!!!!!!splining to astra grids
1414 
1415  rhocn(1)=0.d0
1416  rhocn(iplas+1)=1.d0
1417  dvdro(1)=0.d0
1418  sa(1)=0.d0
1419  g22s(1)=0.d0
1420  g11s(1)=0.d0
1421  fnors(iplas+1)=fvac/(rtor*btor)
1422 
1423 cp n3spl=iplas
1424  n3spl=iplas+1
1425 
1426 cw write(*,*) (b_mint(i), i=1,iplas)
1427 !!!!!!!!!!!!!!!!!!!!! bmint !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1428 
1429  CALL e01baf(n3spl,rhocn,b_mint,rrk,cck,
1430  * n3spl+4,wrk,6*n3spl+16,ifail)
1431 
1432  do i=1,na
1433  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1434  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1435  bmint(i)=cwk(1)
1436  enddo
1437  bmint(na1)=b_mint(iplas)
1438 
1439 cw write(*,*) 'bmint'
1440 
1441 !!!!!!!!!!!!!!!!!!!!! bmaxt !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1442 
1443  CALL e01baf(n3spl,rhocn,b_maxt,rrk,cck,
1444  * n3spl+4,wrk,6*n3spl+16,ifail)
1445 
1446  do i=1,na
1447  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1448  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1449  bmaxt(i)=cwk(1)
1450  enddo
1451  bmaxt(na1)=b_maxt(iplas)
1452 
1453 cw write(*,*) 'bmaxt'
1454 
1455 !!!!!!!!!!!!!!!!!!!!! bdb02 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1456 
1457  CALL e01baf(n3spl,rhocn,b_db02,rrk,cck,
1458  * n3spl+4,wrk,6*n3spl+16,ifail)
1459 
1460  do i=1,na
1461  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1462  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1463  bdb02(i)=cwk(1)
1464  enddo
1465  bdb02(na1)=b_db02(iplas)
1466 
1467 cw write(*,*) 'bdb02'
1468 
1469 !!!!!!!!!!!!!!!!!!!!! b0db2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1470 
1471  CALL e01baf(n3spl,rhocn,b_0db2,rrk,cck,
1472  * n3spl+4,wrk,6*n3spl+16,ifail)
1473 
1474  do i=1,na
1475  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1476  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1477  b0db2(i)=cwk(1)
1478  enddo
1479  b0db2(na1)=b_0db2(iplas)
1480 
1481 cw write(*,*) 'b0db2'
1482 
1483 !!!!!!!!!!!!!!!!!!!!! bdb0 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1484 
1485  CALL e01baf(n3spl,rhocn,b_db0,rrk,cck,
1486  * n3spl+4,wrk,6*n3spl+16,ifail)
1487 
1488  do i=1,na
1489 cp zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1490  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1491  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1492  bdb0(i)=cwk(1)
1493  enddo
1494  bdb0(na1)=b_db0(iplas)
1495 
1496 cw write(*,*) 'bdb0'
1497 
1498 
1499 !!!!!!!!!!!!!!!!!!!!!! slat !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1500  n3spl=iplas
1501 
1502  CALL e01baf(n3spl,rhosn,sa,rrk,cck,
1503  * n3spl+4,wrk,6*n3spl+16,ifail)
1504 
1505  do i=1,na
1506  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1507  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1508  slat(i)=cwk(1)
1509  enddo
1510  slat(na1)=sa(iplas)
1511 
1512 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1513  n3spl=iplas+1
1514 
1515 !!!!!!!!!!!!!!!!!!!!!!! ipol !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1516 
1517  CALL e01baf(n3spl,rhocn,fnors,rrk,cck,
1518  * n3spl+4,wrk,6*n3spl+16,ifail)
1519 
1520  do i=1,na1
1521  zspl=rhowr(i)
1522  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1523  ipol(i)=cwk(1)
1524  enddo
1525  ipol(na1)=fnors(iplas+1)
1526 !!!!!!!!!!!!!!!!!!!!!!! mu !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1527 
1528  CALL e01baf(n3spl,rhocn,amus,rrk,cck,
1529  * n3spl+4,wrk,6*n3spl+16,ifail)
1530 
1531  do i=1,na
1532  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1533  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1534  mu(i)=cwk(1)
1535  enddo
1536  mu(na1)=amus(iplas+1)
1537 
1538 
1539 
1540 !!!!!!!!!!!!!!!! dvdro !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1541  CALL e01baf(n3spl,rhocn,dvdro,rrk,cck,
1542  * n3spl+4,wrk,6*n3spl+16,ifail)
1543 
1544  do i=1,na1
1545  zspl=rhowr(i)
1546  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1547  vr(i)=cwk(1)
1548  enddo
1549 
1550  do i=1,na
1551  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1552  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1553  vrs(i)=cwk(1)
1554  enddo
1555  vrs(na1)=dvdro(iplas)
1556 
1557 
1558 !!!!!!!!!!!!!!!!!! g11 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1559 cw write(*,*) 'g11s',(g11s(j),j=1,iplas)
1560 
1561  CALL e01baf(n3spl,rhocn,g11s,rrk,cck,
1562  * n3spl+4,wrk,6*n3spl+16,ifail)
1563 
1564  do i=1,na
1565  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1566  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1567  g11(i)=cwk(1)
1568  enddo
1569  g11(na1)=g11s(iplas+1)
1570 
1571 
1572 
1573 !!!!!!!!!!!!!!!!!!!! g22 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1574 
1575  CALL e01baf(n3spl,rhocn,g22s,rrk,cck,
1576  * n3spl+4,wrk,6*n3spl+16,ifail)
1577 
1578  do i=1,na
1579  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1580  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1581  g22(i)=cwk(1)
1582 
1583  enddo
1584 
1585  g22(na1)=g22s(iplas+1)
1586 
1587 !!!!!!!!!!!!!!!!!!!!!! g33 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1588 
1589  CALL e01baf(n3spl,rhocn,g33s,rrk,cck,
1590  * n3spl+4,wrk,6*n3spl+16,ifail)
1591 
1592  do i=1,na1
1593  zspl=rhowr(i)
1594  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1595  g33(i)=cwk(1)
1596  enddo
1597 
1598 
1599 !!!!!!!!!!!!!!!!!!!!!!! gradro !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1600 
1601 
1602  CALL e01baf(n3spl,rhocn,gradrs,rrk,cck,
1603  * n3spl+4,wrk,6*n3spl+16,ifail)
1604 
1605  do i=1,na
1606  zspl=(rhowr(i)+rhowr(i+1))*0.5d0
1607  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1608  gradro(i)=cwk(1)
1609  enddo
1610  gradro(na1)=gradrs(iplas+1)
1611 
1612 !!!!!!!!!!!!!!!!!!!!!! droda !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1613 
1614  CALL e01baf(n3spl,rhocn,drodas,rrk,cck,
1615  * n3spl+4,wrk,6*n3spl+16,ifail)
1616 
1617  do i=1,na1
1618  zspl=rhowr(i)
1619  CALL e02bcf(n3spl+4,rrk,cck,zspl,0,cwk,ifail)
1620  droda(i)=cwk(1)
1621  enddo
1622 
1623 
1624 cw write(*,*) 'end'
1625 
1626 
1627  return
1628  end
1629 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1630  subroutine profro(na1,neql,rtor,eqpf,eqff,rho,roc)
1631 c 29 jan 2003
1632  include 'double.inc'
1633  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
1634  include 'dim.inc'
1635  include 'compol.inc'
1636 
1637  real*8 rosn(nrp),ppr(nursp),ffpr(nursp),rhow(nursp)
1638  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
1639  real*8 cwk(4)
1640 
1641  real*8 eqpf(*),eqff(*),rho(*),rtor,roc
1642 
1643  !return
1644 
1645  iplas=neql
1646  nspl=na1+1
1647  do i=1,iplas
1648  rosn(i)=(i-1.d0)/(iplas-1.d0)
1649  enddo
1650 
1651 !! exstrapolation on magn. axis
1652 
1653  rh0=0.d0
1654 
1655  rh1=rho(1)
1656  rh2=rho(2)
1657  rh3=rho(3)
1658 
1659  ppx1=eqpf(1)
1660  ppx2=eqpf(2)
1661  ppx3=eqpf(3)
1662 
1663  ffx1=eqff(1)
1664  ffx2=eqff(2)
1665  ffx3=eqff(3)
1666 
1667  call extrp2(rh0,ppx0, rh1,rh2,rh3, ppx1,ppx2,ppx3)
1668  call extrp2(rh0,ffx0, rh1,rh2,rh3, ffx1,ffx2,ffx3)
1669 
1670  ppr(1)=ppx0
1671  ffpr(1)=ffx0
1672  rhow(1)=0.d0
1673 
1674  do i=2,na1+1
1675  ppr(i)=eqpf(i-1)
1676  ffpr(i)=eqff(i-1)
1677  rhow(i)=rho(i-1)/roc
1678  enddo
1679 
1680  dpdpsi(1)=ppr(1)
1681  dfdpsi(1)=ffpr(1)
1682 
1683  dpdpsi(iplas)=ppr(na1+1)
1684  dfdpsi(iplas)=ffpr(na1+1)
1685 
1686 
1687  CALL e01baf(nspl,rhow,ppr,rrk,cck,
1688  * nspl+4,wrk,6*nspl+16,ifail)
1689 
1690  do i=2,iplas-1
1691 
1692  zrho=rosn(i)
1693  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
1694  dpdpsi(i)=cwk(1)
1695 
1696  enddo
1697 
1698 
1699  CALL e01baf(nspl,rhow,ffpr,rrk,cck,
1700  * nspl+4,wrk,6*nspl+16,ifail)
1701 
1702  do i=2,iplas-1
1703 
1704  zrho=rosn(i)
1705  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
1706  dfdpsi(i)=cwk(1)
1707 
1708  enddo
1709 
1710  do i=1,iplas
1711  dpdpsi(i)=dpdpsi(i)/rtor
1712  dfdpsi(i)=dfdpsi(i)*rtor
1713  enddo
1714 
1715  return
1716  end
1717 
1718 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1719 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1720  subroutine profro_p(na1,neql,rtor,eqpf,rho,roc)
1721 c march 2007
1722  include 'double.inc'
1723  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
1724  include 'dim.inc'
1725  include 'compol.inc'
1726 
1727  real*8 rosn(nrp),rhot(nursp),ppr(nursp),ffpr(nursp),rhow(nursp)
1728  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
1729  real*8 cwk(4)
1730 
1731  real*8 eqpf(*),rho(*),rtor,roc
1732 
1733  iplas=neql
1734  nspl=na1+1
1735  do i=1,iplas
1736  rosn(i)=(i-1.d0)/(iplas-1.d0)
1737  enddo
1738 
1739 !! exstrapolation on magn. axis
1740 
1741  rh0=0.d0
1742 
1743  rh1=rho(1)
1744  rh2=rho(2)
1745  rh3=rho(3)
1746 
1747  ppx1=eqpf(1)
1748  ppx2=eqpf(2)
1749  ppx3=eqpf(3)
1750 
1751  !ffx1=eqff(1)
1752  !ffx2=eqff(2)
1753  !ffx3=eqff(3)
1754 
1755  call extrp2(rh0,ppx0, rh1,rh2,rh3, ppx1,ppx2,ppx3)
1756  !call EXTRP2(rh0,ffx0, rh1,rh2,rh3, ffx1,ffx2,ffx3)
1757 
1758  ppr(1)=ppx0
1759  !ffpr(1)=ffx0
1760  rhow(1)=0.d0
1761 
1762  do i=2,na1+1
1763  ppr(i)=eqpf(i-1)
1764  !ffpr(i)=eqff(i-1)
1765  rhow(i)=rho(i-1)/roc
1766  enddo
1767 
1768  dpdpsi(1)=ppr(1)
1769  !dfdpsi(1)=ffpr(1)
1770 
1771  dpdpsi(iplas)=ppr(na1+1)
1772  !dfdpsi(iplas)=ffpr(na1+1)
1773 
1774 
1775  CALL e01baf(nspl,rhow,ppr,rrk,cck,
1776  * nspl+4,wrk,6*nspl+16,ifail)
1777 
1778  do i=2,iplas-1
1779 
1780  zrho=rosn(i)
1781  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
1782  dpdpsi(i)=cwk(1)
1783 
1784  enddo
1785 
1786 
1787 ! CALL E01BAF(nspl,rhow,ffpr,RRK,CCK,
1788 ! * nspl+4,WRK,6*nspl+16,IFAIL)
1789 !
1790 ! do i=2,iplas-1
1791 !
1792 ! zrho=rosn(i)
1793 ! CALL E02BCF(nspl+4,RRK,CCK,zrho,0,CWk,IFAIL)
1794 ! dfdpsi(i)=cwk(1)
1795 !
1796 ! enddo
1797 
1798  do i=1,iplas
1799  dpdpsi(i)=dpdpsi(i)/rtor
1800  !dfdpsi(i)=dfdpsi(i)*rtor
1801  enddo
1802 
1803  return
1804  end
1805 
1806 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1807 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1808  subroutine profro_p_l(na1,neql,rtor,eqpf,rho,roc)
1809 c march 2007
1810  include 'double.inc'
1811  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
1812  include 'dim.inc'
1813  include 'compol.inc'
1814 
1815  real*8 rosn(nrp),rhot(nursp),ppr(nursp),ffpr(nursp),rhow(nursp)
1816  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
1817  real*8 cwk(4)
1818 
1819  real*8 eqpf(*),rho(*),rtor,roc
1820 
1821  iplas=neql
1822  nspl=na1+1
1823  do i=1,iplas
1824  rosn(i)=(i-1.d0)/(iplas-1.d0)
1825  enddo
1826 
1827 !! exstrapolation on magn. axis
1828 
1829  rh0=0.d0
1830 
1831  rh1=rho(1)
1832  rh2=rho(2)
1833  rh3=rho(3)
1834 
1835  ppx1=eqpf(1)
1836  ppx2=eqpf(2)
1837  ppx3=eqpf(3)
1838 
1839  !ffx1=eqff(1)
1840  !ffx2=eqff(2)
1841  !ffx3=eqff(3)
1842 
1843  call extrp2(rh0,ppx0, rh1,rh2,rh3, ppx1,ppx2,ppx3)
1844  !call EXTRP2(rh0,ffx0, rh1,rh2,rh3, ffx1,ffx2,ffx3)
1845 
1846  ppr(1)=ppr(2)
1847  !ppr(1)=ppx0
1848  !ffpr(1)=ffx0
1849  rhow(1)=0.d0
1850 
1851  do i=2,na1+1
1852  ppr(i)=eqpf(i-1)
1853  !ffpr(i)=eqff(i-1)
1854  rhow(i)=rho(i-1)/roc
1855  enddo
1856 
1857  dpdpsi(1)=ppr(1)
1858  !dfdpsi(1)=ffpr(1)
1859 
1860  dpdpsi(iplas)=ppr(na1+1)
1861  !dfdpsi(iplas)=ffpr(na1+1)
1862 
1863 
1864 !linear interpolation
1865 
1866  do i=2,iplas-1
1867  zrho=rosn(i)
1868 
1869  do ic=1,nspl-1
1870  if(zrho.le.rhow(ic+1) .AnD. zrho.gt.rhow(ic)) then
1871  dpdpsi(i)=( ppr(ic)*(rhow(ic+1)-zrho)
1872  * +ppr(ic+1)*(zrho-rhow(ic)) )
1873  * /(rhow(ic+1)-rhow(ic))
1874  exit
1875  endif
1876  enddo
1877 
1878 
1879  enddo
1880 
1881 
1882 
1883 
1884 
1885 
1886 ! CALL E01BAF(nspl,rhow,ppr,RRK,CCK,
1887 ! * nspl+4,WRK,6*nspl+16,IFAIL)
1888 
1889 ! do i=2,iplas-1
1890 !
1891 ! zrho=rosn(i)
1892 ! CALL E02BCF(nspl+4,RRK,CCK,zrho,0,CWk,IFAIL)
1893 ! dpdpsi(i)=cwk(1)
1894 !
1895 ! enddo
1896 
1897 
1898 ! CALL E01BAF(nspl,rhow,ffpr,RRK,CCK,
1899 ! * nspl+4,WRK,6*nspl+16,IFAIL)
1900 !
1901 ! do i=2,iplas-1
1902 !
1903 ! zrho=rosn(i)
1904 ! CALL E02BCF(nspl+4,RRK,CCK,zrho,0,CWk,IFAIL)
1905 ! dfdpsi(i)=cwk(1)
1906 !
1907 ! enddo
1908 
1909  do i=1,iplas
1910  dpdpsi(i)=dpdpsi(i)/rtor
1911  !dfdpsi(i)=dfdpsi(i)*rtor
1912  enddo
1913 
1914  return
1915  end
1916 
1917 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1918 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1919  subroutine profro_sbc(na1,neql,rtor,btor,cc,Te,cubs,cd,rho,roc)
1920 c 03 jan 2008
1921  include 'double.inc'
1922  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
1923  include 'dim.inc'
1924  include 'compol.inc'
1925  common /com_sigcd/ c_sig(nrp),t_el(nrp),c_bts(nrp),c_driv(nrp)
1926 
1927  real*8 rosn(nrp),funw(nursp),rhow(nursp)
1928  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
1929  real*8 cwk(4)
1930 
1931  real*8 rho(*),rtor,btor,roc,cc(*),te(*),cubs(*),cd(*)
1932 
1933  iplas=neql
1934  nspl=na1+1
1935  do i=1,iplas
1936  rosn(i)=(i-1.d0)/(iplas-1.d0)
1937  enddo
1938 
1939 !! exstrapolation on magn. axis
1940 
1941  rh0=0.d0
1942 
1943  rh1=rho(1)
1944  rh2=rho(2)
1945  rh3=rho(3)
1946 
1947  ccx1=cc(1)
1948  ccx2=cc(2)
1949  ccx3=cc(3)
1950 
1951  cdx1=cd(1)
1952  cdx2=cd(2)
1953  cdx3=cd(3)
1954 
1955  cbx1=cubs(1)
1956  cbx2=cubs(2)
1957  cbx3=cubs(3)
1958 
1959  cex1=te(1)
1960  cex2=te(2)
1961  cex3=te(3)
1962 
1963  call extrp2(rh0,ccx0, rh1,rh2,rh3, ccx1,ccx2,ccx3)
1964  call extrp2(rh0,cdx0, rh1,rh2,rh3, cdx1,cdx2,cdx3)
1965  call extrp2(rh0,cbx0, rh1,rh2,rh3, cbx1,cbx2,cbx3)
1966  call extrp2(rh0,cex0, rh1,rh2,rh3, cex1,cex2,cex3)
1967 
1968 !!!!!!!!!!!!!!!!!!!cc!!!!!!!!!!!!!!!!!!!!!!!
1969 
1970  funw(1)=ccx0
1971  rhow(1)=0.d0
1972 
1973  do i=2,na1+1
1974  funw(i)=cc(i-1)
1975  rhow(i)=rho(i-1)/roc
1976  enddo
1977 
1978  c_sig(1)=funw(1)
1979  c_sig(iplas)=funw(na1+1)
1980 
1981  CALL e01baf(nspl,rhow,funw,rrk,cck,
1982  * nspl+4,wrk,6*nspl+16,ifail)
1983 
1984  do i=2,iplas-1
1985 
1986  zrho=rosn(i)
1987  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
1988  c_sig(i)=cwk(1)
1989 
1990  enddo
1991 
1992  !!!!!!!!!!!!!!!!!!!cc!!!!!!!!!!!!!!!!!!!!!!!
1993 
1994 !!!!!!!!!!!!!!!!!!!cd!!!!!!!!!!!!!!!!!!!!!!!
1995 
1996  funw(1)=cdx0
1997  rhow(1)=0.d0
1998 
1999  do i=2,na1+1
2000  funw(i)=cd(i-1)
2001  rhow(i)=rho(i-1)/roc
2002  enddo
2003 
2004  c_driv(1)=funw(1)
2005  c_driv(iplas)=funw(na1+1)
2006 
2007  CALL e01baf(nspl,rhow,funw,rrk,cck,
2008  * nspl+4,wrk,6*nspl+16,ifail)
2009 
2010  do i=2,iplas-1
2011 
2012  zrho=rosn(i)
2013  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
2014  c_driv(i)=cwk(1)
2015 
2016  enddo
2017 
2018  !!!!!!!!!!!!!!!!!!!cd!!!!!!!!!!!!!!!!!!!!!!!
2019 
2020  !!!!!!!!!!!!!!!!!!!cubs!!!!!!!!!!!!!!!!!!!!!!!
2021 
2022  funw(1)=cbx0
2023  rhow(1)=0.d0
2024 
2025  do i=2,na1+1
2026  funw(i)=cubs(i-1)
2027  rhow(i)=rho(i-1)/roc
2028  enddo
2029 
2030  c_bts(1)=funw(1)
2031  c_bts(iplas)=funw(na1+1)
2032 
2033  CALL e01baf(nspl,rhow,funw,rrk,cck,
2034  * nspl+4,wrk,6*nspl+16,ifail)
2035 
2036  do i=2,iplas-1
2037 
2038  zrho=rosn(i)
2039  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
2040  c_bts(i)=cwk(1)
2041 
2042  enddo
2043 
2044  !!!!!!!!!!!!!!!!!!!cubs!!!!!!!!!!!!!!!!!!!!!!!
2045 
2046  !!!!!!!!!!!!!!!!!!!te!!!!!!!!!!!!!!!!!!!!!!!
2047 
2048  funw(1)=cex0
2049  rhow(1)=0.d0
2050 
2051  do i=2,na1+1
2052  funw(i)=te(i-1)
2053  rhow(i)=rho(i-1)/roc
2054  enddo
2055 
2056  t_el(1)=funw(1)
2057  t_el(iplas)=funw(na1+1)
2058 
2059  CALL e01baf(nspl,rhow,funw,rrk,cck,
2060  * nspl+4,wrk,6*nspl+16,ifail)
2061 
2062  do i=2,iplas-1
2063 
2064  zrho=rosn(i)
2065  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
2066  t_el(i)=cwk(1)
2067 
2068  enddo
2069 
2070  !!!!!!!!!!!!!!!!!!!te!!!!!!!!!!!!!!!!!!!!!!!
2071 
2072  do i=1,iplas
2073  c_bts(i)=c_bts(i)*btor
2074  c_driv(i)=c_driv(i)*btor
2075  t_el(i)=t_el(i)*1.d3
2076  enddo
2077 
2078  return
2079  end
2080 
2081 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2082  subroutine profro_sbc_l(na1,neql,rtor,btor,cc,Te,cubs,cd,rho,roc)
2083 
2084  include 'double.inc'
2085  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
2086  include 'dim.inc'
2087  include 'compol.inc'
2088  common /com_sigcd/ c_sig(nrp),t_el(nrp),c_bts(nrp),c_driv(nrp)
2089 
2090  real*8 rosn(nrp),funw(nursp),rhow(nursp)
2091  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
2092  real*8 cwk(4)
2093 
2094  real*8 rho(*),rtor,btor,roc,cc(*),te(*),cubs(*),cd(*)
2095 
2096  iplas=neql
2097  nspl=na1+1
2098  do i=1,iplas
2099  rosn(i)=(i-1.d0)/(iplas-1.d0)
2100  enddo
2101 
2102 !! exstrapolation on magn. axis
2103 
2104  rh0=0.d0
2105 
2106  rh1=rho(1)
2107  rh2=rho(2)
2108  rh3=rho(3)
2109 
2110  ccx1=cc(1)
2111  ccx2=cc(2)
2112  ccx3=cc(3)
2113 
2114  cdx1=cd(1)
2115  cdx2=cd(2)
2116  cdx3=cd(3)
2117 
2118  cbx1=cubs(1)
2119  cbx2=cubs(2)
2120  cbx3=cubs(3)
2121 
2122  cex1=te(1)
2123  cex2=te(2)
2124  cex3=te(3)
2125 
2126  call extrp2(rh0,ccx0, rh1,rh2,rh3, ccx1,ccx2,ccx3)
2127  call extrp2(rh0,cdx0, rh1,rh2,rh3, cdx1,cdx2,cdx3)
2128  call extrp2(rh0,cbx0, rh1,rh2,rh3, cbx1,cbx2,cbx3)
2129  call extrp2(rh0,cex0, rh1,rh2,rh3, cex1,cex2,cex3)
2130 
2131 !!!!!!!!!!!!!!!!!!!cc!!!!!!!!!!!!!!!!!!!!!!!
2132 
2133  !funw(1)=ccx0
2134  funw(1)=ccx1
2135  rhow(1)=0.d0
2136 
2137  do i=2,na1+1
2138  funw(i)=cc(i-1)
2139  rhow(i)=rho(i-1)/roc
2140  enddo
2141 
2142  c_sig(1)=funw(1)
2143  c_sig(iplas)=funw(na1+1)
2144 
2145 ! CALL E01BAF(nspl,rhow,funw,RRK,CCK,
2146 ! * nspl+4,WRK,6*nspl+16,IFAIL)
2147 !
2148 ! do i=2,iplas-1
2149 !
2150 ! zrho=rosn(i)
2151 ! CALL E02BCF(nspl+4,RRK,CCK,zrho,0,CWk,IFAIL)
2152 ! C_sig(i)=cwk(1)
2153 !
2154 ! enddo
2155 
2156  nspl=na1+1
2157  call linsplin(nspl,funw,rhow, iplas,c_sig,rosn)
2158 
2159  !!!!!!!!!!!!!!!!!!!cc!!!!!!!!!!!!!!!!!!!!!!!
2160 
2161 !!!!!!!!!!!!!!!!!!!cd!!!!!!!!!!!!!!!!!!!!!!!
2162 
2163  !funw(1)=cdx0
2164  funw(1)=cdx1
2165  rhow(1)=0.d0
2166 
2167  do i=2,na1+1
2168  funw(i)=cd(i-1)
2169  rhow(i)=rho(i-1)/roc
2170  enddo
2171 
2172  c_driv(1)=funw(1)
2173  c_driv(iplas)=funw(na1+1)
2174 
2175 ! CALL E01BAF(nspl,rhow,funw,RRK,CCK,
2176 ! * nspl+4,WRK,6*nspl+16,IFAIL)
2177 !
2178 ! do i=2,iplas-1
2179 !
2180 ! zrho=rosn(i)
2181 ! CALL E02BCF(nspl+4,RRK,CCK,zrho,0,CWk,IFAIL)
2182 ! C_driv(i)=cwk(1)
2183 !
2184 ! enddo
2185 
2186  nspl=na1+1
2187  call linsplin(nspl,funw,rhow, iplas,c_driv,rosn)
2188 
2189  !!!!!!!!!!!!!!!!!!!cd!!!!!!!!!!!!!!!!!!!!!!!
2190 
2191  !!!!!!!!!!!!!!!!!!!cubs!!!!!!!!!!!!!!!!!!!!!!!
2192 
2193  !funw(1)=cbx0
2194  funw(1)=cbx1
2195  rhow(1)=0.d0
2196 
2197  do i=2,na1+1
2198  funw(i)=cubs(i-1)
2199  rhow(i)=rho(i-1)/roc
2200  enddo
2201 
2202  c_bts(1)=funw(1)
2203  c_bts(iplas)=funw(na1+1)
2204 
2205 ! CALL E01BAF(nspl,rhow,funw,RRK,CCK,
2206 ! * nspl+4,WRK,6*nspl+16,IFAIL)
2207 !
2208 ! do i=2,iplas-1
2209 !
2210 ! zrho=rosn(i)
2211 ! CALL E02BCF(nspl+4,RRK,CCK,zrho,0,CWk,IFAIL)
2212 ! C_bts(i)=cwk(1)
2213 !
2214 ! enddo
2215 
2216  nspl=na1+1
2217  call linsplin(nspl,funw,rhow, iplas,c_bts,rosn)
2218 
2219  !!!!!!!!!!!!!!!!!!!cubs!!!!!!!!!!!!!!!!!!!!!!!
2220 
2221  !!!!!!!!!!!!!!!!!!!te!!!!!!!!!!!!!!!!!!!!!!!
2222 
2223  !funw(1)=cex0
2224  funw(1)=cex1
2225  rhow(1)=0.d0
2226 
2227  do i=2,na1+1
2228  funw(i)=te(i-1)
2229  rhow(i)=rho(i-1)/roc
2230  enddo
2231 
2232  t_el(1)=funw(1)
2233  t_el(iplas)=funw(na1+1)
2234 
2235 ! CALL E01BAF(nspl,rhow,funw,RRK,CCK,
2236 ! * nspl+4,WRK,6*nspl+16,IFAIL)
2237 !
2238 ! do i=2,iplas-1
2239 !
2240 ! zrho=rosn(i)
2241 ! CALL E02BCF(nspl+4,RRK,CCK,zrho,0,CWk,IFAIL)
2242 ! T_el(i)=cwk(1)
2243 !
2244 ! enddo
2245 
2246  nspl=na1+1
2247  call linsplin(nspl,funw,rhow, iplas,t_el,rosn)
2248 
2249  !!!!!!!!!!!!!!!!!!!te!!!!!!!!!!!!!!!!!!!!!!!
2250 
2251  do i=1,iplas
2252  c_bts(i)=c_bts(i)*btor
2253  c_driv(i)=c_driv(i)*btor
2254  t_el(i)=t_el(i)*1.d3
2255  enddo
2256 
2257  return
2258  end
2259 
2260 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2261 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2262  subroutine profro_cu_l(na1,neql,rtor,btor,cu,rho,roc)
2263 
2264  include 'double.inc'
2265  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
2266  include 'dim.inc'
2267  include 'compol.inc'
2268 
2269  common /com_jb/ bj_av(nrp),curfi_av(nrp)
2270  real*8 rosn(nrp),funw(nursp),rhow(nursp)
2271  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
2272  real*8 cwk(4)
2273 
2274  real*8 rho(*),rtor,btor,roc,cu(*)
2275 
2276  iplas=neql
2277  nspl=na1+1
2278  do i=1,iplas
2279  rosn(i)=(i-1.d0)/(iplas-1.d0)
2280  enddo
2281 
2282 !! exstrapolation on magn. axis
2283 
2284  rh0=0.d0
2285 
2286  rh1=rho(1)
2287  rh2=rho(2)
2288  rh3=rho(3)
2289 
2290  cux1=cu(1)
2291  cux2=cu(2)
2292  cux3=cu(3)
2293 
2294 
2295  call extrp2(rh0,cux0, rh1,rh2,rh3, cux1,cux2,cux3)
2296 
2297 !!!!!!!!!!!!!!!!!!!cu!!!!!!!!!!!!!!!!!!!!!!!
2298 
2299  !funw(1)=cux0
2300  funw(1)=cux1
2301  rhow(1)=0.d0
2302 
2303  do i=2,na1+1
2304  funw(i)=cu(i-1)
2305  rhow(i)=rho(i-1)/roc
2306  enddo
2307 
2308  bj_av(1)=funw(1)
2309  bj_av(iplas)=funw(na1+1)
2310 
2311 ! CALL E01BAF(nspl,rhow,funw,RRK,CCK,
2312 ! * nspl+4,WRK,6*nspl+16,IFAIL)
2313 !
2314 ! do i=2,iplas-1
2315 !
2316 ! zrho=rosn(i)
2317 ! CALL E02BCF(nspl+4,RRK,CCK,zrho,0,CWk,IFAIL)
2318 ! BJ_av(i)=cwk(1)
2319 !
2320 ! enddo
2321 
2322  nspl=na1+1
2323  call linsplin(nspl,funw,rhow, iplas,bj_av,rosn)
2324 
2325 
2326  !!!!!!!!!!!!!!!!!!!cu!!!!!!!!!!!!!!!!!!!!!!!
2327 
2328  do i=1,iplas
2329  bj_av(i)=bj_av(i)*btor*amu0
2330  enddo
2331 
2332  return
2333  end
2334 
2335 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2336 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2337  subroutine profro_cu(na1,neql,rtor,btor,cu,rho,roc)
2338 
2339  include 'double.inc'
2340  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
2341  include 'dim.inc'
2342  include 'compol.inc'
2343 
2344  common /com_jb/ bj_av(nrp),curfi_av(nrp)
2345  real*8 rosn(nrp),funw(nursp),rhow(nursp)
2346  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
2347  real*8 cwk(4)
2348 
2349  real*8 rho(*),rtor,btor,roc,cu(*)
2350 
2351  iplas=neql
2352  nspl=na1+1
2353  do i=1,iplas
2354  rosn(i)=(i-1.d0)/(iplas-1.d0)
2355  enddo
2356 
2357 !! exstrapolation on magn. axis
2358 
2359  rh0=0.d0
2360 
2361  rh1=rho(1)
2362  rh2=rho(2)
2363  rh3=rho(3)
2364 
2365  cux1=cu(1)
2366  cux2=cu(2)
2367  cux3=cu(3)
2368 
2369 
2370  call extrp2(rh0,cux0, rh1,rh2,rh3, cux1,cux2,cux3)
2371 
2372 !!!!!!!!!!!!!!!!!!!cu!!!!!!!!!!!!!!!!!!!!!!!
2373 
2374  funw(1)=cux0
2375  rhow(1)=0.d0
2376 
2377  do i=2,na1+1
2378  funw(i)=cu(i-1)
2379  rhow(i)=rho(i-1)/roc
2380  enddo
2381 
2382  bj_av(1)=funw(1)
2383  bj_av(iplas)=funw(na1+1)
2384 
2385  CALL e01baf(nspl,rhow,funw,rrk,cck,
2386  * nspl+4,wrk,6*nspl+16,ifail)
2387 
2388  do i=2,iplas-1
2389 
2390  zrho=rosn(i)
2391  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
2392  bj_av(i)=cwk(1)
2393 
2394  enddo
2395 
2396  !!!!!!!!!!!!!!!!!!!cu!!!!!!!!!!!!!!!!!!!!!!!
2397 
2398  do i=1,iplas
2399  bj_av(i)=bj_av(i)*btor*amu0
2400  enddo
2401 
2402  return
2403  end
2404 
2405 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2406 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2407  subroutine profro_pres(na1,neql,pres,rho,roc)
2408 
2409  include 'double.inc'
2410  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
2411  include 'dim.inc'
2412  include 'compol.inc'
2413  common /com_pres/ p_rho(0:nrp),dpdro(nrp),dpdfi(1:nrp),romin
2414 
2415  real*8 rocn(nrp),funw(nursp),rhow(nursp)
2416  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
2417  real*8 cwk(4)
2418 
2419  real*8 rho(*),roc,pres(*)
2420  real*8 rokn(nrp)
2421 
2422  iplas=neql
2423  nspl=na1+1
2424 
2425  do i=1,iplas-1
2426  rocn(i)=(i-0.5d0)/(iplas-1.d0)
2427  enddo
2428 
2429  rocn(iplas)=1.d0
2430 
2431  do i=1,iplas
2432  rokn(i)=((i-1.d0)/(iplas-1.d0))
2433  enddo
2434 
2435  rokn(iplas)=1.d0
2436  rokn(1)=0.d0
2437 
2438 !! exstrapolation on magn. axis
2439 
2440  rh0=0.d0
2441 
2442  rh1=rho(1)
2443  rh2=rho(2)
2444  rh3=rho(3)
2445 
2446  px1=pres(1)
2447  px2=pres(2)
2448  px3=pres(3)
2449 
2450 
2451  call extrp2(rh0,px0, rh1,rh2,rh3, px1,px2,px3)
2452 
2453 !!!!!!!!!!!!!!!!!!!p!!!!!!!!!!!!!!!!!!!!!!!
2454 
2455  funw(1)=px0
2456  !funw(1)=pres(1)
2457  rhow(1)=0.d0
2458 
2459  do i=2,na1+1
2460  funw(i)=pres(i-1)
2461  rhow(i)=rho(i-1)/roc
2462  enddo
2463 
2464  p_rho(0)=funw(1)
2465  p_rho(iplas)=funw(na1+1)
2466 
2467  CALL e01baf(nspl,rhow,funw,rrk,cck,
2468  * nspl+4,wrk,6*nspl+16,ifail)
2469 
2470  do i=1,iplas-1
2471 
2472  zrho=rocn(i)
2473  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
2474  p_rho(i)=cwk(1)
2475 
2476  enddo
2477 
2478  do i=2,na1+1
2479 
2480  zrho=rhow(i)
2481  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
2482  dpdfi(i)=cwk(2)/zrho
2483 
2484  enddo
2485  do i=1,iplas
2486 
2487  zrho=rokn(i)
2488  if(zrho.lt.rhow(2)) zrho=rhow(2)
2489  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
2490  dpdro(i)=0.5d0*cwk(2)/zrho
2491 
2492  enddo
2493 
2494  romin=rhow(2)
2495 
2496 !!!!!!!!!!!!!!!!!!!p!!!!!!!!!!!!!!!!!!!!!!!
2497 
2498 
2499 ! do i=1,iplas
2500 ! C_bts(i)=C_bts(i)*btor
2501 ! C_driv(i)=C_driv(i)*btor
2502 ! T_el(i)=T_el(i)*1.d3
2503 ! enddo
2504 
2505  return
2506  end
2507 !!!!!!!!!!!!!!!!!!!p!!!!!!!!!!!!!!!!!!!!!!!
2508  subroutine proffi_pres(na1,neql,pres,rho,roc)
2509 
2510  use spider_spline
2511 
2512  include 'double.inc'
2513  !parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
2514  include 'dim.inc'
2515  include 'compol.inc'
2516  common /com_pres/ p_rho(0:nrp),dpdro(nrp),dpdfi(1:nrp),romin
2517 
2518  real*8 rocn(nrp) !,funw(nursp),rhow(nursp)
2519  !real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
2520  !real*8 cwk(4)
2521 
2522  real*8 rho(*),roc,pres(*)
2523  real*8 rokn(nrp)
2524 
2525  real(8), allocatable :: pxin(:), pyin(:), pxout(:), pyout(:),
2526  + pyoutp(:), pyoutpp(:), py2(:), pwork(:),
2527  + pamat(:,:), pyinnew(:)
2528 
2529  real(8), allocatable :: p05_prim(:),rho05(:),d2pf(:)
2530 
2531  allocate(rho05(na1+1))
2532  allocate(p05_prim(na1+1))
2533  allocate(d2pf(neql))
2534 
2535  do i=2,na1
2536  rho05(i)=0.5d0*( rho(i-1)**2+rho(i)**2 )/rho(na1)**2
2537  enddo
2538  rho05(1)=0.d0
2539  rho05(na1+1)=1.d0
2540 
2541  do i=2,na1
2542  p05_prim(i)=( pres(i)-pres(i-1) )
2543  & /(rho(i)**2-rho(i-1)**2)*rho(na1)**2
2544  enddo
2545 
2546 !! exstrapolation on magn. axis
2547 
2548  rh0=0.d0
2549 
2550  rh1=rho05(2)
2551  rh2=rho05(3)
2552  rh3=rho05(4)
2553 
2554  px1=p05_prim(2)
2555  px2=p05_prim(3)
2556  px3=p05_prim(4)
2557 
2558  call extrp2(rh0,px0, rh1,rh2,rh3, px1,px2,px3)
2559 
2560  p05_prim(1)=px0
2561 
2562 !! exstrapolation to boundry
2563 
2564  rh0=1.d0
2565 
2566  rh1=rho05(na1)
2567  rh2=rho05(na1-1)
2568  rh3=rho05(na1-2)
2569 
2570  px1=p05_prim(na1)
2571  px2=p05_prim(na1-1)
2572  px3=p05_prim(na1-2)
2573 
2574  call extrp2(rh0,px0, rh1,rh2,rh3, px1,px2,px3)
2575 
2576  p05_prim(na1+1)=px0
2577 
2578  iplas=neql
2579  nspl=na1+1
2580 
2581  do i=1,iplas-1
2582  rocn(i)=((i-0.5d0)/(iplas-1.d0))**2
2583  enddo
2584 
2585  rocn(iplas)=1.d0
2586 
2587  do i=1,iplas
2588  rokn(i)=((i-1.d0)/(iplas-1.d0))**2
2589  enddo
2590 
2591  rokn(iplas)=1.d0
2592  rokn(1)=0.d0
2593 
2594 
2595 
2596  knin = nspl
2597  knout = iplas
2598  kopt = 1
2599  ptaus = 1d-6
2600  allocate(pyinnew(knin))
2601  allocate(py2(knin))
2602  allocate(pwork(knin))
2603  mdamat = 7
2604  pbclft = 2.d0
2605  pbcrgt = 2.d0
2606  allocate(pamat(mdamat,knin))
2607  CALL intrptau(rho05,p05_prim,pyinnew,py2,knin,rokn,dpdfi,d2pf,
2608  + pyoutpp,knout,kopt,ptaus,pwork,pamat,mdamat,pbclft,pbcrgt)
2609 
2610 
2611 !!!!!!!!!!!!!!!!!!!p!!!!!!!!!!!!!!!!!!!!!!!
2612 
2613  return
2614  end
2615 
2616 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2617  subroutine proffi_pres_(na1,neql,pres,rho,roc)
2618 
2619  include 'double.inc'
2620  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
2621  include 'dim.inc'
2622  include 'compol.inc'
2623  common /com_pres/ p_rho(0:nrp),dpdro(nrp),dpdfi(1:nrp),romin
2624 
2625  real*8 rocn(nrp),funw(nursp),rhow(nursp)
2626  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
2627  real*8 cwk(4)
2628 
2629  real*8 rho(*),roc,pres(*)
2630  real*8 rokn(nrp)
2631 
2632  iplas=neql
2633  nspl=na1+1
2634 
2635  do i=1,iplas-1
2636  rocn(i)=((i-0.5d0)/(iplas-1.d0))**2
2637  enddo
2638 
2639  rocn(iplas)=1.d0
2640 
2641  do i=1,iplas
2642  rokn(i)=((i-1.d0)/(iplas-1.d0))**2
2643  enddo
2644 
2645  rokn(iplas)=1.d0
2646  rokn(1)=0.d0
2647 
2648 !! exstrapolation on magn. axis
2649 
2650  rh0=0.d0
2651 
2652  rh1=rho(1)
2653  rh2=rho(2)
2654  rh3=rho(3)
2655 
2656  px1=pres(1)
2657  px2=pres(2)
2658  px3=pres(3)
2659 
2660 
2661  call extrp2(rh0,px0, rh1,rh2,rh3, px1,px2,px3)
2662 
2663 !!!!!!!!!!!!!!!!!!!p!!!!!!!!!!!!!!!!!!!!!!!
2664 
2665  funw(1)=px0
2666  rhow(1)=0.d0
2667 
2668  do i=2,na1+1
2669  funw(i)=pres(i-1)
2670  rhow(i)=(rho(i-1)/roc)**2
2671  enddo
2672 
2673  p_rho(0)=funw(1)
2674  p_rho(iplas)=funw(na1+1)
2675 
2676  CALL e01baf(nspl,rhow,funw,rrk,cck,
2677  * nspl+4,wrk,6*nspl+16,ifail)
2678 
2679  do i=1,iplas-1
2680 
2681  zrho=rocn(i)
2682  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
2683  p_rho(i)=cwk(1)
2684 
2685  enddo
2686 
2687  do i=1,iplas
2688 
2689  zrho=rokn(i)
2690  CALL e02bcf(nspl+4,rrk,cck,zrho,0,cwk,ifail)
2691  dpdfi(i)=cwk(2)
2692 
2693  enddo
2694 
2695 !!!!!!!!!!!!!!!!!!!p!!!!!!!!!!!!!!!!!!!!!!!
2696 
2697  return
2698  end
2699 
2700 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2701  subroutine profro_pres_l(na1,neql,pres,rho,roc)
2702 
2703  include 'double.inc'
2704  parameter(nursp=1000)
2705  include 'dim.inc'
2706  include 'compol.inc'
2707  common /com_pres/ p_rho(0:nrp),dpdro(nrp),dpdfi(1:nrp),romin
2708 
2709  real*8 rocn(nrp),funw(nursp),rhow(nursp)
2710  real*8 rho(*),roc,pres(*)
2711 
2712  iplas=neql
2713  nspl=na1+1
2714 
2715  do i=1,iplas-1
2716  rocn(i)=(i-0.5d0)/(iplas-1.d0)
2717  enddo
2718 
2719  rocn(iplas)=1.d0
2720 
2721 !! exstrapolation to magn. axis
2722 
2723  rh0=0.d0
2724 
2725  rh1=rho(1)
2726  rh2=rho(2)
2727  rh3=rho(3)
2728 
2729  px1=pres(1)
2730  px2=pres(2)
2731  px3=pres(3)
2732 
2733 
2734  call extrp2(rh0,px0, rh1,rh2,rh3, px1,px2,px3)
2735 
2736 !!!!!!!!!!!!!!!!!!!p!!!!!!!!!!!!!!!!!!!!!!!
2737 
2738  funw(1)=px0
2739  rhow(1)=0.d0
2740 
2741  do i=2,na1+1
2742  funw(i)=pres(i-1)
2743  rhow(i)=rho(i-1)/roc
2744  enddo
2745 
2746  p_rho(0)=funw(1)
2747  p_rho(iplas)=funw(na1+1)
2748 
2749  do i=1,iplas-1
2750 
2751  zrho=rocn(i)
2752  do j=1,na1
2753  if( zrho.gt.rhow(j) .AnD. zrho.le.rhow(j+1) ) then
2754  pres_rho=( funw(j)*(rhow(j+1)-zrho) +
2755  & funw(j+1)*(zrho-rhow(j)) )
2756  & /( rhow(j+1)-rhow(j) )
2757  go to 100
2758  endif
2759  enddo
2760 
2761  100 p_rho(i)=pres_rho
2762 
2763  enddo
2764 
2765 !!!!!!!!!!!!!!!!!!!p!!!!!!!!!!!!!!!!!!!!!!!
2766 
2767  return
2768  end
2769 
2770 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2771 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2772  subroutine pres_d_psi
2773 
2774  include 'double.inc'
2775  include 'dim.inc'
2776  parameter(nrsp=nrp+1,nrsp4=nrsp+4,nrsp6=nrsp4*6)
2777  include 'compol.inc'
2778  common /com_pres/ p_rho(0:nrp),dpdro(nrp),dpdfi(1:nrp),romin
2779 
2780  real*8 rhocn(nrp+1),funw(nrsp),rhokn(nrp)
2781  real*8 rrk(nrsp4),cck(nrsp4),wrk(nrsp6)
2782  real*8 cwk(4)
2783 
2784 
2785  do i=1,iplas-1
2786  rhocn(i+1)=(i-0.5d0)/(iplas-1.d0)
2787  enddo
2788 
2789  rhocn(iplas+1)=1.d0
2790  rhocn(1)=0.d0
2791 
2792  do i=1,iplas
2793  rhokn(i)=(i-1.d0)/(iplas-1.d0)
2794  enddo
2795 
2796 !! exstrapolation q to magn. axis
2797 
2798  rh0=0.d0
2799 
2800  rh1=0.5d0
2801  rh2=1.5d0
2802  rh3=2.5d0
2803 
2804  qx1=q(1)
2805  qx2=q(2)
2806  qx3=q(3)
2807 
2808  call extrp2(rh0,qx0, rh1,rh2,rh3, qx1,qx2,qx3)
2809 
2810  !dpdpsi(1)=-amu0*qx0*dpdro(1)/flucfm
2811  dpdpsi(1)=-amu0*qx0*dpdfi(1)/flucfm
2812  funw(1)=qx0
2813 
2814  do i=1,iplas
2815  funw(i+1)=q(i)
2816  enddo
2817 
2818 !!!!!!!!!!!!!!!!!!!!!!! q !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2819  n3spl=iplas+1
2820  CALL e01baf(n3spl,rhocn,funw,rrk,cck,
2821  * n3spl+4,wrk,6*n3spl+16,ifail)
2822 
2823  do i=2,iplas
2824  zrho=rhokn(i)
2825  CALL e02bcf(n3spl+4,rrk,cck,zrho,0,cwk,ifail)
2826  qkn=cwk(1)
2827  !dpdpsi(i)=-amu0*qkn*dpdro(i)/flucfm
2828  dpdpsi(i)=-amu0*qkn*dpdfi(i)/flucfm
2829  !if(zrho.lt.romin) dpdpsi(i)=dpdpsi(1)
2830  enddo
2831  !dpdpsi(1)=dpdpsi(2)
2832 
2833 
2834  return
2835  end
2836 
2837 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2838 
2839 !********************************************************************
2840 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2841  subroutine get_rz(rout,zout,ni_p,nj_p)
2842 
2843  include 'double.inc'
2844  include 'dim.inc'
2845  include 'compol.inc'
2846 
2847  real*8 rout(ni_p,nj_p),zout(ni_p,nj_p)
2848 
2849  do i=1,iplas
2850  do j=1,nt
2851  rout(i,j)=r(i,j)
2852  zout(i,j)=z(i,j)
2853  enddo
2854  enddo
2855 
2856  return
2857  end
2858 
2859 
2860  subroutine newgrd_ss( ROC, HRO, NB1, NA1, RHO,btor, HROA )
2861 c----------------------------------------------------------------------|
2862 c define size of the edge cell to be: .6 <= hroa/hro < 1.8
2863 c with a hysteresis of 0.2*hro, so that
2864 c hroa<=0.6*hro jumps to hroa<=1.6*hro
2865 c hroa>=1.8*hro jumps to hroa>=0.8*hro
2866 c write(*,*)na1,time,rho(na1)," -> ",roc
2867 c----------------------------------------------------------------------|
2868  include 'double.inc'
2869  include 'dim.inc'
2870  include 'compol.inc'
2871  dimension rho(*)
2872 
2873  na=na1-1
2874  roc=sqrt(flx_fi(iplas)/(pi*btor))
2875 
2876 c write(*,*)na1,roc,hro*(na+0.1),hro*(na1+0.3),hro*(nb1-0.49)
2877 c write(*,'(1P,6E13.5)')(rho(j),j=na-2,na1+1)
2878 c write(*,'(1P,6E13.5)')(ametr(j),j=na-2,na1+1)
2879  if (roc .gt. hro*(nb1-0.49)) then
2880 c rho_edge gets out of the grid
2881  write(*,*)'>>> WARNING: the allocated grid is too small'
2882  write(*,*)'Time =',time,' Rho_b =',roc,
2883  > ' Rho_max = ',hro*(nb1-0.5)
2884  write(*,*)'Try to increase AWALL'
2885  write(*,*)'If this does not help inspect equilibrium input'
2886  na = nb1-1
2887  elseif (roc .le. hro*(na+0.1)) then ! 0.1 <=> 0.6-0.5
2888 ! if (HROA < 0.6*HRO) then reduce NA
2889 c write(*,*)" <- ",roc,hro*na,rho(na1),hro*na1,na1
2890  do j=na-1,1,-1
2891  na = j
2892 c write(*,*)j,na*hro,roc,(na+1)*hro,roc/hro-(na-0.5)
2893  if (roc .gt. hro*(j+0.1)) goto 10
2894  enddo
2895  elseif (roc .gt. hro*(na1+0.3)) then ! 0.3 <=> 1.8-1.5
2896 ! if (HROA > 1.8*HRO) then increase NA
2897 c write(*,*)" -> ",roc,hro*na,rho(na1),hro*na1,na1
2898  do j=na1,nb1
2899  na = j
2900 c write(*,*)j,na*hro,roc,(na+1)*hro
2901  if (roc .le. hro*(j+1.3)) goto 10
2902  enddo
2903  endif
2904  10 continue
2905 
2906 c hro*(na+0.6) < roc <= hro*(na+1.8)
2907 c write(*,*)ab,abc,time,tstart
2908 c write(*,*)na+1,na+1-na1,hro*(na+0.6),roc,hro*(na+1.8)
2909 
2910  na1 = na+1
2911  do j=1,nb1
2912  rho(j)=(j-0.5)*hro
2913  enddo
2914  hroa = roc-rho(na)
2915  rho(na1) = roc
2916 
2917  return
2918  end
2919 c======================================================================|
2920 
2921  subroutine get_roc(ROC,btor)
2922 
2923 c implicit none
2924  include 'double.inc'
2925  include 'dim.inc'
2926  include 'compol.inc'
2927 
2928  real*8 roc,btor
2929  roc=sqrt(flx_fi(iplas)/(pi*btor))
2930 
2931  return
2932  end
2933 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2934  subroutine get_phi(phi)
2935 
2936 c implicit none
2937  include 'double.inc'
2938  include 'dim.inc'
2939  include 'compol.inc'
2940 
2941  real*8 phi(*)
2942  integer i
2943  do i=1,iplas
2944  phi(i)=flx_fi(i)
2945  enddo
2946 
2947  return
2948  end
2949 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2950 
2951  subroutine pla_volt(dt)
2952  include 'double.inc'
2953  include 'dim.inc'
2954  include 'compol.inc'
2955  common /com_volt/ upls(nrp)
2956  common /com_psisave/ psi_0(nrp)
2957 
2958  call get_psibon(psi_bn1)
2959 
2960  do i=1,iplas
2961  !upls(i)=(psia(i)*psim-psi_0(i))/dt
2962  upls(i)=(psi(i,2)+psi_bn1-psi_0(i))/dt
2963  enddo
2964 
2965  return
2966  end
2967 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2968 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2969 
2970  subroutine savepsi
2971  include 'double.inc'
2972  include 'dim.inc'
2973  include 'compol.inc'
2974  common /com_psisave/ psi_0(nrp)
2975 
2976  call get_psibon(psi_bn1)
2977 
2978  do i=1,iplas
2979  !psi_0(i)=psim*psia(i)
2980  psi_0(i)=psi(i,2)+psi_bn1
2981  enddo
2982 
2983  return
2984  end
2985 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2986 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2987 
2988  subroutine get_psibon(psi_bnd)
2989  include 'double.inc'
2990  include 'dim.inc'
2991  common/selcon/ psi_d(nrp),fi_d(nrp),f_d(nrp),ri_d(nrp),
2992  * ps_pnt(nrp),del_psb,psi_bn1
2993 
2994  psi_bnd=psi_bn1
2995 
2996  return
2997  end
2998 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2999  subroutine flux_state(flx_st)
3000  include 'double.inc'
3001  include 'param.inc'
3002  include 'comblc.inc'
3003 
3004 
3005  sj=0.d0
3006  spsj=0.d0
3007 
3008  do j=1,nj
3009  do i=1,ni
3010 
3011  if(ipr(i,j).eq.1) then
3012  sj=sj+curf(i,j)
3013  spsj=spsj+curf(i,j)*ue(i,j)
3014  endif
3015 
3016  enddo
3017  enddo
3018 
3019  flx_st=spsj/sj
3020 
3021  return
3022  end
3023 
3024  subroutine linsplin(nspl,usp,xsp, n,u,x)
3025  include 'double.inc'
3026  real*8 usp(nspl),xsp(nspl)
3027  real*8 u(n),x(n)
3028 
3029  u(1)=usp(1)
3030  u(n)=usp(nspl)
3031 
3032 
3033  do i=2,n-1
3034  zx=x(i)
3035 
3036  do ic=1,nspl-1
3037  if(zx.le.xsp(ic+1) .AnD. zx.gt.xsp(ic)) then
3038  u(i)=( usp(ic)*(xsp(ic+1)-zx)
3039  * + usp(ic+1)*(zx-xsp(ic)) )
3040  * /(xsp(ic+1)-xsp(ic))
3041  exit
3042  endif
3043  enddo
3044 
3045  enddo
3046 
3047  return
3048  end
3049 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3050 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3051  subroutine prof_ast_sp_l(na1,neql,ArrAS,ArrSP,rho,roc)
3052 
3053 !!!!!ASTRA - SPIDER arrays transformation
3054 !!!!!ArrAS - ASTRA array(input)
3055 !!!!!ArrSP - SPIDER array(output)
3056  include 'double.inc'
3057  include 'dim.inc'
3058  include 'compol.inc'
3059 
3060  parameter(nursp=1000,nursp4=nursp+4,nursp6=nursp4*6)
3061 
3062  real*8 rosn(nrp),funw(nursp),rhow(nursp)
3063  real*8 rrk(nursp4),cck(nursp4),wrk(nursp6)
3064  real*8 cwk(4)
3065 
3066  real*8 rho(*),arras(*),arrsp(*),roc
3067 
3068  iplas=neql
3069  nspl=na1+1
3070  do i=1,iplas
3071  rosn(i)=(i-1.d0)/(iplas-1.d0)
3072  enddo
3073 
3074 !! exstrapolation on magn. axis
3075 
3076  rh0=0.d0
3077 
3078  rh1=rho(1)
3079  rh2=rho(2)
3080  rh3=rho(3)
3081 
3082  cux1=arras(1)
3083  cux2=arras(2)
3084  cux3=arras(3)
3085 
3086 
3087  call extrp2(rh0,cux0, rh1,rh2,rh3, cux1,cux2,cux3)
3088 
3089 !!!!!!!!!!!!!!!!!!!cu!!!!!!!!!!!!!!!!!!!!!!!
3090 
3091  funw(1)=cux0
3092  !funw(1)=cux1
3093  rhow(1)=0.d0
3094 
3095  do i=2,na1+1
3096  funw(i)=arras(i-1)
3097  rhow(i)=rho(i-1)/roc
3098  enddo
3099 
3100  arrsp(1)=funw(1)
3101  arrsp(iplas)=funw(na1+1)
3102 
3103 ! CALL E01BAF(nspl,rhow,funw,RRK,CCK,
3104 ! * nspl+4,WRK,6*nspl+16,IFAIL)
3105 !
3106 ! do i=2,iplas-1
3107 !
3108 ! zrho=rosn(i)
3109 ! CALL E02BCF(nspl+4,RRK,CCK,zrho,0,CWk,IFAIL)
3110 ! BJ_av(i)=cwk(1)
3111 !
3112 ! enddo
3113 
3114  nspl=na1+1
3115  call linsplin(nspl,funw,rhow, iplas,arrsp,rosn)
3116 
3117 
3118  return
3119  end
3120 
3121 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
proffi_pres
subroutine proffi_pres(na1, neql, pres, rho, roc)
Definition: interfac.f:2508
profro_p_l
subroutine profro_p_l(na1, neql, rtor, eqpf, rho, roc)
Definition: interfac.f:1808
profro_p
subroutine profro_p(na1, neql, rtor, eqpf, rho, roc)
Definition: interfac.f:1720
e02bcf
subroutine e02bcf(NCAP7, K, C, X, LEFT, S, IFAIL)
Definition: NAG.f:2761
proffi_pres_
subroutine proffi_pres_(na1, neql, pres, rho, roc)
Definition: interfac.f:2617
profro
subroutine profro(na1, neql, rtor, eqpf, eqff, rho, roc)
Definition: interfac.f:1630
put_ipl
subroutine put_ipl(placur)
Definition: Spider_call.f:243
savepsi
subroutine savepsi
Definition: interfac.f:2970
phys_functions::bp2
real(r8) function bp2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:360
get_phi
subroutine get_phi(phi)
Definition: interfac.f:2934
tab_efit
subroutine tab_efit(tokf, psax, eqdfn, rax, zax, b0, r0)
Definition: B_wrd.f:447
newgrd_ss
subroutine newgrd_ss(ROC, HRO, NB1, NA1, RHO, btor, HROA)
Definition: interfac.f:2860
profro_sbc_l
subroutine profro_sbc_l(na1, neql, rtor, btor, cc, Te, cubs, cd, rho, roc)
Definition: interfac.f:2082
out
subroutine out(rbnd, zbnd, zli3, betpol, bettot, parpla)
Definition: _wrd.f:208
get_roc
subroutine get_roc(ROC, btor)
Definition: interfac.f:2921
helena21_mod::dpdpsi
real(r8) function dpdpsi(psi_n)
Definition: helena21_mod.f90:122
bnd_modul
Definition: bnd_modul.f90:1
e01baf
subroutine e01baf(M, X, Y, K, C, LCK, WRK, LWRK, IFAIL)
Definition: NAG.f:2511
extrp2
subroutine extrp2(X, F, X1, X2, X3, F1, F2, F3)
Definition: com_sub.f:719
grid
subroutine grid
Definition: EQ1_m.f:926
spider_spline
Definition: spider_spline.f:1
b_extrp
subroutine b_extrp(x0, X1, X2, X3, j1, j2, j3, j, yarr, jerr)
Definition: b_extrp.f:1
flux_state
subroutine flux_state(flx_st)
Definition: interfac.f:2999
profro_pres_l
subroutine profro_pres_l(na1, neql, pres, rho, roc)
Definition: interfac.f:2701
pres_d_psi
subroutine pres_d_psi
Definition: interfac.f:2772
profro_pres
subroutine profro_pres(na1, neql, pres, rho, roc)
Definition: interfac.f:2407
astra2spider
subroutine astra2spider(neql, nteta, nbnd, rzbnd, key_dmf,
Definition: interfac.f:1
profro_cu
subroutine profro_cu(na1, neql, rtor, btor, cu, rho, roc)
Definition: interfac.f:2337
profro_sbc
subroutine profro_sbc(na1, neql, rtor, btor, cc, Te, cubs, cd, rho, roc)
Definition: interfac.f:1919
linsplin
subroutine linsplin(nspl, usp, xsp, n, u, x)
Definition: interfac.f:3024
comet
subroutine comet(na1, platok,
Definition: interfac.f:427
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
comet_
subroutine comet_(na1, platok,
Definition: interfac.f:1111
get_rz
subroutine get_rz(rout, zout, ni_p, nj_p)
Definition: interfac.f:2841
get_psibon
subroutine get_psibon(psi_bnd)
Definition: interfac.f:2988
pla_volt
subroutine pla_volt(dt)
Definition: interfac.f:2951
profro_cu_l
subroutine profro_cu_l(na1, neql, rtor, btor, cu, rho, roc)
Definition: interfac.f:2262
prof_ast_sp_l
subroutine prof_ast_sp_l(na1, neql, ArrAS, ArrSP, rho, roc)
Definition: interfac.f:3051
spider2astra
subroutine spider2astra(rout, zout, rtor, btor, rho, roc, na1,
Definition: interfac.f:369
spider_spline::intrptau
subroutine intrptau(PXIN, PYIN, PYINNEW, PY2, KNIN, PXOUT, PYOUT, PYOUTP, PYOUTPP, KNOUT, KOPT, PTAUS, PWORK, PAMAT, MDAMAT, PBCLFT, PBCRGT)
Definition: spider_spline.f:5