ETS-Core  version:0.0.4-46-ge2d8
Core actors for the ETS-6
 All Classes Files Functions Variables Pages
impurity.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
7 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
8 MODULE impurity
9 
10 CONTAINS
11 
12 
13 ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
14 ! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15  SUBROUTINE impurity_imas (EQUILIBRIUM_ITER, CORETRANSPORT_ITER, CORETRANSPORT_NEW , &
16  coreprofiles_old, coreprofiles_iter, coreprofiles_new , &
17  coresources_iter, coresources_new, coresolver_iter, &
18  coreradiation_iter, coreradiation_new, &
19  control_integer, control_double )
20 
21 
22 !-------------------------------------------------------!
23 ! This routine calculate impurity density for !
24 ! different impurity !
25 !-------------------------------------------------------!
26 ! Source: --- !
27 ! Developers: I.M.Ivanova-Stanik !
28 ! Kontacts: irena.ivanova-stanik@ifpilm.pl !
29 ! !
30 ! Comments: for IMAS 20 !
31 ! !
32 ! !
33 !-------------------------------------------------------!
34 
35  use ids_schemas
36  use ids_types, r8 => ids_real, itm_i4 => ids_int
37  use ids_routines
38  use transport_solver
39  use l3interps
40  use imas_constants_module !, only: imas_constants => constants_module
41  use ets_math
42  use get_composition
43  use control_parameters, only: get_control_parameter
44 
45  USE allocate_working_ids
46  USE ets_plasma
47 
48 
49 ! +++++++++ atomic data ++++++++++
50  USE amns_types
51  USE amns_module
52 
53  implicit none
54 
55  TYPE (ids_equilibrium) :: equilibrium_iter
56  TYPE (ids_core_transport) :: coretransport_iter
57  TYPE (ids_core_transport) :: coretransport_new
58  TYPE (ids_core_profiles) :: coreprofiles_old
59  TYPE (ids_core_profiles) :: coreprofiles_iter
60  TYPE (ids_core_profiles) :: coreprofiles_new
61  TYPE (ids_core_sources) :: coresources_iter
62  TYPE (ids_core_sources) :: coresources_new
63  TYPE (ids_transport_solver_numerics) :: coresolver_iter
64  TYPE (ids_radiation) :: coreradiation_iter
65  TYPE (ids_radiation) :: coreradiation_new
66 !
67  TYPE (plasma_composition) :: composition
68 
69 ! +++ Dimensions:
70  INTEGER :: neq !number of radial points (input, determined from EQUILIBRIUM IMAS)
71  INTEGER :: nrho !number of radial points (input, determined from COREPROF IMAS)
72  INTEGER :: nnucl !number of ion species (input, determined from COREPROF IMAS)
73  INTEGER :: nequal !number of equations (input, determined from CORESLVER IMAS)
74  INTEGER :: nsource !number of the source (input, determined from CORESLVER IMAS)
75  INTEGER :: nproces !number of the proces (input, determined from RADIATION_IMAS)
76  INTEGER :: nneut !number of neutrals species (input, determined from COREPROF IMAS)
77 
78  INTEGER :: nion !number of ion species in COREPROFILES
79  INTEGER :: nion_ct !number of ion species in CORETRANSPORT
80  INTEGER :: nion_cs !number of ion species in CORESOURCES
81  INTEGER :: nion_csol !number of ion species in CORESOLVER
82  INTEGER :: nion_lin !number of ion species in RADIATION:Line
83  INTEGER :: nion_brem !number of ion species in RADIATION:Bremst
84  INTEGER :: nion_rec !number of ion species in RADIATION:Recombination
85 
86 
87  INTEGER :: nimp !number of impurity species in COREPROFILES
88  INTEGER :: nimp_ct !number of impurity species in CORETRANSPORT
89  INTEGER :: nimp_csol !number of impurity species in CORESOLVER
90  INTEGER :: nimp_cs !number of impurity species in CORESOURCES
91  INTEGER :: nimp_lin,nimp_brem,nimp_rec !number of impurity species in LIN CORERADIATION
92 
93  INTEGER, ALLOCATABLE :: nzimp(:) !number of ionization state for (:) impurity in COREPROFILES
94  INTEGER, ALLOCATABLE :: nzimp_ct(:) !number of ionization state for (:) impurity in CORETRANSPORT
95  INTEGER, ALLOCATABLE :: nzimp_cs(:) !number of ionization state for (:) impurity in CORESOURCES
96  INTEGER, ALLOCATABLE :: nzimp_csol(:) !number of ionization state for (:) impurity in CORESOLVER
97  INTEGER, ALLOCATABLE :: nzimp_lin(:),nzimp_brem(:),nzimp_rec(:) !number of ionization state for (:) impurity in CORERADIATION
98  INTEGER, ALLOCATABLE :: nz_bnd_type(:) !boundary condition, type, one impurity
99 
100  INTEGER, ALLOCATABLE :: ntype(:) !number of impurity ionization states (input)
101  INTEGER, ALLOCATABLE :: ncomp(:) !max_number of distinct atoms enter the composition-"1" wich is neutral
102  INTEGER, PARAMETER :: nocur = 1 !number of CPO ocurancies in the work flow
103  INTEGER, ALLOCATABLE :: nrho_tr
104  INTEGER, ALLOCATABLE :: nrho_sr(:)
105  INTEGER, ALLOCATABLE :: nrho_lin(:),nrho_brem(:),nrho_rec(:)
106  INTEGER :: nmod !number on the transport model
107  INTEGER :: inmod !index on the transport model
108  INTEGER :: inmod_imp !index on the transport model used for impurity equation
109  INTEGER :: in_tr !index on the transport model used for impurity equation
110 
111  INTEGER :: insource !index on the different source
112  INTEGER :: insource_imp !index on the transport model
113 
114  INTEGER :: inequal !index on the solver numerical transport
115 ! CORERADIATION
116 
117  INTEGER :: inproces !index on the different proces in radiation
118  INTEGER :: inproces_lin !index on the different proces in radiation
119  INTEGER :: inproces_brem !index on the different proces in radiation
120  INTEGER :: inproces_rec !index on the different proces in radiation
121 
122 ! +++ Indexes:
123  INTEGER :: irho, iimp, izimp,ineq !index of impurity component, number of considered impurity components (max ionization state)
124  INTEGER :: ineut, itype, icomp !number of neutrals species (input)
125  INTEGER :: max_nzimp
126  INTEGER :: nzimp2
127  INTEGER :: iion,iion_ct,iion_csol,iion_cs
128  INTEGER :: iion_lin,iion_brem,iion_rec
129  INTEGER :: solver_type !specifies the option for numerical solution
130 
131  REAL (R8) :: b0, b0prime !magnetic field from current time step, [T], previous time steps, [T], time derivative, [T/s]
132  REAL (R8) :: r0
133  REAL (R8), ALLOCATABLE :: rho(:) !toroidal flux coordinate,not normalise,[m] [m]
134  REAL (R8), ALLOCATABLE :: vol(:) !V, [m^3]
135  REAL (R8), ALLOCATABLE :: vpr(:) !V', [m^2]
136  REAL (R8), ALLOCATABLE :: vprm(:) !V' (at previous time step), [m^2]
137  REAL (R8), ALLOCATABLE :: g3(:)
138  REAL (R8), ALLOCATABLE :: ne(:) !electron density [m^-3]
139  REAL (R8), ALLOCATABLE :: te(:) !electron temperature
140  REAL (R8), ALLOCATABLE :: dnz1(:,:) !density gradient, [m^-4]
141  REAL (R8), ALLOCATABLE :: flux(:,:) !ion flux, [1/s]
142  REAL (R8), ALLOCATABLE :: flux_inter(:,:) !ion flux, [1/s]
143  REAL (R8), ALLOCATABLE :: nz1(:,:),aa(:,:) !one impurity density
144  REAL (R8), ALLOCATABLE :: nzm1(:,:) !old one impurity density
145  REAL (R8), ALLOCATABLE :: zn(:) !charge state on the ion
146  REAL (R8), ALLOCATABLE :: am(:) !mas on the ion
147 
148 
149 
150  REAL (R8), ALLOCATABLE :: total_rad_int(:,:) !total radiation for one ionization state
151 ! 26.10.1017 IS
152  REAL (R8), ALLOCATABLE :: tti(:) !ion temperature
153 
154  REAL (R8), ALLOCATABLE :: diff(:,:) !diffusion coefficient for different ionisation, [m^2/s]
155  REAL (R8), ALLOCATABLE :: vcon(:,:) !pinch velocity for different ionisation [m/s]
156  REAL (R8), ALLOCATABLE :: imp_radiation(:,:) !impurity radiation
157  REAL (R8), ALLOCATABLE :: nzsource(:,:) !value of the source term,[m^-3.s^-1]
158  REAL (R8), ALLOCATABLE :: nz_bnd(:,:) !boundary condition, value, one impurity[depends on NZ_BND_TYPE]
159  REAL (R8), ALLOCATABLE :: aneut(:)
160 
161 ! part for radiation
162  REAL (R8), ALLOCATABLE :: lin_rad1(:,:) !profile of lineradiation for one impurity
163  REAL (R8), ALLOCATABLE :: brem_rad1(:,:) !profile of bremst. for one impurity
164  REAL (R8), ALLOCATABLE :: lin_rad(:) !profile of lineradiation for whole impurity
165  REAL (R8), ALLOCATABLE :: brem_rad(:) !profile of bremst. for whole impurity
166  REAL (R8), ALLOCATABLE :: jon_en1(:,:) !profile of jonisation energy for one impurity
167  REAL (R8), ALLOCATABLE :: jon_en(:) !profile of jonisation energy for wholeimpurity
168  REAL (R8), ALLOCATABLE :: rec_los1(:,:) !profile of recombination losses for one impurity
169  REAL (R8), ALLOCATABLE :: rec_los(:) !profile of recombination losses for wholeimpurity
170  REAL (R8), ALLOCATABLE :: qrad(:)
171  REAL (R8), ALLOCATABLE :: se_exp(:)
172 
173  REAL (R8), ALLOCATABLE :: fun(:)
174  REAL (R8), ALLOCATABLE :: fun_in(:)
175  REAL (R8), ALLOCATABLE :: fun_out(:)
176  REAL (R8), ALLOCATABLE :: fun_out_tr(:)
177 
178  REAL (R8) :: amix, tau !mixing factor, time step, [s]
179 
180 
181  REAL (R8) :: time
182 ! LOGICAL, SAVE :: first = .TRUE.
183 
184  REAL (R8), ALLOCATABLE ::r_lin_int(:)
185  REAL (R8), ALLOCATABLE ::r_lin_int1(:,:)
186  REAL (R8), ALLOCATABLE ::r_brem_int(:)
187  REAL (R8), ALLOCATABLE ::r_brem_int1(:,:)
188  REAL (R8), ALLOCATABLE ::r_sum_int(:)
189  REAL (R8), ALLOCATABLE ::r_sum_int1(:,:)
190  REAL (R8), ALLOCATABLE ::e_jon_int(:)
191  REAL (R8), ALLOCATABLE ::e_jon_int1(:,:)
192  REAL (R8), ALLOCATABLE ::e_rec_int(:)
193  REAL (R8), ALLOCATABLE ::e_rec_int1(:,:)
194  REAL (R8), ALLOCATABLE ::e_sum_int(:)
195  REAL (R8), ALLOCATABLE ::e_sum_int1(:,:)
196  REAL (R8), ALLOCATABLE ::sum_lin_int(:)
197  REAL (R8), ALLOCATABLE ::sum_brem_int(:)
198  REAL (R8), ALLOCATABLE ::sum_rad_int(:)
199  REAL (R8), ALLOCATABLE ::sum_jon_int(:)
200  REAL (R8), ALLOCATABLE ::sum_rec_int(:)
201  REAL (R8), ALLOCATABLE ::sum_los_int(:)
202 
203  INTEGER, INTENT(IN) :: control_integer(2) !integer control parameters
204  REAL (R8), INTENT(IN) :: control_double(5) !real control parameters
205  CHARACTER (33) :: filename
206 
207 
208 ! +++ Atomic data:
209 
210  TYPE (amns_handle_type), SAVE :: amns
211  TYPE (amns_handle_rx_type), ALLOCATABLE, SAVE :: &
212  amns_ei(:,:), amns_eip(:,:),amns_rc(:,:), amns_lr(:,:), amns_br(:,:)
213  TYPE (amns_version_type) :: amns_database
214  TYPE (amns_reaction_type) :: ei_rx, eip_rx, rc_rx, lr_rx, br_rx
215  TYPE (amns_reactants_type) :: species
216  TYPE (amns_query_type) :: query
217  TYPE (amns_answer_type) :: answer
218  TYPE (amns_set_type) :: set
219 ! DPC added for testing
220  type (amns_version_type) :: version
221  REAL (R8) :: zn_imp, am_imp
222 
223  CHARACTER (len=80) :: format
224 
225 
226  WRITE (*,*) ' '
227  WRITE (*,*) '===========> IMPURITY started'
228 
229 
230 ! +++ Set dimensions:
231 
232  nrho = SIZE (coreprofiles_iter%profiles_1d(1)%grid%rho_tor_norm)
233  neq = SIZE (equilibrium_iter%time_slice(1)%profiles_1d%rho_tor_norm)
234  IF (ASSOCIATED(coreprofiles_iter%profiles_1d(1)%ion)) THEN
235  nion = SIZE (coreprofiles_iter%profiles_1d(1)%ion)
236  END IF
237 ! +++ Allocate local variables:
238 
239  ALLOCATE (nzimp(nion))
240  ALLOCATE (zn(nion))
241  ALLOCATE (am(nion))
242 
243 ! IF (ASSOCIATED(CORETRANSPORT_ITER(1)%MODEL)) THEN
244  nmod = SIZE (coretransport_iter%MODEL)
245 ! END IF
246 
247 ! IF (ASSOCIATED(CORESOURCES_ITER))
248  nsource = SIZE (coresources_iter%SOURCE)
249 ! END IF
250 
251 ! IF (ASSOCIATED(CORESOLVER_ITER))
252  nequal = SIZE (coresolver_iter%solver_1d(1)%equation)
253 ! END IF
254 
255 ! IF (ASSOCIATED(CORERADIATION_ITER)) THEN
256  nproces = SIZE (coreradiation_iter%PROCESS)
257 ! ELSE
258 ! NPROCES = 0
259 ! END IF
260 
261  write(*,*) 'NRHO=',nrho,'NEQ=',neq,'NION=',nion,'NSOURCE=', nsource, 'NMOD=',nmod, 'NEQUAL=',nequal
262  write(*,*) 'NPROCES=',nproces
263 
264 
265 ! +++++ for output data ++++++++++++++++++
266  ALLOCATE (coreprofiles_new%profiles_1d(1))
267  CALL ids_copy(coreprofiles_iter, coreprofiles_new)
268 
269 
270  ALLOCATE (coresources_new%SOURCE(nsource))
271  CALL ids_copy(coresources_iter, coresources_new)
272 
273 
274  ALLOCATE (coretransport_new%MODEL(nmod))
275  CALL ids_copy(coretransport_iter,coretransport_new)
276 
277  CALL get_prof_composition(coreprofiles_iter,composition)
278 
279  IF (nproces.GT.0.) THEN
280 ! ALLOCATE (CORERADIATION_NEW%PROCESS(NPROCES))
281  CALL ids_copy(coreradiation_iter, coreradiation_new)
282 
283  ELSE
284  nproces=3
285  CALL allocate_coreradiation_ids(1, nproces, nrho, composition, coreradiation_new)
286 
287  coreradiation_new%process(1)%identifier%name="line_radiation"
288  coreradiation_new%process(1)%identifier%index=10
289  coreradiation_new%process(1)%identifier%description="Emission from line radiation"
290 
291  coreradiation_new%process(2)%identifier%name="bremsstrahlung"
292  coreradiation_new%process(2)%identifier%index=8
293  coreradiation_new%process(2)%identifier%description="Emission from bremsstrahlung+recombination"
294 
295  coreradiation_new%process(3)%identifier%name="recombination"
296  coreradiation_new%process(3)%identifier%index=11
297  coreradiation_new%process(3)%identifier%description="Emission from bremsstrahlung+recombination"
298 
299  DO inproces = 1,nproces
300  coreradiation_new%process(inproces)%profiles_1d(1)%emissivity_ion_total(:)=0.0_r8
301  coreradiation_new%process(inproces)%profiles_1d(1)%power_inside_ion_total(:)=0.0_r8
302 
303  DO iion=1, nion
304  nzimp(iion)= SIZE(coreprofiles_iter%profiles_1d(1)%ion(iion)%state)
305  coreradiation_new%process(inproces)%profiles_1d(1)%ion(iion)%multiple_states_flag = coreprofiles_iter%profiles_1d(1)%ion(iion)%multiple_states_flag
306  coreradiation_new%process(inproces)%profiles_1d(1)%grid%rho_tor_norm = coreprofiles_iter%profiles_1d(1)%grid%rho_tor_norm
307  coreradiation_new%process(inproces)%profiles_1d(1)%ion(iion)%element(1) = coreprofiles_iter%profiles_1d(1)%ion(iion)%element(1)
308  coreradiation_new%process(inproces)%profiles_1d(1)%ion(iion)%emissivity(:)=0.0_r8
309  coreradiation_new%process(inproces)%profiles_1d(1)%ion(iion)%power_inside(:)=0.0_r8
310  IF (coreradiation_new%process(inproces)%profiles_1d(1)%ion(iion)%multiple_states_flag.eq.1) then
311 
312  DO izimp=1, nzimp(iion)
313 
314  coreradiation_new%process(inproces)%profiles_1d(1)%ion(iion)%state(izimp)%Z_MIN = izimp
315  coreradiation_new%process(inproces)%profiles_1d(1)%ion(iion)%state(izimp)%Z_MAX = izimp
316  coreradiation_new%process(inproces)%profiles_1d(1)%ion(iion)%state(izimp)%emissivity(:)=0.0_r8
317  coreradiation_new%process(inproces)%profiles_1d(1)%ion(iion)%state(izimp)%power_inside(:)=0.0_r8
318  END DO
319  END IF
320  END DO
321  ENDDO
322  END IF
323 
324 
325  coretransport_new%time(1) = coreprofiles_iter%time(1)
326  coreradiation_new%time(1) = coreprofiles_iter%time(1)
327  coreprofiles_new%time(1) = coreprofiles_iter%time(1)
328  coreradiation_new%time(1) = coreprofiles_iter%time(1)
329 
330 
331 
332 ! ++++++++++++++++++++ radiation +++++++++++++++++++++++++++
333  DO inproces=1,nproces
334  IF (trim(adjustl(coreradiation_new%process(inproces)%identifier%name(1))).eq."line_radiation") THEN
335  inproces_lin = inproces
336  write (*,*) 'INPROCES_LIN=',inproces_lin
337  END IF
338 
339  IF (trim(adjustl(coreradiation_new%process(inproces)%identifier%name(1))).eq."bremsstrahlung") THEN
340  inproces_brem = inproces
341  write (*,*) 'INPROCES_BREM=',inproces_brem
342  END IF
343 
344  IF (trim(adjustl(coreradiation_new%process(inproces)%identifier%name(1))).eq."recombination") THEN
345  inproces_rec = inproces
346  write (*,*) 'INPROCES_REC=',inproces_rec
347  END IF
348  ENDDO
349  ALLOCATE (nrho_lin(inproces_lin))
350  ALLOCATE (nrho_brem(inproces_brem))
351  ALLOCATE (nrho_rec(inproces_rec))
352 
353  nrho_lin(inproces_lin) = SIZE(coreradiation_new%process(inproces_lin)%PROFILES_1D(1)%GRID%rho_tor_norm)
354  write(*,*)'NRHO_LIN(INPROCES_LIN)=',nrho_lin(inproces_lin),'INPROCES_LIN=',inproces_lin
355  nion_lin = SIZE(coreradiation_new%process(inproces_lin)%PROFILES_1D(1)%ion)
356  write(*,*)'nion_lin=',nion_lin
357  IF (nion_lin.NE.nion) THEN
358  write(*,*)'warning different number on the ions between coreprofiles and line radiation'
359  END IF
360 
361  nrho_brem(inproces_brem) = SIZE(coreradiation_new%process(inproces_brem)%PROFILES_1D(1)%GRID%rho_tor_norm)
362  write(*,*)'NRHO_BREM(INPROCES_BREM)=',nrho_brem(inproces_brem),'INPROCES_BREM=',inproces_brem
363  nion_brem = SIZE(coreradiation_new%process(inproces_brem)%PROFILES_1D(1)%ion)
364  IF (nion_brem.NE.nion) THEN
365  write(*,*)'warning different number on the ions between coreprofiles and bremsstrahlung'
366  END IF
367 
368  nrho_rec(inproces_rec) = SIZE(coreradiation_new%process(inproces_rec)%PROFILES_1D(1)%GRID%rho_tor_norm)
369  write(*,*)'NRHO_REC(INPROCES_REC)=',nrho_rec(inproces_rec),'INPROCES_REC=',inproces_rec
370  nion_rec = SIZE(coreradiation_new%process(inproces_rec)%PROFILES_1D(1)%ion)
371  IF (nion_rec.NE.nion) THEN
372  write(*,*)'warning different number on the ions between coreprofiles and recombination'
373  END IF
374 
375  nimp = 0
376  DO iion=1,nion
377  IF (coreprofiles_iter%profiles_1d(1)%ion(iion)%multiple_states_flag.EQ.1) THEN
378  nimp= nimp+1
379  END IF
380  END DO
381 
382  nimp_lin = 0
383  DO iion_lin=1,nion_lin
384  IF (coreradiation_new%process(inproces_lin)%PROFILES_1D(1)%ion(iion_lin)%multiple_states_flag.EQ.1) THEN
385  nimp_lin= nimp_lin + 1
386  END IF
387  END DO
388  write(*,*)'nimp_rad_lin=',nimp_lin
389 
390  IF (nimp_lin.NE.nimp) THEN
391  write(*,*)'warning different number on the impurity between coreprofile and line radiation'
392  END IF
393 
394  ALLOCATE (nzimp_lin(nion_lin))
395 
396  nimp_brem = 0
397  DO iion_brem=1,nion_brem
398  IF (coreradiation_new%process(inproces_brem)%PROFILES_1D(1)%ion(iion_brem)%multiple_states_flag.EQ.1) THEN
399  nimp_brem= nimp_brem+1
400  END IF
401  END DO
402  write(*,*)'nimp_rad_brem=',nimp_brem
403 
404  IF (nimp_brem.NE.nimp) THEN
405  write(*,*)'warning different number on the impurity between coreprofile and brem. radiation'
406  END IF
407 
408  ALLOCATE (nzimp_brem(nion_brem))
409 
410  nimp_rec = 0
411  DO iion_rec = 1,nion_rec
412  IF (coreradiation_new%process(inproces_rec)%PROFILES_1D(1)%ion(iion_rec)%multiple_states_flag.EQ.1) THEN
413  nimp_rec= nimp_rec+1
414  END IF
415  END DO
416  write(*,*)'nimp_rad_rec=',nimp_rec
417  IF (nimp_rec.NE.nimp) THEN
418  write(*,*)'warning different number on the impurity between coreprofile and recombination'
419  END IF
420 
421  ALLOCATE (nzimp_rec(nion_rec))
422 
423  ! Thomas17/9: AMIX = CORESOLVER_ITER%solver_1d(1)%control_parameters%real0d(2)%value
424  call get_control_parameter(coresolver_iter%solver_1d(1)%control_parameters, 'mixing_ratio_kinetic_profiles', amix)
425 
426  ! Thomas17/9: TAU = CORESOLVER_ITER%solver_1d(1)%control_parameters%real0d(1)%value
427  tau = coresolver_iter%time_step_average%data(1)
428 
429  ! Thomas17/9: SOLVER_TYPE = CORESOLVER_ITER%solver_1d(1)%control_parameters%integer0d(1)%value
430  solver_type = coresolver_iter%solver%index
431 
432  write(*,*)'amix=',amix,'tau=',tau,'solver_type',solver_type
433 
434 
435 
436 ! ALLOCATE (NRHO_TR(NMOD))
437  ALLOCATE (nrho_sr(nsource))
438 ! +++ Allocate local variables:
439  ALLOCATE (rho(nrho))
440  ALLOCATE (vol(nrho))
441  ALLOCATE (vpr(nrho))
442  ALLOCATE (vprm(nrho))
443  ALLOCATE (g3(nrho))
444  ALLOCATE (ne(nrho))
445  ALLOCATE (te(nrho))
446  ALLOCATE (fun_out(nrho))
447  ALLOCATE (qrad(nrho))
448  ALLOCATE (se_exp(nrho))
449 
450 ! +++ Copy data to local variables:
451  b0 = equilibrium_iter%vacuum_toroidal_field%B0(1)
452  r0 = equilibrium_iter%vacuum_toroidal_field%R0
453  ne = coreprofiles_iter%profiles_1d(1)%electrons%density
454  te = coreprofiles_iter%profiles_1d(1)%electrons%temperature
455  rho = coreprofiles_iter%profiles_1d(1)%grid%rho_tor_norm
456  b0prime = 0.0_r8
457  qrad = 0.0_r8
458  se_exp = 0.0_r8
459  tti = 0.0_r8
460 
461 ! ++++++++++++ interpolation for equilibrium ++++++++++++++++++++++++++++++++++++
462  CALL l3interp(equilibrium_iter%time_slice(1)%profiles_1d%volume, equilibrium_iter%time_slice(1)%profiles_1d%rho_tor_norm, neq, &
463  vol, rho, nrho)
464 
465  CALL l3deriv(equilibrium_iter%time_slice(1)%profiles_1d%volume, equilibrium_iter%time_slice(1)%profiles_1d%rho_tor_norm, neq, &
466  vpr, rho, nrho)
467 
468  CALL l3deriv(equilibrium_iter%time_slice(1)%profiles_1d%volume, equilibrium_iter%time_slice(1)%profiles_1d%rho_tor_norm, neq, &
469  vprm, rho, nrho)
470  CALL l3interp(equilibrium_iter%time_slice(1)%profiles_1d%gm3, equilibrium_iter%time_slice(1)%profiles_1d%rho_tor_norm, neq, &
471  g3, rho, nrho)
472 ! ++++++++++++ end interpolation for equilibrium ++++++++++++++++++++++++++
473  IF (vpr(1).LE.0.) THEN
474  vpr(1)=0.
475  END IF
476  IF (vprm(1).LE.0.) THEN
477  vprm(1)=0.
478  END IF
479 
480 
481 
482  zn = 0.0_r8
483  am = 0.0_r8
484 
485  DO iion = 1,nion
486  nzimp(iion) = 0
487  END DO
488 
489  nimp=0
490  max_nzimp=0
491  DO iion = 1, nion
492  IF (coreprofiles_iter%profiles_1d(1)%ion(iion)%multiple_states_flag.EQ.1) THEN
493  nimp = nimp+1
494  nzimp(iion)= SIZE(coreprofiles_iter%profiles_1d(1)%ion(iion)%state)
495  IF(max_nzimp.LE.nzimp(iion)) THEN
496  max_nzimp=nzimp(iion)
497  END IF
498  write (*,*) 'nimp=',nimp,'nzimp(iion)= ', nzimp(iion),max_nzimp
499  DO izimp=1,nzimp(iion)
500  IF (coreprofiles_iter%profiles_1d(1)%ion(iion)%state(izimp)%z_min.EQ.coreprofiles_iter%profiles_1d(1)%ion(iion)%state(izimp)%Z_MAX) THEN
501 ! write (*,*) ' this impurity have no bundle and will be solved solved'
502  END IF
503  END DO
504  END IF
505  END DO
506 
507 ! ************************ ATOMOC DATA ****************************
508  WRITE(*,*) 'ITM AMNSPROTO data used (via UAL)'
509  ALLOCATE(amns_ei(0:max_nzimp, nion), amns_rc(0:max_nzimp, nion), &
510  amns_eip(0:max_nzimp,nion), amns_lr(0:max_nzimp, nion), &
511  amns_br(0:max_nzimp, nion))
512  CALL imas_amns_setup(amns, version)
513  query%string = 'version'
514  CALL imas_amns_query(amns,query,answer)
515  WRITE(*,*) 'AMNS data base version = ',trim(answer%string)
516  ei_rx%string = 'EI'
517  eip_rx%string = 'EIP'
518  rc_rx%string = 'RC'
519  lr_rx%string = 'LR'
520  br_rx%string = 'BR'
521  FORMAT = '(''ZN = '',f5.2,'', IS = '',i2,'', RX = '',a,'', SRC = '',a)'
522  query%string = 'source'
523 
524  DO iion=1, nion
525 
526  zn_imp = coreprofiles_old%profiles_1d(1)%ion(iion)%element(1)%Z_N
527  am_imp = coreprofiles_old%profiles_1d(1)%ion(iion)%element(1)%a
528  write(*,*)'AM_imp=',am_imp
529 
530  am_imp = 0
531  DO izimp=0, nzimp(iion)-1
532 ! EI
533  allocate(species%components(4))
534  species%components = &
535  (/ amns_reactant_type(zn_imp, izimp, am_imp, 0), &
536  amns_reactant_type(0, -1, 0, 0), &
537  amns_reactant_type(zn_imp, izimp+1, am_imp, 1), &
538  amns_reactant_type(0, -1, 0, 1) &
539  /)
540  CALL imas_amns_setup_table(amns, ei_rx, species, amns_ei(izimp, iion))
541 
542  deallocate(species%components)
543 
544 ! LR
545  allocate(species%components(2))
546  species%components = &
547  (/ amns_reactant_type(zn_imp, izimp, am_imp, 0), &
548  amns_reactant_type(zn_imp, izimp, am_imp, 1) &
549  /)
550  CALL imas_amns_setup_table(amns, lr_rx, species, amns_lr(izimp, iion))
551  deallocate(species%components)
552  ENDDO
553 
554  DO izimp=1, nzimp(iion)
555 ! RC
556  allocate(species%components(4))
557  species%components = &
558  (/ amns_reactant_type(zn_imp, izimp, am_imp, 0), &
559  amns_reactant_type(0, -1, 0, 0), &
560  amns_reactant_type(zn_imp, izimp-1, am_imp, 1), &
561  amns_reactant_type(0, -1, 0, 1) &
562  /)
563  CALL imas_amns_setup_table(amns, rc_rx, species, amns_rc(izimp, iion))
564  deallocate(species%components)
565 
566  allocate(species%components(2))
567  species%components = &
568  (/ amns_reactant_type(zn_imp, izimp, am_imp, 0), &
569  amns_reactant_type(zn_imp, izimp, am_imp, 1) &
570  /)
571  CALL imas_amns_setup_table(amns, br_rx, species, amns_br(izimp, iion))
572  deallocate(species%components)
573 !
574  ENDDO
575  DO izimp=0,nzimp(iion)
576  !new for potential
577  allocate(species%components(2))
578  species%components = &
579  (/ amns_reactant_type(zn_imp, izimp, am_imp, 0), &
580  amns_reactant_type(zn_imp, izimp, am_imp, 1) &
581  /)
582  CALL imas_amns_setup_table(amns, eip_rx, species, amns_eip(izimp, iion))
583  deallocate(species%components)
584  END DO
585  END DO
586  write(*,*)'after AMNS'
587 ! ++++++++++++++++++++ transport ++++++++++++++
588  ALLOCATE (fun_in(SIZE(coreprofiles_iter%profiles_1d(1)%grid%rho_tor_norm)))
589  in_tr=0
590  DO inmod=1,nmod
591 ! write(*,*)'CORETRANSPORT_ITER%MODEL(INMOD)%identifier%name ',CORETRANSPORT_ITER%MODEL(INMOD)%identifier%name
592 ! write(*,*)'CORETRANSPORT_ITER%MODEL(INMOD)%identifier%index',CORETRANSPORT_ITER%MODEL(INMOD)%identifier%index
593 ! write(*,*)'CORETRANSPORT_ITER%MODEL(INMOD)%identifier%description',CORETRANSPORT_ITER%MODEL(INMOD)%identifier%description
594 ! write(*,*)'CORETRANSPORT_ITER%MODEL(INMOD)%flux_multiplier=',CORETRANSPORT_ITER%MODEL(INMOD)%flux_multiplier
595 ! write(*,*) 'INMOD=',INMOD
596 
597  IF (trim(adjustl(coretransport_iter%MODEL(inmod)%identifier%name(1))).eq."combined".and.&
598  coretransport_iter%MODEL(inmod)%flux_multiplier.eq.0) THEN
599 
600  inmod_imp = inmod
601  in_tr=in_tr+1
602  IF (in_tr.GT.1) THEN
603  write(*,*) 'WARNNING more on 1 combined transport model'
604  END IF
605  END IF
606  ENDDO
607 
608  write (*,*) 'INMOD_IMP=',inmod_imp,'in_TR=',in_tr
609 
610  nrho_tr = SIZE(coretransport_iter%MODEL(inmod_imp)%PROFILES_1D(1)%GRID_D%rho_tor_norm)
611  fun_in = coreprofiles_iter%profiles_1d(1)%grid%rho_tor_norm
612  write(*,*) ' NRHO_TR=', nrho_tr
613  nion_ct = SIZE(coretransport_iter%MODEL(inmod_imp)%PROFILES_1D(1)%ion)
614  write(*,*)'NION_CT=',nion_ct,'nion=',nion
615 
616  IF (nion_ct.NE.nion) THEN
617  write(*,*)'warning different number on the ions between coreprofiles and coretransport'
618  END IF
619 
620  nimp_ct=0
621  DO iion_ct=1,nion_ct
622  IF (coretransport_iter%MODEL(inmod_imp)%PROFILES_1D(1)%ion(iion_ct)%multiple_states_flag.EQ.1) THEN
623  nimp_ct= nimp_ct+1
624  write(*,*)'nimp_ct=',nimp_ct
625  END IF
626  END DO
627 
628  IF (nimp_ct.NE.nimp) THEN
629  write(*,*)'warning different number on the impurity between coreprofiles and coretransport'
630  END IF
631  ALLOCATE(fun_out_tr(nrho_tr))
632  ALLOCATE (nzimp_ct(nion_ct))
633 ! ++++++++++++++++ source ++++++++
634  insource_imp = 0
635  DO insource=1,nsource
636  IF (trim(adjustl(coresources_iter%source(insource)%identifier%name(1))).eq."total") THEN
637  insource_imp = insource
638  write (*,*) 'INSOURCE_IMP=',insource_imp
639  END IF
640  ENDDO
641 
642  nrho_sr(insource_imp) = SIZE(coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%GRID%rho_tor_norm)
643  nion_cs = SIZE(coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%ion)
644  IF (nion_cs.NE.nion) THEN
645  write(*,*)'warning different number on the ions between coreprofiles and coresource'
646  END IF
647 
648 
649  nimp_cs = 0
650  DO iion_cs=1,nion_cs
651  IF (coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%ion(iion_cs)%multiple_states_flag.EQ.1) THEN
652  nimp_cs= nimp_cs+1
653  write(*,*)'nimp_cs=',nimp_cs
654  END IF
655  END DO
656 
657  IF (nimp_cs.NE.nimp) THEN
658  write(*,*)'warning different number on the impurity between coresource and coresource'
659  END IF
660 
661  ALLOCATE (nzimp_cs(nion_cs))
662  ALLOCATE (coresources_new%source(insource_imp)%profiles_1d(1)%electrons%energy(nrho_sr(insource_imp)))
663  coresources_new%source(insource_imp)%profiles_1d(1)%electrons%energy = 0.0_r8
664 
665 
666 
667 ! +++ Start simulation for impurity "petla po ionow"
668 
669 
670  DO iion = 1, nion
671 
672  IF (coreprofiles_iter%profiles_1d(1)%ion(iion)%multiple_states_flag.EQ.1) THEN
673  nzimp2=nzimp(iion)+2
674 ! ALLOCATE (ANEUT(NNEUT))
675  ALLOCATE (nz1(nrho,nzimp2))
676  ALLOCATE (nzm1(nrho,nzimp2))
677  ALLOCATE (dnz1(nrho,nzimp2))
678  ALLOCATE (flux(nrho,nzimp2))
679  ALLOCATE (flux_inter(nrho,nzimp2))
680  ALLOCATE (imp_radiation(nrho,nzimp2))
681  ALLOCATE (nz_bnd(3,nzimp2))
682  ALLOCATE (nz_bnd_type(nzimp2))
683  ALLOCATE (nzsource(nrho,nzimp2))
684  ALLOCATE (diff(nrho,nzimp2))
685  ALLOCATE (vcon(nrho,nzimp2))
686  ALLOCATE (total_rad_int(nrho,nzimp2))
687  nz1 = 0.0_r8
688  nzm1 = 0.0_r8
689  diff = 0.0_r8
690  vcon = 0.0_r8
691  dnz1 = 0.0_r8
692  flux = 0.0_r8
693  flux_inter = 0.0_r8
694  imp_radiation = 0.0_r8
695  nz_bnd = 0.0_r8
696  nz_bnd_type = 0
697  nzsource = 0.0_r8
698  aa = 0.0_r8
699 
700 ! for radiation
701  ALLOCATE (lin_rad1(nrho,nzimp2))
702  ALLOCATE (brem_rad1(nrho,nzimp2))
703  ALLOCATE (jon_en1(nrho,nzimp2))
704  ALLOCATE (rec_los1(nrho,nzimp2))
705 
706  lin_rad1 = 0.0_r8
707  brem_rad1 = 0.0_r8
708  jon_en1 = 0.0_r8
709  rec_los1 = 0.0_r8
710 
711 
712  zn(iion) = coreprofiles_old%profiles_1d(1)%ion(iion)%element(1)%Z_N
713  am(iion) = coreprofiles_old%profiles_1d(1)%ion(iion)%element(1)%a
714 
715 
716 ! +++++++++ for impurity density
717  DO izimp = 1, nzimp(iion)
718  nz1(:,izimp+1) = coreprofiles_iter%profiles_1d(1)%ion(iion)%state(izimp)%density(:)
719  nzm1(:,izimp+1) = coreprofiles_old%profiles_1d(1)%ion(iion)%state(izimp)%density(:)
720  END DO
721 
722 
723 ! ++++++++ transport coefficient for impurity
724  nimp_ct=0
725  nzimp_ct=0
726  DO iion_ct = 1, nion_ct
727  IF (coretransport_iter%MODEL(inmod_imp)%PROFILES_1D(1)%ion(iion_ct)%multiple_states_flag.EQ.1.AND.&
728  zn(iion).EQ.coretransport_iter%MODEL(inmod_imp)%PROFILES_1D(1)%ion(iion_ct)%element(1)%Z_N.AND.&
729  am(iion).EQ.coretransport_iter%MODEL(inmod_imp)%PROFILES_1D(1)%ion(iion_ct)%element(1)%a) THEN
730 
731  nzimp_ct(iion_ct)= SIZE(coretransport_iter%MODEL(inmod_imp)%profiles_1d(1)%ion(iion_ct)%state)
732  nimp_ct=iion_ct
733  write(*,*)nzimp_ct(iion_ct)
734  write(*,*)'irena2'
735  IF (nzimp_ct(iion_ct).NE.nzimp(iion)) THEN
736  write(*,*)'coretransport is no for evrything ionization state'
737  END IF
738 
739  DO izimp=1,nzimp_ct(iion_ct)
740  CALL l3interp(coretransport_iter%model(inmod_imp)%profiles_1d(1)%ion(iion_ct)%state(izimp)%particles%d,&
741  coretransport_iter%MODEL(inmod_imp)%PROFILES_1D(1)%GRID_D%rho_tor, &
742  nrho_tr, fun_out, rho, nrho)
743  diff(:,izimp+1) = fun_out
744 
745  CALL l3interp(coretransport_iter%model(inmod_imp)%profiles_1d(1)%ion(iion_ct)%state(izimp)%particles%v, &
746  coretransport_iter%MODEL(inmod_imp)%PROFILES_1D(1)%GRID_D%rho_tor, &
747  nrho_tr, fun_out, rho, nrho)
748  vcon(:,izimp+1) = fun_out
749 !
750  coretransport_new%model(inmod_imp)%profiles_1d(1)%ion(iion_ct)%state(izimp)%particles%flux= 0.0_r8
751  END DO
752 
753  END IF
754  END DO
755 
756 ! +++++ for boundary condition +++++++++++
757  DO inequal=1, nequal
758 
759  IF (coresolver_iter%solver_1d(1)%equation(inequal)%primary_quantity%ion_index.EQ.iion.AND. &
760  coreprofiles_iter%profiles_1d(1)%ion(iion)%multiple_states_flag.EQ.1) THEN
761  write(*,*) 'nequal=',nequal
762 
763  izimp = coresolver_iter%solver_1d(1)%equation(inequal)%primary_quantity%state_index
764  nz_bnd_type(izimp+1) = coresolver_iter%solver_1d(1)%equation(inequal)%boundary_condition(1)%type%index
765 
766  IF (ASSOCIATED(coresolver_iter%solver_1d(1)%equation(nequal)%boundary_condition(1)%value)) THEN
767  nz_bnd(:,izimp+1) = coresolver_iter%solver_1d(1)%equation(inequal)%boundary_condition(1)%value(:)
768 ! ELSE
769 ! IF(NZ_BND_TYPE(IZIMP+1).EQ.1) THEN
770 ! NZ_BND(1,IZIMP+1)=NZM1(NRHO,IZIMP+1)
771 ! END IF
772  END IF
773  END IF
774  END DO
775 
776 
777 ! +++++++++ connect with neutrals, index ¨1¨ is for neutrals
778  ! INEUT = COREPROFILES_ITER%profiles_1d(1)%ion(IION)%neutral_index
779  ! IF (INEUT.GT.0) THEN
780  ! IF (COREPROFILES_ITER%profiles_1d(1)%neutral(INEUT)%ion_index.EQ.IION) THEN
781  ! DO ITYPE = 1, 2
782  ! FUN_IN=COREPROFILES_ITER%profiles_1d(1)%neutral(INEUT)%state(ITYPE)%density
783  ! NZ1(:,1)= NZ1(:,1)+ FUN_IN
784  ! ENDDO
785  ! NZM1(:,1) = NZ1(:,1)
786  ! END IF
787  ! END IF
788 
789 
790 ! *********** SOURCE **************
791 
792 ! For impurity give source from total
793  nimp_cs=0
794  DO iion_cs = 1, nion_cs
795  IF (coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%ion(iion_cs)%multiple_states_flag.EQ.1.AND. &
796  zn(iion).EQ.coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%ion(iion_cs)%element(1)%Z_N.AND.&
797  am(iion).EQ.coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%ion(iion_cs)%element(1)%a) THEN
798  nzimp_cs(iion_cs)= SIZE(coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%ion(iion_cs)%state)
799  write(*,*)nzimp_cs(iion_cs)
800 
801  IF (nzimp_cs(iion_cs).NE.nzimp(iion)) THEN
802  write(*,*)' WARNING coresource is no for evrything ionization state'
803  END IF
804 
805  nimp_cs=iion_cs
806 
807  DO izimp=1,nzimp_cs(iion_cs)
808  IF (associated(coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%ion(iion_cs)%state(izimp)%particles)) THEN
809  fun_out= 0.0_r8
810  CALL l3interp(coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%ion(iion_cs)%state(izimp)%particles, &
811  coresources_iter%SOURCE(insource_imp)%PROFILES_1D(1)%GRID%rho_tor_norm, &
812  nrho_sr(insource_imp), fun_out, rho, nrho)
813  nzsource(:,izimp+1) = fun_out
814  ELSE
815  nzsource(:,izimp+1)= 0.0_r8
816 ! Write(*,*)'proba'
817  END IF
818 
819  END DO
820  END IF
821  END DO
822 
823 
824 
825 ! DEALLOCATE (FUN_IN)
826 
827  max_nzimp= nzimp(iion)
828  write(*,*)'max_nzimp=',max_nzimp,'nzimp2=',nzimp2
829 
830 
831 
832 !!! +++++++++++++++ tutaj impurity one +++++++++++++++++++
833  CALL impurity_one(te, ne, nz1, nzm1, vpr, vprm, r0, b0, b0prime, diff, flux, flux_inter, rho, &
834  vcon, nrho, nzimp2, nzsource, nz_bnd, nz_bnd_type, &
835  g3, imp_radiation, se_exp, max_nzimp,amix,tau,solver_type, &
836  amns_ei(:,iion), amns_rc(:,iion), amns_lr(:,iion), amns_br(:,iion), &
837  amns_eip(:,iion), lin_rad1,brem_rad1,jon_en1,rec_los1, &
838  control_double, control_integer)
839 
840 !!! +++++++++++++++
841 
842 
843 ! for radiation
844  ALLOCATE (lin_rad(nrho))
845  ALLOCATE (brem_rad(nrho))
846  ALLOCATE (jon_en(nrho))
847  ALLOCATE (rec_los(nrho))
848 ! for source
849 ! ALLOCATE (FUN_OUT_SR(NRHO_SR(INSOURCE_IMP)))
850 
851  DO irho=1,nrho
852  lin_rad(irho) = 0.0_r8
853  brem_rad(irho) = 0.0_r8
854  jon_en(irho) = 0.0_r8
855  rec_los(irho) = 0.0_r8
856  END DO
857 
858  DO izimp = 1,nzimp(iion)
859  IF (nz1(nrho,izimp+1).LE.1.0d-200) nz1(nrho,izimp+1) = 0._r8
860  END DO
861 
862 
863  write (*,*)'irena after ti'
864 ! ++++++++++++++++++++++++++++
865 ! writing in coreprofiles
866  DO irho=1,nrho
867  DO izimp=1,nzimp(iion)
868  qrad(irho) = qrad(irho) + imp_radiation(irho,izimp+1)
869 
870  coreprofiles_new%profiles_1d(1)%ion(iion)%state(izimp)%density(irho) = nz1(irho,izimp+1)
871  coreprofiles_new%profiles_1d(1)%ion(iion)%state(izimp)%Z_MIN = izimp
872  coreprofiles_new%profiles_1d(1)%ion(iion)%state(izimp)%Z_MAX = izimp
873  coreprofiles_new%profiles_1d(1)%ion(iion)%state(izimp)%Z_AVERAGE = izimp
874  coreprofiles_new%profiles_1d(1)%ion(iion)%state(izimp)%Z_SQUARE_AVERAGE = izimp**2
875 
876 ! for radiation for one impurity - PROFILES
877  lin_rad(irho) = lin_rad(irho) + lin_rad1(irho,izimp+1)
878  brem_rad(irho) = brem_rad(irho) + brem_rad1(irho,izimp+1)
879 
880 ! for total radiation by ONE IION
881  END DO
882  CALL l3interp(-qrad,rho, nrho, &
883  coresources_new%source(insource_imp)%profiles_1d(1)%electrons%energy,&
884  coresources_new%source(insource_imp)%PROFILES_1D(1)%GRID%rho_tor_norm, nrho_sr(insource_imp))
885  END DO
886 
887 ! ***********LINE RADIATION **************
888 ! Radiation profile for evryting ionization state and for ion
889  nimp_lin=0
890  DO iion_lin = 1, nion_lin
891  write(*,*)'NION_LIN=',nion_lin
892  IF (coreradiation_new%process(inproces_lin)%profiles_1d(1)%ion(iion_lin)%multiple_states_flag.EQ.1.AND. &
893  zn(iion).EQ.coreradiation_new%process(inproces_lin)%PROFILES_1D(1)%ion(iion_lin)%element(1)%Z_N.AND.&
894  am(iion).EQ.coreradiation_new%process(inproces_lin)%PROFILES_1D(1)%ion(iion_lin)%element(1)%a) THEN
895  write(*,*) 'in radiation'
896  nzimp_lin(iion_lin)= SIZE(coreradiation_new%process(inproces_lin)%PROFILES_1D(1)%ion(iion_lin)%state)
897  write(*,*)'NZIMP_LIN(IION_LIN) = ',nzimp_lin(iion_lin)
898 
899  IF (nzimp_lin(iion_lin).NE.nzimp(iion)) THEN
900  write(*,*)' WARNING coresource is no for evrything ionization state'
901  END IF
902  nimp_lin=iion_lin
903  END IF
904  END DO
905  CALL l3interp(-lin_rad,rho, nrho, &
906  coreradiation_new%process(inproces_lin)%profiles_1d(1)%ion(nimp_lin)%emissivity,&
907  coreradiation_new%process(inproces_lin)%PROFILES_1D(1)%GRID%rho_tor_norm, nrho_lin(inproces_lin))
908 
909  DO izimp=1,nzimp(iion)
910  CALL l3interp(-lin_rad1(:,izimp+1),rho, nrho, &
911  coreradiation_new%process(inproces_lin)%profiles_1d(1)%ion(nimp_lin)%state(izimp)%emissivity,&
912  coreradiation_new%process(inproces_lin)%PROFILES_1D(1)%GRID%rho_tor_norm, nrho_lin(inproces_lin))
913 
914  END DO
915  DO izimp=1,nzimp(iion)
916  total_rad_int(1,izimp)=rho(1)* lin_rad1(1,izimp)
917  DO irho=2,nrho
918  total_rad_int(irho,izimp)= total_rad_int(irho-1,izimp)+ &
919  (vol(irho)-vol(irho-1))*0.5_r8 *(lin_rad1(irho,izimp)+ lin_rad1(irho-1,izimp))
920  END DO
921  END DO
922 
923  DO izimp=1,nzimp(iion)
924  CALL l3interp(-total_rad_int(:,izimp+1),rho, nrho, &
925  coreradiation_new%process(inproces_lin)%profiles_1d(1)%ion(nimp_lin)%state(izimp)%power_inside, &
926  coreradiation_new%process(inproces_lin)%PROFILES_1D(1)%GRID%rho_tor_norm, nrho_lin(inproces_lin))
927 
928  END DO
929 
930 
931  DO izimp=1,nzimp(iion)
932  coreradiation_new%process(inproces_lin)%profiles_1d(1)%ion(nimp_lin)%power_inside= &
933  coreradiation_new%process(inproces_lin)%profiles_1d(1)%ion(nimp_lin)%power_inside + total_rad_int(:,izimp)
934  END DO
935 
936 
937 
938 ! ***********BREM. RADIATION **************
939 ! Radiation profile for evryting ionization state and for ion
940  nimp_brem=0
941  DO iion_brem = 1, nion_brem
942  write(*,*)'NION_BREM=',nion_brem
943  IF (coreradiation_new%process(inproces_brem)%profiles_1d(1)%ion(iion_brem)%multiple_states_flag.EQ.1.AND. &
944  zn(iion).EQ.coreradiation_new%process(inproces_brem)%PROFILES_1D(1)%ion(iion_brem)%element(1)%Z_N.AND.&
945  am(iion).EQ.coreradiation_new%process(inproces_brem)%PROFILES_1D(1)%ion(iion_brem)%element(1)%a) THEN
946  write(*,*) 'in brem radiation'
947  nzimp_brem(iion_brem)= SIZE(coreradiation_new%process(inproces_brem)%PROFILES_1D(1)%ion(iion_brem)%state)
948  write(*,*)'NZIMP_BREM(IION_BREM) = ',nzimp_brem(iion_brem)
949 
950  IF (nzimp_brem(iion_brem).NE.nzimp(iion)) THEN
951  write(*,*)' WARNING coresource is no for evrything ionization state'
952  END IF
953  nimp_brem=iion_brem
954  END IF
955  END DO
956  CALL l3interp(-brem_rad,rho, nrho, &
957  coreradiation_new%process(inproces_brem)%profiles_1d(1)%ion(nimp_brem)%emissivity,&
958  coreradiation_new%process(inproces_brem)%PROFILES_1D(1)%GRID%rho_tor_norm, nrho_brem(inproces_brem))
959 
960  DO izimp=1,nzimp(iion)
961  CALL l3interp(-brem_rad1(:,izimp+1),rho, nrho, &
962  coreradiation_new%process(inproces_brem)%profiles_1d(1)%ion(nimp_brem)%state(izimp)%emissivity,&
963  coreradiation_new%process(inproces_brem)%PROFILES_1D(1)%GRID%rho_tor_norm, nrho_brem(inproces_brem))
964  END DO
965  total_rad_int =0.0
966 
967  DO izimp=1,nzimp(iion)
968  total_rad_int(1,izimp)=rho(1)* brem_rad1(1,izimp)
969  DO irho=2,nrho
970  total_rad_int(irho,izimp)= total_rad_int(irho-1,izimp)+ &
971  (vol(irho)-vol(irho-1))*0.5_r8 *(brem_rad1(irho,izimp)+ brem_rad1(irho-1,izimp))
972  END DO
973  END DO
974 
975 
976  DO izimp=1,nzimp(iion)
977  CALL l3interp(-total_rad_int(:,izimp+1),rho, nrho, &
978  coreradiation_new%process(inproces_brem)%profiles_1d(1)%ion(nimp_brem)%state(izimp)%power_inside, &
979  coreradiation_new%process(inproces_brem)%PROFILES_1D(1)%GRID%rho_tor_norm, nrho_brem(inproces_brem))
980 
981  END DO
982 
983 
984  DO izimp=1,nzimp(iion)
985  coreradiation_new%process(inproces_brem)%profiles_1d(1)%ion(nimp_brem)%power_inside= &
986  coreradiation_new%process(inproces_brem)%profiles_1d(1)%ion(nimp_brem)%power_inside + total_rad_int(:,izimp)
987  END DO
988 
989 
990 ! writing in transport
991  write (*,*)'before write coretransp' , iion
992  DO irho=1,nrho
993  DO izimp=1,nzimp(iion)
994  ALLOCATE (coretransport_new%model(inmod_imp)%profiles_1d(1)%ion(nimp_ct)%state(izimp)%particles%flux(irho))
995  coretransport_new%model(inmod_imp)%profiles_1d(1)%ion(nimp_ct)%state(izimp)%particles%flux(irho) = 0.0_r8
996  END DO
997  END DO
998 
999  DO izimp=1,nzimp(iion)
1000  IF (zn(iion).EQ.coretransport_new%MODEL(inmod_imp)%PROFILES_1D(1)%ion(nimp_ct)%element(1)%Z_N.AND.&
1001  am(iion).EQ.coretransport_new%MODEL(inmod_imp)%PROFILES_1D(1)%ion(nimp_ct)%element(1)%a) THEN
1002  fun_out_tr = 0.0_r8
1003  CALL l3interp(flux(:,izimp+1),rho, nrho, &
1004  fun_out_tr(:),&
1005  coretransport_iter%MODEL(inmod_imp)%PROFILES_1D(1)%GRID_D%rho_tor_norm, nrho_tr)
1006  coretransport_new%model(inmod_imp)%profiles_1d(1)%ion(nimp_ct)%state(izimp)%particles%flux(:) = fun_out_tr(:)
1007  ENDIF
1008  ENDDO
1009 
1010 ! ++++++++++++++++++++++++++++++++++++++++++++++
1011 ! writing in coresource
1012 
1013 
1014 
1015 ! DEALLOCATE (ANEUT)
1016  DEALLOCATE (nz1)
1017  DEALLOCATE (nzm1)
1018  DEALLOCATE (dnz1)
1019  DEALLOCATE (flux)
1020  DEALLOCATE (flux_inter)
1021  DEALLOCATE (imp_radiation)
1022  DEALLOCATE (nz_bnd)
1023  DEALLOCATE (nz_bnd_type)
1024  DEALLOCATE (nzsource)
1025  DEALLOCATE (diff)
1026  DEALLOCATE (vcon)
1027 
1028 
1029 ! for radiation
1030  DEALLOCATE (lin_rad1)
1031  DEALLOCATE (brem_rad1)
1032  DEALLOCATE (jon_en1)
1033  DEALLOCATE (rec_los1)
1034 
1035  DEALLOCATE(total_rad_int)
1036 
1037 
1038  DEALLOCATE (lin_rad)
1039  DEALLOCATE (brem_rad)
1040  DEALLOCATE (jon_en)
1041  DEALLOCATE (rec_los)
1042 
1043 
1044 ! for radiation
1045 
1046 
1047  WRITE(*,*)'END OF IMPURITY',iion
1048 
1049  END IF
1050  END DO
1051 ! +++++++++++++++++++++ finish pentla po jonow
1052  DEALLOCATE (nrho_tr)
1053  DEALLOCATE (nrho_sr)
1054  DEALLOCATE (rho)
1055  DEALLOCATE (vol)
1056  DEALLOCATE (vpr)
1057  DEALLOCATE (vprm)
1058  DEALLOCATE (g3)
1059  DEALLOCATE (ne)
1060  DEALLOCATE (te)
1061  DEALLOCATE (fun_out)
1062  DEALLOCATE (fun_out_tr)
1063  DEALLOCATE (qrad)
1064 ! DEALLOCATE (SE_EXP)
1065 ! DEALLOCATE (FUN)
1066 
1067  DEALLOCATE (nzimp)
1068  DEALLOCATE (zn)
1069  DEALLOCATE (am)
1070 
1071  WRITE (*,*) 'IMPURITY finished <==========='
1072  WRITE (*,*) ' '
1073 
1074  RETURN
1075 
1076  entry impurity_finish
1077 
1078  DO iimp = 1, SIZE(amns_ei, dim=2)
1079  DO izimp = 0, SIZE(amns_ei, dim=1)-1
1080 ! if(allocated(amns_ei(izimp, iimp))) then
1081  CALL imas_amns_finish_table(amns_ei(izimp, iimp))
1082 ! endif
1083  ENDDO
1084  ENDDO
1085  DO iimp = 1, SIZE(amns_rc, dim=2)
1086  DO izimp = 0, SIZE(amns_rc, dim=1)-1
1087 ! if(allocated(amns_rc(izimp, iimp))) then
1088  CALL imas_amns_finish_table(amns_rc(izimp, iimp))
1089 ! endif
1090  ENDDO
1091  ENDDO
1092  DO iimp = 1, SIZE(amns_eip, dim=2)
1093  DO izimp = 0, SIZE(amns_eip, dim=1)-1
1094 ! if(allocated(amns_eip(izimp, iimp))) then
1095  CALL imas_amns_finish_table(amns_eip(izimp, iimp))
1096 ! endif
1097  ENDDO
1098  ENDDO
1099  DO iimp = 1, SIZE(amns_lr, dim=2)
1100  DO izimp = 0, SIZE(amns_lr, dim=1)-1
1101 ! if(allocated(amns_lr(izimp, iimp))) then
1102  CALL imas_amns_finish_table(amns_lr(izimp, iimp))
1103 ! endif
1104  ENDDO
1105  ENDDO
1106  DO iimp = 1, SIZE(amns_br, dim=2)
1107  DO izimp = 0, SIZE(amns_br, dim=1)-1
1108 ! if(allocated(amns_br(izimp, iimp))) then
1109  CALL imas_amns_finish_table(amns_br(izimp, iimp))
1110 ! endif
1111  ENDDO
1112  ENDDO
1113  !WRITE(*,*) 'Calling ITM_AMNS_FINISH'
1114  CALL imas_amns_finish(amns)
1115  !WRITE(*,*) 'deallocating amns_ei'
1116  DEALLOCATE(amns_ei)
1117  !WRITE(*,*) 'deallocating amns_rc'
1118  DEALLOCATE(amns_rc)
1119  !WRITE(*,*) 'deallocating amns_eip'
1120  DEALLOCATE(amns_eip)
1121  !WRITE(*,*) 'deallocating amns_lr'
1122  DEALLOCATE(amns_lr)
1123  !WRITE(*,*) 'deallocating amns_br'
1124  DEALLOCATE(amns_br)
1125  !WRITE(*,*) 'returning'
1126 
1127 
1128 
1129  RETURN
1130 
1131  END SUBROUTINE impurity_imas
1132 
1133  END MODULE impurity
1134 
1135 
1136 SUBROUTINE impurity_one (TE, NE, NZ1, NZM1, VPR, VPRM ,R0 ,BT, BTPRIME ,DIFF,FLUX, FLUX_INTER, RHO,&
1137  vconv, nrho, simp2, nsource, nz_bnd, nz_bnd_type, &
1138  g1,imp_radiation, se_exp, max_nzimp,amix,tau,solver_type, &
1139  amns_ei,amns_rc,amns_lr,amns_br, &
1140  amns_eip,lin_rad1,brem_rad1,jon_en1,rec_los1, &
1141  control_double, control_integer)
1142 
1143 
1144 !--------------------------------------------------------------!
1145 ! This subroutine solves impuriy particle transport !
1146 ! equations for impurity components from 1 to NSTATE, !
1147 ! and provides: density and flux of impurity components !
1148 ! from 1 to NSTATE !
1149 !--------------------------------------------------------------!
1150 ! Source: --- !
1151 ! Developers: I.M.Ivanova-Stanik !
1152 ! Kontacts: irena.ivanova-stanik@ifpilm.pl !
1153 ! !
1154 ! Comments: might change after the IMAS !
1155 ! data stucture is finalized !
1156 ! !
1157 !--------------------------------------------------------------!
1158  use ids_types, r8 => ids_real, itm_i4 => ids_int
1159  use transport_solver_output
1160  use type_transport_solver_numerics
1161  use imas_constants_module!, only: imas_constants => constants_module
1162  use ets_math
1163 
1164 ! +++++++++ atomic data ++++++++++
1165  USE amns_types
1166  USE amns_module
1167 
1168 
1169  IMPLICIT NONE
1170 
1171 
1172 ! +++ Input/Output to numerical solver:
1173  TYPE (solver_io) :: solver !contains all I/O quantities to numerics part
1174  TYPE (transport_solver_numerics) :: numerics
1175 
1176 
1177 ! +++ Internal parameters:
1178  INTEGER :: irho, nrho !radius index, number of radial points
1179  INTEGER :: iimp, simp2,simp,isimp !index of impurity component, number of considered impurity components (max ionization state)
1180  INTEGER :: izimp, nzimp
1181  INTEGER :: nz_bnd_type(simp2) !boundary condition, type
1182  INTEGER :: max_nzimp
1183 
1184  REAL (R8) :: bt, btm, btprime !magnetic field from current time step, [T], previous time steps, [T], time derivative, [T/s]
1185  REAL (R8) :: r0
1186  REAL (R8) :: rho(nrho) !normalised minor radius, [m]
1187  real(r8) :: vpr(nrho) !V', [m^2]
1188  real(r8) :: vprm(nrho) !V' (at previous time step), [m^2]
1189  REAL (R8) :: g1(nrho)
1190  REAL (R8) :: nz1(nrho,simp2) !density from current ans previous time step, [m^-3]
1191  REAL (R8) :: nzm1(nrho,simp2) !density from previous time step,
1192  real(r8) :: flux(nrho,simp2) !ion flux, [1/s]
1193  real(r8) :: flux_inter(nrho,simp2) !ion flux, contributing to convective heat transport [1/s]
1194  REAL (R8) :: ne(nrho),te(nrho) !electron density [m^-3]
1195  REAL (R8) :: diff(nrho,simp2) !diffusion coefficient [m^2/s] and [m/s]
1196  REAL (R8) :: vconv(nrho,simp2) !pinch velocity, [m^2/s] and [m/s]
1197  REAL (R8) :: nsource(nrho,simp2)
1198  REAL (R8) :: nz_bnd(3,simp2) !boundary condition, value, [depends on NZ_BND_TYPE]
1199 
1200 ! REAL (R8) :: REAL_INDEX_T !index for linear interpolation for atomic data
1201  REAL (R8) :: alfa(nrho,simp2) !atom data: ionization coefficient(after interpolation)[m^3/s]
1202  REAL (R8) :: beta(nrho,simp2) !atom data: recombination coefficient (after interpolation) [m^3/s]
1203  REAL (R8) :: gama(nrho,simp2) !atom data: cooling coefficient (after interpolation)
1204  REAL (R8) :: slin(nrho,simp2) !atom data: linear radiation (after interpolation)
1205  REAL (R8) :: imp_radiation(nrho,simp2) !impurity radiation
1206  REAL (R8) :: potential(nrho,1:max_nzimp+2)
1207  REAL (R8) :: amix, tau !mixing factor, time step, [s]
1208 
1209 ! radiation
1210  REAL (R8) :: lin_rad1(nrho,simp2) !profile of lineradiation for one impurity
1211  REAL (R8) :: brem_rad1(nrho,simp2) !profile of bremst. for one impurity
1212  REAL (R8) :: jon_en1(nrho,simp2) !profile of jonisation energy for one impurity
1213  REAL (R8) :: rec_los1(nrho,simp2) !profile of bremst. for one impurity
1214 
1215  REAL (R8) :: se_exp(nrho)
1216 
1217 
1218  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
1219  INTEGER :: ndim !number of equations to be solved
1220  INTEGER :: solver_type !specifies the option for numerical solution
1221  REAL (R8) :: y(nrho) !function at the current amd previous time steps
1222  REAL (R8) :: ym(nrho) !function at the current amd previous time steps
1223  REAL (R8) :: dy(nrho) !derivative of function
1224  REAL (R8) :: a(nrho) !coefficients for numerical solver
1225  REAL (R8) :: b(nrho) !coefficients for numerical solver
1226  REAL (R8) :: c(nrho) !coefficients for numerical solver
1227  REAL (R8) :: d(nrho) !coefficients for numerical solver
1228  REAL (R8) :: e(nrho) !coefficients for numerical solver
1229  REAL (R8) :: f(nrho) !coefficients for numerical solver
1230  REAL (R8) :: g(nrho) !coefficients for numerical solver
1231  REAL (R8) :: h !coefficients for numerical solver
1232  REAL (R8) :: v(2), u(2), w(2) !boundary conditions for numerical solver
1233 
1234  REAL (R8) :: fun1(nrho), intfun1(nrho)
1235  REAL (R8) :: pi,ev
1236 
1237  INTEGER, INTENT(IN) :: control_integer(2) !integer control parameters
1238  REAL (R8), INTENT(IN) :: control_double(5) !real control parameters
1239 
1240  REAL (R8) :: dnz(nrho) !density gradient, [m^-4]
1241 
1242 
1243  REAL (R8) :: time
1244 
1245  REAL :: densityfluxin(nrho)
1246  REAL :: fluxin, fluxout
1247  REAL :: calkatrapez
1248 
1249  REAL ::error
1250 ! ***************
1251  INTEGER :: ifail
1252 
1253 
1254  type (amns_handle_rx_type) :: amns_ei(1:max_nzimp+1), amns_rc(1:max_nzimp+1),&
1255  amns_lr(1:max_nzimp+1), amns_br(1:max_nzimp+1),amns_eip(1:max_nzimp+1)
1256 
1257  write(*,*)'tutaj'
1258 
1259 
1260 ! +++ Set up dimensions:
1261  ndim = 1 !no coupling between density equations
1262  simp = simp2-2
1263  nzimp = simp
1264 
1265 
1266 ! +++ Set up local variables:
1267 
1268 
1269  btm =bt-btprime*tau
1270 
1271  amix = control_double(2)
1272  tau = control_double(1) ! 0.5 - Lackner's method
1273  solver_type = control_integer(1)
1274 
1275 
1276 
1277 
1278 ! +++ Allocate types for interface with numerical solver:
1279  CALL allocate_solver(ndim, nrho, solver, ifail)
1280 
1281  alfa = 0.0_r8
1282  beta = 0.0_r8
1283  gama = 0.0_r8
1284  slin = 0.0_r8
1285  potential = 0.0_r8
1286 !++++ radiation+energy losses
1287  lin_rad1 = 0.0_r8
1288  brem_rad1 = 0.0_r8
1289  jon_en1 = 0.0_r8
1290  rec_los1 = 0.0_r8
1291 
1292 
1293  write(*,*)'on impurity one','nzimp=',nzimp
1294 
1295  do izimp=1, nzimp+1
1296  if(izimp .ne. nzimp+1) then
1297  call imas_amns_rx(amns_ei(izimp),alfa(:,izimp),te,ne)
1298  write(*,*) 'alfa(nrho,izimp)=',alfa(nrho,izimp)
1299  call imas_amns_rx(amns_lr(izimp),slin(:,izimp),te,ne) ! guess: slin == line radiation loss
1300  call imas_amns_rx(amns_eip(izimp+1),potential(:,izimp),te,ne) ! guess: potential= potential of ionization
1301 
1302 
1303  endif
1304 
1305  if(izimp .ne. 1) then
1306  call imas_amns_rx(amns_rc(izimp),beta(:,izimp),te,ne)
1307  call imas_amns_rx(amns_br(izimp),gama(:,izimp),te,ne) ! guess: gama == bremsstrahlung + recombination energy loss
1308 
1309  endif
1310 
1311 
1312  enddo
1313 
1314 ! write(*,'(a,1p,100(e15.6))') 'alfa: ', (alfa(nrho/2,izimp),izimp=1,nzimp+1)
1315 ! write(*,'(a,1p,100(e15.6))') 'beta: ', (beta(nrho/2,izimp),izimp=1,nzimp+1)
1316 ! write(*,'(a,1p,100(e15.6))') 'gama: ', (gama(nrho/2,izimp),izimp=1,nzimp+1)
1317 ! write(*,'(a,1p,100(e15.6))') 'slin: ', (slin(nrho/2,izimp),izimp=1,nzimp+1)
1318 ! write(*,'(a,1p,100(e15.6))') 'potential: ', (potential(nrho/2,izimp),izimp=1,nzimp+1)
1319 
1320 
1321 
1322 
1323 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + !
1324 ! solution of particle transport equation for !
1325 ! individual state of impurity !
1326 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + !
1327 
1328  imp_loop1: DO isimp=2,simp+1 ! "1" - Neutral impurity, "2" - "+1", ...
1329 
1330 ! +++ Set equation to 'predictive' and all coefficients to zero:
1331  flag = 1
1332  y(:) = 0.0d0
1333  dy(:) = 0.0d0
1334  ym(:) = 0.0d0
1335  a(:) = 0.0d0
1336  b(:) = 0.0d0
1337  c(:) = 0.0d0
1338  d(:) = 0.0d0
1339  e(:) = 0.0d0
1340  f(:) = 0.0d0
1341  g(:) = 0.0d0
1342  h = 0.0d0
1343  v(:) = 0.0d0
1344  u(:) = 0.0d0
1345  w(:) = 0.0d0
1346 
1347 
1348 ! +++ Coefficients for ion diffusion equation in form:
1349 !
1350 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
1351 
1352  rho_loop3: DO irho=1,nrho
1353 
1354  y(irho) = nz1(irho,isimp)
1355  ym(irho) = nzm1(irho,isimp)
1356  a(irho) = vpr(irho)
1357  b(irho) = vprm(irho)
1358  c(irho) = 1.d0
1359  d(irho) = vpr(irho)*g1(irho)*diff(irho,isimp)
1360  e(irho) = vpr(irho)*g1(irho)*vconv(irho,isimp) &
1361  - btprime/2.d0/bt*rho(irho)*vpr(irho)
1362  f(irho) = vpr(irho)*nsource(irho,isimp)+vpr(irho)*( nz1(irho,isimp-1)*ne(irho)*alfa(irho,isimp-1) &
1363  + nz1(irho,isimp+1)*ne(irho)*beta(irho,isimp+1) )
1364  g(irho) = vpr(irho)*(ne(irho)*alfa(irho,isimp) + ne(irho)*beta(irho,isimp) )
1365  IF (irho.eq.nrho) then
1366 
1367 
1368 
1369  ENDIF
1370 
1371  END DO rho_loop3
1372 
1373  h = 0.5*tau
1374 ! +++ Boundary conditions for ion diffusion equation in form:
1375 !
1376 ! V*Y' + U*Y =W
1377 !
1378 ! +++ On axis:
1379 ! dNZ/drho(rho=0)=0:
1380  IF(diff(1,isimp).GT.0.d0) THEN
1381  v(1) = -diff(1,isimp)
1382  u(1) = vconv(1,isimp)
1383  w(1) = 0.d0
1384  ELSE
1385  v(1) = 1.d0
1386  u(1) = 0.d0
1387  w(1) = 0.d0
1388  END IF
1389 
1390 
1391 ! +++ At the edge:
1392 ! FIXED NZ
1393 
1394  IF(nz_bnd_type(isimp).EQ.1) THEN
1395  v(2) = 0.d0
1396  u(2) = 1.d0
1397  w(2) = nz_bnd(1,isimp)
1398  ENDIF
1399 
1400 ! FIXED grad_NZ
1401  IF(nz_bnd_type(isimp).EQ.2) THEN
1402  v(2) = 1.d0
1403  u(2) = 0.d0
1404  w(2) = nz_bnd(1,isimp)
1405  ENDIF
1406 
1407 ! FIXED L_NZ
1408  IF(nz_bnd_type(isimp).EQ.3) THEN
1409  v(2) = nz_bnd(1,isimp)
1410  u(2) = 1.d0
1411  w(2) = 0.d0
1412  ENDIF
1413 
1414 ! FIXED Flux_NZ
1415  IF(nz_bnd_type(isimp).EQ.4) THEN
1416  v(2) = -g1(nrho)*diff(nrho,isimp)*vpr(nrho)
1417  u(2) = g1(nrho)*vconv(nrho,isimp)*vpr(nrho)
1418  w(2) = nz_bnd(1,isimp)
1419  ENDIF
1420 
1421 ! Generic boundary condition
1422  IF(nz_bnd_type(isimp).EQ.5) THEN
1423  v(2) = nz_bnd(1,isimp)
1424  u(2) = nz_bnd(2,isimp)
1425  w(2) = nz_bnd(3,isimp)
1426  ENDIF
1427 
1428 
1429 ! +++ Density equation is not solved:
1430  IF(nz_bnd_type(isimp).EQ.0) THEN
1431 
1432  CALL deriv_fun(nrho,rho,y,dy) !temperature gradient
1433 
1434  flag = 0
1435 
1436  rho_loop4: DO irho=1,nrho
1437 
1438  a(irho) = 1.0d0
1439  b(irho) = 1.0d0
1440  c(irho) = 1.0d0
1441  d(irho) = 0.0d0
1442  e(irho) = 0.0d0
1443  f(irho) = 0.0d0
1444  g(irho) = 0.0d0
1445 
1446  END DO rho_loop4
1447 
1448  v(2) = 0.0d0
1449  u(2) = 1.0d0
1450  w(2) = y(nrho)
1451  END IF
1452 
1453  solver_type= 4
1454 ! +++ Defining coefficients for numerical solver:
1455  solver%TYPE = solver_type
1456  solver%EQ_FLAG(ndim) = flag
1457  solver%NDIM = ndim
1458  solver%NRHO = nrho
1459  solver%AMIX = amix
1460 
1461 
1462  rho_loop5: DO irho=1,nrho
1463 
1464  solver%RHO(irho) = rho(irho)
1465  solver%Y(ndim,irho) = y(irho)
1466  solver%DY(ndim,irho) = dy(irho)
1467  solver%YM(ndim,irho) = ym(irho)
1468  solver%A(ndim,irho) = a(irho)
1469  solver%B(ndim,irho) = b(irho)
1470  solver%C(ndim,irho) = c(irho)
1471  solver%D(ndim,irho) = d(irho)
1472  solver%E(ndim,irho) = e(irho)
1473  solver%F(ndim,irho) = f(irho)
1474  solver%G(ndim,irho) = g(irho)
1475 
1476  END DO rho_loop5
1477 
1478 
1479  solver%H = h
1480 
1481  solver%V(ndim,1) = v(1)
1482  solver%U(ndim,1) = u(1)
1483  solver%W(ndim,1) = w(1)
1484  solver%V(ndim,2) = v(2)
1485  solver%U(ndim,2) = u(2)
1486  solver%W(ndim,2) = w(2)
1487 
1488 
1489 ! +++ Solution of density diffusion equation:
1490  CALL solution_interface(solver, ifail)
1491 
1492 
1493 ! +++ New impurity density:
1494  rho_loop6: DO irho=1,nrho
1495 
1496  y(irho) = solver%Y(ndim,irho)
1497  dy(irho) = solver%DY(ndim,irho)
1498 
1499  END DO rho_loop6
1500 
1501 
1502  pi=imas_constants%pi
1503  ev=imas_constants%ev
1504 
1505 
1506 ! +++ New profiles of impurity density flux and integral source:
1507  rho_loop7: DO irho=1,nrho
1508 
1509  nzm1(irho,isimp) = nz1(irho,isimp)
1510  nz1(irho,isimp) = y(irho)
1511  dnz(irho) = dy(irho)
1512  IF (rho(irho).NE.0.0d0) THEN
1513  fun1(irho) = 1.d0/rho(irho)*(vpr(irho)*nsource(irho,isimp) &
1514  + vprm(irho)*nzm1(irho,isimp)/tau &
1515  - nz1(irho,isimp)*vpr(irho)*(1.d0/tau))
1516  ELSE IF (rho(irho).EQ.0.0d0) THEN
1517  fun1(irho) = 4.d0*pi**2*r0* (nsource(irho,isimp) &
1518  + nzm1(irho,isimp)/tau - nz1(irho,isimp)*(1.d0/tau))
1519  END IF
1520 
1521  END DO rho_loop7
1522 
1523 
1524 
1525 
1526  CALL integr_fun(nrho,rho,fun1,intfun1) !Integral source
1527 
1528  rho_loop8: DO irho=1,nrho
1529  flux_inter(irho,isimp) = intfun1(irho) &
1530  + btprime/2.d0/bt*rho(irho)*vpr(irho)*y(irho)
1531 
1532 
1533  flux(irho,isimp) = vpr(irho)*g1(irho)* &
1534  ( y(irho)*vconv(irho,isimp) - dy(irho)*diff(irho,isimp) )
1535 
1536  END DO rho_loop8
1537 
1538 
1539  END DO imp_loop1
1540 
1541 
1542  loop_change: DO isimp=2,simp+1
1543  DO irho=1,nrho
1544 
1545  nzm1(irho,isimp)=nz1(irho,isimp)
1546 
1547  end do
1548  end do loop_change
1549 
1550 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1551 ! solution of particle transport equation for
1552 ! individual state of impurity
1553 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1554 
1555  imp_loop21: DO isimp=simp+1, 2, -1 ! "1" - Neutral impurity, "2" - "+1", ...
1556 
1557 ! +++ Set equation to 'predictive' and all coefficients to zero:
1558  flag = 1
1559  y(:) = 0.0d0
1560  dy(:) = 0.0d0
1561  ym(:) = 0.0d0
1562  a(:) = 0.0d0
1563  b(:) = 0.0d0
1564  c(:) = 0.0d0
1565  d(:) = 0.0d0
1566  e(:) = 0.0d0
1567  f(:) = 0.0d0
1568  g(:) = 0.0d0
1569  h = 0.0d0
1570  v(:) = 0.0d0
1571  u(:) = 0.0d0
1572  w(:) = 0.0d0
1573 
1574 ! +++ Coefficients for ion diffusion equation in form:
1575 !
1576 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
1577 
1578  rho_loop23: DO irho=1,nrho
1579 
1580  y(irho) = nz1(irho,isimp)
1581  ym(irho) = nzm1(irho,isimp)
1582  a(irho) = vpr(irho)
1583  b(irho) = vprm(irho)
1584  c(irho) = 1.d0
1585  d(irho) = vpr(irho)*g1(irho)*diff(irho,isimp)
1586  e(irho) = vpr(irho)*g1(irho)*vconv(irho,isimp) &
1587  - btprime/2.d0/bt*rho(irho)*vpr(irho)
1588  f(irho) = vpr(irho)*nsource(irho,isimp)+vpr(irho)*(nz1(irho,isimp-1)*ne(irho)*alfa(irho,isimp-1) &
1589  + nz1(irho,isimp+1)*ne(irho)*beta(irho,isimp+1) )
1590  g(irho) = vpr(irho)*(ne(irho)*alfa(irho,isimp) + ne(irho)*beta(irho,isimp) )
1591 
1592 
1593  END DO rho_loop23
1594 
1595  h = 0.5*tau
1596 
1597 
1598 
1599 ! +++ Boundary conditions for ion diffusion equation in form:
1600 !
1601 ! V*Y' + U*Y =W
1602 !
1603 ! +++ On axis:
1604  IF(diff(1,isimp).GT.0.d0) THEN
1605  v(1) = -diff(1,isimp)
1606  u(1) = vconv(1,isimp)
1607  w(1) = 0.d0
1608  ELSE
1609  v(1) = 1.d0
1610  u(1) = 0.d0
1611  w(1) = 0.d0
1612  END IF
1613 
1614 ! +++ At the edge:
1615 ! FIXED NZ
1616 !!write(*,*)'NZ_BND_TYPE(ISIMP)=',NZ_BND_TYPE(ISIMP),'odwrotnej'
1617 !!pause
1618  IF(nz_bnd_type(isimp).EQ.1) THEN
1619  v(2) = 0.d0
1620  u(2) = 1.d0
1621  w(2) = nz_bnd(1,isimp)
1622  ENDIF
1623 
1624 ! FIXED grad_NZ
1625  IF(nz_bnd_type(isimp).EQ.2) THEN
1626  v(2) = 1.d0
1627  u(2) = 0.d0
1628  w(2) = nz_bnd(1,isimp)
1629  ENDIF
1630 
1631 ! FIXED L_NZ
1632  IF(nz_bnd_type(isimp).EQ.3) THEN
1633  v(2) = nz_bnd(1,isimp)
1634  u(2) = 1.d0
1635  w(2) = 0.d0
1636  ENDIF
1637 
1638 ! FIXED Flux_NZ
1639  IF(nz_bnd_type(isimp).EQ.4) THEN
1640  v(2) = -g1(nrho)*diff(nrho,isimp)*vpr(nrho)
1641  u(2) = g1(nrho)*vconv(nrho,isimp)*vpr(nrho)
1642  w(2) = nz_bnd(1,isimp)
1643  ENDIF
1644 
1645 ! Generic boundary condition
1646  IF(nz_bnd_type(isimp).EQ.5) THEN
1647  v(2) = nz_bnd(1,isimp)
1648  u(2) = nz_bnd(2,isimp)
1649  w(2) = nz_bnd(3,isimp)
1650  ENDIF
1651 
1652 ! +++ Density equation is not solved:
1653  IF(nz_bnd_type(isimp).EQ.0) THEN
1654 
1655  CALL deriv_fun(nrho,rho,y,dy) !temperature gradient
1656 
1657  flag = 0
1658 
1659  rho_loop24: DO irho=1,nrho
1660 
1661  a(irho) = 1.0d0
1662  b(irho) = 1.0d0
1663  c(irho) = 1.0d0
1664  d(irho) = 0.0d0
1665  e(irho) = 0.0d0
1666  f(irho) = 0.0d0
1667  g(irho) = 0.0d0
1668 
1669  END DO rho_loop24
1670 
1671  v(2) = 0.0d0
1672  u(2) = 1.0d0
1673  w(2) = y(nrho)
1674  END IF
1675 ! +++ Defining coefficients for numerical solver:
1676  solver%TYPE = solver_type
1677  solver%EQ_FLAG(ndim) = flag
1678  solver%NDIM = ndim
1679  solver%NRHO = nrho
1680  solver%AMIX = amix
1681 
1682  rho_loop25: DO irho=1,nrho
1683 
1684  solver%RHO(irho) = rho(irho)
1685  solver%Y(ndim,irho) = y(irho)
1686  solver%DY(ndim,irho) = dy(irho)
1687  solver%YM(ndim,irho) = ym(irho)
1688  solver%A(ndim,irho) = a(irho)
1689  solver%B(ndim,irho) = b(irho)
1690  solver%C(ndim,irho) = c(irho)
1691  solver%D(ndim,irho) = d(irho)
1692  solver%E(ndim,irho) = e(irho)
1693  solver%F(ndim,irho) = f(irho)
1694  solver%G(ndim,irho) = g(irho)
1695 
1696  END DO rho_loop25
1697 
1698  solver%H = h
1699  solver%V(ndim,1) = v(1)
1700  solver%U(ndim,1) = u(1)
1701  solver%W(ndim,1) = w(1)
1702  solver%V(ndim,2) = v(2)
1703  solver%U(ndim,2) = u(2)
1704  solver%W(ndim,2) = w(2)
1705 
1706 ! +++ Solution of density diffusion equation:
1707  CALL solution_interface(solver, ifail)
1708 
1709 
1710 ! +++ New impurity density:
1711  rho_loop26: DO irho=1,nrho
1712 
1713  y(irho) = solver%Y(ndim,irho)
1714  dy(irho) = solver%DY(ndim,irho)
1715 
1716  END DO rho_loop26
1717 
1718 ! +++ New profiles of impurity density flux and integral source:
1719  rho_loop27: DO irho=1,nrho
1720 
1721  nzm1(irho,isimp) = nz1(irho,isimp)
1722  nz1(irho,isimp) = y(irho)
1723  dnz(irho) = dy(irho)
1724  IF (rho(irho).NE.0.0d0) THEN
1725  fun1(irho) = 1.d0/rho(irho)*(vpr(irho)*nsource(irho,isimp) &
1726  + vprm(irho)*nzm1(irho,isimp)/tau &
1727  - nz1(irho,isimp)*vpr(irho)*(1.d0/tau))
1728  ELSE IF (rho(irho).EQ.0.0d0) THEN
1729  fun1(irho) = 4.0d0*pi**2*r0*(nsource(irho,isimp) &
1730  + nzm1(irho,isimp)/tau &
1731  - nz1(irho,isimp)*(1.d0/tau))
1732  END IF
1733 
1734  END DO rho_loop27
1735 
1736 
1737 
1738 
1739  CALL integr_fun(nrho,rho,fun1,intfun1) !Integral source
1740 
1741  rho_loop28: DO irho=1,nrho
1742 
1743 
1744  flux_inter(irho,isimp) = intfun1(irho) &
1745  + btprime/2.d0/bt*rho(irho)*vpr(irho)*y(irho)
1746 
1747 
1748  flux(irho,isimp) = vpr(irho)*g1(irho)* &
1749  ( y(irho)*vconv(irho,isimp) - dy(irho)*diff(irho,isimp) )
1750 
1751  END DO rho_loop28
1752 
1753  END DO imp_loop21
1754 
1755  fluxin = 0.0
1756  fluxout = 0.0
1757 
1758 
1759  DO irho=1,nrho
1760 
1761  densityfluxin(irho) = alfa(irho,1)*vpr(irho)*nz1(irho,1)*ne(irho)-beta(irho,2)*vpr(irho)*nz1(irho,2)*ne(irho)
1762 
1763  ENDDO
1764 
1765 
1766  fluxin = calkatrapez( densityfluxin, nrho-1, rho(nrho)-rho(nrho-1) )
1767 
1768  DO isimp=2, simp+1
1769 
1770  fluxout = fluxout - vpr(nrho)*g1(nrho)*diff(nrho,isimp)*(nz1(nrho,isimp)-nz1(nrho-1,isimp))/(rho(nrho)-rho(nrho-1))
1771 
1772  END DO
1773 
1774 
1775 
1776  do isimp=1,simp+1
1777  DO irho=1,nrho
1778 
1779 
1780  y(irho)=nz1(irho,isimp)
1781 
1782  IF (rho(irho).NE.0.0d0) THEN
1783  fun1(irho) = 1.d0/rho(irho)*(vpr(irho)*nsource(irho,isimp) &
1784  +vprm(irho)*nzm1(irho,isimp)/tau &
1785  -nz1(irho,isimp)*vpr(irho)*(1.d0/tau))
1786  ELSE IF (rho(irho).EQ.0.0d0) THEN
1787  fun1(irho) = 4.d0*pi**2*r0*(nsource(irho,isimp) &
1788  +nzm1(irho,isimp)/tau &
1789  -nz1(irho,isimp)*(1.d0/tau))
1790  END IF
1791 
1792  END DO
1793 
1794  CALL integr_fun(nrho,rho,fun1,intfun1) !Integral source
1795 
1796  CALL deriv_fun(nrho,rho,y,dy)
1797 
1798  DO irho=1,nrho
1799 
1800 
1801  flux_inter(irho,isimp) = intfun1(irho) &
1802  + btprime/2.d0/bt*rho(irho)*vpr(irho)*y(irho)
1803 
1804 
1805  flux(irho,isimp) = vpr(irho)*g1(irho)* &
1806  ( y(irho)*vconv(irho,isimp) - dy(irho)*diff(irho,isimp) )
1807 
1808  END DO
1809 
1810  END DO
1811 
1812 ! ImpurityRadiation & Electron source
1813 
1814  se_exp = 0.0_r8
1815 
1816  DO irho = 1,nrho
1817 
1818  DO isimp = 2,simp+1
1819 
1820 
1821  lin_rad1(irho,isimp) = ne(irho)*nz1(irho,isimp)*slin(irho,isimp)
1822  brem_rad1(irho,isimp) = ne(irho)*nz1(irho,isimp)*gama(irho,isimp)
1823  jon_en1(irho,isimp) = ne(irho)*nz1(irho,isimp)*potential(irho,isimp)*alfa(irho,isimp)
1824  rec_los1(irho,isimp) = ne(irho)*nz1(irho,isimp)*potential(irho,isimp-1)*beta(irho,isimp)
1825  imp_radiation(irho,isimp) = lin_rad1(irho,isimp) + brem_rad1(irho,isimp) + jon_en1(irho,isimp)*ev &
1826  - rec_los1(irho,isimp)*ev
1827 
1828  ENDDO
1829 
1830 
1831 ! SE_EXP(IRHO) = NZ1(IRHO,1)*NE(IRHO)*ALFA(IRHO,1) !Included in NEUTRAL routine
1832  DO isimp = 2,simp
1833  se_exp(irho) = se_exp(irho) &
1834  + nz1(irho,isimp)*ne(irho)*alfa(irho,isimp) - nz1(irho,isimp)*ne(irho)*beta(irho,isimp)
1835  ENDDO
1836  se_exp(irho) = se_exp(irho) - nz1(irho,simp+1)*ne(irho)*beta(irho,simp+1)
1837 
1838 
1839  ENDDO
1840 
1841 ! +++ Deallocate types for interface with numerical solver:
1842  CALL deallocate_solver(solver, ifail)
1843 
1844  RETURN
1845 
1846 
1847 
1848 END SUBROUTINE impurity_one
1849 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1850 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1851 
1852  REAL FUNCTION calkatrapez(wartosc, n, dx)
1853 
1854  INTEGER n ! liczba przedzialow
1855  REAL, DIMENSION(n+1) :: wartosc
1856  REAL (8) :: dx
1857  INTEGER i
1858 
1859  calkatrapez = 0.0
1860  calkatrapez = 0.5*(wartosc(1)+wartosc(n+1))*dx
1861  DO i = 2 , n
1862  calkatrapez = calkatrapez + wartosc(i)*dx
1863  END DO
1864 
1865 END FUNCTION calkatrapez
1866 
1867 
1868 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1869 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1870 
1871 SUBROUTINE integr2(N,X,Y,INTY) !AF 11.Oct.2011 - assumes that Y is zero f0r X.eq.0, just as INTEGR does too...
1872 !-------------------------------------------------------!
1873 ! This subroutine calculates integral of function !
1874 ! Y(X) from X=0 until X=X(N) !
1875 !-------------------------------------------------------!
1876 
1877  use ids_types, r8 => ids_real, itm_i4 => ids_int
1878 
1879  IMPLICIT NONE
1880 
1881  INTEGER :: n ! number of radial points (input)
1882  INTEGER :: i
1883 
1884  REAL (R8) :: x(n), & ! argument array (input)
1885  y(n), & ! function array (input)
1886  inty(n) ! function integral array (output)
1887 
1888  inty(1)=y(1)*x(1)/2.e0_r8
1889  DO i=2,n
1890  inty(i)=inty(i-1)+(y(i-1)+y(i))*(x(i)-x(i-1))/2.e0_r8
1891  END DO
1892 
1893  RETURN
1894 
1895 END SUBROUTINE integr2
1896 
1897 
1898 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1899 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1900 
1901 
transport_solver
Definition: transport_solver.f90:1
solution_interface
subroutine solution_interface(SOLVER, ifail)
Definition: solution_interface.f90:1
impurity
IMPURITY.
Definition: impurity.f90:8
control_parameters::get_control_parameter
Definition: control_parameters.f90:35
type_transport_solver_numerics::solver_io
The SOLVER_IO type includes variables that are used input and output to the nuermical solver of the t...
Definition: type_transport_solver_numerics.f90:63
allocate_working_ids
This module contains routines for allocation/deallocation if IDSs used in ETS.
Definition: allocate_working_ids.f90:13
calkatrapez
REAL function calkatrapez(wartosc, n, dx)
Definition: impurity.f90:1852
type_transport_solver_numerics::allocate_solver
subroutine allocate_solver(NDIM, NRHO, SOLVER, ifail)
Allocate variables in the SOLVER object.
Definition: type_transport_solver_numerics.f90:148
ets_math::deriv_fun
subroutine deriv_fun(N, X, Y, DY)
Definition: ets_math.f90:19
get_composition::get_prof_composition
subroutine get_prof_composition(COREPROFILES, COMPOSITION)
This routine detects composition of CORE_PROFILES IDS.
Definition: get_composition.f90:26
get_composition
This module contains routines for detecting plasma composition in IDSs.
Definition: get_composition.f90:13
integr2
subroutine integr2(N, X, Y, INTY)
Definition: impurity.f90:1871
ets_math::integr_fun
subroutine integr_fun(N, X, Y, INTY)
Definition: ets_math.f90:63
allocate_working_ids::allocate_coreradiation_ids
subroutine allocate_coreradiation_ids(NTIME, NPROCES, NRHO, COMPOSITION, CORERADIATION)
Definition: allocate_working_ids.f90:853
control_parameters
Definition: control_parameters.f90:1
ets_math
Definition: ets_math.f90:11
ets_plasma::plasma_composition
Definition: ets_plasma.f90:94
ets_plasma
The module defines derived types used by ETS6-CoreActor and subroutines to allocate and deallocate in...
Definition: ets_plasma.f90:26
type_transport_solver_numerics::transport_solver_numerics
Definition: type_transport_solver_numerics.f90:130
type_transport_solver_numerics
The module declares types of variables used by numerical solvers.
Definition: type_transport_solver_numerics.f90:12
impurity::impurity_imas
subroutine impurity_imas(EQUILIBRIUM_ITER, CORETRANSPORT_ITER, CORETRANSPORT_NEW, COREPROFILES_OLD, COREPROFILES_ITER, COREPROFILES_NEW, CORESOURCES_ITER, CORESOURCES_NEW, CORESOLVER_ITER, CORERADIATION_ITER, CORERADIATION_NEW, CONTROL_INTEGER, CONTROL_DOUBLE)
Definition: impurity.f90:15
type_transport_solver_numerics::deallocate_solver
subroutine deallocate_solver(SOLVER, ifail)
Dellocate variables in the SOLVER object.
Definition: type_transport_solver_numerics.f90:266
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, G1, IMP_RADIATION, SE_EXP, MAX_NZIMP, AMIX, TAU, SOLVER_TYPE, amns_ei, amns_rc, amns_lr, amns_br, amns_eip, LIN_RAD1, BREM_RAD1, JON_EN1, REC_LOS1, CONTROL_DOUBLE, CONTROL_INTEGER)
Definition: impurity.f90:1136
transport_solver_output
Definition: transport_solver_output.f90:1
Generated by   doxygen 1.8.5