ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
B_MESH.f
Go to the documentation of this file.
1 
2  subroutine qst(qcen,cnor,b0ax,r0ax)
3 
4  include 'double.inc'
5  ! include 'param.inc'
6 
7  parameter(nshp=10,ntp=128,nrp=64)
8  ! include 'comblc.inc'
9 
10  common /complr/ rplr(nrp,ntp),zplr(nrp,ntp),dp(5),
11  * psia(nrp),psi(nrp),psim,psip,rm,zm,
12  * fplr(nrp),qplr(nrp),dflufi(nrp),fluxfi(nrp),
13  * iplas,nt,iplas1,nt1
14 
15  common /comhel/ helinp, helout
16 
17 c*************************************************************************
18 
19 c write(6,*) '"qst" => enter'
20 
21  pi=3.14159265358d0
22 
23  drr=dp(3)
24  drz=dp(4)
25  dzz=dp(5)
26 
27  tg2a=2.d0*drz/(drr-dzz)
28  cos2a=1.d0/sqrt(1.d0+tg2a**2)
29  sin2a=cos2a*tg2a
30 
31  dxx=0.5d0*(drr+dzz)+0.5d0*cos2a*(drr-dzz)+sin2a*drz
32  dyy=0.5d0*(drr+dzz)+0.5d0*cos2a*(dzz-drr)-sin2a*drz
33 
34  1 continue
35 
36  bfcen=qcen*dsqrt(dxx*dyy)
37 
38  fcen=bfcen*rm
39  fvac=b0ax*r0ax
40 
41 c write(6,*) 'fcen, qcen = ', fcen, qcen
42 c pause 'pause: fcen, qcen'
43 
44  dpsi=(psia(iplas)-psia(iplas-1))
45  ps14=-0.25d0*dpsi
46  ffp=tabf(ps14)
47  !ffp=dfdpsi(iplas)
48  fplr(iplas-1)=sqrt(fvac**2-ffp*cnor*(psim-psip)*dpsi)
49 
50  do 200 i=iplas-2,1,-1
51  dpsi=0.5d0*(psia(i+2)-psia(i))*(psim-psip)
52  ffp=tabf(psia(i+1))
53  fplr(i)=sqrt(fplr(i+1)**2-2.d0*ffp*dpsi)
54  200 continue
55 
56  fplr(iplas)=fvac
57 
58 c---------------------------------------
59 c open(1,file='q_rec.pr')
60 
61 c write(1,*) '----------------------'
62 
63  do 16 i=1,iplas1
64 
65  dflufi(i)=0.
66 
67  do 17 j=2,nt1
68 
69  r1=rplr(i,j)
70  r2=rplr(i+1,j)
71  r3=rplr(i+1,j+1)
72  r4=rplr(i,j+1)
73 
74  z1=zplr(i,j)
75  z2=zplr(i+1,j)
76  z3=zplr(i+1,j+1)
77  z4=zplr(i,j+1)
78 
79  r0=(r1+r2+r3+r4)*0.25d0
80 
81  sqpol=funsq(r1,r2,r3,r4,z1,z2,z3,z4)
82 
83  dflufi(i)=dflufi(i)+sqpol*fplr(i)/r0
84 
85  17 continue
86 
87  q2pi=-dflufi(i)/((psia(i+1)-psia(i))*(psim-psip))
88 
89  qplr(i)=0.5d0*q2pi/pi
90 !!!!!! q(i)=q2pi
91 
92 c write(1,*) 'qplr(i), dflufi(i), i ', qplr(i), dflufi(i), i
93 c write(6,*) 'qplr(i), dflufi(i), i ', qplr(i), dflufi(i), i
94 
95 coment if( dflufi(i).lt.0.d0 ) pause 'pause: dflufi(i) < 0'
96 
97  16 continue
98 c...................................................
99 
100  fluxfi(1)=0.d0
101 
102  do i=2,iplas
103  fluxfi(i)=fluxfi(i-1) + dflufi(i-1)
104  enddo
105 
106 c write(1,*) '----------------------'
107 c write(1,*) 'dflufi(i), i=1,iplas =', iplas
108 c write(1,*) (dflufi(i), i=1,iplas)
109 c write(1,*) '----------------------'
110 c write(1,*) '----------------------'
111 c write(1,*) 'qplr(i), i=1,iplas =', iplas
112 c write(1,*) (qplr(i), i=1,iplas)
113 c write(1,*) '----------------------'
114 c write(1,*) '----------------------'
115 c write(1,*) 'fplr(i), i=1,iplas =', iplas
116 c write(1,*) (fplr(i), i=1,iplas)
117 c write(1,*) '----------------------'
118 c write(1,*) '----------------------'
119 c write(1,*) 'fluxfi(i), i=1,iplas =', iplas
120 c write(1,*) (fluxfi(i), i=1,iplas)
121 c write(1,*) '----------------------'
122 c write(1,*) 'psim, psip =', psim, psip
123 c write(1,*) '----------------------'
124 c write(1,*) 'psi(i), i=1,iplas =', iplas
125 c write(1,*) (psi(i), i=1,iplas)
126 c write(1,*) '----------------------'
127 c write(1,*) '----------------------'
128 c write(1,*) 'psia(i), i=1,iplas =', iplas
129 c write(1,*) (psia(i), i=1,iplas)
130 c write(1,*) '----------------------'
131 c..................................................
132 
133  helout = - fluxfi(iplas)*psip
134 
135  do i=1,iplas1
136  helout = helout + 0.5d0*( psi(i)+psi(i+1) )*dflufi(i)
137  enddo
138 
139 c write(1,*) '----------------------'
140 c write(6,*) 'helinp, helout =', helinp, helout
141 c write(1,*) 'helinp, helout =', helinp, helout
142 c write(1,*) '----------------------'
143 c--------------------------------------------------
144 c close(1)
145 c--------------------------------------------------
146 
147 c write(6,*) '"qst" => exit'
148 
149  !!!call PET_DIN()
150 
151  return
152  end
153 
154 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
155 
156  subroutine gridpl
157 
158  include 'double.inc'
159  include 'param.inc'
160  parameter(nshp=10,ntp=128,nrp=64)
161  include 'comblc.inc'
162 
163  common /complr/ rplr(nrp,ntp),zplr(nrp,ntp),dp(5),
164  * psia(nrp),psi(nrp),psim,psip,xm,ym,
165  * fplr(nrp),qplr(nrp),dflufi(nrp),fluxfi(nrp),
166  * iplas,nt,iplas1,nt1
167 
168  real*8 ron(ntp),teta(ntp)
169  real*8 xs(nshp),ys(nshp),fun(nshp)
170 
171  cos(xx) =dcos(xx)
172  sin(xx) =dsin(xx)
173  sqrt(xx)=dsqrt(xx)
174 
175 c*******************************************************************
176 
177 c write(6,*) '"gridpl" => enter'
178 
179  iplas = 61
180  nt = 95
181  iplas1 = iplas-1
182  nt1 = nt-1
183 
184  ixp1=0
185  jxp1=0
186 
187  ixp2=0
188  jxp2=0
189 
190  psim=um
191  psip=up
192 
193  xm=rm
194  ym=zm
195 
196  dtet=2.d0*pi/(nt-2)
197 
198  teta(1)=-dtet
199 
200  do j=2,nt
201  teta(j) = teta(j-1) + dtet
202  enddo
203 
204  teta(1)=teta(nt1)
205  teta(nt)=teta(2)
206 
207  psia(1)=1.d0
208 
209 
210  nsh = 1
211  xs(nsh) = r(imax)
212  ys(nsh) = z(jmax)
213  fun(nsh) = u(imax,jmax)
214 
215  do k=-1,1
216  ii = imax + k
217  do l=-1,1
218 
219  jj = jmax + l
220 
221  if(ii.ne.imax .OR. jj.ne.jmax) then
222  nsh = nsh+1
223  xs(nsh) = r(ii)
224  ys(nsh) = z(jj)
225  fun(nsh) = u(ii,jj)
226  endif
227 
228  enddo
229  enddo
230 
231  call deriv5(xs,ys,fun,nsh,5,dp)
232 
233  ! write(6,*) dp
234 
235  stpx=(r(imax+1)-r(imax))*1.0
236  stpy=(z(jmax+1)-z(jmax))*0.7
237 
238 c----------------------------------------
239  do i=2,iplas
240 
241  u0=1.d0-((i-1)/(iplas-1.)) ! **2
242  psia(i)=u0
243  ur0=up+u0*(um-up)
244 
245  do j=1,nt
246 
247  ttj=teta(j)
248  ro2j=(ur0-um)/( 0.5*dp(3)*cos(ttj)**2
249  + + dp(4)*cos(ttj)*sin(ttj)
250  + + 0.5*dp(5)*sin(ttj)**2 )
251 
252  ron(j)=sqrt(ro2j)
253 
254 
255  rplr(i,j)=rm+ron(j)*cos(teta(j))
256  zplr(i,j)=zm+ron(j)*sin(teta(j))
257 
258 
259  enddo ! j-loop
260 
261  ibeg=i+1
262 
263  if( ron(2) .GT. dmax1(stpx,stpy) ) go to 2799
264 
265  enddo !i-loop
266 
267 c------------------------
268  2799 continue
269 c------------------------
270 
271 c write(6,*)'ibeg',ibeg,u0
272 c pause 'pause: ibeg'
273 
274  do 100 i=ibeg,iplas
275 
276  u0=(iplas-i)/(iplas-1.d0)
277  ! u0=1.d0-((i-1)/(iplas-1.))**2
278 
279  psia(i)=u0
280 
281  call loop95(teta,nt,ron,u0)
282 
283  do 110 j=1,nt
284  rplr(i,j)=rm+ron(j)*cos(teta(j))
285  zplr(i,j)=zm+ron(j)*sin(teta(j))
286  110 continue
287 
288  100 continue
289 c-------------------------
290 
291  do 120 j=1,nt
292  rplr(1,j) = rm
293  zplr(1,j) = zm
294  120 continue
295 
296 
297  do 333 i=1,iplas
298  rplr(i,1)=rplr(i,nt1)
299  zplr(i,1)=zplr(i,nt1)
300 
301  rplr(i,nt)=rplr(i,2)
302  zplr(i,nt)=zplr(i,2)
303  333 continue
304 
305  teta(1)=teta(nt1)-2.d0*pi
306  teta(nt)=teta(2)+2.d0*pi
307 
308  do i=1,iplas
309  do j=1,nt
310  psi(i)=psip+psia(i)*(psim-psip)
311  enddo
312  enddo
313 
314 c write(6,*) '"gridpl" => exit'
315 
316 ! open(1,file='map_bnd.wr')
317 ! write(1,*) nt-2
318 ! do j=2,nt1
319 ! write(1,*) rplr(iplas,j),zplr(iplas,j)
320 ! enddo
321 ! close(1)
322 
323  return
324  end
325 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
326 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
327 
328  subroutine loop95(tetpol,ntet,ro0,u0)
329 
330  include 'double.inc'
331  include 'param.inc'
332  include 'comblc.inc'
333 
334  dimension rxb(nbndp2),zxb(nbndp2)
335 
336  real*8 ut(nip,njp)
337 
338  real*8 roxb(nbndp2),tetxb(nbndp2)
339 
340  real*8 rrk(nbndp4),cck(nbndp4),wrk(nbndp6)
341 
342  real*8 cwk(4),tetpol(1),ro0(1)
343 
344  xzer(x1,x2,v1,v2) = (x1*v2-x2*v1)/(v2-v1)
345 
346  ! write(6,*) '***loop95:enter'
347  ! write(6,*) '***loop95:imax,jmax',imax,jmax
348  ! write(6,*) '***loop95:rmax,zmax',rm,zm
349 
350  if(u0.lt.2.d-5) u0=2.d-5
351 
352  imax1=imax-1
353  jmax1=jmax-1
354 
355  do 10 i=imax1,imax
356  if(rm.le.r(i+1) .AND. rm.gt.r(i)) icell=i
357  10 continue
358 
359  do 20 j=jmax1,jmax
360  if(zm.le.z(j+1) .AND. zm.gt.z(j)) jcell=j
361  20 continue
362 
363  ! write(6,*) 'loop95:icell,jcell',icell,jcell
364 
365  i=icell
366  j=jcell+1
367 
368 cccc u0=0.5
369  do 15 ii=1,ni
370  do 15 jj=1,nj
371  ut(ii,jj)=un(ii,jj)-u0
372  15 continue
373 
374  ig=1
375 
376  888 continue
377 
378  do 885 i=imax,ni1
379 
380  ! write(6,*) 'un(i,j)',un(i,j),un(i+1,j)
381 
382  if(ut(i,j)*ut(i+1,j).le.0.) then
383 
384  ic=i
385  jc=j
386 
387  rxb(ig)=xzer(r(i+1),r(i),ut(i+1,j),ut(i,j))
388  zxb(ig)=z(j)
389 
390  go to 886
391  endif
392 
393  885 continue
394  ! write(6,*) 'loop95: first point was not find'
395  stop
396 
397  886 continue
398 
399  ic1=ic
400  jc1=jc
401 
402 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
403 
404  lin=1
405 
406  100 ig=ig+1
407 
408  i=ic
409  j=jc
410 
411  if(ut(i+1,j)*ut(i,j).le.0. .AND. lin.ne.1) then
412 
413  ic=i
414  jc=j-1
415 
416  lin=3
417 
418  rxb(ig)=xzer(r(i+1),r(i),ut(i+1,j),ut(i,j))
419  zxb(ig)=z(j)
420 
421  elseif(ut(i+1,j+1)*ut(i+1,j).le.0..AND. lin.ne.2) then
422 
423  ic=i+1
424  jc=j
425 
426  lin=4
427 
428  rxb(ig)=r(i+1)
429  zxb(ig)=xzer(z(j+1),z(j),ut(i+1,j+1),ut(i+1,j))
430 
431  elseif(ut(i+1,j+1)*ut(i,j+1).le.0. .AND. lin.ne.3) then
432 
433  ic=i
434  jc=j+1
435 
436  lin=1
437 
438  rxb(ig)=xzer(r(i+1),r(i),ut(i+1,j+1),ut(i,j+1))
439  zxb(ig)=z(j+1)
440 
441  elseif(ut(i,j+1)*ut(i,j).le.0..AND. lin.ne.4) then
442 
443  ic=i-1
444  jc=j
445 
446  lin=2
447 
448  rxb(ig)=r(i)
449  zxb(ig)=xzer(z(j+1),z(j),ut(i,j+1),ut(i,j))
450 
451  endif
452 
453 c write(6,*) 'loop:ic,jc',ic,jc,ig
454 
455  if(jc.eq.jc1 .AND. ic.eq.ic1) then
456 
457  ig=ig+1
458 
459  rxb(ig)=rxb(2 )
460  zxb(ig)=zxb(2 )
461 
462  go to 1212
463 
464  endif
465 
466  go to 100
467 
468  1212 nxb=ig
469 
470  do 200 ig=1,nxb
471 
472  drx=rxb(ig)-rm
473  dzx=zxb(ig)-zm
474 
475  tetp=datan(dzx/drx)
476 
477  if(ig.ne.1) then
478  if(drx.lt.0.) tetp=tetp+pi
479  deltet=tetp-tetxb(ig-1)
480  if(deltet.lt.0.) tetp=tetp+2.*pi
481  endif
482 
483  tetxb(ig)=tetp
484 
485  roxb(ig)=dsqrt(drx**2+dzx**2)
486 
487 c write(6,*) 'ig,ro,tet',ig,roxb(ig),tetxb(ig)
488 
489  200 continue
490 
491 
492  CALL e01baf(nxb,tetxb,roxb,rrk,cck,nxb+4,wrk,6*nxb+16,ifail)
493 
494 c write(6,*) 'ifail=',ifail
495 
496  DO 210 j=1,ntet
497  tetp=tetpol(j)
498  if(tetp.lt.tetxb(1)) then
499  tetp=tetp+2.*pi
500  elseif(tetp.gt.tetxb(nxb)) then
501  tetp=tetp-2.*pi
502  endif
503  CALL e02bcf(nxb+4,rrk,cck,tetp ,0,cwk,ifail)
504  ro0(j)=cwk(1)
505 c write(6,*) 'ifail=',ifail,j
506  210 CONTINUE
507 
508  return
509  end
510 
511 
512 
513 
514 
515 
516 
517 
e02bcf
subroutine e02bcf(NCAP7, K, C, X, LEFT, S, IFAIL)
Definition: NAG.f:2761
funsq
REAL *8 function funsq(R1, R2, R3, R4, Z1, Z2, Z3, Z4)
Definition: scu.f:803
fun
subroutine fun(X, F)
Definition: Ev2.f:10
gridpl
subroutine gridpl
Definition: B_MESH.f:156
e01baf
subroutine e01baf(M, X, Y, K, C, LCK, WRK, LWRK, IFAIL)
Definition: NAG.f:2511
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
loop95
subroutine loop95(tetpol, ntet, ro0, u0)
Definition: B_MESH.f:328
qst
subroutine qst(qcen, cnor, b0ax, r0ax)
Definition: B_MESH.f:2
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121