ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
main_plasma.f90
Go to the documentation of this file.
1 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
15 ! + + + + + + + + + + MAIN PLASMA + + + + + + + + + + + +
16 SUBROUTINE main_plasma &
17 ! (GEOMETRY, PROFILES, TRANSPORT, SOURCES, IMPURITY, EVOLUTION, CONTROL, ifail,failstring)
18  (geometry, profiles, transport, sources, impurity, evolution, control, hyper_diff, j_boun,ifail,failstring) !AF 25.Apr.2016, 22.Aug.2016
19 
20 !-------------------------------------------------------!
21 ! This routine finds the solution for the set of !
22 ! transport equations describing the main plasma, !
23 ! with given sources and transport coefficients !
24 ! for all components. !
25 !-------------------------------------------------------!
26 ! Source: --- !
27 ! Developers: D.Kalupin, R.Stankiewicz !
28 ! Kontacts: D.Kalupin@fz-juelich.de !
29 ! Roman.Stankiewich@gmail.com !
30 ! !
31 ! Comments: equation are derived following !
32 ! Hinton&Hazeltine, Rev. Mod. Phys. !
33 ! vol. 48 (1976), pp.239-308 !
34 ! !
35 !-------------------------------------------------------!
36 
37 
38 ! +++ Declaration of variables:
39  USE ets_plasma
40 
41  USE s4_parameters !AF 12.Oct.2011 - to avoid a few interface blocks
42 
43  IMPLICIT NONE
44 
45  INTEGER :: ifail
46  CHARACTER(LEN=500) :: failstring
47  TYPE (diagnostic) :: diag !contains error messages and warnings
48 
49  INTEGER :: nrho !number of radial points (input)
50  INTEGER :: nion !number of ion species (input)
51 
52  INTEGER :: irho !current radial knot
53  INTEGER :: iion !current ion type
54 
55  REAL (R8) :: conv !convergency at current iteration
56 
57  TYPE (magnetic_geometry) :: geometry !contains all geometry quantities
58  TYPE (plasma_profiles) :: profiles !contains profiles of plasma parameters
59  TYPE (transport_coefficients) :: transport !contains profiles of trasport coefficients
60  TYPE (sources_and_sinks) :: sources !contains profiles of sources
61  TYPE (time_evolution) :: evolution !contains all parameters required by time evolution
62  TYPE (run_control) :: control !contains all parameters required by run
63  TYPE (impurity_profiles) :: impurity !contains profiles of impurities calculated
64 !by separate module
65 !control and iterations loop
66 
67 ! +++ Stabilization scheme !AF 25.Apr.2016, 22.Aug.2016
68  REAL (R8), DIMENSION(2) :: hyper_diff !hyper diffusivity
69  integer :: j_boun
70 
71 ! dy set initial ifail
72  ifail=0
73 
74 ! +++ Set up dimensions:
75  nrho = profiles%NRHO
76  nion = profiles%NION
77 
78 !AF, 11.Oct.2011
79 if (allocated(local_flux_ni_s4)) deallocate(local_flux_ni_s4)
80 if (allocated(local_flux_ni_conv_s4)) deallocate(local_flux_ni_conv_s4)
81 ALLOCATE(local_flux_ni_s4(nion))
82 ALLOCATE(local_flux_ni_conv_s4(nion))
83 local_flux_ni_s4(:) = 0.0_r8
84 local_flux_ni_conv_s4(:) = 0.0_r8
85 !AF - End
86 
87 ! +++ calculation of new current density profile:
88  CALL current &
89  (geometry,profiles,transport,sources,evolution,control,j_boun,ifail,failstring)
90 
91  if (ifail.lt.0) return
92 
93 
94 ! +++ calculation of new ion density profiles:
95  CALL ion_density &
96 ! (GEOMETRY,PROFILES,TRANSPORT,SOURCES,EVOLUTION,CONTROL,ifail,failstring)
97  (geometry,profiles,transport,sources,evolution,control,hyper_diff,ifail,failstring) !AF 25.Apr.2016, 22.Aug.2016
98 
99  if (ifail.lt.0) return
100 
101 ! +++ calculation of new electron density profile:
102  IF (control%QUASI_NEUT.GT.0.5) &
103  CALL electron_density &
104 ! (GEOMETRY,PROFILES,TRANSPORT,SOURCES,EVOLUTION,CONTROL,ifail,failstring)
105  (geometry,profiles,transport,sources,evolution,control,hyper_diff,ifail,failstring) !AF 25.Apr.2016, 22.Aug.2016
106 
107  if (ifail.lt.0) return
108 
109 ! +++ calculation of electron/ion density profile from quasi-neutrality:
110  CALL quasi_neutrality &
111  (geometry,profiles,impurity,control, ifail, failstring)
112 
113  if (ifail.lt.0) return
114 
115 
116 ! +++ calculation of new ion temperature profiles:
117  CALL temperatures &
118 ! (GEOMETRY,PROFILES,TRANSPORT,SOURCES,EVOLUTION,CONTROL,ifail,failstring)
119  (geometry,profiles,transport,sources,evolution,control,hyper_diff,ifail,failstring) !AF 25.Apr.2016, 22.Aug.2016
120 
121  if (ifail.lt.0) return
122 
123 ! +++ calculation of new toroidal rotation profiles:
124  CALL rotation &
125  (geometry,profiles,transport,sources,evolution,control,ifail,failstring)
126 
127  if (ifail.lt.0) return
128 
129 
130 !AF, 11.Oct.2011
131 DEALLOCATE(local_flux_ni_s4)
132 DEALLOCATE(local_flux_ni_conv_s4)
133 !AF - End
134 
135  RETURN
136 
137 
138 
139 END SUBROUTINE main_plasma
140 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
141 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
142 
143 
144 
145 
146 
147 
148 
149 
150 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
151 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
152 !-------------------------------------------------------!
153 ! !
154 !___________ SOLUTION OF TRANSPORT EQUATIONS: _________!
155 ! !
156 !-------------------------------------------------------!
157 ! These subroutines define generic numerical !
158 ! coefficients and boundary conditions, required by !
159 ! standardized interface to numerical solver !
160 !-------------------------------------------------------!
161 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
162 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
163 
164 
165 
166 
167 ! +++ CURRENT TRANSPORT EQUATION +++
168 
169 
170 
171 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
204 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
205 SUBROUTINE current &
206  (geometry,profiles,transport,sources,evolution,control,j_boun,ifail,failstring)
207 
208 !-------------------------------------------------------!
209 ! This subroutine solves current diffusion equation !
210 ! and provides the flux function, density of !
211 ! parallel current, density of toroidal current, !
212 ! safety factor power density due to Ohmic heating !
213 ! and parallel electric field !
214 !-------------------------------------------------------!
215 
216  USE itm_constants
217  USE ets_plasma
218  USE type_solver
219 
220  IMPLICIT NONE
221 
222  INTEGER, INTENT (INOUT) :: ifail
223  CHARACTER(LEN=500) :: failstring
224 
225 ! +++ External parameters:
226  TYPE (magnetic_geometry) :: geometry !contains all geometry quantities
227  TYPE (plasma_profiles) :: profiles !contains profiles of plasma parameters
228  TYPE (transport_coefficients) :: transport !contains profiles of trasport coefficients
229  TYPE (sources_and_sinks) :: sources !contains profiles of sources
230  TYPE (time_evolution) :: evolution !contains all parameters required by time evolution
231  TYPE (run_control) :: control !contains all parameters required by run
232 !evolution and convergence loop
233 ! +++ Input/Output to numerical solver:
234  TYPE (numerics) :: solver !contains all I/O quantities to numerics part
235 
236 ! +++ Internal parameters:
237  INTEGER :: irho, nrho !radius index, number of radial points
238  INTEGER :: nion !number of considered ion components
239 
240  REAL (R8) :: bt, btm, btprime !magnetic field from current time step, [T], previous time steps, [T], time derivative, [T/s]
241  REAL (R8) :: rgeo !major radius, [m]
242  REAL (R8) :: rho(profiles%nrho) !normalised minor radius, [m]
243  REAL (R8) :: vpr(profiles%nrho) !V', [m^2]
244  REAL (R8) :: g3(profiles%nrho)
245  REAL (R8) :: fdia(profiles%nrho) !diamagnetic function, [T*m]
246 
247  REAL (R8) :: psi(profiles%nrho), psim(profiles%nrho) !flux function from current ans previous time step, [V*s]
248  REAL (R8) :: dpsi(profiles%nrho), dpsim(profiles%nrho) !AF - 25.Sep.2014
249 
250  REAL (R8) :: psi_bnd(2,3) !boundary condition, value, [depend on PSI_BND_TYPE]
251  INTEGER :: psi_bnd_type(2) !boundary condition, type
252 
253  REAL (R8) :: sigma(profiles%nrho) !plasma parallel conductivity, [(Ohm*m)^-1]
254 
255  REAL (R8) :: qsf(profiles%nrho) !safety factor
256  REAL (R8) :: curr_tor(profiles%nrho) !current density, toroidal, [A/m^2]
257  REAL (R8) :: curr_par(profiles%nrho) !current density, parallel, [A/m^2]
258  REAL (R8) :: curr_ni_exp(profiles%nrho) !total non inductive current, PSI independent component, [A/m^2]
259  REAL (R8) :: curr_ni_imp(profiles%nrho) !total non inductive current, component proportional to PSI, [A/m^2/V/s]
260  REAL (R8) :: qoh(profiles%nrho) !Ohmic heating, [W/m^3]
261  REAL (R8) :: qoh_tot !Ohmic heating, [W]
262  REAL (R8) :: curr_tot, curr_ni !Total Ohmic and non-inductive currents [A]
263  REAL (R8) :: e_par(profiles%nrho) !parallel electric field,, [V/m]
264 
265  REAL (R8) :: fun1(profiles%nrho), dfun1(profiles%nrho)
266  REAL (R8) :: fun2(profiles%nrho)
267  REAL (R8) :: fun3(profiles%nrho), intfun3(profiles%nrho)
268  REAL (R8) :: fun4(profiles%nrho), dfun4(profiles%nrho)
269  REAL (R8) :: fun5(profiles%nrho), dfun5(profiles%nrho)
270  REAL (R8) :: fun7(profiles%nrho), intfun7(profiles%nrho)
271  REAL (R8) :: fun71(profiles%nrho), dfun71(profiles%nrho)
272  REAL (R8) :: fun72(profiles%nrho), dfun72(profiles%nrho)
273 
274 
275  REAL (R8) :: amix, tau !mixing factor, time step, [s]
276 
277  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
278  INTEGER :: ndim !number of equations to be solved
279  INTEGER :: solver_type !specifies the option for numerical solution
280  REAL (R8) :: y(profiles%nrho) !function at the current amd previous time steps
281  REAL (R8) :: ym(profiles%nrho) !function at the current amd previous time steps
282  REAL (R8) :: dy(profiles%nrho) !derivative of function
283  REAL (R8) :: dym(profiles%nrho) !derivative of function at previous time step !AF - 25.Sep.2014
284  REAL (R8) :: a(profiles%nrho) !coefficients for numerical solver
285  REAL (R8) :: b(profiles%nrho) !coefficients for numerical solver
286  REAL (R8) :: c(profiles%nrho) !coefficients for numerical solver
287  REAL (R8) :: d(profiles%nrho) !coefficients for numerical solver
288  REAL (R8) :: e(profiles%nrho) !coefficients for numerical solver
289  REAL (R8) :: f(profiles%nrho) !coefficients for numerical solver
290  REAL (R8) :: g(profiles%nrho) !coefficients for numerical solver
291  REAL (R8) :: h !coefficients for numerical solver
292  REAL (R8) :: v(2), u(2), w(2) !boundary conditions for numerical solver
293 
294  REAL (R8) :: dpc1(profiles%nrho), dpc2(profiles%nrho), solution_norm
295 
296  INTEGER :: dy_from_solver !AF, 4.Oct.2011
297  INTEGER :: fix_axis !AF, 5.Oct.2011
298  integer :: j_boun
299 
300  dy_from_solver = 0 !AF 4.Oct.2011 - to more easily switch between the 2 ways of getting DPSI
301  fix_axis = 1 !AF, 5.Oct.2011 - to more easily switch on/off the axis corrections
302 
303 ! +++ Set up dimensions:
304  ndim = 1
305  nrho = profiles%NRHO
306  nion = profiles%NION
307 
308 
309 ! +++ Allocate types for interface with numerical solver:
310  CALL allocate_numerics(ndim, nrho, solver, ifail)
311 
312 
313 ! +++ Set equation to 'predictive' and all coefficients to zero:
314  flag = 1
315  y(:) = 0.0e0_r8
316  dy(:) = 0.0e0_r8
317  ym(:) = 0.0e0_r8
318  dym(:) = 0.0e0_r8 !AF - 25.Sep.2014
319  a(:) = 0.0e0_r8
320  b(:) = 0.0e0_r8
321  c(:) = 0.0e0_r8
322  d(:) = 0.0e0_r8
323  e(:) = 0.0e0_r8
324  f(:) = 0.0e0_r8
325  g(:) = 0.0e0_r8
326  h = 0.0e0_r8
327  v(:) = 0.0e0_r8
328  u(:) = 0.0e0_r8
329  w(:) = 0.0e0_r8
330 
331 
332 ! +++ Set up local variables:
333  amix = control%AMIX
334  tau = control%TAU
335  solver_type = control%SOLVER_TYPE
336 
337  bt = geometry%BGEO
338  btm = evolution%BTM
339  btprime = (bt-btm)/tau
340  rgeo = geometry%RGEO
341 
342 
343 
344 ! +++ Set up boundary condition values:
345  psi_bnd_type(2) = profiles%PSI_BND_TYPE
346  psi_bnd(2,1) = profiles%PSI_BND(1)
347  psi_bnd(2,2) = profiles%PSI_BND(2)
348  psi_bnd(2,3) = profiles%PSI_BND(3)
349 
350 
351 
352 ! +++ Set up profiles:
353  rho_loop4: DO irho=1,nrho
354  rho(irho) = geometry%RHO(irho)
355  vpr(irho) = geometry%VPR(irho)
356  fdia(irho) = geometry%FDIA(irho)
357  g3(irho) = geometry%G3(irho)
358 
359  psi(irho) = profiles%PSI(irho)
360  dpsi(irho) = profiles%DPSI(irho)
361  psim(irho) = evolution%PSIM(irho)
362  dpsim(irho) = evolution%DPSIM(irho)
363  qsf(irho) = profiles%QSF(irho)
364 
365  sigma(irho) = transport%SIGMA(irho)
366 
367  curr_tor(irho) = profiles%CURR_TOR(irho)
368  curr_par(irho) = profiles%CURR_PAR(irho)
369 
370  curr_ni_exp(irho) = 0.e0_r8
371  curr_ni_imp(irho) = 0.e0_r8
372 
373  curr_ni_exp(irho) = curr_ni_exp(irho) + sources%CURR_EXP(irho)
374  curr_ni_imp(irho) = curr_ni_imp(irho) + sources%CURR_IMP(irho)
375 
376  fun1(irho) = sigma(irho)*itm_mu0*(rho(irho)/fdia(irho))**2
377  fun2(irho) = vpr(irho)*g3(irho)/4.e0_r8/itm_pi**2
378 
379  END DO rho_loop4
380 
381  CALL derivn(nrho,rho,fun1,dfun1) !Derivation of function FUN1
382 
383 
384 
385 
386 ! +++ Coefficients for for current diffusion
387 ! equation in form:
388 !
389 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
390 
391  rho_loop5: DO irho=1,nrho
392  y(irho) = psi(irho)
393  dy(irho) = dpsi(irho) !AF - 25.Sep.2014
394  ym(irho) = psim(irho)
395  dym(irho)= dpsim(irho) !AF - 25.Sep.2014
396 
397  a(irho) = sigma(irho)
398  b(irho) = sigma(irho)
399  c(irho) = itm_mu0*bt*rho(irho)/fdia(irho)**2
400  d(irho) = vpr(irho)/4.e0_r8/itm_pi**2*g3(irho)/fdia(irho)
401  e(irho) = -fun1(irho)*btprime/2.e0_r8
402  IF (rho(irho).NE.0.e0_r8) THEN
403  f(irho) = vpr(irho)/2.e0_r8/itm_pi/rho(irho)*curr_ni_exp(irho)
404  g(irho) = vpr(irho)/2.e0_r8/itm_pi/rho(irho)*curr_ni_imp(irho) &
405  + btprime/2.e0_r8*sigma(irho)*dfun1(irho)
406  ELSE
407  f(irho) = 1.e0_r8/2.e0_r8/itm_pi*curr_ni_exp(irho)
408  g(irho) = 1.e0_r8/2.e0_r8/itm_pi*curr_ni_imp(irho) &
409  + btprime/2.e0_r8*sigma(irho)*dfun1(irho)
410 
411  ENDIF
412  END DO rho_loop5
413 
414  h = tau
415 
416 
417 
418 ! +++ Boundary conditions for current diffusion
419 ! equation in form:
420 !
421 ! V*Y' + U*Y =W
422 !
423 ! +++ On axis:
424 ! dpsi/drho(rho=0)=0
425  v(1) = 1.e0_r8
426  u(1) = 0.e0_r8
427  w(1) = 0.e0_r8
428 
429 ! +++ At the edge:
430 
431 ! FIXED psi
432  IF(psi_bnd_type(2).EQ.1) THEN
433  v(2) = 0.e0_r8
434  u(2) = 1.e0_r8
435  w(2) = psi_bnd(2,1)
436  ENDIF
437 
438 ! FIXED total current
439  IF(psi_bnd_type(2).EQ.2) THEN
440  v(2) = 1.e0_r8
441  u(2) = 0.e0_r8
442 ! dpc test (this could be wrong!
443  w(2) = - itm_mu0/fun2(nrho)*psi_bnd(2,1)
444  write(*,*) 'Current equivalent current for B.C. ', &
445  -(y(nrho)-y(nrho-1))/(rho(nrho)-rho(nrho-1))*fun2(nrho)/itm_mu0
446 ! cpd
447  ENDIF
448 
449 ! FIXED loop voltage
450  IF(psi_bnd_type(2).EQ.3) THEN
451  v(2) = 0.e0_r8
452  u(2) = 1.e0_r8
453  w(2) = tau*psi_bnd(2,1)+psim(nrho)
454  ENDIF
455 
456 ! Generic boundary condition
457  IF(psi_bnd_type(2).EQ.4) THEN
458  v(2) = psi_bnd(2,1)
459  u(2) = psi_bnd(2,2)
460  w(2) = psi_bnd(2,3)
461  ENDIF
462 
463 
464 
465 ! +++ Current equation is not solved: !Interpretative value of safety factor should be given
466  IF(psi_bnd_type(2).EQ.0) THEN
467 
468  rho_loop6: DO irho =1,nrho
469  IF (qsf(irho).NE.0.e0_r8) THEN
470  dy(irho) = 2.e0_r8*itm_pi*bt*rho(irho)/qsf(irho)
471  fun3(irho) = 2.e0_r8*itm_pi*bt/qsf(irho)
472  END IF
473  END DO rho_loop6
474 
475  CALL integr(nrho,rho,fun3,y)
476 
477  flag = 0
478 
479  rho_loop7: DO irho=1,nrho
480  a(irho) = 1.0e0_r8
481  b(irho) = 1.0e0_r8
482  c(irho) = 1.0e0_r8
483  d(irho) = 0.0e0_r8
484  e(irho) = 0.0e0_r8
485  f(irho) = 0.0e0_r8
486  g(irho) = 0.0e0_r8
487  END DO rho_loop7
488 
489  v(2) = 0.0e0_r8
490  u(2) = 1.0e0_r8
491  w(2) = y(nrho)
492  END IF
493 
494 
495 
496 
497 
498 
499 ! +++ Defining coefficients for numerical solver:
500  solver%TYPE = solver_type
501  solver%EQ_FLAG(ndim) = flag
502  solver%NDIM = ndim
503  solver%NRHO = nrho
504  solver%AMIX = amix
505  solver%DERIVATIVE_FLAG(1) = 2
506 
507  rho_loop8: DO irho=1,nrho
508 
509  solver%RHO(irho) = rho(irho)
510 
511  solver%Y(ndim,irho) = y(irho)
512  solver%DY(ndim,irho) = dy(irho)
513  solver%YM(ndim,irho) = ym(irho)
514 
515  solver%A(ndim,irho) = a(irho)
516  solver%B(ndim,irho) = b(irho)
517  solver%C(ndim,irho) = c(irho)
518  solver%D(ndim,irho) = d(irho)
519  solver%E(ndim,irho) = e(irho)
520  solver%F(ndim,irho) = f(irho)
521  solver%G(ndim,irho) = g(irho)
522 
523  END DO rho_loop8
524 
525  solver%H = h
526 
527  solver%V(ndim,1) = v(1)
528  solver%U(ndim,1) = u(1)
529  solver%W(ndim,1) = w(1)
530  solver%V(ndim,2) = v(2)
531  solver%U(ndim,2) = u(2)
532  solver%W(ndim,2) = w(2)
533 
534 
535 
536 ! +++ Solution of current diffusion equation:
537  CALL solution_interface(solver, ifail)
538 
539 ! dy check for nans in solution, return if true)
540  IF (any( isnan(solver%y))) then
541  failstring='Error in the current diffusion equeation: nans in the solution vector'
542  ifail=-1
543  return
544  end if
545 
546  solution_norm=max( &
547  maxval(abs( solver%c(ndim,:)*solver%a(ndim,:)*solver%y(ndim,:) )), &
548  maxval(abs( solver%c(ndim,:)*solver%b(ndim,:)*solver%ym(ndim,:) )), &
549  maxval(abs( solver%h*solver%d(ndim,:)*solver%dy(ndim,:) )), &
550  maxval(abs( solver%h*solver%e(ndim,:)*solver%y(ndim,:) )), &
551  maxval(abs( solver%h*solver%c(ndim,:)*solver%f(ndim,:) )), &
552  maxval(abs( solver%h*solver%c(ndim,:)*solver%g(ndim,:)*solver%y(ndim,:) )))
553 
554  dpc1 = solver%c(ndim,:)*solver%a(ndim,:)*solver%y(ndim,:) &
555  - solver%c(ndim,:)*solver%b(ndim,:)*solver%ym(ndim,:) &
556  - solver%h*solver%d(ndim,:)*solver%dy(ndim,:) &
557  + solver%h*solver%e(ndim,:)*solver%y(ndim,:) &
558  - solver%h*solver%c(ndim,:)*solver%f(ndim,:) &
559  + solver%h*solver%c(ndim,:)*solver%g(ndim,:)*solver%y(ndim,:)
560 
561 ! write(*,*) 'CURRENT: SOLUTION_NORM = ', solution_norm
562 ! write(*,*) 'CURRENT: RHS-LHS = ', DPC1
563 ! write(*,*) 'CURRENT: norm(RHS-LHS) = ', DPC1/solution_norm
564 
565  IF (psi_bnd_type(2).NE.0) THEN
566 ! +++ New magnetic flux function and current density:
567  intfun7 = 0.d0
568 
569  y(:) = solver%Y(ndim,:)
570 ! the following is experimental
571  IF (fix_axis.EQ.1) THEN !AF 5.Oct.2011
572  call f_par_axis(nrho,rho,y)
573  END IF
574 ! end of the experiment
575  rho_loop9: DO irho=1,nrho
576  IF (dy_from_solver.EQ.0) THEN !AF 4.Oct.2011
577 !dpc
578  dpc2(irho) = (c(irho)*((a(irho)*y(irho)-b(irho)*ym(irho))/h &
579  + g(irho)*y(irho) - f(irho)) + e(irho)*y(irho)) / d(irho)
580 !cpd
581  fun7(irho) = c(irho)*((a(irho)*y(irho)-b(irho)*ym(irho))/h &
582  + g(irho)*y(irho) - f(irho))
583  IF (rho(irho).NE.0.e0_r8) THEN
584  intfun7(irho) = intfun7(irho-1) &
585  + (fun7(irho-1)+fun7(irho))*(rho(irho)-rho(irho-1))/2.d0
586  END IF
587 !DPC: try and prevent division by 0
588  IF(d(irho).NE.0.0_r8) THEN
589  dy(irho) = intfun7(irho)/d(irho) + e(irho)/d(irho)*y(irho)
590  ELSE
591  dy(irho) = 0.0_r8
592  WRITE(*,*) 'Warning: DY(',irho,') set to 0'
593  END IF
594  ELSE !AF 4.Oct.2011
595 !dpc-temp
596  dy(irho) = solver%DY(ndim,irho) !AF 4.Oct.2011
597 !dpc-end
598  END IF !AF 4.Oct.2011
599 
600  END DO rho_loop9
601 
602  IF (dy_from_solver.EQ.0) THEN !AF 5.Oct.2011 - this is useless if the solver DY is being used already
603 ! write(*,*) '### ', DY(NRHO), ' replaced by ',SOLVER%DY(NDIM,NRHO)
604  dy(nrho) = solver%DY(ndim,nrho)
605  END IF
606 
607  rho_loop10: DO irho=1,nrho
608  fun4(irho) = fun2(irho)*dy(irho)
609  fun5(irho) = fun2(irho)*dy(irho)*rgeo*bt/fdia(irho)
610  END DO rho_loop10
611 
612  CALL derivn(nrho, rho, fun4, dfun4) !Derivation of function FUN4
613  CALL derivn(nrho, rho, fun5, dfun5) !Derivation of function FUN5
614 
615 ! +++ New profiles of plasma parameters obtained
616 ! from current diffusion equation:
617  rho_loop11: DO irho=1,nrho
618  psi(irho) = y(irho)
619  if(rho(irho).ne.0.e0_r8) then
620  qsf(irho) = 2.e0_r8*itm_pi*bt*rho(irho)/dy(irho)
621 
622  IF(vpr(irho).NE.0.0_r8) THEN
623  curr_tor(irho) = - 2.e0_r8*itm_pi*rgeo/itm_mu0/vpr(irho)*dfun4(irho)
624  curr_par(irho) = - 2.e0_r8*itm_pi/rgeo/itm_mu0/vpr(irho) &
625  *(fdia(irho)/bt)**2*dfun5(irho)
626  END IF
627  endif
628  END DO rho_loop11
629 
630 
631  if (j_boun.eq.1) then
632 ! artifisially set curr_tor and curr_par near the boundary
633 ! to avoid jump, use 5 points
634 
635  do irho=nrho-4,nrho
636  curr_tor(irho)=curr_tor(nrho-5)
637  curr_par(irho)=curr_par(nrho-5)
638  enddo
639  end if
640 
641  IF (psi(1).NE.0.0_r8) psi = psi - psi(1)
642 
643  IF(rho(1).EQ.0.e0_r8) THEN
644  IF (fix_axis.EQ.1) THEN !AF 5.Oct.2011
645  CALL f_axis(nrho, rho, qsf)
646  CALL f_axis(nrho, rho, curr_par)
647  CALL f_axis(nrho, rho, curr_tor)
648  END IF
649  END IF
650 
651  dpsi = dy !AF - 25.Sep.2014
652 
653  END IF
654 
655 !dy comment all modification for the current density calculations
656 ! until full solution will be implemented
657 
658 ! +++ New magnetic flux function if equation is not solved:
659  IF (psi_bnd_type(2).EQ.0) THEN
660  psi = 0.d0
661 !dy curr_tor =0.0D0
662 !dy curr_par =0.0d0
663  rho_loop12: DO irho=1,nrho
664  fun7(irho) = 2.e0_r8*itm_pi*bt/qsf(irho)*rho(irho)
665 !dy fun71(irho) =fun2(irho)*fun7(irho)
666 !dy fun72(irho) =fun71(irho)*RGEO*BT/FDIA(IRHO)
667  IF (rho(irho).NE.0.e0_r8) THEN
668 ! integration
669  psi(irho) = psi(irho-1) &
670  + (fun7(irho-1)+fun7(irho))*(rho(irho)-rho(irho-1))/2.d0
671  END IF
672  END DO rho_loop12
673 
674 !dy call derivn(nrho,rho,fun71,dfun71)
675 !dy call derivn(nrho,rho,fun72,dfun72)
676  do irho=1,nrho
677  if (vpr(irho).ne.0.0d0) then
678 !dy curr_tor(irho)=-2.e0_R8*ITM_PI*RGEO/ITM_MU0/VPR(IRHO)*dfun71(irho)
679 !dy curr_par(irho)=-2.e0_R8*ITM_PI/RGEO/ITM_MU0/VPR(IRHO) &
680 !dy *(FDIA(IRHO)/BT)**2*dfun72(irho)
681  end if
682  end do
683 
684  IF(rho(1).EQ.0.e0_r8) THEN
685  IF (fix_axis.EQ.1) THEN !AF 5.Oct.2011
686  CALL f_axis(nrho, rho, psi)
687 !dy CALL F_AXIS (NRHO, RHO, curr_tor)
688 !dy CALL F_AXIS (NRHO, RHO, curr_par)
689  END IF
690  END IF
691 
692  CALL derivn(nrho,rho,psi,dpsi) !AF - 25.Sep.2014
693 
694  END IF
695 
696 
697 
698 
699 ! +++ Radial electric field and ohmic heating:
700  rho_loop13: DO irho=1,nrho
701  if (sigma(irho).eq.0.0) then
702  e_par(irho)=0.0
703  else
704  e_par(irho) = ( curr_par(irho) - curr_ni_exp(irho) &
705  - curr_ni_imp(irho)*y(irho) ) / sigma(irho)
706  end if
707  qoh(irho) = sigma(irho)*e_par(irho)**2 * control%ohmic_heating_multiplier
708  END DO rho_loop13
709 
710 
711 !+++ Current diagnostics:
712  fun7 = vpr * curr_par / 2.0e0_r8 / itm_pi * bt / fdia**2
713  CALL integr2(nrho, rho, fun7, intfun7)
714  curr_tot = intfun7(nrho) * fdia(nrho)
715  fun7 = vpr * (curr_ni_exp + curr_ni_imp * psi) / 2.0e0_r8 / itm_pi * bt / fdia**2
716  CALL integr2(nrho, rho, fun7, intfun7)
717  curr_ni = intfun7(nrho) * fdia(nrho)
718 
719  fun7 = qoh * vpr
720  CALL integr2(nrho, rho, fun7, intfun7)
721  qoh_tot = intfun7(nrho)
722 
723 
724  WRITE (*,*) '*******************'
725  WRITE (*,*) '*******************'
726  WRITE (*,*) '***** QOH =', qoh_tot/1.d6, '[MW]'
727  WRITE (*,*) '***** CURR =', curr_tot/1.d6,'[MA]'
728  WRITE (*,*) '***** CURR_NI =', curr_ni/1.d6, '[MA]'
729  WRITE (*,*) '***** F_NI =', curr_ni/curr_tot
730  WRITE (*,*) '*******************'
731  WRITE (*,*) '*******************'
732 
733 
734 
735 ! +++ Return new profiles to the work flow:
736 
737  profiles%PSI = psi
738  profiles%DPSI = dpsi !AF - 25.Sep.2014
739  profiles%QSF = qsf
740  profiles%CURR_TOR = curr_tor
741  profiles%CURR_PAR = curr_par
742  profiles%E_PAR = e_par
743  sources%QOH = qoh
744  profiles%QOH = qoh ! added by request of AF; not sure if units are OK
745  profiles%INT_QOH = intfun7
746  profiles%JNI = curr_ni_exp + curr_ni_imp * psi
747  profiles%JOH = curr_par - curr_ni_exp - curr_ni_imp * psi !AF - 25.Set.2014 - this is how it's done in 4.10a
748  profiles%SIGMA = sigma
749 
750 
751 ! +++ Deallocate types for interface with numerical solver:
752  CALL deallocate_numerics(solver, ifail)
753 
754 
755 
756  RETURN
757 
758 
759 
760 END SUBROUTINE current
761 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
762 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
763 
764 
765 
766 
767 
768 
769 
770 
771 
772 
773 ! +++ PARTICLE TRANSPORT EQUATIONS +++
774 
775 
776 
777 
778 
779 
780 
781 
782 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
810 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
811 SUBROUTINE ion_density &
812 ! (GEOMETRY,PROFILES,TRANSPORT,SOURCES,EVOLUTION,CONTROL,ifail,failstring)
813  (geometry,profiles,transport,sources,evolution,control,hyper_diff,ifail,failstring) !AF 25.Apr.2016, 22.Aug.2016
814 
815 !-------------------------------------------------------!
816 ! This subroutine solves ion particle transport !
817 ! equations for ion components from 1 to NION, !
818 ! and provides: density and flux of ion components !
819 ! from 1 to NION !
820 !-------------------------------------------------------!
821 
822  USE ets_plasma
823  USE type_solver
824 !!! dpc --- missing USE PEDESTAL
825 
826  USE s4_parameters !AF 12.Oct.2011 - to avoid a few interface blocks
827 
828  IMPLICIT NONE
829 
830  INTEGER, INTENT (INOUT) :: ifail
831  CHARACTER(LEN=500) :: failstring
832 
833 ! +++ External parameters:
834  TYPE (magnetic_geometry) :: geometry !contains all geometry quantities
835  TYPE (plasma_profiles) :: profiles !contains profiles of plasma parameters
836  TYPE (transport_coefficients) :: transport !contains profiles of trasport coefficients
837  TYPE (sources_and_sinks) :: sources !contains profiles of sources
838  TYPE (time_evolution) :: evolution !contains all parameters required by time evolution
839  TYPE (run_control) :: control !contains all parameters required by run
840 !evolution and convergence loop
841 ! +++ Input/Output to numerical solver:
842  TYPE (numerics) :: solver !contains all I/O quantities to numerics part
843 
844 ! +++ Internal parameters:
845  INTEGER :: irho, nrho !radius index, number of radial points
846  INTEGER :: iion, nion !index of ion component, number of considered ion components
847  INTEGER :: imodel, nmodel !model index, number of transport modules used
848 
849  REAL (R8) :: bt, btm, btprime !magnetic field from current time step, [T], previous time steps, [T], time derivative, [T/s]
850  REAL (R8) :: rho(profiles%nrho) !normalised minor radius, [m]
851  REAL (R8) :: vpr(profiles%nrho) !V', [m^2]
852  REAL (R8) :: vprm(profiles%nrho) !V' (at previous time step), [m^2]
853  REAL (R8) :: g1(profiles%nrho)
854 
855  REAL (R8) :: ni(profiles%nrho),nim(profiles%nrho) !density from current ans previous time step, [m^-3]
856  REAL (R8) :: dni(profiles%nrho) !density gradient, [m^-4]
857  REAL (R8) :: dnim(profiles%nrho) !density gradient at previous time step, [m^-4] !AF - 25.Sep.2014
858  REAL (R8) :: flux(profiles%nrho) !ion flux, [1/s]
859  REAL (R8) :: int_source(profiles%nrho) !integral of source [1/s]
860  REAL (R8) :: flux_ni_conv(profiles%nrho) !ion flux, contributing to convective heat transport [1/s]
861 
862  REAL (R8) :: diff(profiles%nrho),vconv(profiles%nrho) !diffusion coefficient and pinch velocity, [m^2/s] and [m/s]
863  REAL (R8) :: diff_mod(transport%nrho,3) !diffusion coefficient due to particular transport model [m^2/s]
864  REAL (R8) :: vconv_mod(transport%nrho,3) !pinch velocity due to particular transport model [m/s]
865  REAL (R8) :: c1(3) !coefficient in front of convective term in energy balance equation
866 
867  REAL (R8) :: si_exp(profiles%nrho) !explicit particle source, [1/(m^3*s)]
868  REAL (R8) :: si_imp(profiles%nrho) !implicit particle source (proportional to density), [1/s]
869 
870  REAL (R8) :: ni_bnd(2,3) !boundary condition, value, [depends on NI_BND_TYPE]
871  INTEGER :: ni_bnd_type(2) !boundary condition, type
872  REAL (R8) :: rho_bnd !boundary location
873  INTEGER :: nrho_bnd !boundary location index
874 
875  REAL (R8) :: amix, tau !mixing factor, time step, [s]
876 
877  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
878  INTEGER :: ndim !number of equations to be solved
879  INTEGER :: solver_type !specifies the option for numerical solution
880  REAL (R8) :: y(profiles%nrho) !function at the current amd previous time steps
881  REAL (R8) :: ym(profiles%nrho) !function at the current amd previous time steps
882  REAL (R8) :: dy(profiles%nrho) !derivative of function
883  REAL (R8) :: dym(profiles%nrho) !derivative of function at previous time step !AF - 25.Sep.2014
884  REAL (R8) :: a(profiles%nrho) !coefficients for numerical solver
885  REAL (R8) :: b(profiles%nrho) !coefficients for numerical solver
886  REAL (R8) :: c(profiles%nrho) !coefficients for numerical solver
887  REAL (R8) :: d(profiles%nrho) !coefficients for numerical solver
888  REAL (R8) :: e(profiles%nrho) !coefficients for numerical solver
889  REAL (R8) :: f(profiles%nrho) !coefficients for numerical solver
890  REAL (R8) :: g(profiles%nrho) !coefficients for numerical solver
891  REAL (R8) :: h !coefficients for numerical solver
892  REAL (R8) :: v(2), u(2), w(2) !boundary conditions for numerical solver
893 
894  REAL (R8) :: fun1(profiles%nrho), intfun1(profiles%nrho)
895 
896  !AF, 11.Oct.2011 - auxiliary variables to be used in the axis boundary condition of solver 4
897  REAL (R8), ALLOCATABLE :: local_int_source_s4(:) !auxiliary variable with the integral of sources divided by the metrics coefficients Vprime and G1
898  REAL (R8) :: local_fun1_s4, local_intfun1_s4
899  !AF - End
900 
901 ! +++ Stabilization scheme !AF 25.Apr.2016, 22.Aug.2016
902  REAL (R8), DIMENSION(2) :: hyper_diff !hyper diffusivity
903  REAL (R8) :: hyper_diff_exp !explicit hyper diffusivity [m^2/s]
904  REAL (R8) :: hyper_diff_imp !implicit hyper diffusivity
905  REAL (R8) :: diff_hyper !additional diffusivity to remove profile oscillations
906 !dy limit diffusion coefficients
907  REAL (R8), parameter :: nine_diff_limit=1.0e4
908 !dy
909  integer :: irho_err
910 
911 
912  hyper_diff_exp = hyper_diff(1) !AF 22.Aug.2016
913  hyper_diff_imp = hyper_diff(2) !AF 22.Aug.2016
914 
915 ! +++ Set up dimensions:
916  ndim = 1 !no coupling between density equations
917  nrho = profiles%NRHO
918  nion = profiles%NION
919  nmodel = 3
920 
921 !AF, 11.Oct.2011
922 ALLOCATE(local_int_source_s4(nion))
923 local_int_source_s4(:) = 0.0_r8
924 !AF - End
925 
926 ! +++ Allocate types for interface with numerical solver:
927  CALL allocate_numerics(ndim, nrho, solver, ifail)
928 
929 
930 
931 ! +++ Set up local variables:
932  amix = control%AMIX
933  tau = control%TAU
934  solver_type = control%SOLVER_TYPE
935 
936  bt = geometry%BGEO
937  btm = evolution%BTM
938  btprime = (bt-btm)/tau
939 
940 
941  rho_loop1: DO irho=1,nrho
942  rho(irho) = geometry%RHO(irho)
943  vpr(irho) = geometry%VPR(irho)
944  vprm(irho) = evolution%VPRM(irho)
945  g1(irho) = geometry%G1(irho)
946  END DO rho_loop1
947 
948 
949 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
950 ! solution of particle transport equation for
951 ! individual ion species
952 
953  ion_loop1: DO iion=1,nion
954 
955 ! +++ Set equation to 'predictive' and all coefficients to zero:
956  flag = 1
957  y(:) = 0.0e0_r8
958  dy(:) = 0.0e0_r8
959  ym(:) = 0.0e0_r8
960  dym(:) = 0.0e0_r8 !AF - 25.Sep.2014
961  a(:) = 0.0e0_r8
962  b(:) = 0.0e0_r8
963  c(:) = 0.0e0_r8
964  d(:) = 0.0e0_r8
965  e(:) = 0.0e0_r8
966  f(:) = 0.0e0_r8
967  g(:) = 0.0e0_r8
968  h = 0.0e0_r8
969  v(:) = 0.0e0_r8
970  u(:) = 0.0e0_r8
971  w(:) = 0.0e0_r8
972 
973 
974 
975 ! +++ Set up boundary conditions for particular ion type:
976  ni_bnd_type(2) = profiles%NI_BND_TYPE(iion)
977  ni_bnd(2,1) = profiles%NI_BND(1,iion)
978  ni_bnd(2,2) = profiles%NI_BND(2,iion)
979  ni_bnd(2,3) = profiles%NI_BND(3,iion)
980 
981  rho_bnd = profiles%NI_BND_RHO(iion)
982 
983 
984 ! +++ Set up local variables for particular ion type:
985  rho_loop2: DO irho = 1,nrho
986  ni(irho) = profiles%NI(irho,iion)
987  dni(irho) = profiles%DNI(irho,iion) !AF - 25.Sep.2014
988  nim(irho) = evolution%NIM(irho,iion)
989  dnim(irho) = evolution%DNIM(irho,iion) !AF - 25.Sep.2014
990 
991  diff(irho) = 0.e0_r8
992  vconv(irho) = 0.e0_r8
993 
994  model_loop1: DO imodel = 1, nmodel
995  c1(imodel) = transport%C1(imodel)
996 
997  diff_mod(irho,imodel) = transport%DIFF_NI(irho,iion,imodel)
998  vconv_mod(irho,imodel) = transport%VCONV_NI(irho,iion,imodel)
999 
1000  diff(irho) = diff(irho) + diff_mod(irho,imodel)
1001  vconv(irho) = vconv(irho) + vconv_mod(irho,imodel)
1002  END DO model_loop1
1003 
1004  si_exp(irho) = 0.e0_r8
1005  si_imp(irho) = 0.e0_r8
1006 
1007  si_exp(irho) = si_exp(irho) + sources%SI_EXP(irho,iion)
1008  si_imp(irho) = si_imp(irho) + sources%SI_IMP(irho,iion)
1009 
1010  END DO rho_loop2
1011 
1012 
1013 
1014 ! +++ Coefficients for ion diffusion equation in form:
1015 !
1016 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
1017 
1018  diff_hyper = hyper_diff_exp + hyper_diff_imp*maxval(diff) !AF 25.Apr.2016, 22.Aug.2016
1019 
1020  rho_loop3: DO irho=1,nrho
1021  y(irho) = ni(irho)
1022  dy(irho) = dni(irho) !AF - 25.Sep.2014
1023  ym(irho) = nim(irho)
1024  dym(irho) = dnim(irho) !AF - 25.Sep.2014
1025 
1026  a(irho) = vpr(irho)
1027  b(irho) = vprm(irho)
1028  c(irho) = 1.e0_r8
1029 ! D(IRHO) = VPR(IRHO)*G1(IRHO)*DIFF(IRHO)
1030  d(irho) = vpr(irho)*g1(irho)*(diff(irho)+diff_hyper) !AF 25.Apr.2016
1031 ! E(IRHO) = VPR(IRHO)*G1(IRHO)*VCONV(IRHO) &
1032  e(irho) = vpr(irho)*g1(irho)*(vconv(irho)+diff_hyper*dnim(irho)/nim(irho)) & !AF 25.Apr.2016
1033  - btprime/2.e0_r8/bt*rho(irho)*vpr(irho)
1034  f(irho) = vpr(irho)*si_exp(irho)
1035  g(irho) = vpr(irho)*si_imp(irho)
1036  END DO rho_loop3
1037  h = tau
1038 
1039 
1040 
1041 ! +++ Boundary conditions for ion diffusion equation in form:
1042 !
1043 ! V*Y' + U*Y =W
1044 !
1045 ! +++ On axis:
1046 ! dNi/drho(rho=0)=0:
1047  IF (solver_type.NE.4) THEN !AF 11.Oct.2011
1048  v(1) = 1.e0_r8
1049  u(1) = 0.e0_r8
1050  ELSE !AF 11.Oct.2011 - Zero flux instead of zero gradient at the axis for solver 4
1051 ! IF (DIFF(1).GT.1.0E-6) THEN !AF 19.Mar.2012 - To avoid problems with the axis boundary condition
1052  IF ((diff(1)+diff_hyper).GT.1.0e-6) THEN !AF 19.Mar.2012 - To avoid problems with the axis boundary condition !AF 25.Apr.2016
1053 ! V(1) = -DIFF(1)
1054  v(1) = -diff(1)-diff_hyper !AF 25.Apr.2016
1055  ELSE
1056  v(1) = -1.0e-6
1057  ENDIF
1058 ! U(1) = VCONV(1) !AF 25.Apr.2016
1059  u(1) = vconv(1)+diff_hyper*dnim(1)/nim(1) !AF 25.Apr.2016
1060  END IF !AF 11.Oct.2011
1061  w(1) = 0.e0_r8
1062 
1063 ! +++ At the edge:
1064 ! FIXED Ni
1065  IF(ni_bnd_type(2).EQ.1) THEN
1066  v(2) = 0.e0_r8
1067  u(2) = 1.e0_r8
1068  w(2) = ni_bnd(2,1)
1069  ENDIF
1070 
1071 ! FIXED grad_Ni
1072  IF(ni_bnd_type(2).EQ.2) THEN
1073  v(2) = 1.e0_r8
1074  u(2) = 0.e0_r8
1075  w(2) = -ni_bnd(2,1)
1076  ENDIF
1077 
1078 ! FIXED L_Ni
1079  IF(ni_bnd_type(2).EQ.3) THEN
1080  v(2) = ni_bnd(2,1)
1081  u(2) = 1.e0_r8
1082  w(2) = 0.e0_r8
1083  ENDIF
1084 
1085 ! FIXED Flux_Ni
1086  IF(ni_bnd_type(2).EQ.4) THEN
1087 ! V(2) = -G(NRHO)*DIFF(NRHO)
1088 ! U(2) = G(NRHO)*VCONV(NRHO)
1089  v(2) = -vpr(nrho)*g1(nrho)*diff(nrho)
1090  u(2) = vpr(nrho)*g1(nrho)*vconv(nrho)
1091  w(2) = ni_bnd(2,1)
1092  ENDIF
1093 
1094 ! Generic boundary condition
1095  IF(ni_bnd_type(2).EQ.5) THEN
1096  v(2) = ni_bnd(2,1)
1097  u(2) = ni_bnd(2,2)
1098  w(2) = ni_bnd(2,3)
1099  ENDIF
1100 
1101 
1102 
1103 ! +++ Density equation is not solved:
1104  IF(ni_bnd_type(2).EQ.0) THEN
1105 
1106  CALL derivn(nrho,rho,y,dy)
1107 
1108  flag = 0
1109 
1110  rho_loop4: DO irho=1,nrho
1111  a(irho) = 1.0e0_r8
1112  b(irho) = 1.0e0_r8
1113  c(irho) = 1.0e0_r8
1114  d(irho) = 0.0e0_r8
1115  e(irho) = 0.0e0_r8
1116  f(irho) = 0.0e0_r8
1117  g(irho) = 0.0e0_r8
1118  END DO rho_loop4
1119 
1120  v(2) = 0.0e0_r8
1121  u(2) = 1.0e0_r8
1122  w(2) = y(nrho)
1123  END IF
1124 
1125 
1126 
1127 ! +++ Defining coefficients for numerical solver:
1128  solver%TYPE = solver_type
1129  solver%EQ_FLAG(ndim) = flag
1130  solver%NDIM = ndim
1131  solver%NRHO = nrho
1132  solver%AMIX = amix
1133  solver%DERIVATIVE_FLAG(1) = 0
1134 
1135  rho_loop5: DO irho=1,nrho
1136 
1137  solver%RHO(irho) = rho(irho)
1138 
1139  solver%Y(ndim,irho) = y(irho)
1140  solver%DY(ndim,irho) = dy(irho)
1141  solver%YM(ndim,irho) = ym(irho)
1142 
1143  solver%A(ndim,irho) = a(irho)
1144  solver%B(ndim,irho) = b(irho)
1145  solver%C(ndim,irho) = c(irho)
1146  solver%D(ndim,irho) = d(irho)
1147  solver%E(ndim,irho) = e(irho)
1148  solver%F(ndim,irho) = f(irho)
1149  solver%G(ndim,irho) = g(irho)
1150 
1151  END DO rho_loop5
1152 
1153  solver%H = h
1154 
1155  solver%V(ndim,1) = v(1)
1156  solver%U(ndim,1) = u(1)
1157  solver%W(ndim,1) = w(1)
1158  solver%V(ndim,2) = v(2)
1159  solver%U(ndim,2) = u(2)
1160  solver%W(ndim,2) = w(2)
1161 
1162 
1163 
1164 ! +++ Solution of density diffusion equation:
1165  CALL solution_interface(solver, ifail)
1166 
1167  ! dy check for nans in solution, return if true)
1168  IF (any( isnan(solver%y))) then
1169  failstring='Error in the ion density equation, nans in the solution, stop'
1170  ifail=-1
1171  return
1172  end if
1173 
1174 ! +++ New ion density:
1175  rho_loop6: DO irho=1,nrho
1176  y(irho) = solver%Y(ndim,irho)
1177  dy(irho) = solver%DY(ndim,irho)
1178  END DO rho_loop6
1179 
1180  IF(ni_bnd_type(2).EQ.0) THEN
1181  y(:) = profiles%NI(:,iion)
1182  CALL derivn(nrho,rho,y,dy)
1183  END IF
1184 
1185 
1186 ! +++ New profiles of ion density flux and integral source:
1187  rho_loop7: DO irho=1,nrho
1188 ! NI(IRHO) = Y(IRHO)
1189 ! DNI(IRHO) = DY(IRHO)
1190 ! IF (RHO(IRHO).NE.0.E0_R8) THEN
1191 ! FUN1(IRHO) = 1.e0_R8/RHO(IRHO)*(VPR(IRHO)*SI_EXP(IRHO) & !AF 11.Oct.2011 - the division by rho was needed because the routine INTEGR below actually integrates the function multiplied by rho
1192 ! +VPRM(IRHO)*NIM(IRHO)/TAU &
1193 ! -NI(IRHO)*VPR(IRHO)*(1.e0_R8/TAU+SI_IMP(IRHO)))
1194 ! ELSE
1195 ! FUN1(IRHO) = 1.e0_R8*(SI_EXP(IRHO)+NIM(IRHO)/TAU-NI(IRHO) & !AF 11.Oct.2011 - this is only OK with solver_test since Vprime.eq.rho in that case - get it fixed in the trunk!
1196 ! *(1.e0_R8/TAU+SI_IMP(IRHO)))
1197 ! ENDIF
1198 
1199 
1200  ni(irho) = y(irho)
1201  dni(irho) = dy(irho)
1202  if (ni(irho).le.0.0) then
1203  if (irho.eq.1) then
1204  failstring='Error in the density equation: on-axis ion density is negative, stop'
1205  write(*,*) 'value of on-axis density for ion', iion, 'is', ni(1)
1206  write(*,*) 'check ni profile in the file err_tmp'
1207  open (22,file='err_tmp')
1208  do irho_err=1,nrho
1209  write(22,*) irho_err, y(irho_err)
1210  enddo
1211  close(22)
1212  ifail=-1
1213  return
1214  else
1215  write(*,*) 'warning, density for ion ',iion,' and irho ',irho,'is negative, set it to the value at the previous irho'
1216  ni(irho)=ni(irho-1)
1217  dni(irho)=dni(irho-1)
1218  end if
1219  end if
1220  fun1(irho) = (vpr(irho)*si_exp(irho) & !AF 11.Oct.2011 - the division by rho was needed because the routine INTEGR below actually integrates the function multiplied by rho
1221  +vprm(irho)*nim(irho)/tau &
1222  -ni(irho)*vpr(irho)*(1.e0_r8/tau+si_imp(irho)))
1223  END DO rho_loop7
1224 
1225 ! CALL INTEGR(NRHO,RHO,FUN1,INTFUN1) !Integral source !AF 11.Oct.2011
1226  CALL integr2(nrho,rho,fun1,intfun1) !Integral source !AF 11.Oct.2011 - this routine simply integrates the function, not the function times rho as INTEGR was doing
1227 
1228 
1229 !AF 12.Oct.2011 -
1230  local_fun1_s4 = (si_exp(1) & !stripping Vprime by assuming that VPR and VPRM are approximately the same, which they are if, as usually, rho(1).eq.0
1231  +nim(1)/tau &
1232  -ni(1)*(1.e0_r8/tau+si_imp(1)))/g1(1) !stripping G1
1233  local_intfun1_s4 = local_fun1_s4*rho(1)/2.e0_r8 !stripped integral for the axis - from the INTEGR2 routine, first point only
1234 !AF - End
1235 
1236 
1237  rho_loop8: DO irho=1,nrho
1238  int_source(irho) = intfun1(irho) &
1239  + btprime/2.e0_r8/bt*rho(irho)*vpr(irho)*ni(irho)
1240  flux(irho) = vpr(irho)*g1(irho)* &
1241  ( y(irho)*vconv(irho) - dy(irho)*diff(irho) )
1242 
1243 
1244 
1245 ! +++ Contribution to ion energy transport:
1246  flux_ni_conv(irho) = 0.e0_r8
1247 
1248  model_loop2: DO imodel = 1, nmodel
1249  flux_ni_conv(irho) = flux_ni_conv(irho) &
1250  + c1(imodel)*vpr(irho)*g1(irho) &
1251  *( y(irho)*vconv_mod(irho,imodel) &
1252  -dy(irho)*diff_mod(irho,imodel) )
1253  END DO model_loop2
1254 
1255 
1256 
1257 ! +++ If equation is not solved, flux is determined
1258 ! by the integral of sources and transport coefficients
1259 ! are updated with effective values:
1260  IF (ni_bnd_type(2).EQ.0) THEN
1261  diff(irho) = 1.d-6
1262  flux(irho) = int_source(irho)
1263  flux_ni_conv(irho) = 1.5e0_r8*int_source(irho)
1264  IF ((vpr(irho)*g1(irho).NE.0.0_r8).and.(dy(irho).ne.0.0_r8)) &
1265  diff(irho) = - flux(irho) / dy(irho) / (vpr(irho)*g1(irho))
1266  if (abs(diff(irho)).ge.nine_diff_limit) &
1267  diff(irho)=sign(nine_diff_limit,diff(irho))
1268  vconv(irho) = 0.0_r8
1269  IF (diff(irho).LE.1.d-6 .AND. vpr(irho)*g1(irho) .NE. 0 ) THEN
1270  diff(irho) = 1.d-6
1271  vconv(irho) = (flux(irho) / (max(abs(vpr(irho)), 1.e-6)*g1(irho)) + dy(irho)*diff(irho)) /max(abs(y(irho)), 1.e-6)
1272  END IF
1273  END IF
1274 
1275 
1276 
1277 ! +++ Return new ion density and flux profiles to the work flow:
1278  profiles%NI(irho,iion) = ni(irho)
1279  profiles%DNI(irho,iion) = dni(irho) !AF - 25.Sep.2014
1280  profiles%DIFF_NI(irho,iion) = diff(irho)
1281  profiles%VCONV_NI(irho,iion) = vconv(irho)
1282  profiles%FLUX_NI(irho,iion) = flux(irho)
1283  profiles%FLUX_NI_CONV(irho,iion) = flux_ni_conv(irho)
1284 !dy 2017-10-06 PROFILES%SOURCE_NI(IRHO,IION) = SI_EXP(IRHO) + SI_IMP(IRHO) * NI(IRHO)
1285  profiles%SOURCE_NI(irho,iion) = si_exp(irho) - si_imp(irho) * ni(irho)
1286  profiles%INT_SOURCE_NI(irho,iion) = int_source(irho)
1287 
1288 
1289  END DO rho_loop8
1290 
1291  fun1=profiles%SOURCE_NI(:,iion)*vpr
1292  CALL integr2(nrho,rho,fun1,intfun1)
1293  profiles%INT_SOURCE_NI(:,iion) = intfun1
1294 
1295  !AF 11.Oct.2011 - local flux information for the axis boundary conditions in temperature and rotation equations
1296 
1297  local_int_source_s4(iion) = local_intfun1_s4 & !AF 11.Oct.2011 - LOCAL_INT_SOURCE_S4 is also stripped from Vprime and G1 above, with a reasonable approximation in the case of Vprime
1298  + btprime/2.e0_r8/bt*rho(1)*ni(1)/g1(1)
1299  local_flux_ni_s4(iion) = ( y(1)*vconv(1) - dy(1)*diff(1) )
1300 
1301  local_flux_ni_conv_s4(iion) = 0.e0_r8
1302  model_loop3: DO imodel = 1, nmodel
1303  local_flux_ni_conv_s4(iion) = local_flux_ni_conv_s4(iion) &
1304  + c1(imodel) &
1305  *( y(1)*vconv_mod(1,imodel) &
1306  -dy(1)*diff_mod(1,imodel) )
1307  END DO model_loop3
1308 
1309  ! +++ If equation is not solved, flux is determined
1310  ! by the integral of sources:
1311  ! IF (NI_BND_TYPE(2).EQ.0) THEN
1312  IF (1.EQ.0) THEN ! AF - 7.Jul.2010, same post-treatment as in predictive mode, but using the analytical solution instead of the numerical one
1313  local_flux_ni_s4(iion) = local_int_source_s4(iion)
1314  local_flux_ni_conv_s4(iion) = 1.5e0_r8*local_int_source_s4(iion)
1315  END IF
1316  !AF - End
1317 
1318  END DO ion_loop1
1319 
1320 
1321 
1322 ! +++ Deallocate types for interface with numerical solver:
1323  CALL deallocate_numerics(solver, ifail)
1324 
1325 !AF, 11.Oct.2011
1326 DEALLOCATE(local_int_source_s4)
1327 !AF - End
1328 
1329  RETURN
1330 
1331 
1332 
1333 END SUBROUTINE ion_density
1334 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1335 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1336 
1337 
1338 
1339 
1340 
1341 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1342 SUBROUTINE electron_density &
1343  (geometry,profiles,transport,sources,evolution,control,hyper_diff,ifail,failstring) !AF 25.Apr.2016, 22.Aug.2016
1344 
1345 !-------------------------------------------------------!
1346 ! This subroutine solves electron transport !
1347 ! equation and provides: !
1348 ! density and flux of electrons. !
1349 !-------------------------------------------------------!
1350 
1351  USE ets_plasma
1352  USE type_solver
1353 
1354  USE s4_parameters
1355 
1356  IMPLICIT NONE
1357 
1358  INTEGER, INTENT (INOUT) :: ifail
1359  CHARACTER(LEN=500) :: failstring
1360 
1361 ! +++ External parameters:
1362  TYPE (magnetic_geometry) :: geometry !contains all geometry quantities
1363  TYPE (plasma_profiles) :: profiles !contains profiles of plasma parameters
1364  TYPE (transport_coefficients) :: transport !contains profiles of trasport coefficients
1365  TYPE (sources_and_sinks) :: sources !contains profiles of sources
1366  TYPE (time_evolution) :: evolution !contains all parameters required by time evolution
1367  TYPE (run_control) :: control !contains all parameters required by run
1368 
1369 !evolution and convergence loop
1370 ! +++ Input/Output to numerical solver:
1371  TYPE (numerics) :: solver !contains all I/O quantities to numerics part
1372 
1373 ! +++ Internal parameters:
1374  INTEGER :: irho, nrho !radius index, number of radial points
1375  INTEGER :: imodel, nmodel !model index, number of transport modules used
1376 
1377  REAL (R8) :: bt, btm, btprime !magnetic field from current time step, [T], previous time steps, [T], time derivative, [T/s]
1378  REAL (R8) :: rho(profiles%nrho) !normalised minor radius, [m]
1379  REAL (R8) :: vpr(profiles%nrho) !V', [m^2]
1380  REAL (R8) :: vprm(profiles%nrho) !V' (at previous time step), [m^2]
1381  REAL (R8) :: g1(profiles%nrho)
1382 
1383  REAL (R8) :: ne(profiles%nrho),nem(profiles%nrho) !density from current ans previous time step, [m^-3]
1384  REAL (R8) :: dne(profiles%nrho) !density gradient, [m^-4]
1385  REAL (R8) :: dnem(profiles%nrho) !density gradient at previous time step, [m^-4] !AF - 25.Sep.2014
1386  REAL (R8) :: flux(profiles%nrho) !ion flux, [1/s]
1387  REAL (R8) :: int_source(profiles%nrho) !integral of source [1/s]
1388  REAL (R8) :: flux_ne_conv(profiles%nrho) !ion flux, contributing to convective heat transport [1/s]
1389 
1390  REAL (R8) :: diff(profiles%nrho),vconv(profiles%nrho) !diffusion coefficient and pinch velocity, [m^2/s] and [m/s]
1391  REAL (R8) :: diff_mod(transport%nrho,3) !diffusion coefficient due to particular transport model [m^2/s]
1392  REAL (R8) :: vconv_mod(transport%nrho,3) !pinch velocity due to particular transport model [m/s]
1393  REAL (R8) :: c1(3) !coefficient in front of convective term in energy balance equation
1394 
1395  REAL (R8) :: se_exp(profiles%nrho) !explicit particle source, [1/(m^3*s)]
1396  REAL (R8) :: se_imp(profiles%nrho) !implicit particle source (proportional to density), [1/s]
1397 
1398  REAL (R8) :: ne_bnd(2,3) !boundary condition, value, [depends on NI_BND_TYPE]
1399  INTEGER :: ne_bnd_type(2) !boundary condition, type
1400 
1401  REAL (R8) :: amix, tau !mixing factor, time step, [s]
1402 
1403  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
1404  INTEGER :: ndim !number of equations to be solved
1405  INTEGER :: solver_type !specifies the option for numerical solution
1406  REAL (R8) :: y(profiles%nrho) !function at the current amd previous time steps
1407  REAL (R8) :: ym(profiles%nrho) !function at the current amd previous time steps
1408  REAL (R8) :: dy(profiles%nrho) !derivative of function
1409  REAL (R8) :: dym(profiles%nrho) !derivative of function at previous time step !AF - 25.Sep.2014
1410  REAL (R8) :: a(profiles%nrho) !coefficients for numerical solver
1411  REAL (R8) :: b(profiles%nrho) !coefficients for numerical solver
1412  REAL (R8) :: c(profiles%nrho) !coefficients for numerical solver
1413  REAL (R8) :: d(profiles%nrho) !coefficients for numerical solver
1414  REAL (R8) :: e(profiles%nrho) !coefficients for numerical solver
1415  REAL (R8) :: f(profiles%nrho) !coefficients for numerical solver
1416  REAL (R8) :: g(profiles%nrho) !coefficients for numerical solver
1417  REAL (R8) :: h !coefficients for numerical solver
1418  REAL (R8) :: v(2), u(2), w(2) !boundary conditions for numerical solver
1419 
1420  REAL (R8) :: fun1(profiles%nrho), intfun1(profiles%nrho)
1421 
1422  REAL (R8) :: local_int_source_s4 !auxiliary variable with the integral of sources divided by the metrics coefficients Vprime and G1
1423  REAL (R8) :: local_fun1_s4, local_intfun1_s4
1424 
1425 ! +++ Stabilization scheme !AF 25.Apr.2016, 22.Aug.2016
1426  REAL (R8), DIMENSION(2) :: hyper_diff !hyper diffusivity
1427  REAL (R8) :: hyper_diff_exp !explicit hyper diffusivity [m^2/s]
1428  REAL (R8) :: hyper_diff_imp !implicit hyper diffusivity
1429  REAL (R8) :: diff_hyper !additional diffusivity to remove profile oscillations
1430 !dy limit diffusion coefficients
1431  REAL (R8), parameter :: nine_diff_limit=1.0e4
1432 !dy
1433  integer :: irho_err
1434 
1435 
1436  hyper_diff_exp = hyper_diff(1) !AF 22.Aug.2016
1437  hyper_diff_imp = hyper_diff(2) !AF 22.Aug.2016
1438 
1439 ! +++ Set up dimensions:
1440  ndim = 1 !no coupling between density equations
1441  nrho = profiles%NRHO
1442  nmodel = 3
1443 
1444 
1445 ! +++ Allocate types for interface with numerical solver:
1446  CALL allocate_numerics(ndim, nrho, solver, ifail)
1447 
1448 
1449 
1450 ! +++ Set up local variables:
1451  amix = control%AMIX
1452  tau = control%TAU
1453  solver_type = control%SOLVER_TYPE
1454 
1455  bt = geometry%BGEO
1456  btm = evolution%BTM
1457  btprime = (bt-btm)/tau
1458 
1459 
1460  rho_loop1: DO irho=1,nrho
1461  rho(irho) = geometry%RHO(irho)
1462  vpr(irho) = geometry%VPR(irho)
1463  vprm(irho) = evolution%VPRM(irho)
1464  g1(irho) = geometry%G1(irho)
1465  END DO rho_loop1
1466 
1467 
1468 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1469 ! solution of particle transport equation
1470 
1471 ! +++ Set equation to 'predictive' and all coefficients to zero:
1472  flag = 1
1473  y(:) = 0.0e0_r8
1474  dy(:) = 0.0e0_r8
1475  ym(:) = 0.0e0_r8
1476  dym(:) = 0.0e0_r8 !AF - 25.Sep.2014
1477  a(:) = 0.0e0_r8
1478  b(:) = 0.0e0_r8
1479  c(:) = 0.0e0_r8
1480  d(:) = 0.0e0_r8
1481  e(:) = 0.0e0_r8
1482  f(:) = 0.0e0_r8
1483  g(:) = 0.0e0_r8
1484  h = 0.0e0_r8
1485  v(:) = 0.0e0_r8
1486  u(:) = 0.0e0_r8
1487  w(:) = 0.0e0_r8
1488 
1489 
1490 
1491 ! +++ Set up boundary conditions for particular ion type:
1492  ne_bnd_type(2) = profiles%NE_BND_TYPE
1493  ne_bnd(2,1) = profiles%NE_BND(1)
1494  ne_bnd(2,2) = profiles%NE_BND(2)
1495  ne_bnd(2,3) = profiles%NE_BND(3)
1496 
1497 
1498 
1499 ! +++ Set up local variables for particular ion type:
1500  rho_loop2: DO irho = 1,nrho
1501  ne(irho) = profiles%NE(irho)
1502  dne(irho) = profiles%DNE(irho) !AF - 25.Sep.2014
1503  nem(irho) = evolution%NEM(irho)
1504  dnem(irho) = evolution%DNEM(irho) !AF - 25.Sep.2014
1505 
1506  diff(irho) = 0.e0_r8
1507  vconv(irho) = 0.e0_r8
1508 
1509  model_loop1: DO imodel = 1, nmodel
1510  c1(imodel) = transport%C1(imodel)
1511 
1512  diff_mod(irho,imodel) = transport%DIFF_NE(irho,imodel)
1513  vconv_mod(irho,imodel) = transport%VCONV_NE(irho,imodel)
1514 
1515  diff(irho) = diff(irho) + diff_mod(irho,imodel)
1516  vconv(irho) = vconv(irho) + vconv_mod(irho,imodel)
1517  END DO model_loop1
1518 
1519  se_exp(irho) = sources%SE_EXP(irho)
1520  se_imp(irho) = sources%SE_IMP(irho)
1521 
1522  END DO rho_loop2
1523 
1524 
1525 
1526 ! +++ Coefficients for electron diffusion equation in form:
1527 !
1528 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
1529 
1530  diff_hyper = hyper_diff_exp + hyper_diff_imp*maxval(diff) !AF 25.Apr.2016, 22.Aug.2016
1531 
1532  rho_loop3: DO irho=1,nrho
1533  y(irho) = ne(irho)
1534  dy(irho) = dne(irho) !AF - 25.Sep.2014
1535  ym(irho) = nem(irho)
1536  dym(irho) = dnem(irho) !AF - 25.Sep.2014
1537 
1538  a(irho) = vpr(irho)
1539  b(irho) = vprm(irho)
1540  c(irho) = 1.e0_r8
1541 ! D(IRHO) = VPR(IRHO)*G1(IRHO)*DIFF(IRHO)
1542  d(irho) = vpr(irho)*g1(irho)*(diff(irho)+diff_hyper) !AF 25.Apr.2016
1543 ! E(IRHO) = VPR(IRHO)*G1(IRHO)*VCONV(IRHO) &
1544  e(irho) = vpr(irho)*g1(irho)*(vconv(irho)+diff_hyper*dnem(irho)/nem(irho)) & !AF 25.Apr.2016
1545  - btprime/2.e0_r8/bt*rho(irho)*vpr(irho)
1546  f(irho) = vpr(irho)*se_exp(irho)
1547  g(irho) = vpr(irho)*se_imp(irho)
1548  END DO rho_loop3
1549  h = tau
1550 
1551 
1552 
1553 ! +++ Boundary conditions for electron diffusion equation in form:
1554 !
1555 ! V*Y' + U*Y =W
1556 !
1557 ! +++ On axis:
1558 ! dNi/drho(rho=0)=0: !AF 25.Apr.2016 - this is Ne, not Ni
1559  IF (solver_type.NE.4) THEN
1560  v(1) = 1.e0_r8
1561  u(1) = 0.e0_r8
1562  ELSE
1563 ! IF (DIFF(1).GT.1.0E-6) THEN
1564  IF ((diff(1)+diff_hyper).GT.1.0e-6) THEN !AF 25.Apr.2016
1565 ! V(1) = -DIFF(1)
1566  v(1) = -diff(1)-diff_hyper !AF 25.Apr.2016
1567  ELSE
1568  v(1) = -1.0e-6
1569  ENDIF
1570 ! U(1) = VCONV(1)
1571  u(1) = vconv(1)+diff_hyper*dnem(1)/nem(1) !AF 25.Apr.2016
1572  END IF
1573  w(1) = 0.e0_r8
1574 
1575 ! +++ At the edge:
1576 ! FIXED Ne
1577  IF(ne_bnd_type(2).EQ.1) THEN
1578  v(2) = 0.e0_r8
1579  u(2) = 1.e0_r8
1580  w(2) = ne_bnd(2,1)
1581  ENDIF
1582 
1583 ! FIXED grad_Ne
1584  IF(ne_bnd_type(2).EQ.2) THEN
1585  v(2) = 1.e0_r8
1586  u(2) = 0.e0_r8
1587  w(2) = -ne_bnd(2,1)
1588  ENDIF
1589 
1590 ! FIXED L_Ne
1591  IF(ne_bnd_type(2).EQ.3) THEN
1592  v(2) = ne_bnd(2,1)
1593  u(2) = 1.e0_r8
1594  w(2) = 0.e0_r8
1595  ENDIF
1596 
1597 ! FIXED Flux_Ne
1598  IF(ne_bnd_type(2).EQ.4) THEN
1599  v(2) = -vpr(nrho)*g1(nrho)*diff(nrho)
1600  u(2) = vpr(nrho)*g1(nrho)*vconv(nrho)
1601  w(2) = ne_bnd(2,1)
1602  ENDIF
1603 
1604 ! Generic boundary condition
1605  IF(ne_bnd_type(2).EQ.5) THEN
1606  v(2) = ne_bnd(2,1)
1607  u(2) = ne_bnd(2,2)
1608  w(2) = ne_bnd(2,3)
1609  ENDIF
1610 
1611 
1612 
1613 ! +++ Density equation is not solved:
1614  IF(ne_bnd_type(2).EQ.0) THEN
1615 
1616  CALL derivn(nrho,rho,y,dy)
1617 
1618  flag = 0
1619 
1620  rho_loop4: DO irho=1,nrho
1621  a(irho) = 1.0e0_r8
1622  b(irho) = 1.0e0_r8
1623  c(irho) = 1.0e0_r8
1624  d(irho) = 0.0e0_r8
1625  e(irho) = 0.0e0_r8
1626  f(irho) = 0.0e0_r8
1627  g(irho) = 0.0e0_r8
1628  END DO rho_loop4
1629 
1630  v(2) = 0.0e0_r8
1631  u(2) = 1.0e0_r8
1632  w(2) = y(nrho)
1633  END IF
1634 
1635 
1636 
1637 ! +++ Defining coefficients for numerical solver:
1638  solver%TYPE = solver_type
1639  solver%EQ_FLAG(ndim) = flag
1640  solver%NDIM = ndim
1641  solver%NRHO = nrho
1642  solver%AMIX = amix
1643  solver%DERIVATIVE_FLAG(1) = 0
1644 
1645  rho_loop5: DO irho=1,nrho
1646 
1647  solver%RHO(irho) = rho(irho)
1648 
1649  solver%Y(ndim,irho) = y(irho)
1650  solver%DY(ndim,irho) = dy(irho)
1651  solver%YM(ndim,irho) = ym(irho)
1652 
1653  solver%A(ndim,irho) = a(irho)
1654  solver%B(ndim,irho) = b(irho)
1655  solver%C(ndim,irho) = c(irho)
1656  solver%D(ndim,irho) = d(irho)
1657  solver%E(ndim,irho) = e(irho)
1658  solver%F(ndim,irho) = f(irho)
1659  solver%G(ndim,irho) = g(irho)
1660 
1661  END DO rho_loop5
1662 
1663  solver%H = h
1664 
1665  solver%V(ndim,1) = v(1)
1666  solver%U(ndim,1) = u(1)
1667  solver%W(ndim,1) = w(1)
1668  solver%V(ndim,2) = v(2)
1669  solver%U(ndim,2) = u(2)
1670  solver%W(ndim,2) = w(2)
1671 
1672 
1673 
1674 ! +++ Solution of density diffusion equation:
1675  CALL solution_interface(solver, ifail)
1676 
1677  ! dy check for nans in solution, return if true)
1678  IF (any( isnan(solver%y))) then
1679  failstring='Error in the electron density equation: nans in the solution, stop'
1680  ifail=-1
1681  return
1682  end if
1683 
1684 ! +++ New electron density:
1685  rho_loop6: DO irho=1,nrho
1686  y(irho) = solver%Y(ndim,irho)
1687  dy(irho) = solver%DY(ndim,irho)
1688  END DO rho_loop6
1689 
1690  IF(ne_bnd_type(2).EQ.0) THEN
1691  y(:) = profiles%NE(:)
1692  CALL derivn(nrho,rho,y,dy)
1693  END IF
1694 
1695 
1696 ! +++ New profiles of electron density flux and integral source:
1697  rho_loop7: DO irho=1,nrho
1698 
1699  ne(irho) = y(irho)
1700  dne(irho) = dy(irho)
1701  if (ne(irho).le.0.0) then
1702  if (irho.eq.1) then
1703  failstring='Error in the density equation: on-axis electron density is negative, stop'
1704  write(*,*) 'value of on-axis electron density is', ne(1)
1705  write(*,*) 'check ne profile in the file err_tmp'
1706  open (22,file='err_tmp')
1707  do irho_err=1,nrho
1708  write(22,*) irho_err, y(irho_err)
1709  enddo
1710  close(22)
1711  ifail=-1
1712  return
1713  else
1714  write(*,*) 'warning, electron density for ', irho,' is negative, set it to the value at the previous irho'
1715  ne(irho)=ne(irho-1)
1716  dne(irho)=dne(irho-1)
1717  end if
1718  end if
1719  fun1(irho) = (vpr(irho)*se_exp(irho) + vprm(irho)*nem(irho)/tau &
1720  -ne(irho)*vpr(irho)*(1.e0_r8/tau+se_imp(irho)))
1721  END DO rho_loop7
1722 
1723  CALL integr2(nrho,rho,fun1,intfun1)
1724 
1725 
1726  local_fun1_s4 = (se_exp(1) + nem(1)/tau - ne(1)*(1.e0_r8/tau+se_imp(1)))/g1(1)
1727  local_intfun1_s4 = local_fun1_s4*rho(1)/2.e0_r8
1728 
1729 
1730  rho_loop8: DO irho=1,nrho
1731  int_source(irho) = intfun1(irho) + btprime/2.e0_r8/bt*rho(irho)*vpr(irho)*ne(irho)
1732  flux(irho) = vpr(irho)*g1(irho)* ( y(irho)*vconv(irho) - dy(irho)*diff(irho) )
1733 
1734 
1735 
1736 ! +++ Contribution to electron energy transport:
1737  flux_ne_conv(irho) = 0.e0_r8
1738 
1739  model_loop2: DO imodel = 1, nmodel
1740  flux_ne_conv(irho) = flux_ne_conv(irho) + c1(imodel)*vpr(irho)*g1(irho) &
1741  *( y(irho)*vconv_mod(irho,imodel) &
1742  -dy(irho)*diff_mod(irho,imodel) )
1743  END DO model_loop2
1744 
1745 
1746 
1747 ! +++ If equation is not solved, flux is determined
1748 ! by the integral of sources:
1749  IF (ne_bnd_type(2).EQ.0) THEN
1750  diff(irho) = 1.d-6
1751  flux(irho) = int_source(irho)
1752  flux_ne_conv(irho) = 1.5e0_r8*int_source(irho)
1753  IF (vpr(irho)*g1(irho).NE.0.0_r8) &
1754  diff(irho) = - flux(irho) / dy(irho) / (vpr(irho)*g1(irho))
1755  if (abs(diff(irho)).ge.nine_diff_limit) &
1756  diff(irho)=sign(nine_diff_limit,diff(irho))
1757  vconv(irho) = 0.0_r8
1758  IF (diff(irho).LE.1.d-6) THEN
1759  diff(irho) = 1.d-6
1760  vconv(irho) = (flux(irho) / (max(abs(vpr(irho)), 1.e-6)*g1(irho)) + dy(irho)*diff(irho)) / y(irho)
1761  END IF
1762  END IF
1763 
1764 
1765 
1766 ! +++ Return new ion density and flux profiles to the work flow:
1767  profiles%NE(irho) = ne(irho)
1768  profiles%DNE(irho) = dne(irho) !AF - 25.Sep.2014
1769  profiles%DIFF_NE(irho) = diff(irho)
1770  profiles%VCONV_NE(irho) = vconv(irho)
1771  profiles%FLUX_NE(irho) = flux(irho)
1772  profiles%FLUX_NE_CONV(irho) = flux_ne_conv(irho)
1773 !dy 2017-10-06 PROFILES%SOURCE_NE(IRHO) = SE_EXP(IRHO) + SE_IMP(IRHO) * NE(IRHO)
1774  profiles%SOURCE_NE(irho) = se_exp(irho) - se_imp(irho) * ne(irho)
1775  profiles%INT_SOURCE_NE(irho) = int_source(irho)
1776 
1777 
1778  END DO rho_loop8
1779 
1780  fun1 = profiles%SOURCE_NE(:)*vpr
1781  CALL integr2(nrho,rho,fun1,intfun1)
1782  profiles%INT_SOURCE_NE(:) = intfun1
1783 
1784 
1785  local_int_source_s4 = local_intfun1_s4 + btprime/2.e0_r8/bt*rho(1)*ne(1)/g1(1)
1786  local_flux_ne_s4 = ( y(1)*vconv(1) - dy(1)*diff(1) )
1787 
1788  local_flux_ne_conv_s4 = 0.e0_r8
1789  model_loop3: DO imodel = 1, nmodel
1790  local_flux_ne_conv_s4 = local_flux_ne_conv_s4 + c1(imodel) &
1791  *( y(1)*vconv_mod(1,imodel) - dy(1)*diff_mod(1,imodel) )
1792  END DO model_loop3
1793 
1794  IF (1.EQ.0) THEN
1795  local_flux_ne_s4 = local_int_source_s4
1796  local_flux_ne_conv_s4 = 1.5e0_r8*local_int_source_s4
1797  END IF
1798 
1799 
1800 
1801 
1802 ! +++ Deallocate types for interface with numerical solver:
1803  CALL deallocate_numerics(solver, ifail)
1804 
1805 
1806  RETURN
1807 
1808 
1809 
1810 END SUBROUTINE electron_density
1811 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1812 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1813 
1814 
1815 
1816 
1817 
1818 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1837 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
1838 SUBROUTINE quasi_neutrality &
1839  (geometry,profiles,impurity,control, ifail, failstring)
1840 
1841 !DIAG)
1842 
1843 !-------------------------------------------------------!
1844 ! This subroutine calculates electron density, !
1845 ! electron flux plasma effective charge and !
1846 ! convective contribution to electron energy !
1847 ! transport from density and flux of background !
1848 ! ions (all ion components computed by the ETS) !
1849 ! and impurity ions (all ion components computed !
1850 ! by separate impurity routine) !
1851 ! using quasi-neutrality condition !
1852 !-------------------------------------------------------!
1853 
1854  USE ets_plasma
1855 
1856  USE s4_parameters !AF 12.Oct.2011 - to avoid a few interface blocks
1857 
1858  IMPLICIT NONE
1859 
1860  INTEGER, INTENT (INOUT) :: ifail
1861  CHARACTER(LEN=500) :: failstring
1862 
1863  TYPE (magnetic_geometry) :: geometry !contains all geometry quantities
1864  TYPE (plasma_profiles) :: profiles !contains profiles of plasma parameters
1865  TYPE (impurity_profiles) :: impurity !contains profiles of impurities calculated
1866  TYPE (run_control) :: control !contains all parameters required by run
1867 
1868 ! INTEGER, INTENT (INOUT) :: ifail
1869 ! CHARACTER(LEN=500) :: failstring
1870  TYPE (diagnostic) :: diag !contains error messages and warnings
1871 
1872 ! +++ Internal parameters:
1873  INTEGER :: iion, nion !index and total number of ions
1874  INTEGER :: iimp, nimp !index and total number of impurity
1875  INTEGER :: izimp,nzimp !index and total number of impurity ionization states
1876  INTEGER :: nion_qn !number of ions computed from QN
1877 
1878  INTEGER :: irho, nrho
1879  REAL (R8) :: rho(profiles%nrho) !normalised minor radius, [m]
1880 
1881  REAL (R8) :: zion(profiles%nion) !charge of background ions
1882  REAL (R8) :: zion2(profiles%nion) !charge of background ions, squared
1883 
1884  REAL (R8) :: ni(profiles%nrho,profiles%nion) !density of background ions, [m^-3]
1885  REAL (R8) :: dni(profiles%nrho,profiles%nion) !AF - 25.Sep.2014 - derivative of density of background ions, [m^-3]
1886  REAL (R8) :: flux_ni(profiles%nrho,profiles%nion) !flux of background ions [s-1]
1887  REAL (R8) :: flux_ni_conv(profiles%nrho,profiles%nion) !ion flux, contributing to convective heat transport [1/s]
1888 
1889  REAL (R8) :: zimp(impurity%nrho,impurity%nimp,impurity%nzimp) !charge of impurity ions
1890  REAL (R8) :: zimp2(impurity%nrho,impurity%nimp,impurity%nzimp) !charge of impurity ions, squared
1891 
1892  REAL (R8) :: nz(impurity%nrho,impurity%nimp,impurity%nzimp) !density of impurity ions, [m^-3]
1893  REAL (R8) :: flux_nz(impurity%nrho,impurity%nimp,impurity%nzimp) !flux ofimpurity ions, [s-1]
1894 
1895  REAL (R8) :: ne(profiles%nrho) !density of electrons, [m^-3]
1896  REAL (R8) :: dne(profiles%nrho) !AF - 25.Sep.2014 - derivative of density of electrons, [m^-4]
1897  REAL (R8) :: flux_ne(profiles%nrho) !flux of electrons [s-1]
1898  REAL (R8) :: flux_ne_conv(profiles%nrho) !electron flux, contributing to convective heat transport [1/s]
1899 
1900  REAL (R8) :: zeff(profiles%nrho) !plasma effective charge
1901  REAL (R8) :: nz2(profiles%nrho) !density*Z^2 [m^-3]
1902 
1903  REAL (R8) :: znion(profiles%nrho) !total ion_density*Z [m^-3]
1904  REAL (R8) :: zflux(profiles%nrho) !total ion_flux*Z [s-1]
1905  REAL (R8) :: ztot, ftot
1906  REAL (R8) :: fion(profiles%nion)
1907 
1908  !AF 12.Oct.2011
1909  REAL (R8), ALLOCATABLE :: local_int_source_s4(:) !auxiliary variable with the integral of sources divided by the metrics coefficients Vprime and G1
1910  !AF - End
1911 
1912  !AF - 25.Sep.2014
1913  REAL (R8) :: aux1(profiles%nrho)
1914  REAL (R8) :: aux2(profiles%nrho)
1915  !AF - End
1916 
1917  REAL(R8) :: ni_tot(profiles%nrho)
1918  REAL(R8) :: flux_ni_tot(profiles%nrho)
1919  REAL(R8) :: contrib_2_energy_flux_ni_tot(profiles%nrho)
1920 
1921  REAL(R8) :: ni_z2(profiles%nrho)
1922  REAL(R8) :: ion_fraction(profiles%nion)
1923  REAL(R8) :: ni_qn, zion_av
1924 
1925 !dy
1926  integer :: irho_err
1927 
1928 ! +++ Set up dimensions:
1929  nrho = profiles%NRHO
1930  nion = profiles%NION
1931  nimp = impurity%NIMP
1932  nzimp = impurity%NZIMP
1933 
1934 
1935 ! +++ Set up local variables:
1936 
1937  rho = geometry%RHO
1938 
1939  ne = profiles%NE
1940  dne = profiles%DNE !AF, 25.Sep.2014
1941  flux_ne = profiles%FLUX_NE
1942  flux_ne_conv = profiles%FLUX_NE_CONV
1943 
1944 
1945  ni = profiles%NI
1946  dni = profiles%DNI !AF - 25.Sep.2014
1947  flux_ni = profiles%FLUX_NI
1948  flux_ni_conv = profiles%FLUX_NI_CONV
1949  zion = profiles%ZION
1950  zion2 = profiles%ZION2
1951 
1952 
1953  nz = impurity%NZ
1954  flux_nz = impurity%FLUX_NZ
1955  zimp = impurity%ZIMP
1956  zimp2 = impurity%ZIMP2
1957 
1958 
1959 
1960 ! +++ Quasineutrality condition:
1961 
1962 
1963  quasi_neutralty_condition: IF (control%QUASI_NEUT.EQ.0) THEN
1964 ! ELECTRON density & flux obtained from QN
1965  profiles%NE = 0.0_r8
1966  profiles%DNE = 0.0_r8 !AF
1967  profiles%FLUX_NE = 0.0_r8
1968  profiles%CONTRIB_2_ENERGY_FLUX_NE = 0.0_r8
1969 
1970  ion_contribution: DO iion=1,nion
1971  profiles%NE(:) = profiles%NE(:) + profiles%ZION(iion)*(profiles%NI(:,iion) + profiles%NI_FAST(:,iion))
1972  profiles%FLUX_NE(:) = profiles%FLUX_NE(:) + profiles%ZION(iion)*profiles%FLUX_NI(:,iion)
1973  profiles%CONTRIB_2_ENERGY_FLUX_NE(:) = profiles%CONTRIB_2_ENERGY_FLUX_NE(:) &
1974  + profiles%ZION(iion)*profiles%CONTRIB_2_ENERGY_FLUX_NI(:,iion)
1975  END DO ion_contribution
1976 
1977  impurity_contribution: DO iimp=1,nimp
1978  impurity_charge_states: DO izimp=1,nzimp
1979  profiles%NE(:) = profiles%NE(:) + impurity%ZIMP(:,iimp,izimp)*impurity%NZ(:,iimp,izimp)
1980  profiles%FLUX_NE(:) = profiles%FLUX_NE(:) + impurity%ZIMP(:,iimp,izimp)*impurity%FLUX_NZ(:,iimp,izimp)
1981  END DO impurity_charge_states
1982  ENDDO impurity_contribution
1983 
1984  profiles%NE(:) = profiles%NE(:) - profiles%NE_FAST(:)
1985 
1986 !
1987 ! PROFILES%ZEFF(IRHO) = NZ2(IRHO)/NE(IRHO)
1988 !
1989 
1990  !AF - 25.Sep.2014
1991  CALL derivn(nrho,rho,ne,profiles%DNE)
1992  !PROFILES%DNE = DNE
1993  !AF - End
1994 
1995  local_flux_ne_conv_s4 = 0.0_r8
1996  ion_loop3: DO iion = 1, nion
1997  local_flux_ne_conv_s4 = local_flux_ne_conv_s4 + profiles%ZION(iion)*local_flux_ni_conv_s4(iion)
1998  END DO ion_loop3
1999 
2000  ELSE IF (control%QUASI_NEUT.EQ.1.OR.control%QUASI_NEUT.EQ.2) THEN
2001 ! TOTAL ION density & flux obtained from QN
2002 
2003  ni_tot = profiles%NE + profiles%NE_FAST
2004  flux_ni_tot = profiles%FLUX_NE
2005  contrib_2_energy_flux_ni_tot = profiles%CONTRIB_2_ENERGY_FLUX_NE
2006 
2007 
2008  subtract_predicted_ions: DO iion=1,nion
2009  ni_tot(:) = ni_tot(:) - profiles%ZION(iion)*profiles%NI_FAST(:,iion)
2010  IF (profiles%NI_BND_TYPE(iion).GT.0.5) THEN
2011  ni_tot(:) = ni_tot(:) - profiles%ZION(iion)*profiles%NI(:,iion)
2012  flux_ni_tot(:) = flux_ni_tot(:) - profiles%ZION(iion)*profiles%FLUX_NI(:,iion)
2013  contrib_2_energy_flux_ni_tot(:) = contrib_2_energy_flux_ni_tot(:) - profiles%ZION(iion)&
2014  *profiles%CONTRIB_2_ENERGY_FLUX_NI(:,iion)
2015  END IF
2016  ENDDO subtract_predicted_ions
2017 
2018  subtract_impurity: DO iimp=1,nimp
2019  impurity_charge_states1: DO izimp=1,nzimp
2020  ni_tot(:) = ni_tot(:) - impurity%ZIMP(:,iimp,izimp)*impurity%NZ(:,iimp,izimp)
2021  flux_ni_tot(:) = flux_ni_tot(:) - impurity%ZIMP(:,iimp,izimp)*impurity%FLUX_NZ(:,iimp,izimp)
2022  END DO impurity_charge_states1
2023  ENDDO subtract_impurity
2024 
2025 
2026 ! Separation for individual ion species
2027  nion_qn = 0
2028  ni_qn = 0.0_r8
2029  zion_av = 0.0_r8
2030 
2031  total_density_for_qn_ions: DO iion = 1, nion
2032  IF (profiles%NI_BND_TYPE(iion).EQ.0) THEN
2033  ni_qn = ni_qn + profiles%NI_BND(1,iion)
2034  nion_qn = nion_qn + 1
2035  END IF
2036  END DO total_density_for_qn_ions
2037  ion_fraction(:) = profiles%NI_BND(1,:)/ni_qn
2038 
2039  ion_gharge_density: DO iion = 1, nion
2040  IF (profiles%NI_BND_TYPE(iion).EQ.0) THEN
2041  zion_av = zion_av + (ion_fraction(iion)*profiles%ZION(iion))
2042  END IF
2043  END DO ion_gharge_density
2044 
2045 
2046  proportional_to_boundary_value: IF (control%QUASI_NEUT.EQ.1) THEN
2047  DO iion = 1, nion
2048  IF (profiles%NI_BND_TYPE(iion).EQ.0) THEN
2049  profiles%NI(:,iion) = ni_tot(:) /zion_av*ion_fraction(iion)
2050  profiles%FLUX_NI(:,iion) = flux_ni_tot(:) /zion_av*ion_fraction(iion)
2051  profiles%CONTRIB_2_ENERGY_FLUX_NI(:,iion) = contrib_2_energy_flux_ni_tot(:)/zion_av*ion_fraction(iion)
2052  END IF
2053  END DO
2054  END IF proportional_to_boundary_value
2055 
2056  inversely_proportional_to_ion_gharge_density: IF (control%QUASI_NEUT.EQ.2) THEN
2057  DO iion = 1, nion
2058  IF (profiles%NI_BND_TYPE(iion).EQ.0) THEN
2059  profiles%NI(:,iion) = ni_tot(:) /nion_qn/profiles%ZION(iion)
2060  profiles%FLUX_NI(:,iion) = flux_ni_tot(:) /nion_qn/profiles%ZION(iion)
2061  profiles%CONTRIB_2_ENERGY_FLUX_NI(:,iion) = contrib_2_energy_flux_ni_tot(:)/nion_qn/profiles%ZION(iion)
2062  END IF
2063  END DO
2064  END IF inversely_proportional_to_ion_gharge_density
2065 
2066 
2067  !AF - 25.Sep.2014
2068  DO iion = 1, nion
2069  aux1 = profiles%NI(:,iion)
2070  aux2 = 0.e0_r8
2071  CALL derivn(nrho, geometry%RHO, aux1, aux2)
2072  profiles%DNI(:,iion) = aux2
2073  END DO
2074  !AF - End
2075 
2076  local_flux_ni_conv_s4 = 0.0_r8
2077  DO iion = 1, nion
2078  IF (profiles%NI_BND_TYPE(iion).EQ.0) THEN
2079  IF (control%QUASI_NEUT.EQ.1) &
2080  !LOCAL_FLUX_NI_CONV_S4 = LOCAL_FLUX_NE_CONV_S4/ZTOT*FION(IION)
2081  local_flux_ni_conv_s4 = local_flux_ne_conv_s4/zion_av*ion_fraction(iion)
2082  IF (control%QUASI_NEUT.EQ.2) &
2083  local_flux_ni_conv_s4 = local_flux_ne_conv_s4/nion/zion(iion)
2084  END IF
2085  END DO
2086 
2087 
2088  END IF quasi_neutralty_condition
2089 
2090 ! Plasma effective charge:
2091  ni_z2 = 0.0_r8
2092  ion_contribution_2_zeff: DO iion = 1, nion
2093  ni_z2(:) = ni_z2(:) + profiles%ZION2(iion)*(profiles%NI(:,iion) + profiles%NI_FAST(:,iion))
2094  END DO ion_contribution_2_zeff
2095  impurity_contribution_2_zeff: DO iimp = 1, nimp
2096  DO izimp =1, nzimp
2097  ni_z2(:) = ni_z2(:) + impurity%ZIMP2(:,iimp,izimp)*impurity%NZ(:,iimp,izimp)
2098  END DO
2099  END DO impurity_contribution_2_zeff
2100 
2101 
2102  profiles%ZEFF = ni_z2 / (profiles%NE + profiles%NE_FAST)
2103 
2104 !dy case when ne is not computed at all
2105 !dy just check that
2106  if (control%QUASI_NEUT.EQ.-1) then
2107  do irho=1,nrho
2108  if (ne(irho).le.0.0) then
2109  if (irho.eq.1) then
2110  failstring='Error in the density equation: on-axis electron density is negative, stop'
2111  write(*,*) 'value of on-axis electron density is', ne(1)
2112  write(*,*) 'check ne profile in the file err_tmp'
2113  open (22,file='err_tmp')
2114  do irho_err=1,nrho
2115  write(22,*) irho_err, ne(irho_err)
2116  enddo
2117  close(22)
2118  ifail=-1
2119  return
2120  else
2121  write(*,*) 'warning, electron density for ', irho,' is negative, set it to the value at the previous irho'
2122  ne(irho)=ne(irho-1)
2123  end if
2124  end if
2125  enddo
2126  profiles%ne=ne
2127  end if
2128 
2129 
2130 
2131 ! DO IION = 1, NION
2132 ! IF (PROFILES%NI_BND_TYPE(IION).EQ.0) THEN
2133 ! FTOT = FTOT + PROFILES%NI_BND(1,IION)
2134 ! END IF
2135 ! END DO
2136 ! FION(:) = PROFILES%NI_BND(1,:)/FTOT
2137 
2138 ! DO IION = 1, NION
2139 ! IF (PROFILES%NI_BND_TYPE(IION).EQ.0) THEN
2140 ! ZTOT = ZTOT + (FION(IION)*ZION(IION))
2141 ! NI(:,IION) = 0.e0_R8
2142 ! FLUX_NI(:,IION) = 0.e0_R8
2143 ! FLUX_NI_CONV(:,IION) = 0.e0_R8
2144 ! END IF
2145 ! END DO
2146 
2147 
2148 ! RHO_LOOP3: DO IRHO=1,NRHO
2149 
2150 ! total ion density and flux multiplied by charge
2151 ! IMPURITY_LOOP3: DO IIMP=1,NIMP
2152 ! IMPURITY_CHARGE_LOOP3: DO IZIMP=1,NZIMP
2153 ! ZNION(IRHO) = ZNION(IRHO) - ZIMP(IRHO,IIMP,IZIMP)*NZ(IRHO,IIMP,IZIMP)
2154 ! ZFLUX(IRHO) = ZFLUX(IRHO) - ZIMP(IRHO,IIMP,IZIMP)*FLUX_NZ(IRHO,IIMP,IZIMP)
2155 ! NZ2(IRHO) = NZ2(IRHO) + ZIMP2(IRHO,IIMP,IZIMP)*NZ(IRHO,IIMP,IZIMP)
2156 ! END DO IMPURITY_CHARGE_LOOP3
2157 ! ENDDO IMPURITY_LOOP3
2158 
2159 ! ION_LOOP4: DO IION=1,NION
2160 ! IF (PROFILES%NI_BND_TYPE(IION).GT.0.5) THEN
2161 ! ZNION(IRHO) = ZNION(IRHO) - ZION(IION)*NI(IRHO,IION)
2162 ! ZFLUX(IRHO) = ZFLUX(IRHO) - ZION(IION)*FLUX_NI(IRHO,IION)
2163 ! FLUX_NE_CONV(IRHO) = FLUX_NE_CONV(IRHO) - ZION(IION)*FLUX_NI_CONV(IRHO,IION)
2164 ! END IF
2165 ! ENDDO ION_LOOP4
2166 
2167 
2168 ! +++ Separation for individual species:
2169 
2170 ! proportional to boundary value
2171 ! IF (CONTROL%QUASI_NEUT.EQ.1) THEN
2172 ! DO IION = 1, NION
2173 ! IF (PROFILES%NI_BND_TYPE(IION).EQ.0) THEN
2174 ! NI(IRHO,IION) = ZNION(IRHO)/ZTOT*FION(IION)
2175 ! FLUX_NI(IRHO,IION) = ZFLUX(IRHO)/ZTOT*FION(IION)
2176 ! FLUX_NI_CONV(IRHO,IION) = FLUX_NE_CONV(IRHO)/ZTOT*FION(IION)
2177 ! END IF
2178 ! END DO
2179 ! END IF
2180 
2181 ! proportional to charge
2182 ! IF (CONTROL%QUASI_NEUT.EQ.2) THEN
2183 ! DO IION = 1, NION
2184 ! IF (PROFILES%NI_BND_TYPE(IION).EQ.0) THEN
2185 ! NI(IRHO,IION) = ZNION(IRHO)/NION/ZION(IION)
2186 ! FLUX_NI(IRHO,IION) = ZFLUX(IRHO)/NION/ZION(IION)
2187 ! FLUX_NI_CONV(IRHO,IION) = FLUX_NE_CONV(IRHO)/NION/ZION(IION)
2188 ! END IF
2189 ! END DO
2190 ! END IF
2191 
2192 ! DO IION = 1, NION
2193 ! NZ2(IRHO) = NZ2(IRHO) + ZION2(IION)*NI(IRHO,IION)
2194 !
2195 ! +++ Return ion density, flux and effective charge:
2196 ! IF (PROFILES%NI_BND_TYPE(IION).EQ.0) THEN
2197 ! IF (NI(IRHO,IION).LE.0.0_R8) THEN
2198 ! WRITE(*,*)'=================================='
2199 ! WRITE(*,*)'=================================='
2200 ! WRITE(*,*)'== Quasi Neutrality is broken! =='
2201 ! WRITE(*,*)'== If you are solving electrons =='
2202 ! WRITE(*,*)'== and ions simultaneously, =='
2203 ! WRITE(*,*)'== please check your transport =='
2204 ! WRITE(*,*)'== and source terms. =='
2205 ! WRITE(*,*)'=================================='
2206 ! WRITE(*,*)'=================================='
2207 ! ifail = -2 !none acceptable error!
2208 ! failstring = TRIM(failstring)//'; Quasi Neutrality is broken!'
2209 ! END IF
2210 ! PROFILES%NI(IRHO,IION) = NI(IRHO,IION)
2211 ! PROFILES%FLUX_NI(IRHO,IION) = FLUX_NI(IRHO,IION)
2212 ! PROFILES%FLUX_NI_CONV(IRHO,IION) = FLUX_NI_CONV(IRHO,IION)
2213 ! END IF
2214 ! END DO
2215 
2216 ! PROFILES%ZEFF(IRHO) = NZ2(IRHO)/NE(IRHO)
2217 ! END DO RHO_LOOP3
2218 
2219 
2220  RETURN
2221 
2222 END SUBROUTINE quasi_neutrality
2223 
2224 
2225 
2226 
2227 
2228 
2229 ! +++ HEAT TRANSPORT EQUATIONS +++
2230 
2231 
2232 
2233 
2234 
2235 
2236 
2237 
2238 
2239 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
2289 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
2290 SUBROUTINE temperatures &
2291 ! (GEOMETRY,PROFILES,TRANSPORT,SOURCES,EVOLUTION,CONTROL,ifail,failstring)
2292  (geometry,profiles,transport,sources,evolution,control,hyper_diff,ifail,failstring) !AF 25.Apr.2016, 22.Aug.2016
2293 
2294 !-------------------------------------------------------!
2295 ! This subroutine solves transport equations !
2296 ! for ion components from 1 to NION and electrons, !
2297 ! and provides: temperatures, heat fluxes and !
2298 ! its convective and conductive components !
2299 !-------------------------------------------------------!
2300 
2301  USE ets_plasma
2302  USE type_solver
2303  USE itm_constants
2304  USE plasma_collisionality
2305 
2306 
2307  USE s4_parameters !AF 12.Oct.2011 - to avoid a few interface blocks
2308 
2309  IMPLICIT NONE
2310 
2311  INTEGER, INTENT (INOUT) :: ifail
2312  CHARACTER(LEN=500) :: failstring
2313 
2314 ! +++ External parameters:
2315  TYPE (magnetic_geometry) :: geometry !contains all geometry quantities
2316  TYPE (plasma_profiles) :: profiles !contains profiles of plasma parameters
2317  TYPE (transport_coefficients) :: transport !contains profiles of trasport coefficients
2318  TYPE (sources_and_sinks) :: sources !contains profiles of sources
2319  TYPE (time_evolution) :: evolution !contains all parameters required by time evolution
2320  TYPE (collisionality) :: collisions !contains all parameters computed by module for plasma collision
2321  TYPE (run_control) :: control !contains all parameters required by run
2322 !evolution and convergence loop
2323 ! +++ Input/Output to numerical solver:
2324  TYPE (numerics) :: solver !contains all I/O quantities to numerics part
2325 
2326 ! +++ Internal parameters:
2327  INTEGER :: irho, nrho !radius index, number of radial points
2328  INTEGER :: iion, nion !index of ion component, number of considered ion components
2329  INTEGER :: zion !second index of ion component
2330 
2331  REAL (R8) :: bt, btm, btprime !magnetic field, from current time step,[T], previous time steps, [T], time derivative, [T/s]
2332  REAL (R8) :: rho(profiles%nrho) !normalised minor radius, [m]
2333  REAL (R8) :: vpr(profiles%nrho) !V', [m^2]
2334  REAL (R8) :: dvpr(profiles%nrho) !V'', [m]
2335  REAL (R8) :: vprm(profiles%nrho) !V' (at previous time step), [m^2]
2336  REAL (R8) :: g1(profiles%nrho)
2337 
2338  REAL (R8) :: ti(profiles%nrho), ni(profiles%nrho) !ion temperature, [eV], ion density, [m^-3]
2339  REAL (R8) :: tim(profiles%nrho), nim(profiles%nrho) !ion temperature (at time-1),[eV], ion density (at time-1), [m^-3]
2340  REAL (R8) :: dti(profiles%nrho), dtim(profiles%nrho) !AF - 25.Sep.2014
2341  REAL (R8) :: dni(profiles%nrho), dnim(profiles%nrho) !AF - 25.Sep.2014
2342 
2343  REAL (R8) :: flux_ti(profiles%nrho) !total ion heat flux, [W]
2344  REAL (R8) :: flux_ti_cond(profiles%nrho) !conductive ion heat flux, [W]
2345  REAL (R8) :: flux_ti_conv(profiles%nrho) !convective ion heat flux, [W]
2346  REAL (R8) :: int_source(profiles%nrho) !integral ion heat source [W]
2347  REAL (R8) :: flux_ni(profiles%nrho) !ion particle flux contributing to energy transport, [1/s]
2348 
2349  REAL (R8) :: qi_exp(profiles%nrho) !explicit external ion heat source, [eV/(m^3*s)]
2350  REAL (R8) :: qi_imp(profiles%nrho) !implicit external ion heat source, (proportional Ti) [1/(m^3*s)]
2351  REAL (R8) :: vei(profiles%nrho) !energy exchange between electrons and ions, [1/(m^3*s)]
2352  REAL (R8) :: vzi(profiles%nrho) !energy exchange to all other ion components, [1/(m^3*s)]
2353  REAL (R8) :: vii(profiles%nrho,profiles%nion) !energy exchange between different ion components, [1/(m^3*s)]
2354  REAL (R8) :: qei(profiles%nrho) !ion heat source due to interaction with electrons, [eV/(m^3*s)]
2355  REAL (R8) :: qzi(profiles%nrho) !ion heat source due to interaction with all other ion comp., [eV/(m^3*s)]
2356  REAL (R8) :: qgi(profiles%nrho) !flow-work term in energy exchange, [eV/(m^3*s)]
2357 
2358  REAL (R8) :: ti_bnd(2,3) !boundary condition, value, [depends on TI_BND_TYPE]
2359  INTEGER :: ti_bnd_type(2,profiles%nion) !boundary condition, type
2360 
2361  REAL (R8) :: te(profiles%nrho), tem(profiles%nrho) !electron temperature from current ans previous time step, [eV]
2362  REAL (R8) :: ne(profiles%nrho), nem(profiles%nrho) !electron density from current and previous time step, [m^-3]
2363  REAL (R8) :: dte(profiles%nrho), dtem(profiles%nrho) !AF - 25.Sep.2014
2364  REAL (R8) :: dne(profiles%nrho), dnem(profiles%nrho) !AF - 25.Sep.2014
2365 
2366  REAL (R8) :: flux_te(profiles%nrho) !total electron heat flux, [W]
2367  REAL (R8) :: flux_te_cond(profiles%nrho) !non convective electron heat flux, [W]
2368  REAL (R8) :: flux_te_conv(profiles%nrho) !non convective electron heat flux, [W]
2369  REAL (R8) :: flux_ne(profiles%nrho) !electron particle flux, [1/s]
2370 
2371  REAL (R8) :: qe_exp(profiles%nrho) !explicit external electron heat source, [eV/(m^3*s)]
2372  REAL (R8) :: qe_imp(profiles%nrho) !implicit external electron heat source, (proportional Te) [1/(m^3*s)]
2373  REAL (R8) :: vie(profiles%nrho) !energy exchange between electrons and ions, [1/(m^3*s)]
2374  REAL (R8) :: qie(profiles%nrho) !ion heat source due to interaction with electrons, [eV/(m^3*s)]
2375 
2376  REAL (R8) :: te_bnd(2,3) !boundary condition, value, [depends on TE_BND_TYPE]
2377  INTEGER :: te_bnd_type(2) !boundary condition, type
2378 
2379 ! REAL (R8) :: DIFF(PROFILES%NRHO), VCONV(PROFILES%NRHO) !heat diffusion coefficient, [m^2/s], heat pinch velocity, [m/s]
2380  REAL (R8) :: diff_te(profiles%nrho), vconv_te(profiles%nrho) !heat diffusion coefficient, [m^2/s], heat pinch velocity, [m/s]
2381  REAL (R8) :: diff_ti(profiles%nrho), vconv_ti(profiles%nrho) !heat diffusion coefficient, [m^2/s], heat pinch velocity, [m/s]
2382 
2383  REAL (R8) :: amix, tau !mixing factor, time step, [s]
2384 
2385  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
2386  INTEGER :: idim, ndim !index and total number of equations to be solved
2387  INTEGER :: solver_type !specifies the option for numerical solution
2388  REAL (R8) :: y(profiles%nrho) !function at the current amd previous time steps
2389  REAL (R8) :: ym(profiles%nrho) !function at the current amd previous time steps
2390  REAL (R8) :: dy(profiles%nrho) !derivative of function
2391  REAL (R8) :: dym(profiles%nrho) !derivative of function at previous time step !AF - 25.Sep.2014
2392  REAL (R8) :: a(profiles%nrho) !coefficients for numerical solver
2393  REAL (R8) :: b(profiles%nrho) !coefficients for numerical solver
2394  REAL (R8) :: c(profiles%nrho) !coefficients for numerical solver
2395  REAL (R8) :: d(profiles%nrho) !coefficients for numerical solver
2396  REAL (R8) :: e(profiles%nrho) !coefficients for numerical solver
2397  REAL (R8) :: f(profiles%nrho) !coefficients for numerical solver
2398  REAL (R8) :: g(profiles%nrho) !coefficients for numerical solver
2399  REAL (R8) :: h !coefficients for numerical solver
2400  REAL (R8) :: v(2), u(2), w(2) !boundary conditions for numerical solver
2401 
2402  REAL (R8) :: fun1(profiles%nrho), intfun1(profiles%nrho)
2403  REAL (R8) :: fun2(profiles%nrho), intfun2(profiles%nrho)
2404 
2405  real (r8), parameter :: twothird = 2.0_r8/3.0_r8, &
2406  fivethird = 5.0_r8/3.0_r8
2407 ! real (r8), parameter :: twothird = 0.66667, &
2408 ! fivethird = 1.66667
2409 
2410 ! +++ Stabilization scheme !AF 25.Apr.2016, 22.Aug.2016
2411  REAL (R8), DIMENSION(2) :: hyper_diff !hyper diffusivity
2412  REAL (R8) :: hyper_diff_exp !explicit hyper diffusivity [m^2/s]
2413  REAL (R8) :: hyper_diff_imp !implicit hyper diffusivity
2414  REAL (R8) :: diff_hyper !additional diffusivity to remove profile oscillations
2415 
2416 !dy limit diffusion coefficients
2417  REAL (R8), parameter :: tite_diff_limit=1.0e4
2418 !dy
2419  integer :: irho_err
2420 
2421 
2422  hyper_diff_exp = hyper_diff(1) !AF 22.Aug.2016
2423  hyper_diff_imp = hyper_diff(2) !AF 22.Aug.2016
2424 
2425 ! +++ Set up dimensions:
2426  nrho = profiles%NRHO
2427  nion = profiles%NION
2428  ndim = nion + 1
2429 
2430 
2431 ! +++ Allocate types for interface with PLASMA_COLLISIONS:
2432  CALL allocate_collisionality(nrho, nion, collisions, ifail)
2433 
2434 
2435 ! +++ Allocate types for interface with numerical solver:
2436  CALL allocate_numerics(ndim, nrho, solver, ifail)
2437 
2438 
2439 
2440 ! +++ Set up local variables:
2441  amix = control%AMIX
2442  tau = control%TAU
2443  solver_type = control%SOLVER_TYPE
2444 
2445  bt = geometry%BGEO
2446  btm = evolution%BTM
2447  btprime = (bt-btm)/tau
2448 
2449 
2450  rho_loop1: DO irho=1,nrho
2451  rho(irho) = geometry%RHO(irho)
2452  vpr(irho) = geometry%VPR(irho)
2453  vprm(irho) = evolution%VPRM(irho)
2454  g1(irho) = geometry%G1(irho)
2455  END DO rho_loop1
2456 
2457  CALL derivn(nrho,rho,vpr,dvpr) !Derivation of V'
2458 
2459 
2460 
2461 ! +++ Energy exchange terms due to collisions
2462 ! (defined from previous iteration):
2463  CALL plasma_collisions(geometry,profiles,collisions,ifail)
2464 
2465  vie = collisions%VIE !DPC 2009-01-19
2466 
2467 
2468 ! +++ Parameters for numerical solver
2469 ! (common for all components):
2470  solver%TYPE = solver_type
2471  solver%NDIM = ndim
2472  solver%NRHO = nrho
2473  solver%AMIX = amix
2474  solver%DERIVATIVE_FLAG(1) = 0
2475 
2476  rho_loop2: DO irho=1,nrho
2477  solver%RHO(irho) = rho(irho)
2478  END DO rho_loop2
2479 
2480 
2481 
2482 !-------------------------------------------------------!
2483 ! !
2484 ! ION HEAT TRANSPORT: !
2485 ! !
2486 !-------------------------------------------------------!
2487 
2488  ion_loop1: DO iion=1,nion
2489 
2490 ! +++ Set equation to 'predictive' and all coefficients to zero:
2491  flag = 1
2492  y(:) = 0.0e0_r8
2493  dy(:) = 0.0e0_r8
2494  ym(:) = 0.0e0_r8
2495  dym(:) = 0.0e0_r8 !AF - 25.Sep.2014
2496  a(:) = 0.0e0_r8
2497  b(:) = 0.0e0_r8
2498  c(:) = 0.0e0_r8
2499  d(:) = 0.0e0_r8
2500  e(:) = 0.0e0_r8
2501  f(:) = 0.0e0_r8
2502  g(:) = 0.0e0_r8
2503  h = 0.0e0_r8
2504  v(:) = 0.0e0_r8
2505  u(:) = 0.0e0_r8
2506  w(:) = 0.0e0_r8
2507 
2508 
2509 
2510 ! +++ Set up boundary conditions for particular ion type:
2511  ti_bnd_type(2,iion) = profiles%TI_BND_TYPE(iion)
2512  ti_bnd(2,1) = profiles%TI_BND(1,iion)
2513  ti_bnd(2,2) = profiles%TI_BND(2,iion)
2514  ti_bnd(2,3) = profiles%TI_BND(3,iion)
2515 
2516 
2517 
2518 ! +++ Set up local variables for particular ion type:
2519  rho_loop3: DO irho = 1,nrho
2520  ti(irho) = profiles%TI(irho,iion)
2521  dti(irho) = profiles%DTI(irho,iion) !AF - 25.Sep.2014
2522  tim(irho) = evolution%TIM(irho,iion)
2523  dtim(irho) = evolution%DTIM(irho,iion) !AF - 25.Sep.2014
2524  ni(irho) = profiles%NI(irho,iion)
2525  dni(irho) = profiles%DNI(irho,iion) !AF - 25.Sep.2014
2526  nim(irho) = evolution%NIM(irho,iion)
2527  dnim(irho) = evolution%DNIM(irho,iion) !AF - 25.Sep.2014
2528 
2529  flux_ti(irho) = profiles%FLUX_TI(irho,iion)
2530  flux_ti_cond(irho) = profiles%FLUX_TI_COND(irho,iion)
2531  flux_ti_conv(irho) = profiles%FLUX_TI_CONV(irho,iion)
2532  flux_ni(irho) = profiles%FLUX_NI_CONV(irho,iion)
2533 
2534  diff_ti(irho) = transport%DIFF_TI(irho,iion)
2535  vconv_ti(irho) = transport%VCONV_TI(irho,iion)
2536  qgi(irho) = transport%QGI(irho,iion)
2537 
2538  qi_exp(irho) = sources%QI_EXP(irho,iion)
2539  qi_imp(irho) = sources%QI_IMP(irho,iion)
2540 
2541  vei(irho) = collisions%VEI(irho,iion)
2542  qei(irho) = collisions%QEI(irho,iion)
2543  vzi(irho) = collisions%VZI(irho,iion)
2544  qzi(irho) = collisions%QZI(irho,iion)
2545 
2546  ion_loop2: DO zion = 1,nion
2547 !!DPC was VII(IRHO,IION) = COLLISIONS%VII(IRHO,IION,ZION) !!! check if the following is what was intended
2548  vii(irho,zion) = collisions%VII(irho,iion,zion)
2549  END DO ion_loop2
2550 
2551 
2552 
2553  END DO rho_loop3
2554 
2555 
2556 
2557 ! +++ Coefficients for ion heat transport equation in form:
2558 !
2559 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
2560 
2561  diff_hyper = hyper_diff_exp + hyper_diff_imp*maxval(diff_ti) !AF 25.Apr.2016, 22.Aug.2016
2562 
2563  rho_loop4: DO irho=nrho,1,-1
2564  y(irho) = ti(irho)
2565  dy(irho) = dti(irho) !AF - 25.Sep.2014
2566  ym(irho) = tim(irho)
2567  dym(irho) = dtim(irho) !AF - 25.Sep.2014
2568 
2569  a(irho) = 1.5e0_r8*vpr(irho)*ni(irho)
2570 !DPC temporary "fix" for NaNs
2571  if(vpr(irho).le.0.0_r8.and.irho.eq.1) then
2572 ! write(*,*) B(2:5)
2573  b(irho)=0.0_r8
2574  else
2575  b(irho) = 1.5e0_r8*vprm(irho)**fivethird/vpr(irho)**twothird*nim(irho)
2576  endif
2577 !DPC end of temporary fix
2578  c(irho) = 1.e0_r8
2579 ! D(IRHO) = VPR(IRHO)*G1(IRHO)*NI(IRHO)*DIFF_TI(IRHO)
2580  d(irho) = vpr(irho)*g1(irho)*ni(irho)*(diff_ti(irho)+diff_hyper) !AF 25.Apr.2016
2581  ! E(IRHO) = VPR(IRHO)*G1(IRHO)*NI(IRHO)*VCONV_TI(IRHO) &
2582  e(irho) = vpr(irho)*g1(irho)*ni(irho)*(vconv_ti(irho)+diff_hyper*dtim(irho)/tim(irho)) & !AF 25.Apr.2016
2583  + flux_ni(irho) &
2584  - 1.5e0_r8*btprime/2.e0_r8/bt*rho(irho)*ni(irho)*vpr(irho)
2585  f(irho) = vpr(irho) &
2586  * (qi_exp(irho) + qei(irho) + qzi(irho) + qgi(irho))
2587  g(irho) = vpr(irho)*(qi_imp(irho) + vei(irho) + vzi(irho)) &
2588  - btprime/2.e0_r8/bt*rho(irho)*ni(irho)*dvpr(irho)
2589  END DO rho_loop4
2590 
2591  h = tau
2592 
2593 
2594 
2595 ! +++ Boundary conditions for ion heat transport equation in form:
2596 !
2597 ! V*Y' + U*Y =W
2598 
2599 ! +++ On axis:
2600 
2601 ! dTi/drho(rho=0)=0
2602  IF (solver_type.NE.4) THEN !AF 11.Oct.2011
2603  v(1) = 1.e0_r8
2604  u(1) = 0.e0_r8
2605  ELSE !AF 11.Oct.2011 - Zero flux instead of zero gradient at the axis for solver 4
2606 ! IF (DIFF_TI(1).GT.1.0E-6) THEN !AF 19.Mar.2012 - To avoid problems with the axis boundary condition
2607  IF ((diff_ti(1)+diff_hyper).GT.1.0e-6) THEN !AF 19.Mar.2012 - To avoid problems with the axis boundary condition !AF 25.Apr.2016
2608  ! V(1) = -DIFF_TI(1)*NI(1)
2609  v(1) = -(diff_ti(1)+diff_hyper)*ni(1) !AF 25.Apr.2016
2610  ELSE
2611  v(1) = -1.0e-6*ni(1)
2612  ENDIF
2613 ! U(1) = VCONV_TI(1)*NI(1)+LOCAL_FLUX_NI_CONV_S4(IION)
2614  u(1) = (vconv_ti(1)+diff_hyper*dtim(1)/tim(1))*ni(1)+local_flux_ni_conv_s4(iion) !AF 25.Apr.2016
2615  END IF !AF 11.Oct.2011
2616  w(1) = 0.e0_r8
2617 
2618 ! +++ At the edge:
2619 
2620 ! FIXED Ti
2621  IF(ti_bnd_type(2,iion).EQ.1) THEN
2622  v(2) = 0.e0_r8
2623  u(2) = 1.e0_r8
2624  w(2) = ti_bnd(2,1)
2625  ENDIF
2626 
2627 ! FIXED grad_Ti
2628  IF(ti_bnd_type(2,iion).EQ.2) THEN
2629  v(2) = 1.e0_r8
2630  u(2) = 0.e0_r8
2631  w(2) = -ti_bnd(2,1)
2632  ENDIF
2633 
2634 ! FIXED L_Ti
2635  IF(ti_bnd_type(2,iion).EQ.3) THEN
2636  v(2) = ti_bnd(2,1)
2637  u(2) = 1.e0_r8
2638  w(2) = 0.e0_r8
2639  ENDIF
2640 
2641 ! FIXED Flux_Ti
2642  IF(ti_bnd_type(2,iion).EQ.4) THEN
2643  v(2) = -vpr(nrho)*g1(nrho)*diff_ti(nrho)*ni(nrho)
2644  u(2) = vpr(nrho)*g1(nrho)*vconv_ti(nrho)*ni(nrho)+flux_ni(nrho)
2645  w(2) = ti_bnd(2,1)
2646  ENDIF
2647 
2648 ! Generic boundary condition
2649  IF(ti_bnd_type(2,iion).EQ.5) THEN
2650  v(2) = ti_bnd(2,1)
2651  u(2) = ti_bnd(2,2)
2652  w(2) = ti_bnd(2,3)
2653  ENDIF
2654 
2655 
2656 
2657 ! +++ Temperature equation is not solved:
2658  IF(ti_bnd_type(2,iion).EQ.0) THEN
2659 
2660  CALL derivn(nrho,rho,y,dy) !temperature gradient
2661 
2662  flag = 0
2663 
2664  rho_loop5: DO irho=1,nrho
2665  a(irho) = 1.0e0_r8
2666  b(irho) = 1.0e0_r8
2667  c(irho) = 1.0e0_r8
2668  d(irho) = 0.0e0_r8
2669  e(irho) = 0.0e0_r8
2670  f(irho) = 0.0e0_r8
2671  g(irho) = 0.0e0_r8
2672  END DO rho_loop5
2673 
2674  v(2) = 0.0e0_r8
2675  u(2) = 1.0e0_r8
2676  w(2) = y(nrho)
2677  END IF
2678 
2679 
2680 
2681 ! +++ Defining coefficients for numerical solver:
2682  solver%EQ_FLAG(iion) = flag
2683 
2684  rho_loop6: DO irho=1,nrho
2685  solver%Y(iion,irho) = y(irho)
2686  solver%DY(iion,irho) = dy(irho)
2687  solver%YM(iion,irho) = ym(irho)
2688 
2689  solver%A(iion,irho) = a(irho)
2690  solver%B(iion,irho) = b(irho)
2691  solver%C(iion,irho) = c(irho)
2692  solver%D(iion,irho) = d(irho)
2693  solver%E(iion,irho) = e(irho)
2694  solver%F(iion,irho) = f(irho)
2695  solver%G(iion,irho) = g(irho)
2696 
2697  solver%CM1(ndim,iion,irho) = vie(irho)
2698  solver%CM1(iion,ndim,irho) = vie(irho)
2699 
2700  ion_loop3: DO zion = 1,nion
2701  solver%CM1(iion,zion,irho) = vii(irho,zion)
2702  END DO ion_loop3
2703 
2704  END DO rho_loop6
2705 
2706  solver%H = h
2707 
2708  solver%V(iion,1) = v(1)
2709  solver%U(iion,1) = u(1)
2710  solver%W(iion,1) = w(1)
2711  solver%V(iion,2) = v(2)
2712  solver%U(iion,2) = u(2)
2713  solver%W(iion,2) = w(2)
2714 
2715 
2716 
2717  END DO ion_loop1
2718 
2719 
2720 
2721 !-------------------------------------------------------!
2722 ! !
2723 ! ELECTRON HEAT TRANSPORT: !
2724 ! !
2725 !-------------------------------------------------------!
2726 
2727 ! +++ Set equation to 'predictive' and all coefficients to zero:
2728  flag = 1
2729  y(:) = 0.0e0_r8
2730  dy(:) = 0.0e0_r8
2731  ym(:) = 0.0e0_r8
2732  dym(:) = 0.0e0_r8 !AF - 25.Sep.2014
2733  a(:) = 0.0e0_r8
2734  b(:) = 0.0e0_r8
2735  c(:) = 0.0e0_r8
2736  d(:) = 0.0e0_r8
2737  e(:) = 0.0e0_r8
2738  f(:) = 0.0e0_r8
2739  g(:) = 0.0e0_r8
2740  h = 0.0e0_r8
2741  v(:) = 0.0e0_r8
2742  u(:) = 0.0e0_r8
2743  w(:) = 0.0e0_r8
2744 
2745 
2746 
2747 ! +++ Set up local variables for electron heat transport equation:
2748  rho_loop7: DO irho=1,nrho
2749  te(irho) = profiles%TE(irho)
2750  dte(irho) = profiles%DTE(irho) !AF - 25.Sep.2014
2751  ne(irho) = profiles%NE(irho)
2752  dne(irho) = profiles%DNE(irho) !AF - 25.Sep.2014
2753  tem(irho) = evolution%TEM(irho)
2754  dtem(irho) = evolution%DTEM(irho) !AF - 25.Sep.2014
2755  nem(irho) = evolution%NEM(irho)
2756  dnem(irho) = evolution%DNEM(irho) !AF - 25.Sep.2014
2757 
2758  flux_ne(irho) = profiles%FLUX_NE_CONV(irho)
2759 
2760  diff_te(irho) = transport%DIFF_TE(irho)
2761  vconv_te(irho) = transport%VCONV_TE(irho)
2762  qgi(irho) = 0.e0_r8
2763  ion_loop4: DO iion = 1, nion
2764  qgi(irho) = qgi(irho) + transport%QGI(irho,iion)
2765  END DO ion_loop4
2766 
2767 
2768 
2769  qe_exp(irho) = sources%QOH(irho) / itm_ev
2770  qe_imp(irho) = 0.e0_r8
2771 
2772  qe_exp(irho) = qe_exp(irho) + sources%QE_EXP(irho)
2773  qe_imp(irho) = qe_imp(irho) + sources%QE_IMP(irho)
2774 
2775  qie(irho) = collisions%QIE(irho)
2776  vie(irho) = collisions%VIE(irho)
2777 
2778  END DO rho_loop7
2779 
2780 
2781 
2782 ! +++ Set up boundary conditions for electron heat transport equation:
2783  te_bnd_type(2) = profiles%TE_BND_TYPE
2784  te_bnd(2,1) = profiles%TE_BND(1)
2785  te_bnd(2,2) = profiles%TE_BND(2)
2786  te_bnd(2,3) = profiles%TE_BND(3)
2787 
2788 
2789 
2790 ! +++ Coefficients for electron heat transport equation in form:
2791 !
2792 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
2793 
2794  diff_hyper = hyper_diff_exp + hyper_diff_imp*maxval(diff_te) !AF 25.Apr.2016, 22.Aug.2016
2795 
2796  rho_loop8: DO irho=nrho,1,-1
2797  y(irho) = te(irho)
2798  dy(irho) = dte(irho) !AF - 25.Sep.2014
2799  ym(irho) = tem(irho)
2800  dym(irho) = dtem(irho) !AF - 25.Sep.2014
2801 
2802  a(irho) = 1.5e0_r8*vpr(irho)*ne(irho)
2803 !DPC temporary "fix" for NaNs
2804  if(vpr(irho).le.0.0_r8.and.irho.eq.1) then
2805 ! write(*,*) B(2:5)
2806  b(irho)=0.0_r8
2807  else
2808  b(irho) = 1.5e0_r8*vprm(irho)**fivethird/vpr(irho)**twothird*nem(irho)
2809  endif
2810 !DPC end of temporary fix
2811  c(irho) = 1.e0_r8
2812 ! D(IRHO) = VPR(IRHO)*G1(IRHO)*NE(IRHO)*DIFF_TE(IRHO)
2813  d(irho) = vpr(irho)*g1(irho)*ne(irho)*(diff_te(irho)+diff_hyper) !AF 25.Apr.2016
2814 ! E(IRHO) = VPR(IRHO)*G1(IRHO)*NE(IRHO)*VCONV_TE(IRHO) &
2815  e(irho) = vpr(irho)*g1(irho)*ne(irho)*(vconv_te(irho)+diff_hyper*dtem(irho)/tem(irho)) & !AF 25.Apr.2016
2816  + flux_ne(irho) &
2817  - 1.5e0_r8*btprime/2.e0_r8/bt*rho(irho)*ne(irho)*vpr(irho)
2818  f(irho) = vpr(irho) * ( qe_exp(irho) + qie(irho) - qgi(irho) )
2819  g(irho) = vpr(irho) * (qe_imp(irho) + vie(irho)) &
2820  - btprime/2.e0_r8/bt*rho(irho)*ne(irho)*dvpr(irho)
2821  END DO rho_loop8
2822 
2823  h = tau
2824 
2825 
2826 ! +++ Boundary conditions for electron heat
2827 ! transport equation in form:
2828 !
2829 ! V*Y' + U*Y =W
2830 
2831 ! +++ On axis
2832 ! dTe/drho(rho=0)=0:
2833  IF (solver_type.NE.4) THEN !AF 11.Oct.2011
2834  v(1) = 1.e0_r8
2835  u(1) = 0.e0_r8
2836  ELSE !AF 11.Oct.2011 - Zero flux instead of zero gradient at the axis for solver 4
2837 ! IF (DIFF_TE(1).GT.1.0E-6) THEN !AF 19.Mar.2012 - To avoid problems with the axis boundary condition
2838  IF ((diff_te(1)+diff_hyper).GT.1.0e-6) THEN !AF 19.Mar.2012 - To avoid problems with the axis boundary condition !AF 25.Apr.2016
2839 ! V(1) = -DIFF_TE(1)*NE(1)
2840  v(1) = -(diff_te(1)+diff_hyper)*ne(1) !AF 25.Apr.2016
2841  ELSE
2842  v(1) = -1.0e-6*ne(1)
2843  ENDIF
2844 ! U(1) = VCONV_TE(1)*NE(1)+LOCAL_FLUX_NE_CONV_S4
2845  u(1) = (vconv_te(1)+diff_hyper*dtem(1)/tem(1))*ne(1)+local_flux_ne_conv_s4 !AF 25.Apr.2016
2846  END IF !AF 11.Oct.2011
2847  w(1) = 0.e0_r8
2848 
2849 ! +++ At the edge:
2850 
2851 ! FIXED Te
2852  IF(te_bnd_type(2).EQ.1) THEN
2853  v(2) = 0.e0_r8
2854  u(2) = 1.e0_r8
2855  w(2) = te_bnd(2,1)
2856  ENDIF
2857 
2858 ! FIXED grad_Te
2859  IF(te_bnd_type(2).EQ.2) THEN
2860  v(2) = 1.e0_r8
2861  u(2) = 0.e0_r8
2862  w(2) = -te_bnd(2,1)
2863  ENDIF
2864 
2865 ! FIXED L_Te
2866  IF(te_bnd_type(2).EQ.3) THEN
2867  v(2) = te_bnd(2,1)
2868  u(2) = 1.e0_r8
2869  w(2) = 0.e0_r8
2870  ENDIF
2871 
2872 ! FIXED Flux_Te
2873  IF(te_bnd_type(2).EQ.4) THEN
2874 ! V(2) = -G(NRHO)*DIFF(NRHO)*NE(NRHO)
2875 ! U(2) = G(NRHO)*VCONV(NRHO)*NE(NRHO)+FLUX_NE(NRHO)
2876  v(2) = -vpr(nrho)*g1(nrho)*diff_te(nrho)*ne(nrho)
2877  u(2) = vpr(nrho)*g1(nrho)*vconv_te(nrho)*ne(nrho)+flux_ne(nrho)
2878  w(2) = te_bnd(2,1)
2879  ENDIF
2880 
2881 ! Generic boundary condition
2882  IF(te_bnd_type(2).EQ.5) THEN
2883  v(2) = te_bnd(2,1)
2884  u(2) = te_bnd(2,2)
2885  w(2) = te_bnd(2,3)
2886  ENDIF
2887 
2888 
2889 
2890 ! +++ Temperature equation is not solved:
2891  IF(te_bnd_type(2).EQ.0) THEN
2892 
2893  CALL derivn(nrho,rho,y,dy) !temperature gradient
2894 
2895  flag = 0
2896 
2897  rho_loop9: DO irho=1,nrho
2898  a(irho) = 1.0e0_r8
2899  b(irho) = 1.0e0_r8
2900  c(irho) = 1.0e0_r8
2901  d(irho) = 0.0e0_r8
2902  e(irho) = 0.0e0_r8
2903  f(irho) = 0.0e0_r8
2904  g(irho) = 0.0e0_r8
2905  END DO rho_loop9
2906 
2907  v(2) = 0.0e0_r8
2908  u(2) = 1.0e0_r8
2909  w(2) = y(nrho)
2910  END IF
2911 
2912 
2913 
2914 ! +++ Defining coefficients for numerical solver:
2915  solver%TYPE = solver_type
2916  solver%EQ_FLAG(ndim) = flag
2917  solver%NDIM = ndim
2918  solver%NRHO = nrho
2919  solver%AMIX = amix
2920  solver%DERIVATIVE_FLAG(1) = 0
2921 
2922  rho_loop10: DO irho=1,nrho
2923 
2924  solver%RHO(irho) = rho(irho)
2925 
2926  solver%Y(ndim,irho) = y(irho)
2927  solver%DY(ndim,irho) = dy(irho)
2928  solver%YM(ndim,irho) = ym(irho)
2929 
2930  solver%A(ndim,irho) = a(irho)
2931  solver%B(ndim,irho) = b(irho)
2932  solver%C(ndim,irho) = c(irho)
2933  solver%D(ndim,irho) = d(irho)
2934  solver%E(ndim,irho) = e(irho)
2935  solver%F(ndim,irho) = f(irho)
2936  solver%G(ndim,irho) = g(irho)
2937 
2938  solver%CM1(ndim,ndim,irho) = 0.e0_r8
2939 
2940  END DO rho_loop10
2941 
2942  solver%H = h
2943 
2944  solver%V(ndim,1) = v(1)
2945  solver%U(ndim,1) = u(1)
2946  solver%W(ndim,1) = w(1)
2947  solver%V(ndim,2) = v(2)
2948  solver%U(ndim,2) = u(2)
2949  solver%W(ndim,2) = w(2)
2950 
2951 
2952 
2953 ! +++ Solution of heat transport equations:
2954  CALL solution_interface(solver, ifail)
2955 
2956  ! dy check for nans in solution, return if true)
2957  IF (any( isnan(solver%y))) then
2958  failstring='Error in the temperature equation: nans in the solution, stop'
2959  ifail=-1
2960  return
2961  end if
2962 
2963 
2964 ! +++ New temperatures and heat fluxes:
2965 !-------------------------------------------------------!
2966 ! IONS: !
2967 !-------------------------------------------------------!
2968  ion_loop5: DO iion = 1,nion
2969 
2970 ! dpc 2011-08-11: I think we need most of the following
2971  ti_bnd_type(2,iion) = profiles%TI_BND_TYPE(iion)
2972  ti_bnd(2,1) = profiles%TI_BND(1,iion)
2973  ti_bnd(2,2) = profiles%TI_BND(2,iion)
2974  ti_bnd(2,3) = profiles%TI_BND(3,iion)
2975 ! dpc end
2976 
2977  rho_loop11: DO irho=1,nrho
2978 
2979 ! dpc 2011-08-11: I think we need most of the following
2980  ti(irho) = profiles%TI(irho,iion)
2981  tim(irho) = evolution%TIM(irho,iion)
2982  ni(irho) = profiles%NI(irho,iion)
2983  nim(irho) = evolution%NIM(irho,iion)
2984  flux_ti(irho) = profiles%FLUX_TI(irho,iion)
2985  flux_ti_cond(irho) = profiles%FLUX_TI_COND(irho,iion)
2986  flux_ti_conv(irho) = profiles%FLUX_TI_CONV(irho,iion)
2987  flux_ni(irho) = profiles%FLUX_NI_CONV(irho,iion)
2988  diff_ti(irho) = transport%DIFF_TI(irho,iion)
2989  vconv_ti(irho) = transport%VCONV_TI(irho,iion)
2990  qgi(irho) = transport%QGI(irho,iion)
2991  qi_exp(irho) = sources%QI_EXP(irho,iion)
2992  qi_imp(irho) = sources%QI_IMP(irho,iion)
2993  vei(irho) = collisions%VEI(irho,iion)
2994  qei(irho) = collisions%QEI(irho,iion)
2995  vzi(irho) = collisions%VZI(irho,iion)
2996  qzi(irho) = collisions%QZI(irho,iion)
2997 ! dpc end
2998 
2999  y(irho) = solver%Y(iion,irho)
3000  dy(irho) = solver%DY(iion,irho)
3001 
3002  IF(ti_bnd_type(2,iion).EQ.0) THEN
3003  y(:) = profiles%TI(:,iion)
3004  CALL derivn(nrho,rho,y,dy)
3005  END IF
3006 
3007  ti(irho) = y(irho)
3008  dti(irho) = dy(irho) !AF - 25.Sep.2014
3009  if (ti(irho).lt.0.0) then
3010  if (irho.eq.1) then
3011  failstring='Error in the temperature equation: on-axis ion temperature is negative, stop'
3012  write(*,*) 'value of on-axis temperature for ion', iion, 'is', ti(1)
3013  write(*,*) 'check ti profile in the file err_tmp'
3014  open (22,file='err_tmp')
3015  do irho_err=1,nrho
3016  write(22,*) irho_err, solver%y(iion,irho_err)
3017  enddo
3018  close(22)
3019  ifail=-1
3020  return
3021  else
3022  write(*,*) 'warning, temperatuire for ion ',iion,' and irho ',irho,'is negative, set it to the value at the previous irho'
3023  ti(irho)=ti(irho-1)
3024  dti(irho)=dti(irho-1)
3025  end if
3026  end if
3027  IF (rho(irho).NE.0.e0_r8) THEN
3028  fun1(irho) = vpr(irho)/rho(irho)* &
3029  ((1.5e0_r8*nim(irho)*tim(irho)/tau*(vprm(irho)/vpr(irho))**fivethird &
3030  + qi_exp(irho) + qei(irho) + qzi(irho) + qgi(irho)) &
3031  - (1.5e0_r8*ni(irho)/tau + qi_imp(irho) + vei(irho) + vzi(irho) &
3032  - btprime/2.e0_r8/bt*rho(irho)*ni(irho)*dvpr(irho)) * y(irho))
3033  ELSE
3034  fun1(irho) = ((1.5e0_r8*nim(irho)*tim(irho)/tau &
3035  + qi_exp(irho) + qei(irho) + qzi(irho) + qgi(irho)) &
3036  - (1.5e0_r8*ni(irho)/tau + qi_imp(irho) + vei(irho) + vzi(irho) &
3037  - btprime/2.e0_r8/bt*rho(irho)*ni(irho)*dvpr(irho)) * y(irho))
3038  ENDIF
3039  END DO rho_loop11
3040 
3041  CALL integr(nrho,rho,fun1,intfun1) !Integral source
3042 
3043  rho_loop12: DO irho=1,nrho
3044  flux_ti_conv(irho) = y(irho)*flux_ni(irho)
3045 
3046  flux_ti_cond(irho) = vpr(irho)*g1(irho)*ni(irho) &
3047  * ( y(irho)*vconv_ti(irho) - dy(irho)*diff_ti(irho) )
3048 
3049  flux_ti(irho) = flux_ti_conv(irho) + flux_ti_cond(irho)
3050 
3051  int_source(irho) = intfun1(irho) &
3052  + y(irho) * 1.5e0_r8*btprime/2.e0_r8/bt*rho(irho)*ni(irho)*vpr(irho)
3053 
3054 
3055 
3056 ! +++ If equation is not solved, total and conductive ion heat flux
3057 ! are determined from the integral of sources:
3058  IF (ti_bnd_type(2,iion).EQ.0) THEN
3059 
3060  diff_ti(irho) = 1.d-6
3061  flux_ti(irho) = int_source(irho)
3062  flux_ti_cond(irho) = int_source(irho) - flux_ni(irho)*y(irho)
3063 
3064  IF ((vpr(irho)*g1(irho).NE.0.0_r8)) &
3065 !dy limit also DY if less than sqrt(epsilon(1.0))
3066  diff_ti(irho) = - flux_ti_cond(irho) / sign(max(abs(dy(irho)),sqrt(epsilon(1.0))),dy(irho)) / (vpr(irho)*g1(irho)*ni(irho))
3067 !dy further limit diff_ti
3068  if (abs(diff_ti(irho)).ge.tite_diff_limit) &
3069  diff_ti(irho)=sign(tite_diff_limit,diff_ti(irho))
3070  vconv_ti(irho) = 0.0_r8
3071  IF (diff_ti(irho).LE.1.d-6) THEN
3072  diff_ti(irho) = 1.d-6
3073  vconv_ti(irho) = (flux_ti_cond(irho) / (max(abs(vpr(irho)), 1.e-6)*g1(irho)*ni(irho)) + dy(irho)*diff_ti(irho)) /max(abs(y(irho)), 1.e-6)
3074  END IF
3075  END IF
3076 
3077 
3078 
3079 ! +++ Return ion temperature and flux:
3080  profiles%TI(irho,iion) = ti(irho)
3081  profiles%DTI(irho,iion) = dti(irho) !AF, 25.Sep.2014
3082  profiles%DIFF_TI(irho,iion) = diff_ti(irho)
3083  profiles%VCONV_TI(irho,iion) = vconv_ti(irho)
3084  profiles%FLUX_TI_CONV(irho,iion) = flux_ti_conv(irho)
3085  profiles%FLUX_TI_COND(irho,iion) = flux_ti_cond(irho)
3086  profiles%FLUX_TI(irho,iion) = flux_ti(irho)
3087  profiles%SOURCE_TI(irho,iion) = qi_exp(irho) + qei(irho) + qzi(irho) + qgi(irho) &
3088  - (qi_imp(irho) + vei(irho) + vzi(irho)) * ti(irho)
3089  profiles%INT_SOURCE_TI(irho,iion) = int_source(irho)
3090  profiles%QEI_OUT(irho)=profiles%QEI_OUT(irho)+qei(irho)
3091 
3092  END DO rho_loop12
3093 
3094  fun1=profiles%SOURCE_TI(:,iion)*vpr
3095  CALL integr2(nrho,rho,fun1,intfun1)
3096  profiles%INT_SOURCE_TI(:,iion) = intfun1
3097 
3098  END DO ion_loop5
3099 
3100 
3101 
3102 !-------------------------------------------------------!
3103 ! ELECTRONS: !
3104 !-------------------------------------------------------!
3105  rho_loop13: DO irho=1,nrho
3106 
3107 ! dpc 2011-08-11: I think we need most of the following
3108  qgi(irho) = 0.e0_r8
3109  ion_loop6: DO iion = 1, nion
3110  qgi(irho) = qgi(irho) + transport%QGI(irho,iion)
3111  END DO ion_loop6
3112 ! dpc end
3113 
3114  y(irho) = solver%Y(ndim,irho)
3115  dy(irho) = solver%DY(ndim,irho)
3116 
3117  IF(te_bnd_type(2).EQ.0) THEN
3118  y(:) = profiles%TE(:)
3119  CALL derivn(nrho,rho,y,dy)
3120  END IF
3121 
3122  te(irho) = y(irho)
3123  dte(irho) = dy(irho) !AF - 25.Set.2014
3124 
3125  if (te(irho).lt.0.0) then
3126  if (irho.eq.1) then
3127  failstring='Error in the temperature equation: on-axis electron temperature is negative, stop'
3128  write(*,*) 'value of on-axis electron temperature is', te(1)
3129  write(*,*) 'check te profile in the file err_tmp'
3130  open (22,file='err_tmp')
3131  rewind(22)
3132  do irho_err=1,nrho
3133  write(22,*) irho_err, solver%y(ndim,irho_err)
3134  enddo
3135  close(22)
3136  ifail=-1
3137  return
3138  else
3139  write(*,*) 'warning, electron temperatuire for irho ',irho,'is negative, set it to the value at the previous irho'
3140  te(irho)=te(irho-1)
3141  dte(irho)=dte(irho-1)
3142  end if
3143  end if
3144 
3145  IF (rho(irho).NE.0.e0_r8) THEN
3146  fun2(irho) = vpr(irho)/rho(irho)* &
3147  ( 1.5e0_r8*nem(irho)*tem(irho)/tau*(vprm(irho)/vpr(irho))**fivethird &
3148  + qe_exp(irho) + qie(irho) - qgi(irho) &
3149  - y(irho) * ( 1.5e0_r8*ne(irho)/tau + qe_imp(irho) + vie(irho) &
3150  - btprime/2.e0_r8/bt*rho(irho)*ne(irho)*dvpr(irho)/vpr(irho) ) )
3151  ELSE
3152  fun2(irho) = ( 1.5e0_r8*nem(irho)*tem(irho)/tau &
3153  + qe_exp(irho) + qie(irho) - qgi(irho) &
3154  - y(irho) * ( 1.5e0_r8*ne(irho)/tau + qe_imp(irho) + vie(irho) &
3155  - btprime/2.e0_r8/bt*ne(irho)*dvpr(irho) ) )
3156  ENDIF
3157  END DO rho_loop13
3158 
3159  CALL integr(nrho,rho,fun2,intfun2) !Integral source
3160 
3161  rho_loop14: DO irho=1,nrho
3162  flux_te_conv(irho) = y(irho)*flux_ne(irho)
3163 
3164  flux_te_cond(irho) = vpr(irho)*g1(irho)*ne(irho) &
3165  * ( y(irho)*vconv_te(irho) - dy(irho)*diff_te(irho) )
3166 
3167  flux_te(irho) = flux_te_conv(irho) + flux_te_cond(irho)
3168 
3169  int_source(irho) = intfun2(irho) &
3170  + y(irho) * 1.5e0_r8*btprime/2.e0_r8/bt*rho(irho)*ne(irho)*vpr(irho)
3171 
3172 
3173 ! +++ If equation is not solved, conductive component of electron heat flux
3174 ! is determined from the integral of sources:
3175  IF (te_bnd_type(2).EQ.0) THEN
3176 
3177  diff_te(irho) = 1.d-6
3178  flux_te(irho) = int_source(irho)
3179 
3180  flux_te_cond(irho) = int_source(irho) - flux_ne(irho)*y(irho)
3181 
3182  IF (vpr(irho)*g1(irho).NE.0.0_r8) &
3183 !dy limit also DY if less than sqrt(epsilon(1.0))
3184  diff_te(irho) = - flux_te_cond(irho) / sign(max(abs(dy(irho)),sqrt(epsilon(1.0))),dy(irho)) / (vpr(irho)*g1(irho)*ne(irho))
3185 !dy further limit diff_ti
3186  if (abs(diff_te(irho)).ge.tite_diff_limit) &
3187  diff_te(irho)=sign(tite_diff_limit,diff_te(irho))
3188  vconv_te(irho) = 0.0_r8
3189  IF (diff_te(irho).LE.1.d-6) THEN
3190  diff_te(irho) = 1.d-6
3191  vconv_te(irho) = (flux_te_cond(irho) / (max(abs(vpr(irho)), 1.e-6)*g1(irho)*ne(irho)) + dy(irho)*diff_te(irho)) / y(irho)
3192  END IF
3193  END IF
3194 
3195 
3196 
3197 ! +++ Return electron temperature and flux:
3198  profiles%TE(irho) = te(irho)
3199  profiles%DTE(irho) = dte(irho) !AF, 25.Sep.2014
3200  profiles%DIFF_TE(irho) = diff_te(irho)
3201  profiles%VCONV_TE(irho) = vconv_te(irho)
3202  profiles%FLUX_TE(irho) = flux_te(irho)
3203  profiles%FLUX_TE_CONV(irho) = flux_te_conv(irho)
3204  profiles%FLUX_TE_COND(irho) = flux_te_cond(irho)
3205  profiles%SOURCE_TE(irho) = qe_exp(irho) + qie(irho) - qgi(irho) &
3206  - (qe_imp(irho) + vie(irho)) * te(irho)
3207  profiles%INT_SOURCE_TE(irho) = int_source(irho)
3208  END DO rho_loop14
3209 
3210  fun1=profiles%SOURCE_TE*vpr
3211  CALL integr2(nrho,rho,fun1,intfun1)
3212  profiles%INT_SOURCE_TE = intfun1
3213 
3214 
3215 ! +++ Deallocate types for interface with numerical solver:
3216  CALL deallocate_numerics(solver, ifail)
3217 
3218 
3219 ! +++ Deallocate types for interface with PLASMA_COLLISIONS:
3220  CALL deallocate_collisionality(collisions, ifail)
3221 
3222  RETURN
3223 
3224 
3225 
3226 END SUBROUTINE temperatures
3227 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3228 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3229 
3230 
3231 
3232 
3233 
3234 
3235 
3236 
3237 ! +++ ROTATION TRANSPORT EQUATIONS +++
3238 
3239 
3240 
3241 
3242 
3243 
3244 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3280 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3281 SUBROUTINE rotation &
3282  (geometry,profiles,transport,sources,evolution,control,ifail,failstring)
3283 
3284 !-------------------------------------------------------!
3285 ! This subroutine solves the momentum transport !
3286 ! equations for ion components fron 1 to NION, !
3287 ! and provides: ion toroidal rotation velocity,ion !
3288 ! angular velocity, ion momentum (total and !
3289 ! individual per ion component), ion momentum flux !
3290 ! (total and individual per ion component), !
3291 !-------------------------------------------------------!
3292 
3293  USE ets_plasma
3294  USE type_solver
3295  USE itm_constants
3296  USE plasma_collisionality
3297 
3298 
3299  USE s4_parameters !AF 12.Oct.2011 - to avoid a few interface blocks
3300 
3301  IMPLICIT NONE
3302 
3303  INTEGER, INTENT (INOUT) :: ifail
3304  CHARACTER(LEN=500) :: failstring
3305 
3306 ! +++ External parameters:
3307  TYPE (magnetic_geometry) :: geometry !contains all geometry quantities
3308  TYPE (plasma_profiles) :: profiles !contains profiles of plasma parameters
3309  TYPE (transport_coefficients) :: transport !contains profiles of trasport coefficients
3310  TYPE (sources_and_sinks) :: sources !contains profiles of sources
3311  TYPE (collisionality) :: collisions !contains all parameters computed by module for plasma collision
3312  TYPE (time_evolution) :: evolution !contains all parameters required by time evolution
3313  TYPE (run_control) :: control !contains all parameters required by run
3314 !evolution and convergence loop
3315 ! +++ Input/Output to numerical solver:
3316  TYPE (numerics) :: solver !contains all I/O quantities to numerics part
3317 
3318 ! +++ Internal parameters:
3319  INTEGER :: irho, nrho !radius index, number of radial points
3320  INTEGER :: iion, nion !index of ion component, number of considered ion components
3321  INTEGER :: zion !second index of ion component
3322 
3323  REAL (R8) :: aion !ion mass number [proton mass]
3324  REAL (R8) :: mion !ion mass [kg]
3325 
3326  REAL (R8) :: bt, btm, btprime !magnetic field, from current time step, [T], previous time steps, [T], time derivative, [T/s]
3327  REAL (R8) :: rho(profiles%nrho) !normalised minor radius, [m]
3328  REAL (R8) :: vpr(profiles%nrho) !V', [m^2]
3329  REAL (R8) :: vprm(profiles%nrho) !V' (at previous time step), [m^2]
3330  REAL (R8) :: g1(profiles%nrho)
3331  REAL (R8) :: g2(profiles%nrho), g2m(profiles%nrho)
3332 
3333  REAL (R8) :: vtor(profiles%nrho) !ion toroidal velocity [m/s]
3334  REAL (R8) :: dvtor(profiles%nrho) !derivative of ion toroidal velocity [m/s] !AF - 25.Sep.2014
3335  REAL (R8) :: vtorm(profiles%nrho) !ion toroidal velocity (at time-1) [m/s]
3336  REAL (R8) :: dvtorm(profiles%nrho) !derivative of ion toroidal velocity at previous time !AF - 25.Sep.2014
3337  REAL (R8) :: wtor(profiles%nrho) !ion angular velocity [1/s]
3338  REAL (R8) :: mtor(profiles%nrho) !ion toroidal momentum [kg*m/s]
3339  REAL (R8) :: mtor_tot(profiles%nrho) !total ion toroidal momentum [kg*m/s]
3340 
3341  REAL (R8) :: diff(profiles%nrho),vconv(profiles%nrho) !velocity diffusion coefficient, [m^2/s], and pinch velocity, [m/s]
3342 
3343  REAL (R8) :: flux_mtor(profiles%nrho) !ion toroidal momentum flux [kg*m^2/s^2]
3344  REAL (R8) :: flux_mtor_conv(profiles%nrho) !convective contribution to ion toroidal momentum flux [kg*m^2/s^2]
3345  REAL (R8) :: flux_mtor_cond(profiles%nrho) !conductive contribution ion toroidal momentum flux [kg*m^2/s^2]
3346  REAL (R8) :: flux_mtor_tot(profiles%nrho) !total ion toroidal momentum flux [kg*m^2/s^2]
3347  REAL (R8) :: int_source(profiles%nrho) !integral of momentum sources [kg*m^2/s^2]
3348  REAL (R8) :: flux_ni(profiles%nrho) !ion particle flux [1/s]
3349 
3350  REAL (R8) :: vtor_bnd(2,3) !boundary condition, value [depends on VTOR_BND_TYPE]
3351  INTEGER :: vtor_bnd_type(2) !boundary condition, type
3352 
3353  REAL (R8) :: ti(profiles%nrho) !ion temperature, [eV]
3354  REAL (R8) :: ni(profiles%nrho) !ion density [m^-3]
3355  REAL (R8) :: nim(profiles%nrho) !ion density (at time-1) [m^-3]
3356  REAL (R8) :: dni(profiles%nrho), dnim(profiles%nrho) !AF - 25.Sep.2014
3357 
3358  REAL (R8) :: ui_exp(profiles%nrho) !external explicit ion momentum source [kg/m/s^2]
3359  REAL (R8) :: ui_imp(profiles%nrho) !external implicit ion momentum source [kg/m^2/s]
3360  REAL (R8) :: uzi(profiles%nrho) !ion momentum source due to interaction with ions [kg/m/s^2]
3361  REAL (R8) :: wzi(profiles%nrho) !friction between different ion components (normalized frequency) [kg/m^2/s]
3362  REAL (R8) :: wii(profiles%nrho,profiles%nion) !momentum exchange between different ion components, [kg/(m^2*s)]
3363 
3364  REAL (R8) :: amix, tau !mixing factor, time step, [s]
3365 
3366  INTEGER :: flag !flag for equation: 0 - interpretative (not solved), 1 - predictive (solved)
3367  INTEGER :: idim, ndim !index and total number of equations to be solved
3368  INTEGER :: solver_type !specifies the option for numerical solution
3369  REAL (R8) :: y(profiles%nrho), ym(profiles%nrho) !function at the current amd previous time steps
3370  REAL (R8) :: dy(profiles%nrho) !derivative of function
3371  REAL (R8) :: dym(profiles%nrho) !derivative of function at previous time step !AF - 25.Sep.2014
3372  REAL (R8) :: a(profiles%nrho), b(profiles%nrho) !coefficients for numerical solver
3373  REAL (R8) :: c(profiles%nrho), d(profiles%nrho) !coefficients for numerical solver
3374  REAL (R8) :: e(profiles%nrho), f(profiles%nrho) !coefficients for numerical solver
3375  REAL (R8) :: g(profiles%nrho), h !coefficients for numerical solver
3376  REAL (R8) :: v(2), u(2), w(2) !boundary conditions for numerical solver
3377 
3378  REAL (R8) :: fun1(profiles%nrho), intfun1(profiles%nrho)
3379 
3380 ! +++ Set up dimensions:
3381  nrho = profiles%NRHO
3382  nion = profiles%NION
3383  ndim = nion
3384 
3385 
3386 ! +++ Allocate types for interface with PLASMA_COLLISIONS:
3387  CALL allocate_collisionality(nrho, nion, collisions, ifail)
3388 
3389 
3390 ! +++ Allocate types for interface with numerical solver:
3391  CALL allocate_numerics(ndim, nrho, solver, ifail)
3392 
3393 
3394 
3395 ! +++ Set up local variables:
3396  amix = control%AMIX
3397  tau = control%TAU
3398  solver_type = control%SOLVER_TYPE
3399 
3400  bt = geometry%BGEO
3401  btm = evolution%BTM
3402  btprime = (bt-btm)/tau
3403 
3404  rho_loop1: DO irho=1,nrho
3405  rho(irho) = geometry%RHO(irho)
3406  vpr(irho) = geometry%VPR(irho)
3407  vprm(irho) = evolution%VPRM(irho)
3408  g1(irho) = geometry%G1(irho)
3409  g2(irho) = geometry%G2(irho)
3410  g2m(irho) = evolution%G2M(irho)
3411 
3412  mtor_tot(irho) = 0.e0_r8
3413  flux_mtor_tot(irho) = 0.e0_r8
3414  END DO rho_loop1
3415 
3416 
3417 
3418 ! +++ Exchange terms (defined on previous iteration):
3419 ! CALL PLASMA_COLLISIONS (GEOMETRY,PROFILES,COLLISIONS,ifail)
3420 
3421 
3422 
3423 ! +++ Numeric types common for all ions:
3424  solver%TYPE = solver_type
3425  solver%NDIM = ndim
3426  solver%NRHO = nrho
3427  solver%AMIX = amix
3428  solver%DERIVATIVE_FLAG(1) = 0
3429 
3430 
3431 
3432 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3433 ! solution of rotation transport equation for
3434 ! individual ion species
3435 
3436  ion_loop1: DO iion=1,nion
3437 
3438  idim = iion
3439 
3440 ! +++ Set equation to 'predictive' and all coefficients to zero:
3441  flag = 1
3442  y(:) = 0.0e0_r8
3443  dy(:) = 0.0e0_r8
3444  ym(:) = 0.0e0_r8
3445  dym(:) = 0.0e0_r8 !AF - 25.Sep.2014
3446  a(:) = 0.0e0_r8
3447  b(:) = 0.0e0_r8
3448  c(:) = 0.0e0_r8
3449  d(:) = 0.0e0_r8
3450  e(:) = 0.0e0_r8
3451  f(:) = 0.0e0_r8
3452  g(:) = 0.0e0_r8
3453  h = 0.0e0_r8
3454  v(:) = 0.0e0_r8
3455  u(:) = 0.0e0_r8
3456  w(:) = 0.0e0_r8
3457 
3458 
3459 
3460 ! +++ Set up ion mass:
3461  mion = profiles%MION(iion)!*MP
3462 
3463 
3464 
3465 ! +++ Set up boundary conditions for particular ion type:
3466  vtor_bnd_type(2) = profiles%VTOR_BND_TYPE(iion)
3467  vtor_bnd(2,1) = profiles%VTOR_BND(1,iion)
3468  vtor_bnd(2,2) = profiles%VTOR_BND(2,iion)
3469  vtor_bnd(2,3) = profiles%VTOR_BND(3,iion)
3470 
3471 
3472 
3473 ! +++ Set up local variables for particular ion type:
3474  rho_loop2: DO irho = 1,nrho
3475  vtor(irho) = profiles%VTOR(irho,iion)
3476  dvtor(irho) = profiles%DVTOR(irho,iion) !AF - 25.Sep.2014
3477  vtorm(irho) = evolution%VTORM(irho,iion)
3478  dvtorm(irho) = evolution%DVTORM(irho,iion) !AF - 25.Sep.2014
3479  wtor(irho) = profiles%WTOR(irho,iion)
3480  ni(irho) = profiles%NI(irho,iion)
3481  dni(irho) = profiles%DNI(irho,iion) !AF - 25.Sep.2014
3482  nim(irho) = evolution%NIM(irho,iion)
3483  dnim(irho) = evolution%DNIM(irho,iion) !AF - 25.Sep.2014
3484  mtor(irho) = profiles%MTOR(irho,iion)
3485 
3486  flux_mtor(irho) = profiles%FLUX_MTOR(irho,iion)
3487  flux_ni(irho) = profiles%FLUX_NI(irho,iion)
3488 
3489  diff(irho) = transport%DIFF_VTOR(irho,iion)
3490  vconv(irho) = transport%VCONV_VTOR(irho,iion)
3491 
3492  ui_exp(irho) = sources%UI_EXP(irho,iion)
3493  ui_imp(irho) = sources%UI_IMP(irho,iion)
3494 
3495  wzi(irho) = collisions%WZI(irho,iion)
3496  uzi(irho) = collisions%UZI(irho,iion)
3497 
3498  ion_loop2: DO zion = 1,nion
3499 !!DPC was WII(IRHO,IION) = COLLISIONS%WII(IRHO,IION,ZION) !!! check if the following is what was intended
3500  wii(irho,zion) = collisions%WII(irho,iion,zion)
3501  END DO ion_loop2
3502 
3503  END DO rho_loop2
3504 
3505 
3506 ! +++ Coefficients for rotation transport equation in form:
3507 !
3508 ! (A*Y-B*Y(t-1))/H + 1/C * (-D*Y' + E*Y) = F - G*Y
3509 
3510  rho_loop3: DO irho=1,nrho
3511  y(irho) = vtor(irho)
3512  dy(irho) = dvtor(irho) !AF - 25.Sep.2014
3513  ym(irho) = vtorm(irho)
3514  dym(irho) = dvtorm(irho) !AF - 25.Sep.2014
3515 
3516  a(irho) = vpr(irho)*g2(irho)*ni(irho)*mion
3517  b(irho) = vprm(irho)*g2m(irho)*nim(irho)*mion
3518  c(irho) = 1.e0_r8
3519  d(irho) = vpr(irho)*g1(irho)*ni(irho)*mion*diff(irho)*g2(irho) !AF, 14.May.2011 - multipication by G2, which in analytics is 1
3520  e(irho) = (vpr(irho)*g1(irho)*ni(irho)*vconv(irho) &
3521  + flux_ni(irho) &
3522  - btprime/2.e0_r8/bt*rho(irho)*ni(irho)*vpr(irho)) &
3523  * g2(irho)*mion
3524  f(irho) = vpr(irho)*(ui_exp(irho) + uzi(irho))
3525  g(irho) = vpr(irho)*(ui_imp(irho) + wzi(irho))
3526  END DO rho_loop3
3527 
3528  h = tau
3529 
3530 
3531 
3532 ! +++ Boundary conditions for numerical solver in form:
3533 !
3534 ! V*Y' + U*Y =W
3535 
3536 ! +++ On axis
3537 ! dVtor,i/drho(rho=0)=0: ! AF - 24.Jun.2010, replaces "! dTi/drho(rho=0)=0:"
3538  IF (solver_type.NE.4) THEN !AF 11.Oct.2011
3539  v(1) = 1.e0_r8
3540  u(1) = 0.e0_r8
3541  ELSE !AF 11.Oct.2011 - Zero flux instead of zero gradient at the axis for solver 4
3542  IF (diff(1).GT.1.0e-6) THEN !AF 19.Mar.2012 - To avoid problems with the axis boundary condition
3543  v(1) = -diff(1)*ni(1)
3544  ELSE
3545  v(1) = -1.0e-6*ni(1)
3546  ENDIF
3547  u(1) = vconv(1)*ni(1) & !AF 19.Jul.2011 - Roman found this one
3548  +local_flux_ni_s4(iion) !AF 19.Jul.2011 - Roman found this one
3549  END IF !AF 11.Oct.2011
3550  w(1) = 0.e0_r8
3551 
3552 ! +++ At the edge:
3553 ! FIXED Vtor,i
3554  IF(vtor_bnd_type(2).EQ.1) THEN
3555  v(2) = 0.e0_r8
3556  u(2) = 1.e0_r8
3557  w(2) = vtor_bnd(2,1)
3558  ENDIF
3559 
3560 ! FIXED grad_Vtor,i
3561  IF(vtor_bnd_type(2).EQ.2) THEN
3562  v(2) = 1.e0_r8
3563  u(2) = 0.e0_r8
3564  w(2) = -vtor_bnd(2,1)
3565  ENDIF
3566 
3567 ! FIXED L_Vtor,i
3568  IF(vtor_bnd_type(2).EQ.3) THEN
3569  v(2) = 1.e0_r8
3570  u(2) = 1.e0_r8/vtor_bnd(2,1)
3571  w(2) = 0.e0_r8
3572  ENDIF
3573 
3574 ! FIXED Flux_Mtor,i
3575  IF(vtor_bnd_type(2).EQ.4) THEN
3576  v(2) = -vpr(nrho)*g1(nrho)*g2(nrho)*diff(nrho)*ni(nrho)*mion
3577  u(2) = vpr(nrho)*g1(nrho)*g2(nrho)*vconv(nrho)*ni(nrho)*mion &
3578  +g2(nrho)*flux_ni(nrho)*mion
3579  w(2) = vtor_bnd(2,1)
3580  ENDIF
3581 
3582 ! Generic boundary condition
3583  IF(vtor_bnd_type(2).EQ.5) THEN
3584  v(2) = vtor_bnd(2,1)
3585  u(2) = vtor_bnd(2,2)
3586  w(2) = vtor_bnd(2,3)
3587  ENDIF
3588 
3589 
3590 
3591 ! +++ Rotation equation is not solved:
3592  IF(vtor_bnd_type(2).EQ.0) THEN
3593 
3594  CALL derivn(nrho,rho,y,dy)
3595 
3596  flag = 0
3597 
3598  rho_loop4: DO irho=1,nrho
3599  a(irho) = 1.0e0_r8
3600  b(irho) = 1.0e0_r8
3601  c(irho) = 1.0e0_r8
3602  d(irho) = 0.0e0_r8
3603  e(irho) = 0.0e0_r8
3604  f(irho) = 0.0e0_r8
3605  g(irho) = 0.0e0_r8
3606  END DO rho_loop4
3607 
3608  v(2) = 0.0e0_r8
3609  u(2) = 1.0e0_r8
3610  w(2) = y(nrho)
3611  END IF
3612 
3613 
3614 
3615 ! +++ Defining coefficients for numerical solver:
3616  solver%EQ_FLAG(idim) = flag
3617 
3618  rho_loop5: DO irho=1,nrho
3619  solver%RHO(irho) = rho(irho)
3620 
3621  solver%Y(idim,irho) = y(irho)
3622  solver%DY(idim,irho) = dy(irho)
3623  solver%YM(idim,irho) = ym(irho)
3624 
3625  solver%A(idim,irho) = a(irho)
3626  solver%B(idim,irho) = b(irho)
3627  solver%C(idim,irho) = c(irho)
3628  solver%D(idim,irho) = d(irho)
3629  solver%E(idim,irho) = e(irho)
3630  solver%F(idim,irho) = f(irho)
3631  solver%G(idim,irho) = g(irho)
3632 
3633  solver%CM1(idim,idim,irho) = 0.e0_r8
3634 
3635  ion_loop3: DO zion = 1,nion
3636  solver%CM1(iion,zion,irho)= wii(irho,zion)
3637  END DO ion_loop3
3638 
3639  END DO rho_loop5
3640 
3641  solver%H = h
3642 
3643  solver%V(idim,1) = v(1)
3644  solver%U(idim,1) = u(1)
3645  solver%W(idim,1) = w(1)
3646  solver%V(idim,2) = v(2)
3647  solver%U(idim,2) = u(2)
3648  solver%W(idim,2) = w(2)
3649 
3650 
3651 
3652  END DO ion_loop1
3653 
3654 
3655 ! +++ Solution of rotation transport equation:
3656  CALL solution_interface(solver, ifail)
3657 
3658  ! dy check for nans in solution, return if true)
3659  IF (any( isnan(solver%y))) then
3660  failstring='Error in the vtor equation: nans in the solution, stop'
3661  ifail=-1
3662  return
3663  end if
3664 
3665 
3666  ion_loop4: DO iion = 1,nion
3667 
3668 
3669  idim = iion
3670 
3671 ! dpc 2011-08-11: I think we need most of the following
3672  mion = profiles%MION(iion)!*MP
3673  vtor_bnd_type(2) = profiles%VTOR_BND_TYPE(iion)
3674  vtor_bnd(2,1) = profiles%VTOR_BND(1,iion)
3675  vtor_bnd(2,2) = profiles%VTOR_BND(2,iion)
3676  vtor_bnd(2,3) = profiles%VTOR_BND(3,iion)
3677 ! dpc end
3678 
3679  rho_loop6: DO irho=1,nrho
3680 
3681 ! dpc 2011-08-11: I think we need most of the following
3682  vtor(irho) = profiles%VTOR(irho,iion)
3683  vtorm(irho) = evolution%VTORM(irho,iion)
3684  wtor(irho) = profiles%WTOR(irho,iion)
3685  ni(irho) = profiles%NI(irho,iion)
3686  nim(irho) = evolution%NIM(irho,iion)
3687  mtor(irho) = profiles%MTOR(irho,iion)
3688  flux_mtor(irho) = profiles%FLUX_MTOR(irho,iion)
3689  flux_ni(irho) = profiles%FLUX_NI(irho,iion)
3690  diff(irho) = transport%DIFF_VTOR(irho,iion)
3691  vconv(irho) = transport%VCONV_VTOR(irho,iion)
3692  ui_exp(irho) = sources%UI_EXP(irho,iion)
3693  ui_imp(irho) = sources%UI_IMP(irho,iion)
3694  wzi(irho) = collisions%WZI(irho,iion)
3695  uzi(irho) = collisions%UZI(irho,iion)
3696 ! dpc end
3697 
3698 ! +++ New solution and its derivative:
3699  y(irho) = solver%Y(idim,irho)
3700  dy(irho) = solver%DY(idim,irho)
3701 
3702  IF(vtor_bnd_type(2).EQ.0) THEN
3703  y(:) = profiles%VTOR(:,iion)
3704  CALL derivn(nrho,rho,y,dy)
3705  END IF
3706 
3707 ! +++ New rotation velocity and momentum flux:
3708  vtor(irho) = y(irho)
3709  dvtor(irho) = dy(irho) !AF - 25.Set.2014
3710 !dy 2017-10-06 WTOR(IRHO) = Y(IRHO)/G2(NRHO)
3711  wtor(irho) = y(irho)/g2(irho)
3712 
3713  mtor(irho) = g2(irho)*ni(irho)*mion*y(irho)
3714  mtor_tot(irho) = mtor_tot(irho) + mtor(irho)
3715 
3716  IF (rho(irho).NE.0.e0_r8) THEN
3717  fun1(irho) = vpr(irho)/rho(irho) * ( ui_exp(irho) + uzi(irho) &
3718  + ( wzi(irho) + g2(irho)*mion*ni(irho)/tau- ui_imp(irho)) * y(irho) )
3719  ELSE
3720  fun1(irho) = ( ui_exp(irho) + uzi(irho)+ ( wzi(irho) &
3721  + g2(irho)*mion*ni(irho)/tau- ui_imp(irho)) * y(irho) )
3722  ENDIF
3723  END DO rho_loop6
3724 
3725  CALL integr(nrho,rho,fun1,intfun1) !Integral source
3726 
3727  rho_loop7: DO irho=1,nrho
3728  flux_mtor_conv(irho) = g2(irho)*mion*flux_ni(irho)*y(irho)
3729 
3730  flux_mtor_cond(irho) = vpr(irho)*g1(irho)*g2(irho)*mion*ni(irho) &
3731  * ( y(irho)*vconv(irho) - dy(irho)*diff(irho) )
3732 
3733  flux_mtor(irho) = flux_mtor_conv(irho) + flux_mtor_cond(irho)
3734 
3735  int_source(irho) = intfun1(irho) &
3736  + vpr(irho)*g2(irho)*btprime/2.e0_r8/bt*rho(irho)*mion*ni(irho)*y(irho)
3737 
3738  IF (vtor_bnd_type(2).EQ.0) THEN !if equation is not solved, conductive component of electron heat flux is determined from the integral of sources
3739 
3740  diff(irho) = 1.d-6
3741  flux_mtor(irho) = int_source(irho)
3742 
3743  flux_mtor_cond(irho)= int_source(irho) - g2(irho)*mion*flux_ni(irho)*y(irho)
3744 
3745 
3746  IF ((vpr(irho)*g1(irho)*g2(irho)*mion*ni(irho).NE.0.0_r8).and.(dy(irho).ne.0.0_r8)) &
3747  diff(irho) = - flux_mtor_cond(irho) / dy(irho) / (vpr(irho)*g1(irho)*g2(irho)*mion*ni(irho))
3748  vconv(irho) = 0.0_r8
3749  IF (diff(irho).LE.1.d-6) THEN
3750  diff(irho) = 1.d-6
3751  vconv(irho) = (flux_mtor_cond(irho) / (max(abs(vpr(irho)), 1.e-6)*g1(irho)*g2(irho)*mion*ni(irho)) + dy(irho)*diff(irho)) / max(abs(y(irho)), 1.e-6)
3752  END IF
3753 
3754  END IF
3755 
3756  flux_mtor_tot(irho) = flux_mtor_tot(irho) + flux_mtor(irho)
3757 
3758 
3759 
3760 ! +++ Return new profiles to the work flow:
3761  profiles%VTOR(irho,iion) = vtor(irho)
3762  profiles%DVTOR(irho,iion) = dvtor(irho) !AF, 25.Sep.2014
3763 !dy 2019-08-31 make sure profiles from db are written to the coreprof
3764 ! way around, address afterwords for consistent fix
3765 ! PROFILES%WTOR(IRHO,IION) = WTOR(IRHO)
3766  profiles%MTOR(irho,iion) = mtor(irho)
3767  profiles%DIFF_VTOR(irho,iion) = diff(irho)
3768  profiles%VCONV_VTOR(irho,iion) = vconv(irho)
3769  profiles%FLUX_MTOR(irho,iion) = flux_mtor(irho)
3770  profiles%FLUX_MTOR_CONV(irho,iion) = flux_mtor_conv(irho)
3771  profiles%FLUX_MTOR_COND(irho,iion) = flux_mtor_cond(irho)
3772  profiles%SOURCE_MTOR(irho,iion) = ui_exp(irho) + uzi(irho) - (ui_imp(irho) + wzi(irho)) * vtor(irho)
3773  profiles%INT_SOURCE_MTOR(irho,iion) = int_source(irho)
3774  END DO rho_loop7
3775 
3776  fun1=profiles%SOURCE_MTOR(:,iion)*vpr
3777  CALL integr2(nrho,rho,fun1,intfun1)
3778  profiles%INT_SOURCE_MTOR(:,iion) = intfun1
3779 
3780 
3781  END DO ion_loop4
3782 
3783 
3784 
3785  rho_loop8: DO irho=1,nrho
3786  profiles%MTOR_TOT(irho) = mtor_tot(irho)
3787  profiles%FLUX_MTOR_TOT(irho) = flux_mtor_tot(irho)
3788  END DO rho_loop8
3789 
3790 
3791 
3792 ! +++ Deallocate types for interface with numerical solver:
3793  CALL deallocate_numerics(solver, ifail)
3794 
3795 
3796 
3797 ! +++ Deallocate types for interface with PLASMA_COLLISIONS:
3798  CALL deallocate_collisionality(collisions, ifail)
3799 
3800 
3801 
3802  RETURN
3803 
3804 
3805 
3806 END SUBROUTINE rotation
3807 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3808 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3809 
3810 
3811 
3812 
3813 
3814 
3815 
3816 
3817 
3818 
3819 
3820 
3821 
3822 
3823 
3824 
3825 
3826 
3827 
3828 
3829 
3830 
3831 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3832 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3833 !-------------------------------------------------------!
3834 ! !
3835 !______________ MATHEMATICAL SUBROUTINES: _____________!
3836 ! !
3837 !-------------------------------------------------------!
3838 ! These subroutines have been extracted from RITM code, !
3839 ! and consist of derivation and integration routines !
3840 !-------------------------------------------------------!
3841 
3842 
3843 
3844 
3845 
3846 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3855 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3856 SUBROUTINE derivn(N,X,Y,DY1)
3857 
3858 !-------------------------------------------------------!
3859 ! These subroutines calculate first and second !
3860 ! derivatives, DY1 and DY2, of function Y respect !
3861 ! to argument X !
3862 !-------------------------------------------------------!
3863 
3864  use itm_types
3865  IMPLICIT NONE
3866 
3867  INTEGER :: n ! number of radial points (input)
3868  INTEGER :: i
3869 
3870  REAL (R8) :: x(n), & ! argument array (input)
3871  y(n), & ! function array (input)
3872  dy1(n) ! function derivative array (output)
3873  REAL (R8) :: h(n),dy2(n)
3874 
3875  REAL (R8) :: ddy !AF 6.Oct.2011
3876 
3877  DO i=1,n-1
3878  h(i)=x(i+1)-x(i)
3879  END DO
3880 
3881  DO i=2,n-1
3882  dy1(i)=((y(i+1)-y(i))*h(i-1)/h(i)+(y(i)-y(i-1))*h(i)/h(i-1)) &
3883  /(h(i)+h(i-1))
3884 ! DY2(I)=2.e0_R8*((Y(I-1)-Y(I))/H(I-1)+(Y(I+1)-Y(I))/H(I)) & !AF 6.Oct.2011
3885 ! /(H(I)+H(I-1))
3886  END DO
3887 
3888 ! DY1(1)=DY1(2)-DY2(2)*H(1) !AF 6.Oct.2011
3889 ! DY1(N)=DY1(N-1)+DY2(N-1)*H(N-1) !AF 6.Oct.2011
3890 
3891  ddy = 2.e0_r8*((y(1)-y(2))/h(1)+(y(3)-y(2))/h(2))/(h(2)+h(1))
3892  dy1(1) = dy1(2)-ddy*h(1)
3893  ddy = 2.e0_r8*((y(n-2)-y(n-1))/h(n-2)+(y(n)-y(n-1))/h(n-1))/(h(n-1)+h(n-2))
3894  dy1(n) = dy1(n-1)+ddy*h(n-1)
3895 
3896  RETURN
3897 END SUBROUTINE derivn
3898 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3899 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3900 
3901 
3902 
3903 
3904 
3905 
3906 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3915 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3916 SUBROUTINE integr(N,X,Y,INTY)
3917 !-------------------------------------------------------!
3918 ! This subroutine calculates integral of function !
3919 ! Y(X)*X from X=0 until X=X(N) !
3920 !-------------------------------------------------------!
3921 
3922  use itm_types
3923 
3924  IMPLICIT NONE
3925 
3926  INTEGER :: n ! number of radial points (input)
3927  INTEGER :: i
3928 
3929  REAL (R8) :: x(n), & ! argument array (input)
3930  y(n), & ! function array (input)
3931  inty(n) ! function integral array (output)
3932 
3933  inty(1)=y(1)*x(1)**2/2.e0_r8
3934  DO i=2,n
3935  inty(i)=inty(i-1)+(y(i-1)*x(i-1)+y(i)*x(i))*(x(i)-x(i-1))/2.e0_r8
3936  END DO
3937 
3938  RETURN
3939 END SUBROUTINE integr
3940 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3941 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3942 
3943 
3944 SUBROUTINE integr2(N,X,Y,INTY) !AF 11.Oct.2011 - assumes that Y is zero f0r X.eq.0, just as INTEGR does too...
3945 !-------------------------------------------------------!
3946 ! This subroutine calculates integral of function !
3947 ! Y(X) from X=0 until X=X(N) !
3948 !-------------------------------------------------------!
3949 
3950  use itm_types
3951 
3952  IMPLICIT NONE
3953 
3954  INTEGER :: n ! number of radial points (input)
3955  INTEGER :: i
3956 
3957  REAL (R8) :: x(n), & ! argument array (input)
3958  y(n), & ! function array (input)
3959  inty(n) ! function integral array (output)
3960 
3961  inty(1)=y(1)*x(1)/2.e0_r8
3962  DO i=2,n
3963  inty(i)=inty(i-1)+(y(i-1)+y(i))*(x(i)-x(i-1))/2.e0_r8
3964  END DO
3965 
3966  RETURN
3967 
3968 END SUBROUTINE integr2
3969 
3970 
3971 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3977 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
3978 SUBROUTINE f_axis(n, r, f)
3979  !-------------------------------------------------------!
3980  ! !
3981  ! This subroutine finds !
3982  ! f(r_1=0) from f(r_2), f(r_3) and f(r_4) !
3983  ! !
3984  !-------------------------------------------------------!
3985 
3986  IMPLICIT NONE
3987 
3988  INTEGER n, i
3989  REAL *8 h(n), r(n), f(n)
3990 
3991 
3992  DO i=1,3
3993  h(i)=r(i+1)-r(i)
3994  END DO
3995 
3996  f(1) = ((f(2)*r(4)/h(2)+f(4)*r(2)/h(3))*r(3) &
3997  -f(3)*(r(2)/h(2)+r(4)/h(3))*r(2)*r(4)/r(3)) &
3998  /(r(4)-r(2))
3999 
4000 
4001  RETURN
4002 
4003 
4004 END SUBROUTINE f_axis
4005 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
4006 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
4007 
4008 
4009 
4010 
4011 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
4017 ! + + + + + + + + + + + + + + + + + + + + + + + + + + + +
4018 
4019 SUBROUTINE f_par_axis(n, r, f)
4020  use itm_types
4021  implicit none
4022  integer n
4023  real(R8) r(n), f(n), d1, d2
4024  if(n.lt.3) stop 'n too small in F_par_AXIS'
4025  d1=r(2)-r(1)
4026  if(d1.le.0.0_r8) stop 'd1 <= 0 in F_par_AXIS'
4027  d2=r(3)-r(2)
4028  if(d2.le.0.0_r8) stop 'd2 <= 0 in F_par_AXIS'
4029  f(1)=f(2)-d1**2*(f(3)-f(2))/(d2**2+2*d1*d2)
4030  return
4031 end SUBROUTINE f_par_axis
solution_interface
subroutine solution_interface(SOLVER, ifail)
INTERFACE TO NUMERICAL SOLVER.
Definition: solution_interface.f90:12
impurity
IMPURITY.
Definition: impurity.F90:8
profiles
subroutine profiles(p0, rbphi, dp0, drbphi, a)
Definition: profiles.f90:1
solver
Definition: solver.f90:1
type_solver
Definition: type_solver.f90:1
electron_density
subroutine electron_density(GEOMETRY, PROFILES, TRANSPORT, SOURCES, EVOLUTION, CONTROL, HYPER_DIFF, ifail, failstring)
Definition: main_plasma.f90:1342
derivn
subroutine derivn(N, X, Y, DY1)
These subroutines calculate first and second derivatives, DY1 and DY2, of function Y respect to argum...
Definition: main_plasma.f90:3856
ets_plasma::sources_and_sinks
Definition: ets_plasma.f90:294
main_plasma
subroutine main_plasma
Main plasma.
Definition: main_plasma.f90:16
f_axis
subroutine f_axis(n, r, f)
This subroutine finds f(r_1=0) from f(r_2), f(r_3) and f(r_4)
Definition: main_plasma.f90:3978
temperatures
subroutine temperatures
HEAT TRANSPORT EQUATIONS.
Definition: main_plasma.f90:2290
ets_plasma::impurity_profiles
Definition: ets_plasma.f90:373
ets_plasma::time_evolution
Definition: ets_plasma.f90:395
integr2
subroutine integr2(N, X, Y, INTY)
Definition: main_plasma.f90:3944
type_solver::allocate_numerics
subroutine allocate_numerics(NDIM, NRHO, SOLVER, ifail)
Definition: type_solver.f90:66
ets_plasma::diagnostic
Definition: ets_plasma.f90:463
flux
subroutine flux(psitok, rk, zk, nk)
Definition: EQ1_m.f:786
plasma_collisionality
This routine calculates the collision frquencies and various exchange terms determined by collisions...
Definition: plasma_collisionality.f90:10
ion_density
subroutine ion_density
PARTICLE TRANSPORT EQUATIONS.
Definition: main_plasma.f90:811
ets_plasma::plasma_profiles
Definition: ets_plasma.f90:80
rotation
subroutine rotation(GEOMETRY, PROFILES, TRANSPORT, SOURCES, EVOLUTION, CONTROL, ifail, failstring)
ROTATION TRANSPORT EQUATIONS.
Definition: main_plasma.f90:3281
current
subroutine current(GEOMETRY, PROFILES, TRANSPORT, SOURCES, EVOLUTION, CONTROL, j_boun, ifail, failstring)
CURRENT TRANSPORT EQUATION.
Definition: main_plasma.f90:205
helena21_mod::fdia
real(r8) function fdia(psi_n)
Definition: helena21_mod.f90:139
integr
subroutine integr(N, X, Y, INTY)
This subroutine calculates integral of function Y(X)*X from X=0 until X=X(N)
Definition: main_plasma.f90:3916
type_solver::deallocate_numerics
subroutine deallocate_numerics(SOLVER, ifail)
Deallocate plasma profiles needed by the transport solver.
Definition: type_solver.f90:140
ets_plasma::run_control
Definition: ets_plasma.f90:426
plasma_collisionality::plasma_collisions
subroutine plasma_collisions(GEOMETRY, PROFILES, COLLISIONS, ifail)
Definition: plasma_collisionality.f90:17
quasi_neutrality
subroutine quasi_neutrality(GEOMETRY, PROFILES, IMPURITY, CONTROL, ifail, failstring)
QUASI NEUTRALITY.
Definition: main_plasma.f90:1838
ets_plasma::deallocate_collisionality
subroutine deallocate_collisionality(COLLISIONS, ifail)
Deallocate plasma profiles needed by the transport solver.
Definition: ets_plasma.f90:1519
type_solver::numerics
Definition: type_solver.f90:22
ets_plasma::collisionality
Definition: ets_plasma.f90:343
ets_plasma::allocate_collisionality
subroutine allocate_collisionality(NRHO, NION, COLLISIONS, ifail)
Allocate profiles of sources needed by the transport solver???
Definition: ets_plasma.f90:1456
ets_plasma
The module declares types of variables used in ETS (transport code)
Definition: ets_plasma.f90:8
ets_plasma::magnetic_geometry
Definition: ets_plasma.f90:44
evolution
subroutine evolution(T, R_in, R_out, El, Tr_l, Tr_U, Ip)
Definition: ets_workflow.F90:1871
f_par_axis
subroutine f_par_axis(n, r, f)
This subroutine finds f(r_1=0) from f(r_2), f(r_3) d/dr f(r_1)=0.
Definition: main_plasma.f90:4019
ets_plasma::transport_coefficients
Definition: ets_plasma.f90:237
s4_parameters
Definition: s4_parameters.f90:1