ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
impurity_one.F90
Go to the documentation of this file.
1 ! +++ IMPURITY TRANSPORT EQUATIONS FOR ONE IMPURITY +++
2 
3 
4 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
5 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
6 #ifdef GOT_AMNSPROTO
7 SUBROUTINE impurity_one(TE,NE,NZ1,NZM1,VPR,VPRM,R0,BT,BTPRIME,DIFF,FLUX,FLUX_INTER,rho,&
8  vconv,nrho,simp2,nsource,nz_bnd,nz_bnd_type,control_double,control_integer,g1,imp_radiation,se_exp,&
9  max_nzimp,amns_ei,amns_rc,amns_lr,amns_br,amns_eip,lin_rad1,brem_rad1,jon_en1,rec_los1)
10 #else
11 SUBROUTINE impurity_one(TE,NE,NZ1,NZM1,VPR,VPRM,R0,BT,BTPRIME,DIFF,FLUX,FLUX_INTER,rho,&
12  vconv,nrho,simp2,nsource,nz_bnd,nz_bnd_type,control_double,control_integer,g1,imp_radiation,se_exp,&
13  max_nzimp)
14 #endif
15 
16 !--------------------------------------------------------------!
17 ! This subroutine solves impuriy particle transport !
18 ! equations for impurity components from 1 to NSTATE, !
19 ! and provides: density and flux of impurity components !
20 ! from 1 to NSTATE !
21 !--------------------------------------------------------------!
22 ! Source: --- !
23 ! Developers: R.Stankiewicz,I.M.Ivanova-Stanik !
24 ! Kontacts: irena@ifpilm.waw.pl !
25 ! !
26 ! Comments: might change after the ITM !
27 ! data stucture is finalized !
28 ! !
29 !--------------------------------------------------------------!
30 
31  USE itm_types
32  USE type_solver
33  USE itm_constants
34 #ifdef GOT_AMNSPROTO
35  use amns_types
36  use amns_module
37 #endif
38 
39 
40  IMPLICIT NONE
41 
42 
43 ! +++ Input/Output to numerical solver:
44  TYPE (numerics) :: solver !contains all I/O quantities to numerics part
45 
46 
47 ! +++ Internal parameters:
48  INTEGER :: irho, nrho !radius index, number of radial points
49  INTEGER :: iimp, simp2,simp,isimp !index of impurity component, number of considered impurity components (max ionization state)
50  INTEGER :: izimp, nzimp
51  INTEGER :: nz_bnd_type(simp2) !boundary condition, type
52  INTEGER :: index_t !index for linear interpolation for atomic data
53  INTEGER :: max_nzimp
54 
55  REAL (R8) :: bt, btm, btprime !magnetic field from current time step, [T], previous time steps, [T], time derivative, [T/s]
56  REAL (R8) :: r0
57  REAL (R8) :: rho(nrho) !normalised minor radius, [m]
58  real(r8) :: vpr(nrho) !V', [m^2]
59  real(r8) :: vprm(nrho) !V' (at previous time step), [m^2]
60  REAL (R8) :: g1(nrho)
61  REAL (R8) :: nz1(nrho,simp2) !density from current ans previous time step, [m^-3]
62  REAL (R8) :: nzm1(nrho,simp2) !density from previous time step,
63  real(r8) :: flux(nrho,simp2) !ion flux, [1/s]
64  real(r8) :: flux_inter(nrho,simp2) !ion flux, contributing to convective heat transport [1/s]
65  REAL (R8) :: ne(nrho),te(nrho) !electron density [m^-3]
66  REAL (R8) :: diff(nrho,simp2) !diffusion coefficient [m^2/s] and [m/s]
67  REAL (R8) :: vconv(nrho,simp2) !pinch velocity, [m^2/s] and [m/s]
68  REAL (R8) :: nsource(nrho,simp2)
69  REAL (R8) :: nz_bnd(3,simp2) !boundary condition, value, [depends on NZ_BND_TYPE]
70 
71  REAL (R8) :: real_index_t !index for linear interpolation for atomic data
72  REAL (R8) :: alfa(nrho,simp2) !atom data: ionization coefficient(after interpolation)[m^3/s]
73  REAL (R8) :: beta(nrho,simp2) !atom data: recombination coefficient (after interpolation) [m^3/s]
74  REAL (R8) :: gama(nrho,simp2) !atom data: cooling coefficient (after interpolation)
75  REAL (R8) :: slin(nrho,simp2) !atom data: linear radiation (after interpolation)
76  REAL (R8) :: imp_radiation(nrho,simp2) !impurity radiation
77  REAL (R8) :: alfa_nz(max_nzimp+2,500) !atom data: ionization coefficient [m^3/s]
78  REAL (R8) :: beta_nz(max_nzimp+2,500) !atom data: recombination coefficient [m^3/s]
79  REAL (R8) :: slin_nz(max_nzimp+2,500) !atom data: linear radiation
80  REAL (R8) :: gama_nz(max_nzimp+2,500) !atom data: cooling coefficient
81  REAL (R8) :: potential(nrho,1:max_nzimp+2)
82  REAL (R8) :: amix, tau !mixing factor, time step, [s]
83 
84 ! radiation
85  REAL (R8) :: lin_rad1(nrho,simp2) !profile of lineradiation for one impurity
86  REAL (R8) :: brem_rad1(nrho,simp2) !profile of bremst. for one impurity
87  REAL (R8) :: jon_en1(nrho,simp2) !profile of jonisation energy for one impurity
88  REAL (R8) :: rec_los1(nrho,simp2) !profile of bremst. for one impurity
89 
90  REAL (R8) :: se_exp(nrho)
91 
92  INTEGER :: solut_method
93  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
94  INTEGER :: ndim !number of equations to be solved
95  INTEGER :: solver_type !specifies the option for numerical solution
96  REAL (R8) :: y(nrho) !function at the current amd previous time steps
97  REAL (R8) :: ym(nrho) !function at the current amd previous time steps
98  REAL (R8) :: dy(nrho) !derivative of function
99  REAL (R8) :: a(nrho) !coefficients for numerical solver
100  REAL (R8) :: b(nrho) !coefficients for numerical solver
101  REAL (R8) :: c(nrho) !coefficients for numerical solver
102  REAL (R8) :: d(nrho) !coefficients for numerical solver
103  REAL (R8) :: e(nrho) !coefficients for numerical solver
104  REAL (R8) :: f(nrho) !coefficients for numerical solver
105  REAL (R8) :: g(nrho) !coefficients for numerical solver
106  REAL (R8) :: h !coefficients for numerical solver
107  REAL (R8) :: v(2), u(2), w(2) !boundary conditions for numerical solver
108 
109  REAL (R8) :: fun1(nrho), intfun1(nrho)
110 
111 
112  INTEGER, INTENT(IN) :: control_integer(2) !integer control parameters
113  REAL (R8), INTENT(IN) :: control_double(5) !real control parameters
114 
115  REAL (R8) :: dnz(nrho) !density gradient, [m^-4]
116  REAL (R8) :: int_source(nrho) !integral of source [1/s]
117  REAL (R8) :: nzp(nrho,simp2) !impurity density (Alternated Direction method)
118  REAL (R8) :: a1p(nrho), b1p(nrho) !coefficients for numerical solver (AD method)
119  REAL (R8) :: c1p(nrho), d1pp(nrho) !coefficients for numerical solver (AD method)
120  REAL (R8) :: vp(2), up(2), wp(2) !coefficients for numerical solver (AD method)
121  REAL (R8) :: drho, drhop, drhom
122  REAL (R8) :: a2p(simp2), b2p(simp2) !coefficients for numerical solver (AD method)
123  REAL (R8) :: c2p(simp2), d2pp(simp2) !coefficients for numerical solver (AD method)
124 
125  REAL (R8) :: time
126 
127  REAL :: densityfluxin(nrho)
128  REAL :: fluxin, fluxout
129  REAL :: calkatrapez
130 ! declaration for cubint
131  INTEGER ::ntab
132  REAL ::xtab(nrho)
133  REAL ::ftab(nrho)
134  INTEGER ::ia_in
135  INTEGER ::ib_in
136  REAL ::result
137  REAL ::error
138 ! ***************
139  INTEGER :: ifail
140 
141 #ifdef GOT_AMNSPROTO
142  type (amns_handle_rx_type) :: amns_ei(1:max_nzimp+1), amns_rc(1:max_nzimp+1), &
143  amns_lr(1:max_nzimp+1), amns_br(1:max_nzimp+1),amns_eip(1:max_nzimp+1)
144 #endif
145 
146 ! +++ Set up dimensions:
147  ndim = 1 !no coupling between density equations
148  simp = simp2-2
149  nzimp = simp
150 
151 ! +++ Allocate types for interface with numerical solver:
152  CALL allocate_numerics(ndim, nrho, solver, ifail)
153 
154 
155 ! +++ Set up local variables:
156 
157 !irena
158  solut_method = 1 !temporary fix
159 
160  amix = control_double(2)
161  tau = control_double(1) ! 0.5 - Lackner's method
162  solver_type = control_integer(1)
163 
164  btm =bt-btprime*tau
165 
166 ! +++ AtomData Coefficients: Linear Interpolation
167 
168  alfa = 0.0d0
169  beta = 0.0d0
170  gama = 0.0d0
171  slin = 0.0d0
172  potential = 0.0d0
173 !++++ radiation+energy losses
174  lin_rad1 = 0.0d0
175  brem_rad1 = 0.0d0
176  jon_en1 = 0.0d0
177  rec_los1 = 0.0d0
178 
179 
180 #ifdef GOT_AMNSPROTO
181  do izimp=1, nzimp+1
182 
183 
184 
185 
186  if(izimp .ne. nzimp+1) then
187  call itm_amns_rx(amns_ei(izimp),alfa(:,izimp),te,ne)
188 !DPC write(*,*) 'alfa(nrho,izimp)=',alfa(nrho,izimp)
189  call itm_amns_rx(amns_lr(izimp),slin(:,izimp),te,ne) ! guess: slin == line radiation loss
190  call itm_amns_rx(amns_eip(izimp+1),potential(:,izimp),te,ne) ! guess: potential= potential of ionization
191 
192 
193  endif
194 
195  if(izimp .ne. 1) then
196  call itm_amns_rx(amns_rc(izimp),beta(:,izimp),te,ne)
197  call itm_amns_rx(amns_br(izimp),gama(:,izimp),te,ne) ! guess: gama == bremsstrahlung + recombination energy loss
198 
199  endif
200 
201 
202  enddo
203 
204 
205 
206 !!! stop
207 #else
208 
209  call read_data(alfa_nz,beta_nz,gama_nz,slin_nz,max_nzimp,simp)
210  slin_nz = slin_nz * 1.e-30 ! DPC addition
211  gama_nz = gama_nz * 1.e-20 ! DPC addition --- this one is just a guess!
212 
213 
214 
215  rho_loopa: DO irho=1, nrho
216 
217  index_t = int( 499.0/5.0*log10(te(irho)) ) + 1
218  real_index_t = 499.0/5.0*log10(te(irho)) + 1.0
219 
220  imp_loopa: DO isimp=1,simp+1
221 
222  alfa(irho,isimp) = ( alfa_nz(isimp,index_t+1) - alfa_nz(isimp,index_t) )*(real_index_t - index_t) + alfa_nz(isimp,index_t)
223  beta(irho,isimp) = ( beta_nz(isimp,index_t+1) - beta_nz(isimp,index_t) )*(real_index_t - index_t) + beta_nz(isimp,index_t)
224  slin(irho,isimp) = ( slin_nz(isimp,index_t+1) - slin_nz(isimp,index_t) )*(real_index_t - index_t) + slin_nz(isimp,index_t)
225  gama(irho,isimp) = ( gama_nz(isimp,index_t+1) - gama_nz(isimp,index_t) )*(real_index_t - index_t) + gama_nz(isimp,index_t)
226 
227  END DO imp_loopa
228 
229  END DO rho_loopa
230 #endif
231 
232 
233 !DPC write(*,'(a,1p,100(e15.6))') 'alfa: ', (alfa(nrho/2,izimp),izimp=1,nzimp+1)
234 !DPC write(*,'(a,1p,100(e15.6))') 'beta: ', (beta(nrho/2,izimp),izimp=1,nzimp+1)
235 !DPC write(*,'(a,1p,100(e15.6))') 'gama: ', (gama(nrho/2,izimp),izimp=1,nzimp+1)
236 !DPC write(*,'(a,1p,100(e15.6))') 'slin: ', (slin(nrho/2,izimp),izimp=1,nzimp+1)
237 !DPC write(*,'(a,1p,100(e15.6))') 'potential: ', (potential(nrho/2,izimp),izimp=1,nzimp+1)
238 
239 !DPC write(*,'(a,1p,100(e25.16))') 'BOUNDARY atomic physics', &
240 !DPC te(nrho), ne(nrho), alfa(nrho,1:nzimp+1), slin(nrho,1:nzimp+1), beta(nrho,1:nzimp+1), gama(nrho,1:nzimp+1)
241 
242  SELECT CASE (solut_method)
243 
244 
245  CASE(1)
246 
247 
248 
249 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
250 ! solution of particle transport equation for
251 ! individual state of impurity
252 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
253 
254  imp_loop1: DO isimp=2,simp+1 ! "1" - Neutral impurity, "2" - "+1", ...
255 
256 ! +++ Set equation to 'predictive' and all coefficients to zero:
257  flag = 1
258  y(:) = 0.0d0
259  dy(:) = 0.0d0
260  ym(:) = 0.0d0
261  a(:) = 0.0d0
262  b(:) = 0.0d0
263  c(:) = 0.0d0
264  d(:) = 0.0d0
265  e(:) = 0.0d0
266  f(:) = 0.0d0
267  g(:) = 0.0d0
268  h = 0.0d0
269  v(:) = 0.0d0
270  u(:) = 0.0d0
271  w(:) = 0.0d0
272 
273 
274 ! +++ Coefficients for ion diffusion equation in form:
275 !
276 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
277 
278  rho_loop3: DO irho=1,nrho
279 
280  y(irho) = nz1(irho,isimp)
281  ym(irho) = nzm1(irho,isimp)
282  a(irho) = vpr(irho)
283  b(irho) = vprm(irho)
284  c(irho) = 1.d0
285  d(irho) = vpr(irho)*g1(irho)*diff(irho,isimp)
286  e(irho) = vpr(irho)*g1(irho)*vconv(irho,isimp) &
287  - btprime/2.d0/bt*rho(irho)*vpr(irho)
288  f(irho) = vpr(irho)*nsource(irho,isimp)+vpr(irho)*( nz1(irho,isimp-1)*ne(irho)*alfa(irho,isimp-1) &
289  + nz1(irho,isimp+1)*ne(irho)*beta(irho,isimp+1) )
290  g(irho) = vpr(irho)*(ne(irho)*alfa(irho,isimp) + ne(irho)*beta(irho,isimp) )
291  IF (irho.eq.nrho) then
292 
293 
294 
295  ENDIF
296 
297  END DO rho_loop3
298 
299 
300  h = 0.5*tau
301 ! +++ Boundary conditions for ion diffusion equation in form:
302 !
303 ! V*Y' + U*Y =W
304 !
305 ! +++ On axis:
306 ! dNZ/drho(rho=0)=0:
307  IF(diff(1,isimp).GT.0.d0) THEN
308  v(1) = -diff(1,isimp)
309  u(1) = vconv(1,isimp)
310  w(1) = 0.d0
311  ELSE
312  v(1) = 1.d0
313  u(1) = 0.d0
314  w(1) = 0.d0
315  END IF
316 
317 
318 ! +++ At the edge:
319 ! FIXED NZ
320 
321  IF(nz_bnd_type(isimp).EQ.1) THEN
322  v(2) = 0.d0
323  u(2) = 1.d0
324  w(2) = nz_bnd(1,isimp)
325  ENDIF
326 
327 ! FIXED grad_NZ
328  IF(nz_bnd_type(isimp).EQ.2) THEN
329  v(2) = 1.d0
330  u(2) = 0.d0
331  w(2) = nz_bnd(1,isimp)
332  ENDIF
333 
334 ! FIXED L_NZ
335  IF(nz_bnd_type(isimp).EQ.3) THEN
336  v(2) = nz_bnd(1,isimp)
337  u(2) = 1.d0
338  w(2) = 0.d0
339  ENDIF
340 
341 ! FIXED Flux_NZ
342  IF(nz_bnd_type(isimp).EQ.4) THEN
343  v(2) = -g1(nrho)*diff(nrho,isimp)*vpr(nrho)
344  u(2) = g1(nrho)*vconv(nrho,isimp)*vpr(nrho)
345  w(2) = nz_bnd(1,isimp)
346  ENDIF
347 
348 ! Generic boundary condition
349  IF(nz_bnd_type(isimp).EQ.5) THEN
350  v(2) = nz_bnd(1,isimp)
351  u(2) = nz_bnd(2,isimp)
352  w(2) = nz_bnd(3,isimp)
353  ENDIF
354 
355 ! +++ Density equation is not solved:
356  IF(nz_bnd_type(isimp).EQ.0) THEN
357 
358  CALL derivn(nrho,rho,y,dy) !temperature gradient
359 
360  flag = 0
361 
362  rho_loop4: DO irho=1,nrho
363 
364  a(irho) = 1.0d0
365  b(irho) = 1.0d0
366  c(irho) = 1.0d0
367  d(irho) = 0.0d0
368  e(irho) = 0.0d0
369  f(irho) = 0.0d0
370  g(irho) = 0.0d0
371 
372  END DO rho_loop4
373 
374  v(2) = 0.0d0
375  u(2) = 1.0d0
376  w(2) = y(nrho)
377  END IF
378 
379 
380 
381 ! +++ Defining coefficients for numerical solver:
382  solver%TYPE = solver_type
383  solver%EQ_FLAG(ndim) = flag
384  solver%NDIM = ndim
385  solver%NRHO = nrho
386  solver%AMIX = amix
387 
388 
389  rho_loop5: DO irho=1,nrho
390 
391  solver%RHO(irho) = rho(irho)
392  solver%Y(ndim,irho) = y(irho)
393  solver%DY(ndim,irho) = dy(irho)
394  solver%YM(ndim,irho) = ym(irho)
395  solver%A(ndim,irho) = a(irho)
396  solver%B(ndim,irho) = b(irho)
397  solver%C(ndim,irho) = c(irho)
398  solver%D(ndim,irho) = d(irho)
399  solver%E(ndim,irho) = e(irho)
400  solver%F(ndim,irho) = f(irho)
401  solver%G(ndim,irho) = g(irho)
402 
403  END DO rho_loop5
404 
405  solver%H = h
406 
407  solver%V(ndim,1) = v(1)
408  solver%U(ndim,1) = u(1)
409  solver%W(ndim,1) = w(1)
410  solver%V(ndim,2) = v(2)
411  solver%U(ndim,2) = u(2)
412  solver%W(ndim,2) = w(2)
413 
414 
415 ! +++ Solution of density diffusion equation:
416  CALL solution_interface(solver, ifail)
417 
418 
419 ! +++ New impurity density:
420  rho_loop6: DO irho=1,nrho
421 
422  y(irho) = solver%Y(ndim,irho)
423  dy(irho) = solver%DY(ndim,irho)
424 
425  END DO rho_loop6
426 
427 
428 
429 ! +++ New profiles of impurity density flux and integral source:
430  rho_loop7: DO irho=1,nrho
431 
432  nzm1(irho,isimp) = nz1(irho,isimp)
433  nz1(irho,isimp) = y(irho)
434  dnz(irho) = dy(irho)
435  IF (rho(irho).NE.0.0d0) THEN
436  fun1(irho) = 1.d0/rho(irho)*(vpr(irho)*nsource(irho,isimp) &
437  + vprm(irho)*nzm1(irho,isimp)/tau &
438  - nz1(irho,isimp)*vpr(irho)*(1.d0/tau))
439  ELSE IF (rho(irho).EQ.0.0d0) THEN
440  fun1(irho) = 4.d0*itm_pi**2*r0* (nsource(irho,isimp) &
441  + nzm1(irho,isimp)/tau - nz1(irho,isimp)*(1.d0/tau))
442  END IF
443 
444 
445 
446  END DO rho_loop7
447 
448 
449 
450  CALL integr(nrho,rho,fun1,intfun1) !Integral source
451 
452  rho_loop8: DO irho=1,nrho
453  flux_inter(irho,isimp) = intfun1(irho) &
454  + btprime/2.d0/bt*rho(irho)*vpr(irho)*y(irho)
455 
456 
457  flux(irho,isimp) = vpr(irho)*g1(irho)* &
458  ( y(irho)*vconv(irho,isimp) - dy(irho)*diff(irho,isimp) )
459 
460  END DO rho_loop8
461 
462 
463  END DO imp_loop1
464 
465 
466  loop_change: DO isimp=2,simp+1
467  DO irho=1,nrho
468 
469  nzm1(irho,isimp)=nz1(irho,isimp)
470 
471  end do
472  end do loop_change
473 
474 
475 
476 
477 
478 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
479 ! solution of particle transport equation for
480 ! individual state of impurity
481 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
482 
483  imp_loop21: DO isimp=simp+1, 2, -1 ! "1" - Neutral impurity, "2" - "+1", ...
484 
485 
486 
487 
488 ! +++ Set equation to 'predictive' and all coefficients to zero:
489  flag = 1
490  y(:) = 0.0d0
491  dy(:) = 0.0d0
492  ym(:) = 0.0d0
493  a(:) = 0.0d0
494  b(:) = 0.0d0
495  c(:) = 0.0d0
496  d(:) = 0.0d0
497  e(:) = 0.0d0
498  f(:) = 0.0d0
499  g(:) = 0.0d0
500  h = 0.0d0
501  v(:) = 0.0d0
502  u(:) = 0.0d0
503  w(:) = 0.0d0
504 
505 ! +++ Coefficients for ion diffusion equation in form:
506 !
507 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
508 
509  rho_loop23: DO irho=1,nrho
510 
511  y(irho) = nz1(irho,isimp)
512  ym(irho) = nzm1(irho,isimp)
513  a(irho) = vpr(irho)
514  b(irho) = vprm(irho)
515  c(irho) = 1.d0
516  d(irho) = vpr(irho)*g1(irho)*diff(irho,isimp)
517  e(irho) = vpr(irho)*g1(irho)*vconv(irho,isimp) &
518  - btprime/2.d0/bt*rho(irho)*vpr(irho)
519  f(irho) = vpr(irho)*nsource(irho,isimp)+vpr(irho)*(nz1(irho,isimp-1)*ne(irho)*alfa(irho,isimp-1) &
520  + nz1(irho,isimp+1)*ne(irho)*beta(irho,isimp+1) )
521  g(irho) = vpr(irho)*(ne(irho)*alfa(irho,isimp) + ne(irho)*beta(irho,isimp) )
522 
523 
524 
525 
526  END DO rho_loop23
527 
528  h = 0.5*tau
529 
530 
531 
532 
533 ! +++ Boundary conditions for ion diffusion equation in form:
534 !
535 ! V*Y' + U*Y =W
536 !
537 ! +++ On axis:
538  IF(diff(1,isimp).GT.0.d0) THEN
539  v(1) = -diff(1,isimp)
540  u(1) = vconv(1,isimp)
541  w(1) = 0.d0
542  ELSE
543  v(1) = 1.d0
544  u(1) = 0.d0
545  w(1) = 0.d0
546  END IF
547 
548 ! +++ At the edge:
549 ! FIXED NZ
550 !!write(*,*)'NZ_BND_TYPE(ISIMP)=',NZ_BND_TYPE(ISIMP),'odwrotnej'
551 !!pause
552  IF(nz_bnd_type(isimp).EQ.1) THEN
553  v(2) = 0.d0
554  u(2) = 1.d0
555  w(2) = nz_bnd(1,isimp)
556  ENDIF
557 
558 ! FIXED grad_NZ
559  IF(nz_bnd_type(isimp).EQ.2) THEN
560  v(2) = 1.d0
561  u(2) = 0.d0
562  w(2) = nz_bnd(1,isimp)
563  ENDIF
564 
565 ! FIXED L_NZ
566  IF(nz_bnd_type(isimp).EQ.3) THEN
567  v(2) = nz_bnd(1,isimp)
568  u(2) = 1.d0
569  w(2) = 0.d0
570  ENDIF
571 
572 ! FIXED Flux_NZ
573  IF(nz_bnd_type(isimp).EQ.4) THEN
574  v(2) = -g1(nrho)*diff(nrho,isimp)*vpr(nrho)
575  u(2) = g1(nrho)*vconv(nrho,isimp)*vpr(nrho)
576  w(2) = nz_bnd(1,isimp)
577  ENDIF
578 
579 ! Generic boundary condition
580  IF(nz_bnd_type(isimp).EQ.5) THEN
581  v(2) = nz_bnd(1,isimp)
582  u(2) = nz_bnd(2,isimp)
583  w(2) = nz_bnd(3,isimp)
584  ENDIF
585 
586 
587 
588 ! +++ Density equation is not solved:
589  IF(nz_bnd_type(isimp).EQ.0) THEN
590 
591  CALL derivn(nrho,rho,y,dy) !temperature gradient
592 
593  flag = 0
594 
595  rho_loop24: DO irho=1,nrho
596 
597  a(irho) = 1.0d0
598  b(irho) = 1.0d0
599  c(irho) = 1.0d0
600  d(irho) = 0.0d0
601  e(irho) = 0.0d0
602  f(irho) = 0.0d0
603  g(irho) = 0.0d0
604 
605  END DO rho_loop24
606 
607  v(2) = 0.0d0
608  u(2) = 1.0d0
609  w(2) = y(nrho)
610  END IF
611 
612 
613 ! +++ Defining coefficients for numerical solver:
614  solver%TYPE = solver_type
615  solver%EQ_FLAG(ndim) = flag
616  solver%NDIM = ndim
617  solver%NRHO = nrho
618  solver%AMIX = amix
619 
620  rho_loop25: DO irho=1,nrho
621 
622  solver%RHO(irho) = rho(irho)
623  solver%Y(ndim,irho) = y(irho)
624  solver%DY(ndim,irho) = dy(irho)
625  solver%YM(ndim,irho) = ym(irho)
626  solver%A(ndim,irho) = a(irho)
627  solver%B(ndim,irho) = b(irho)
628  solver%C(ndim,irho) = c(irho)
629  solver%D(ndim,irho) = d(irho)
630  solver%E(ndim,irho) = e(irho)
631  solver%F(ndim,irho) = f(irho)
632  solver%G(ndim,irho) = g(irho)
633 
634  END DO rho_loop25
635 
636  solver%H = h
637  solver%V(ndim,1) = v(1)
638  solver%U(ndim,1) = u(1)
639  solver%W(ndim,1) = w(1)
640  solver%V(ndim,2) = v(2)
641  solver%U(ndim,2) = u(2)
642  solver%W(ndim,2) = w(2)
643 
644 
645 
646 ! +++ Solution of density diffusion equation:
647  CALL solution_interface(solver, ifail)
648 
649 
650 ! +++ New impurity density:
651  rho_loop26: DO irho=1,nrho
652 
653  y(irho) = solver%Y(ndim,irho)
654  dy(irho) = solver%DY(ndim,irho)
655 
656  END DO rho_loop26
657 
658 ! +++ New profiles of impurity density flux and integral source:
659  rho_loop27: DO irho=1,nrho
660 
661  nzm1(irho,isimp) = nz1(irho,isimp)
662  nz1(irho,isimp) = y(irho)
663  dnz(irho) = dy(irho)
664  IF (rho(irho).NE.0.0d0) THEN
665  fun1(irho) = 1.d0/rho(irho)*(vpr(irho)*nsource(irho,isimp) &
666  + vprm(irho)*nzm1(irho,isimp)/tau &
667  - nz1(irho,isimp)*vpr(irho)*(1.d0/tau))
668  ELSE IF (rho(irho).EQ.0.0d0) THEN
669  fun1(irho) = 4.0d0*itm_pi**2*r0*(nsource(irho,isimp) &
670  + nzm1(irho,isimp)/tau &
671  - nz1(irho,isimp)*(1.d0/tau))
672  END IF
673 
674  END DO rho_loop27
675 
676 
677 
678 
679  CALL integr(nrho,rho,fun1,intfun1) !Integral source
680 
681  rho_loop28: DO irho=1,nrho
682 
683 
684  flux_inter(irho,isimp) = intfun1(irho) &
685  + btprime/2.d0/bt*rho(irho)*vpr(irho)*y(irho)
686 
687 
688  flux(irho,isimp) = vpr(irho)*g1(irho)* &
689  ( y(irho)*vconv(irho,isimp) - dy(irho)*diff(irho,isimp) )
690 
691  END DO rho_loop28
692 
693  END DO imp_loop21
694 
695 
696 
697  CASE(2)
698 
699 ! + + Metoda naprzemiennych kierunkow + +
700 
701 ! + + + + + + + + + + + + + + + + + + + + +
702 
703 
704 ! Definition of Progonka Coefficients - SpaceLoop
705 
706 
707  DO isimp=2,simp+1 ! Impurity State Loop
708 
709  nzp = 0.0
710 
711 ! +++ Boundary conditions for ion diffusion equation in form:
712 !
713 ! V*Y' + U*Y =W
714 !
715 ! +++ On axis:
716  IF(diff(1,isimp).GT.0.d0) THEN
717  v(1) = -diff(1,isimp)
718  u(1) = vconv(1,isimp)
719  w(1) = 0.d0
720  ELSE
721  v(1) = 1.d0
722  u(1) = 0.d0
723  w(1) = 0.d0
724  END IF
725 
726 ! +++ At the edge:
727 ! FIXED NZ
728  IF(nz_bnd_type(isimp).EQ.1) THEN
729  v(2) = 0.d0
730  u(2) = 1.d0
731  w(2) = nz_bnd(1,isimp)
732  ENDIF
733 
734 ! FIXED grad_NZ
735  IF(nz_bnd_type(isimp).EQ.2) THEN
736  v(2) = 1.d0
737  u(2) = 0.d0
738  w(2) = nz_bnd(1,isimp)
739  ENDIF
740 
741 ! FIXED L_NZ
742  IF(nz_bnd_type(isimp).EQ.3) THEN
743  v(2) = nz_bnd(1,isimp)
744  u(2) = 1.d0
745  w(2) = 0.d0
746  ENDIF
747 
748 ! FIXED Flux_NZ
749  IF(nz_bnd_type(isimp).EQ.4) THEN
750  v(2) = -g(nrho)*diff(nrho,isimp)
751  u(2) = g(nrho)*vconv(nrho,isimp)
752  w(2) = nz_bnd(1,isimp)
753  ENDIF
754 
755 ! Generic boundary condition
756  IF(nz_bnd_type(isimp).EQ.5) THEN
757  v(2) = nz_bnd(1,isimp)
758  u(2) = nz_bnd(2,isimp)
759  w(2) = nz_bnd(3,isimp)
760  ENDIF
761 
762 
763 ! For progonka:
764 
765  vp(1) = u(1) - v(1)/(rho(2)-rho(1))
766  up(1) = v(1)/(rho(2)-rho(1))
767  wp(1) = w(1)
768 
769  vp(2) = u(2) + v(2)/(rho(nrho)-rho(nrho-1))
770  up(2) = - v(2)/(rho(nrho)-rho(nrho-1))
771  wp(2) = w(2)
772 
773 
774 
775  ! +++ Coefficients for ion diffusion equation in form:
776 !
777 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = 0
778 
779  rho_loop33: DO irho=1,nrho
780 
781 
782  a(irho) = vpr(irho)
783  b(irho) = vprm(irho)
784  c(irho) = 1.d0
785  d(irho) = vpr(irho)*g1(irho)*diff(irho,isimp)
786  e(irho) = vpr(irho)*g1(irho)*vconv(irho,isimp) &
787  - btprime/2.d0/bt*rho(irho)*vpr(irho)
788  END DO rho_loop33
789  h = tau
790 
791  a1p = 0.0
792  b1p = 0.0
793  c1p = 0.0
794  d1pp = 0.0
795 
796 
797  DO irho=2,nrho-1 ! SPACE LOOP
798 
799  drho = 0.5*(rho(irho+1)-rho(irho-1))
800 
801  drhop = rho(irho+1)-rho(irho)
802 
803  drhom = rho(irho)-rho(irho-1)
804 
805 
806  a1p(irho) = h/drho/c(irho)*( 0.5*(d(irho+1)+d(irho))/drhop ) - &
807  h/drho/c(irho)*( 0.5*(e(irho)+e(irho)) )
808 
809  b1p(irho) = a(irho) + &
810  h/drho/c(irho)*( 0.5*(d(irho+1)+d(irho))/drhop + 0.5*(d(irho)+d(irho-1))/drhom ) - &
811  h/drho/c(irho)*( 0.25*(e(irho+1)-e(irho-1)) )
812 
813  c1p(irho) = h/drho/c(irho)*( 0.5*(d(irho)+d(irho-1))/drhom ) + &
814  h/drho/c(irho)*( 0.5*(e(irho)+e(irho-1)) )
815 
816  d1pp(irho) = b(irho)*nz1(irho,isimp)
817 
818  END DO ! SPACE LOOP
819 
820 ! Finding Ni(NP) - Solotion of Set of Algrbraic Equation by Progonka
821 !
822 ! -Ap*Ni(+1,t+dt) + Bp*Ni(0,t+dt) - Cp*Ni(-1,t+dt) = Dpp
823 
824 
825  CALL unhide_iimp(nrho, simp,isimp, a1p, b1p, c1p, d1pp, vp, up, wp, nzp)
826 
827 
828 ! Przepisanie macierzy
829 
830  nz1(:,isimp) = nzp(:,isimp)
831 
832  END DO !IONS LOOP =========
833 
834 
835 
836  DO irho=2, nrho-1 ! SPACE LOOP
837 
838 
839 ! Definition of Progonka Coefficients - IonizationState
840 
841 
842  nzp = 0.0
843 
844  a2p = 0.0
845  b2p = 0.0
846  c2p = 0.0
847  d2pp = 0.0
848 
849 
850 
851  DO isimp=2,simp+1 ! ION LOOP
852 
853 
854  a2p(isimp) = a(irho)*ne(irho)*beta(irho,isimp+1)*tau
855 
856  b2p(isimp) = a(irho)*ne(irho)*( alfa(irho,isimp)+beta(irho,isimp) )*tau + a(irho)
857 
858  c2p(isimp) = a(irho)*ne(irho)*alfa(irho,isimp-1)*tau
859 
860  d2pp(isimp) = a(irho)*nz1(irho,isimp)
861 
862  END DO ! ION LOOP
863 
864 
865 ! Finding Ni(NP) - Solotion of Set of Algrbraic Equation by Progonka
866 !
867 ! -Ap*Ni(+1,t+dt) + Bp*Ni(0,t+dt) - Cp*Ni(-1,t+dt) = dp
868 
869 
870  CALL unhide_nrho(nrho, simp, irho, a2p, b2p, c2p, d2pp, nz1(irho,1), nzp)
871 
872 ! Rewriting matrix for Next Time Step
873 
874  nz1(irho,:) = nzp(irho,:)
875 
876 
877  ENDDO ! SPACE LOOP
878 
879 
880 
881  END SELECT
882 
883 
884 
885  fluxin = 0.0
886  fluxout = 0.0
887 
888 
889  DO irho=1,nrho
890 
891  densityfluxin(irho) = alfa(irho,1)*vpr(irho)*nz1(irho,1)*ne(irho)-beta(irho,2)*vpr(irho)*nz1(irho,2)*ne(irho)
892 
893  ENDDO
894 
895 
896  fluxin = calkatrapez( densityfluxin, nrho-1, rho(nrho)-rho(nrho-1) )
897 
898  DO isimp=2, simp+1
899 
900  fluxout = fluxout - vpr(nrho)*g1(nrho)*diff(nrho,isimp)*(nz1(nrho,isimp)-nz1(nrho-1,isimp))/(rho(nrho)-rho(nrho-1))
901 
902  END DO
903 
904 
905 
906  !WRITE(*,*) ALFA(2,2),FluxIn, FluxOut, NRHO
907 ! pause
908 
909 ! Flux after "progonka"
910 
911 
912  do isimp=1,simp+1
913  DO irho=1,nrho
914 
915 
916  y(irho)=nz1(irho,isimp)
917 
918  IF (rho(irho).NE.0.0d0) THEN
919  fun1(irho) = 1.d0/rho(irho)*(vpr(irho)*nsource(irho,isimp) &
920  +vprm(irho)*nzm1(irho,isimp)/tau &
921  -nz1(irho,isimp)*vpr(irho)*(1.d0/tau))
922  ELSE IF (rho(irho).EQ.0.0d0) THEN
923  fun1(irho) = 4.d0*itm_pi**2*r0*(nsource(irho,isimp) &
924  +nzm1(irho,isimp)/tau &
925  -nz1(irho,isimp)*(1.d0/tau))
926  END IF
927 
928  END DO
929 
930  CALL integr(nrho,rho,fun1,intfun1) !Integral source
931 
932  CALL derivn(nrho,rho,y,dy)
933 
934  DO irho=1,nrho
935 
936 
937  flux_inter(irho,isimp) = intfun1(irho) &
938  + btprime/2.d0/bt*rho(irho)*vpr(irho)*y(irho)
939 
940 
941  flux(irho,isimp) = vpr(irho)*g1(irho)* &
942  ( y(irho)*vconv(irho,isimp) - dy(irho)*diff(irho,isimp) )
943 
944  END DO
945 
946  END DO
947 
948 
949 ! ImpurityRadiation & Electron source
950 
951  se_exp = 0.0_r8
952 
953  DO irho = 1,nrho
954 
955  DO isimp = 2,simp+1
956 
957 
958  lin_rad1(irho,isimp) = ne(irho)*nz1(irho,isimp)*slin(irho,isimp)
959  brem_rad1(irho,isimp) = ne(irho)*nz1(irho,isimp)*gama(irho,isimp)
960  jon_en1(irho,isimp) = ne(irho)*nz1(irho,isimp)*potential(irho,isimp)*alfa(irho,isimp)
961  rec_los1(irho,isimp) = ne(irho)*nz1(irho,isimp)*potential(irho,isimp-1)*beta(irho,isimp)
962  imp_radiation(irho,isimp) = lin_rad1(irho,isimp) + brem_rad1(irho,isimp) + jon_en1(irho,isimp)*itm_ev &
963  - rec_los1(irho,isimp)*itm_ev
964 
965  ENDDO
966 
967 
968 
969 ! SE_EXP(IRHO) = NZ1(IRHO,1)*NE(IRHO)*ALFA(IRHO,1) !Included in NEUTRAL routine
970  DO isimp = 2,simp
971  se_exp(irho) = se_exp(irho) &
972  + nz1(irho,isimp)*ne(irho)*alfa(irho,isimp) - nz1(irho,isimp)*ne(irho)*beta(irho,isimp)
973  ENDDO
974  se_exp(irho) = se_exp(irho) - nz1(irho,simp+1)*ne(irho)*beta(irho,simp+1)
975 
976 
977  ENDDO
978 
979 
980 
981 ! +++ Deallocate types for interface with numerical solver:
982  CALL deallocate_numerics(solver, ifail)
983 
984 
985  RETURN
986 
987 
988 
989 END SUBROUTINE impurity_one
990 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
991 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
992 
993 
994 
995 
996 
997 
998 
999 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1000 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1001 !-------------------------------------------------------!
1002 ! !
1003 !______________ MATHEMATICAL SUBROUTINES: _____________!
1004 ! !
1005 !-------------------------------------------------------!
1006 ! These subroutines have been extracted from RITM code, !
1007 ! and consist of derivation and integration routines !
1008 !-------------------------------------------------------!
1009 
1010 
1011 
1012 
1013 
1014 
1015 
1016 
1017 
1018 
1019 
1020 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1021 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1022 
1023 
1024 SUBROUTINE unhide_iimp(NRHO, SIMP, ISIMP, As, Bs, Cs, ds, Vbnd, Ubnd, Wbnd, NZPs)
1025 
1026 ! Wyznaczenie Ni(NP) - Solotion of Set of Algrbraic Equation by Progonka:
1027 ! -As*Ni2(-1) + Bs*Ni2(0) - Cs*Ni2(-1) = dp
1028 ! Boundary Condition
1029 ! Ni2(RHO=0,t) = -b1/a2*n(RHO=HRHO(),t)+c1/a1
1030 ! Ni2(L-HRHO(),t) = -aN/bN*n(L,t)+cN/bN
1031 
1032  INTEGER, PARAMETER :: dp = kind(1.0d0)! Double precision
1033 
1034  INTEGER nrho, simp, isimp
1035 
1036  REAL (DP), DIMENSION(NRHO, SIMP+2) :: nzps ! Auxiliary Impurity density ( Position, IonizationState)
1037 
1038  REAL (DP), DIMENSION(NRHO) :: as, bs, cs, ds
1039 
1040  REAL (DP), DIMENSION(2) :: vbnd, ubnd, wbnd
1041 
1042  REAL (DP), DIMENSION(NRHO) :: alfas, betas
1043 
1044  INTEGER irho
1045 
1046 
1047 
1048  REAL (DP) a1_bound, b1_bound, c1_bound
1049 
1050  REAL (DP) an_bound, bn_bound, cn_bound
1051 
1052 
1053  a1_bound = vbnd(1)
1054  b1_bound = ubnd(1)
1055  c1_bound = wbnd(1)
1056 
1057  an_bound = vbnd(2)
1058  bn_bound = ubnd(2)
1059  cn_bound = wbnd(2)
1060 
1061 
1062  alfas(1) = -b1_bound/a1_bound
1063  betas(1) = c1_bound/a1_bound
1064 
1065 
1066 
1067 
1068 
1069  DO irho=2,nrho-1
1070 
1071  alfas(irho) = as(irho)/( bs(irho) - cs(irho)*alfas(irho-1) )
1072 
1073  betas(irho) = ( cs(irho)*betas(irho-1) + ds(irho) )/( bs(irho) - cs(irho)*alfas(irho-1) )
1074 
1075  END DO
1076 
1077 
1078 
1079  nzps(nrho,isimp) = ( cn_bound - betas(nrho-1)*bn_bound )/( an_bound - alfas(nrho-1)*bn_bound )
1080 
1081 
1082  DO irho=nrho-1, 2, -1
1083 
1084  nzps(irho,isimp) = alfas(irho)*nzps(irho+1,isimp) + betas(irho)
1085 
1086  END DO
1087 
1088 
1089 
1090  nzps(1,isimp)=alfas(1)*nzps(2,isimp) + betas(1)
1091 
1092 
1093 END SUBROUTINE unhide_iimp
1094 
1095 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1096 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1097 
1098 
1099 
1100 
1101 
1102 SUBROUTINE unhide_nrho(NRHO, SIMP, IRHO, As, Bs, Cs, ds, Nbnd, NZPs)
1103 
1104 ! Wyznaczenie Ni(NP) - Solotion of Set of Algrbraic Equation by Progonka:
1105 ! -As*Ni2(-1) + Bs*Ni2(0) - Cs*Ni2(-1) = dp
1106 ! Boundary Condition
1107 ! Ni2(RHO=0,t) = -b1/b2*n(RHO=HRHO(),t)+c1/a1
1108 ! Ni2(L-HRHO(),t) = -aN/bN*n(L,t)+cN/bN
1109 
1110 
1111  INTEGER, PARAMETER :: dp = kind(1.0d0)! Double precision
1112 
1113  INTEGER nrho, simp, irho
1114 
1115  REAL (DP) :: nbnd
1116 
1117  REAL (DP), DIMENSION(NRHO,SIMP+2) :: nzps ! Auxiliary Impurity density ( Position, IonizationState)
1118 
1119  REAL (DP), DIMENSION(SIMP+2) :: as, bs, cs, ds
1120 
1121  REAL (DP), DIMENSION(SIMP+2) :: alfas, betas
1122 
1123  INTEGER isimp
1124 
1125 
1126  REAL a1_bound, b1_bound, c1_bound
1127 
1128  REAL an_bound, bn_bound, cn_bound
1129 
1130  a1_bound = 1.0
1131  b1_bound = 0.0
1132  c1_bound = nbnd
1133 
1134 
1135 
1136  an_bound = 1.0
1137  bn_bound = 0.0
1138  cn_bound = 0.0
1139 
1140 
1141  alfas(1) = -b1_bound/a1_bound
1142  betas(1) = c1_bound/a1_bound
1143 
1144 
1145 
1146 
1147 
1148  DO isimp=2,simp+1
1149 
1150  alfas(isimp) = as(isimp)/( bs(isimp) - cs(isimp)*alfas(isimp-1) )
1151 
1152  betas(isimp) = ( cs(isimp)*betas(isimp-1) + ds(isimp) )/( bs(isimp) - cs(isimp)*alfas(isimp-1) )
1153 
1154  END DO
1155 
1156 
1157  nzps(irho,simp+2) = ( cn_bound - betas(simp+1)*bn_bound )/( an_bound - alfas(simp+1)*bn_bound )
1158 
1159 
1160 
1161  DO isimp=simp+1, 2, -1
1162 
1163  nzps(irho,isimp) = alfas(isimp)*nzps(irho,isimp+1) + betas(isimp)
1164 
1165  END DO
1166 
1167 
1168  nzps(irho,1) = alfas(1)*nzps(irho,2) + betas(1)
1169 
1170 
1171 END SUBROUTINE unhide_nrho
1172 
1173 
1174 REAL FUNCTION calkatrapez(wartosc, n, dx)
1175 
1176  INTEGER n ! liczba przedzialow
1177  REAL, DIMENSION(n+1) :: wartosc
1178  REAL (8) :: dx
1179  INTEGER i
1180 
1181  calkatrapez = 0.0
1182  calkatrapez = 0.5*(wartosc(1)+wartosc(n+1))*dx
1183  DO i = 2 , n
1184  calkatrapez = calkatrapez + wartosc(i)*dx
1185  END DO
1186 
1187 END FUNCTION calkatrapez
1188 
1189 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1190 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1191 SUBROUTINE read_data (ALFA_NZ,BETA_NZ,Gama_nz,SLIN_NZ,max_nzimp,simp,potential)
1192 
1193 ! This subroutine reads atom data coeficients for
1194 ! impurity module
1195 
1196 
1197 
1198  INTEGER, PARAMETER :: dp = kind(1.0d0) ! Double precision
1199  INTEGER isimp, jt,simp,max_nzimp
1200 
1201  CHARACTER*70, DIMENSION(3,74) :: atomname
1202 
1203  REAL (DP) :: alfa_nz(max_nzimp+2,500)
1204  REAL (DP) :: beta_nz(max_nzimp+2,500)
1205  REAL (DP) :: gama_nz(max_nzimp+2,500)
1206  REAL (DP) :: slin_nz(max_nzimp+2,500)
1207  REAL (DP) :: potential(max_nzimp+2)
1208 
1209  atomname(1,1)='./src/atomdata/HTABLE'
1210  atomname(2,1)='./src/atomdata/HPROM_OLD'
1211  atomname(3,1)='./src/atomdata/HJON'
1212 
1213  atomname(1,2)='./src/atomdata/HETABLE1'
1214  atomname(2,2)='./src/atomdata/HEPROM1'
1215  atomname(3,2)='./src/atomdata/HEJON'
1216 
1217  atomname(1,3)='./src/atomdata/LITABLE1'
1218  atomname(2,3)='./src/atomdata/LIPROM1'
1219  atomname(3,3)='./src/atomdata/LIJON'
1220 
1221  atomname(1,4)='./src/atomdata/BETABLE1'
1222  atomname(2,4)='./src/atomdata/BEPROM1'
1223  atomname(3,4)='./src/atomdata/BEJON'
1224 
1225  atomname(1,5)='./src/atomdata/BTABLE1'
1226  atomname(2,5)='./src/atomdata/BPROM1'
1227  atomname(3,5)='./src/atomdata/BJON'
1228 
1229  atomname(1,6)='./src/atomdata/CTABLE1'
1230  atomname(2,6)='./src/atomdata/CPROM1'
1231  atomname(3,6)='./src/atomdata/CJON'
1232 
1233 
1234  atomname(1,7)='./src/atomdata/NTABLE1'
1235  atomname(2,7)='./src/atomdata/NPROM1'
1236  atomname(3,7)='./src/atomdata/NJON'
1237 
1238  atomname(1,8)='./src/atomdata/OTABLE1'
1239  atomname(2,8)='./src/atomdata/OPROM1'
1240  atomname(3,8)='./src/atomdata/OJON'
1241 
1242  atomname(1,10)='./src/atomdata/NETABLE1'
1243  atomname(2,10)='./src/atomdata/NEPROM1'
1244  atomname(3,10)='./src/atomdata/NEJON'
1245 
1246  atomname(1,14)='./src/atomdata/SITABLE1'
1247  atomname(2,14)='./src/atomdata/SIPROM1'
1248  atomname(3,14)='./src/atomdata/SIJON'
1249 
1250  atomname(1,18)='./src/atomdata/ARTABLE1'
1251  atomname(2,18)='./src/atomdata/ARPROM1'
1252  atomname(3,18)='./src/atomdata/ARJON'
1253 
1254  atomname(1,26)='./src/atomdata/FETABLE1'
1255  atomname(2,26)='./src/atomdata/FEPROM1'
1256  atomname(3,26)='./src/atomdata/FEJON'
1257 
1258  atomname(1,28)='./src/atomdata/NITABLE1'
1259  atomname(2,28)='./src/atomdata/NIPROM1'
1260 
1261  atomname(1,42)='./src/atomdata/MOTABLE1'
1262  atomname(2,42)='./src/atomdata/MOPROM1'
1263  atomname(3,42)='./src/atomdata/MOJON'
1264 
1265  atomname(1,74)='./src/atomdata/WTABLE1'
1266  atomname(2,74)='./src/atomdata/WPROM1'
1267  atomname(3,74)='./src/atomdata/WJON'
1268 
1269 
1270  alfa_nz = 0.0; beta_nz = 0.0; gama_nz = 0.0
1271  potential = 0.0
1272 
1273 
1274  OPEN(23,file=atomname(1,simp),status='OLD')
1275 
1276 
1277 
1278  DO isimp=1,simp; READ(23,*) (alfa_nz(isimp,jt),beta_nz(isimp+1,jt),gama_nz(isimp+1,jt), jt=1, 500); ENDDO
1279 
1280  CLOSE(23)
1281 
1282  IF (simp.EQ.1.)then
1283  go to 40
1284  endif
1285 
1286  OPEN(33,file=atomname(2,simp),status='OLD')
1287 
1288  READ(33,*) ((slin_nz(isimp,jt),jt=1, 500 ), isimp = 1,simp+1 )
1289 
1290  CLOSE(33)
1291 
1292 40 OPEN(43,file=atomname(3,simp),status='OLD')
1293 
1294  READ(43,*) (potential(isimp), isimp = 1,simp )
1295 
1296  CLOSE(43)
1297 
1298 
1299 END SUBROUTINE read_data
1300 
1301 
solution_interface
subroutine solution_interface(SOLVER, ifail)
INTERFACE TO NUMERICAL SOLVER.
Definition: solution_interface.f90:12
solver
Definition: solver.f90:1
unhide_nrho
subroutine unhide_nrho(NRHO, SIMP, IRHO, As, Bs, Cs, ds, Nbnd, NZPs)
Definition: impurity_one.F90:1102
calkatrapez
REAL function calkatrapez(wartosc, n, dx)
Definition: impurity_one.F90:1174
type_solver
Definition: type_solver.f90:1
derivn
subroutine derivn(N, X, Y, DY1)
These subroutines calculate first and second derivatives, DY1 and DY2, of function Y respect to argum...
Definition: main_plasma.f90:3856
type_solver::allocate_numerics
subroutine allocate_numerics(NDIM, NRHO, SOLVER, ifail)
Definition: type_solver.f90:66
flux
subroutine flux(psitok, rk, zk, nk)
Definition: EQ1_m.f:786
impurity_one
subroutine impurity_one(TE, NE, NZ1, NZM1, VPR, VPRM, R0, BT, BTPRIME, DIFF, FLUX, FLUX_INTER, rho, VCONV, NRHO, SIMP2, NSOURCE, NZ_BND, NZ_BND_TYPE, control_double, CONTROL_INTEGER, G1, IMP_RADIATION, SE_EXP, max_nzimp)
Definition: impurity_one.F90:11
integr
subroutine integr(N, X, Y, INTY)
This subroutine calculates integral of function Y(X)*X from X=0 until X=X(N)
Definition: main_plasma.f90:3916
read_data
subroutine read_data(ALFA_NZ, BETA_NZ, Gama_nz, SLIN_NZ, max_nzimp, simp, potential)
Definition: impurity_one.F90:1191
type_solver::deallocate_numerics
subroutine deallocate_numerics(SOLVER, ifail)
Deallocate plasma profiles needed by the transport solver.
Definition: type_solver.f90:140
unhide_iimp
subroutine unhide_iimp(NRHO, SIMP, ISIMP, As, Bs, Cs, ds, Vbnd, Ubnd, Wbnd, NZPs)
Definition: impurity_one.F90:1024
type_solver::numerics
Definition: type_solver.f90:22