ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
•All Classes Files Functions Variables Pages
main_neutrals.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
2 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
8  SUBROUTINE main_neutrals &
9  (geometry,profiles,neutrals,sources,control,ifail)
10 
11 !-------------------------------------------------------!
12 ! This subroutine calculates densities, fluxes and !
13 ! temperatures of neutrals, source of main ions !
14 ! and sinks of energy from electrons and ions !
15 ! due to interaction with main neutrals4 !
16 !-------------------------------------------------------!
17 ! Source: RITM code !
18 ! Developers: M.Tokar, D.Kalupin !
19 ! Contacts: M.Tokar@fz-juelich.de !
20 ! D.Kalupin@fz-juelich.de !
21 ! !
22 ! Comments: -- !
23 ! !
24 !-------------------------------------------------------!
25 
26 
27 ! +++ Declaration of variables:
28  USE itm_constants
29  USE ets_plasma
30  USE type_neutrals
31  USE type_solver
32 
33 
34  IMPLICIT NONE
35 
36  INTEGER, INTENT (INOUT) :: ifail
37 
38 ! +++ External parameters:
39  TYPE (magnetic_geometry) :: geometry !contains all geometry quantities
40  TYPE (plasma_profiles) :: profiles !contains profiles of plasma parameters
41  TYPE (sources_and_sinks) :: sources !contains profiles of sources
42  TYPE (neutral_profiles) :: neutrals !contains profiles of incoming neutrals
43  TYPE (run_control) :: control !contains all control parameters
44 
45 ! +++ Input/Output to numerical solver:
46  TYPE (numerics) :: solver !contains all I/O quantities to numerics part
47 
48 ! +++ Internal parameters:
49  INTEGER :: nrho !number of radial points (input)
50  INTEGER :: nion !number of ion species (input)
51 
52  INTEGER :: ineu !neutral type index
53  INTEGER :: irho !current radial knot
54  INTEGER :: iion !current ion type
55 
56  REAL (R8) :: rho(profiles%nrho) !normalised minor radius, [m]
57  REAL (R8) :: vpr(profiles%nrho) !V', [m^2]
58 
59  REAL (R8) :: aion !ion mass number [proton mass]
60  REAL (R8) :: ti(profiles%nrho) !ion temperature, [eV]
61  REAL (R8) :: ni(profiles%nrho) !ion density, [m^-3]
62  REAL (R8) :: te(profiles%nrho) !electron temperature, [eV]
63  REAL (R8) :: ne(profiles%nrho) !electron density, [m^-3]
64 
65  REAL (R8) :: t01 !temperature of primary neutrals, [eV]
66  REAL (R8) :: v01 !velocity^2 of primary neutrals, [m^2/s^2]
67  REAL (R8) :: n01(profiles%nrho) !density of primary neutrals, [m^-3]
68  REAL (R8) :: dn01(profiles%nrho) !gradient of primary neutrals density, [m^-4]
69  REAL (R8) :: flux_n01(profiles%nrho) !flux density of primary neutrals, [1/s]
70 
71  INTEGER :: neu_bnd_type !boundary condition type (specified for cold neutrals)
72  REAL (R8) :: neu_bnd(3) !boundary condition value (specified for cold neutrals)
73 
74  REAL (R8) :: t02(profiles%nrho) !temperature of hot neutrals, [eV]
75  REAL (R8) :: v02(profiles%nrho) !velocity^2 of hot neutrals, [m^2/s^2]
76  REAL (R8) :: n02(profiles%nrho) !density of hot neutrals, [m^-3]
77  REAL (R8) :: p02(profiles%nrho) !pressure of hot neutrals, [eV/m^3]
78  REAL (R8) :: dp02(profiles%nrho) !pressure gradient of hot neutrals, [eV/m^4]
79  REAL (R8) :: flux_n02(profiles%nrho) !flux density of hot neutrals, [1/s]
80  REAL (R8) :: source_n02(profiles%nrho) !source of hot neutrals, [1//m^3/s]
81 
82  REAL (R8) :: si(profiles%nrho) !source density of ions, [1/(m^3*s)]
83  REAL (R8) :: qe(profiles%nrho) !energy source/loss density for electrons, [eV/(m^3*s)]
84  REAL (R8) :: qi(profiles%nrho) !energy source/loss density for ions, [eV/(m^3*s)]
85 
86  REAL (R8) :: rec !recycling coefficient of hot neutrals
87  REAL (R8) :: c_ion !ionization rate coefficient [m^3/s]
88  REAL (R8) :: c_cx !charge-exchange rate coefficient [m^3/s]
89  REAL (R8) :: c_rec !recombination rate coefficient [m^3/s]
90  REAL (R8) :: s_rec(profiles%nrho) !recombination source [1/m^3/s]
91  REAL (R8) :: nu_ion(profiles%nrho) !ionization rate coefficient [1/s]
92  REAL (R8) :: nu_cx(profiles%nrho) !charge-exchange rate coefficient [1/s]
93  REAL (R8) :: nu_rec(profiles%nrho) !ionization rate coefficient [1/s]
94  REAL (R8) :: nu_n(profiles%nrho) !ionization rate coefficient [1/s]
95 
96  INTEGER :: iter !iterration index
97  REAL (R8) :: tn02(profiles%nrho) !
98  INTEGER :: i0 !index of R0 point
99  REAL (R8) :: r0 !point with zero flux and T02 = T_i
100  REAL (R8) :: t020 !temperature of neutrals in point R0
101  REAL (R8) :: source_n020 !source of neutrals in point R0
102  REAL (R8) :: hla, tx !
103 
104  REAL (R8) :: amix, tau !mixing factor, time step, [s]
105 
106  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
107  INTEGER :: ndim !number of equations to be solved
108  INTEGER :: solver_type !specifies the option for numerical solution
109  REAL (R8) :: y(profiles%nrho) !function at the current amd previous time steps
110  REAL (R8) :: ym(profiles%nrho) !function at the current amd previous time steps
111  REAL (R8) :: dy(profiles%nrho) !derivative of function
112  REAL (R8) :: a(profiles%nrho) !coefficients for numerical solver
113  REAL (R8) :: b(profiles%nrho) !coefficients for numerical solver
114  REAL (R8) :: c(profiles%nrho) !coefficients for numerical solver
115  REAL (R8) :: d(profiles%nrho) !coefficients for numerical solver
116  REAL (R8) :: e(profiles%nrho) !coefficients for numerical solver
117  REAL (R8) :: f(profiles%nrho) !coefficients for numerical solver
118  REAL (R8) :: g(profiles%nrho) !coefficients for numerical solver
119  REAL (R8) :: h !coefficients for numerical solver
120  REAL (R8) :: v(2), u(2), w(2) !boundary conditions for numerical solver
121 
122  REAL (R8) :: fun1(profiles%nrho), dfun1(profiles%nrho)
123 
124 
125 ! +++ Set up dimensions:
126  nrho = profiles%NRHO
127  nion = profiles%NION
128  ndim = 1
129 
130 
131 ! +++ Allocate types for interface with numerical solver:
132  CALL allocate_numerics(ndim, nrho, solver, ifail)
133 
134 
135 
136 ! +++ Set up local variables:
137  amix = 1.d0
138  solver_type = control%SOLVER_TYPE
139 
140 
141 
142 ! +++ Set up local variables:
143  rho_loop1: DO irho=1,nrho
144  rho(irho) = geometry%RHO(irho)
145  vpr(irho) = geometry%VPR(irho)
146 
147  te(irho) = profiles%TE(irho)
148  ne(irho) = profiles%NE(irho)
149 
150  qe(irho) = 0.d0
151  END DO rho_loop1
152 
153 !!! dpc: HACK
154  IF(rho(1).EQ.0.0_r8 .AND. vpr(1).EQ.0.0_r8) THEN
155  vpr(1) = vpr(2) / 1e6_r8
156  ENDIF
157 
158 
159  rec = neutrals%COEF_RECYCLE
160 
161 
162  ion_loop1: DO iion=1,nion
163  aion = neutrals%MION(iion)
164  t01 = neutrals%T0(nrho,iion,1)
165  n01 = neutrals%N0(nrho,iion,1)
166 
167 
168 
169 ! +++ Set up boundary conditions for particular ion type:
170  neu_bnd_type = neutrals%NEU_BND_TYPE(iion)
171  neu_bnd(1) = neutrals%NEU_BND(1,iion)
172  neu_bnd(2) = neutrals%NEU_BND(2,iion)
173  neu_bnd(3) = neutrals%NEU_BND(3,iion)
174 
175 
176 
177 ! +++ Set up profiles:
178  rho_loop2: DO irho = 1,nrho
179  si(irho) = 0.d0
180  qi(irho) = 0.d0
181 
182  ti(irho) = profiles%TI(irho,iion)
183  ni(irho) = profiles%NI(irho,iion)
184 
185 
186 ! +++ Frequencies of elementary processes:
187  CALL const_atom(te(irho), ti(irho), c_ion, c_cx, c_rec)
188 
189  nu_ion(irho) = c_ion * ne(irho)
190  nu_cx(irho) = c_cx * ni(irho)
191  nu_n(irho) = nu_ion(irho) + nu_cx(irho)
192  s_rec(irho) = c_rec * ne(irho) * ni(irho)
193  END DO rho_loop2
194 
195 
196 
197 !-------------------------------------------------------!
198 ! Primary neutrals !
199 !-------------------------------------------------------!
200 
201  v01 = 9.58d7 * t01 / aion
202 
203 
204 
205 ! +++ Set equation to 'predictive' and all coefficients to zero:
206  flag = 1
207  y(:) = 0.0d0
208  dy(:) = 0.0d0
209  ym(:) = 0.0d0
210  a(:) = 0.0d0
211  b(:) = 0.0d0
212  c(:) = 0.0d0
213  d(:) = 0.0d0
214  e(:) = 0.0d0
215  f(:) = 0.0d0
216  g(:) = 0.0d0
217  h = 0.0d0
218  v(:) = 0.0d0
219  u(:) = 0.0d0
220  w(:) = 0.0d0
221 
222 
223 
224 ! +++ Coefficients for for continuity equation in form:
225 !
226 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
227 
228  rho_loop3: DO irho = 1,nrho
229  y(irho) = n01(irho)
230 
231  c(irho) = vpr(irho)
232  d(irho) = vpr(irho) / nu_n(irho)
233  g(irho) = nu_n(irho) / v01
234  END DO rho_loop3
235  h = 1.d0
236 
237 
238 ! +++ Boundary conditions for ion diffusion equation in form:
239 !
240 ! V*Y' + U*Y =W
241 !
242 ! +++ On axis:
243 ! dN0/drho(rho=0)=0:
244  v(1) = 1.d0
245  u(1) = 0.d0
246  w(1) = 0.d0
247 
248 ! +++ At the edge:
249 ! FIXED N0
250  IF(neu_bnd_type.EQ.1) THEN
251  v(2) = 0.d0
252  u(2) = 1.d0
253  w(2) = neu_bnd(1)
254  END IF
255 
256 
257 ! FIXED FLUX_N0
258  IF(neu_bnd_type.EQ.2) THEN
259  v(2) = 1.d0
260  u(2) = dsqrt(v01)/nu_n(nrho)
261  w(2) = 2.d0 * neu_bnd(1) / dsqrt(v01)
262  END IF
263 
264 ! Generic boundary condition
265  IF(neu_bnd_type.EQ.3) THEN
266  v(2) = neu_bnd(1)
267  u(2) = neu_bnd(2)
268  w(2) = neu_bnd(3)
269  ENDIF
270 
271 
272 
273 ! +++ Defining coefficients for numerical solver:
274  solver%TYPE = solver_type
275  solver%EQ_FLAG(ndim) = flag
276  solver%NDIM = ndim
277  solver%NRHO = nrho
278  solver%AMIX = amix
279 
280 
281  rho_loop4: DO irho=1,nrho
282 
283  solver%RHO(irho) = rho(irho)
284 
285  solver%Y(ndim,irho) = y(irho)
286  solver%DY(ndim,irho) = dy(irho)
287  solver%YM(ndim,irho) = ym(irho)
288 
289  solver%A(ndim,irho) = a(irho)
290  solver%B(ndim,irho) = b(irho)
291  solver%C(ndim,irho) = c(irho)
292  solver%D(ndim,irho) = d(irho)
293  solver%E(ndim,irho) = e(irho)
294  solver%F(ndim,irho) = f(irho)
295  solver%G(ndim,irho) = g(irho)
296 
297  END DO rho_loop4
298 
299  solver%H = h
300 
301  solver%V(ndim,1) = v(1)
302  solver%U(ndim,1) = u(1)
303  solver%W(ndim,1) = w(1)
304  solver%V(ndim,2) = v(2)
305  solver%U(ndim,2) = u(2)
306  solver%W(ndim,2) = w(2)
307 
308 
309 
310 ! +++ Solution of density diffusion equation:
311  CALL solution_interface(solver, ifail)
312 
313 
314 
315  rho_loop5: DO irho=1,nrho
316 ! +++ New solution:
317  y(irho) = solver%Y(ndim,irho)
318  dy(irho) = solver%DY(ndim,irho)
319 
320 
321 
322 ! +++ New profiles of neutral density flux and flux:
323  n01(irho) = y(irho)
324  dn01(irho) = dy(irho)
325 
326 
327  flux_n01(irho) = - dn01(irho) / nu_n(irho) * v01 * vpr(irho)
328 
329 
330 
331 !-------------------------------------------------------!
332 ! Secondary neutrals !
333 !-------------------------------------------------------!
334 
335 
336 ! +++ Source and initial temperature
337  source_n02(irho) = s_rec(irho) + nu_cx(irho) * n01(irho)
338  t02(irho) = ti(irho)
339  END DO rho_loop5
340 
341 
342 
343  iter_loop1: DO iter = 1,5 ! iteration loop for secondary ion temperature T02
344 
345 
346 ! +++ Set equation to 'predictive' and all coefficients to zero:
347  flag = 1
348  y(:) = 0.0d0
349  dy(:) = 0.0d0
350  ym(:) = 0.0d0
351  a(:) = 0.0d0
352  b(:) = 0.0d0
353  c(:) = 0.0d0
354  d(:) = 0.0d0
355  e(:) = 0.0d0
356  f(:) = 0.0d0
357  g(:) = 0.0d0
358  h = 0.0d0
359  v(:) = 0.0d0
360  u(:) = 0.0d0
361  w(:) = 0.0d0
362 
363 
364 
365  rho_loop6: DO irho = 1,nrho
366 
367  v02(irho) = 9.58d7 * t02(irho) / aion
368 
369 
370 ! +++ Coefficients for for continuity equation in form:
371 !
372 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
373 
374  y(irho) = p02(irho)
375 
376  c(irho) = vpr(irho)
377  d(irho) = vpr(irho) / nu_n(irho)
378  f(irho) = source_n02(irho) / v02(irho) * t02(irho)*1.6d-12
379  g(irho) = nu_ion(irho) / v02(irho)
380  END DO rho_loop6
381  h = 1.d0
382 
383 
384 ! +++ Boundary conditions:
385 ! +++ On axis:
386 ! dP0/drho(rho=0)=0:
387  v(1) = 1.d0
388  u(1) = 0.d0
389  w(1) = 0.d0
390 
391 ! +++ At the edge:
392  v(2) = (1.d0+rec) / nu_n(nrho)
393  u(2) = (1.d0-rec) / dsqrt(v02(nrho))
394  w(2) = 0.d0
395 
396 
397 
398 
399 ! +++ Defining coefficients for numerical solver:
400  solver%TYPE = solver_type
401  solver%EQ_FLAG(ndim) = flag
402  solver%NDIM = ndim
403  solver%NRHO = nrho
404  solver%AMIX = amix
405 
406  rho_loop7: DO irho=1,nrho
407 
408  solver%RHO(irho) = rho(irho)
409 
410  solver%Y(ndim,irho) = y(irho)
411  solver%DY(ndim,irho) = dy(irho)
412  solver%YM(ndim,irho) = ym(irho)
413 
414  solver%A(ndim,irho) = a(irho)
415  solver%B(ndim,irho) = b(irho)
416  solver%C(ndim,irho) = c(irho)
417  solver%D(ndim,irho) = d(irho)
418  solver%E(ndim,irho) = e(irho)
419  solver%F(ndim,irho) = f(irho)
420  solver%G(ndim,irho) = g(irho)
421 
422  END DO rho_loop7
423 
424  solver%H = h
425 
426  solver%V(ndim,1) = v(1)
427  solver%U(ndim,1) = u(1)
428  solver%W(ndim,1) = w(1)
429  solver%V(ndim,2) = v(2)
430  solver%U(ndim,2) = u(2)
431  solver%W(ndim,2) = w(2)
432 
433 
434 
435 ! +++ Solution of density diffusion equation:
436  CALL solution_interface(solver, ifail)
437 
438 
439 
440  rho_loop8: DO irho=1,nrho
441 ! +++ New solution:
442  y(irho) = solver%Y(ndim,irho)
443  dy(irho) = solver%DY(ndim,irho)
444 
445 
446 ! +++ New neutral pressure and its derivative:
447  p02(irho) = y(irho)
448  dp02(irho) = dy(irho)
449 
450 
451 ! +++ New neutral density and flux:
452  n02(irho) = p02(irho) / (1.6d-12*t02(irho))
453  flux_n02(irho) = - dp02(irho) * v02(irho) / (1.6d-12*t02(irho)) &
454  / nu_n(irho) * vpr(irho)
455  END DO rho_loop8
456 
457 !!! DPC: HACK
458  flux_n02(1) = 0.0_r8
459 
460 ! +++ temperature
461 ! point R0 with zero flux and T02 = T_i
462  rho_loop9: DO irho = nrho, 2, -1
463  i0=irho
464  IF (flux_n02(irho) * flux_n02(irho-1).LE.0.d0) goto 10
465  END DO rho_loop9
466 
467  10 CONTINUE
468  r0 = ( flux_n02(i0)/vpr(i0)*rho(i0-1) &
469  - flux_n02(i0-1)/vpr(i0-1)*rho(i0) ) &
470  / (flux_n02(i0)/vpr(i0) - flux_n02(i0-1)/vpr(i0-1))
471  t020 = ( ti(i0) * (r0 - rho(i0-1)) &
472  + ti(i0-1)*(rho(i0)-r0)) &
473  / (rho(i0)-rho(i0-1))
474  source_n020 = (source_n02(i0)*(r0-rho(i0-1)) &
475  + source_n02(i0-1) * (rho(i0)-r0)) &
476  / (rho(i0)-rho(i0-1))
477 
478 
479 
480 ! T02 towards edge from R0
481  hla = (r0-rho(i0)) * (source_n020+source_n02(i0)) &
482  / flux_n02(i0) * vpr(i0)
483  tx = (t020+ti(i0)) / 2.d0
484 
485  CALL step_t02(hla, tx, t020, tn02(i0))
486 
487  rho_loop10: DO irho = i0,nrho-1
488  hla = (rho(irho)-rho(irho+1)) &
489  * (source_n02(irho)+source_n02(irho+1)) &
490  / (flux_n02(irho)/vpr(irho)+flux_n02(irho+1)/vpr(irho+1))
491  tx = (ti(irho)+ti(irho+1)) / 2.d0
492 
493  CALL step_t02(hla, tx, tn02(irho),tn02(irho+1))
494 
495  END DO rho_loop10
496 
497 
498 
499 ! T02 towards core from R0
500 !!! DPC: another hack
501  IF(i0.GT.2) THEN
502  hla = (r0-rho(i0-1)) * (source_n020+source_n02(i0-1)) &
503  / flux_n02(i0-1) * vpr(i0-1)
504  tx = (t020+ti(i0-1)) / 2.d0
505 
506  CALL step_t02(hla, tx, t020, tn02(i0-1))
507 
508  rho_loop11: DO irho = i0-1, 2, -1
509  hla = (rho(irho)-rho(irho-1)) &
510  * (source_n02(irho)+source_n02(irho-1)) &
511  / (flux_n02(irho)/vpr(irho)+flux_n02(irho-1)/vpr(irho-1))
512  tx = (ti(irho)+ti(irho-1))/2.d0
513 
514  tn02(irho-1) = ti(irho-1)
515 
516  IF (hla.LT.0.d0) CALL step_t02(hla, tx, tn02(irho), tn02(irho-1))
517 
518  END DO rho_loop11
519  ENDIF
520 
521 
522 ! new approximation for T02
523  amix = 1.d0
524 
525  rho_loop12: DO irho = 2, nrho
526  t02(irho) = (1.d0-amix)*t02(irho) + amix*tn02(irho)
527  END DO rho_loop12
528  t02(1) = t02(2)
529 
530  END DO iter_loop1
531 
532 
533 
534  rho_loop13: DO irho = 1, nrho
535 ! +++ Particle and heat sources due to neutrals:
536  si(irho) = (n01(irho)+n02(irho)) * nu_ion(irho) - s_rec(irho)
537  qe(irho) = qe(irho) - 30.d0 * (n01(irho)+n02(irho)) &
538  * nu_ion(irho) - 1.5d0*te(irho) * s_rec(irho)
539  qi(irho) = 1.5d0 * (t01*n01(irho) + t02(irho)*n02(irho)) &
540  * nu_n(irho) - 1.5d0*ti(irho) &
541  * ((n01(irho)+n02(irho)) * nu_cx(irho) + s_rec(irho))
542 
543 
544 ! +++ Output neutral particle and heat sources to the work flow:
545  neutrals%N0(irho,iion,1) = n01(irho)
546  neutrals%FLUX_N0(irho,iion,1) = flux_n01(irho)
547  neutrals%T0(irho,iion,2) = t01
548 
549  neutrals%N0(irho,iion,2) = n02(irho)
550  neutrals%FLUX_N0(irho,iion,2) = flux_n02(irho)
551  neutrals%T0(irho,iion,2) = t02(irho)
552 
553 
554  sources%SI_EXP(irho,iion) = si(irho)
555  sources%QI_EXP(irho,iion) = qi(irho) * itm_ev
556  sources%QE_EXP(irho) = qe(irho) * itm_ev
557 
558  END DO rho_loop13
559 
560 
561 
562  END DO ion_loop1
563 
564 
565  RETURN
566 
567 
568 
569  END SUBROUTINE main_neutrals
570 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
571 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
572 
573 
574 
575 
576 
577 
578 
579 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
580 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
586  SUBROUTINE step_t02(FUN, TX, T0, T1)
587 
588  USE itm_constants
589 
590  IMPLICIT NONE
591 
592  REAL (R8) :: fun, tx, t0, t1
593 
594 
595  t1 = tx + (t0-tx) * dexp(fun)
596 
597 
598  RETURN
599 
600 
601  END SUBROUTINE step_t02
602 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
603 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
604 
605 
606 
607 
608 
609 
610 
611 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
612 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
618  SUBROUTINE const_atom (TE, TI, C_ION, C_CX, C_REC)
619 !-------------------------------------------------------!
620 ! The temperature dependence of the rate !
621 ! coefficients for atom ionization, !
622 ! charge-exchange and recombination !
623 !-------------------------------------------------------!
624 
625  USE itm_constants
626 
627  IMPLICIT NONE
628 
629  REAL (R8) :: te, ti ! eectron and ion temperatur [eV]
630  REAL (R8) :: c_ion ! ionization rate coefficient [m^3/s]
631  REAL (R8) :: c_cx ! charge-exchange rate coefficient [m^3/s]
632  REAL (R8) :: c_rec ! ionization rate coefficient [m^3/s]
633 
634  REAL (R8) :: fun1, fun2
635 
636 
637  c_ion = 0.73d-14 * dsqrt(te) * dexp(-13.6d0/te) &
638  / (1.d0+0.01d0*te)
639  c_cx = 1.d-14*ti**0.3d0
640  fun1 = 13.6d0/te
641  fun2 = 0.d0
642 
643  IF (fun1.GT.0.d0.AND.fun1.LE.1.d0) &
644  fun2 = fun1 *dexp(-fun1) * ( -0.58d0 -dlog(fun1) &
645  +fun1 -0.25d0*fun1**2 +0.055d0*fun1**3 )
646 
647  IF (fun1.GT.1.d0) &
648  fun2 = (fun1**2 +2.33d0*fun1 +0.25d0) &
649  / (fun1**2 +3.33d0*fun1 +1.68d0)
650 
651  c_rec = 5.2d-20*dsqrt(fun1)*fun2
652 
653 
654  RETURN
655 
656 
657  END SUBROUTINE const_atom
658 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
659 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
660 
661 
662 
663 
solution_interface
subroutine solution_interface(SOLVER, ifail)
INTERFACE TO NUMERICAL SOLVER.
Definition: solution_interface.f90:12
profiles
subroutine profiles(p0, rbphi, dp0, drbphi, a)
Definition: profiles.f90:1
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
fun
subroutine fun(X, F)
Definition: Ev2.f:10
ets_plasma::sources_and_sinks
Definition: ets_plasma.f90:294
const_atom
subroutine const_atom(TE, TI, C_ION, C_CX, C_REC)
??
Definition: main_neutrals.f90:618
type_solver::allocate_numerics
subroutine allocate_numerics(NDIM, NRHO, SOLVER, ifail)
Definition: type_solver.f90:66
step_t02
subroutine step_t02(FUN, TX, T0, T1)
??
Definition: main_neutrals.f90:586
ets_plasma::plasma_profiles
Definition: ets_plasma.f90:80
ets_plasma::run_control
Definition: ets_plasma.f90:426
type_neutrals::neutral_profiles
Definition: type_neutrals.f90:32
type_solver::numerics
Definition: type_solver.f90:22
ets_plasma
The module declares types of variables used in ETS (transport code)
Definition: ets_plasma.f90:8
main_neutrals
subroutine main_neutrals(GEOMETRY, PROFILES, NEUTRALS, SOURCES, CONTROL, ifail)
??
Definition: main_neutrals.f90:8
ets_plasma::magnetic_geometry
Definition: ets_plasma.f90:44
type_neutrals
??
Definition: type_neutrals.f90:9