ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
neutrals.F90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
2 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3  MODULE neutrals
4 
5  CONTAINS
6 
7 
8 
9 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
10 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
11 
12  SUBROUTINE neutrals_ets (COREIMPUR_ITER, EQUILIBRIUM_ITER, COREPROF_ITER, &
13  coreneutrals_old, coreneutrals_iter, &
14  coresource_new, coreneutrals_new, &
15  control_integer, control_double)
16 
17 
18 
19 !-------------------------------------------------------!
20 ! This routine provides interface between CPO !
21 ! derived types and internal ETS derived types. !
22 !-------------------------------------------------------!
23 ! Source: --- !
24 ! Developers: R. Stankiewicz, I. Ivanova-Stanik !
25 ! Contacts: romsta@ifpilm.waw.pl !
26 ! !
27 ! Comments: input data for NEUTRALS MODULE !
28 ! !
29 ! !
30 ! D.Kalupin: Some corrections are introduced !
31 !-------------------------------------------------------!
32 
33 
34  USE euitm_schemas
35  USE ets_plasma
36  USE type_solver
37  USE copy_structures
38  USE deallocate_structures
39  USE itm_types
40  USE itm_constants
41  USE allocate_deallocate
42 
43 #ifdef GOT_AMNSPROTO
44  USE amns_types
45  USE amns_module
46  USE euitm_routines
47 #endif
48 
49 
50  IMPLICIT NONE
51 
52 ! +++ Input/Output to numerical solver:
53  TYPE (numerics) :: solver !contains all I/O quantities to numerics part
54 
55 ! +++ CPO derived types:
56  TYPE (type_coreneutrals),POINTER :: coreneutrals_old(:) !input CPO with neutrals
57  TYPE (type_coreneutrals),POINTER :: coreneutrals_iter(:) !input CPO with neutrals
58  TYPE (type_coreneutrals),POINTER :: coreneutrals_new(:) !input CPO with neutrals
59  TYPE (type_equilibrium), POINTER :: equilibrium_iter(:) !input CPO with geometry quantities from previous ITERration
60  TYPE (type_coreprof), POINTER :: coreprof_iter(:) !input CPO with plasma profiles
61  TYPE (type_coreimpur), POINTER :: coreimpur_iter(:) !input CPO with impurities
62  TYPE (type_coresource), POINTER :: coresource_new(:) !output CPO for neutral radiation
63 
64 
65 ! +++ Dimensions:
66  INTEGER :: neq !number of radial points (input, determined from EQUILIBRIUM CPO)
67  INTEGER :: nrho !number of radial points (input, determined from COREPROF CPO)
68  INTEGER :: nion !number of ion species (input, determined from COREPROF CPO)
69  INTEGER :: nimp !number of impurity species (input, determined from COREIMPUR CPO)
70  INTEGER, ALLOCATABLE :: nzimp(:)
71  INTEGER, ALLOCATABLE :: nn_bnd_type(:,:) !boundary condition, type, one neutral
72  INTEGER :: nneut !number of neutrals species (input)
73  INTEGER :: nnucl !number of neutrals species (input)
74  INTEGER, ALLOCATABLE :: ntype(:) !number of impurity ionization states (input)
75  INTEGER, ALLOCATABLE :: ncomp(:) !max_number of distinct atoms enter the composition-"1" wich is neutral
76  INTEGER, PARAMETER :: nocur = 1 !number of CPO ocurancies in the work flow
77  INTEGER :: zmax, tmax
78 
79 ! +++ Indexes:
80  INTEGER :: irho !
81  INTEGER :: iion !
82  INTEGER :: iimp !
83  INTEGER :: icomp !
84  INTEGER :: itype !
85  INTEGER :: ineut !
86  INTEGER :: inucl !
87  INTEGER :: iz !
88  INTEGER :: hot_neutrals, cold_neutrals
89 
90 
91 
92  REAL (R8), ALLOCATABLE :: rho(:) !toroidal flux coordinate,not normalise,[m]
93  REAL (R8), ALLOCATABLE :: aneut(:) !masa of the neutral
94  REAL (R8), ALLOCATABLE :: vpr(:) !V', [m^2]
95  REAL (R8), ALLOCATABLE :: vprm(:) !V' (at previous time step), [m^2]
96  REAL (R8), ALLOCATABLE :: g3(:)
97  REAL (R8), ALLOCATABLE :: ne(:) !electron density [m^-3]
98  REAL (R8), ALLOCATABLE :: te(:) !electron temperature
99  REAL (R8), ALLOCATABLE :: ni(:,:) !ion density of "1"-ionization state [m^-3]
100  REAL (R8), ALLOCATABLE :: ti(:,:) !ion temperature of "1"-ionization state
101  REAL (R8), ALLOCATABLE :: dn0(:,:,:) !density gradient, [m^-4]
102  REAL (R8), ALLOCATABLE :: flux(:,:,:) !ion flux, [1/s]
103  REAL (R8), ALLOCATABLE :: flux_inter(:,:,:) !ion flux, [1/s]
104  REAL (R8), ALLOCATABLE :: n0(:,:,:) !one neutral density
105  REAL (R8), ALLOCATABLE :: n0m(:,:,:) !old one neutral density
106  REAL (R8), ALLOCATABLE :: t0(:,:,:) !temperature one neutrals, [eV]
107  REAL (R8), ALLOCATABLE :: vt(:,:,:) !toroidal velosicty for one neutrals,
108  REAL (R8), ALLOCATABLE :: vr(:,:,:) !radial velosicty for one neutrals
109  REAL (R8), ALLOCATABLE :: vp(:,:,:) !poloidal velosicty for one neutrals
110  REAL (R8), ALLOCATABLE :: vconv(:,:,:) !pinch velocity for different ionisation [m/s]
111  REAL (R8), ALLOCATABLE :: diff(:,:,:)
112  REAL (R8), ALLOCATABLE :: nn_bnd(:,:,:) !boundary condition, value, one neutral[depends on NZ_BND_TYPE]
113  REAL (R8), ALLOCATABLE :: nnsource(:,:,:) !value of the source term,[m^-3.s^-1]
114  REAL (R8) :: amix, tau !mixing factor, time step, [s]
115 
116 ! dimension for solver
117  INTEGER, INTENT(IN) :: control_integer(1) !integer control parameters
118  REAL (R8), INTENT(IN) :: control_double(2) !real control parameters
119  INTEGER :: solut_method
120  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
121  INTEGER :: ndim !number of equations to be solved
122  INTEGER :: solver_type !specifies the option for numerical solution
123  REAL (R8), ALLOCATABLE :: y(:) !function at the current amd previous time steps
124  REAL (R8), ALLOCATABLE :: ym(:) !function at the current amd previous time steps
125  REAL (R8), ALLOCATABLE :: dy(:) !derivative of function
126  REAL (R8), ALLOCATABLE :: a(:) !coefficients for numerical solver
127  REAL (R8), ALLOCATABLE :: b(:) !coefficients for numerical solver
128  REAL (R8), ALLOCATABLE :: c(:) !coefficients for numerical solver
129  REAL (R8), ALLOCATABLE :: d(:) !coefficients for numerical solver
130  REAL (R8), ALLOCATABLE :: e(:) !coefficients for numerical solver
131  REAL (R8), ALLOCATABLE :: f(:) !coefficients for numerical solver
132  REAL (R8), ALLOCATABLE :: g(:) !coefficients for numerical solver
133  REAL (R8) :: h !coefficients for numerical solver
134  REAL (R8) :: v(2), u(2), w(2) !boundary conditions for numerical solver
135  REAL (R8), ALLOCATABLE :: fun(:), intfun(:)
136 
137 
138 ! +++ parameters for atomic data
139  REAL (R8), ALLOCATABLE :: ionizat(:,:) !atom data: ionization coefficient (after interpolation) [m^3/s]
140  REAL (R8), ALLOCATABLE :: recomb(:,:) !atom data: ionization coefficient (after interpolation) [m^3/s]
141  REAL (R8), ALLOCATABLE :: potential(:,:) !atom data: potential on jonisation [eV]
142  REAL (R8), ALLOCATABLE :: chargeexch(:,:) !atom data: charge-exchange rate coefficient [m^3/s]
143  REAL (R8), ALLOCATABLE :: si_exp(:,:)
144  REAL (R8), ALLOCATABLE :: ssi_exp(:,:)
145  REAL (R8), ALLOCATABLE :: qi_exp(:,:)
146  REAL (R8), ALLOCATABLE :: sz_exp(:,:,:)
147  REAL (R8), ALLOCATABLE :: qz_exp(:,:,:)
148  REAL (R8), ALLOCATABLE :: qe_exp(:)
149  REAL (R8), ALLOCATABLE :: qe_exp_tot(:)
150 
151  REAL (R8) :: time !Time
152 
153  LOGICAL, SAVE :: first = .true.
154  INTEGER :: ifail
155 
156 
157 #ifdef GOT_AMNSPROTO
158  TYPE (amns_handle_type), SAVE :: amns
159  TYPE (amns_handle_rx_type), ALLOCATABLE, SAVE :: amns_ei(:,:), amns_eip(:,:), amns_rc(:,:)
160  TYPE (amns_handle_rx_type), ALLOCATABLE, SAVE :: amns_cx(:,:), amns_lr(:,:), amns_br(:,:)
161  TYPE (amns_version_type) :: amns_database
162  TYPE (amns_reaction_type) :: ei_rx, eip_rx, rc_rx, cx_rx, lr_rx, br_rx
163  TYPE (amns_reactants_type) :: species
164  TYPE (amns_query_type) :: query
165  TYPE (amns_answer_type) :: answer
166  TYPE (amns_set_type) :: set
167  REAL (R8) :: zn_neut, mi_neut
168 #endif
169  CHARACTER (len=80) :: format
170 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
171 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
172 
173 
174 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
175 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
176 ! +++ Set dimensions:
177  neq = SIZE (equilibrium_iter(1)%profiles_1d%rho_tor )
178  nrho = SIZE (coreneutrals_iter(1)%rho_tor )
179  CALL get_comp_dimensions(coreneutrals_iter(1)%COMPOSITIONS, nnucl, nion, nimp, nzimp, nneut, ntype, ncomp)
180  ndim = 1
181 
182 
183 ! Maximum charge state of Ions and Impurities:
184 ! ZMAX = 0
185 ! DO IION = 1,NION
186 ! ZMAX = MAX(INT(CORENEUTRALS_ITER(1)%COMPOSITIONS%IONS(IION)%ZION),ZMAX)
187 ! END DO
188 ! ZMAX = MAX(MAXVAL(NZIMP),ZMAX)
189  tmax = maxval(ntype)
190 
191  zmax = 1
192 
193 
194 
195 
196 ! +++ Allocate output CPOs:
197  CALL allocate_coreneutrals_cpo(nocur, nrho, nnucl, nion, nimp, nzimp, nneut, ntype, ncomp, coreneutrals_new)
198  CALL copy_cpo(coreneutrals_iter(1), coreneutrals_new(1))
199  CALL allocate_coresource_cpo(nocur, nrho, nnucl, nion, nimp, nzimp, nneut, ntype, ncomp, coresource_new)
200  CALL deallocate_cpo(coresource_new(1)%COMPOSITIONS)
201  CALL copy_cpo(coreneutrals_iter(1)%COMPOSITIONS, coresource_new(1)%COMPOSITIONS)
202 
203 
204 
205 
206 ! +++ Copy time and grids to output CPOs:
207  coreneutrals_new(1)%datainfo%cocos = 13
208  coreneutrals_new(1)%time = coreprof_iter(1)%time
209  coreneutrals_new(1)%rho_tor = coreprof_iter(1)%rho_tor
210  coreneutrals_new(1)%rho_tor_norm = coreprof_iter(1)%rho_tor_norm
211 
212  coresource_new(1)%datainfo%cocos = 13
213  coresource_new(1)%time = coreprof_iter(1)%time
214  coresource_new(1)%VALUES(1)%rho_tor = coreprof_iter(1)%rho_tor
215  coresource_new(1)%VALUES(1)%rho_tor_norm = coreprof_iter(1)%rho_tor_norm
216 
217 
218 
219 
220 
221 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
222 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
223 ! +++ Allocate types for interface with numerical solver:
224  CALL allocate_numerics(ndim, nrho, solver, ifail)
225 
226 
227 ! +++ Allocate internal variables:
228  ALLOCATE ( rho(nrho) )
229  ALLOCATE ( vpr(nrho) )
230  ALLOCATE ( vprm(nrho) )
231  ALLOCATE ( g3(nrho) )
232  ALLOCATE ( ne(nrho) )
233  ALLOCATE ( te(nrho) )
234  ALLOCATE ( ti(nrho,nion) )
235  ALLOCATE ( ni(nrho,nneut) )
236 
237  ALLOCATE ( qe_exp(nrho) )
238  ALLOCATE ( qe_exp_tot(nrho) )
239  ALLOCATE ( si_exp(nrho,nion) )
240  ALLOCATE ( qi_exp(nrho,nion) )
241  ALLOCATE ( sz_exp(nrho,nimp,zmax) )
242  ALLOCATE ( qz_exp(nrho,nimp,zmax) )
243  ALLOCATE ( ssi_exp(nrho,nion) )
244 
245  ALLOCATE ( y(nrho) )
246  ALLOCATE ( ym(nrho) )
247  ALLOCATE ( dy(nrho) )
248  ALLOCATE ( a(nrho) )
249  ALLOCATE ( b(nrho) )
250  ALLOCATE ( c(nrho) )
251  ALLOCATE ( d(nrho) )
252  ALLOCATE ( e(nrho) )
253  ALLOCATE ( f(nrho) )
254  ALLOCATE ( g(nrho) )
255 
256  ALLOCATE ( fun(nrho) )
257  ALLOCATE ( intfun(nrho) )
258 
259  ALLOCATE ( aneut(nneut) )
260 
261  ALLOCATE ( n0(nrho,nneut,tmax) )
262  ALLOCATE ( n0m(nrho,nneut,tmax) )
263  ALLOCATE ( vconv(nrho,nneut,tmax) )
264  ALLOCATE ( diff(nrho,nneut,tmax) )
265  ALLOCATE ( dn0(nrho,nneut,tmax) )
266  ALLOCATE ( flux(nrho,nneut,tmax) )
267  ALLOCATE ( nn_bnd(3, nneut,tmax) )
268  ALLOCATE (nn_bnd_type( nneut,tmax) )
269  ALLOCATE ( t0(nrho,nneut,tmax) )
270  ALLOCATE ( vt(nrho,nneut,tmax) )
271  ALLOCATE ( vp(nrho,nneut,tmax) )
272  ALLOCATE ( vr(nrho,nneut,tmax) )
273  ALLOCATE ( nnsource(nrho,nneut,tmax) )
274 
275  ALLOCATE ( ionizat(nrho,0:zmax) )
276  ALLOCATE ( recomb(nrho,0:zmax) )
277  ALLOCATE ( potential(nrho,0:zmax) )
278  ALLOCATE ( chargeexch(nrho,0:zmax) )
279 
280 
281 
282 
283 ! +++ Fill internal variables:
284  si_exp = 0._r8
285  ssi_exp = 0._r8
286  sz_exp = 0._r8
287  qi_exp = 0._r8
288  qe_exp = 0._r8
289  qe_exp_tot = 0._r8
290 
291  n0 = 0._r8
292  n0m = 0._r8
293  dn0 = 0._r8
294  vconv = 0._r8
295  nn_bnd = 0._r8
296  nn_bnd_type = 0
297  t0 = 0._r8
298  vt = 0._r8
299  vr = 0._r8
300  vp = 0._r8
301  nnsource = 0._r8
302  diff = 0._r8
303 
304  rho = coreprof_iter(1)%RHO_TOR
305  ne = coreprof_iter(1)%NE%VALUE
306  te = coreprof_iter(1)%TE%VALUE
307 
308  DO iion = 1,nion
309  ti(:,iion) = coreprof_iter(1)%TI%VALUE(:,iion)
310  ENDDO
311 
312  amix = control_double(2)
313  tau = control_double(1)
314  solver_type = control_integer(1)
315 
316 
317  CALL l3deriv(equilibrium_iter(1)%profiles_1d%volume, equilibrium_iter(1)%profiles_1d%rho_tor, neq, &
318  vpr, rho, nrho)
319 
320  CALL l3deriv(equilibrium_iter(1)%profiles_1d%volume, equilibrium_iter(1)%profiles_1d%rho_tor, neq, &
321  vprm, rho, nrho)
322 
323  CALL l3interp(equilibrium_iter(1)%profiles_1d%gm3, equilibrium_iter(1)%profiles_1d%rho_tor, neq, &
324  g3, rho, nrho)
325 
326  DO ineut=1, nneut
327  inucl = coreneutrals_old(1)%compositions%NEUTRALSCOMP(ineut)%NEUTCOMP(1)%nucindex
328  aneut(ineut) = coreneutrals_old(1)%compositions%nuclei(inucl)%amn
329 
330 ! CORENEUTRALS_ITER(1)%PROFILES(INEUT)%NEUTRALTYPE(2)%t0%value(:) = COREPROF_ITER(1)%TI%VALUE(:,1)
331  coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(1)%t0%value(:) = coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(1)%t0%boundary%value(1)
332 
333  DO itype= 1, ntype(ineut)
334  CALL l3interp(coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%n0%value, coreneutrals_iter(1)%rho_tor, SIZE(coreneutrals_iter(1)%rho_tor), &
335  n0(:,ineut,itype), rho, nrho)
336  CALL l3interp(coreneutrals_old(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%n0%value, coreneutrals_iter(1)%rho_tor, SIZE(coreneutrals_iter(1)%rho_tor), &
337  n0m(:,ineut,itype), rho, nrho)
338  CALL l3interp(coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%t0%value, coreneutrals_iter(1)%rho_tor, SIZE(coreneutrals_iter(1)%rho_tor), &
339  t0(:,ineut,itype), rho, nrho)
340  CALL l3interp(coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%v0%toroidal%value, coreneutrals_iter(1)%rho_tor, SIZE(coreneutrals_iter(1)%rho_tor), &
341  vt(:,ineut,itype), rho, nrho)
342  CALL l3interp(coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%v0%poloidal%value, coreneutrals_iter(1)%rho_tor, SIZE(coreneutrals_iter(1)%rho_tor), &
343  vp(:,ineut,itype), rho, nrho)
344  CALL l3interp(coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%v0%radial%value, coreneutrals_iter(1)%rho_tor, SIZE(coreneutrals_iter(1)%rho_tor), &
345  vr(:,ineut,itype), rho, nrho)
346 
347  nn_bnd(:,ineut,itype) = coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%n0%boundary%value(:)
348  nn_bnd_type(ineut,itype) = coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%n0%boundary%type
349  END DO
350 
351 
352  DO iion=1,nion
353  IF(abs(aneut(ineut) - coreprof_iter(1)%compositions%nuclei(coreprof_iter(1)%compositions%IONS(iion)%nucindex)%amn).LE. 0.25) THEN
354  ni(:,ineut) = coreprof_iter(1)%NI%VALUE(:,iion)
355  coreneutrals_iter(1)%PROFILES(ineut)%NEUTRALTYPE(2)%t0%value(:) = coreprof_iter(1)%TI%VALUE(:,iion)
356  ENDIF
357  ENDDO
358 
359  DO iimp=1,nimp
360  IF(abs(aneut(ineut) - coreimpur_iter(1)%compositions%nuclei(coreimpur_iter(1)%compositions%IMPURITIES(iimp)%nucindex)%amn).LE. 0.25) THEN
361  CALL l3interp(coreimpur_iter(1)%IMPURITY(iimp)%NZ(:,1), coreimpur_iter(1)%rho_tor, SIZE(coreimpur_iter(1)%rho_tor), &
362  ni(:,ineut), rho, nrho)
363  ENDIF
364  END DO
365  END DO
366 
367 
368 
369 
370 
371 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
372 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
373 ! +++ Build atomic data:
374  IF(first) THEN
375 #ifdef GOT_AMNSPROTO
376  WRITE(*,*) 'ITM AMNSPROTO data used (via UAL)'
377  ALLOCATE(amns_ei(0:zmax, nneut), amns_rc(0:zmax, nneut), &
378  amns_eip(0:zmax,nneut), amns_lr(0:zmax, nneut), &
379  amns_br(0:zmax, nneut), amns_cx(0:zmax, nneut))
380  CALL itm_amns_setup(amns)
381  query%string = 'version'
382  CALL itm_amns_query(amns,query,answer)
383  WRITE(*,*) 'AMNS data base version = ',trim(answer%string)
384  ei_rx%string = 'EI'
385  eip_rx%string = 'EIP'
386  rc_rx%string = 'RC'
387  lr_rx%string = 'LR'
388  br_rx%string = 'BR'
389  cx_rx%string = 'CX'
390  FORMAT = '(''ZN = '',f5.2,'', IS = '',i2,'', RX = '',a,'', SRC = '',a)'
391  query%string = 'source'
392 
393  DO ineut=1, nneut
394 
395  inucl = coreneutrals_old(1)%compositions%NEUTRALSCOMP(ineut)%NEUTCOMP(1)%nucindex
396  zn_neut = coreneutrals_old(1)%compositions%nuclei(inucl)%zn
397  mi_neut = coreneutrals_old(1)%compositions%nuclei(inucl)%amn
398  mi_neut = 0
399 
400  DO iz = 0, zmax-1
401  IF (zn_neut .GE. iz) THEN
402 ! EI
403  allocate(species%components(4))
404  species%components = &
405  (/ amns_reactant_type(zn_neut, iz, mi_neut, 0), &
406  amns_reactant_type(0, -1, 0, 0), &
407  amns_reactant_type(zn_neut, iz+1, mi_neut, 1), &
408  amns_reactant_type(0, -1, 0, 1) &
409  /)
410  CALL itm_amns_setup_table(amns, ei_rx, species, amns_ei(iz, ineut))
411  deallocate(species%components)
412  END IF
413  END DO
414 
415 
416  DO iz = 0, zmax
417  IF (zn_neut .GE. iz) THEN
418 ! EIP
419  allocate(species%components(2))
420  species%components = &
421  (/ amns_reactant_type(zn_neut, iz, mi_neut, 0), &
422  amns_reactant_type(zn_neut, iz, mi_neut, 1) &
423  /)
424  CALL itm_amns_setup_table(amns, eip_rx, species, amns_eip(iz, ineut))
425  deallocate(species%components)
426 
427  END IF
428  ENDDO
429 
430  DO iz = 1, zmax
431  IF (zn_neut .GE. iz) THEN
432 ! CX
433  allocate(species%components(4))
434  species%components = &
435  (/ amns_reactant_type(zn_neut, iz, mi_neut, 0), &
436  amns_reactant_type(1, 0, 0, 0), &
437  amns_reactant_type(zn_neut, iz-1, mi_neut, 1), &
438  amns_reactant_type(1, 1, 0, 1) &
439  /)
440  CALL itm_amns_setup_table(amns, cx_rx, species, amns_cx(iz, ineut))
441  deallocate(species%components)
442 ! RC
443  allocate(species%components(4))
444  species%components = &
445  (/ amns_reactant_type(zn_neut, iz, mi_neut, 0), &
446  amns_reactant_type(0, -1, 0, 0), &
447  amns_reactant_type(zn_neut, iz-1, mi_neut, 1), &
448  amns_reactant_type(0, -1, 0, 1) &
449  /)
450  CALL itm_amns_setup_table(amns, rc_rx, species, amns_rc(iz, ineut))
451  deallocate(species%components)
452  END IF
453  ENDDO
454 
455  END DO
456 #endif
457  first=.false.
458  ENDIF
459 
460 
461 
462 
463 
464 
465 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
466 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
467 ! +++ Calculate neutrals:
468  neutral_type_loop1: DO ineut=1,nneut
469 
470 
471 
472 
473 
474 
475 
476 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
477 ! + + + + start calculation for neutrals + + + + + + + +
478 
479  potential = 0.0_r8
480  ionizat = 0.0_r8
481  recomb = 0.0_r8
482  chargeexch = 0.0_r8
483 
484 
485 #ifdef GOT_AMNSPROTO
486  DO iz=0,zmax-1
487  CALL itm_amns_rx(amns_ei(iz,ineut), ionizat(:,iz), te, ne)
488  END DO
489  DO iz=1,zmax
490  CALL itm_amns_rx(amns_eip(iz,ineut), potential(:,iz), te, ne)
491  CALL itm_amns_rx(amns_rc(iz,ineut), recomb(:,iz), te, ne)
492  CALL itm_amns_rx(amns_cx(iz,ineut), chargeexch(:,iz), te, ne)
493  END DO
494 #endif
495 
496 
497 
498 
499 
500 
501 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
502 ! solution of particle transport equation for +
503 ! individual state of neutrals for low energy +
504 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
505  cold_neutrals = 0
506  IF(ASSOCIATED(coreneutrals_old(1)%compositions%neutralscomp(ineut)%type))THEN
507  DO itype = 1, ntype(ineut)
508  IF (coreneutrals_old(1)%compositions%neutralscomp(ineut)%type(itype)%flag .EQ. 0) &
509  cold_neutrals = 1
510  END DO
511  END IF
512 
513  IF (cold_neutrals .EQ. 0) goto 123
514 
515 
516 
517 ! ***************** ITYPE = 1, cold neutrals
518 
519  itype = 1
520 
521 
522  DO irho =1 , nrho
523  nnsource(irho,ineut,itype) = 0.0
524  diff(irho,ineut,itype) = (9.56d7*t0(irho,ineut,itype)/aneut(ineut)) &
525  / (3.0*(ne(irho)*ionizat(irho,0)+chargeexch(irho,1)*ni(irho,ineut)))
526  END DO
527 
528 ! +++ Set equation t0 'predictive' and all coefficients to zero:
529  flag = 1
530  y(:) = 0.0d0
531  dy(:) = 0.0d0
532  ym(:) = 0.0d0
533  a(:) = 0.0d0
534  b(:) = 0.0d0
535  c(:) = 0.0d0
536  d(:) = 0.0d0
537  e(:) = 0.0d0
538  f(:) = 0.0d0
539  g(:) = 0.0d0
540  h = 0.0d0
541  v(:) = 0.0d0
542  u(:) = 0.0d0
543  w(:) = 0.0d0
544 
545 
546 ! +++ Coefficients for ion diffusion equation in form:
547 !
548 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
549 
550  DO irho=1,nrho
551  y(irho) = n0(irho,ineut,itype)
552  ym(irho) = n0m(irho,ineut,itype)
553  a(irho) = vpr(irho)
554  b(irho) = vprm(irho)
555  c(irho) = 1.d0
556  d(irho) = vpr(irho)*g3(irho)*diff(irho,ineut,itype)
557  e(irho) = vpr(irho)*g3(irho)*vconv(irho,ineut,itype)
558  f(irho) = vpr(irho)*nnsource(irho,ineut,itype)
559  g(irho) = vpr(irho)*(ne(irho)*ionizat(irho,0)+ni(irho,ineut)*chargeexch(irho,1))
560  END DO
561  h = tau
562 
563 
564 
565 ! +++ Boundary conditions for ion diffusion equation in form:
566 !
567 ! V*Y' + U*Y =W
568 !
569 ! +++ On axis:
570 ! dNN/drho(rho=0)=0:
571  v(1) = 1.d0
572  u(1) = 0.d0
573  w(1) = 0.d0
574 
575 ! +++ At the edge:
576 ! FIXED NN
577  IF(nn_bnd_type(ineut,itype).EQ.1) THEN
578  v(2) = 0.d0
579  u(2) = 1.d0
580  w(2) = nn_bnd(1,ineut,itype)
581  ENDIF
582 
583 
584 ! FIXED grad_NN
585  IF(nn_bnd_type(ineut,itype).EQ.2) THEN
586  v(2) = 1.d0
587  u(2) = 0.d0
588  w(2) = nn_bnd(1,ineut,itype)
589  ENDIF
590 
591 ! FIXED L_NN
592  IF(nn_bnd_type(ineut,itype).EQ.3) THEN
593  v(2) = nn_bnd(1,ineut,itype)
594  u(2) = 1.d0
595  w(2) = 0.d0
596  ENDIF
597 
598 ! FIXED Flux_NN
599  IF(nn_bnd_type(ineut,itype).EQ.4) THEN
600  v(2) = -g3(nrho)*diff(nrho,ineut,itype)*vpr(nrho)
601  u(2) = g3(nrho)*vconv(nrho,ineut,itype)*vpr(nrho)
602  w(2) = nn_bnd(1,ineut,itype)
603  ENDIF
604 
605 
606 ! Generic boundary condition
607  IF(nn_bnd_type(ineut,itype).EQ.5) THEN
608  v(2) = nn_bnd(1,ineut,itype)
609  u(2) = nn_bnd(2,ineut,itype)
610  w(2) = nn_bnd(3,ineut,itype)
611  ENDIF
612 
613 
614 
615 ! +++ Density equation is not solved:
616  IF (nn_bnd_type(ineut,itype).EQ.0) THEN
617  CALL l3deriv(y, rho, nrho, dy, rho, nrho)
618 
619  flag = 0
620 
621  rho_loop4: DO irho=1,nrho
622  a(irho) = 1.0d0
623  b(irho) = 1.0d0
624  c(irho) = 1.0d0
625  d(irho) = 0.0d0
626  e(irho) = 0.0d0
627  f(irho) = 0.0d0
628  g(irho) = 0.0d0
629 
630  END DO rho_loop4
631 
632  v(2) = 0.0d0
633  u(2) = 1.0d0
634  w(2) = y(nrho)
635  END IF
636 
637 
638 ! +++ Defining coefficients for numerical solver:
639  solver%TYPE = solver_type
640  solver%EQ_FLAG(ndim) = flag
641  solver%NDIM = ndim
642  solver%NRHO = nrho
643  solver%AMIX = amix
644 
645 
646  rho_loop5: DO irho=1,nrho
647 
648  solver%RHO(irho) = rho(irho)
649  solver%Y(ndim,irho) = y(irho)
650  solver%DY(ndim,irho) = dy(irho)
651  solver%YM(ndim,irho) = ym(irho)
652  solver%A(ndim,irho) = a(irho)
653  solver%B(ndim,irho) = b(irho)
654  solver%C(ndim,irho) = c(irho)
655  solver%D(ndim,irho) = d(irho)
656  solver%E(ndim,irho) = e(irho)
657  solver%F(ndim,irho) = f(irho)
658  solver%G(ndim,irho) = g(irho)
659 
660  END DO rho_loop5
661 
662  solver%H = h
663 
664  solver%V(ndim,1) = v(1)
665  solver%U(ndim,1) = u(1)
666  solver%W(ndim,1) = w(1)
667  solver%V(ndim,2) = v(2)
668  solver%U(ndim,2) = u(2)
669  solver%W(ndim,2) = w(2)
670 
671 
672 ! +++ Solution of density diffusion equation:
673  CALL solution_interface(solver, ifail)
674 
675 
676 ! +++ New neutrals density:
677  rho_loop6: DO irho=1,nrho
678  y(irho) = solver%Y(ndim,irho)
679  dy(irho) = solver%DY(ndim,irho)
680 
681  END DO rho_loop6
682 
683 
684 ! +++ New profiles of neutrals density flux and integral source:
685  rho_loop7: DO irho=1,nrho
686 
687  n0m(irho,ineut,itype) = n0(irho,ineut,itype)
688 
689  n0(irho,ineut,itype) = y(irho)
690  dn0(irho,ineut,itype) = dy(irho)
691  fun(irho) = 1.d0/rho(irho)*(vpr(irho)*nnsource(irho,ineut,itype) &
692  + vprm(irho)*n0m(irho,ineut,itype)/tau &
693  - n0(irho,ineut,itype)*vpr(irho)*(1.d0/tau))
694  flux(irho,ineut,itype) = vpr(irho)*g3(irho)* &
695  ( y(irho)*vconv(irho,ineut,itype) - dy(irho)*diff(irho,ineut,itype))
696 
697  END DO rho_loop7
698 
699 
700 
701  123 CONTINUE
702 
703 
704  hot_neutrals = 0
705  IF(ASSOCIATED(coreneutrals_old(1)%compositions%neutralscomp(ineut)%type))THEN
706  DO itype = 1, ntype(ineut)
707  IF (coreneutrals_old(1)%compositions%neutralscomp(ineut)%type(itype)%flag .EQ. 1) &
708  hot_neutrals = 1
709  END DO
710  END IF
711 
712  IF (hot_neutrals .EQ. 0) goto 124
713 
714 !********************* ITYPE = 2, hot neutrals
715  itype=2
716 
717  DO irho=1,nrho
718  nnsource(irho,ineut,itype)=0.0
719  END DO
720 
721  DO irho=1,nrho
722 ! DIFF(IRHO,INEUT,ITYPE)=(9.56D7*T0(IRHO,INEUT,ITYPE)/ANEUT(INEUT))/(3.0*(ne(IRHO)*IONIZAT(IRHO,0)+CHARGEEXCH(IRHO,1)*ni(IRHO,INEUT)))
723  diff(irho,ineut,itype)=(9.56d7*t0(irho,ineut,itype)/aneut(ineut))/(3.0*(ne(irho)*ionizat(irho,0)+chargeexch(irho,1)*ni(irho,ineut)))
724  END DO
725 
726 ! +++ Set equation t0 'predictive' and all coefficients to zero:
727  flag = 1
728  y(:) = 0.0d0
729  dy(:) = 0.0d0
730  ym(:) = 0.0d0
731  a(:) = 0.0d0
732  b(:) = 0.0d0
733  c(:) = 0.0d0
734  d(:) = 0.0d0
735  e(:) = 0.0d0
736  f(:) = 0.0d0
737  g(:) = 0.0d0
738  h = 0.0d0
739  v(:) = 0.0d0
740  u(:) = 0.0d0
741  w(:) = 0.0d0
742 
743 
744 ! +++ Coefficients for ion diffusion equation in form:
745 !
746 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
747 
748  rho_loop8: DO irho=1,nrho
749  y(irho) = n0(irho,ineut,itype)
750  ym(irho) = n0m(irho,ineut,itype)
751  a(irho) = vpr(irho)
752  b(irho) = vprm(irho)
753  c(irho) = 1.d0
754  d(irho) = vpr(irho)*g3(irho)*diff(irho,ineut,itype)
755  e(irho) = vpr(irho)*g3(irho)*vconv(irho,ineut,itype)
756  f(irho) = vpr(irho)*nnsource(irho,ineut,itype)+ &
757  vpr(irho)*ni(irho,ineut)*chargeexch(irho,1)*n0(irho,ineut,1)+ &
758  vpr(irho)*ni(irho,ineut)*ne(irho)*recomb(irho,1)
759  g(irho) = vpr(irho)*(ni(irho,ineut)*chargeexch(irho,1)+ne(irho)*ionizat(irho,0))
760 
761  END DO rho_loop8
762 
763  h = tau
764 
765 
766 ! +++ Boundary conditions for ion diffusion equation in form:
767 !
768 ! V*Y' + U*Y =W
769 !
770 ! +++ On axis:
771 ! dNN/drho(rho=0)=0:
772  v(1) = 1.d0
773  u(1) = 0.d0
774  w(1) = 0.d0
775 
776 ! +++ At the edge:
777 ! FIXED NN
778  IF(nn_bnd_type(ineut,itype).EQ.1) THEN
779  v(2) = 0.d0
780  u(2) = 1.d0
781  w(2) = nn_bnd(1,ineut,itype)
782  ENDIF
783 
784 ! FIXED grad_NN
785  IF(nn_bnd_type(ineut,itype).EQ.2) THEN
786  v(2) = 1.d0
787  u(2) = 0.d0
788  w(2) = nn_bnd(1,ineut,itype)
789  ENDIF
790 
791 ! FIXED L_NN
792  IF(nn_bnd_type(ineut,itype).EQ.3) THEN
793  v(2) = nn_bnd(1,ineut,itype)
794  u(2) = 1.d0
795  w(2) = 0.d0
796  ENDIF
797 
798 ! FIXED Flux_NN
799  IF(nn_bnd_type(ineut,itype).EQ.4) THEN
800  v(2) = -g3(nrho)*diff(nrho,ineut,itype)*vpr(nrho)
801  u(2) = g3(nrho)*vconv(nrho,ineut,itype)*vpr(nrho)
802  w(2) = nn_bnd(1,ineut,itype)
803  ENDIF
804 
805 
806 ! Generic boundary condition
807  IF(nn_bnd_type(ineut,itype).EQ.5) THEN
808  v(2) = nn_bnd(1,ineut,itype)
809  u(2) = nn_bnd(2,ineut,itype)
810  w(2) = nn_bnd(3,ineut,itype)
811  ENDIF
812 
813 
814 
815 ! +++ Density equation is not solved:
816  IF(nn_bnd_type(ineut,itype).EQ.0) THEN
817 
818  CALL l3deriv(y, rho, nrho, dy, rho, nrho)
819 
820  flag = 0
821 
822  rho_loop9: DO irho=1,nrho
823 
824  a(irho) = 1.0d0
825  b(irho) = 1.0d0
826  c(irho) = 1.0d0
827  d(irho) = 0.0d0
828  e(irho) = 0.0d0
829  f(irho) = 0.0d0
830  g(irho) = 0.0d0
831 
832  END DO rho_loop9
833 
834  v(2) = 0.0d0
835  u(2) = 1.0d0
836  w(2) = y(nrho)
837  END IF
838 
839 
840 ! +++ Defining coefficients for numerical solver:
841  solver%TYPE = solver_type
842  solver%EQ_FLAG(ndim) = flag
843  solver%NDIM = ndim
844  solver%NRHO = nrho
845  solver%AMIX = amix
846 
847 
848  rho_loop10: DO irho=1,nrho
849 
850  solver%RHO(irho) = rho(irho)
851  solver%Y(ndim,irho) = y(irho)
852  solver%DY(ndim,irho) = dy(irho)
853  solver%YM(ndim,irho) = ym(irho)
854  solver%A(ndim,irho) = a(irho)
855  solver%B(ndim,irho) = b(irho)
856  solver%C(ndim,irho) = c(irho)
857  solver%D(ndim,irho) = d(irho)
858  solver%E(ndim,irho) = e(irho)
859  solver%F(ndim,irho) = f(irho)
860  solver%G(ndim,irho) = g(irho)
861 
862  END DO rho_loop10
863 
864  solver%H = h
865 
866  solver%V(ndim,1) = v(1)
867  solver%U(ndim,1) = u(1)
868  solver%W(ndim,1) = w(1)
869  solver%V(ndim,2) = v(2)
870  solver%U(ndim,2) = u(2)
871  solver%W(ndim,2) = w(2)
872 
873 
874 ! +++ Solution of density diffusion equation:
875  CALL solution_interface(solver, ifail)
876 
877 
878 ! +++ New impurity density:
879  rho_loop11: DO irho=1,nrho
880 
881  y(irho) = solver%Y(ndim,irho)
882  dy(irho) = solver%DY(ndim,irho)
883 
884  END DO rho_loop11
885 
886 
887 ! +++ New profiles of impurity density flux and integral source:
888  rho_loop12: DO irho=1,nrho
889 
890  n0m(irho,ineut,itype) = n0(irho,ineut,itype)
891  n0(irho,ineut,itype) = y(irho)
892  dn0(irho,ineut,itype) = dy(irho)
893  fun(irho) = 1.d0/rho(irho)*(vpr(irho)*nnsource(irho,ineut,itype) &
894  + vprm(irho)*n0m(irho,ineut,itype)/tau &
895  - n0(irho,ineut,itype)*vpr(irho)*(1.d0/tau))
896  flux(irho,ineut,itype) = vpr(irho)*g3(irho)* &
897  ( y(irho)*vconv(irho,ineut,itype) - dy(irho)*diff(irho,ineut,itype))
898 
899  END DO rho_loop12
900 
901  DO irho=1,nrho
902 
903 
904  DO iion=1,nion
905  IF(abs(aneut(ineut)-coreprof_iter(1)%compositions%nuclei(coreprof_iter(1)%compositions%IONS(iion)%nucindex)%amn).LE.0.25) THEN
906  DO itype=1,ntype(ineut)
907  si_exp(irho,iion) = si_exp(irho,iion)+ne(irho)*(n0(irho,ineut,itype)*ionizat(irho,0)-recomb(irho,1)*ni(irho,ineut))
908  END DO
909  ENDIF
910  END DO
911 
912  IF (ASSOCIATED(coresource_new(1)%VALUES(1)%sz)) THEN
913 
914  DO iimp=1,nimp
915  IF(abs(aneut(ineut)-coreimpur_iter(1)%compositions%nuclei(coreimpur_iter(1)%compositions%IMPURITIES(iimp)%nucindex)%amn).LE.0.25) THEN
916  DO itype=1,ntype(ineut)
917  sz_exp(irho,iimp,1) = sz_exp(irho,iimp,1)+ne(irho)*(n0(irho,ineut,itype)*ionizat(irho,0)-recomb(irho,1)*ni(irho,ineut))
918  END DO
919  ENDIF
920  END DO
921 
922  ENDIF
923 
924  END DO
925 
926 
927 ! +++ Energy source/sink caused by neutrals:
928  DO irho=1,nrho
929 
930 
931  DO iion=1,nion
932  IF(abs(aneut(ineut)-coreprof_iter(1)%compositions%nuclei(coreprof_iter(1)%compositions%IONS(iion)%nucindex)%amn).LE.0.25) THEN
933  DO itype=1,ntype(ineut)
934 
935  qe_exp(irho) = qe_exp(irho) - 1.5_r8 * itm_ev * &
936  (ionizat(irho,0) * potential(irho,1) * n0(irho,ineut,itype) * ne(irho))
937 ! losses by recombination are included in impurity part
938  qi_exp(irho,iion) = qi_exp(irho,iion) + 1.5_r8 * itm_ev * &
939  (chargeexch(irho,1) * (t0(irho,ineut,itype)-ti(irho,iion)) * n0(irho,ineut,itype) * ni(irho,ineut))
940 
941  END DO
942  ENDIF
943 
944  qe_exp_tot(irho) = qe_exp_tot(irho) + qe_exp(irho)
945 
946 
947  END DO
948 
949  IF (ASSOCIATED(coresource_new(1)%VALUES(1)%sz)) THEN
950 
951  DO iimp=1,nimp
952  IF(abs(aneut(ineut)-coreimpur_iter(1)%compositions%nuclei(coreimpur_iter(1)%compositions%IMPURITIES(iimp)%nucindex)%amn).LE.0.25) THEN
953  DO itype=1,ntype(ineut)
954 
955  qe_exp(irho) = qe_exp(irho) - 1.5_r8 * itm_ev * &
956  ( ionizat(irho,0) * potential(irho,1) * n0(irho,ineut,itype) * ne(irho))
957 ! losses by recombination are included in impurity part
958 
959  qz_exp(irho,iimp,1)= qz_exp(irho,iimp,1) + 1.5_r8 * itm_ev * &
960  (chargeexch(irho,1) * (t0(irho,ineut,itype)-ti(irho,1)) * n0(irho,ineut,itype) * ni(irho,ineut))
961 !irena Split
962 ! QZ_EXP(IRHO,IIMP,1)= QZ_EXP(IRHO,IIMP,1) + 1.5_R8 * ITM_EV * &
963 ! ( CHARGEEXCH(IRHO,1) * (T0(IRHO,INEUT,ITYPE)-TI(IRHO)) * N0(IRHO,INEUT,ITYPE) * NI(IRHO,INEUT))
964 
965  END DO
966  ENDIF
967 
968  qe_exp_tot(irho) = qe_exp_tot(irho) + qe_exp(irho)
969 
970  END DO
971  END IF
972 
973 
974 
975  124 CONTINUE
976 
977 
978  END DO
979 
980 ! Split irena calculate the radiation only ionization losses
981 
982 
983 
984 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
985 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
986 
987 
988 ! +++ Save output to CPOs:
989 
990  loop_irho20: DO irho=1,nrho
991  DO itype=1,ntype(ineut)
992  coreneutrals_new(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%n0%value(irho) = n0(irho,ineut,itype)
993  coreneutrals_new(1)%PROFILES(ineut)%NEUTRALTYPE(itype)%n0%flux(irho) = flux(irho,ineut,itype)
994  END DO
995 
996 
997  DO iion=1,nion
998  IF(abs(aneut(ineut)-coreprof_iter(1)%compositions%nuclei(coreprof_iter(1)%compositions%IONS(iion)%nucindex)%amn).LE.0.25) THEN
999  coresource_new(1)%VALUES(1)%si%exp(irho,iion) = si_exp(irho,iion)
1000  coresource_new(1)%VALUES(1)%qi%exp(irho,iion) = qi_exp(irho,iion)
1001 ! irena split
1002  coresource_new(1)%VALUES(1)%qe%exp(irho) = qe_exp(irho)
1003 ! irena split
1004  coresource_new(1)%VALUES(1)%se%exp(irho) = coresource_new(1)%VALUES(1)%se%exp(irho) + &
1005  coresource_new(1)%VALUES(1)%si%exp(irho,iion)
1006  END IF
1007  END DO
1008 
1009 
1010  IF (ASSOCIATED(coresource_new(1)%VALUES(1)%sz)) THEN
1011  DO iimp=1,nimp
1012  IF(abs(aneut(ineut)-coreimpur_iter(1)%compositions%nuclei(coreimpur_iter(1)%compositions%IMPURITIES(iimp)%nucindex)%amn).LE.0.25) THEN
1013  coresource_new(1)%VALUES(1)%sz(iimp)%exp(irho,1) = sz_exp(irho,iimp,1)
1014  coresource_new(1)%VALUES(1)%se%exp(irho) = coresource_new(1)%VALUES(1)%se%exp(irho) + &
1015  coresource_new(1)%VALUES(1)%sz(iimp)%exp(irho,1)
1016  END IF
1017  END DO
1018  END IF
1019 
1020 
1021  coresource_new(1)%VALUES(1)%qe%exp(irho) = qe_exp_tot(irho)
1022 
1023 
1024 
1025  END DO loop_irho20
1026 
1027 
1028 
1029 
1030  END DO neutral_type_loop1
1031 
1032 ! Irena Split jonization losses by neutrals
1033 
1034  DO ineut=1, nneut
1035  coreneutrals_new(1)%PROFILES(ineut)%prad0 = 0.d0
1036  END DO
1037  DO irho=1,nrho
1038  DO ineut=1, nneut
1039  DO itype=1,ntype(ineut)
1040  coreneutrals_new(1)%PROFILES(ineut)%prad0(irho) = coreneutrals_new(1)%PROFILES(ineut)%prad0(irho) + &
1041  itm_ev * ionizat(irho,0) * potential(irho,1) * n0(irho,ineut,itype) * ne(irho)
1042  END DO
1043  END DO
1044  END DO
1045 
1046 ! Finished loop over neutral types
1047 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1048 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1049 
1050 
1051 
1052 
1053  ALLOCATE (coresource_new(1)%VALUES(1)%sourceid%id(1))
1054  ALLOCATE (coresource_new(1)%VALUES(1)%sourceid%description(1))
1055  coresource_new(1)%VALUES(1)%sourceid%id = 'gaspuff'
1056  coresource_new(1)%VALUES(1)%sourceid%flag = 21
1057  coresource_new(1)%VALUES(1)%sourceid%description = 'Gas puff'
1058 
1059  DEALLOCATE ( rho )
1060  DEALLOCATE ( vpr )
1061  DEALLOCATE ( vprm )
1062  DEALLOCATE ( g3 )
1063  DEALLOCATE ( ne )
1064  DEALLOCATE ( te )
1065  DEALLOCATE ( ti )
1066  DEALLOCATE ( ni )
1067 
1068  DEALLOCATE ( qe_exp )
1069  DEALLOCATE ( qe_exp_tot )
1070  DEALLOCATE ( si_exp )
1071  DEALLOCATE ( qi_exp )
1072  DEALLOCATE ( sz_exp )
1073  DEALLOCATE ( qz_exp )
1074  DEALLOCATE (ssi_exp )
1075 
1076  DEALLOCATE ( y )
1077  DEALLOCATE ( ym )
1078  DEALLOCATE ( dy )
1079  DEALLOCATE ( a )
1080  DEALLOCATE ( b )
1081  DEALLOCATE ( c )
1082  DEALLOCATE ( d )
1083  DEALLOCATE ( e )
1084  DEALLOCATE ( f )
1085  DEALLOCATE ( g )
1086 
1087  DEALLOCATE ( fun )
1088  DEALLOCATE (intfun )
1089  DEALLOCATE ( aneut )
1090 
1091  DEALLOCATE ( n0 )
1092  DEALLOCATE ( n0m )
1093  DEALLOCATE ( vconv )
1094  DEALLOCATE ( diff )
1095  DEALLOCATE ( dn0 )
1096  DEALLOCATE ( flux )
1097  DEALLOCATE ( nn_bnd )
1098  DEALLOCATE (nn_bnd_type )
1099  DEALLOCATE ( t0 )
1100  DEALLOCATE ( vt )
1101  DEALLOCATE ( vp )
1102  DEALLOCATE ( vr )
1103  DEALLOCATE ( nnsource )
1104 
1105  DEALLOCATE ( ionizat )
1106  DEALLOCATE ( recomb )
1107  DEALLOCATE ( potential )
1108  DEALLOCATE ( chargeexch )
1109 
1110  WRITE (*,*) 'NEUTRALS finished <==========='
1111  WRITE (*,*) ' '
1112 
1113 
1114  RETURN
1115 
1116 
1117  END SUBROUTINE neutrals_ets
1118 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1119 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1120 
1121 
1122 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1123 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1124  SUBROUTINE writeoutneutrals (ITIME_OUT,CORENEUTRALS)
1125 
1126 
1127 ! This subroutine stores the results of computations
1128 ! into files
1129 
1130  USE euitm_schemas
1131  USE itm_types
1132  IMPLICIT NONE
1133 
1134  TYPE (type_coreneutrals), POINTER :: coreneutrals(:) !input CPO with neutrals
1135 
1136 ! +++ Internal parameters:
1137  INTEGER :: nrho !number of radial points (input, determined from COREPROF CPO)
1138  INTEGER :: nneut !number of neutrals species (input, determined from CORENEUTRALSCPO)
1139  integer :: itime_out
1140  integer :: ineut !index of impurity component, number of considered impurity components (max ionization state)
1141  integer :: irho
1142  INTEGER :: ntype,itype !type part for one neutrals
1143  REAL (R8), ALLOCATABLE :: rho(:) !toroidal flux coordinate,not normalise,[m]
1144 
1145 
1146  CHARACTER (33) filename
1147 
1148 
1149  nrho = SIZE (coreneutrals(1)%RHO_TOR, dim=1)
1150  nneut = SIZE (coreneutrals(1)%profiles,dim=1)
1151 
1152  ALLOCATE (rho(nrho))
1153 
1154  rho = coreneutrals(1)%RHO_TOR
1155 
1156  IF (nneut.GE.2) THEN
1157  WRITE(*,*)'in programe Write neutrals'
1158  WRITE(*,*)'CORENEUTRALS(1)%profiles(1)%neutraltype(1)%n0%value(NRHO)=',coreneutrals(1)%profiles(1)%neutraltype(1)%n0%value(nrho)
1159  ENDIF
1160 
1161  DO ineut=1,nneut
1162 
1163  ntype = SIZE(coreneutrals(1)%profiles(ineut)%neutraltype)
1164 
1165  WRITE(*,*)'ntype=',ntype,'irena'
1166 
1167  WRITE(filename,'(a,i1.1,a,i7.7,a)') 'eq_ets_data/OUTNE',ineut,'/NEU',itime_out,'.DAT'
1168  OPEN (unit=20, file=filename)
1169 
1170  DO irho = 1, nrho
1171  WRITE (20,'(100(1x,e14.7))')rho(irho),(coreneutrals(1)%profiles(ineut)%neutraltype(itype)%n0%value(irho),itype=1,ntype),&
1172  coreneutrals(1)%profiles(ineut)%prad0(irho)
1173  END DO
1174 
1175  CLOSE (20)
1176 
1177  END DO
1178 
1179  DEALLOCATE (rho)
1180 
1181 
1182  RETURN
1183  END SUBROUTINE writeoutneutrals
1184 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1185 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1186 
1187 
1188  END MODULE neutrals
l3deriv
subroutine l3deriv(y_in, x_in, nr_in, dydx_out, x_out, nr_out)
Definition: l3interp.f90:59
solution_interface
subroutine solution_interface(SOLVER, ifail)
INTERFACE TO NUMERICAL SOLVER.
Definition: solution_interface.f90:12
solver
Definition: solver.f90:1
neutrals::writeoutneutrals
subroutine writeoutneutrals(ITIME_OUT, CORENEUTRALS)
Definition: neutrals.F90:1124
type_solver
Definition: type_solver.f90:1
fun
subroutine fun(X, F)
Definition: Ev2.f:10
allocate_deallocate::get_comp_dimensions
subroutine get_comp_dimensions(COMPOSITIONS, NNUCL, NION, NIMP, NZIMP, NNEUT, NTYPE, NCOMP)
Definition: allocate_deallocate.f90:1845
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
l3interp
subroutine l3interp(y_in, x_in, nr_in, y_out, x_out, nr_out)
Definition: l3interp.f90:1
allocate_deallocate
This module contains routines for allocation/deallocation if CPOs used in ETS.
Definition: allocate_deallocate.f90:8
allocate_deallocate::allocate_coresource_cpo
subroutine allocate_coresource_cpo(NSLICE, NRHO, NNUCL, NION, NIMP, NZIMP, NNEUT, NTYPE, NCOMP, CORESOURCE)
This routine allocates CORESOURCE CPO.
Definition: allocate_deallocate.f90:460
neutrals::neutrals_ets
subroutine neutrals_ets(COREIMPUR_ITER, EQUILIBRIUM_ITER, COREPROF_ITER, CORENEUTRALS_OLD, CORENEUTRALS_ITER, CORESOURCE_NEW, CORENEUTRALS_NEW, CONTROL_INTEGER, CONTROL_DOUBLE)
Definition: neutrals.F90:12
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
allocate_deallocate::allocate_coreneutrals_cpo
subroutine allocate_coreneutrals_cpo(NSLICE, NRHO, NNUCL, NION, NIMP, NZIMP, NNEUT, NTYPE, NCOMP, CORENEUTRALS)
Definition: allocate_deallocate.f90:971