ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
curfit_f.f
Go to the documentation of this file.
1  subroutine prefit(rk,zk,ncpfc,NECON,WECON,rax,zax,alp_b,psi_bnd)
2 
3  use bnd_modul
4 
5  include 'double.inc'
6  include 'parcur.inc'
7  include 'prm.inc'
8  include 'comevl.inc'
9  common/com_flag/kastr
10  common/com_snf/ksnf
11  common/com_xwx/ kxwx
12 
13  real*8 rk(*),zk(*)
14  dimension wecon(*)
15  integer necon(*),ncpfc
16 
17 ! if(.not.allocated(rbtab))
18 ! %allocate( rbtab(nbtab), zbtab(nbtab) )
19 
20 ! ksnf=0 !if ksnf=1,then kxwx must be =1
21 
22 c -----------------------------
23  alph=alp_b
24  r_ax=rax
25  z_ax=zax
26 
27 c -----------------------------
28  if(kastr.eq.0) then
29 
30  write(fname,'(a,a)') path(1:kname),'bonfit.dat'
31  open(1,file=fname,form='formatted')
32  !open(1,file='bonfit.dat')
33 
34  read(1,*) wwl,wdk,wsig
35 
36  read(1,*) (d_wght(ik),ik=1,nequi)
37 
38  read(1,*) lfit
39  if(lfit.ne.0) then
40  do l=1,lfit
41  read(1,*) rfit(l),zfit(l)
42  enddo
43  endif
44  read(1,*) lpre
45  if(lpre.ne.0) then
46  do l=1,lpre
47  read(1,*) rpre(l),zpre(l)
48  enddo
49  endif
50  read(1,*) rx_p,zx_p
51  read(1,*) psi_bnd
52  if(kxwx.eq.2) then
53  read(1,*) rx2_p,zx2_p
54  endif
55  close(1)
56 
57 !!!!!!! C'*C
58  if(ksnf.eq.1) then
59  write(fname,'(a,a)') path(1:kname),'matC.wr'
60  open(1,file=fname,form='formatted')
61  read(1,*) c_wght
62  read(1,*) n_cc,m_cc
63  read(1,*) (g_wght(i),i=1,n_cc)
64  read(1,*) ((cpir(i,j),i=1,n_cc),j=1,m_cc)
65  close(1)
66 
67  do i=1,n_cc
68  do j=1,m_cc
69  cpir(i,j)=cpir(i,j)*g_wght(i)
70  enddo
71  enddo
72 
73  do i=1,m_cc
74  do j=1,m_cc
75  cxcs=0.d0
76  do k=1,n_cc
77  cxcs=cxcs+cpir(k,i)*cpir(k,j)
78  enddo
79  cxc(i,j)=cxcs
80  enddo
81  enddo
82 
83  endif
84 
85 !temporary
86 ! do l=1,Lfit
87 ! rfit(l)=rfit(l)+0.02d0
88 ! enddo
89 ! do l=1,Lpre
90 ! rpre(l)=rpre(l)+0.02d0
91 ! enddo
92 ! rx_p=rx_p+0.02d0
93 !temporary
94 
95  elseif(kastr.eq.1) then
96 
97  wwl= 1.d0
98  wdk= 1.d0
99  wsig=1.d-2
100 
101  do ik=1,nequi
102  d_wght(ik)=1.d0
103  enddo
104 
105 ! write(fname,'(a,a)') path(1:kname),'tab_bnd.dat'
106 ! open(1,file=fname,form='formatted')
107 ! !open(1,file='tab_bnd.dat')
108 !
109 ! read(1,*) Lfit
110  lfit=nbnd
111 
112  do l=1,lfit
113 ! read(1,*) rfit(l),zfit(l)
114  rfit(l)=rbtab(l)
115  zfit(l)=zbtab(l)
116  enddo
117 
118  lpre=0
119 
120  rx_p=rax !dummy
121  zx_p=zax !dummy
122 
123  endif
124 
125  do l=1,lfit
126  w_wght(l)=wwl*2.d0
127  enddo
128  s_wght=wsig
129 
130  call grimat(rk,zk,ncpfc,necon, wecon)
131 
132  if(kxwx.eq.2)then
133  call grimat2x(rk,zk,ncpfc,necon, wecon)
134  endif
135 ! do iq=1,NEQUI
136 ! curref(iq)=PFCEQW(iq)
137 ! enddo
138 
139 ! deallocate( rbtab, zbtab )
140  return
141  end
142 
143 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
144 
145  subroutine curfit(rk,zk,nk,NECON,WECON,psi_bnd)
146 
147  include 'double.inc'
148  include 'parcur.inc'
149  include 'prm.inc'
150  include 'comevl.inc'
151  common/com_snf/ksnf
152  common/com_xwx/ kxwx
153 
154 c -----------------------------
155  real*8 rk(*),zk(*)
156  dimension wecon(*)
157  integer necon(*),nk
158 
159 
160  if(ksnf.eq.1) then
161  call psi_cpn
162  call cursol_snf(nequi,pfceqw,psi_bnd)
163  else
164  if(kxwx.eq.1)then
165  call psi_cpn
166  call cursol(nequi,pfceqw,psi_bnd)
167  elseif(kxwx.eq.2)then
168  call psi_cpn
169  call psi_cpn_2x
170  call cursol_2x(nequi,pfceqw,psi_bnd)
171  elseif(kxwx.eq.0)then
172  call psi_cpn_wx
173  call cursol_wx(nequi,pfceqw,psi_bnd)
174  endif
175  endif
176  write(*,*) 'coil currents'
177 
178  do ik=1,nequi
179  write(*,*) pfceqw(ik),ik
180  enddo
181 
182  return
183  end
184 
185 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
186  subroutine curfit_l(rk,zk,nk,NECON,WECON,psi_bnd)
187 
188  include 'double.inc'
189  include 'parcur.inc'
190  include 'prm.inc'
191  include 'comevl.inc'
192 
193 
194 c -----------------------------
195  real*8 rk(*),zk(*)
196  dimension wecon(*)
197  integer necon(*),nk
198 
199  call psi_cpn
200 
201  call cursol_(nequi,pfceqw,psi_bnd)
202 
203  write(*,*) 'coil currents'
204 
205  do ik=1,nequi
206  write(*,*) pfceqw(ik),ik
207  enddo
208 
209  return
210  end
211 
212 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
213 
214  subroutine curfit_(rk,zk,nk,NECON,WECON,psi_bnd)
215 
216  include 'double.inc'
217  include 'dim.inc'
218  include 'parcur.inc'
219  include 'prm.inc'
220  include 'comevl.inc'
221 
222 c -----------------------------
223 
224  real*8 rk(*),zk(*)
225  dimension wecon(*)
226  integer necon(*),nk
227 
228  call bonpsi
229 
230  call cursol_(nequi,pfceqw,psi_bnd)
231 
232  write(*,*) 'coil currents'
233 
234  do ik=1,nequi
235  write(*,*) pfceqw(ik),ik
236  enddo
237 
238  return
239  end
240 
241 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
242 
243  subroutine grimat(rk,zk,ncpfc,NECON, WECON )
244 
245  include 'double.inc'
246  include 'parcur.inc'
247  include 'prm.inc'
248  include 'comevl.inc'
249 
250  real*8 rk(*),zk(*), wecon(*)
251  integer nk
252  integer necon(*)
253 
254  pi = 3.14159265359d0
255 
256  !ncpfc=nloc(npfc)
257 
258  do 1010 iq=1,nequi
259 
260  do 1876 l=1,lfit
261  gindk(iq,l)=0.d0
262  1876 continue
263 
264  do l=1,lpre
265  gindp(iq,l)=0.d0
266  enddo
267 
268  gx_r(iq)=0.d0
269  gx_z(iq)=0.d0
270 
271  gx_rz(iq)=0.d0
272  gx_zz(iq)=0.d0
273  gx_rr(iq)=0.d0
274 
275  if(lfit.ne.0) then
276  do 110 l=1,lfit
277 
278  do 112 ik=1,ncpfc
279  if( necon(ik) .eq. iq ) then
280  zgindk=greeni(rfit(l),zfit(l),rk(ik),zk(ik))/pi
281  gindk(iq,l)=gindk(iq,l)+zgindk*wecon(ik)
282  endif
283  112 continue
284 
285  110 continue
286  endif
287 
288  if(lpre.ne.0) then
289  do l=1,lpre
290 
291  do ik=1,ncpfc
292  if( necon(ik) .eq. iq ) then
293  zgindp=greeni(rpre(l),zpre(l),rk(ik),zk(ik))/pi
294  gindp(iq,l)=gindp(iq,l)+zgindp*wecon(ik)
295  endif
296  enddo
297 
298  enddo
299  endif
300 
301  do ik=1,ncpfc
302  if( necon(ik) .eq. iq ) then
303 
304  call grgren(rk(ik),zk(ik),rx_p,zx_p, dgdr,dgdz)
305  gx_r(iq)=gx_r(iq)+dgdr*wecon(ik)/pi
306  gx_z(iq)=gx_z(iq)+dgdz*wecon(ik)/pi
307 
308 !!for snowflake
309  call d2gren(rk(ik),zk(ik),rx_p,zx_p, d2gdrz,d2gdzz,d2gdrr)
310  gx_rz(iq)=gx_rz(iq)+d2gdrz*wecon(ik)/pi
311  gx_zz(iq)=gx_zz(iq)+d2gdzz*wecon(ik)/pi
312  gx_rr(iq)=gx_rr(iq)+d2gdrr*wecon(ik)/pi
313 !!for snowflake
314 
315  endif
316  enddo
317 
318  1010 continue
319 
320  return
321  end
322 
323 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
324 
325  subroutine grimat2x(rk,zk,ncpfc,NECON, WECON )
326 
327  include 'double.inc'
328  include 'parcur.inc'
329  include 'prm.inc'
330  include 'comevl.inc'
331 
332  real*8 rk(*),zk(*), wecon(*)
333  integer nk
334  integer necon(*)
335 
336  pi = 3.14159265359d0
337 
338  do 1010 iq=1,nequi
339 
340  gx2_r(iq)=0.d0
341  gx2_z(iq)=0.d0
342 
343  do ik=1,ncpfc
344  if( necon(ik) .eq. iq ) then
345 
346  call grgren(rk(ik),zk(ik),rx2_p,zx2_p, dgdr,dgdz)
347  gx2_r(iq)=gx2_r(iq)+dgdr*wecon(ik)/pi
348  gx2_z(iq)=gx2_z(iq)+dgdz*wecon(ik)/pi
349 
350 !!for snowflake
351 ! call d2GREN(rk(ik),zk(ik),Rx_p,Zx_p, d2Gdrz,d2Gdzz,d2Gdrr)
352 ! Gx_rz(iq)=Gx_rz(iq)+d2Gdrz*wecon(ik)/pi
353 ! Gx_zz(iq)=Gx_zz(iq)+d2Gdzz*wecon(ik)/pi
354 ! Gx_rr(iq)=Gx_rr(iq)+d2Gdrr*wecon(ik)/pi
355 !!for snowflake
356 
357  endif
358  enddo
359 
360  1010 continue
361 
362  return
363  end
364 
365 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
366  subroutine precal(rk,zk,nk,NECON, WECON )
367 
368  include 'double.inc'
369  include 'parcur.inc'
370  include 'prm.inc'
371  include 'comevl.inc'
372  parameter(nshp=10)
373  include 'param_.inc'
374  include 'comblc.inc'
375  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
376 
377  real*8 rk(*),zk(*), wecon(*)
378  integer nk
379  integer necon(*)
380 
381  !ncpfc=nloc(npfc)
382  ncpfc=nk
383 
384  do 1010 iq=1,nequi
385 
386  ginda(iq)=0.d0
387  gindx(iq)=0.d0
388 
389  do ik=1,ncpfc
390  if( necon(ik) .eq. iq ) then
391  zginda=greeni(r_ax,z_ax,rk(ik),zk(ik))/pi
392  ginda(iq)=ginda(iq)+zginda*wecon(ik)
393  zgindx=greeni(rx_p,zx_p,rk(ik),zk(ik))/pi
394  gindx(iq)=gindx(iq)+zgindx*wecon(ik)
395  endif
396  enddo
397 
398  1010 continue
399 
400  r0=r_ax
401  z0=z_ax
402 
403  ic=(r0-rmin)/dr(1)+1
404  jc=(z0-zmin)/dz(1)+1
405 
406 c---definition of the nearest knote
407 
408  sdmin=rmax
409 
410  do 320 k=0,1
411  rr=r(ic+k)
412  do 325 l=0,1
413  zz=z(jc+l)
414  dlx=dsqrt( (rr-r_ax)**2+(zz-z_ax)**2 )
415  if(dlx.lt.sdmin) then
416  sdmin=dlx
417  ik=ic+k
418  jk=jc+l
419  endif
420  325 continue
421  320 continue
422 
423  nsh=1
424  xs(nsh)=r(ik)
425  ys(nsh)=z(jk)
426  fun(nsh)=ui(ik,jk)
427 
428  do 400 k=-1,1
429 
430  i= ik+k
431 
432  do 410 l=-1,1
433 
434  j= jk+l
435 
436  if(k.eq.0 .AND. l.eq.0 ) go to 410
437  nsh=nsh+1
438  xs(nsh)=r(i)
439  ys(nsh)=z(j)
440  fun(nsh)=ui(i,j)
441 
442  410 continue
443  400 continue
444 
445  rra=xs(1)
446  zza=ys(1)
447 
448  call deriv5(xs,ys,fun,nsh,5,dp)
449 
450  psip_a=fun(1)+ dp(1)*(r_ax-rra) + dp(2)*(z_ax-zza)
451  + + 0.5d0*dp(3)*(r_ax-rra)*(r_ax-rra)
452  + + dp(4)*(r_ax-rra)*(z_ax-zza)
453  + + 0.5d0*dp(5)*(z_ax-zza)*(z_ax-zza)
454 
455  return
456  end
457 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
458  subroutine precal_wx(rk,zk,nk,NECON, WECON )
459 
460  include 'double.inc'
461  include 'parcur.inc'
462  include 'prm.inc'
463  include 'comevl.inc'
464  parameter(nshp=10)
465  include 'param_.inc'
466  include 'comblc.inc'
467  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
468 
469  real*8 rk(*),zk(*), wecon(*)
470  integer nk
471  integer necon(*)
472 
473  !ncpfc=nloc(npfc)
474  ncpfc=nk
475 
476  do 1010 iq=1,nequi
477 
478  ginda(iq)=0.d0
479  gindx(iq)=0.d0
480 
481  do ik=1,ncpfc
482  if( necon(ik) .eq. iq ) then
483  zginda=greeni(r_ax,z_ax,rk(ik),zk(ik))/pi
484  ginda(iq)=ginda(iq)+zginda*wecon(ik)
485  zgindx=greeni(rx0,zx0,rk(ik),zk(ik))/pi
486  gindx(iq)=gindx(iq)+zgindx*wecon(ik)
487  endif
488  enddo
489 
490  1010 continue
491 
492  r0=r_ax
493  z0=z_ax
494 
495  ic=(r0-rmin)/dr(1)+1
496  jc=(z0-zmin)/dz(1)+1
497 
498 c---definition of the nearest knote
499 
500  sdmin=rmax
501 
502  do 320 k=0,1
503  rr=r(ic+k)
504  do 325 l=0,1
505  zz=z(jc+l)
506  dlx=dsqrt( (rr-r_ax)**2+(zz-z_ax)**2 )
507  if(dlx.lt.sdmin) then
508  sdmin=dlx
509  ik=ic+k
510  jk=jc+l
511  endif
512  325 continue
513  320 continue
514 
515  nsh=1
516  xs(nsh)=r(ik)
517  ys(nsh)=z(jk)
518  fun(nsh)=ui(ik,jk)
519 
520  do 400 k=-1,1
521 
522  i= ik+k
523 
524  do 410 l=-1,1
525 
526  j= jk+l
527 
528  if(k.eq.0 .AND. l.eq.0 ) go to 410
529  nsh=nsh+1
530  xs(nsh)=r(i)
531  ys(nsh)=z(j)
532  fun(nsh)=ui(i,j)
533 
534  410 continue
535  400 continue
536 
537  rra=xs(1)
538  zza=ys(1)
539 
540  call deriv5(xs,ys,fun,nsh,5,dp)
541 
542  psip_a=fun(1)+ dp(1)*(r_ax-rra) + dp(2)*(z_ax-zza)
543  + + 0.5d0*dp(3)*(r_ax-rra)*(r_ax-rra)
544  + + dp(4)*(r_ax-rra)*(z_ax-zza)
545  + + 0.5d0*dp(5)*(z_ax-zza)*(z_ax-zza)
546 
547  return
548  end
549 
550 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
551 
552  subroutine cursol(NEQUI,PFCEQW,psi_bnd)
553 
554  include 'double.inc'
555  include 'parcur.inc'
556 
557  dimension a(ncf_p,ncf_p),x(ncf_p),y(ncf_p),ip(ncf_p)
558  dimension psictr(ncf_p)
559  dimension pfceqw(*)
560 
561 ! a(j,k)*I(k)+a(j,NEQUI+l)*Lam(l)+a(j,NEQUI+Lpre+1)*Psi_b=y(j)
562 
563  pi=3.14159265359d0
564  amu0=0.4d0*pi
565  !amu0=1.0d0
566 
567 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
568 ! equation (.)*d(I_j)=0
569 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
570 
571 ! a(j,k)*I(k) ,j-number of equation
572 
573  do k=1,nequi
574  do j=1,nequi
575  asum=0.d0
576  do l=1,lfit
577  asum=asum+gindk(k,l)*gindk(j,l)*w_wght(l)
578  enddo
579  a(j,k)=asum
580  enddo
581  enddo
582 
583  do j=1,nequi
584  a(j,j)=a(j,j)+d_wght(j)*s_wght
585  enddo
586 
587 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
588 
589  do j=1,nequi
590  do l=1,lpre
591  a(j,nequi+l)=gindp(j,l)
592  enddo
593  enddo
594 
595 ! a(j,NEQUI+Lpre+1)*Lam_r(l) ,j-number of equation
596 
597  do j=1,nequi
598  a(j,nequi+lpre+1)=gx_r(j)
599  enddo
600 
601 ! a(j,NEQUI+Lpre+2)*Lam_z(l) ,j-number of equation
602 
603  do j=1,nequi
604  a(j,nequi+lpre+2)=gx_z(j)
605  enddo
606 
607 
608 
609 
610 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
611 
612  do j=1,nequi
613  a(j,nequi+lpre+3)=(1.d0-alph)*ginda(j)+alph*gindx(j)
614  enddo
615 
616 
617 
618 
619 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
620 
621  do j=1,nequi
622  asum=0.d0
623  do l=1,lfit
624  asum=asum+w_wght(l)*gindk(j,l)
625  enddo
626  a(j,nequi+lpre+4)=-asum
627  enddo
628 
629 ! right hand side
630 
631  do j=1,nequi
632  asum=0.d0
633  do l=1,lfit
634  asum=asum+w_wght(l)*gindk(j,l)*psifit(l)
635  enddo
636  y(j)=-asum+d_wght(j)*s_wght*curref(j)*amu0
637  enddo
638 
639 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
640 ! equation (.)*d(Lam_l)=0
641 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
642 
643 ! a(j,k)*I(k) ,j-number of equation
644 
645  do l=1,lpre
646  j=nequi+l
647  do k=1,nequi
648  a(j,k)=gindp(k,l)
649  enddo
650  enddo
651 
652 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
653 
654  do l=1,lpre
655  j=nequi+l
656  do ll=1,lpre
657  k=nequi+ll
658  a(j,k)=0.d0
659  enddo
660  enddo
661 
662 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
663 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
664 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
665 
666  do l=1,lpre
667  j=nequi+l
668  k=nequi+lpre+1
669  a(j,k)=0.d0
670  k=nequi+lpre+2
671  a(j,k)=0.d0
672  k=nequi+lpre+3
673  a(j,k)=0.d0
674  enddo
675 
676 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
677 
678  do l=1,lpre
679  j=nequi+l
680  k=nequi+lpre+4
681  a(j,k)=-1.d0
682  enddo
683 ! right hand side
684  do l=1,lpre
685  j=nequi+l
686  y(j)=-psipre(l)
687  enddo
688 
689 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
690 ! equation (.)*d(Lam_r)=0
691 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
692 
693 ! a(j,k)*I(k) ,j-number of equation
694 
695  j=nequi+lpre+1
696 
697  do k=1,nequi
698  a(j,k)=gx_r(k)
699  enddo
700 
701 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
702 
703  do ll=1,lpre
704  k=nequi+ll
705  a(j,k)=0.d0
706  enddo
707 
708 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
709 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
710 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
711 
712  k=nequi+lpre+1
713  a(j,k)=0.d0
714  k=nequi+lpre+2
715  a(j,k)=0.d0
716  k=nequi+lpre+3
717  a(j,k)=0.d0
718 
719 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
720 
721  k=nequi+lpre+4
722  a(j,k)=0.d0
723 ! right hand side
724  y(j)=-dpsx_r
725 
726 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
727 ! equation (.)*d(Lam_z)=0
728 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
729 
730 ! a(j,k)*I(k) ,j-number of equation
731 
732  j=nequi+lpre+2
733 
734  do k=1,nequi
735  a(j,k)=gx_z(k)
736  enddo
737 
738 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
739 
740  do ll=1,lpre
741  k=nequi+ll
742  a(j,k)=0.d0
743  enddo
744 
745 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
746 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
747 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
748 
749  k=nequi+lpre+1
750  a(j,k)=0.d0
751  k=nequi+lpre+2
752  a(j,k)=0.d0
753  k=nequi+lpre+3
754  a(j,k)=0.d0
755 
756 ! a(j,NEQUI+Lpre+3)*Psi_b ,j-number of equation
757 
758  k=nequi+lpre+4
759  a(j,k)=0.d0
760 ! right hand side
761  y(j)=-dpsx_z
762 
763 
764 
765 
766 
767 
768 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
769 ! equation (.)*d(Lamx)=0
770 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
771 
772 ! a(j,k)*I(k) ,j-number of equation
773 
774  j=nequi+lpre+3
775 
776  do k=1,nequi
777  a(j,k)=(1.d0-alph)*ginda(k)+alph*gindx(k)
778  enddo
779 
780 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
781 
782  do ll=1,lpre
783  k=nequi+ll
784  a(j,k)=0.d0
785  enddo
786 
787 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
788 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
789 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
790 
791  k=nequi+lpre+1
792  a(j,k)=0.d0
793  k=nequi+lpre+2
794  a(j,k)=0.d0
795  k=nequi+lpre+3
796  a(j,k)=0.d0
797 
798 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
799 
800  k=nequi+lpre+4
801  a(j,k)=-1.d0
802 
803 ! right hand side
804  y(j)=-(1.d0-alph)*psip_a-alph*psip_x
805 
806 
807 
808 
809 
810 
811 
812 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
813 ! equation (.)*d(-psi_b)=0
814 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
815 
816 ! a(j,k)*I(k) ,j-number of equation
817 
818  j=nequi+lpre+4
819 
820  do k=1,nequi
821  asum=0.d0
822  do l=1,lfit
823  asum=asum+w_wght(l)*gindk(k,l)
824  enddo
825  a(j,k)=asum
826  enddo
827 
828 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
829 
830  do l=1,lpre
831  k=nequi+l
832  a(j,k)=1.d0
833  enddo
834 
835 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
836  a(j,nequi+lpre+1)=0.d0
837 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
838  a(j,nequi+lpre+2)=0.d0
839 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
840  a(j,nequi+lpre+3)=1.d0
841 
842 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
843 
844  asum=0.d0
845  do l=1,lfit
846  asum=asum+w_wght(l)
847  enddo
848  a(j,nequi+lpre+4)=-asum
849 !!!
850 ! right hand side
851  asum=0.d0
852  do l=1,lfit
853  asum=asum+w_wght(l)*psifit(l)
854  enddo
855  y(j)=-asum
856 
857 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
858  call ge(nequi+lpre+4,ncf_p,a,y,x,ip)
859 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
860 
861  do iq=1,nequi
862  pfceqw(iq)=x(iq)
863  enddo
864 
865 !!!!!!checking!!!!!
866  do l=1,lfit
867  psictr(l)=psifit(l)
868  do iq=1,nequi
869  psictr(l)=psictr(l)+gindk(iq,l)*pfceqw(iq)
870  enddo
871  enddo
872 
873  gr_xp=dpsx_r
874  gz_xp=dpsx_z
875  ps_xp=psip_x
876  ps_ma=psip_a
877  do iq=1,nequi
878  gr_xp=gr_xp+gx_r(iq)*pfceqw(iq)
879  gz_xp=gz_xp+gx_z(iq)*pfceqw(iq)
880 
881  ps_xp=ps_xp+gindx(iq)*pfceqw(iq)
882  ps_ma=ps_ma+ginda(iq)*pfceqw(iq)
883  enddo
884  write(*,*) 'cursol:'
885  write(*,*) 'gr_xp=',gr_xp
886  write(*,*) 'gz_xp=',gz_xp
887 
888  do l=1,lpre
889  psictr(lfit+l)=psipre(l)
890  do iq=1,nequi
891  psictr(lfit+l)=psictr(lfit+l)+gindp(iq,l)*pfceqw(iq)
892  enddo
893  enddo
894 
895  do iq=1,nequi
896  pfceqw(iq)=x(iq)/amu0
897  enddo
898 
899  psi_bnd=x(nequi+lpre+4)
900 
901  alp_out=(ps_ma-psi_bnd)/(ps_ma-ps_xp)
902 
903  return
904  end
905 
906 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
907  subroutine cursol_snf(NEQUI,PFCEQW,psi_bnd)
908 
909  include 'double.inc'
910  include 'parcur.inc'
911 
912  dimension a(ncf_p,ncf_p),x(ncf_p),y(ncf_p),ip(ncf_p)
913  dimension psictr(ncf_p)
914  dimension pfceqw(*)
915 
916 ! a(j,k)*I(k)+a(j,NEQUI+l)*Lam(l)+a(j,NEQUI+Lpre+1)*Psi_b=y(j)
917 
918  pi=3.14159265359d0
919  amu0=0.4d0*pi
920  !amu0=1.0d0
921 
922 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
923 ! equation (.)*d(I_j)=0
924 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
925 
926 ! a(j,k)*I(k) ,j-number of equation
927 
928  do k=1,nequi
929  do j=1,nequi
930  asum=0.d0
931  do l=1,lfit
932  asum=asum+gindk(k,l)*gindk(j,l)*w_wght(l)
933  enddo
934  a(j,k)=asum
935  enddo
936  enddo
937 
938  do j=1,nequi
939  a(j,j)=a(j,j)+d_wght(j)*s_wght
940  enddo
941 
942 !!!!!!!!!!!! + C'*C*I
943 
944  do k=1,nequi
945  do j=1,nequi
946  a(j,k)=a(j,k)+cxc(j,k)*c_wght
947  enddo
948  enddo
949 
950 !!!!!!!!!!!!
951 
952 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
953 
954  do j=1,nequi
955  do l=1,lpre
956  a(j,nequi+l)=gindp(j,l)
957  enddo
958  enddo
959 
960 ! a(j,NEQUI+Lpre+1)*Lam_r(l) ,j-number of equation
961 
962  do j=1,nequi
963  a(j,nequi+lpre+1)=gx_r(j)
964  enddo
965 
966 ! a(j,NEQUI+Lpre+2)*Lam_z(l) ,j-number of equation
967 
968  do j=1,nequi
969  a(j,nequi+lpre+2)=gx_z(j)
970  enddo
971 
972 
973 
974 
975 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
976 
977  do j=1,nequi
978  a(j,nequi+lpre+3)=(1.d0-alph)*ginda(j)+alph*gindx(j)
979  enddo
980 
981 
982 
983 
984 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
985 
986  do j=1,nequi
987  asum=0.d0
988  do l=1,lfit
989  asum=asum+w_wght(l)*gindk(j,l)
990  enddo
991  a(j,nequi+lpre+4)=-asum
992  enddo
993 
994 ! a(j,NEQUI+Lpre+5)*Lam_rz(l) ,j-number of equation
995 
996  do j=1,nequi
997  a(j,nequi+lpre+5)=gx_rz(j)
998  enddo
999 
1000 ! a(j,NEQUI+Lpre+6)*Lam_zz(l) ,j-number of equation
1001 
1002  do j=1,nequi
1003  a(j,nequi+lpre+6)=gx_zz(j)
1004  enddo
1005 
1006 ! right hand side
1007 
1008  do j=1,nequi
1009  asum=0.d0
1010  do l=1,lfit
1011  asum=asum+w_wght(l)*gindk(j,l)*psifit(l)
1012  enddo
1013  y(j)=-asum+d_wght(j)*s_wght*curref(j)*amu0
1014  enddo
1015 
1016 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1017 ! equation (.)*d(Lam_l)=0
1018 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1019 
1020 ! a(j,k)*I(k) ,j-number of equation
1021 
1022  do l=1,lpre
1023  j=nequi+l
1024  do k=1,nequi
1025  a(j,k)=gindp(k,l)
1026  enddo
1027  enddo
1028 
1029 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1030 
1031  do l=1,lpre
1032  j=nequi+l
1033  do ll=1,lpre
1034  k=nequi+ll
1035  a(j,k)=0.d0
1036  enddo
1037  enddo
1038 
1039 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1040 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1041 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1042 
1043  do l=1,lpre
1044  j=nequi+l
1045  k=nequi+lpre+1
1046  a(j,k)=0.d0
1047  k=nequi+lpre+2
1048  a(j,k)=0.d0
1049  k=nequi+lpre+3
1050  a(j,k)=0.d0
1051  enddo
1052 
1053 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1054 
1055  do l=1,lpre
1056  j=nequi+l
1057  k=nequi+lpre+4
1058  a(j,k)=-1.d0
1059  enddo
1060 
1061 
1062 ! a(j,NEQUI+Lpre+5)*Lam_rz ,j-number of equation
1063 ! a(j,NEQUI+Lpre+6)*Lam_zz ,j-number of equation
1064 
1065  do l=1,lpre
1066  j=nequi+l
1067  k=nequi+lpre+5
1068  a(j,k)=0.d0
1069  k=nequi+lpre+6
1070  a(j,k)=0.d0
1071  enddo
1072 
1073 ! right hand side
1074  do l=1,lpre
1075  j=nequi+l
1076  y(j)=-psipre(l)
1077  enddo
1078 
1079 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1080 ! equation (.)*d(Lam_r)=0
1081 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1082 
1083 ! a(j,k)*I(k) ,j-number of equation
1084 
1085  j=nequi+lpre+1
1086 
1087  do k=1,nequi
1088  a(j,k)=gx_r(k)
1089  enddo
1090 
1091 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1092 
1093  do ll=1,lpre
1094  k=nequi+ll
1095  a(j,k)=0.d0
1096  enddo
1097 
1098 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1099 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1100 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1101 
1102  k=nequi+lpre+1
1103  a(j,k)=0.d0
1104  k=nequi+lpre+2
1105  a(j,k)=0.d0
1106  k=nequi+lpre+3
1107  a(j,k)=0.d0
1108 
1109 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1110 
1111  k=nequi+lpre+4
1112  a(j,k)=0.d0
1113 
1114 ! a(j,NEQUI+Lpre+5)*Lam_rz ,j-number of equation
1115 ! a(j,NEQUI+Lpre+6)*Lam_zz ,j-number of equation
1116 
1117  k=nequi+lpre+5
1118  a(j,k)=0.d0
1119  k=nequi+lpre+6
1120  a(j,k)=0.d0
1121 
1122 ! right hand side
1123  y(j)=-dpsx_r
1124 
1125 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1126 ! equation (.)*d(Lam_z)=0
1127 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1128 
1129 ! a(j,k)*I(k) ,j-number of equation
1130 
1131  j=nequi+lpre+2
1132 
1133  do k=1,nequi
1134  a(j,k)=gx_z(k)
1135  enddo
1136 
1137 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1138 
1139  do ll=1,lpre
1140  k=nequi+ll
1141  a(j,k)=0.d0
1142  enddo
1143 
1144 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1145 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1146 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1147 
1148  k=nequi+lpre+1
1149  a(j,k)=0.d0
1150  k=nequi+lpre+2
1151  a(j,k)=0.d0
1152  k=nequi+lpre+3
1153  a(j,k)=0.d0
1154 
1155 ! a(j,NEQUI+Lpre+3)*Psi_b ,j-number of equation
1156 
1157  k=nequi+lpre+4
1158  a(j,k)=0.d0
1159 
1160 ! a(j,NEQUI+Lpre+5)*Lam_rz ,j-number of equation
1161 ! a(j,NEQUI+Lpre+6)*Lam_zz ,j-number of equation
1162 
1163  k=nequi+lpre+5
1164  a(j,k)=0.d0
1165  k=nequi+lpre+6
1166  a(j,k)=0.d0
1167 
1168 ! right hand side
1169  y(j)=-dpsx_z
1170 
1171 
1172 
1173 
1174 
1175 
1176 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1177 ! equation (.)*d(Lamx)=0
1178 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1179 
1180 ! a(j,k)*I(k) ,j-number of equation
1181 
1182  j=nequi+lpre+3
1183 
1184  do k=1,nequi
1185  a(j,k)=(1.d0-alph)*ginda(k)+alph*gindx(k)
1186  enddo
1187 
1188 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1189 
1190  do ll=1,lpre
1191  k=nequi+ll
1192  a(j,k)=0.d0
1193  enddo
1194 
1195 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1196 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1197 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1198 
1199  k=nequi+lpre+1
1200  a(j,k)=0.d0
1201  k=nequi+lpre+2
1202  a(j,k)=0.d0
1203  k=nequi+lpre+3
1204  a(j,k)=0.d0
1205 
1206 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1207 
1208  k=nequi+lpre+4
1209  a(j,k)=-1.d0
1210 
1211 ! a(j,NEQUI+Lpre+5)*Lam_rz ,j-number of equation
1212 ! a(j,NEQUI+Lpre+6)*Lam_zz ,j-number of equation
1213 
1214  k=nequi+lpre+5
1215  a(j,k)=0.d0
1216  k=nequi+lpre+6
1217  a(j,k)=0.d0
1218 
1219 ! right hand side
1220  y(j)=-(1.d0-alph)*psip_a-alph*psip_x
1221 
1222 
1223 
1224 
1225 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1226 ! equation (.)*d(-psi_b)=0
1227 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1228 
1229 ! a(j,k)*I(k) ,j-number of equation
1230 
1231  j=nequi+lpre+4
1232 
1233  do k=1,nequi
1234  asum=0.d0
1235  do l=1,lfit
1236  asum=asum+w_wght(l)*gindk(k,l)
1237  enddo
1238  a(j,k)=asum
1239  enddo
1240 
1241 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1242 
1243  do l=1,lpre
1244  k=nequi+l
1245  a(j,k)=1.d0
1246  enddo
1247 
1248 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1249  a(j,nequi+lpre+1)=0.d0
1250 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1251  a(j,nequi+lpre+2)=0.d0
1252 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1253  a(j,nequi+lpre+3)=1.d0
1254 
1255 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1256 
1257  asum=0.d0
1258  do l=1,lfit
1259  asum=asum+w_wght(l)
1260  enddo
1261  a(j,nequi+lpre+4)=-asum
1262 
1263 ! a(j,NEQUI+Lpre+5)*Lam_rz ,j-number of equation
1264  a(j,nequi+lpre+5)=0.d0
1265 ! a(j,NEQUI+Lpre+2)*Lam_zz ,j-number of equation
1266  a(j,nequi+lpre+6)=0.d0
1267 
1268 !!!
1269 ! right hand side
1270  asum=0.d0
1271  do l=1,lfit
1272  asum=asum+w_wght(l)*psifit(l)
1273  enddo
1274  y(j)=-asum
1275 
1276 
1277 
1278 
1279 
1280 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1281 ! equation (.)*d(Lam_rz)=0
1282 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1283 
1284 ! a(j,k)*I(k) ,j-number of equation
1285 
1286  j=nequi+lpre+5
1287 
1288  do k=1,nequi
1289  a(j,k)=gx_rz(k)
1290  enddo
1291 
1292 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1293 
1294  do ll=1,lpre
1295  k=nequi+ll
1296  a(j,k)=0.d0
1297  enddo
1298 
1299 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1300 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1301 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1302 
1303  k=nequi+lpre+1
1304  a(j,k)=0.d0
1305  k=nequi+lpre+2
1306  a(j,k)=0.d0
1307  k=nequi+lpre+3
1308  a(j,k)=0.d0
1309 
1310 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1311 
1312  k=nequi+lpre+4
1313  a(j,k)=0.d0
1314 ! a(j,NEQUI+Lpre+5)*Lam_rz ,j-number of equation
1315  k=nequi+lpre+5
1316  a(j,k)=0.d0
1317 ! a(j,NEQUI+Lpre+6)*Lam_zz ,j-number of equation
1318  k=nequi+lpre+6
1319  a(j,k)=0.d0
1320 ! right hand side
1321  y(j)=-dpsx_rz
1322 
1323 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1324 ! equation (.)*d(Lam_zz)=0
1325 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1326 
1327 ! a(j,k)*I(k) ,j-number of equation
1328 
1329  j=nequi+lpre+6
1330 
1331  do k=1,nequi
1332  a(j,k)=gx_zz(k)
1333  enddo
1334 
1335 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1336 
1337  do ll=1,lpre
1338  k=nequi+ll
1339  a(j,k)=0.d0
1340  enddo
1341 
1342 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1343 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1344 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1345 
1346  k=nequi+lpre+1
1347  a(j,k)=0.d0
1348  k=nequi+lpre+2
1349  a(j,k)=0.d0
1350  k=nequi+lpre+3
1351  a(j,k)=0.d0
1352 
1353 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1354 
1355  k=nequi+lpre+4
1356  a(j,k)=0.d0
1357 ! a(j,NEQUI+Lpre+5)*Lam_rz ,j-number of equation
1358  k=nequi+lpre+5
1359 ! a(j,NEQUI+Lpre+6)*Lam_zz ,j-number of equation
1360  a(j,k)=0.d0
1361  k=nequi+lpre+6
1362  a(j,k)=0.d0
1363 ! right hand side
1364  y(j)=-dpsx_zz
1365 
1366 
1367 
1368 
1369 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1370  call ge(nequi+lpre+6,ncf_p,a,y,x,ip)
1371 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1372 
1373  do iq=1,nequi
1374  pfceqw(iq)=x(iq)
1375  enddo
1376 
1377 !!!!!!checking!!!!!
1378  do l=1,lfit
1379  psictr(l)=psifit(l)
1380  do iq=1,nequi
1381  psictr(l)=psictr(l)+gindk(iq,l)*pfceqw(iq)
1382  enddo
1383  enddo
1384 
1385  gr_xp=dpsx_r
1386  gz_xp=dpsx_z
1387  ps_xp=psip_x
1388  ps_ma=psip_a
1389  d2rz_xp=dpsx_rz
1390  d2zz_xp=dpsx_zz
1391  d2rr_xp=dpsx_rr
1392 
1393  do iq=1,nequi
1394 
1395  gr_xp=gr_xp+gx_r(iq)*pfceqw(iq)
1396  gz_xp=gz_xp+gx_z(iq)*pfceqw(iq)
1397 
1398  ps_xp=ps_xp+gindx(iq)*pfceqw(iq)
1399  ps_ma=ps_ma+ginda(iq)*pfceqw(iq)
1400 
1401  d2rz_xp=d2rz_xp+gx_rz(iq)*pfceqw(iq)
1402  d2zz_xp=d2zz_xp+gx_zz(iq)*pfceqw(iq)
1403  d2rr_xp=d2rr_xp+gx_rr(iq)*pfceqw(iq)
1404 
1405  enddo
1406 
1407  write(*,*) 'cursol:'
1408  write(*,*) 'd2rz_xp=',d2rz_xp
1409  write(*,*) 'd2zz_xp=',d2zz_xp
1410  write(*,*) 'd2rr_xp=',d2rr_xp
1411  write(*,*) 'gr_xp=',gr_xp
1412  write(*,*) 'gz_xp=',gz_xp
1413 
1414  do l=1,lpre
1415  psictr(lfit+l)=psipre(l)
1416  do iq=1,nequi
1417  psictr(lfit+l)=psictr(lfit+l)+gindp(iq,l)*pfceqw(iq)
1418  enddo
1419  enddo
1420 
1421  do iq=1,nequi
1422  pfceqw(iq)=x(iq)/amu0
1423  enddo
1424 
1425  psi_bnd=x(nequi+lpre+4)
1426 
1427  alp_out=(ps_ma-psi_bnd)/(ps_ma-ps_xp)
1428 
1429  return
1430  end
1431 
1432 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1433 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1434  subroutine cursol_2x(NEQUI,PFCEQW,psi_bnd)
1435 
1436  include 'double.inc'
1437  include 'parcur.inc'
1438 
1439  dimension a(ncf_p,ncf_p),x(ncf_p),y(ncf_p),ip(ncf_p)
1440  dimension psictr(ncf_p)
1441  dimension pfceqw(*)
1442 
1443 ! a(j,k)*I(k)+a(j,NEQUI+l)*Lam(l)+a(j,NEQUI+Lpre+1)*Psi_b=y(j)
1444 
1445  pi=3.14159265359d0
1446  amu0=0.4d0*pi
1447  !amu0=1.0d0
1448 
1449 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1450 ! equation (.)*d(I_j)=0
1451 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1452 
1453 ! a(j,k)*I(k) ,j-number of equation
1454 
1455  do k=1,nequi
1456  do j=1,nequi
1457  asum=0.d0
1458  do l=1,lfit
1459  asum=asum+gindk(k,l)*gindk(j,l)*w_wght(l)
1460  enddo
1461  a(j,k)=asum
1462  enddo
1463  enddo
1464 
1465  do j=1,nequi
1466  a(j,j)=a(j,j)+d_wght(j)*s_wght
1467  enddo
1468 
1469 !!!!!!!!!!!! + C'*C*I
1470 
1471  do k=1,nequi
1472  do j=1,nequi
1473  a(j,k)=a(j,k)+cxc(j,k)*c_wght
1474  enddo
1475  enddo
1476 
1477 !!!!!!!!!!!!
1478 
1479 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1480 
1481  do j=1,nequi
1482  do l=1,lpre
1483  a(j,nequi+l)=gindp(j,l)
1484  enddo
1485  enddo
1486 
1487 ! a(j,NEQUI+Lpre+1)*Lam_r(l) ,j-number of equation
1488 
1489  do j=1,nequi
1490  a(j,nequi+lpre+1)=gx_r(j)
1491  enddo
1492 
1493 ! a(j,NEQUI+Lpre+2)*Lam_z(l) ,j-number of equation
1494 
1495  do j=1,nequi
1496  a(j,nequi+lpre+2)=gx_z(j)
1497  enddo
1498 
1499 
1500 
1501 
1502 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1503 
1504  do j=1,nequi
1505  a(j,nequi+lpre+3)=(1.d0-alph)*ginda(j)+alph*gindx(j)
1506  enddo
1507 
1508 
1509 
1510 
1511 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1512 
1513  do j=1,nequi
1514  asum=0.d0
1515  do l=1,lfit
1516  asum=asum+w_wght(l)*gindk(j,l)
1517  enddo
1518  a(j,nequi+lpre+4)=-asum
1519  enddo
1520 
1521 ! a(j,NEQUI+Lpre+5)*Lam_r2(l) ,j-number of equation
1522 
1523  do j=1,nequi
1524  a(j,nequi+lpre+5)=gx2_r(j)
1525  enddo
1526 
1527 ! a(j,NEQUI+Lpre+6)*Lam_z2(l) ,j-number of equation
1528 
1529  do j=1,nequi
1530  a(j,nequi+lpre+6)=gx2_z(j)
1531  enddo
1532 
1533 ! right hand side
1534 
1535  do j=1,nequi
1536  asum=0.d0
1537  do l=1,lfit
1538  asum=asum+w_wght(l)*gindk(j,l)*psifit(l)
1539  enddo
1540  y(j)=-asum+d_wght(j)*s_wght*curref(j)*amu0
1541  enddo
1542 
1543 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1544 ! equation (.)*d(Lam_l)=0
1545 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1546 
1547 ! a(j,k)*I(k) ,j-number of equation
1548 
1549  do l=1,lpre
1550  j=nequi+l
1551  do k=1,nequi
1552  a(j,k)=gindp(k,l)
1553  enddo
1554  enddo
1555 
1556 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1557 
1558  do l=1,lpre
1559  j=nequi+l
1560  do ll=1,lpre
1561  k=nequi+ll
1562  a(j,k)=0.d0
1563  enddo
1564  enddo
1565 
1566 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1567 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1568 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1569 
1570  do l=1,lpre
1571  j=nequi+l
1572  k=nequi+lpre+1
1573  a(j,k)=0.d0
1574  k=nequi+lpre+2
1575  a(j,k)=0.d0
1576  k=nequi+lpre+3
1577  a(j,k)=0.d0
1578  enddo
1579 
1580 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1581 
1582  do l=1,lpre
1583  j=nequi+l
1584  k=nequi+lpre+4
1585  a(j,k)=-1.d0
1586  enddo
1587 
1588 
1589 ! a(j,NEQUI+Lpre+5)*Lam_r2 ,j-number of equation
1590 ! a(j,NEQUI+Lpre+6)*Lam_z2 ,j-number of equation
1591 
1592  do l=1,lpre
1593  j=nequi+l
1594  k=nequi+lpre+5
1595  a(j,k)=0.d0
1596  k=nequi+lpre+6
1597  a(j,k)=0.d0
1598  enddo
1599 
1600 ! right hand side
1601  do l=1,lpre
1602  j=nequi+l
1603  y(j)=-psipre(l)
1604  enddo
1605 
1606 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1607 ! equation (.)*d(Lam_r)=0
1608 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1609 
1610 ! a(j,k)*I(k) ,j-number of equation
1611 
1612  j=nequi+lpre+1
1613 
1614  do k=1,nequi
1615  a(j,k)=gx_r(k)
1616  enddo
1617 
1618 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1619 
1620  do ll=1,lpre
1621  k=nequi+ll
1622  a(j,k)=0.d0
1623  enddo
1624 
1625 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1626 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1627 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1628 
1629  k=nequi+lpre+1
1630  a(j,k)=0.d0
1631  k=nequi+lpre+2
1632  a(j,k)=0.d0
1633  k=nequi+lpre+3
1634  a(j,k)=0.d0
1635 
1636 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1637 
1638  k=nequi+lpre+4
1639  a(j,k)=0.d0
1640 
1641 ! a(j,NEQUI+Lpre+5)*Lam_r2 ,j-number of equation
1642 ! a(j,NEQUI+Lpre+6)*Lam_z2 ,j-number of equation
1643 
1644  k=nequi+lpre+5
1645  a(j,k)=0.d0
1646  k=nequi+lpre+6
1647  a(j,k)=0.d0
1648 
1649 ! right hand side
1650  y(j)=-dpsx_r
1651 
1652 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1653 ! equation (.)*d(Lam_z)=0
1654 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1655 
1656 ! a(j,k)*I(k) ,j-number of equation
1657 
1658  j=nequi+lpre+2
1659 
1660  do k=1,nequi
1661  a(j,k)=gx_z(k)
1662  enddo
1663 
1664 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1665 
1666  do ll=1,lpre
1667  k=nequi+ll
1668  a(j,k)=0.d0
1669  enddo
1670 
1671 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1672 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1673 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1674 
1675  k=nequi+lpre+1
1676  a(j,k)=0.d0
1677  k=nequi+lpre+2
1678  a(j,k)=0.d0
1679  k=nequi+lpre+3
1680  a(j,k)=0.d0
1681 
1682 ! a(j,NEQUI+Lpre+3)*Psi_b ,j-number of equation
1683 
1684  k=nequi+lpre+4
1685  a(j,k)=0.d0
1686 
1687 ! a(j,NEQUI+Lpre+5)*Lam_r2 ,j-number of equation
1688 ! a(j,NEQUI+Lpre+6)*Lam_z2 ,j-number of equation
1689 
1690  k=nequi+lpre+5
1691  a(j,k)=0.d0
1692  k=nequi+lpre+6
1693  a(j,k)=0.d0
1694 
1695 ! right hand side
1696  y(j)=-dpsx_z
1697 
1698 
1699 
1700 
1701 
1702 
1703 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1704 ! equation (.)*d(Lamx)=0
1705 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1706 
1707 ! a(j,k)*I(k) ,j-number of equation
1708 
1709  j=nequi+lpre+3
1710 
1711  do k=1,nequi
1712  a(j,k)=(1.d0-alph)*ginda(k)+alph*gindx(k)
1713  enddo
1714 
1715 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1716 
1717  do ll=1,lpre
1718  k=nequi+ll
1719  a(j,k)=0.d0
1720  enddo
1721 
1722 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1723 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1724 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1725 
1726  k=nequi+lpre+1
1727  a(j,k)=0.d0
1728  k=nequi+lpre+2
1729  a(j,k)=0.d0
1730  k=nequi+lpre+3
1731  a(j,k)=0.d0
1732 
1733 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1734 
1735  k=nequi+lpre+4
1736  a(j,k)=-1.d0
1737 
1738 ! a(j,NEQUI+Lpre+5)*Lam_r2 ,j-number of equation
1739 ! a(j,NEQUI+Lpre+6)*Lam_z2 ,j-number of equation
1740 
1741  k=nequi+lpre+5
1742  a(j,k)=0.d0
1743  k=nequi+lpre+6
1744  a(j,k)=0.d0
1745 
1746 ! right hand side
1747  y(j)=-(1.d0-alph)*psip_a-alph*psip_x
1748 
1749 
1750 
1751 
1752 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1753 ! equation (.)*d(-psi_b)=0
1754 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1755 
1756 ! a(j,k)*I(k) ,j-number of equation
1757 
1758  j=nequi+lpre+4
1759 
1760  do k=1,nequi
1761  asum=0.d0
1762  do l=1,lfit
1763  asum=asum+w_wght(l)*gindk(k,l)
1764  enddo
1765  a(j,k)=asum
1766  enddo
1767 
1768 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1769 
1770  do l=1,lpre
1771  k=nequi+l
1772  a(j,k)=1.d0
1773  enddo
1774 
1775 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1776  a(j,nequi+lpre+1)=0.d0
1777 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1778  a(j,nequi+lpre+2)=0.d0
1779 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1780  a(j,nequi+lpre+3)=1.d0
1781 
1782 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1783 
1784  asum=0.d0
1785  do l=1,lfit
1786  asum=asum+w_wght(l)
1787  enddo
1788  a(j,nequi+lpre+4)=-asum
1789 
1790 ! a(j,NEQUI+Lpre+5)*Lam_r2 ,j-number of equation
1791  a(j,nequi+lpre+5)=0.d0
1792 ! a(j,NEQUI+Lpre+2)*Lam_z2 ,j-number of equation
1793  a(j,nequi+lpre+6)=0.d0
1794 
1795 !!!
1796 ! right hand side
1797  asum=0.d0
1798  do l=1,lfit
1799  asum=asum+w_wght(l)*psifit(l)
1800  enddo
1801  y(j)=-asum
1802 
1803 
1804 
1805 
1806 
1807 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1808 ! equation (.)*d(Lam_r2)=0
1809 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1810 
1811 ! a(j,k)*I(k) ,j-number of equation
1812 
1813  j=nequi+lpre+5
1814 
1815  do k=1,nequi
1816  a(j,k)=gx2_r(k)
1817  enddo
1818 
1819 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1820 
1821  do ll=1,lpre
1822  k=nequi+ll
1823  a(j,k)=0.d0
1824  enddo
1825 
1826 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1827 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1828 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1829 
1830  k=nequi+lpre+1
1831  a(j,k)=0.d0
1832  k=nequi+lpre+2
1833  a(j,k)=0.d0
1834  k=nequi+lpre+3
1835  a(j,k)=0.d0
1836 
1837 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1838 
1839  k=nequi+lpre+4
1840  a(j,k)=0.d0
1841 ! a(j,NEQUI+Lpre+5)*Lam_r2 ,j-number of equation
1842  k=nequi+lpre+5
1843  a(j,k)=0.d0
1844 ! a(j,NEQUI+Lpre+6)*Lam_z2 ,j-number of equation
1845  k=nequi+lpre+6
1846  a(j,k)=0.d0
1847 ! right hand side
1848  y(j)=-dpsx2_r
1849 
1850 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1851 ! equation (.)*d(Lam_z2)=0
1852 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1853 
1854 ! a(j,k)*I(k) ,j-number of equation
1855 
1856  j=nequi+lpre+6
1857 
1858  do k=1,nequi
1859  a(j,k)=gx2_z(k)
1860  enddo
1861 
1862 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1863 
1864  do ll=1,lpre
1865  k=nequi+ll
1866  a(j,k)=0.d0
1867  enddo
1868 
1869 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
1870 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
1871 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
1872 
1873  k=nequi+lpre+1
1874  a(j,k)=0.d0
1875  k=nequi+lpre+2
1876  a(j,k)=0.d0
1877  k=nequi+lpre+3
1878  a(j,k)=0.d0
1879 
1880 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
1881 
1882  k=nequi+lpre+4
1883  a(j,k)=0.d0
1884 ! a(j,NEQUI+Lpre+5)*Lam_r2 ,j-number of equation
1885  k=nequi+lpre+5
1886 ! a(j,NEQUI+Lpre+6)*Lam_z2 ,j-number of equation
1887  a(j,k)=0.d0
1888  k=nequi+lpre+6
1889  a(j,k)=0.d0
1890 ! right hand side
1891  y(j)=-dpsx2_z
1892 
1893 
1894 
1895 
1896 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1897  call ge(nequi+lpre+6,ncf_p,a,y,x,ip)
1898 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1899 
1900  do iq=1,nequi
1901  pfceqw(iq)=x(iq)
1902  enddo
1903 
1904 !!!!!!checking!!!!!
1905  do l=1,lfit
1906  psictr(l)=psifit(l)
1907  do iq=1,nequi
1908  psictr(l)=psictr(l)+gindk(iq,l)*pfceqw(iq)
1909  enddo
1910  enddo
1911 
1912  gr_xp=dpsx_r
1913  gz_xp=dpsx_z
1914  ps_xp=psip_x
1915  ps_ma=psip_a
1916  gr_xp2=dpsx2_r
1917  gz_xp2=dpsx2_z
1918 
1919  do iq=1,nequi
1920 
1921  gr_xp=gr_xp+gx_r(iq)*pfceqw(iq)
1922  gz_xp=gz_xp+gx_z(iq)*pfceqw(iq)
1923 
1924  ps_xp=ps_xp+gindx(iq)*pfceqw(iq)
1925  ps_ma=ps_ma+ginda(iq)*pfceqw(iq)
1926 
1927  gr_xp2=gr_xp2+gx2_r(iq)*pfceqw(iq)
1928  gz_xp2=gz_xp2+gx2_z(iq)*pfceqw(iq)
1929 
1930  enddo
1931 
1932  write(*,*) 'cursol:'
1933  write(*,*) 'gr_xp=',gr_xp
1934  write(*,*) 'gz_xp=',gz_xp
1935  write(*,*) 'gr_xp2=',gr_xp2
1936  write(*,*) 'gz_xp2=',gz_xp2
1937 
1938  do l=1,lpre
1939  psictr(lfit+l)=psipre(l)
1940  do iq=1,nequi
1941  psictr(lfit+l)=psictr(lfit+l)+gindp(iq,l)*pfceqw(iq)
1942  enddo
1943  enddo
1944 
1945  do iq=1,nequi
1946  pfceqw(iq)=x(iq)/amu0
1947  enddo
1948 
1949  psi_bnd=x(nequi+lpre+4)
1950 
1951  alp_out=(ps_ma-psi_bnd)/(ps_ma-ps_xp)
1952 
1953  return
1954  end
1955 
1956 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1957 
1958  subroutine cursol_wx(NEQUI,PFCEQW,psi_bnd)
1959 
1960  include 'double.inc'
1961  include 'parcur.inc'
1962 
1963  dimension a(ncf_p,ncf_p),x(ncf_p),y(ncf_p),ip(ncf_p)
1964  dimension psictr(ncf_p)
1965  dimension pfceqw(*)
1966 
1967 ! a(j,k)*I(k)+a(j,NEQUI+l)*Lam(l)+a(j,NEQUI+Lpre+1)*Psi_b=y(j)
1968 
1969  pi=3.14159265359d0
1970  amu0=0.4d0*pi
1971  !amu0=1.0d0
1972 
1973 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1974 ! equation (.)*d(I_j)=0
1975 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1976 
1977 ! a(j,k)*I(k) ,j-number of equation
1978 
1979  do k=1,nequi
1980  do j=1,nequi
1981  asum=0.d0
1982  do l=1,lfit
1983  asum=asum+gindk(k,l)*gindk(j,l)*w_wght(l)
1984  enddo
1985  a(j,k)=asum
1986  enddo
1987  enddo
1988 
1989  do j=1,nequi
1990  a(j,j)=a(j,j)+d_wght(j)*s_wght
1991  enddo
1992 
1993 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
1994 
1995  do j=1,nequi
1996  do l=1,lpre
1997  a(j,nequi+l)=gindp(j,l)
1998  enddo
1999  enddo
2000 
2001 
2002 ! *a(j,NEQUI+Lpre+1)*Lam_r(l) ,j-number of equation
2003 !
2004 ! do j=1,NEQUI
2005 ! a(j,NEQUI+Lpre+1)=Gx_r(j)
2006 ! enddo
2007 !
2008 ! *a(j,NEQUI+Lpre+2)*Lam_z(l) ,j-number of equation
2009 !
2010 ! do j=1,NEQUI
2011 ! a(j,NEQUI+Lpre+2)=Gx_z(j)
2012 ! enddo
2013 
2014 
2015 
2016 ! a(j,NEQUI+Lpre+1)*Lamx ,j-number of equation
2017 
2018  do j=1,nequi
2019  a(j,nequi+lpre+1)=(1.d0-alph)*ginda(j)+alph*gindx(j)
2020  enddo
2021 
2022 
2023 
2024 
2025 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
2026 
2027  do j=1,nequi
2028  asum=0.d0
2029  do l=1,lfit
2030  asum=asum+w_wght(l)*gindk(j,l)
2031  enddo
2032  a(j,nequi+lpre+2)=-asum
2033  enddo
2034 
2035 ! right hand side
2036 
2037  do j=1,nequi
2038  asum=0.d0
2039  do l=1,lfit
2040  asum=asum+w_wght(l)*gindk(j,l)*psifit(l)
2041  enddo
2042  y(j)=-asum+d_wght(j)*s_wght*curref(j)*amu0
2043  enddo
2044 
2045 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2046 ! equation (.)*d(Lam_l)=0
2047 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2048 
2049 ! a(j,k)*I(k) ,j-number of equation
2050 
2051  do l=1,lpre
2052  j=nequi+l
2053  do k=1,nequi
2054  a(j,k)=gindp(k,l)
2055  enddo
2056  enddo
2057 
2058 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
2059 
2060  do l=1,lpre
2061  j=nequi+l
2062  do ll=1,lpre
2063  k=nequi+ll
2064  a(j,k)=0.d0
2065  enddo
2066  enddo
2067 
2068 ! *a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
2069 ! *a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
2070 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
2071 
2072  do l=1,lpre
2073  j=nequi+l
2074  k=nequi+lpre+1
2075  a(j,k)=0.d0
2076 ! k=NEQUI+Lpre+2
2077 ! a(j,k)=0.d0
2078 ! k=NEQUI+Lpre+3
2079 ! a(j,k)=0.d0
2080  enddo
2081 
2082 ! a(j,NEQUI+Lpre+2)*Psi_b ,j-number of equation
2083 
2084  do l=1,lpre
2085  j=nequi+l
2086  k=nequi+lpre+2
2087  a(j,k)=-1.d0
2088  enddo
2089 ! right hand side
2090  do l=1,lpre
2091  j=nequi+l
2092  y(j)=-psipre(l)
2093  enddo
2094 
2095 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2096 ! equation (.)*d(Lam_r)=0
2097 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2098 
2099 ! a(j,k)*I(k) ,j-number of equation
2100 !
2101 ! j=NEQUI+Lpre+1
2102 !
2103 ! do k=1,NEQUI
2104 ! a(j,k)=Gx_r(k)
2105 ! enddo
2106 !
2107 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
2108 !
2109 ! do ll=1,Lpre
2110 ! k=NEQUI+ll
2111 ! a(j,k)=0.d0
2112 ! enddo
2113 !
2114 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
2115 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
2116 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
2117 !
2118 ! k=NEQUI+Lpre+1
2119 ! a(j,k)=0.d0
2120 ! k=NEQUI+Lpre+2
2121 ! a(j,k)=0.d0
2122 ! k=NEQUI+Lpre+3
2123 ! a(j,k)=0.d0
2124 !
2125 ! a(j,NEQUI+Lpre+4)*Psi_b ,j-number of equation
2126 !
2127 ! k=NEQUI+Lpre+4
2128 ! a(j,k)=0.d0
2129 ! right hand side
2130 ! y(j)=-dpsx_r
2131 !
2132 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2133 ! equation (.)*d(Lam_z)=0
2134 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2135 !
2136 ! a(j,k)*I(k) ,j-number of equation
2137 !
2138 ! j=NEQUI+Lpre+2
2139 !
2140 ! do k=1,NEQUI
2141 ! a(j,k)=Gx_z(k)
2142 ! enddo
2143 !
2144 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
2145 !
2146 ! do ll=1,Lpre
2147 ! k=NEQUI+ll
2148 ! a(j,k)=0.d0
2149 ! enddo
2150 !
2151 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
2152 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
2153 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
2154 !
2155 ! k=NEQUI+Lpre+1
2156 ! a(j,k)=0.d0
2157 ! k=NEQUI+Lpre+2
2158 ! a(j,k)=0.d0
2159 ! k=NEQUI+Lpre+3
2160 ! a(j,k)=0.d0
2161 !
2162 ! a(j,NEQUI+Lpre+3)*Psi_b ,j-number of equation
2163 !
2164 ! k=NEQUI+Lpre+4
2165 ! a(j,k)=0.d0
2166 ! right hand side
2167 ! y(j)=-dpsx_z
2168 
2169 
2170 
2171 
2172 
2173 
2174 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2175 ! equation (.)*d(Lamx)=0
2176 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2177 
2178 ! a(j,k)*I(k) ,j-number of equation
2179 
2180  j=nequi+lpre+1
2181 
2182  do k=1,nequi
2183  a(j,k)=(1.d0-alph)*ginda(k)+alph*gindx(k)
2184  enddo
2185 
2186 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
2187 
2188  do ll=1,lpre
2189  k=nequi+ll
2190  a(j,k)=0.d0
2191  enddo
2192 
2193 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
2194 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
2195 ! a(j,NEQUI+Lpre+3)*Lamx ,j-number of equation
2196 
2197  k=nequi+lpre+1
2198  a(j,k)=0.d0
2199 ! k=NEQUI+Lpre+2
2200 ! a(j,k)=0.d0
2201 ! k=NEQUI+Lpre+3
2202 ! a(j,k)=0.d0
2203 
2204 ! a(j,NEQUI+Lpre+2)*Psi_b ,j-number of equation
2205 
2206  k=nequi+lpre+2
2207  a(j,k)=-1.d0
2208 
2209 ! right hand side
2210  y(j)=-(1.d0-alph)*psip_a-alph*psip_x
2211 
2212 
2213 
2214 
2215 
2216 
2217 
2218 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2219 ! equation (.)*d(-psi_b)=0
2220 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2221 
2222 ! a(j,k)*I(k) ,j-number of equation
2223 
2224  j=nequi+lpre+2
2225 
2226  do k=1,nequi
2227  asum=0.d0
2228  do l=1,lfit
2229  asum=asum+w_wght(l)*gindk(k,l)
2230  enddo
2231  a(j,k)=asum
2232  enddo
2233 
2234 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
2235 
2236  do l=1,lpre
2237  k=nequi+l
2238  a(j,k)=1.d0
2239  enddo
2240 
2241 ! a(j,NEQUI+Lpre+1)*Lam_r ,j-number of equation
2242 ! a(j,NEQUI+Lpre+1)=0.d0
2243 ! a(j,NEQUI+Lpre+2)*Lam_z ,j-number of equation
2244 ! a(j,NEQUI+Lpre+2)=0.d0
2245 ! a(j,NEQUI+Lpre+1)*Lamx ,j-number of equation
2246  a(j,nequi+lpre+1)=1.d0
2247 
2248 ! a(j,NEQUI+Lpre+2)*Psi_b ,j-number of equation
2249 
2250  asum=0.d0
2251  do l=1,lfit
2252  asum=asum+w_wght(l)
2253  enddo
2254  a(j,nequi+lpre+2)=-asum
2255 !!!
2256 ! right hand side
2257  asum=0.d0
2258  do l=1,lfit
2259  asum=asum+w_wght(l)*psifit(l)
2260  enddo
2261  y(j)=-asum
2262 
2263 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2264  call ge(nequi+lpre+2,ncf_p,a,y,x,ip)
2265 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2266 
2267  do iq=1,nequi
2268  pfceqw(iq)=x(iq)
2269  enddo
2270 
2271 !!!!!!checking!!!!!
2272  do l=1,lfit
2273  psictr(l)=psifit(l)
2274  do iq=1,nequi
2275  psictr(l)=psictr(l)+gindk(iq,l)*pfceqw(iq)
2276  enddo
2277  enddo
2278 
2279  gr_xp=dpsx_r
2280  gz_xp=dpsx_z
2281  ps_xp=psip_x
2282  ps_ma=psip_a
2283  do iq=1,nequi
2284  gr_xp=gr_xp+gx_r(iq)*pfceqw(iq)
2285  gz_xp=gz_xp+gx_z(iq)*pfceqw(iq)
2286 
2287  ps_xp=ps_xp+gindx(iq)*pfceqw(iq)
2288  ps_ma=ps_ma+ginda(iq)*pfceqw(iq)
2289  enddo
2290 
2291  do l=1,lpre
2292  psictr(lfit+l)=psipre(l)
2293  do iq=1,nequi
2294  psictr(lfit+l)=psictr(lfit+l)+gindp(iq,l)*pfceqw(iq)
2295  enddo
2296  enddo
2297 
2298  do iq=1,nequi
2299  pfceqw(iq)=x(iq)/amu0
2300  enddo
2301 
2302  psi_bnd=x(nequi+lpre+2)
2303 
2304  alp_out=(ps_ma-psi_bnd)/(ps_ma-ps_xp)
2305 
2306  return
2307  end
2308 
2309 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2310 
2311  !subroutine psi_cpn
2312  subroutine psi_cpn2
2313 
2314  include 'double.inc'
2315  parameter(nshp=10)
2316  include 'parcur.inc'
2317  include 'param.inc'
2318  include 'comblc.inc'
2319  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
2320 
2321  do l=1,lfit
2322 
2323  r0=rfit(l)
2324  z0=zfit(l)
2325 
2326  ic=(r0-rmin)/dr(1)+1
2327  jc=(z0-zmin)/dz(1)+1
2328 
2329  psifit(l)=blin(ic,jc,r0,z0)
2330 
2331  enddo
2332 
2333  do l=1,lpre
2334 
2335  r0=rpre(l)
2336  z0=zpre(l)
2337 
2338  ic=(r0-rmin)/dr(1)+1
2339  jc=(z0-zmin)/dz(1)+1
2340 
2341  psipre(l)=blin(ic,jc,r0,z0)
2342 
2343  enddo
2344 
2345  r0=rx_p
2346  z0=zx_p
2347 
2348  ic=(r0-rmin)/dr(1)+1
2349  jc=(z0-zmin)/dz(1)+1
2350 
2351 
2352 c---definition of the nearest knote
2353 
2354  sdmin=rmax
2355 
2356  do 320 k=0,1
2357  rr=r(ic+k)
2358  do 325 l=0,1
2359  zz=z(jc+l)
2360  dlx=dsqrt( (rr-rx_p)**2+(zz-zx_p)**2 )
2361  if(dlx.lt.sdmin) then
2362  sdmin=dlx
2363  ik=ic+k
2364  jk=jc+l
2365  endif
2366  325 continue
2367  320 continue
2368 
2369  nsh=1
2370  xs(nsh)=r(ik)
2371  ys(nsh)=z(jk)
2372  fun(nsh)=ui(ik,jk)
2373 
2374  do 400 k=-1,1
2375 
2376  i= ik+k
2377 
2378  do 410 l=-1,1
2379 
2380  j= jk+l
2381 
2382  if(k.eq.0 .AND. l.eq.0 ) go to 410
2383  nsh=nsh+1
2384  xs(nsh)=r(i)
2385  ys(nsh)=z(j)
2386  fun(nsh)=ui(i,j)
2387 
2388  410 continue
2389  400 continue
2390 
2391  rrx=xs(1)
2392  zzx=ys(1)
2393 
2394  call deriv5(xs,ys,fun,nsh,5,dp)
2395 
2396  dpsx_r = dp(1) + dp(3)*(rx_p-rrx) + dp(4)*(zx_p-zzx)
2397  dpsx_z = dp(2) + dp(5)*(zx_p-zzx) + dp(4)*(rx_p-rrx)
2398 
2399  dpsx_rz = dp(4)
2400  dpsx_zz = dp(5)
2401  dpsx_rr = dp(3)
2402 
2403  psip_x=fun(1)+ dp(1)*(rx_p-rrx) + dp(2)*(zx_p-zzx)
2404  + + 0.5d0*dp(3)*(rx_p-rrx)*(rx_p-rrx)
2405  + + dp(4)*(rx_p-rrx)*(zx_p-zzx)
2406  + + 0.5d0*dp(5)*(zx_p-zzx)*(zx_p-zzx)
2407 
2408  return
2409  end
2410 
2411 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2412 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2413 
2414  !subroutine psi_cpn
2415  subroutine psi_cpn_wx
2416 
2417  include 'double.inc'
2418  parameter(nshp=10)
2419  include 'parcur.inc'
2420  include 'param.inc'
2421  include 'comblc.inc'
2422  real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
2423 
2424  do l=1,lfit
2425 
2426  r0=rfit(l)
2427  z0=zfit(l)
2428 
2429  ic=(r0-rmin)/dr(1)+1
2430  jc=(z0-zmin)/dz(1)+1
2431 
2432  psifit(l)=blin(ic,jc,r0,z0)
2433 
2434  enddo
2435 
2436  do l=1,lpre
2437 
2438  r0=rpre(l)
2439  z0=zpre(l)
2440 
2441  ic=(r0-rmin)/dr(1)+1
2442  jc=(z0-zmin)/dz(1)+1
2443 
2444  psipre(l)=blin(ic,jc,r0,z0)
2445 
2446  enddo
2447 
2448  r0=rx0
2449  z0=zx0
2450 
2451  ic=(r0-rmin)/dr(1)+1
2452  jc=(z0-zmin)/dz(1)+1
2453 
2454 
2455 c---definition of the nearest knote
2456 
2457  sdmin=rmax
2458 
2459  do 320 k=0,1
2460  rr=r(ic+k)
2461  do 325 l=0,1
2462  zz=z(jc+l)
2463  dlx=dsqrt( (rr-rx0)**2+(zz-zx0)**2 )
2464  if(dlx.lt.sdmin) then
2465  sdmin=dlx
2466  ik=ic+k
2467  jk=jc+l
2468  endif
2469  325 continue
2470  320 continue
2471 
2472  nsh=1
2473  xs(nsh)=r(ik)
2474  ys(nsh)=z(jk)
2475  fun(nsh)=ui(ik,jk)
2476 
2477  do 400 k=-1,1
2478 
2479  i= ik+k
2480 
2481  do 410 l=-1,1
2482 
2483  j= jk+l
2484 
2485  if(k.eq.0 .AND. l.eq.0 ) go to 410
2486  nsh=nsh+1
2487  xs(nsh)=r(i)
2488  ys(nsh)=z(j)
2489  fun(nsh)=ui(i,j)
2490 
2491  410 continue
2492  400 continue
2493 
2494  rrx=xs(1)
2495  zzx=ys(1)
2496 
2497  call deriv5(xs,ys,fun,nsh,5,dp)
2498 
2499  dpsx_r = dp(1) + dp(3)*(rx0-rrx) + dp(4)*(zx0-zzx)
2500  dpsx_z = dp(2) + dp(5)*(zx0-zzx) + dp(4)*(rx0-rrx)
2501 
2502  dpsx_rz = dp(4)
2503  dpsx_zz = dp(5)
2504  dpsx_rr = dp(3)
2505 
2506  psip_x=fun(1)+ dp(1)*(rx0-rrx) + dp(2)*(zx0-zzx)
2507  + + 0.5d0*dp(3)*(rx0-rrx)*(rx0-rrx)
2508  + + dp(4)*(rx0-rrx)*(zx0-zzx)
2509  + + 0.5d0*dp(5)*(zx0-zzx)*(zx0-zzx)
2510 
2511  return
2512  end
2513 
2514 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2515 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2516 
2517  subroutine psi_cpn
2518  !subroutine psi_cpn_3
2519 
2520  include 'double.inc'
2521 ! parameter(nshp=10)
2522  parameter(nshp=26)
2523  include 'parcur.inc'
2524  include 'param.inc'
2525  include 'comblc.inc'
2526 ! real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
2527  real*8 xs(nshp),ys(nshp),fun(nshp),dp(9)
2528 
2529  do l=1,lfit
2530 
2531  r0=rfit(l)
2532  z0=zfit(l)
2533 
2534  ic=(r0-rmin)/dr(1)+1
2535  jc=(z0-zmin)/dz(1)+1
2536 
2537  psifit(l)=blin(ic,jc,r0,z0)
2538 
2539  enddo
2540 
2541  do l=1,lpre
2542 
2543  r0=rpre(l)
2544  z0=zpre(l)
2545 
2546  ic=(r0-rmin)/dr(1)+1
2547  jc=(z0-zmin)/dz(1)+1
2548 
2549  psipre(l)=blin(ic,jc,r0,z0)
2550 
2551  enddo
2552 
2553  r0=rx_p
2554  z0=zx_p
2555 
2556  ic=(r0-rmin)/dr(1)+1
2557  jc=(z0-zmin)/dz(1)+1
2558 
2559 
2560 c---definition of the nearest node
2561 
2562  sdmin=rmax
2563 
2564  do 320 k=0,1
2565  rr=r(ic+k)
2566  do 325 l=0,1
2567  zz=z(jc+l)
2568  dlx=dsqrt( (rr-rx_p)**2+(zz-zx_p)**2 )
2569  if(dlx.lt.sdmin) then
2570  sdmin=dlx
2571  ik=ic+k
2572  jk=jc+l
2573  endif
2574  325 continue
2575  320 continue
2576 
2577  nsh=1
2578  xs(nsh)=r(ik)
2579  ys(nsh)=z(jk)
2580  fun(nsh)=ui(ik,jk)
2581 
2582 ! do 400 k=-1,1
2583  do 400 k=-2,2
2584 
2585  i= ik+k
2586 
2587 ! do 410 l=-1,1
2588  do 410 l=-2,2
2589 
2590  j= jk+l
2591 
2592  if(k.eq.0 .AND. l.eq.0 ) go to 410
2593  nsh=nsh+1
2594  xs(nsh)=r(i)
2595  ys(nsh)=z(j)
2596  fun(nsh)=ui(i,j)
2597 
2598  410 continue
2599  400 continue
2600 
2601  rrx=xs(1)
2602  zzx=ys(1)
2603 
2604 ! call deriv5(xs,ys,fun,nsh,5,dp)
2605  call deriv9(xs,ys,fun,nsh,9,dp)
2606 
2607  dpsx_r = dp(1) + dp(3)*(rx_p-rrx) + dp(4)*(zx_p-zzx)
2608  + +dp(6)*(rx_p-rrx)**2*3+dp(8)*(zx_p-zzx)**2
2609  + +dp(7)*(rx_p-rrx)*(zx_p-zzx)*2
2610  dpsx_z = dp(2) + dp(5)*(zx_p-zzx) + dp(4)*(rx_p-rrx)
2611  + +dp(7)*(rx_p-rrx)**2+dp(9)*(zx_p-zzx)**2*3
2612  + +dp(8)*(rx_p-rrx)*(zx_p-zzx)*2
2613  dpsx_rz = dp(4)
2614  + +dp(7)*(rx_p-rrx)*2+dp(8)*(zx_p-zzx)*2
2615  dpsx_zz = dp(5)
2616  + +dp(8)*(rx_p-rrx)*2+dp(9)*(zx_p-zzx)*6
2617  dpsx_rr = dp(3)
2618  + +dp(7)*(zx_p-zzx)*2+dp(6)*(rx_p-rrx)*6
2619 
2620  psip_x=fun(1)+ dp(1)*(rx_p-rrx) + dp(2)*(zx_p-zzx)
2621  + + 0.5d0*dp(3)*(rx_p-rrx)*(rx_p-rrx)
2622  + + dp(4)*(rx_p-rrx)*(zx_p-zzx)
2623  + + 0.5d0*dp(5)*(zx_p-zzx)*(zx_p-zzx)
2624  + + dp(6)*(rx_p-rrx)*(rx_p-rrx)*(rx_p-rrx)
2625  + + dp(7)*(rx_p-rrx)*(rx_p-rrx)*(zx_p-zzx)
2626  + + dp(8)*(rx_p-rrx)*(zx_p-zzx)*(zx_p-zzx)
2627  + + dp(9)*(zx_p-zzx)*(zx_p-zzx)*(zx_p-zzx)
2628 
2629  return
2630  end
2631 
2632 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2633 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2634 
2635  subroutine psi_cpn_2x
2636 
2637  include 'double.inc'
2638 ! parameter(nshp=10)
2639  parameter(nshp=26)
2640  include 'parcur.inc'
2641  include 'param.inc'
2642  include 'comblc.inc'
2643 ! real*8 xs(nshp),ys(nshp),fun(nshp),dp(5)
2644  real*8 xs(nshp),ys(nshp),fun(nshp),dp(9)
2645 
2646 
2647  r0=rx2_p
2648  z0=zx2_p
2649 
2650  ic=(r0-rmin)/dr(1)+1
2651  jc=(z0-zmin)/dz(1)+1
2652 
2653 
2654 c---definition of the nearest node
2655 
2656  sdmin=rmax
2657 
2658  do 320 k=0,1
2659  rr=r(ic+k)
2660  do 325 l=0,1
2661  zz=z(jc+l)
2662  dlx=dsqrt( (rr-rx2_p)**2+(zz-zx2_p)**2 )
2663  if(dlx.lt.sdmin) then
2664  sdmin=dlx
2665  ik=ic+k
2666  jk=jc+l
2667  endif
2668  325 continue
2669  320 continue
2670 
2671  nsh=1
2672  xs(nsh)=r(ik)
2673  ys(nsh)=z(jk)
2674  fun(nsh)=ui(ik,jk)
2675 
2676 ! do 400 k=-1,1
2677  do 400 k=-2,2
2678 
2679  i= ik+k
2680 
2681 ! do 410 l=-1,1
2682  do 410 l=-2,2
2683 
2684  j= jk+l
2685 
2686  if(k.eq.0 .AND. l.eq.0 ) go to 410
2687  nsh=nsh+1
2688  xs(nsh)=r(i)
2689  ys(nsh)=z(j)
2690  fun(nsh)=ui(i,j)
2691 
2692  410 continue
2693  400 continue
2694 
2695  rrx=xs(1)
2696  zzx=ys(1)
2697 
2698 ! call deriv5(xs,ys,fun,nsh,5,dp)
2699  call deriv9(xs,ys,fun,nsh,9,dp)
2700 
2701  dpsx2_r = dp(1) + dp(3)*(rx2_p-rrx) + dp(4)*(zx2_p-zzx)
2702  + +dp(6)*(rx2_p-rrx)**2*3+dp(8)*(zx2_p-zzx)**2
2703  + +dp(7)*(rx2_p-rrx)*(zx2_p-zzx)*2
2704  dpsx2_z = dp(2) + dp(5)*(zx2_p-zzx) + dp(4)*(rx2_p-rrx)
2705  + +dp(7)*(rx2_p-rrx)**2+dp(9)*(zx2_p-zzx)**2*3
2706  + +dp(8)*(rx2_p-rrx)*(zx2_p-zzx)*2
2707  dpsx2_rz = dp(4)
2708  + +dp(7)*(rx2_p-rrx)*2+dp(8)*(zx2_p-zzx)*2
2709  dpsx2_zz = dp(5)
2710  + +dp(8)*(rx2_p-rrx)*2+dp(9)*(zx2_p-zzx)*6
2711  dpsx2_rr = dp(3)
2712  + +dp(7)*(zx2_p-zzx)*2+dp(6)*(rx2_p-rrx)*6
2713 
2714  psip_x2= fun(1)+ dp(1)*(rx2_p-rrx) + dp(2)*(zx2_p-zzx)
2715  + + 0.5d0*dp(3)*(rx2_p-rrx)*(rx2_p-rrx)
2716  + + dp(4)*(rx2_p-rrx)*(zx2_p-zzx)
2717  + + 0.5d0*dp(5)*(zx2_p-zzx)*(zx2_p-zzx)
2718  + + dp(6)*(rx2_p-rrx)*(rx2_p-rrx)*(rx2_p-rrx)
2719  + + dp(7)*(rx2_p-rrx)*(rx2_p-rrx)*(zx2_p-zzx)
2720  + + dp(8)*(rx2_p-rrx)*(zx2_p-zzx)*(zx2_p-zzx)
2721  + + dp(9)*(zx2_p-zzx)*(zx2_p-zzx)*(zx2_p-zzx)
2722 
2723  return
2724  end
2725 
2726 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2727 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2728 
2729  subroutine cursol_(NEQUI,PFCEQW,psi_bnd)
2730 
2731  include 'double.inc'
2732  include 'parcur.inc'
2733 
2734  dimension a(ncf_p,ncf_p),x(ncf_p),y(ncf_p),ip(ncf_p)
2735  dimension psictr(ncf_p)
2736  dimension pfceqw(*)
2737 
2738 
2739 ! a(j,k)*I(k)+a(j,NEQUI+l)*Lam(l)+a(j,NEQUI+Lpre+1)*Psi_b=y(j)
2740 
2741  pi=3.14159265359d0
2742  amu0=0.4d0*pi
2743 
2744 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2745 ! equation (.)*d(I_j)=0
2746 !
2747 ! a(j,k)*I(k) ,j-number of equation
2748 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2749 
2750  do k=1,nequi
2751  do j=1,nequi
2752  asum=0.d0
2753  do l=1,lfit
2754  asum=asum+gindk(k,l)*gindk(j,l)*w_wght(l)
2755  enddo
2756  a(j,k)=asum
2757  enddo
2758  enddo
2759 
2760  do j=1,nequi
2761  a(j,j)=a(j,j)+d_wght(j)*s_wght
2762  enddo
2763 
2764 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
2765 
2766  do j=1,nequi
2767  do l=1,lpre
2768  a(j,nequi+l)=gindp(j,l)
2769  enddo
2770  enddo
2771 
2772 ! a(j,NEQUI+Lpre+1)*Psi_b ,j-number of equation
2773 
2774  do j=1,nequi
2775  asum=0.d0
2776  do l=1,lfit
2777  asum=asum+w_wght(l)*gindk(j,l)
2778  enddo
2779  a(j,nequi+lpre+1)=-asum
2780  enddo
2781 
2782  do j=1,nequi
2783  asum=0.d0
2784  do l=1,lfit
2785  asum=asum+w_wght(l)*gindk(j,l)*psifit(l)
2786  enddo
2787  y(j)=-asum+d_wght(j)*s_wght*curref(j)*amu0
2788  enddo
2789 
2790 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2791 !
2792 ! equation (.)*d(Lam_l)=0
2793 !
2794 ! a(j,k)*I(k) ,j-number of equation
2795 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2796 
2797  do l=1,lpre
2798  j=nequi+l
2799  do k=1,nequi
2800  a(j,k)=gindp(k,l)
2801  enddo
2802  enddo
2803 
2804 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
2805 
2806  do l=1,lpre
2807  j=nequi+l
2808  do ll=1,lpre
2809  k=nequi+ll
2810  a(j,k)=0.d0
2811  enddo
2812  enddo
2813 
2814 ! a(j,NEQUI+Lpre+1)*Psi_b ,j-number of equation
2815 
2816  do l=1,lpre
2817  j=nequi+l
2818  k=nequi+lpre+1
2819  a(j,k)=-1.d0
2820  enddo
2821 
2822  do l=1,lpre
2823  j=nequi+l
2824  y(j)=-psipre(l)
2825  enddo
2826 
2827 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2828 ! equation (.)*d(-psi_b)=0
2829 !
2830 ! a(j,k)*I(k) ,j-number of equation
2831 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2832 
2833  do k=1,nequi
2834  asum=0.d0
2835  do l=1,lfit
2836  asum=asum+w_wght(l)*gindk(k,l)
2837  enddo
2838  a(nequi+lpre+1,k)=asum
2839  enddo
2840 
2841 ! a(j,NEQUI+l)*Lam(l) ,j-number of equation
2842 
2843  do l=1,lpre
2844  k=nequi+l
2845  a(nequi+lpre+1,k)=1.d0
2846  enddo
2847 
2848 ! a(j,NEQUI+Lpre+1)*Psi_b ,j-number of equation
2849 
2850  asum=0.d0
2851  do l=1,lfit
2852  asum=asum+w_wght(l)
2853  enddo
2854  a(nequi+lpre+1,nequi+lpre+1)=-asum
2855 !!!
2856  asum=0.d0
2857  do l=1,lfit
2858  asum=asum+w_wght(l)*psifit(l)
2859  enddo
2860  y(nequi+lpre+1)=-asum
2861 
2862 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2863  call ge(nequi+lpre+1,ncf_p,a,y,x,ip)
2864 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2865 
2866  do iq=1,nequi
2867  pfceqw(iq)=x(iq)
2868  enddo
2869 
2870 
2871  do l=1,lfit
2872  psictr(l)=psifit(l)
2873  do iq=1,nequi
2874  psictr(l)=psictr(l)+gindk(iq,l)*pfceqw(iq)
2875  enddo
2876  enddo
2877 
2878  do l=1,lpre
2879  psictr(lfit+l)=psipre(l)
2880  do iq=1,nequi
2881  psictr(lfit+l)=psictr(lfit+l)+gindp(iq,l)*pfceqw(iq)
2882  enddo
2883  enddo
2884 
2885  do iq=1,nequi
2886  pfceqw(iq)=x(iq)/amu0
2887  enddo
2888 
2889  psi_bnd=x(nequi+lpre+1)
2890 
2891  return
2892  end
2893 
2894 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2895 
2896  subroutine bonpsi
2897 
2898  include 'double.inc'
2899  include 'dim.inc'
2900  include 'parcur.inc'
2901  include 'compol.inc'
2902 
2903  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
2904  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
2905 
2906  dimension binadg(ntp,ntp),dgdn(ntp)
2907 
2908  i=iplas !!!!!!!!!!!!!!!!!!!!!!1
2909 
2910  do 10 j=2,nt1
2911 
2912  a1=a13(i-1,j-1)
2913  a2=a34(i-1,j-1)+a12(i-1,j)
2914  a3=a24(i-1,j)
2915 
2916  g1=psi(i-1,j-1)
2917  g2=psi(i-1,j)
2918  g3=psi(i-1,j+1)
2919 
2920  dltk=(dlt(i,j-1)+dlt(i,j))*0.5d0
2921  dgdnl=a1*g1+a2*g2+a3*g3
2922  dgdn(j)=dgdnl/dltk
2923 
2924  10 continue
2925 
2926  dgdn(1)=dgdn(nt1)
2927  dgdn(nt)=dgdn(2)
2928 
2929 
2930  do l=1,lfit
2931 
2932  rr=rfit(l)
2933  zz=zfit(l)
2934 
2935  do jb=1,nt1
2936 
2937  r0=r(iplas,jb)
2938  z0=z(iplas,jb)
2939 
2940  r1=r(iplas,jb+1)
2941  z1=z(iplas,jb+1)
2942 
2943  call bint(rr,zz,r0,z0,r1,z1,fint,1)
2944 
2945  binadg(jb,l)=fint
2946 
2947  enddo
2948  enddo
2949 
2950  do l=1,lfit
2951 
2952  psb=0.d0
2953 
2954  do jb=2,nt1
2955 
2956  psb=psb+binadg(jb,l)*(dgdn(jb)+dgdn(jb+1))*0.5d0
2957 
2958  enddo
2959 
2960  psifit(l)=-psb
2961 
2962  enddo
2963 
2964  do l=1,lpre
2965 
2966  rr=rpre(l)
2967  zz=zpre(l)
2968 
2969  do jb=1,nt1
2970 
2971  r0=r(iplas,jb)
2972  z0=z(iplas,jb)
2973 
2974  r1=r(iplas,jb+1)
2975  z1=z(iplas,jb+1)
2976 
2977  call bint(rr,zz,r0,z0,r1,z1,fint,1)
2978 
2979  binadg(jb,l)=fint
2980 
2981  enddo
2982  enddo
2983 
2984  do l=1,lpre
2985 
2986  psb=0.d0
2987 
2988  do jb=2,nt1
2989 
2990  psb=psb+binadg(jb,l)*(dgdn(jb)+dgdn(jb+1))*0.5d0
2991 
2992  enddo
2993 
2994  psipre(l)=-psb
2995 
2996  enddo
2997 
2998  return
2999  end
3000 
3001 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
greeni
REAL *8 function greeni(R, Z, RP, ZP)
Definition: scu.f:484
fun
subroutine fun(X, F)
Definition: Ev2.f:10
deriv9
subroutine deriv9(X, Y, F, M, N, U)
Definition: scu.f:338
bnd_modul
Definition: bnd_modul.f90:1
bint
subroutine bint(X, Y, R0, Z0, r1, z1, F, I)
Definition: scu.f:557
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
ge
subroutine ge(N, NZ, A, X, Y, IP)
Definition: scu.f:420
curfit
subroutine curfit(iopt, m, x, y, w, xb, xe, k, s, nest, n, t, c, fp, wrk, lwrk, iwrk, ier)
Definition: curfit.f:1
grgren
subroutine grgren(R0, Z0, R, Z, dGdr, dGdz)
Definition: scu.f:519
blin
real *8 function blin(i, j, r0, z0)
Definition: EQ1_m.f:835
d2gren
subroutine d2gren(R0, Z0, R, Z, d2Gdrz, d2Gdzz, d2Gdrr)
Definition: scu.f:696