ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
_fluxt.f
Go to the documentation of this file.
1  subroutine flux_p(psitok,rk,zk,nk)
2 
3  include 'double.inc'
4  include 'parevo.inc'
5  parameter(nkp=njlim)
6  include 'dim.inc'
7  include 'compol.inc'
8  include 'compol_add.inc'
9  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
10  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
11 
12  real*8 psitok(*),rk(*),zk(*)
13  dimension dg_dn(ntp)
14 
15  sqrt(arg)=dsqrt(arg)
16 
17  i=iplas !!!!!!!!!!!!!!!!!!!!!!
18 
19  do j=2,nt1
20 
21  il=numlin(i,j,nr,nt)
22 
23  a1=a13(i-1,j-1)
24  a2=a34(i-1,j-1)+a12(i-1,j)
25  a3=a24(i-1,j)
26  a4=a23(i-1,j-1)
27  a6=a23(i-1,j)
28 
29  a5=-(a1+a2+a3+a4+a6)
30 
31 
32  g1=psi(i-1,j-1)
33  g2=psi(i-1,j)
34  g3=psi(i-1,j+1)
35 
36  g4=psi(i,j-1)
37  g5=psi(i,j)
38  g6=psi(i,j+1)
39 
40  dgdnl=a1*g1+a2*g2+a3*g3+a4*g4+a5*g5+a6*g6-right(il)
41  dltk=(dlt(i,j-1)+dlt(i,j))*0.5d0
42  dg_dn(j)=dgdnl/dltk
43 
44  enddo
45 
46  dg_dn(1)=dg_dn(nt1)
47  dg_dn(nt)=dg_dn(2)
48 
49 
50 
51  do ik=1,nk
52 
53  psitok(ik)=0.d0
54 
55  rr=rk(ik)
56  zz=zk(ik)
57 
58  do j=2,nt1
59 
60  r0=r(i,j)
61  z0=z(i,j)
62 
63  r1=r(i,j+1)
64  z1=z(i,j+1)
65 
66  call bint(rr,zz,r0,z0,r1,z1,fint,1)
67 
68  psitok(ik)=psitok(ik)-fint*(dg_dn(j)+dg_dn(j+1))*0.5d0
69 
70  enddo
71 
72  enddo
73 
74  return
75  end
76 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
77  subroutine flux_g(psitok,rk,zk,nk)
78 
79  include 'double.inc'
80  include 'parevo.inc'
81  parameter(nkp=njlim)
82  include 'dim.inc'
83  include 'compol.inc'
84  include 'compol_add.inc'
85  common/comaaa/ a12(nrp,ntp),a23(nrp,ntp),a34(nrp,ntp),
86  + a14(nrp,ntp),a13(nrp,ntp),a24(nrp,ntp)
87 
88  real*8 psitok(*),rk(*),zk(*)
89  dimension dg_dn(ntp)
90 
91  sqrt(arg)=dsqrt(arg)
92 
93  i=iplas+2 !!!!!!!!!!!!!!!!!!!!!!
94 
95  do j=2,nt1
96 
97  a1=a13(i-1,j-1)
98  a2=a34(i-1,j-1)+a12(i-1,j)
99  a3=a24(i-1,j)
100  a4=a23(i-1,j-1)
101  a6=a23(i-1,j)
102 
103  a5=-(a1+a2+a3+a4+a6)
104 
105 
106  g1=g(i-1,j-1)
107  g2=g(i-1,j)
108  g3=g(i-1,j+1)
109 
110  g4=g(i,j-1)
111  g5=g(i,j)
112  g6=g(i,j+1)
113 
114  dgdnl=a1*g1+a2*g2+a3*g3+a4*g4+a5*g5+a6*g6
115  dltk=(dlt(i,j-1)+dlt(i,j))*0.5d0
116  dg_dn(j)=dgdnl/dltk
117 
118  enddo
119 
120  dg_dn(1)=dg_dn(nt1)
121  dg_dn(nt)=dg_dn(2)
122 
123 
124 
125  do ik=1,nk
126 
127  psitok(ik)=0.d0
128 
129  rr=rk(ik)
130  zz=zk(ik)
131 
132  do j=2,nt1
133 
134  r0=r(i,j)
135  z0=z(i,j)
136  g0=g(i,j)
137 
138  r1=r(i,j+1)
139  z1=z(i,j+1)
140  g1=g(i,j+1)
141 
142  r_05=0.5d0*(r0+r1)
143  z_05=0.5d0*(z0+z1)
144 
145  d_r=r1-r0
146  d_z=z1-z0
147 
148  call greng(rr,zz, r_05,z_05, fgreen,dgdr,dgdz)
149 
150  dgr_dn=(dgdr*d_z-dgdz*d_r)/r_05
151 
152  fint=
153  * 0.5d0*( (g0+g1)*dgr_dn - (dg_dn(j)+dg_dn(j+1))*fgreen*dlt(i,j) )
154  * /pi
155 
156 
157  psitok(ik)=psitok(ik)+fint
158 
159  enddo
160 
161  enddo
162 
163  return
164  end
165 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
166 
167 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
168 
169  subroutine flux_r(psitok,ncequi)
170 
171  include 'double.inc'
172 
173  include 'prm.inc'
174  parameter(nekp=npfc0+nplim)
175  parameter(nkp=njlim)
176  include 'dimpl1.inc'
177  include 'dimpl2.inc'
178  include 'parrc1.inc'
179  include 'comrec.inc'
180  include 'compol.inc'
181  include 'compol_add.inc'
182 
183  common /comext/ zaindk(nip,njp,nekp)
184  real *4 zaindk
185  real*8 psitok(*)
186 
187 
188  ddx=x(2)-x(1)
189  ddy=y(2)-y(1)
190 
191  sqcen=0.d0
192  do j=2,nt1
193  sqcen=sqcen+sq1(1,j)+sq4(1,j)
194  enddo
195 
196  r0=rm
197  z0=zm
198 
199  ic=(r0-x(1))/ddx+1
200  jc=(z0-y(1))/ddy+1
201 
202  r1=x(ic)
203  r2=x(ic+1)
204 
205  z1=y(jc)
206  z2=y(jc+1)
207 
208  do k=1,ncequi
209 
210  u1=zaindk(ic,jc,k)
211  u2=zaindk(ic+1,jc,k)
212  u3=zaindk(ic+1,jc+1,k)
213  u4=zaindk(ic,jc+1,k)
214 
215  psitok(k)=blin_(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4)*cur(1,2)*sqcen
216 
217  enddo
218 
219  do i=2,iplas
220  do j=2,nt1
221 
222  if(i.ne.iplas) then
223  sqk=sq1(i,j)+sq2(i-1,j)+sq3(i-1,j-1)+sq4(i,j-1)
224  else
225  sqk=sq2(i-1,j)+sq3(i-1,j-1)
226  endif
227 
228  r0=r(i,j)
229  z0=z(i,j)
230 
231  ic=(r0-x(1))/ddx+1
232  jc=(z0-y(1))/ddy+1
233 
234  r1=x(ic)
235  r2=x(ic+1)
236 
237  z1=y(jc)
238  z2=y(jc+1)
239 
240 
241  do k=1,ncequi
242 
243  u1=zaindk(ic,jc,k)
244  u2=zaindk(ic+1,jc,k)
245  u3=zaindk(ic+1,jc+1,k)
246  u4=zaindk(ic,jc+1,k)
247 
248  psitok(k)=psitok(k)+blin_(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4 )
249  * *cur(i,j)*sqk
250 
251  enddo
252 
253  enddo
254  enddo
255 
256  return
257  end
258 
259  subroutine numcel(rrk,zzk,icell,jcell)
260 
261  include 'double.inc'
262  include 'parevo.inc'
263  parameter(nkp=njlim)
264  include 'dim.inc'
265  include 'compol.inc'
266  include 'compol_add.inc'
267 
268  sqrt(arg)=dsqrt(arg)
269 
270  dr0=rrk-rm
271  dz0=zzk-zm
272 
273  ro0=sqrt(dr0**2+dz0**2)
274 
275  tetp=dacos(dr0/ro0)
276  if(dz0.lt.0.d0) then
277  tet0=2.d0*pi-tetp
278  else
279  tet0=tetp
280  endif
281 
282  if(tet0.lt.teta(1)) tet0=tet0+2.d0*pi
283  if(tet0.gt.teta(nt)) tet0=tet0-2.d0*pi
284 
285  if(tet0.lt.teta(1)) tet0=tet0+2.d0*pi
286  if(tet0.gt.teta(nt)) tet0=tet0-2.d0*pi
287 
288  if(tet0.lt.teta(1)) pause 'numcel:tet0<teta(1)'
289  if(tet0.gt.teta(nt)) pause 'numcel:tet0>teta(nt)'
290 
291  do j=1,nt1
292 
293  if(tet0.ge.teta(j) .AND. tet0.lt.teta(j+1)) jc=j
294 
295  enddo
296 
297  jcell=jc
298 
299  drvb=rrk-r(nr,jc)
300  dzvb=zzk-z(nr,jc)
301 
302  drcb=r(nr,jc+1)-r(nr,jc)
303  dzcb=z(nr,jc+1)-z(nr,jc)
304 
305  vecpro=drcb*dzvb-drvb*dzcb
306 
307 
308  if(vecpro.le.0.d0) then
309 
310  icell=nr
311 
312  else
313 
314  do i=nr,2,-1
315 
316  drv=rrk-r(i,jc)
317  dzv=zzk-z(i,jc)
318 
319  drc=r(i,jc+1)-r(i,jc)
320  dzc=z(i,jc+1)-z(i,jc)
321 
322  vecpro=drc*dzv-drv*dzc
323 
324  if(vecpro.gt.0.d0) then
325  ic=i-1
326  else
327  go to 100
328  endif
329 
330  enddo
331 
332  100 icell=ic
333 
334  endif
335 
336  return
337  end
338 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
339 
340  subroutine f_flux(psitok,rk,zk,nk)
341 
342  include 'double.inc'
343  include 'parevo.inc'
344  parameter(nkp=njlim)
345  include 'dim.inc'
346  include 'compol.inc'
347  include 'compol_add.inc'
348 
349  real*8 psitok(*),rk(*),zk(*)
350 
351  sqrt(arg)=dsqrt(arg)
352 
353  ik_out=0
354  do 100 ik=1,nk
355  ! write(6,*) 'ik:',ik
356  psitok(ik)=0.d0
357 
358  rr=rk(ik)
359  zz=zk(ik)
360 
361  do 110 j=2,nt1
362 
363  psitok(ik)=psitok(ik)
364  * -pinadg(ik,j)*(dgdn(j)+dgdn(j+1))*0.5d0
365 
366  110 continue
367 
368  if(iprcon(ik).eq.0) then !!!
369 
370  call numcel(rr,zz,ic,jc)
371 
372  if(ic.ne.nr) then
373  write(*,*) 'flux:conductor is in the box',ic,nr
374  write(*,*) 'program is terminated'
375  write(*,*) 'you need to consult program provider'
376  stop
377  endif
378 
379  ik_out=ik_out+1
380 
381 ! do 110 j=2,nt1
382 !
383 ! psitok(ik)=psitok(ik)
384 ! * -pinadg(ik_out,j)*(dgdn(j)+dgdn(j+1))*0.5d0
385 !
386 ! 110 continue
387 
388 
389  else !!!!!!!!!!!!!!!!!!!!!
390 
391  call numcel(rr,zz,ic,jc)
392  write(*,*) 'ik ic jc',ik,ic,jc
393 
394  if(ic.eq.nr) then
395  write(*,*) 'flux:conductor is out the box',ic,nr
396  write(*,*) 'program is terminated'
397  write(*,*) 'you need to consult program provider'
398  stop
399  endif
400 
401  r1=r(ic,jc)
402  r2=r(ic+1,jc)
403  r3=r(ic+1,jc+1)
404  r4=r(ic,jc+1)
405 
406  z1=z(ic,jc)
407  z2=z(ic+1,jc)
408  z3=z(ic+1,jc+1)
409  z4=z(ic,jc+1)
410 
411  !u1=psii(ic,jc)
412  !u2=psii(ic+1,jc)
413  !u3=psii(ic+1,jc+1)
414  !u4=psii(ic,jc+1)
415 
416  u1=g(ic,jc)
417  u2=g(ic+1,jc)
418  u3=g(ic+1,jc+1)
419  u4=g(ic,jc+1)
420 
421  call blic_d(rr,zz,r1,r2,r3,r4,z1,z2,z3,z4,
422  * u1,u2,u3,u4,u0,dudr,dudz)
423 
424  !psitok(ik)=u0
425  psitok(ik)=psitok(ik)+u0
426 
427  endif !!!!!!!!!!!!!!!!!!!!
428 
429  100 continue
430 
431  if(ik_out.ne.nk_out) then
432  write(*,*) .ne.'ik_outnk_out',ik_out,nk_out
433  write(*,*) 'program is terminated'
434  write(*,*) 'you need to consult program provider'
435  stop
436  endif
437 
438  return
439  end
440 
441 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
442 
443  real*8 function blin_tr(tet0,ro0,
444  * tet1,tet2,ro1,ro2,ro3,ro4,u1,u2,u3,u4)
445 
446  include 'double.inc'
447 
448  ro14=(ro1*(tet2-tet0)+ro4*(tet0-tet1))/(tet2-tet1)
449  ro23=(ro2*(tet2-tet0)+ro3*(tet0-tet1))/(tet2-tet1)
450 
451  u14=(u1*(tet2-tet0)+u4*(tet0-tet1))/(tet2-tet1)
452  u23=(u2*(tet2-tet0)+u3*(tet0-tet1))/(tet2-tet1)
453 
454  blin_tr=(u14*(ro23-ro0)+u23*(ro0-ro14))/(ro23-ro14)
455 
456  return
457  end
458 
459 c====================================================================
460 
f_flux
subroutine f_flux(psitok, rk, zk, nk)
Definition: _fluxt.f:340
blin_tr
real *8 function blin_tr(tet0, ro0,
Definition: _fluxt.f:443
drv
subroutine drv(ZIX1, ZIY1, ZIX2, ZIY2)
Definition: ppplib.f:5671
flux_p
subroutine flux_p(psitok, rk, zk, nk)
Definition: _fluxt.f:1
flux_r
subroutine flux_r(psitok, ncequi)
Definition: _fluxt.f:169
numcel
subroutine numcel(rrk, zzk, icell, jcell)
Definition: _fluxt.f:259
blic_d
subroutine blic_d(r0, z0, r1, r2, r3, r4, z1, z2, z3, z4,
Definition: _extp.f:144
numlin
function numlin(i, j, nr, nt)
Definition: com_sub.f:1040
bint
subroutine bint(X, Y, R0, Z0, r1, z1, F, I)
Definition: scu.f:557
flux_g
subroutine flux_g(psitok, rk, zk, nk)
Definition: _fluxt.f:77
blin_
real *8 function blin_(r0, z0, r1, r2, z1, z2, u1, u2, u3, u4)
Definition: scu.f:857
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
greng
subroutine greng(R0, Z0, R, Z, fgreen, dGdr, dGdz)
Definition: scu.f:658