ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
B_metric_O.f
Go to the documentation of this file.
1 !!!this file contains the following routines:
2 !!!
3 !!!!!!!! metric
4 !!!!!!!! geom
5 !!!!!!!! funsq
6 !!!!!!!! numlin
7 !!!!!!!! matcof
8 !!!!!!!! matpla
9 !!!!!!!! numlin
10 !!!!!!!! prog1d
11 !!!!!!!! procof
12 !!!!!!!! extrpc
13 !!!!!!!!!!!!!!!!
14 
15 ch4astra !subroutine metric
16  subroutine metrix
17 
18  !-----------------
19 
20  include 'double.inc'
21  include 'dim.inc'
22  include 'compol.inc'
23 
24  sqrt(xx)=dsqrt(xx)
25 
26  !write(6,*) 'metric:enter'
27 
28  do 10 j=1,nt
29  do 10 i=1,iplas1
30 
31  dlr(i,j)=sqrt( (r(i+1,j)-r(i,j))**2+(z(i+1,j)-z(i,j))**2 )
32  st(i,j)=dlr(i,j)*(r(i+1,j)+r(i,j))*0.5d0
33 
34  10 continue
35  !write(6,*) 'metric:10 '
36 
37  do 20 j=1,nt1
38  do 20 i=2,iplas
39 
40  dlt(i,j)=sqrt( (r(i,j+1)-r(i,j))**2+(z(i,j+1)-z(i,j))**2 )
41  sr(i,j)=dlt(i,j)*(r(i,j+1)+r(i,j))*0.5d0
42 
43  20 continue
44 
45  !write(6,*) 'metric:20 '
46 
47 
48  do 35 j=1,nt1
49  do 30 i=1,iplas1
50 
51  r1=r(i,j)
52  r2=r(i+1,j)
53  r3=r(i+1,j+1)
54  r4=r(i,j+1)
55  r0=(r1+r2+r3+r4)*0.25d0
56 
57  r12=(r1+r2)*0.5d0
58  r23=(r3+r2)*0.5d0
59  r34=(r3+r4)*0.5d0
60  r14=(r1+r4)*0.5d0
61 
62  z1=z(i,j)
63  z2=z(i+1,j)
64  z3=z(i+1,j+1)
65  z4=z(i,j+1)
66  z0=(z1+z2+z3+z4)*0.25d0
67 
68  z12=(z1+z2)*0.5d0
69  z23=(z3+z2)*0.5d0
70  z34=(z3+z4)*0.5d0
71  z14=(z1+z4)*0.5d0
72 
73  dl12=dlr(i,j)
74  dl14=dlt(i,j)
75  dl23=dlt(i+1,j)
76  dl34=dlr(i,j+1)
77 
78  sq_1=funsq(r1,r12,r0,r14,z1,z12,z0,z14)
79  sq_2=funsq(r2,r23,r0,r12,z2,z23,z0,z12)
80  sq_3=funsq(r3,r34,r0,r23,z3,z34,z0,z23)
81  sq_4=funsq(r4,r14,r0,r34,z4,z14,z0,z34)
82 
83  !S(i,j)=funsq(r1,r2,r3,r4,z1,z2,z3,z4)
84 
85  !sq1(i,j)=0.25d0*S(i,j)
86  !sq2(i,j)=0.25d0*S(i,j)
87  !sq3(i,j)=0.25d0*S(i,j)
88  !sq4(i,j)=0.25d0*S(i,j)
89 
90  !vol1(i,j)=funsq(r1,r2,r4,r4,z1,z2,z4,z4)*(r1+r2+r4)/6.d0
91  !vol2(i,j)=funsq(r1,r2,r3,r3,z1,z2,z3,z3)*(r1+r2+r3)/6.d0
92  !vol3(i,j)=funsq(r2,r3,r4,r4,z2,z3,z4,z4)*(r2+r3+r4)/6.d0
93  !vol4(i,j)=funsq(r1,r3,r4,r4,z1,z3,z4,z4)*(r1+r3+r4)/6.d0
94 
95  !vol1(i,j)=funsq(r1,r2,r4,r4,z1,z2,z4,z4)*r0*0.5d0
96  !vol2(i,j)=funsq(r1,r2,r3,r3,z1,z2,z3,z3)*r0*0.5d0
97  !vol3(i,j)=funsq(r2,r3,r4,r4,z2,z3,z4,z4)*r0*0.5d0
98  !vol4(i,j)=funsq(r1,r3,r4,r4,z1,z3,z4,z4)*r0*0.5d0
99 
100  !vol(i,j)=(vol1(i,j)+vol2(i,j)+vol3(i,j)+vol4(i,j))
101 
102  !vol1(i,j)=0.25d0*vol(i,j)
103  !vol2(i,j)=0.25d0*vol(i,j)
104  !vol3(i,j)=0.25d0*vol(i,j)
105  !vol4(i,j)=0.25d0*vol(i,j)
106 
107  dr12=r2-r1
108  dz12=z2-z1
109  !dl12= (dr12**2+dz12**2)
110 
111  dr14=r4-r1
112  dz14=z4-z1
113  !dl14= (dr14**2+dz14**2)
114 
115  dr23=r3-r2
116  dz23=z3-z2
117  !dl23= (dr23**2+dz23**2)
118 
119  dr34=r3-r4
120  dz34=z3-z4
121  !dl34= (dr34**2+dz34**2)
122 
123 
124  dlm2=dl12*dl23
125  cos_2=(dr12*dr23+dz12*dz23)/dlm2
126  sin_2=1.d0-cos_2**2
127  vol_2=dlm2*sqrt(sin_2)*(r1+r2+r4)/12.d0
128 
129  dlm3=dl34*dl23
130  cos_3=(dr34*dr23+dz34*dz23)/dlm3
131  sin_3=1.d0-cos_3**2
132  vol_3=dlm3*sqrt(sin_3)*(r2+r3+r4)/12.d0
133 
134  if(i.ne.1) then
135 
136  dlm4=dl34*dl14
137  cos_4=(dr34*dr14+dz34*dz14)/dlm4
138  sin_4=1.d0-cos_4**2
139  vol_4=dlm4*sqrt(sin_4)*(r1+r3+r4)/12.d0
140 
141  dlm1=dl12*dl14
142  cos_1=(dr12*dr14+dz12*dz14)/dlm1
143  sin_1=1.d0-cos_1**2
144  vol_1=dlm1*sqrt(sin_1)*(r1+r2+r4)/12.d0
145 
146  else
147  !cos1(i,j)=0.d0
148  cos_1=0.d0 ! DPC
149  sin_1=1.d0
150  vol_1=0.d0
151  !
152  !cos2(i,j)=0.d0
153  !sin2(i,j)=1.d0
154  !
155  !cos3(i,j)=0.d0
156  !sin3(i,j)=1.d0
157  !
158  !cos4(i,j)=0.d0
159  cos_4=0.d0 ! DPC
160  sin_4=1.d0
161  vol_4=0.d0
162  endif
163 
164  cos2(i,j)=cos_2
165  sin2(i,j)=sin_2
166  vol2(i,j)=vol_2
167 
168  cos3(i,j)=cos_3
169  sin3(i,j)=sin_3
170  vol3(i,j)=vol_3
171 
172  cos4(i,j)=cos_4
173  sin4(i,j)=sin_4
174  vol4(i,j)=vol_4
175 
176  cos1(i,j)=cos_1
177  sin1(i,j)=sin_1
178  vol1(i,j)=vol_1
179 
180  vol(i,j)=vol_1+vol_2+vol_3+vol_4
181 
182  sq1(i,j)=sq_1
183  sq2(i,j)=sq_2
184  sq3(i,j)=sq_3
185  sq4(i,j)=sq_4
186 
187  s(i,j)=sq_1+sq_2+sq_3+sq_4
188 
189  30 continue
190  35 continue
191 
192  return
193  end
194 
195 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
196 
197  subroutine geom_b
198 !!!!cylindr geometry
199 
200  include 'double.inc'
201  include 'dim.inc'
202  include 'compol.inc'
203 
204  sqrt(xx)=dsqrt(xx)
205  !write(6,*) 'metric:enter'
206 
207  do 10 i=1,nr1
208  do 10 j=1,nt
209 
210  dlr(i,j)=sqrt( (r(i+1,j)-r(i,j))**2+(z(i+1,j)-z(i,j))**2 )
211  st(i,j)=dlr(i,j)*rm
212 
213  10 continue
214 
215  !write(6,*) 'metric:10 '
216 
217  do 20 i=2,nr
218  do 20 j=1,nt1
219 
220  dlt(i,j)=sqrt( (r(i,j+1)-r(i,j))**2+(z(i,j+1)-z(i,j))**2 )
221  sr(i,j)=dlt(i,j)*rm
222 
223  20 continue
224 
225  !write(6,*) 'metric:20 '
226 
227  do 30 i=1,nr1
228  do 35 j=1,nt1
229 
230  r1=r(i,j)
231  r2=r(i+1,j)
232  r3=r(i+1,j+1)
233  r4=r(i,j+1)
234  r0=(r1+r2+r3+r4)*0.25d0
235 
236  r12=(r1+r2)*0.5d0
237  r23=(r3+r2)*0.5d0
238  r34=(r3+r4)*0.5d0
239  r14=(r1+r4)*0.5d0
240 
241  z1=z(i,j)
242  z2=z(i+1,j)
243  z3=z(i+1,j+1)
244  z4=z(i,j+1)
245  z0=(z1+z2+z3+z4)*0.25d0
246 
247  z12=(z1+z2)*0.5d0
248  z23=(z3+z2)*0.5d0
249  z34=(z3+z4)*0.5d0
250  z14=(z1+z4)*0.5d0
251 
252  sq1(i,j)=funsq(r1,r12,r0,r14,z1,z12,z0,z14)
253  sq2(i,j)=funsq(r2,r23,r0,r12,z2,z23,z0,z12)
254  sq3(i,j)=funsq(r3,r34,r0,r23,z3,z34,z0,z23)
255  sq4(i,j)=funsq(r4,r14,r0,r34,z4,z14,z0,z34)
256 
257  s(i,j)=funsq(r1,r2,r3,r4,z1,z2,z3,z4)
258 
259  vol1(i,j)=funsq(r1,r2,r4,r4,z1,z2,z4,z4)*(rm+rm+rm)/6.d0
260  vol2(i,j)=funsq(r1,r2,r3,r3,z1,z2,z3,z3)*(rm+rm+rm)/6.d0
261  vol3(i,j)=funsq(r2,r3,r4,r4,z2,z3,z4,z4)*(rm+rm+rm)/6.d0
262  vol4(i,j)=funsq(r1,r3,r4,r4,z1,z3,z4,z4)*(rm+rm+rm)/6.d0
263 
264  vol(i,j)=(vol1(i,j)+vol2(i,j)+vol3(i,j)+vol4(i,j))
265 
266  dr12=r2-r1
267  dz12=z2-z1
268 
269  dl12= (dr12**2+dz12**2)
270 
271  dr14=r4-r1
272  dz14=z4-z1
273  dl14= (dr14**2+dz14**2)
274 
275  dr23=r3-r2
276  dz23=z3-z2
277  dl23= (dr23**2+dz23**2)
278 
279  dr34=r3-r4
280  dz34=z3-z4
281  dl34= (dr34**2+dz34**2)
282 
283 
284  cos2(i,j)=(dr12*dr23+dz12*dz23)/sqrt(dl12*dl23)
285  sin2(i,j)=1.-cos2(i,j)**2
286 
287  cos3(i,j)=(dr34*dr23+dz34*dz23)/sqrt(dl34*dl23)
288  sin3(i,j)=1.-cos3(i,j)**2
289 
290  if(i.eq.1) go to 35
291 
292  cos4(i,j)=(dr34*dr14+dz34*dz14)/sqrt(dl34*dl14)
293  sin4(i,j)=1.-cos4(i,j)**2
294 
295  cos1(i,j)=(dr12*dr14+dz12*dz14)/sqrt(dl12*dl14)
296  sin1(i,j)=1.-cos1(i,j)**2
297 
298 
299 
300  cos1(i,j)=0.d0
301  sin1(i,j)=1.d0
302 
303  cos2(i,j)=0.d0
304  sin2(i,j)=1.d0
305 
306  cos3(i,j)=0.d0
307  sin3(i,j)=1.d0
308 
309  cos4(i,j)=0.d0
310  sin4(i,j)=1.d0
311 
312  35 continue
313 
314  30 continue
315 
316  return
317  end
318 
319 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
320 !
321 ! REAL FUNCTION FUNSQ*8(R1,R2,R3,R4,Z1,Z2,Z3,Z4)
322 ! include 'double.inc'
323 !
324 ! R13=R3-R1
325 ! R24=R4-R2
326 ! Z13=Z3-Z1
327 ! Z24=Z4-Z2
328 !
329 ! ZS=R13*Z24-R24*Z13
330 !
331 ! IF(ZS.LT.0.) WRITE(6,*) '***ERROR:SQUARE<0',zs
332 !
333 ! FUNSQ=ZS*0.5D0
334 !
335 ! RETURN
336 ! END
337 !
338 
339  subroutine matrot
340 
341  include 'double.inc'
342  include 'dim.inc'
343  include 'compol.inc'
344  common/com_rot/ ar1(nrp,ntp),ar2(nrp,ntp),ar3(nrp,ntp),
345  + ar4(nrp,ntp),ar5(nrp,ntp),ar6(nrp,ntp)
346  dimension b_tet(nrp,ntp),b_ro(nrp,ntp),rot_b(nrp,ntp)
347  dimension z_err(nrp,ntp)
348 
349  do i=2,iplas1
350  do j=2,nt1
351 
352 !!!! B_tet(i-1,j-1)
353 
354  ar1(i,j)= cos2(i-1,j-1)*vol2(i-1,j-1)/(sr(i,j-1)*sin2(i-1,j-1))
355 
356 !!!! B_tet(i-1,j)
357 
358  ar2(i,j)=( -cos2(i-1,j)/sr(i,j)-1.d0/st(i-1,j) )
359  & *vol2(i-1,j)/sin2(i-1,j) +
360  & ( -1.d0/st(i-1,j) )*vol1(i-1,j)/sin1(i-1,j)+
361  & ( cos3(i-1,j-1)/sr(i,j-1)-1.d0/st(i-1,j) )
362  & *vol3(i-1,j-1)/sin3(i-1,j-1) +
363  & ( -1.d0/st(i-1,j) )*vol4(i-1,j-1)/sin4(i-1,j-1)
364 
365 !!!! B_tet(i-1,j+1)
366 
367  ar3(i,j)=-cos3(i-1,j)*vol3(i-1,j)/(sr(i,j)*sin3(i-1,j))
368 
369 !!!! B_tet(i,j-1)
370 
371  ar4(i,j)= cos1(i,j-1)*vol1(i,j-1)/(sr(i,j-1)*sin1(i,j-1))
372 
373 !!!! B_tet(i,j)
374 
375  ar5(i,j)=( -cos1(i,j)/sr(i,j)+1.d0/st(i,j) )
376  & *vol1(i,j)/sin1(i,j) +
377  & ( 1.d0/st(i,j) )*vol2(i,j)/sin2(i,j)+
378  & ( cos4(i,j-1)/sr(i,j-1)+1.d0/st(i,j) )
379  & *vol4(i,j-1)/sin4(i,j-1) +
380  & ( 1.d0/st(i,j) )*vol3(i,j-1)/sin3(i,j-1)
381 
382 
383 !!!! B_tet(i,j+1)
384 
385  ar6(i,j)=-cos4(i,j)*vol4(i,j)/(sr(i,j)*sin4(i,j))
386 
387 
388  enddo
389  enddo
390 
391 
392  do i=2,iplas1
393  do j=1,nt1
394 
395  b_ro(i,j)=(psi(i,j+1)-psi(i,j))/sr(i,j)
396 
397  enddo
398  enddo
399 
400  do i=1,iplas1
401  do j=1,nt
402 
403  b_tet(i,j)=(psi(i+1,j)-psi(i,j))/st(i,j)
404 
405  enddo
406  enddo
407 
408  do i=2,iplas1
409  do j=2,nt1
410 
411  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
412 
413  bt1=b_tet(i-1,j-1)
414  bt2=b_tet(i-1,j)
415  bt3=b_tet(i-1,j+1)
416  bt4=b_tet(i,j-1)
417  bt5=b_tet(i,j)
418  bt6=b_tet(i,j+1)
419 
420  rot_b(i,j)=-( ar1(i,j)*bt1+ar2(i,j)*bt2+ar3(i,j)*bt3+
421  & ar4(i,j)*bt4+ar5(i,j)*bt5+ar6(i,j)*bt6 )/sqk
422 
423  z_err(i,j)=rot_b(i,j)-cur(i,j)
424 
425  enddo
426  enddo
427 
428  return
429  end
430 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
431 
432  subroutine matcof
433 
434  include 'double.inc'
435  include 'dim.inc'
436  include 'compol.inc'
437 
438  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
439  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
440 
441 
442  do 10 j=1,nt1
443  do 10 i=1,iplas1
444 
445  s12=st(i,j)
446  s14=sr(i,j)
447  s34=st(i,j+1)
448  s23=sr(i+1,j)
449 
450  if(i.ne.1) then
451 
452  a12(i,j)=( (-1.d0/s12 + cos1(i,j)/s14)*vol1(i,j)/sin1(i,j) +
453  + (-1.d0/s12 - cos2(i,j)/s23)*vol2(i,j)/sin2(i,j) )/s12
454 
455  a23(i,j)=( (-1.d0/s23 - cos2(i,j)/s12)*vol2(i,j)/sin2(i,j) +
456  + (-1.d0/s23 + cos3(i,j)/s34)*vol3(i,j)/sin3(i,j) )/s23
457 
458  a34(i,j)=( (-1.d0/s34 + cos3(i,j)/s23)*vol3(i,j)/sin3(i,j) +
459  + (-1.d0/s34 - cos4(i,j)/s14)*vol4(i,j)/sin4(i,j) )/s34
460 
461  a14(i,j)=( (-1.d0/s14 + cos1(i,j)/s12)*vol1(i,j)/sin1(i,j) +
462  + (-1.d0/s14 - cos4(i,j)/s34)*vol4(i,j)/sin4(i,j) )/s14
463 
464  a13(i,j)= (cos2(i,j)/(s12*s23))*vol2(i,j)/sin2(i,j) +
465  + (cos4(i,j)/(s34*s14))*vol4(i,j)/sin4(i,j)
466 
467  a24(i,j)= (-cos1(i,j)/(s12*s14))*vol1(i,j)/sin1(i,j) +
468  + (-cos3(i,j)/(s34*s23))*vol3(i,j)/sin3(i,j)
469 
470  else
471 
472  a12(i,j)=( (-1.d0/s12 - cos2(i,j)/s23)*vol2(i,j)/sin2(i,j) )/s12
473 
474  a23(i,j)=( (-1.d0/s23 - cos2(i,j)/s12)*vol2(i,j)/sin2(i,j) +
475  + (-1.d0/s23 + cos3(i,j)/s34)*vol3(i,j)/sin3(i,j) )/s23
476 
477  a34(i,j)=( (-1.d0/s34 + cos3(i,j)/s23)*vol3(i,j)/sin3(i,j) )/s34
478 
479  a13(i,j)= (cos2(i,j)/(s12*s23))*vol2(i,j)/sin2(i,j)
480 
481  a24(i,j)= (-cos3(i,j)/(s34*s23))*vol3(i,j)/sin3(i,j)
482 
483  endif
484 
485  10 continue
486 
487  return
488  end
489 
490 
491  subroutine procof(icq,cur_mu)
492 
493  include 'double.inc'
494  include 'dim.inc'
495  include 'compol.inc'
496 
497  if(ngav.eq.0) cur_mu=tok*amu0
498  call psib_pla(pspl_av)
499 
500  if(ngav.eq.1) then
501  call procof_i(icq)
502  elseif(ngav.eq.2) then
503  call procof_fl(icq)
504  elseif(ngav.eq.-10 .AND. erru.lt.5.d-3) then
505  call promat_j(iplas,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psim,iter,
506  * kstep,fvac,psi_eav,pspl_av,psibon0,cur_mu )
507  elseif(ngav.eq.-1 .AND. erru.lt.5.d-3) then
508  call promat(iplas,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psim,iter,
509  * kstep,fvac,psi_eav,pspl_av,psibon0,cur_mu )
510  elseif(ngav.eq.-2 .AND. erru.lt.5.d-3) then
511  call promat(iplas,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psim,iter,
512  * kstep,fvac,psi_eav,pspl_av,psibon0,cur_mu )
513  elseif(ngav.eq.-3 .AND. erru.lt.5.d-3) then
514  cur_mu=tok*amu0
515  call promat_i(iplas,dfdpsi,dpdpsi,f,flx_fi,psia,q,rm,psim,iter,
516  * kstep,fvac,psi_eav,pspl_av,psibon0,cur_mu )
517  endif
518 
519  return
520  end
521 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
522 
523  subroutine procof_fl(icq)
524 
525 !!!flux conserving
526 
527  include 'double.inc'
528  include 'dim.inc'
529  include 'compol.inc'
530 
531  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
532  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
533 
534  common/compsf/ psf(nrp),sqtor(nrp)
535 
536  dimension amn(nrp),a0(nrp),apl(nrp),capp(nrp)
537  dimension avrc(nrp),delsc(nrp),avrk(nrp),delsk(nrp),delv(nrp)
538  dimension bmn(nrp),b0(nrp),bpl(nrp),wrk1(nrp),wrk2(nrp)
539  dimension rhs(nrp),delv3(nrp)
540  dimension dfdpsn(nrp),delf(nrp)
541  dimension qs(nrp),qsn(nrp),psfn(nrp),dpsdas(nrp)
542 
543  sqrt(xx)=dsqrt(xx)
544 
545  ! write(6,*)' procof_fl:enter'
546  ! write(6,*)'psiax,,psip',psiax,psip
547 
548  do i=1,iplas
549 
550  dfdpsn(i)=dfdpsi(i)
551  if(iter.eq.1) dfdpsn(i)=0.d0
552  qs(i)=q(i)
553  psfn(i)=psia(i)
554  dpsdas(i)=dpsda(i)
555 
556  enddo
557 
558  964 continue
559 
560  do i=1,iplas
561 
562  qsn(i)=qs(i)
563 
564  enddo
565 
566  do 30 i=1,iplas1
567 
568  zdelsc=0.d0
569  zavrc=0.d0
570 
571  do 35 j=2,nt1
572 
573  r1=r(i,j)
574  r2=r(i+1,j)
575  r3=r(i+1,j+1)
576  r4=r(i,j+1)
577  r0=(r1+r2+r3+r4)*0.25d0
578 
579  zdelsc=zdelsc+s(i,j)
580  zavrc=zavrc+s(i,j)/r0
581 
582  35 continue
583 
584  avrc(i)=zavrc/zdelsc
585  delsc(i)=zdelsc
586 
587  30 continue
588 
589  sqtor(1)=0.d0
590 
591  do 33 i=2,iplas
592  sqtor(i)=sqtor(i-1)+delsc(i-1)
593  33 continue
594 
595  do 40 i=2,iplas
596 
597  zdelv=0.d0
598  zdelv3=0.d0
599  zdelsk=0.d0
600  zavrk=0.d0
601  zapro=0.d0
602 
603  do 45 j=2,nt1
604 
605  if(i.ne.iplas) then
606  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
607  r0=r(i,j)
608  else
609  sqk=sq2(i-1,j)+sq3(i-1,j-1)
610  r0=(3.d0*r(i,j)+r(i-1,j))*0.25d0
611  endif
612 
613  zdelsk=zdelsk+sqk
614  zdelv=zdelv +sqk*r0
615  zdelv3=zdelv3 +sqk*r0**3
616  zavrk=zavrk+sqk/r0
617  zapro=zapro
618 
619  45 continue
620 
621  avrk(i)=zavrk/zdelsk
622  delsk(i)=zdelsk
623  delv(i)=zdelv
624  delv3(i)=zdelv3
625  rhs(i-1)=0.d0
626 
627  40 continue
628 
629  do 10 i=2,iplas1
630 
631  samn=0.d0
632  sa0 =0.d0
633  sapl=0.d0
634 
635  do 110 j=2,nt1
636 
637  a1=a13(i-1,j-1)
638  a2=a34(i-1,j-1)+a12(i-1,j)
639  a3=a24(i-1,j)
640  a4=a23(i-1,j-1)+a14(i,j-1)
641  a6=a14(i,j)+a23(i-1,j)
642  a7=a24(i,j-1)
643  a8=a34(i,j-1)+a12(i,j)
644  a9=a13(i,j)
645  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
646 
647  zamn=a1+a2+a3
648  za0 =a4+a5+a6
649  zapl=a7+a8+a9
650 
651  samn=samn+zamn
652  sa0 =sa0 +za0
653  sapl=sapl+zapl
654 
655  110 continue
656 
657  qmn=-qs(i-1)/(avrc(i-1)*delsc(i-1))
658  qpl=-qs(i) / (avrc(i)*delsc(i))
659  qk=(qmn*dpsdas(i-1)+qpl*dpsdas(i))/(dpsdas(i-1)+dpsdas(i))
660 
661  afmn=qk*qmn*avrk(i)*delsk(i)
662  afpl=qk*qpl*avrk(i)*delsk(i)
663 
664  af0=-(afmn+afpl)
665 
666  amn(i)=samn - afmn
667  a0(i) =sa0 - af0
668  apl(i)=sapl - afpl
669 
670  capp(i)=sapl*delsc(i)
671 
672  if(i.eq.iplas-1) then
673 
674  cappl=sapl*delsc(i)
675  capmn=samn*delsc(i-1)
676 
677  endif
678 
679  10 continue
680 
681  npro=iplas-2
682 
683  do 100 i=1,npro
684 
685  if(i.ne.1) then
686  bmn(i)=amn(i+1)
687  else
688  bmn(i)=0.d0
689  endif
690 
691  b0(i)=-a0(i+1)
692 
693  if(i.ne.npro) then
694  bpl(i)=apl(i+1)
695  else
696  bpl(i)=0.d0
697  endif
698 
699  rhs(i)=rhs(i)-dpdpsi(i+1)*delv(i+1)-dwdpsi(i+1)*delv3(i+1)
700 
701  !write(6,*) i
702  !write(6,*) bmn(i),b0(i),bpl(i)
703 
704  100 continue
705 
706  rhs(1)=rhs(1)+amn(2)*psiax
707  rhs(npro)=rhs(npro)+apl(iplas-1)*psip
708 
709  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
710 
711  call prog1d(npro,psf(2),bpl,bmn,b0,rhs,wrk1,wrk2)
712 
713  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
714 
715  psf(1)=psiax
716  psf(iplas)=psip
717 
718  psfn(1)=1.d0
719  psfn(iplas)=0.d0
720 
721  do i=2,iplas1
722 
723  psfn(i)=(psf(i)-psip)/(psiax-psip)
724 
725  enddo
726 
727  do 300 i=1,iplas1
728 
729  f(i)=-qs(i)*(psf(i+1)-psf(i))/(delsc(i)*avrc(i))
730 
731  300 continue
732 
733  do 400 i=2,iplas1
734 
735  qmn=-qs(i-1)/(avrc(i-1)*delsc(i-1))
736  qpl=-qs(i)/(avrc(i)*delsc(i))
737  qk=(qmn*dpsdas(i-1)+qpl*dpsdas(i))/(dpsdas(i-1)+dpsdas(i))
738 
739  dfdpsi(i)=qk*(f(i)-f(i-1))
740 
741  400 continue
742 
743 !!!!!!!!!!///////////////////////////////////////////
744 
745 !!! equation in central node
746 
747  ac0=0.d0
748  skcen=0.d0
749 
750  do 20 j=2,nt1
751 
752  acj=a12(1,j)+a24(1,j)+a34(1,j-1)+a13(1,j-1)
753  ac0=ac0-acj
754 
755  skcen=skcen+sq1(1,j)+sq4(1,j-1)
756 
757  20 continue
758 
759  dfdpsi(1)=
760  * ( (ac0*(psiax-psf(2)))/skcen-rm*dpdpsi(1)-dwdpsi(1)*rm**3 )*rm
761 
762 !!! equation on bound node!
763 !!!!!!!!!!!!!!!!!!!!!!!!!!!
764 
765  i=iplas
766 
767  samn=0.d0
768  sa0 =0.d0
769 
770  do 50 j=2,nt1
771 
772  a1=a13(i-1,j-1)
773  a2=a34(i-1,j-1)+a12(i-1,j)
774  a3=a24(i-1,j)
775 
776  zamn=a1+a2+a3
777  za0 =-zamn
778  samn=samn+zamn
779  sa0 =sa0 +za0
780 
781  50 continue
782 
783  !cappa=cappl+delsk(i)*(cappl-capmn)/delsk(i-1)
784 
785  call extrpc( sqtor(i),capp(i),
786  * sqtor(i),sqtor(i-1),sqtor(i-2),sqtor(i-3),
787  * capp(i-1) ,capp(i-2) ,capp(i-3) )
788 
789 ! call extrpc( psia(i),capp(i),
790 ! * psia(i),psia(i-1),psia(i-2),psia(i-3),
791 ! * capp(i-1) ,capp(i-2) ,capp(i-3) )
792 
793  cappa=capp(i)
794 
795  !!!!!!!!!!!!!!!!!!!!!!!!
796 
797  if(icq.eq.0) then
798 
799  dpsids=(
800  * samn*(psip-psf(i-1))+avrk(i)*delsk(i)*dfdpsi(i)
801  * +delv(i)*dpdpsi(i)+delv3(i)*dwdpsi(i)
802  * ) /cappa
803 
804  f(i)=sqrt( f(i-1)**2+dfdpsi(i)*(psip-psf(i-1)) )
805 
806  q(i)=-f(i)*avrk(i)/dpsids
807 
808  write(6,*) 'q(iplas)=',q(i)
809 
810 ! pause ' '
811 
812  endif
813 
814  !q(i)=q(i-1)+delsk(i)*(q(i-1)-q(i-2))/delsk(i-1)
815 
816  !!!!!!!!!!!!!!!!!!!!!!!!
817 
818  qn=-q(i)/avrk(i)
819  qmn=-q(i-1)/avrc(i-1)
820 
821  !Q14= 0.5d0*(Qn+Qmn)
822 
823  q14=-0.5d0*(q(i)+q(i-1))
824 
825  avr14=avrk(i)
826 
827 ! f(i)=(samn*(psip-psf(i-1))+delv(i)*dpdpsi(i)-Q14*avr14*f(i-1))
828 ! * /( cappa/Qn-Q14*avr14 )
829 
830 
831  f(i)=
832  * ( samn*(psip-psf(i-1))
833  * +delv(i)*dpdpsi(i)+delv3(i)*dwdpsi(i)-q14*f(i-1) )
834  * /( cappa/qn-q14 )
835 
836 ! dfdpsi(i)=Q14*(f(i)-f(i-1))/delsk(i)
837  dfdpsi(i)=(q14/avr14)*(f(i)-f(i-1))/delsk(i)
838 
839  psiai=psia(i)-0.25d0*(psia(i)- psia(i-1))
840 
841  call extrp2( psiai,dfdpsi(i),
842  * psia(i-3),psia(i-2),psia(i-1),
843  * dfdpsi(i-3) ,dfdpsi(i-2) ,dfdpsi(i-1) )
844 
845 
846  f(i)=dsqrt(f(i-1)**2+dfdpsi(i)*(psip- psf(i-1)))
847 
848 !!!!!!!!!!!!!>>>>>>>>
849 
850  !dfdpsi(i)=0.d0
851 
852 !!!!!!!!!!!!!>>>>>>>>
853 
854  if(itin.le.2) then
855 
856  vrh=1.00d0
857 
858  else
859 
860  vrh=1.00d0
861 
862  endif
863 
864  do i=1,iplas
865 
866  dfdpsi(i)=vrh*dfdpsi(i)+(1.d0-vrh)*dfdpsn(i)
867 
868  enddo
869 
870  do i=1,iplas1
871 
872  delf(i)=dfdpsi(i)-dfdpsn(i)
873 
874  enddo
875  delf(iplas)=0.d0
876 
877  !if(itin.eq.1) wgt=1.0d0
878 
879  ! do ki=1,2
880  do i=2,iplas1
881 
882  wgt=1.0d0
883 
884  dmon=(delf(i+1)-delf(i))/(delf(i)-delf(i-1))
885 
886  !if(dmon.lt.-0.5d0.and.dmon.gt.-2.0d0) wgt=0.5d0
887  if(dmon.lt.0.d0) wgt=0.5d0
888 
889  delf(i)=wgt*delf(i)+(1.d0-wgt)*
890  * (delf(i-1)*(psf(i+1)-psf(i))+delf(i+1)*(psf(i)-psf(i-1)))/
891  * (psf(i+1)-psf(i-1))
892 
893  enddo
894  ! enddo
895 
896 
897  do i=2,iplas1
898 
899  dfdpsi(i)=dfdpsn(i)+delf(i)
900 
901  enddo
902 
903 
904  !if(icq.ne.0) then
905 
906 
907  ! if(ngav.eq.2) then
908  !open(1,file='ddps1.pr')
909  !do 567 i=1,iplas
910  !567 !write(1,*) dfdpsi(i),f(i),i
911  ! write(1,*) 'q(iplas)=',q(iplas)
912  ! close(1)
913  ! stop
914  ! endif
915 
916 
917  return
918  end
919 
920 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
921 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
922 
923  subroutine procof_i(icq)
924 
925  include 'double.inc'
926  include 'dim.inc'
927  include 'compol.inc'
928 
929  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
930  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
931 
932  common/compsf/ psf(nrp),sqtor(nrp)
933 
934  dimension amn(nrp),a0(nrp),apl(nrp),qcapp(nrp)
935  dimension avrc(nrp),delsc(nrp),avrk(nrp),delsk(nrp),delv(nrp)
936  dimension bmn(nrp),b0(nrp),bpl(nrp),wrk1(nrp),wrk2(nrp)
937  dimension rhs(nrp),delv3(nrp)
938  dimension dfdpsn(nrp),vw(nrp)
939  dimension dfdpsw(nrp),delf(nrp)
940 
941  sqrt(xx)=dsqrt(xx)
942  ! write(6,*) 'procof_I:enter'
943  do i=1,iplas
944 
945  dfdpsw(i)=dfdpsi(i)
946  dfdpsn(i)=dfdpsi(i)
947  psf(i)=psia(i)*(psim-psip)+psip
948 
949  enddo
950 
951 c!!!!!!! dS(i+1/2) and <1/R>(i+1/2) !!!!!!!!!!!
952 
953  do 30 i=1,iplas1
954 
955  zdelsc=0.d0
956  zavrc=0.d0
957 
958  do 35 j=2,nt1
959 
960  r1=r(i,j)
961  r2=r(i+1,j)
962  r3=r(i+1,j+1)
963  r4=r(i,j+1)
964 
965  r0=(r1+r2+r3+r4)*0.25d0
966 
967  zdelsc=zdelsc+s(i,j)
968  zavrc=zavrc+s(i,j)/r0
969 
970  35 continue
971 
972  avrc(i)=zavrc/zdelsc
973  delsc(i)=zdelsc
974 
975  30 continue
976 
977 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
978 
979  sqtor(1)=0.d0
980 
981  do 33 i=2,iplas
982 
983  sqtor(i)=sqtor(i-1)+delsc(i-1)
984 
985  33 continue
986 
987 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
988 
989 c!!!!!!! dS(i), dV(i) and <1/R>(i) !!!!!!!!!!!
990 
991  do 40 i=2,iplas
992 
993  zdelv=0.d0
994  zdelv3=0.d0
995  zdelsk=0.d0
996  zavrk=0.d0
997 
998  do 45 j=2,nt1
999 
1000  if(i.ne.iplas) then
1001 
1002  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
1003  r0=r(i,j)
1004 
1005  else
1006 
1007  sqk=sq2(i-1,j)+sq3(i-1,j-1)
1008  r0=(3.d0*r(i,j)+r(i-1,j))*0.25d0
1009 
1010  endif
1011 
1012  zdelsk=zdelsk+sqk
1013  zdelv=zdelv +sqk*r0
1014  zdelv3=zdelv3 +sqk*r0**3
1015  zavrk=zavrk+sqk/r0
1016 
1017  45 continue
1018 
1019  avrk(i)=zavrk/zdelsk
1020  delsk(i)=zdelsk
1021  delv(i)=zdelv
1022  delv3(i)=zdelv3
1023  rhs(i-1)=0.d0 !right hand side for averaged equation
1024 
1025  40 continue
1026 
1027  do 10 i=2,iplas1
1028 
1029  samn=0.d0
1030  sa0 =0.d0
1031  sapl=0.d0
1032 
1033  do 110 j=2,nt1
1034 
1035  a1=a13(i-1,j-1)
1036  a2=a34(i-1,j-1)+a12(i-1,j)
1037  a3=a24(i-1,j)
1038  a4=a23(i-1,j-1)+a14(i,j-1)
1039  a6=a14(i,j)+a23(i-1,j)
1040  a7=a24(i,j-1)
1041  a8=a34(i,j-1)+a12(i,j)
1042  a9=a13(i,j)
1043  a5=-(a1+a2+a3+a4+a6+a7+a8+a9)
1044 
1045  zamn=a1+a2+a3
1046  za0 =a4+a5+a6
1047  zapl=a7+a8+a9
1048 
1049  samn=samn+zamn
1050  sa0 =sa0 +za0
1051  sapl=sapl+zapl
1052 
1053  110 continue
1054 
1055  qmn=-q(i-1)/(avrc(i-1)*delsc(i-1))
1056  qpl=-q(i) / (avrc(i)*delsc(i))
1057  qk=(qmn*dpsda(i-1)+qpl*dpsda(i))/(dpsda(i-1)+dpsda(i))
1058 
1059  afmn=qk*qmn*avrk(i)*delsk(i)
1060  afpl=qk*qpl*avrk(i)*delsk(i)
1061 
1062  af0=-(afmn+afpl)
1063 
1064  amn(i)=samn !- afmn
1065  a0(i) =sa0 !- af0
1066  apl(i)=sapl !- afpl
1067 
1068  qcapp(i)=qpl/sapl
1069 
1070  if(i.eq.iplas-1) then
1071 
1072  qcappl=qpl/sapl
1073  qcapmn=qmn/samn
1074 
1075  endif
1076 
1077  10 continue
1078 
1079  if(itin.eq.1) then
1080 
1081  errdf=0.d0
1082 
1083  do i=2,iplas1
1084 
1085  !plcu=apl(i)*(psf(i+1)-psf(i))
1086  !write(6,*) 'tok',plcu,i
1087  !pause ' pause'
1088 
1089  zff=amn(i)*psf(i-1)+a0(i)*psf(i)+apl(i)*psf(i+1)
1090 
1091  zdfdps=( zff-dpdpsi(i)*delv(i)-dwdpsi(i)*delv3(i) )
1092  * /(avrk(i)*delsk(i))
1093 
1094  !znvz=zff-dpdpsi(i)*delv(i)-dfdpsn(i)*avrk(i)*delsk(i)
1095 
1096  dfdpsi(i)=zdfdps
1097  deldf=dabs(zdfdps-dfdpsn(i))
1098 
1099  !write(6,*) 'deldf',deldf,dfdpsn(i),i
1100  !write(6,*) 'nev.L(psi)',znvz,zff
1101  !pause ' pause'
1102 
1103  errdf=dmax1(deldf,errdf)
1104  dfdpsn(i)=dfdpsi(i)
1105  dfdpsi(i)=0.d0 !!!!!!!!!!!!!<<<<<<<<<<<<
1106 
1107  enddo
1108 
1109  !write(6,*) 'errdf',errdf
1110  !pause ' pause'
1111  endif
1112  !! return
1113 
1114  !!!!!!!!!!
1115 ! i=iplas
1116 ! !!!!!!!!!!
1117 !
1118 ! zff=tok-apl(i-1)*(psf(i)-psf(i-1))
1119 !
1120 ! dfdpsi(i)=(zff-dpdpsi(i)*delv(i))/(avrk(i)*delsk(i))
1121 ! deldf=dabs(dfdpsi(i)-dfdpsn(i))
1122 ! !write(6,*) 'deldf*',deldf,i
1123 ! !write(6,*) '*',dfdpsi(i),dfdpsn(i)
1124 ! !pause ' pause'
1125 !
1126 !
1127 !
1128 ! Qcappa=Qcappl+delsk(i)*(Qcappl-Qcapmn)/delsk(i-1)
1129 !
1130 ! call extrpc( sqtor(i),Qcapp(i),
1131 ! * sqtor(i),sqtor(i-1), sqtor(i-2), sqtor(i-3),
1132 ! * Qcapp(i-1) ,Qcapp(i-2) ,Qcapp(i-3) )
1133 !
1134 ! call extrpc( psia(i),Qcapp(i),
1135 ! * psia(i),psia(i-1), psia(i-2), psia(i-3),
1136 ! * Qcapp(i-1) ,Qcapp(i-2) ,Qcapp(i-3) )
1137 !
1138 ! !ff2n=f(i-1)**2+dfdpsn(i)*(psf(i)-psf(i-1))
1139 !
1140 ! !f2n=(tok*Qcappa)**2
1141 ! f2n=(tok*Qcapp(i))**2
1142 !
1143 ! !write(6,*) 'fnold fn',ff2n,f2n
1144 !
1145 ! f(i)=dsqrt(f2n)
1146 ! f(i-1)=dsqrt( f2n-dfdpsi(i)*(psf(i)-psf(i-1)) )
1147 !
1148  !!amu0=0.4d0*pi
1149 
1150  736 continue
1151 
1152  errdf=0.d0
1153 
1154  f2n1=(amu0*tok*qcapp(iplas-1))**2
1155  f(iplas-1)=dsqrt(f2n1)
1156 
1157  do i=iplas-1,2,-1
1158 
1159  f(i-1)=dsqrt( f(i)**2-dfdpsi(i)*(psf(i+1)-psf(i-1)) )
1160 
1161  enddo
1162 
1163  psf(iplas)=0.d0
1164 
1165  do i=iplas1,1,-1
1166 
1167  dftor=f(i)*delsc(i)*avrc(i)
1168  delpsi=dftor/q(i)
1169  psf(i)=psf(i+1)+delpsi
1170 
1171  enddo
1172 
1173  do i=2,iplas1
1174 
1175  !plcu=apl(i)*(psf(i+1)-psf(i))
1176  !write(6,*) 'tok',plcu,i
1177  !pause ' pause'
1178 
1179  zff=amn(i)*psf(i-1)+a0(i)*psf(i)+apl(i)*psf(i+1)
1180 
1181  dfdpsi(i)=( zff-dpdpsi(i)*delv(i)-dwdpsi(i)*delv3(i) )
1182  * /(avrk(i)*delsk(i))
1183  deldf=dabs(dfdpsi(i)-dfdpsn(i))
1184 
1185  !write(6,*) 'deldf',deldf,errdf,i
1186  !pause ' pause'
1187 
1188  errdf=dmax1(deldf,errdf)
1189  dfdpsn(i)=dfdpsi(i)
1190 
1191  enddo
1192 
1193  !write(6,*) 'errdf',errdf
1194  !pause ' pause'
1195 
1196  if(errdf.gt.1.d-7) go to 736
1197 
1198  i=iplas !!!<<<<<<<<<<<<<
1199 
1200  call extrp2( psia(i),dfdpsi(i),
1201  * psia(i-3),psia(i-2),psia(i-1),
1202  * dfdpsi(i-3) ,dfdpsi(i-2) ,dfdpsi(i-1) )
1203 
1204 !!!!!!!!!!///////////////////////////////////////////
1205 
1206 !!! equation in central node
1207 
1208  ac0=0.d0
1209  skcen=0.d0
1210 
1211  do 20 j=2,nt1
1212 
1213  acj=a12(1,j)+a24(1,j)+a34(1,j-1)+a13(1,j-1)
1214  ac0=ac0-acj
1215 
1216  skcen=skcen+sq1(1,j)+sq4(1,j-1)
1217 
1218  20 continue
1219 
1220  dfdpsi(1)=( ac0*(psf(1)-psf(2))/skcen
1221  * -rm*dpdpsi(1)-dwdpsi(1)*rm**3 )*rm
1222 
1223 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1224  do i=1,iplas1
1225 
1226  delf(i)=dfdpsi(i)-dfdpsw(i)
1227 
1228  enddo
1229  delf(iplas)=0.d0
1230 
1231  !if(itin.eq.1) wgt=1.0d0
1232 
1233  ! do ki=1,2
1234  do i=2,iplas1
1235 
1236  wgt=1.0d0
1237 
1238  dmon=(delf(i+1)-delf(i))/(delf(i)-delf(i-1))
1239 
1240  !if(dmon.lt.-0.5d0.and.dmon.gt.-2.0d0) wgt=0.5d0
1241  if(dmon.lt.0.d0) wgt=0.5d0
1242 
1243  delf(i)=wgt*delf(i)+(1.d0-wgt)*
1244  * (delf(i-1)*(psf(i+1)-psf(i))+delf(i+1)*(psf(i)-psf(i-1)))/
1245  * (psf(i+1)-psf(i-1))
1246 
1247  enddo
1248  ! enddo
1249 
1250 
1251  do i=2,iplas1
1252 
1253  dfdpsi(i)=dfdpsw(i)+delf(i)
1254 
1255  enddo
1256 
1257 
1258 
1259 
1260  ! write(6,*) 'procof:exit'
1261 
1262  return
1263  end
1264 
1265 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1266 
1267 
1268 
1269 
1270 
1271 
1272 
1273 
1274 
1275  subroutine metcof_
1276 
1277  include 'double.inc'
1278  include 'dim.inc'
1279  include 'compol.inc'
1280 
1281 
1282  dimension ax(nrp),gr_a(nrp)
1283  dimension delsc(nrp),delv(nrp),dvda(nrp)
1284 
1285 
1286 
1287  sqrt(xx)=dsqrt(xx)
1288 
1289  do i=1,iplas
1290  ax(i)=(i-1.d0)/(iplas-1.d0)
1291  enddo
1292 
1293 
1294  i=1
1295 
1296  zdelv=0.d0
1297  sss=0.d0
1298  sss0=0.d0
1299 
1300  dadr=0.d0
1301  dadz=0.d0
1302 
1303  u0=ax(1)
1304 
1305  do j=2,nt1
1306 
1307  u1=ax(2)
1308  u2=ax(2)
1309 
1310  r1=r(2,j)
1311  r2=r(2,j+1)
1312 
1313  z1=z(2,j)
1314  z2=z(2,j+1)
1315 
1316  dadr=dadr+0.5d0*(u1+u2)*(z2-z1)
1317  dadz=dadz-0.5d0*(u1+u2)*(r2-r1)
1318 
1319  sss=sss+funsq(r1,r2,0.d0,0.d0,z1,z2,0.d0,0.d0)
1320  sqk=sq1(i,j)+sq4(i,j-1)
1321  zdelv=zdelv +sqk*r(1,2)
1322 
1323  enddo
1324 
1325  dadr=dadr/sss
1326  dadz=dadz/sss
1327 
1328  gr_a(1)=dadr**2+dadz**2
1329  delv(1)=zdelv
1330 
1331 
1332  do i=1,iplas
1333 
1334  zdelv=0.d0
1335  zavgr=0.d0
1336 
1337  do j=2,nt1
1338 
1339  r1=r(i-1,j)
1340  r2=r(i,j-1)
1341  r3=r(i+1,j)
1342  r4=r(i,j+1)
1343 
1344  z1=z(i-1,j)
1345  z2=z(i,j-1)
1346  z3=z(i+1,j)
1347  z4=z(i,j+1)
1348 
1349  u1=ax(i-1)
1350  u2=ax(i)
1351  u3=ax(i+1)
1352  u4=ax(i)
1353 
1354  if(i.eq.iplas) then
1355 
1356  r3=r(i,j)
1357  z3=z(i,j)
1358  u3=ax(i)
1359 
1360  endif
1361 
1362  sss=funsq(r1,r2,r3,r4,z1,z2,z3,z4)
1363 
1364  dadr= 0.5d0*((u1+u2)*(z2-z1)+(u2+u3)*(z3-z2)+
1365  * (u3+u4)*(z4-z3)+(u4+u1)*(z1-z4))/sss
1366 
1367  dadz=-0.5d0*((u1+u2)*(r2-r1)+(u2+u3)*(r3-r2)+
1368  * (u3+u4)*(r4-r3)+(u4+u1)*(r1-r4))/sss
1369 
1370  gr2=dadr**2+dadz**2
1371 
1372  if(i.ne.iplas) then
1373  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
1374  r0=r(i,j)
1375  else
1376  sqk=sq2(i-1,j)+sq3(i-1,j-1)
1377  r0=(3.d0*r(i,j)+r(i-1,j))*0.25d0
1378  endif
1379 
1380  zdelv=zdelv +sqk*r0
1381  zavgr=zavgr+gr2*sqk*r0
1382 
1383  enddo
1384  gr_a(i)=zavgr/zdelv
1385  delv(i)=zdelv
1386  enddo
1387 
1388  do i=1,iplas-1
1389  dvda(i)=(delv(i+1)-delv(i))/(ax(i+1)-ax(i))
1390  enddo
1391 
1392  return
1393  end
1394 
1395 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1396 
1397 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1398 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1399 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1400 
1401  subroutine qunew(qs,psfn)
1402 
1403  include 'double.inc'
1404  include 'dim.inc'
1405  include 'compol.inc'
1406 
1407  dimension qs(1),psfn(1)
1408 
1409  qs(1)=q(1)
1410  qs(iplas)=q(iplas)
1411 
1412  do i=2,iplas1
1413 
1414  psx=0.5d0*(psfn(i)+psfn(i+1))
1415 
1416  do is=1,iplas1
1417 
1418  q0=q(is)
1419  qpl=q(is+1)
1420 
1421  ps0=0.5d0*(psia(is)+psia(is+1))
1422 
1423  if(is.ne.iplas1) then
1424 
1425  pspl=0.5d0*(psia(is+1)+psia(is+2))
1426 
1427  else
1428 
1429  pspl=1.d0
1430 
1431  endif
1432 
1433  if(psx.le.ps0 .AND. psx.ge.pspl) then
1434 
1435  qs(i)=(q0*(pspl-psx)+qpl*(psx-ps0))/(pspl-ps0)
1436 
1437  go to 10
1438 
1439  endif
1440 
1441  enddo
1442 
1443  10 continue
1444 
1445  enddo
1446 
1447  return
1448  end
1449 
1450 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1451 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1452 
1453 !**********************************************************************
1454 ! PROGONKA - USUAL
1455 !
1456 ! A(J)*U(J+1)-C(J)*U(J)+B(J)*U(J-1)=-F(J) , J=1,M
1457 !
1458 ! B(1)=0 , A(M)=0
1459 !
1460 !**********************************************************************
1461 
1462 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1463 
1464 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
extrpc
subroutine extrpc(X, F, X0, X1, X2, X3, F1, F2, F3)
Definition: com_sub.f:699
metcof_
subroutine metcof_
Definition: B_metric_O.f:1275
matrot
subroutine matrot
Definition: B_metric_O.f:339
funsq
REAL *8 function funsq(R1, R2, R3, R4, Z1, Z2, Z3, Z4)
Definition: scu.f:803
qunew
subroutine qunew(qs, psfn)
Definition: B_metric_O.f:1401
procof
subroutine procof(icq, cur_mu)
Definition: B_metric_O.f:491
metrix
subroutine metrix
Definition: B_metric_O.f:16
prog1d
subroutine prog1d(M, U, A, B, C, F, ALF, BET)
Definition: scu.f:822
promat
subroutine promat(n, dfdpsi, dpdpsi, f, flx_fi, psia, q, rm, psim, iter, kstep, fvac, psi_eav, pspl_av, psibon0, cur_mu)
Definition: exp_promat_e_I_f0.f:1
helena21_mod::dpdpsi
real(r8) function dpdpsi(psi_n)
Definition: helena21_mod.f90:122
procof_fl
subroutine procof_fl(icq)
Definition: B_metric_O.f:523
extrp2
subroutine extrp2(X, F, X1, X2, X3, F1, F2, F3)
Definition: com_sub.f:719
matcof
subroutine matcof
Definition: B_metric_O.f:432
psib_pla
subroutine psib_pla(pspl_av)
Definition: B_rig.f:499
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
promat_i
subroutine promat_i(n, dfdpsi, dpdpsi, f, flx_fi, psia, q, rm, psim, iter, kstep, fvac, psi_eav, pspl_av, psibon0, tok)
Definition: exp_promat_e_I_f0.f:428
promat_j
subroutine promat_j(n, dfdpsi, dpdpsi, f, flx_fi, psia, q, rm, psim, iter, kstep, fvac, psi_eav, pspl_av, psibon0, cur_mu)
Definition: exp_promat_e_I_f0.f:1088
procof_i
subroutine procof_i(icq)
Definition: B_metric_O.f:923
geom_b
subroutine geom_b
Definition: B_metric_O.f:197