ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
neo.f90
Go to the documentation of this file.
1 
14 
15 MODULE spec_kind_mod
16  !-------------------------------------------------------------------------------
17  !For REAL variables:
18  ! Use p=6, r=35 for single precision on 32 bit machine
19  ! Use p=12,r=100 for double precision on 32 bit machine
20  ! for single precision on 64 bit machine
21  !Parameters for SELECTED_REAL_KIND:
22  ! p -number of digits
23  ! r -range of exponent from -r to r
24  !-------------------------------------------------------------------------------
25  implicit none
26 
27  INTEGER :: iecrit=0,iindice, kcrit
28 
29  !INTEGER, PARAMETER :: rspec=kind(1.0D0),ispec=kind(iecrit)
30  INTEGER, PARAMETER ::rspec = SELECTED_REAL_KIND(p=12,r=100),ispec = 8
31 
32 END MODULE spec_kind_mod
33 
77 
78 MODULE nclass_itm
79  !-------------------------------------------------------------------------------
80  !NCLASS-Calculates NeoCLASSical transport properties
81  !
82  !NCLASS_ITM is an F90 module to calculate neoclassical transport properties in
83  ! for an axisymmetric toroidal plasma based on CRONOS coupling
84  !
85  !References:
86  !
87  ! W.A.Houlberg 1/2002
88  !
89  !Contains PUBLIC routines:
90  !
91  ! NCLASS -neoclassical properties on a single surface
92  !
93  !Contains PRIVATE routines:
94  !
95  ! NCLASS_FLOW -flow velocities within a surface
96  ! -called from NCLASS
97  ! NCLASS_INIT -initialize arrays
98  ! -called from NCLASS
99  ! NCLASS_K -velocity-dependent viscositite
100  ! -called from NCLASS_MU
101  ! NCLASS_MN -friction coefficients
102  ! -called from NCLASS
103  ! NCLASS_MU -fluid viscosities
104  ! -called from NCLASS
105  ! NCLASS_NU -collision frequencies
106  ! -called from NCLASS_K
107  ! NCLASS_TAU -collision times
108  ! -called from NCLASS
109  ! NCLASS_BACKSUB -LU matrix back substitution
110  ! -called from NCLASS_FLOW
111  ! NCLASS_DECOMP -LU matrix decomposition
112  ! -called from NCLASS_FLOW
113  ! NCLASS_ERF -error function
114  ! -called from NCLASS_NU
115  !
116  !Comments: cvs
117  !
118  !-------------------------------------------------------------------------------
119  USE spec_kind_mod
120  IMPLICIT NONE
121 
122  !-------------------------------------------------------------------------------
123  ! Private procedures
124  !-------------------------------------------------------------------------------
125  PRIVATE :: &
126  & nclass_flow, &
127  & nclass_init, &
128  & nclass_k, &
129  & nclass_mn, &
130  & nclass_mu, &
131  & nclass_nu, &
132  & nclass_tau, &
133  & nclass_backsub, &
134  & nclass_decomp, &
135  & nclass_erf, &
136  & normalisation
137 
138  !-------------------------------------------------------------------------------
139  ! Private data
140  !-------------------------------------------------------------------------------
141  !Logical switches
142  LOGICAL :: &
143  & l_banana_nc,l_pfirsch_nc,l_potato_nc
144 
145  !Options
146  INTEGER, PRIVATE :: &
147  & k_order_nc
148 
149  !Constants
150  REAL, PRIVATE, SAVE :: &
151  & cden_nc,cpotb_nc,cpotl_nc,errsum
152 
153  !Terms summed over species
154  REAL, PRIVATE, SAVE :: &
155  & petap_nc,pjbbs_nc,pjbex_nc,pjboh_nc
156 
157  !Dimensions
158  INTEGER, PRIVATE, SAVE :: &
159  & mf_nc,mi_nc,ms_nc,mz_nc
160 
161  !Electron isotope identification
162  INTEGER, PRIVATE, SAVE :: &
163  & imel_nc
164 
165  !Isotope arrays
166  REAL(KIND=rspec), PRIVATE, SAVE, ALLOCATABLE :: &
167  & vti_nc(:),amntii_nc(:,:), &
168  & calmi_nc(:,:,:),calnii_nc(:,:,:,:), &
169  & capmii_nc(:,:,:,:),capnii_nc(:,:,:,:)
170 
171  !Isotope and charge arrays
172  REAL(KIND=rspec), PRIVATE, SAVE, ALLOCATABLE :: &
173  & grppiz_nc(:,:)
174 
175  !Species arrays
176  INTEGER, PRIVATE, SAVE, ALLOCATABLE :: &
177  & jms_nc(:),jzs_nc(:)
178 
179  REAL(KIND=rspec), PRIVATE, SAVE, ALLOCATABLE :: &
180  & sqzs_nc(:),xis_nc(:),ymus_nc(:,:,:), &
181  & bsjbps_nc(:),bsjbts_nc(:), &
182  & upars_nc(:,:,:),uthetas_nc(:,:,:), &
183  & gfls_nc(:,:),dpss_nc(:,:),dtss_nc(:,:), &
184  & qfls_nc(:,:),chipss_nc(:,:),chitss_nc(:,:), &
185  & tauss_nc(:,:)
186  !
187  !Physical constants, mathematical constants, conversion factors
188  !
189  REAL(KIND=rspec), PARAMETER :: &
190  z_coulomb=1.6022e-19_rspec, & !Coulomb charge [coul]
191  z_epsilon0=8.8542e-12_rspec, & !permittivity of free space [F/m]
192  z_j7kv=1.6022e-16_rspec, & !energy conversion factor [J/keV]
193  z_pi=3.141592654_rspec, & !pi [-]
194  z_protonmass=1.6726e-27_rspec, & !proton mass [kg]
195  z_electronmass=9.1095e-31_rspec, & !electron mass [kg]
196  z_mu0 = 1.2566e-06_rspec !Permeability of free space : (Henry/meter)
197 
198  REAL(KIND=rspec), PRIVATE, PARAMETER :: &
199  one=1.0_rspec, & !REAL 1
200  zero=0.0_rspec !REAL 0
201 
202 CONTAINS
203 
204 SUBROUTINE nclass(m_i , m_z ,p_b2 , p_bm2 , p_eb, &
205  & p_fhat , p_fm ,p_ft , p_grbm2 , p_grphi, &
206  & p_gr2phi, p_ngrth ,amu_i , grt_i ,temp_i, &
207  & den_iz , fex_iz ,grp_iz ,ipr , iflag, &
208  & l_banana,l_pfirsch,l_potato,k_order, &
209  & c_den,c_potb,c_potl, &
210  & p_etap,p_jbbs,p_jbex,p_jboh, &
211  & m_s,jm_s,jz_s, &
212  & bsjbp_s,bsjbt_s, &
213  & gfl_s,dn_s,vnnt_s,vneb_s,vnex_s,dp_ss,dt_ss, &
214  & upar_s,utheta_s, &
215  & qfl_s,chi_s,vqnt_s,vqeb_s,vqex_s, &
216  & chip_ss,chit_ss, &
217  & calm_i,caln_ii,capm_ii,capn_ii,ymu_s, &
218  & sqz_s,xi_s,tau_ss)
219  !-------------------------------------------------------------------------------
220  !NCLASS calculates the neoclassical transport properties of a multiple
221  ! species axisymmetric plasma using k_order parallel and radial force
222  ! balance equations for each species
223  !References:
224  ! S.P.Hirshman, D.J.Sigmar, Nucl Fusion 21 (1981) 1079
225  ! W.A.Houlberg, K.C.Shaing, S.P.Hirshman, M.C.Zarnstorff, Phys Plasmas 4 (1997)
226  ! 3230
227  ! W.A.Houlberg 1/2002
228  !Input:
229  ! m_i -number of isotopes (> 1) [-]
230  ! m_z -highest charge state [-]
231  ! p_b2 -<B**2> [T**2]
232  ! p_bm2 -<1/B**2> [/T**2]
233  ! p_eb -<E.B> [V*T/m]
234  ! p_fhat -mu_0*F/(dPsi/dr) [rho/m]
235  ! p_fm(m) -poloidal moments of drift factor for PS [/m**2]
236  ! p_ft -trapped fraction [-]
237  ! p_grbm2 -<grad(rho)**2/B**2> [rho**2/m**2/T**2]
238  ! p_grphi -potential gradient Phi' [V/rho]
239  ! p_gr2phi -second potential gradient Psi'(Phi'/Psi')' [V/rho**2]
240  ! p_ngrth -<n.grad(Theta)> [/m]
241  ! amu_i(i) -atomic mass number [-]
242  ! grt_i(i) -temperature gradient [keV/rho]
243  ! temp_i(i) -temperature [keV]
244  ! den_iz(i,z) -density [/m**3]
245  ! fex_iz(3,i,z) -moments of external parallel force [T*n/m**3]
246  ! grp_iz(i,z) -pressure gradient [keV/m**3/rho]
247  !Output:
248  ! iflag -error and warning flag [-]
249  ! =-1 warning
250  ! =0 no warnings or errors
251  ! =1 error
252  !Optional input:
253  ! L_BANANA -option to include banana viscosity [logical]
254  ! L_PFIRSCH -option to include Pfirsch-Schluter viscosity [logical]
255  ! L_POTATO -option to include potato orbits [logical]
256  ! K_ORDER -order of v moments to be solved [-]
257  ! =2 u and q (default)
258  ! =3 u, q, and u2
259  ! =else error
260  ! C_DEN-density cutoff below which species is ignored (default 1.e10) [/m**3]
261  ! C_POTB-kappa(0)*Bt(0)/[2*q(0)**2] [T]
262  ! C_POTL-q(0)*R(0) [m]
263  !Optional output:
264  !
265  ! Terms summed over species
266  !
267  ! P_ETAP-parallel electrical resistivity [Ohm*m]
268  ! P_JBBS-<J_bs.B> [A*T/m**2]
269  ! P_JBEX-<J_ex.B> current response to fex_iz [A*T/m**2]
270  ! P_JBOH-<J_OH.B> Ohmic current [A*T/m**2]
271  !
272  ! Species mapping
273  !
274  ! M_S -number of species [ms>1]
275  ! JM_S(s) -isotope number of s [-]
276  ! JZ_S(s) -charge state of s [-]
277  !
278  ! Bootstrap current and electrical resistivity
279  !
280  ! BSJBP_S(s) -<J_bs.B> driven by unit p'/p of s [A*T*rho/m**2]
281  ! BSJBT_S(s) -<J_bs.B> driven by unit T'/T of s [A*T*rho/m**2]
282  !
283  ! Continuity equation
284  !
285  ! GFL_S(m,s)-radial particle flux comps of s [rho/m**3/s]
286  ! m=1, banana-plateau, p' and T'
287  ! m=2, Pfirsch-Schluter
288  ! m=3, classical
289  ! m=4, banana-plateau, <E.B>
290  ! m=5, banana-plateau, external parallel force fex_iz
291  ! DN_S(s) -diffusion coefficients (diag comp) [rho**2/s]
292  ! VNNT_S(s) -convection velocity (off diag p',T' comps) [rho/s]
293  ! VNEB_S(s) -<E.B> particle convection velocity [rho/s]
294  ! VNEX_S(s) -external force particle convection velocity [rho/s]
295  ! DP_SS(s1,s2)-diffusion coefficient of s2 on p'/p of s1 [rho**2/s]
296  ! DT_SS(s1,s2)-diffusion coefficient of s2 on T'/T of s1 [rho**2/s]
297  !
298  !Momentum equation
299  !
300  ! UPAR_S(3,m,s)-parallel flow of s from force m [T*m/s]
301  ! m=1, p', T', Phi'
302  ! m=2, <E.B>
303  ! m=3, fex_iz
304  ! UTHETA_S(3,m,s)-poloidal flow of s from force m [m/s/T]
305  ! m=1, p', T'
306  ! m=2, <E.B>
307  ! m=3, fex_iz
308  !
309  !Energy equation
310  !
311  ! QFL_S(m,s)-radial heat conduction flux comps of s [W*rho/m**3]
312  ! m=1, banana-plateau, p' and T'
313  ! m=2, Pfirsch-Schluter
314  ! m=3, classical
315  ! m=4, banana-plateau, <E.B>
316  ! m=5, banana-plateau, external parallel force fex_iz
317  ! CHI_S(s) -conduction coefficients (diag comp) [rho**2/s]
318  ! VQNT_S(s) -conduction velocity (off diag p',T' comps) [rho/s]
319  ! VQEB_S(s) -<E.B> heat convection velocity [rho/s]
320  ! VQEX_S(s) -external force heat convection velocity [rho/s]
321  ! CHIP_SS(s1,s2)-heat cond coefficient of s2 on p'/p of s1 [rho**2/s]
322  ! CHIT_SS(s1,s2)-heat cond coefficient of s2 on T'/T of s1 [rho**2/s]
323  !
324  !Friction coefficients
325  !
326  ! CAPM_II(K_ORDER,K_ORDER,m_i,m_i)-test particle (tp) friction matrix [-]
327  ! CAPN_II(K_ORDER,K_ORDER,m_i,m_i)-field particle (fp) friction matrix [-]
328  ! CALM_I(K_ORDER,K_ORDER,m_i)-tp eff friction matrix [kg/m**3/s]
329  ! CALN_II(K_ORDER,K_ORDER,m_i,m_i)-fp eff friction matrix [kg/m**3/s]
330  !
331  !Viscosity coefficients
332  !
333  ! YMU_S(s)-normalized viscosity for s [kg/m**3/s]
334  !
335  !Miscellaneous
336  !
337  ! SQZ_S(s) -orbit squeezing factor for s [-]
338  ! XI_S(s) -charge weighted density factor of s [-]
339  ! TAU_SS(s1,s2) -90 degree scattering time [s]
340  !-------------------------------------------------------------------------------
341  !Declaration of input variables
342 implicit none
343 INTEGER, INTENT(IN) :: &
344  & m_i,m_z,ipr
345 
346 REAL(KIND=rspec), INTENT(IN) :: &
347  & p_b2,p_bm2,p_eb,p_fhat,p_fm(:),p_ft,p_grbm2,p_grphi,p_gr2phi, &
348  & p_ngrth, &
349  & amu_i(:),grt_i(:),temp_i(:), &
350  & den_iz(:,:),grp_iz(:,:),fex_iz(:,:,:)
351 
352 !Declaration of output variables
353 INTEGER, INTENT(OUT) :: &
354  & iflag
355 
356 !Declaration of optional input variables
357 LOGICAL, INTENT(IN), OPTIONAL :: &
358  & L_BANANA,L_PFIRSCH,L_POTATO
359 
360 INTEGER, INTENT(IN), OPTIONAL :: &
361  & K_ORDER
362 
363 REAL(KIND=rspec), INTENT(IN), OPTIONAL :: &
364  & C_DEN,C_POTB,C_POTL
365 
366 !Declaration of optional output variables
367 REAL(KIND=rspec), INTENT(OUT), OPTIONAL :: &
368  & P_ETAP,P_JBBS,P_JBEX,P_JBOH
369 
370 INTEGER, INTENT(OUT), OPTIONAL :: &
371  & M_S,JM_S(:),JZ_S(:)
372 
373 REAL(KIND=rspec), INTENT(OUT), OPTIONAL :: &
374  & BSJBP_S(:),BSJBT_S(:)
375 
376 REAL(KIND=rspec), INTENT(OUT), OPTIONAL :: &
377  & DP_SS(:,:),DT_SS(:,:), &
378  & GFL_S(:,:),DN_S(:),VNEB_S(:),VNEX_S(:),VNNT_S(:)
379 
380 REAL(KIND=rspec), INTENT(OUT), OPTIONAL :: &
381  & UPAR_S(:,:,:),UTHETA_S(:,:,:)
382 
383 REAL(KIND=rspec), INTENT(OUT), OPTIONAL :: &
384  & CHIP_SS(:,:),CHIT_SS(:,:), &
385  & QFL_S(:,:),CHI_S(:),VQEB_S(:),VQEX_S(:),VQNT_S(:)
386 
387 REAL(KIND=rspec), INTENT(OUT), OPTIONAL :: &
388  & CAPM_II(:,:,:,:),CAPN_II(:,:,:,:),CALM_I(:,:,:),CALN_II(:,:,:,:)
389 
390 REAL(KIND=rspec), INTENT(OUT), OPTIONAL :: &
391  & YMU_S(:,:,:)
392 
393 REAL(KIND=rspec), INTENT(OUT), OPTIONAL :: &
394  & SQZ_S(:),XI_S(:),TAU_SS(:,:)
395 
396 !Declaration of local variables
397 INTEGER :: &
398  & i,iza,k,l
399 
400 REAL(KIND=rspec) :: &
401  & dent, rval
402 
403 REAL(KIND=rspec), ALLOCATABLE :: &
404  & denz2(:)
405 
406 !-------------------------------------------------------------------------------
407 !Initialization
408 !-------------------------------------------------------------------------------
409 !Warning and error flag
410 iflag = 0
411 if (ipr .le. 10) then
412  ! !,*) ' -- dans fortran 90 --',m_i
413 endif
414 ! write(*,*) 'iflag=',iflag,' ipr=',ipr
415 ! write(*,*) 'M_S',M_S
416 ! write(*,*) 'JM_S',JM_S
417 ! CALL mexErrMsgTxt('stop')
418 
419 !Number of velocity moments
420 IF(present(k_order)) THEN
421 
422  IF(k_order < 2 .OR. k_order > 3) THEN
423 
424  iflag=1
425  write(*,*) ' K_ORDER =',k_order
426  write(*,*) 'NCLASS(1)/ERROR: K_ORDER must be 2 or 3'
427  goto 9999
428 
429  ELSE
430 
431  k_order_nc=k_order
432 
433  ENDIF
434 
435 ELSE
436 
437  k_order_nc=2
438 
439 ENDIF
440 
441 !Banana contribution to viscosity
442 !On by default, unless user turns it off or specifies p_ft=0
443 ! write(*,*) 'p_ft=',p_ft
444 IF(abs(p_ft) > 0.0_rspec) THEN
445 
446  l_banana_nc=.true.
447 
448 ELSE
449 
450  l_banana_nc=.false.
451 
452 ENDIF
453 
454 IF(present(l_banana) .AND. (.NOT. l_banana)) l_banana_nc=.false.
455 
456 !Pfirsch-Schluter contribution to viscosity
457 !On by default, unless user turns it off or specifies p_fm=0
458 
459 IF(abs(sum(p_fm(:))) > 0.0_rspec) THEN
460 
461  l_pfirsch_nc=.true.
462  mf_nc=SIZE(p_fm)
463 
464 ELSE
465 
466  l_pfirsch_nc=.false.
467  mf_nc=0
468 
469 ENDIF
470 
471 IF(present(l_pfirsch) .AND. (.NOT. l_pfirsch)) &
472  & l_pfirsch_nc=.false.
473 
474 !Potato orbit contribution to viscosity
475 !Off by default, unless user turns it on and provides constants
476 IF(present(l_potato) .AND. &
477  & present(c_potb) .AND. &
478  & present(c_potl)) THEN
479 
480  IF(l_potato) THEN
481 
482  l_potato_nc=.true.
483  cpotb_nc=c_potb
484  cpotl_nc=c_potl
485 
486  ELSE
487 
488  l_potato_nc=.false.
489  cpotb_nc=0.0_rspec
490  cpotl_nc=0.0_rspec
491 
492  ENDIF
493 
494 ELSE
495 
496  l_potato_nc=.false.
497  cpotb_nc=0.0_rspec
498  cpotl_nc=0.0_rspec
499 
500 ENDIF
501 
502 !No potato viscosity if there are no banana or Pfirsch-Schluter terms
503 IF( (.NOT. l_banana_nc) .AND. (.NOT. l_pfirsch_nc)) THEN
504 
505  l_potato_nc=.false.
506  cpotb_nc=0.0_rspec
507  cpotl_nc=0.0_rspec
508 
509 ENDIF
510 
511 !Cutoff density
512 IF(present(c_den)) THEN
513 
514  !Use input value
515  cden_nc=c_den
516 
517 ELSE
518 
519  !Use default value
520  cden_nc=1.0e10_rspec
521 
522 ENDIF
523 
524 !Trapped fraction between 0 and 1 inclusive
525 IF((p_ft < 0.0_rspec) .OR. (p_ft > 1.0_rspec)) THEN
526 
527  iflag=1
528  write(*,*) 'NCLASS(2)/ERROR: must have 0 <= p_ft <= 1'
529  goto 9999
530 
531 ENDIF
532 
533 !Allocate private data
534 ! write(*,*) 'lancement de nclass_init',cden_nc
535 CALL nclass_init(m_i,m_z,amu_i,den_iz,iflag)
536 
537 !Check messages
538 IF(iflag /= 0) THEN
539 
540  ! write(*,*) 'NCLASS(3)/'
541  IF(iflag > 0) goto 9999
542 
543 ENDIF
544 
545 !-------------------------------------------------------------------------------
546 !Set various plasma properties
547 !-------------------------------------------------------------------------------
548 !Friction coefficients
549 ! write(*,*) 'lancement de mn',temp_i
550 CALL nclass_mn(amu_i,temp_i)
551 
552 !Thermal velocities
553 vti_nc(1:mi_nc)=sqrt(2.0_rspec*z_j7kv*temp_i(1:mi_nc) &
554  & /amu_i(1:mi_nc)/z_protonmass)
555 ! write(*,*) 'mi_nc=',mi_nc
556 ! write(*,*) 'Te=',temp_i(1:mi_nc)
557 !Collision times
558 ! write(*,*) 'lancement de collision'
559 CALL nclass_tau(amu_i,temp_i,den_iz)
560 
561 !Reduced friction coefficients
562 calmi_nc(:,:,:)=0.0_rspec
563 
564 DO i=1,mi_nc !Over isotopes
565 
566  DO k=1,k_order_nc !Over first moment
567 
568  DO l=1,k_order_nc !Over second moment
569 
570  !Test particle component
571  calmi_nc(k,l,i)= &
572  & sum(amntii_nc(i,1:mi_nc)*capmii_nc(k,l,i,1:mi_nc))
573  ! write(*,*) 'i=',i,' k=',k,' l=',l
574  ! write(*,*) 'amntii_nc(i,1:mi_nc)=',amntii_nc(i,1)
575  !Field particle component
576  calnii_nc(k,l,i,1:mi_nc)= &
577  & amntii_nc(i,1:mi_nc)*capnii_nc(k,l,i,1:mi_nc)
578  if (k .eq. 1 .and. l .lt. 3 .and. i .eq. 1) then
579  ! write(*,*) 'calmi_nc',calmi_nc(k,l,i),l
580  endif
581 
582  END DO !Over second moment
583 
584  END DO !Over first moment
585 
586 END DO !Over isotopes
587 
588 !Species charge state density factor, total pressure, squeezing
589 dent=0.0_rspec
590 ALLOCATE(denz2(1:(mi_nc+ms_nc)))
591 ! write(*,*) 'ms_nc=',ms_nc,' mi_nc=',mi_nc
592 DO i=1,mi_nc !Over isotopes
593 
594  denz2(i)=0.0_rspec
595 
596  DO iza=1,mz_nc !Over charge states
597 
598  IF(den_iz(i,iza) > cden_nc) THEN
599 
600  denz2(i)=denz2(i)+den_iz(i,iza)*iza**2
601  dent=dent+den_iz(i,iza)*temp_i(i)
602 
603  ENDIF
604 
605  ENDDO !Over charge states
606 
607 ENDDO !Over isotopes
608 
609 DO i=1,ms_nc !Over species
610 
611  iza=abs(jzs_nc(i))
612  xis_nc(i)=den_iz(jms_nc(i),iza)*REAL(jzs_nc(i),rspec)**2 &
613  & /denz2(jms_nc(i))
614  sqzs_nc(i)=1.0_rspec+p_fhat**2/p_b2*amu_i(jms_nc(i)) &
615  & *z_protonmass*p_gr2phi &
616  & /(z_coulomb*REAL(jzs_nc(i),rspec))
617  if(sqzs_nc(i) > 10.0_rspec) sqzs_nc(i)=10.0_rspec
618  IF(sqzs_nc(i) < 0.5_rspec) sqzs_nc(i)=0.5_rspec
619 
620 ENDDO !Over species
621 
622 DEALLOCATE(denz2)
623 
624 !Normalized viscosities
625 CALL nclass_mu(p_fm,p_ft,p_ngrth,amu_i,temp_i,den_iz,ipr)
626 !Add potential gradient to pressure gradient
627 DO i=1,ms_nc !Over species
628  iza=abs(jzs_nc(i))
629  grppiz_nc(jms_nc(i),iza)=grp_iz(jms_nc(i),iza) &
630  & +p_grphi*den_iz(jms_nc(i),iza) &
631  & *REAL(jzs_nc(i),rspec)*z_coulomb/z_j7kv
632  rval = p_grphi*den_iz(jms_nc(i),iza) &
633  & *REAL(jzs_nc(i),rspec)*z_coulomb/z_j7kv
634 
635 END DO !Over species
636 !-------------------------------------------------------------------------------
637 !Solve for parallel flows within a surface
638 !-------------------------------------------------------------------------------
639 iflag=0
640 ! write(*,*) 'lancement de nclass_flow ipr=',ipr
641 ! write(*,*) 'Te=',temp_i
642 ! write(*,*) 'grt=',grt_i
643 ! write(*,*) 'ne(1,1)=',den_iz(1,1)
644 ! write(*,*) 'pfhat=',p_fhat,'b2=',p_b2
645 ! write(*,*) 'flow, gradient Te',grt_i
646 ! write(*,*) 'flow, fex',fex_iz(1,1,1)
647 ! write(*,*) 'flow, den',den_iz(1,1)
648 
649 CALL nclass_flow(p_b2,p_bm2,p_eb,p_fhat,p_grbm2,grt_i,temp_i, &
650  & den_iz,fex_iz,iflag,ipr)
651 ! CALL mexErrMsgTxt('stop')
652 
653 
654 !Check messages
655 IF(iflag /= 0) THEN
656 
657  write(*,*) 'NCLASS(4)/'
658  IF(iflag > 0) goto 9999
659 
660 ENDIF
661 
662 !-------------------------------------------------------------------------------
663 !Calculate poloidal velocity from parallel velocity
664 !-------------------------------------------------------------------------------
665 DO i=1,ms_nc !Over species
666 
667  iza=iabs(jzs_nc(i))
668 
669  DO k=1,k_order_nc !Over moments
670 
671  DO l=1,3 !Over forces
672 
673  uthetas_nc(k,l,i)=upars_nc(k,l,i)/p_b2
674  ! if (l .eq. 2 .and. ipr .eq. 1 .and. k .eq. 1) then
675  ! write(*,*) 'uthetas_nc(1,2,i)=',uthetas_nc(1,2,i)
676  ! write(*,*) 'upars_nc(1,2,i)=',upars_nc(1,2,i)
677  ! write(*,*) 'i=',i,' p_b2=',p_b2
678  ! endif
679 
680  IF(l == 1) THEN
681 
682  IF(k == 1) THEN
683 
684  !Add p' and Phi'
685  uthetas_nc(k,l,i)=uthetas_nc(k,l,i) &
686  & +p_fhat*grppiz_nc(jms_nc(i),iza)*z_j7kv &
687  & /(z_coulomb*REAL(jzs_nc(i),rspec) &
688  & *den_iz(jms_nc(i),iza))/p_b2
689 
690  ELSEIF(k == 2) THEN
691 
692  !Add T'
693  uthetas_nc(k,l,i)=uthetas_nc(k,l,i) &
694  & +p_fhat*z_j7kv*grt_i(jms_nc(i)) &
695  & /(REAL(jzs_nc(i),rspec)*z_coulomb*p_b2)
696 
697  ENDIF
698 
699  ENDIF
700  if (ipr .eq. 1 .and. l .eq. 2) then
701  ! write(*,*) 'i=',i
702  ! write(*,*) 'uthetas_nc(1,2,i)=',uthetas_nc(1,2,i)
703  ! if (i .eq. ms_nc) stop
704  endif
705 
706  ENDDO !Over forces
707 
708  ENDDO !Over moments
709 
710 ENDDO !Over species
711 
712 !-------------------------------------------------------------------------------
713 !Output
714 !-------------------------------------------------------------------------------
715 !Species mapping
716 IF(present(m_s)) m_s=ms_nc
717 IF(present(jm_s)) jm_s(1:ms_nc)=jms_nc(1:ms_nc)
718 IF(present(jz_s)) jz_s(1:ms_nc)=jzs_nc(1:ms_nc)
719 
720 !Electrical resistivity and currents
721 IF(present(p_etap)) p_etap=petap_nc
722 ! write(*,*) 'transfert boot :',pjbbs_nc
723 IF(present(p_jbbs)) p_jbbs=pjbbs_nc
724 IF(present(p_jbex)) p_jbex=pjbex_nc
725 IF(present(p_jboh)) p_jboh=pjboh_nc
726 IF(present(bsjbp_s)) bsjbp_s(1:ms_nc)=bsjbps_nc(1:ms_nc)
727 IF(present(bsjbt_s)) bsjbt_s(1:ms_nc)=bsjbts_nc(1:ms_nc)
728 
729 !Continuity equation
730 IF(present(dp_ss)) dp_ss(1:ms_nc,1:ms_nc)=dpss_nc(1:ms_nc,1:ms_nc)
731 IF(present(dt_ss)) dt_ss(1:ms_nc,1:ms_nc)=dtss_nc(1:ms_nc,1:ms_nc)
732 IF(present(gfl_s)) gfl_s(1:5,1:ms_nc)=gfls_nc(1:5,1:ms_nc)
733 
734 IF(present(dn_s)) THEN
735 
736  DO i=1,ms_nc !Over species
737 
738  dn_s(i)=dpss_nc(i,i)
739 
740  ENDDO !Over species
741 
742 ENDIF
743 
744 IF(present(vnnt_s)) THEN
745 
746  DO i=1,ms_nc !Over species
747 
748  iza=abs(jzs_nc(i))
749  vnnt_s(i)=(sum(gfls_nc(1:3,i)) &
750  & +dn_s(i)*(grp_iz(jms_nc(i),iza) &
751  & -den_iz(jms_nc(i),iza)*grt_i(jms_nc(i))) &
752  & /temp_i(jms_nc(i)))/den_iz(jms_nc(i),iza)
753 
754  ENDDO !Over species
755 
756 ENDIF
757 
758 IF(present(vneb_s)) THEN
759 
760  DO i=1,ms_nc !Over species
761 
762  iza=abs(jzs_nc(i))
763  vneb_s(i)=gfls_nc(4,i)/den_iz(jms_nc(i),iza)
764 
765  ENDDO !Over species
766 
767 ENDIF
768 
769 IF(present(vnex_s)) THEN
770 
771  DO i=1,ms_nc !Over species
772 
773  iza=abs(jzs_nc(i))
774  vnex_s(i)=gfls_nc(5,i)/den_iz(jms_nc(i),iza)
775 
776  ENDDO !Over species
777 
778 ENDIF
779 
780 !Momentum equation
781 
782 IF(present(upar_s)) upar_s(1:k_order_nc,1:3,1:ms_nc) &
783  & =upars_nc(1:k_order_nc,1:3,1:ms_nc)
784 IF(present(utheta_s)) utheta_s(1:k_order_nc,1:3,1:ms_nc) &
785  & =uthetas_nc(1:k_order_nc,1:3,1:ms_nc)
786 
787 !Energy equation
788 IF(present(qfl_s)) qfl_s(1:5,1:ms_nc)=qfls_nc(1:5,1:ms_nc)
789 ! write(6,*) 'qfls_nc(1,1)=',qfls_nc(1,1)
790 IF(present(chi_s)) THEN
791 
792  DO i=1,ms_nc !Over species
793 
794  chi_s(i)=chipss_nc(i,i)+chitss_nc(i,i)
795 
796  ENDDO !Over species
797 
798 ENDIF
799 
800 IF(present(vqnt_s)) THEN
801 
802  DO i=1,ms_nc !Over species
803 
804  iza=abs(jzs_nc(i))
805  vqnt_s(i)=(sum(qfls_nc(1:3,i))/z_j7kv/den_iz(jms_nc(i),iza) &
806  & +chi_s(i)*grt_i(jms_nc(i)))/temp_i(jms_nc(i))
807 
808  ENDDO !Over species
809 
810 ENDIF
811 
812 IF(present(vqeb_s)) THEN
813 
814  DO i=1,ms_nc !Over species
815 
816  iza=abs(jzs_nc(i))
817  vqeb_s(i)=qfls_nc(4,i)/den_iz(jms_nc(i),iza)/temp_i(jms_nc(i)) &
818  & /z_j7kv
819 
820  ENDDO !Over species
821 
822 ENDIF
823 
824 IF(present(vqex_s)) THEN
825 
826  DO i=1,ms_nc !Over species
827 
828  iza=abs(jzs_nc(i))
829  vqex_s(i)=qfls_nc(5,i)/den_iz(jms_nc(i),iza)/temp_i(jms_nc(i)) &
830  & /z_j7kv
831 
832  ENDDO !Over species
833 
834 ENDIF
835 
836 IF(present(chip_ss)) chip_ss(1:ms_nc,1:ms_nc) &
837  & =chipss_nc(1:ms_nc,1:ms_nc)
838 IF(present(chit_ss)) chit_ss(1:ms_nc,1:ms_nc) &
839  & =chitss_nc(1:ms_nc,1:ms_nc)
840 
841 !Friction coefficients
842 IF(present(calm_i)) &
843  & calm_i(1:k_order_nc,1:k_order_nc,1:mi_nc) &
844  & =calmi_nc(1:k_order_nc,1:k_order_nc,1:mi_nc)
845 IF(present(caln_ii)) &
846  & caln_ii(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc) &
847  & =calnii_nc(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc)
848 IF(present(capm_ii)) &
849  & capm_ii(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc) &
850  & =capmii_nc(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc)
851 IF(present(capn_ii)) &
852  & capn_ii(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc) &
853  & =capnii_nc(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc)
854 
855 !Viscosity coefficients
856 IF(present(ymu_s)) ymu_s(1:k_order_nc,1:k_order_nc,1:ms_nc) &
857  & =ymus_nc(1:k_order_nc,1:k_order_nc,1:ms_nc)
858 
859 !Miscellaneous
860 IF(present(sqz_s)) sqz_s(1:ms_nc)=sqzs_nc(1:ms_nc)
861 IF(present(xi_s)) xi_s(1:ms_nc)=xis_nc(1:ms_nc)
862 IF(present(tau_ss)) tau_ss(1:ms_nc,1:ms_nc) &
863  & =tauss_nc(1:ms_nc,1:ms_nc)
864 
865 9999 CONTINUE
866 
867 END SUBROUTINE nclass
868 
869 SUBROUTINE nclass_flow(p_b2,p_bm2,p_eb,p_fhat,p_grbm2,grt_i, &
870  & temp_i,den_iz,fex_iz,iflag,ipr)
871  !-------------------------------------------------------------------------------
872  !NCLASS_FLOW calculates the neoclassical flows u0, u1=q/p (and and optionally
873  ! u2) plus some of the other transport properties
874  !References:
875  ! S.P.Hirshman, D.J.Sigmar, Nucl Fusion 21 (1981) 1079
876  ! W.A.Houlberg, K.C.Shaing, S.P.Hirshman, M.C.Zarnstorff, Phys Plasmas 4 (1997)
877  ! 3230
878  ! W.A.Houlberg 1/2002
879  !Input:
880  ! p_b2 -<B**2> [T**2]
881  ! p_bm2 -<1/B**2> [/T**2]
882  ! p_eb -<E.B> [V*T/m]
883  ! p_fhat -mu_0*F/(dPsi/dr) [rho/m]
884  ! p_grbm2 -<grad(rho)**2/B**2> [rho**2/m**2/T**2]
885  ! grt_i(i) -temperature gradient [keV/rho]
886  ! temp_i(i) -temperature [keV]
887  ! den_iz(i,z) -density [/m**3]
888  ! fex_iz(3,i,z) -moments of external parallel force [T*n/m**3]
889  !Output:
890  ! iflag -error and warning flag [-]
891  ! =-1 warning
892  ! =0 no warnings or errors
893  ! =1 error
894  !-------------------------------------------------------------------------------
895 implicit none
896  !Declaration of input variables
897 REAL(KIND=rspec), INTENT(IN) :: &
898  & p_b2,p_bm2,p_eb,p_fhat,p_grbm2, &
899  & grt_i(:),temp_i(:),den_iz(:,:),fex_iz(:,:,:)
900 
901 INTEGER, INTENT(in) :: &
902  & ipr
903 
904 !Declaration of output variables
905 INTEGER, INTENT(OUT) :: &
906  & iflag
907 
908 
909 !Declaration of local variables
910 INTEGER :: &
911  & i,im,iz,iza,j,jm,k,l,l1,m,m1
912 
913 INTEGER :: &
914  & indxk(1:k_order_nc),indxki(1:k_order_nc*mi_nc)
915 
916 REAL(KIND=rspec) :: &
917  & cbp,cbpa,cbpaq,cc,ccl,ccla,cclaq, cclb,cclbq,cps,cpsa,cpsaq, &
918  & cpsb,cpsbq,d,denzc
919 
920 REAL(KIND=rspec) :: &
921  & akk(1:k_order_nc,1:k_order_nc), &
922  & xk(1:k_order_nc), &
923  & akiki(1:k_order_nc*mi_nc,1:k_order_nc*mi_nc), &
924  & akikis(1:k_order_nc*mi_nc,1:k_order_nc*mi_nc), &
925  & xki3(1:k_order_nc*mi_nc,3), &
926  & crk6i(1:k_order_nc,6,mi_nc), &
927  & rk6s(1:k_order_nc,6,ms_nc), &
928  & srcth(1:k_order_nc,ms_nc),srcthp(ms_nc),srctht(ms_nc), &
929  & crhatp(1:k_order_nc,ms_nc,mi_nc), &
930  & crhatt(1:k_order_nc,ms_nc,mi_nc), &
931  & rhatp(1:k_order_nc,ms_nc,ms_nc), &
932  & rhatt(1:k_order_nc,ms_nc,ms_nc), &
933  & uaip(1:k_order_nc,ms_nc,ms_nc), &
934  & uait(1:k_order_nc,ms_nc,ms_nc), &
935  & xabp(k_order_nc*mi_nc,ms_nc), &
936  & xabt(k_order_nc*mi_nc,ms_nc)
937 !
938 ! LAPACK decomposition
939 !
940 REAL(KIND=rspec) :: &
941  & alp(1:k_order_nc,1:k_order_nc), &
942  & alps(1:k_order_nc,1:k_order_nc), &
943  & vecm(1:k_order_nc), &
944  & alp2(1:k_order_nc*mi_nc,1:k_order_nc*mi_nc), &
945  & alp2s(1:k_order_nc*mi_nc,1:k_order_nc*mi_nc), &
946  & vecm2(1:k_order_nc*mi_nc), &
947  & ama(1:k_order_nc), &
948  & ama2(1:k_order_nc*mi_nc), &
949  & sumdgesv(1:k_order_nc*mi_nc)
950 REAL(KIND=rspec) :: &
951  & ipvt(1:k_order_nc), &
952  & rcond,z(1:k_order_nc)
953 
954 INTEGER :: &
955  & indl2(1:k_order_nc*mi_nc), &
956  & indl(1:k_order_nc),n,job,lda,ichoix
957 
958 !-------------------------------------------------------------------------------
959 !Initialization
960 !-------------------------------------------------------------------------------
961 !Local arrays
962 ! CALL mexErrMsgTxt('stop')
963 
964 akk(:,:)=0.0_rspec
965 xk(:)=0.0_rspec
966 akiki(:,:)=0.0_rspec
967 xki3(:,:)=0.0_rspec
968 crk6i(:,:,:)=0.0_rspec
969 rk6s(:,:,:)=0.0_rspec
970 srcth(:,:)=0.0_rspec
971 srcthp(:)=0.0_rspec
972 srctht(:)=0.0_rspec
973 crhatp(:,:,:)=0.0_rspec
974 crhatt(:,:,:)=0.0_rspec
975 rhatp(:,:,:)=0.0_rspec
976 rhatt(:,:,:)=0.0_rspec
977 uaip(:,:,:)=0.0_rspec
978 uait(:,:,:)=0.0_rspec
979 xabp(:,:)=0.0_rspec
980 xabt(:,:)=0.0_rspec
981 
982 !Private data
983 petap_nc=0.0_rspec
984 pjbbs_nc=0.0_rspec
985 pjbex_nc=0.0_rspec
986 pjboh_nc=0.0_rspec
987 gfls_nc(:,:)=0.0_rspec
988 qfls_nc(:,:)=0.0_rspec
989 
990 !Local constants
991 cc=(p_fhat/z_coulomb)*z_j7kv
992 cbp=p_fhat/p_b2/z_coulomb
993 cps=(p_fhat/z_coulomb)*(1.0/p_b2-p_bm2)
994 ccl=(p_grbm2/z_coulomb)/p_fhat
995 ! write(*,*) 'cc=',cc,p_fhat,z_coulomb,z_j7kv
996 
997 !-------------------------------------------------------------------------------
998 !Calculate species information for reduced charge state formalism
999 !-------------------------------------------------------------------------------
1000 DO i=1,ms_nc !Over species
1001 
1002  !Isotopic and charge state indices
1003  im=jms_nc(i)
1004  iz=jzs_nc(i)
1005  iza=iabs(iz)
1006 
1007  !Response matrix
1008  akk(:,:)=xis_nc(i)*calmi_nc(:,:,im)-ymus_nc(:,:,i)
1009 
1010  !Get lu decomposition of response matrix
1011  iflag=0
1012 
1013  ! CALL NCLASS_DECOMP(akk,k_order_nc,indxk,d,iflag,ipr)
1014  ! write(*,*) 'fin nclass_decomp i=',i,im,iz,iza
1015 
1016  !Check messages
1017  ! IF(iflag == 1) THEN
1018 
1019  ! write(*,*) 'NCLASS_FLOW(1)/'
1020  ! GOTO 9999
1021 
1022  ! ENDIF
1023 
1024  !Source and response from lambda terms involving isotopic flows
1025  ichoix = 0
1026  if (ichoix .eq. 0) then
1027  DO k=1,k_order_nc !Over velocity moments
1028  ! write(*,*) '---------------------',k,'----------------'
1029  rk6s(k,k,i)= xis_nc(i)
1030  alp = xis_nc(i)*calmi_nc(:,:,im)-ymus_nc(:,:,i)
1031  alps = alp
1032  ama = rk6s(:,k,i)
1033  ! write(*,*) 'B=',ama(1),ama(2),ama(3)
1034  call normalisation(alp,ama,k_order_nc,vecm)
1035  ! write(*,*) 'BN=',ama(1),ama(2),ama(3)
1036  ! if (vecm(1) .le. 1e-6) then
1037  ! write(*,*) 'vecm=',vecm(1),vecm(2),vecm(3)
1038  ! endif
1039  call dgesv(k_order_nc,1,alp,k_order_nc, &
1040  & indl,ama,k_order_nc,iflag)
1041 
1042  ! write(*,*) 'BS=',ama(1),ama(2),ama(3)
1043  ! write(*,*) 'A1=',alps(1,1)*ama(1)+alps(1,2)*ama(2) &
1044  ! & +alps(1,3)*ama(3)
1045  ! write(*,*) 'A2=',alps(2,1)*ama(1)+alps(2,2)*ama(2) &
1046  ! & +alps(2,3)*ama(3)
1047  ! write(*,*) 'A3=',alps(3,1)*ama(1)+alps(3,2)*ama(2) &
1048  ! & +alps(3,3)*ama(3)
1049  ! CALL NCLASS_BACKSUB(akk,k_order_nc,indxk,rk6s(:,k,i),ipr)
1050  crk6i(:,k,im)=crk6i(:,k,im)+xis_nc(i)*ama
1051  rk6s(:,k,i) = ama
1052  ! write(*,*) 'rk6s(:,k,i)=',rk6s(:,k,i)
1053  END DO !Over velocity moments
1054  endif
1055  ! stop
1056  !Source and response from poloidal source (p' and T') terms
1057  srcth(1,i)=(cc/iz)*grppiz_nc(im,iza)/den_iz(im,iza)
1058  srcth(2,i)=(cc/iz)*grt_i(im)
1059 
1060  DO k=1,k_order_nc !Over velocity moments
1061 
1062  rk6s(k,4,i)=sum(srcth(:,i)*ymus_nc(k,:,i))
1063 
1064  END DO !Over velocity moments
1065 
1066  alp = xis_nc(i)*calmi_nc(:,:,im)-ymus_nc(:,:,i)
1067  ama = rk6s(:,4,i)
1068  ! write(*,*) 'ama=',ama
1069  call normalisation(alp,ama,k_order_nc,vecm)
1070  call dgesv(k_order_nc,1,alp,k_order_nc, &
1071  & indl,ama,k_order_nc,iflag)
1072  ! CALL NCLASS_BACKSUB(akk,k_order_nc,indxk,rk6s(:,4,i),ipr)
1073 
1074  crk6i(:,4,im) = crk6i(:,4,im)+xis_nc(i)*ama
1075  rk6s(:,4,i) = ama
1076  !Source and response from unit p'/p and T'/T terms for decomposition of fluxes
1077  srcthp(i) =-(cc/iz)*temp_i(im)
1078  srctht(i) =-(cc/iz)*temp_i(im)
1079  rhatp(:,i,i) = srcthp(i)*ymus_nc(:,1,i)
1080  rhatt(:,i,i) = srctht(i)*ymus_nc(:,2,i)
1081 
1082  alp = xis_nc(i)*calmi_nc(:,:,im)-ymus_nc(:,:,i)
1083  ama = rhatp(1:k_order_nc,i,i)
1084  call normalisation(alp,ama,k_order_nc,vecm)
1085  call dgesv(k_order_nc,1,alp,k_order_nc, &
1086  & indl,ama,k_order_nc,iflag)
1087  rhatp(1:k_order_nc,i,i) = ama
1088  crhatp(:,i,im)=crhatp(:,i,im)+xis_nc(i)*ama
1089  ! CALL NCLASS_BACKSUB(akk,k_order_nc,indxk, &
1090  ! & rhatp(1:k_order_nc,i,i),ipr)
1091  alp = xis_nc(i)*calmi_nc(:,:,im)-ymus_nc(:,:,i)
1092  ama = rhatt(1:k_order_nc,i,i)
1093  call normalisation(alp,ama,k_order_nc,vecm)
1094  call dgesv(k_order_nc,1,alp,k_order_nc, &
1095  & indl,ama,k_order_nc,iflag)
1096  ! CALL NCLASS_BACKSUB(akk,k_order_nc,indxk, &
1097  ! & rhatt(1:k_order_nc,i,i),ipr)
1098  rhatt(1:k_order_nc,i,i) = ama
1099  crhatt(:,i,im)=crhatt(:,i,im)+xis_nc(i)*ama
1100 
1101  !Source and response from parallel electric field terms for resistivity
1102  rk6s(1,5,i) =-iz*z_coulomb*den_iz(im,iza)
1103  rk6s(2,5,i) = 0.0
1104  alp = xis_nc(i)*calmi_nc(:,:,im)-ymus_nc(:,:,i)
1105  ama = rk6s(1:k_order_nc,5,i)
1106  call normalisation(alp,ama,k_order_nc,vecm)
1107  call dgesv(k_order_nc,1,alp,k_order_nc, &
1108  & indl,ama,k_order_nc,iflag)
1109  ! CALL NCLASS_BACKSUB(akk,k_order_nc,indxk,rk6s(1:k_order_nc,5,i),ipr)
1110  crk6i(:,5,im) = crk6i(:,5,im)+xis_nc(i)*ama
1111  rk6s(1:k_order_nc,5,i) = ama
1112  !Source and response from external force terms
1113  rk6s(1:k_order_nc,6,i)=-fex_iz(1:k_order_nc,im,iza)
1114  alp = xis_nc(i)*calmi_nc(:,:,im)-ymus_nc(:,:,i)
1115  ama = rk6s(1:k_order_nc,6,i)
1116  call normalisation(alp,ama,k_order_nc,vecm)
1117  call dgesv(k_order_nc,1,alp,k_order_nc, &
1118  & indl,ama,k_order_nc,iflag)
1119  ! CALL NCLASS_BACKSUB(akk,k_order_nc,indxk,rk6s(1:k_order_nc,6,i),ipr)
1120  rk6s(1:k_order_nc,6,i) = ama
1121  crk6i(:,6,im)=crk6i(:,6,im)+xis_nc(i)*ama
1122 
1123 END DO !Over species
1124 
1125 !-------------------------------------------------------------------------------
1126 !Load coefficient matrix and source terms for isotopic flows
1127 !-------------------------------------------------------------------------------
1128 DO im=1,mi_nc !Over isotopes 1
1129 
1130  DO m=1,k_order_nc !Over velocity moments 1
1131 
1132  m1 = im+(m-1)*mi_nc
1133  !Diagonal coefficients
1134  akiki(m1,m1) = 1.0_rspec
1135  !p' and T' force terms
1136  xki3(m1,1) = crk6i(m,4,im)
1137  ! write(*,*) 'xki3(m1,1)=',xki3(m1,1),m1,m,im
1138  !Unit p'/p and T'/T
1139  xabp(m1,:) = crhatp(m,:,im)
1140  xabt(m1,:) = crhatt(m,:,im)
1141 
1142  xki3(m1,2) = crk6i(m,5,im)
1143  !External force
1144  xki3(m1,3) = crk6i(m,6,im)
1145  ! Field particle friction
1146  DO jm=1,mi_nc !Over isotopes 2
1147 
1148  DO l=1,k_order_nc !Over velocity moments 2
1149 
1150  l1=jm+(l-1)*mi_nc
1151  ! if (ipr .eq. 1 .and. l1 .eq. 2) then
1152 
1153  ! write(*,*) ' avant somme, akiki(m1,2)=',akiki(m1,l1)
1154 
1155  ! endif
1156  akiki(m1,l1)=akiki(m1,l1) &
1157  & +sum(calnii_nc(1:k_order_nc,l,im,jm)*crk6i(m,1:k_order_nc,im))
1158 
1159  ENDDO !Over velocity moments 2
1160 
1161  ENDDO !Over isotopes 2
1162 
1163  ENDDO !Over velocity moments 1
1164 
1165 ENDDO !Over isotopes 1
1166 
1167 !-------------------------------------------------------------------------------
1168 !Get lu decomposition of coefficient matrix
1169 !-------------------------------------------------------------------------------
1170 iflag=0
1171 akikis = akiki
1172 ! CALL NCLASS_DECOMP(akiki,k_order_nc*mi_nc,indxki,d,iflag,ipr)
1173 
1174 !Check messages
1175 ! IF(iflag == 1) THEN
1176 
1177 ! write(*,*) 'NCLASS_FLOW(2)/'
1178 ! GOTO 9999
1179 
1180 ! ENDIF
1181 
1182 !-------------------------------------------------------------------------------
1183 !Evaluate isotopic flows from back substitution for each source
1184 !-------------------------------------------------------------------------------
1185 !xki3(1,k) to xki3(mi_nc,k) are the isotopic velocities
1186 !xki3(mi_nc+1,k) to xki3(2*mi_nc,k) are the isotopic heat flows
1187 !xki3(2*mi_nc+1,k) to xki3(3*mi_nc,k) are the u2 flows
1188 !Evaluate species flows
1189 DO k=1,3 !Over forces
1190  alp2 = akikis
1191  ama2 = xki3(1:k_order_nc*mi_nc,k)
1192  call normalisation(alp2,ama2,k_order_nc*mi_nc,vecm2)
1193  ! if (k .eq. 2) then
1194  ! write(*,*) 'ama2=',ama2
1195  ! endif
1196  ! stop
1197  call dgesv(k_order_nc*mi_nc,1,alp2,k_order_nc*mi_nc, &
1198  & indl2,ama2,k_order_nc*mi_nc,iflag)
1199  errsum = 0.0
1200  do l=1,k_order_nc*mi_nc
1201  sumdgesv(l) = 0.0
1202  do j=1,k_order_nc*mi_nc
1203  sumdgesv(l) = sumdgesv(l)+akikis(l,j)*ama2(j)
1204  enddo
1205  ! write(*,*) 'A',l,'=',sumdgesv(l),'/ ',xki3(l,k)
1206  if (xki3(l,k) /= 0) then
1207  errsum = abs(sumdgesv(l)-xki3(l,k))/xki3(l,k) + errsum
1208  endif
1209  enddo
1210  ! write(*,*) '----'
1211  ! write(*,*) 'errsum=',errsum
1212  xki3(1:k_order_nc*mi_nc,k) = ama2
1213  ! write(*,*) 'ama2=',ama2
1214 
1215 ENDDO !Over forces
1216 ! stop
1217 
1218 DO i=1,ms_nc !Over species
1219  alp2 = akikis
1220  ama2 = xabp(1:k_order_nc*mi_nc,i)
1221  call normalisation(alp2,ama2,k_order_nc*mi_nc,vecm2)
1222  call dgesv(k_order_nc*mi_nc,1,alp2,k_order_nc*mi_nc, &
1223  & indl2,ama2,k_order_nc*mi_nc,iflag)
1224  xabp(1:k_order_nc*mi_nc,i) = ama2
1225  ! CALL NCLASS_BACKSUB(akiki,k_order_nc*mi_nc,indxki, &
1226  ! & xabp(1:k_order_nc*mi_nc,i),ipr)
1227  alp2 = akikis
1228  ama2 = xabt(1:k_order_nc*mi_nc,i)
1229  call normalisation(alp2,ama2,k_order_nc*mi_nc,vecm2)
1230  call dgesv(k_order_nc*mi_nc,1,alp2,k_order_nc*mi_nc, &
1231  & indl2,ama2,k_order_nc*mi_nc,iflag)
1232  xabt(1:k_order_nc*mi_nc,i) = ama2
1233  ! CALL NCLASS_BACKSUB(akiki,k_order_nc*mi_nc,indxki, &
1234  ! & xabt(1:k_order_nc*mi_nc,i),ipr)
1235 
1236 ENDDO !Over species
1237 ! stop
1238 DO i=1,ms_nc !Over species 1
1239 
1240  im=jms_nc(i)
1241 
1242  !Force contributions
1243  DO m=1,3 !Over forces
1244 
1245  DO k=1,k_order_nc !Over velocity moments
1246 
1247  IF(m == 2) THEN
1248 
1249  !Add normalization to <E.B>
1250  upars_nc(k,m,i)=rk6s(k,5,i)
1251 
1252  ELSE
1253 
1254  upars_nc(k,m,i)=rk6s(k,m+3,i)
1255 
1256  ENDIF
1257 
1258  ENDDO !Over velocity moments
1259 
1260  !Response contributions
1261  DO jm=1,mi_nc !Over isotopes
1262 
1263  xk(:)=0.0_rspec
1264 
1265  DO l=1,k_order_nc !Over velocity moments
1266 
1267  l1=jm+(l-1)*mi_nc
1268  xk(:)=xk(:)-calnii_nc(1:k_order_nc,l,im,jm)*xki3(l1,m)
1269 
1270  ENDDO !Over velocity moments
1271 
1272  DO l=1,k_order_nc !Over velocity moments
1273 
1274  upars_nc(:,m,i)=upars_nc(:,m,i)+xk(l)*rk6s(:,l,i)
1275 
1276  ENDDO !Over velocity moments
1277 
1278  ENDDO !Over isotopes
1279 
1280  ENDDO !Over forces
1281  ! stop
1282  !Unit p'/p and T'/T
1283  DO j=1,ms_nc !Over species 2
1284 
1285  uaip(:,j,i)=rhatp(:,j,i)
1286  uait(:,j,i)=rhatt(:,j,i)
1287 
1288  !Response contributions
1289  DO jm=1,mi_nc !Over isotopes
1290 
1291  xk(:)=0.0_rspec
1292 
1293  DO l=1,k_order_nc !Over velocity moments
1294 
1295  l1=jm+(l-1)*mi_nc
1296  xk(:)=xk(:)-calnii_nc(1:k_order_nc,l,im,jm)*xabp(l1,j)
1297 
1298  ENDDO !Over velocity moments
1299 
1300  DO l=1,k_order_nc !Over velocity moments
1301 
1302  uaip(:,j,i)=uaip(:,j,i)+xk(l)*rk6s(:,l,i)
1303 
1304  ENDDO !Over velocity moments
1305 
1306  xk(:)=0.0_rspec
1307 
1308  DO l=1,k_order_nc !Over velocity moments
1309  l1=jm+(l-1)*mi_nc
1310 
1311  xk(:)=xk(:)-calnii_nc(1:k_order_nc,l,im,jm)*xabt(l1,j)
1312 
1313  ENDDO !Over velocity moments
1314 
1315  DO l=1,k_order_nc !Over velocity moments
1316 
1317  uait(:,j,i)=uait(:,j,i)+xk(l)*rk6s(:,l,i)
1318 
1319  ENDDO !Over velocity moments
1320 
1321  ENDDO !Over isotopes
1322 
1323  ENDDO !Over species 2
1324 
1325 ENDDO !Over species 1
1326 
1327 !-------------------------------------------------------------------------------
1328 !Evaluate currents and fluxes from flows
1329 !-------------------------------------------------------------------------------
1330 bsjbps_nc(:) = 0.0_rspec
1331 bsjbts_nc(:) = 0.0_rspec
1332 dpss_nc(:,:) = 0.0_rspec
1333 dtss_nc(:,:) = 0.0_rspec
1334 chipss_nc(:,:) = 0.0_rspec
1335 chitss_nc(:,:) = 0.0_rspec
1336 
1337 DO i=1,ms_nc !Over species 1
1338 
1339  im=jms_nc(i)
1340  iz=jzs_nc(i)
1341  iza=iabs(iz)
1342  !Currents
1343 
1344  denzc=den_iz(im,iza)*iz*z_coulomb
1345  pjbbs_nc=pjbbs_nc+denzc*upars_nc(1,1,i)
1346 
1347  pjboh_nc=pjboh_nc+denzc*upars_nc(1,2,i)
1348 
1349  pjbex_nc=pjbex_nc+denzc*upars_nc(1,3,i)
1350  !Unit p'/p and T'/T
1351  bsjbps_nc(:)=bsjbps_nc(:)+denzc*uaip(1,:,i)
1352  bsjbts_nc(:)=bsjbts_nc(:)+denzc*uait(1,:,i)
1353  !Fluxes
1354  !Banana-Plateau
1355  cbpa=cbp/iz
1356  cbpaq=cbpa*(z_j7kv*temp_i(im))
1357  !Unit p'/p and T'/T
1358  dpss_nc(i,i)=dpss_nc(i,i)-cbpa*ymus_nc(1,1,i)*srcthp(i)
1359  dtss_nc(i,i)=dtss_nc(i,i)-cbpa*ymus_nc(1,2,i)*srctht(i)
1360  chipss_nc(i,i)=chipss_nc(i,i)-cbpaq*ymus_nc(2,1,i)*srcthp(i)
1361  chitss_nc(i,i)=chitss_nc(i,i)-cbpaq*ymus_nc(2,2,i)*srctht(i)
1362  !p' and T'
1363  gfls_nc(1,i)=gfls_nc(1,i)-cbpa &
1364  & *sum(ymus_nc(1,:,i)*(upars_nc(:,1,i)+srcth(:,i)))
1365  qfls_nc(1,i)=qfls_nc(1,i)-cbpaq &
1366  & *sum(ymus_nc(2,:,i)*(upars_nc(:,1,i)+srcth(:,i)))
1367 
1368  gfls_nc(4,i)=gfls_nc(4,i)-cbpa &
1369  & *sum(ymus_nc(1,:,i)*upars_nc(:,2,i))
1370  qfls_nc(4,i)=qfls_nc(4,i)-cbpaq &
1371  & *sum(ymus_nc(2,:,i)*upars_nc(:,2,i))
1372  ! write(6,*) 'qfls_nc(4,i)=',qfls_nc(4,i),i
1373  ! write(6,*) 'flow, upars_nc(1,2,i)=',upars_nc(1,2,i)
1374  ! External force
1375  gfls_nc(5,i)=gfls_nc(5,i)-cbpa &
1376  & *sum(ymus_nc(1,:,i)*upars_nc(:,3,i))
1377  qfls_nc(5,i)=qfls_nc(5,i)-cbpaq &
1378  & *sum(ymus_nc(2,:,i)*upars_nc(:,3,i))
1379 
1380  ! write(6,*) 'qfls_nc(5,i)=',qfls_nc(5,i)
1381  DO k=1,k_order_nc !Over velocity moments
1382 
1383  !Unit p'/p and T'/T
1384  dpss_nc(:,i)=dpss_nc(:,i)-cbpa*ymus_nc(1,k,i)*uaip(k,:,i)
1385  dtss_nc(:,i)=dtss_nc(:,i)-cbpa*ymus_nc(1,k,i)*uait(k,:,i)
1386  chipss_nc(:,i)=chipss_nc(:,i)-cbpaq*ymus_nc(2,k,i)*uaip(k,:,i)
1387  chitss_nc(:,i)=chitss_nc(:,i)-cbpaq*ymus_nc(2,k,i)*uait(k,:,i)
1388  ! if (ipr .eq. 2) then
1389  ! write(*,*) 'i=',i,'k=',k,' ymus_nc(2,k,i)=',ymus_nc(2,k,i)
1390  ! write(*,*) 'chipss_nc(i,i)=',chipss_nc(i,i)
1391  ! write(*,*) 'ymus_nc(2,k,i)=',ymus_nc(2,k,i)
1392  ! if (i .eq. ms_nc) stop
1393  ! endif
1394 
1395  END DO !Over velocity moments
1396 
1397  !Pfirsch-Schluter and classical
1398  !Test particle
1399  cpsa=cps*(xis_nc(i)/iz)
1400  cpsaq=cpsa*(z_j7kv*temp_i(im))
1401  ccla=ccl*(xis_nc(i)/iz)
1402  cclaq=ccla*(z_j7kv*temp_i(im))
1403  !Pfirsch-Schluter
1404  gfls_nc(2,i)=gfls_nc(2,i)-cpsa &
1405  & *sum(calmi_nc(1,:,im)*srcth(:,i))
1406  qfls_nc(2,i)=qfls_nc(2,i)-cpsaq &
1407  & *sum(calmi_nc(2,:,im)*srcth(:,i))
1408  !Classical
1409  gfls_nc(3,i)=gfls_nc(3,i)+ccla &
1410  & *sum(calmi_nc(1,:,im)*srcth(:,i))
1411  qfls_nc(3,i)=qfls_nc(3,i)+cclaq &
1412  & *sum(calmi_nc(2,:,im)*srcth(:,i))
1413  !Unit p'/p and T'/T
1414  dpss_nc(i,i)=dpss_nc(i,i)-(cpsa-ccla)*calmi_nc(1,1,im)*srcthp(i)
1415  dtss_nc(i,i)=dtss_nc(i,i)-(cpsa-ccla)*calmi_nc(1,2,im)*srctht(i)
1416  ! if (i.eq.1) then
1417  ! write(*,*) 'dpss_nc(1,1)=',dpss_nc(1,1),cpsa-ccla
1418  ! write(*,*) 'dtss_nc(1,1)=',dtss_nc(1,1),calmi_nc(1,1,im)
1419  ! write(*,*) 'calmi_nc(1,2,im)=',calmi_nc(1,2,im),im
1420  ! endif
1421  chipss_nc(i,i)=chipss_nc(i,i)-(cpsaq-cclaq)*calmi_nc(2,1,im) &
1422  & *srcthp(i)
1423  chitss_nc(i,i)=chitss_nc(i,i)-(cpsaq-cclaq)*calmi_nc(2,2,im) &
1424  & *srctht(i)
1425 
1426  ! if (ipr .eq. 2) then
1427  ! write(*,*) 'i=',i,'chipss_nc(i,i)=',chipss_nc(i,i)
1428  ! write(*,*) 'i=',i,'chitss_nc(i,i)=',chitss_nc(i,i)
1429  ! if (i .eq. ms_nc) stop
1430  ! endif
1431  !Field particle
1432  DO j=1,ms_nc !Over species 2
1433 
1434  jm=jms_nc(j)
1435  cpsb=cpsa*xis_nc(j)
1436  cpsbq=cpsb*(z_j7kv*temp_i(im))
1437  cclb=ccla*xis_nc(j)
1438  cclbq=cclb*(z_j7kv*temp_i(im))
1439  !Pfirsch-Schluter
1440  gfls_nc(2,i)=gfls_nc(2,i)-cpsb &
1441  & *sum(calnii_nc(1,1:k_order_nc,im,jm)*srcth(:,j))
1442  qfls_nc(2,i)=qfls_nc(2,i)-cpsbq &
1443  & *sum(calnii_nc(2,1:k_order_nc,im,jm)*srcth(:,j))
1444  !Classical
1445  gfls_nc(3,i)=gfls_nc(3,i)+cclb &
1446  & *sum(calnii_nc(1,1:k_order_nc,im,jm)*srcth(:,j))
1447  qfls_nc(3,i)=qfls_nc(3,i)+cclbq &
1448  & *sum(calnii_nc(2,1:k_order_nc,im,jm)*srcth(:,j))
1449  !Unit p'/p and T'/T
1450  dpss_nc(j,i)=dpss_nc(j,i) &
1451  & -(cpsb-cclb)*calnii_nc(1,1,im,jm)*srcthp(j)
1452  dtss_nc(j,i)=dtss_nc(j,i) &
1453  & -(cpsb-cclb)*calnii_nc(1,2,im,jm)*srctht(j)
1454  chipss_nc(j,i)=chipss_nc(j,i) &
1455  & -(cpsbq-cclbq)*calnii_nc(2,1,im,jm)*srcthp(j)
1456  chitss_nc(j,i)=chitss_nc(j,i) &
1457  & -(cpsbq-cclbq)*calnii_nc(2,2,im,jm)*srctht(j)
1458 
1459  ENDDO !Over species 2
1460 
1461 ENDDO !Over species 1
1462 
1463 !Electrical resistivity
1464 petap_nc=1.0_rspec/pjboh_nc
1465 ! if (ipr .eq. 83) then
1466 ! write(*,*) 'ipr=',ipr,'petap_nc =',petap_nc,1.0_rspec,pjboh_nc
1467 ! endif
1468 !Add <E.B> normalization
1469 pjboh_nc=p_eb*pjboh_nc
1470 upars_nc(:,2,:)=p_eb*upars_nc(:,2,:)
1471 qfls_nc(4,:) = p_eb*qfls_nc(4,:)
1472 gfls_nc(4,:) = p_eb*gfls_nc(4,:)
1473 
1474 !-------------------------------------------------------------------------------
1475 !Convert full coefficient matrices to diffusivities and conductivities
1476 !-------------------------------------------------------------------------------
1477 DO i=1,ms_nc !Over species
1478 
1479  im=jms_nc(i)
1480  iza=iabs(jzs_nc(i))
1481  dpss_nc(:,i)=dpss_nc(:,i)/den_iz(im,iza)
1482  dtss_nc(:,i)=dtss_nc(:,i)/den_iz(im,iza)
1483  chipss_nc(:,i)=chipss_nc(:,i)/den_iz(im,iza)/temp_i(im)/z_j7kv
1484  chitss_nc(:,i)=chitss_nc(:,i)/den_iz(im,iza)/temp_i(im)/z_j7kv
1485 
1486 ENDDO !Over species
1487 
1488 9999 CONTINUE
1489 
1490 END SUBROUTINE nclass_flow
1491 
1492 SUBROUTINE nclass_init(m_i,m_z,amu_i,den_iz,iflag)
1493  !-------------------------------------------------------------------------------
1494  !NCLASS_INIT initializes the species information and allocates arrays
1495  !References:
1496  ! W.A.Houlberg 1/2002
1497  !Input:
1498  ! m_i -number of isotopes (> 1) [-]
1499  ! m_z -highest charge state [-]
1500  ! amu_i(i) -atomic mass number [-]
1501  ! den_iz(i,z) -density [/m**3]
1502  !Output:
1503  ! iflag -error and warning flag [-]
1504  ! =-1 warning
1505  ! =0 no warnings or errors
1506  ! =1 error
1507  !-------------------------------------------------------------------------------
1508 implicit none
1509  !Declaration of input variables
1510 INTEGER, INTENT(IN) :: &
1511  & m_i,m_z
1512 
1513 REAL(KIND=rspec), INTENT(IN) :: &
1514  & amu_i(:), &
1515  & den_iz(:,:)
1516 
1517 !Declaration of output variables
1518 INTEGER, INTENT(OUT) :: &
1519  & iflag
1520 
1521 
1522 !Declaration of local variables
1523 INTEGER :: &
1524  & i,j,k, &
1525  & nset1,nset2
1526 
1527 !-------------------------------------------------------------------------------
1528 !Get number of significant isotopes and charge states, and max charge
1529 !-------------------------------------------------------------------------------
1530 mi_nc=0
1531 mz_nc=0
1532 ms_nc=0
1533 
1534 DO i=1,m_i !Over isotopes
1535 
1536  DO j=1,m_z !Over charge states
1537  IF(den_iz(i,j) >= cden_nc) THEN
1538 
1539  ms_nc=ms_nc+1
1540  IF(i > mi_nc) mi_nc=i
1541  IF(j > mz_nc) mz_nc=j
1542 
1543  ENDIF
1544 
1545  ! write(*,*) 'i=',i,' j=',j,' m_i=',m_i,' m_z=',m_z
1546  ! write(*,*) 'ms_nc=',ms_nc,' mi_nc=',mi_nc
1547  ! write(*,*) 'cden_nc =',cden_nc,' den_iz=',den_iz(i,j)
1548 
1549  END DO !Over charge states
1550 
1551 END DO !Over isotopes
1552 
1553 !At least two species
1554 IF(mi_nc < 2) THEN
1555 
1556  iflag=1
1557  write(*,*) 'NCLASS_INIT(1)/ERROR:m_i must be >= 2'
1558  write(*,*) 'mi_nc =',mi_nc
1559  goto 9999
1560 
1561 ENDIF
1562 
1563 !Highest charge state at least 1
1564 IF(mz_nc < 1) THEN
1565 
1566  iflag=1
1567  write(*,*) 'NCLASS_INIT(2)/ERROR:m_z must be >= 1'
1568  goto 9999
1569 
1570 ENDIF
1571 
1572 !-------------------------------------------------------------------------------
1573 !Allocate or reallocate isotope info
1574 !-------------------------------------------------------------------------------
1575 IF(ALLOCATED(calmi_nc)) THEN
1576 
1577  nset1=SIZE(calmi_nc,1)
1578  nset2=SIZE(calmi_nc,3)
1579 
1580  !If storage requirements have changed, reallocate
1581  IF(nset1 /= k_order_nc .OR. nset2 /= mi_nc) THEN
1582 
1583  !Reallocate
1584  DEALLOCATE(vti_nc, &
1585  & amntii_nc, &
1586  & calmi_nc, &
1587  & calnii_nc, &
1588  & capmii_nc, &
1589  & capnii_nc)
1590  ALLOCATE(vti_nc(1:mi_nc), &
1591  & amntii_nc(1:mi_nc,1:mi_nc), &
1592  & calmi_nc(1:k_order_nc,1:k_order_nc,1:mi_nc), &
1593  & calnii_nc(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc), &
1594  & capmii_nc(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc), &
1595  & capnii_nc(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc))
1596 
1597  ENDIF
1598 
1599 ELSE
1600 
1601  !Allocate
1602  ALLOCATE(vti_nc(1:mi_nc), &
1603  & amntii_nc(1:mi_nc,1:mi_nc), &
1604  & calmi_nc(1:k_order_nc,1:k_order_nc,1:mi_nc), &
1605  & calnii_nc(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc), &
1606  & capmii_nc(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc), &
1607  & capnii_nc(1:k_order_nc,1:k_order_nc,1:mi_nc,1:mi_nc))
1608 
1609 ENDIF
1610 
1611 !-------------------------------------------------------------------------------
1612 !Allocate or reallocate isotope/charge state info
1613 !-------------------------------------------------------------------------------
1614 IF(ALLOCATED(grppiz_nc)) THEN
1615 
1616  nset1=SIZE(grppiz_nc,1)
1617  nset2=SIZE(grppiz_nc,2)
1618 
1619  !If storage requirements have changed, reallocate
1620  IF(nset1 /= mi_nc .OR. nset2 /= mz_nc) THEN
1621 
1622  !Reallocate
1623  DEALLOCATE(grppiz_nc)
1624  ALLOCATE(grppiz_nc(1:mi_nc,1:mz_nc))
1625 
1626  ENDIF
1627 
1628 ELSE
1629 
1630  !Allocate
1631  ALLOCATE(grppiz_nc(1:mi_nc,1:mz_nc))
1632 
1633 ENDIF
1634 
1635 !-------------------------------------------------------------------------------
1636 !Allocate or reallocate species info
1637 !-------------------------------------------------------------------------------
1638 IF(ALLOCATED(jms_nc)) THEN
1639 
1640  nset1=SIZE(jms_nc)
1641 
1642  !If storage requirements have changed, reallocate
1643  IF(nset1 /= ms_nc) THEN
1644 
1645  !Reallocate
1646  DEALLOCATE(jms_nc, &
1647  & jzs_nc, &
1648  & sqzs_nc, &
1649  & xis_nc, &
1650  & bsjbps_nc, &
1651  & bsjbts_nc, &
1652  & gfls_nc, &
1653  & dpss_nc, &
1654  & dtss_nc, &
1655  & qfls_nc, &
1656  & chipss_nc, &
1657  & chitss_nc, &
1658  & tauss_nc)
1659 
1660  ALLOCATE(jms_nc(1:ms_nc), &
1661  & jzs_nc(1:ms_nc), &
1662  & sqzs_nc(1:ms_nc), &
1663  & xis_nc(1:ms_nc), &
1664  & bsjbps_nc(1:ms_nc), &
1665  & bsjbts_nc(1:ms_nc), &
1666  & gfls_nc(5,1:ms_nc), &
1667  & dpss_nc(1:ms_nc,1:ms_nc), &
1668  & dtss_nc(1:ms_nc,1:ms_nc), &
1669  & qfls_nc(5,1:ms_nc), &
1670  & chipss_nc(1:ms_nc,1:ms_nc), &
1671  & chitss_nc(1:ms_nc,1:ms_nc), &
1672  & tauss_nc(1:ms_nc,1:ms_nc))
1673 
1674  ENDIF
1675 
1676 ELSE
1677 
1678  !Allocate
1679  ALLOCATE(jms_nc(1:ms_nc), &
1680  & jzs_nc(1:ms_nc), &
1681  & sqzs_nc(1:ms_nc), &
1682  & xis_nc(1:ms_nc), &
1683  & bsjbps_nc(1:ms_nc), &
1684  & bsjbts_nc(1:ms_nc), &
1685  & gfls_nc(5,1:ms_nc), &
1686  & dpss_nc(1:ms_nc,1:ms_nc), &
1687  & dtss_nc(1:ms_nc,1:ms_nc), &
1688  & qfls_nc(5,1:ms_nc), &
1689  & chipss_nc(1:ms_nc,1:ms_nc), &
1690  & chitss_nc(1:ms_nc,1:ms_nc), &
1691  & tauss_nc(1:ms_nc,1:ms_nc))
1692 
1693 
1694 ENDIF
1695 
1696 !-------------------------------------------------------------------------------
1697 !Allocate or reallocate species/moments info
1698 !-------------------------------------------------------------------------------
1699 IF(ALLOCATED(ymus_nc)) THEN
1700 
1701  nset1=SIZE(ymus_nc,1)
1702  nset2=SIZE(ymus_nc,3)
1703 
1704  !If storage requirements have changed, reallocate
1705  IF(nset1 /= ms_nc .OR. nset2 /= k_order_nc) THEN
1706 
1707  !Reallocate
1708  DEALLOCATE(ymus_nc, &
1709  & upars_nc, &
1710  & uthetas_nc)
1711 
1712  ALLOCATE(ymus_nc(1:k_order_nc,1:k_order_nc,1:ms_nc), &
1713  & upars_nc(1:k_order_nc,1:5,1:ms_nc), &
1714  & uthetas_nc(1:k_order_nc,1:5,1:ms_nc))
1715 
1716  ENDIF
1717 
1718 ELSE
1719 
1720  !Allocate
1721  ALLOCATE(ymus_nc(1:k_order_nc,1:k_order_nc,1:ms_nc), &
1722  & upars_nc(1:k_order_nc,1:5,1:ms_nc), &
1723  & uthetas_nc(1:k_order_nc,1:5,1:ms_nc))
1724 
1725 ENDIF
1726 
1727 !-------------------------------------------------------------------------------
1728 !Find significant charge states and mapping
1729 !-------------------------------------------------------------------------------
1730 k=0
1731 imel_nc=0
1732 
1733 DO i=1,mi_nc !Over isotopes
1734 
1735  DO j=1,mz_nc !Over charge states
1736 
1737  IF(den_iz(i,j) >= cden_nc) THEN
1738 
1739  ! Set isotope number and charge state for this species
1740  k=k+1
1741  jms_nc(k)=i
1742 
1743  IF(amu_i(i) < 0.5) THEN
1744 
1745  !Electrons
1746  jzs_nc(k)=-j
1747  imel_nc=i
1748 
1749  ELSE
1750 
1751  !Ions
1752  jzs_nc(k)=j
1753 
1754  ENDIF
1755 
1756  ENDIF
1757 
1758  ENDDO !Over charge states
1759 
1760 ENDDO !Over isotopes
1761 
1762 9999 CONTINUE
1763 
1764 END SUBROUTINE nclass_init
1765 
1766 SUBROUTINE nclass_k(p_fm,p_ft,p_ngrth,x,amu_i,temp_i,ykb_s,ykp_s, &
1767  & ykpo_s,ykpop_s)
1768  !-------------------------------------------------------------------------------
1769  !NCLASS_K calculates the velocity-dependent neoclassical viscosity coefficients
1770  !References:
1771  ! S.P.Hirshman, D.J.Sigmar, Nucl Fusion 21 (1981) 1079
1772  ! C.E.Kessel, Nucl Fusion 34 (1994) 1221
1773  ! K.C.Shaing, M.Yokoyama, M.Wakatani, C.T.Hsu, Phys Plasmas 3 (1996) 965
1774  ! W.A.Houlberg, K.C.Shaing, S.P.Hirshman, M.C.Zarnstorff, Phys Plasmas 4 (1997)
1775  ! 3230
1776  ! W.A.Houlberg 1/2002
1777  !Input:
1778  ! p_fm(:) -poloidal moments of drift factor for PS [/m**2]
1779  ! p_ft -trapped fraction [-]
1780  ! p_ngrth -<n.grad(Theta)> [/m]
1781  ! x -normalized velocity v/(2kT/m)**0.5 [-]
1782  ! amu_i(i) -atomic mass number [-]
1783  ! temp_i(i) -temperature [keV]
1784  !Output:
1785  ! ykb_s(s) -banana viscosity for s [/s]
1786  ! ykp_s(s) -Pfirsch-Schluter viscosity for s [/s]
1787  ! ykpo_s(s) -potato viscosity for s [/s]
1788  ! ykpop_s(s) -potato-plateau viscosity for s [/s]
1789  !-------------------------------------------------------------------------------
1790 implicit none
1791  !Declaration of input variables
1792 REAL(KIND=rspec), INTENT(IN) :: &
1793  & p_fm(:),p_ft,p_ngrth,x, &
1794  & amu_i(:),temp_i(:)
1795 
1796 !Declaration of output variables
1797 REAL(KIND=rspec) :: &
1798  & ykb_s(:),ykp_s(:),ykpo_s(:),ykpop_s(:)
1799 
1800 !Declaration of local variables
1801 INTEGER :: &
1802  & i,im,iz
1803 
1804 REAL(KIND=rspec) :: &
1805  & c1,c2,c3,c4
1806 
1807 REAL(KIND=rspec) :: &
1808  & ynud_s(ms_nc),ynut_s(ms_nc),ynutis(mf_nc,ms_nc)
1809 
1810 !-------------------------------------------------------------------------------
1811 !Initialization
1812 !-------------------------------------------------------------------------------
1813 ykb_s(:)=0.0_rspec
1814 ykp_s(:)=0.0_rspec
1815 ykpo_s(:)=0.0_rspec
1816 ykpop_s(:)=0.0_rspec
1817 
1818 !Get collisional frequencies
1819 CALL nclass_nu(p_ngrth,x,temp_i,ynud_s,ynut_s,ynutis)
1820 
1821 !-------------------------------------------------------------------------------
1822 !Set velocity dependent viscosities (K's)
1823 !-------------------------------------------------------------------------------
1824 c1=1.5_rspec*x**2
1825 c2=3.0_rspec*2.19_rspec/(2.0_rspec**1.5)*x**(1.0/3.0)
1826 
1827 IF(l_potato_nc) c3=3.0_rspec*z_pi/(64.0_rspec*2.0_rspec**0.33333) &
1828  & /abs(cpotl_nc)
1829 
1830 DO i=1,ms_nc !Over species
1831 
1832  im=jms_nc(i)
1833  iz=jzs_nc(i)
1834 
1835  !Banana
1836  IF(l_banana_nc) THEN
1837 
1838  !Provide cutoff to eliminate failure at unity trapped fraction
1839  !At A=>1 viscosity will go over to Pfirsch-Schluter value
1840  c4=1.0_rspec-p_ft
1841  IF(c4 < 1.0e-3_rspec) c4=1.0e-3_rspec
1842  ykb_s(i)=p_ft/c4/sqzs_nc(i)**1.5*ynud_s(i)
1843 
1844  ENDIF
1845 
1846  !Pfirsch-Schuter
1847  IF(l_pfirsch_nc) THEN
1848 
1849  ykp_s(i)=ykp_s(i)+c1*vti_nc(im)**2 &
1850  & *sum(p_fm(:)*ynutis(:,i))/ynut_s(i)
1851 
1852  ENDIF
1853 
1854  !Potato
1855  IF(l_potato_nc) THEN
1856 
1857  c4=abs(amu_i(im)*z_protonmass*vti_nc(im) &
1858  & /(iz*z_coulomb*cpotb_nc*cpotl_nc))
1859  ykpo_s(i)=c2*c4**(0.333333_rspec)*ynud_s(i) &
1860  & /sqzs_nc(i)**(1.666667_rspec)
1861  ykpop_s(i)=c3*vti_nc(im)*c4**(1.333333_rspec)
1862 
1863  ENDIF
1864 
1865 ENDDO !Over species
1866 
1867 END SUBROUTINE nclass_k
1868 
1869 SUBROUTINE nclass_mn(amu_i,temp_i)
1870  !-------------------------------------------------------------------------------
1871  !NCLASS_MN calculates neoclassical friction coefficients
1872  !References:
1873  ! S.P.Hirshman, D.J.Sigmar, Nucl Fusion 21 (1981) 1079
1874  ! S.P.Hirshman, Phys Fluids 20 (1977) 589
1875  ! W.A.Houlberg, K.C.Shaing, S.P.Hirshman, M.C.Zarnstorff, Phys Plasmas 4 (1997)
1876  ! 3230
1877  ! W.A.Houlberg 1/2002
1878  !Input:
1879  ! amu_i(i) -atomic mass number [-]
1880  ! temp_i(i) -temperature [keV]
1881  !Comments:
1882  ! The k_order*korder matrix of test particle (M) and field particle (N)
1883  ! coefficients of the collision operator use the Laguerre polynomials of
1884  ! order 3/2 as basis functions for each isotopic combination
1885  ! The indices on the M and N matrices are one greater than the notation in the
1886  ! H&S review article so as to avoid 0 as an index
1887  !-------------------------------------------------------------------------------
1888 implicit none
1889  !Declaration of input variables
1890 REAL(KIND=rspec), INTENT(IN) :: &
1891  & amu_i(:),temp_i(:)
1892 
1893 !Declaration of local variables
1894 INTEGER :: &
1895  & im,jm
1896 
1897 REAL(KIND=rspec) :: &
1898  & xab,xab2,xmab,xtab,yab32,yab52,yab72,yab92
1899 
1900 !-------------------------------------------------------------------------------
1901 !Calculate friction coefficients
1902 !-------------------------------------------------------------------------------
1903 DO im=1,mi_nc !Over isotopes a
1904 
1905  DO jm=1,mi_nc !Over isotopes b
1906 
1907  !Mass ratio
1908  xmab=amu_i(im)/amu_i(jm)
1909  !Temperature ratio
1910  xtab=temp_i(im)/temp_i(jm)
1911  !Thermal velocity ratio, vtb/vta
1912  xab=sqrt(xmab/xtab)
1913 
1914  xab2=xab**2
1915  yab32=(1.0_rspec+xab2)*sqrt(1.0_rspec+xab2)
1916  yab52=(1.0_rspec+xab2)*yab32
1917 
1918  IF(k_order_nc == 3) THEN
1919 
1920  yab72=(1.0_rspec+xab2)*yab52
1921  yab92=(1.0_rspec+xab2)*yab72
1922 
1923  ENDIF
1924 
1925  !Test particle coefficients, M
1926  !Eqn 4.11 (HS81) for M00
1927  capmii_nc(1,1,im,jm)=-(1.0_rspec+xmab)/yab32
1928  !Eqn 4.12 (HS81) for M01
1929  capmii_nc(1,2,im,jm)=1.5_rspec*(1.0_rspec+xmab)/yab52
1930  !Eqn 4.8 (HS81) for M10
1931  capmii_nc(2,1,im,jm)=capmii_nc(1,2,im,jm)
1932  !Eqn 4.13 (HS81) for M11
1933  capmii_nc(2,2,im,jm)=-(3.25_rspec+xab2 &
1934  & *(4.0_rspec+7.5_rspec*xab2))/yab52
1935 
1936  IF(k_order_nc == 3) THEN
1937 
1938  !Eqn 4.15 (HS81) for M02
1939  capmii_nc(1,3,im,jm)=-1.875_rspec*(1.0_rspec+xmab)/yab72
1940  !Eqn 4.16 (HS81) for M12
1941  capmii_nc(2,3,im,jm)=(4.3125_rspec+xab2 &
1942  & *(6.0_rspec+15.75_rspec*xab2))/yab72
1943  !Eqn 4.8 (HS81) for M20
1944  capmii_nc(3,1,im,jm)=capmii_nc(1,3,im,jm)
1945  !Eqn 4.8 (HS81) for M21
1946  capmii_nc(3,2,im,jm)=capmii_nc(2,3,im,jm)
1947  !Eqn 5.21 (HS81) for M22
1948  capmii_nc(3,3,im,jm)=-(433.0/64.0+xab2 &
1949  & *(17.0_rspec+xab2 &
1950  & *(57.375_rspec+xab2 &
1951  & *(28.0_rspec+xab2 &
1952  & *21.875_rspec))))/yab92
1953 
1954  ENDIF
1955 
1956  !Field particle coefficients, N
1957  !Momentum conservation, Eqn 4.11 (HS81) for N00
1958  capnii_nc(1,1,im,jm)=-capmii_nc(1,1,im,jm)
1959  !Eqn 4.9 and 4.12 (HS81) for N01
1960  capnii_nc(1,2,im,jm)=-xab2*capmii_nc(1,2,im,jm)
1961  !Momentum conservation, Eqn 4.12 (HS81) for N10
1962  capnii_nc(2,1,im,jm)=-capmii_nc(2,1,im,jm)
1963  !Eqn 4.14 (HS81) for N11 - corrected rhs
1964  capnii_nc(2,2,im,jm)=6.75_rspec*sqrt(xtab)*xab2/yab52
1965 
1966  IF(k_order_nc == 3) THEN
1967 
1968  !Eqn 4.15 for N02 (HS81) - corrected rhs by Ta/Tb
1969  capnii_nc(1,3,im,jm)=-xab2**2*capmii_nc(1,3,im,jm)
1970  !Eqn 4.17 for N12 (HS81)
1971  capnii_nc(2,3,im,jm)=-14.0625_rspec*xtab*xab2**2/yab72
1972  !Momentum conservation for N20 (HS81)
1973  capnii_nc(3,1,im,jm)=-capmii_nc(3,1,im,jm)
1974  !Eqn 4.9 and 4.17 for N21 (HS81)
1975  capnii_nc(3,2,im,jm)=-14.0625_rspec*xab2**2/yab72
1976  !Eqn 5.22 for N22 (HS81)
1977  capnii_nc(3,3,im,jm)=2625.0/64.0*xtab*xab2**2/yab92
1978 
1979  ENDIF
1980 
1981  ENDDO !Over isotopes b
1982 
1983 ENDDO !Over isotopes a
1984 
1985 END SUBROUTINE nclass_mn
1986 
1987 SUBROUTINE nclass_mu(p_fm,p_ft,p_ngrth,amu_i,temp_i,den_iz,ipr)
1988  !-------------------------------------------------------------------------------
1989  !NCLASS_MU calculates the matrix of neoclassical fluid viscosities
1990  !References:
1991  ! K.C.Shaing, M.Yokoyama, M.Wakatani, C.T.Hsu, Phys Plasmas 3 (1996) 965
1992  ! W.A.Houlberg, K.C.Shaing, S.P.Hirshman, M.C.Zarnstorff, Phys Plasmas 4 (1997)
1993  ! 3230
1994  ! W.A.Houlberg 1/2002
1995  !Input:
1996  ! p_fm(:) -poloidal moments of drift factor for PS [/m**2]
1997  ! p_ft -trapped fraction [-]
1998  ! p_ngrth -<n.grad(Theta)> [/m]
1999  ! amu_i(i) -atomic mass number [-]
2000  ! temp_i(i) -temperature [keV]
2001  ! den_iz(i,z) -density [/m**3]
2002  !Comments:
2003  ! Integrates the velocity-dependent banana and Pfirsch-Schluter contributions
2004  ! over velocity space
2005  !-------------------------------------------------------------------------------
2006 implicit none
2007  !Declaration of input variables
2008 REAL(KIND=rspec), INTENT(IN) :: &
2009  & p_ft,p_ngrth,p_fm(:), &
2010  & amu_i(:),temp_i(:), &
2011  & den_iz(:,:)
2012 
2013 !Declaration of local variables
2014 !Parameters
2015 INTEGER, PARAMETER :: &
2016  & mpnts=13
2017 
2018 REAL(KIND=rspec), PARAMETER :: &
2019  & bmax=3.2
2020 
2021 !Saved
2022 INTEGER, SAVE :: &
2023  & init=0
2024 
2025 REAL(KIND=rspec), SAVE :: &
2026  & c1,h, &
2027  & x(mpnts),w(mpnts,5)
2028 
2029 !Other
2030 INTEGER :: &
2031  & i,im,iza,k,l,m,ipr
2032 
2033 REAL(KIND=rspec) :: &
2034  & c2,ewt,expmx2,x2,x4,x6,x8,x10,x12,xk,xx, &
2035  & dum(3), &
2036  & ykb_s(ms_nc),ykp_s(ms_nc),ykpo_s(ms_nc),ykpop_s(ms_nc), &
2037  & ymubs(1:k_order_nc,1:k_order_nc,ms_nc), &
2038  & ymubps(1:k_order_nc,1:k_order_nc,ms_nc), &
2039  & ymupps(1:k_order_nc,1:k_order_nc,ms_nc), &
2040  & ymupos(1:k_order_nc,1:k_order_nc,ms_nc)
2041 
2042 !-------------------------------------------------------------------------------
2043 !Initialization
2044 !-------------------------------------------------------------------------------
2045 IF(init == 0) THEN
2046 
2047  !Set integration points and weights
2048  h=bmax/(mpnts-1)
2049  x(1)=0.0_rspec
2050  w(:,:)=0.0_rspec
2051 
2052  DO m=2,mpnts !Over velocity nodes for integration
2053 
2054  x(m)=h*(m-1)
2055  x2=x(m)*x(m)
2056  expmx2=exp(-x2)
2057  x4=x2*x2
2058  w(m,1)=x4*expmx2
2059  x6=x4*x2
2060  w(m,2)=x6*expmx2
2061  x8=x4*x4
2062  w(m,3)=x8*expmx2
2063  x10=x4*x6
2064  w(m,4)=x10*expmx2
2065  x12=x6*x6
2066  w(m,5)=x12*expmx2
2067 
2068  ENDDO !Over velocity nodes for integration
2069 
2070  c1=8.0_rspec/3.0_rspec/sqrt(z_pi)*h
2071  init=1
2072 
2073 ENDIF
2074 
2075 !Null out viscosity arrays
2076 ymus_nc(:,:,:)=0.0_rspec
2077 ymubs(:,:,:)=0.0_rspec
2078 ymubps(:,:,:)=0.0_rspec
2079 ymupos(:,:,:)=0.0_rspec
2080 ymupps(:,:,:)=0.0_rspec
2081 
2082 !-------------------------------------------------------------------------------
2083 !Integrate over velocity space for fluid viscosities
2084 !-------------------------------------------------------------------------------
2085 IF(l_banana_nc .OR. l_pfirsch_nc) THEN
2086 
2087  DO m=2,mpnts !Over velocity grid
2088 
2089  IF(m == mpnts) THEN
2090 
2091  !Use half weight for end point
2092  ewt=0.5_rspec
2093 
2094  ELSE
2095 
2096  !Use full weight for intervening points
2097  ewt=1.0_rspec
2098 
2099  ENDIF
2100 
2101  xx=x(m)
2102  !Get velocity-dependent k values
2103  CALL nclass_k(p_fm,p_ft,p_ngrth,xx,amu_i,temp_i,ykb_s,ykp_s, &
2104  & ykpo_s,ykpop_s)
2105 
2106  DO i=1,ms_nc !Over species
2107 
2108  im=jms_nc(i)
2109  iza=iabs(jzs_nc(i))
2110  c2=c1*ewt*den_iz(im,iza)*amu_i(im)*z_protonmass
2111  dum(1)=c2*w(m,1)
2112  dum(2)=c2*(w(m,2)-2.5_rspec*w(m,1))
2113  dum(3)=c2*(w(m,3)-5.0_rspec*w(m,2)+6.25_rspec*w(m,1))
2114 
2115  IF(.NOT. l_banana_nc) THEN
2116 
2117  !Only Pfirsch-Schluter
2118  xk=ykp_s(i)
2119 
2120  ELSEIF(.NOT. l_pfirsch_nc) THEN
2121 
2122  !Only banana
2123  xk=ykb_s(i)
2124 
2125  ELSE
2126 
2127  !Both banana and Pfirsch-Schluter
2128  xk=ykb_s(i)*ykp_s(i)/(ykb_s(i)+ykp_s(i))
2129 
2130  ENDIF
2131 
2132  ymubs(1,1,i)=ymubs(1,1,i)+ykb_s(i)*dum(1)
2133  ymubs(1,2,i)=ymubs(1,2,i)+ykb_s(i)*dum(2)
2134  ymubs(2,2,i)=ymubs(2,2,i)+ykb_s(i)*dum(3)
2135  ymubps(1,1,i)=ymubps(1,1,i)+xk*dum(1)
2136  ymubps(1,2,i)=ymubps(1,2,i)+xk*dum(2)
2137  ymubps(2,2,i)=ymubps(2,2,i)+xk*dum(3)
2138 
2139  IF(l_potato_nc) THEN
2140 
2141  ymupos(1,1,i)=ymupos(1,1,i)+ykpo_s(i)*dum(1)
2142  ymupos(1,2,i)=ymupos(1,2,i)+ykpo_s(i)*dum(2)
2143  ymupos(2,2,i)=ymupos(2,2,i)+ykpo_s(i)*dum(3)
2144  xk=ykpo_s(i)*ykpop_s(i)/(ykpo_s(i)+ykpop_s(i))
2145  ymupps(1,1,i)=ymupps(1,1,i)+xk*dum(1)
2146  ymupps(1,2,i)=ymupps(1,2,i)+xk*dum(2)
2147  ymupps(2,2,i)=ymupps(2,2,i)+xk*dum(3)
2148 
2149  ENDIF
2150 
2151  IF(k_order_nc == 3) THEN
2152 
2153  dum(1)=c2*(0.5_rspec*w(m,3)-3.5_rspec*w(m,2) &
2154  & +4.375_rspec*w(m,1))
2155  dum(2)=c2*(0.5_rspec*w(m,4)-4.75_rspec*w(m,3) &
2156  & +13.125_rspec*w(m,2)-10.9375_rspec*w(m,1))
2157  dum(3)=c2*(0.25_rspec*w(m,5)-3.5_rspec*w(m,4) &
2158  & +16.625_rspec*w(m,3)-30.625*w(m,2) &
2159  & +(1225.0/64.0)*w(m,1))
2160  ymubs(1,3,i)=ymubs(1,3,i)+ykb_s(i)*dum(1)
2161  ymubs(2,3,i)=ymubs(2,3,i)+ykb_s(i)*dum(2)
2162  ymubs(3,3,i)=ymubs(3,3,i)+ykb_s(i)*dum(3)
2163  xk=ykb_s(i)*ykp_s(i)/(ykb_s(i)+ykp_s(i))
2164  ymubps(1,3,i)=ymubps(1,3,i)+xk*dum(1)
2165  ymubps(2,3,i)=ymubps(2,3,i)+xk*dum(2)
2166  ymubps(3,3,i)=ymubps(3,3,i)+xk*dum(3)
2167 
2168  IF(l_potato_nc) THEN
2169 
2170  ymupos(1,3,i)=ymupos(1,3,i)+ykpo_s(i)*dum(1)
2171  ymupos(2,3,i)=ymupos(2,3,i)+ykpo_s(i)*dum(2)
2172  ymupos(3,3,i)=ymupos(3,3,i)+ykpo_s(i)*dum(3)
2173  xk=ykpo_s(i)*ykpop_s(i)/(ykpo_s(i)+ykpop_s(i))
2174  ymupps(1,3,i)=ymupps(1,3,i)+xk*dum(1)
2175  ymupps(2,3,i)=ymupps(2,3,i)+xk*dum(2)
2176  ymupps(3,3,i)=ymupps(3,3,i)+xk*dum(3)
2177 
2178  ENDIF
2179 
2180  ENDIF
2181 
2182  ENDDO !Over species
2183 
2184  ENDDO !Over velocity grid
2185 
2186  !Load net viscosity
2187  DO i=1,ms_nc !Over species
2188 
2189  DO l=1,k_order_nc !Over rows
2190 
2191  DO k=1,l !Over columns in upper part of matrix
2192 
2193  IF(l_potato_nc) THEN
2194 
2195  !Banana Pfirsch-Schluter plus potato potato-plateau
2196  ymus_nc(k,l,i)=(ymupos(k,l,i)**3*ymupps(k,l,i) &
2197  & +ymubs(k,l,i)**3*ymubps(k,l,i)) &
2198  & /(ymupos(k,l,i)**3+ymubs(k,l,i)**3)
2199  ELSE
2200 
2201  !Banana Pfirsch-Schluter
2202  ymus_nc(k,l,i)=ymubps(k,l,i)
2203 
2204  ENDIF
2205  ! if (ipr .eq. 2 .and. k .eq. 2) then
2206  ! write(*,*) 'ymubps(2,l,i)=',ymubps(k,l,i),l,i
2207  ! write(*,*) 'ymupos(2,l,i)=',ymupos(k,l,i)
2208  ! write(*,*) 'ymupps(2,l,i)=',ymupps(k,l,i)
2209  ! write(*,*) 'ymubs(2,l,i)=',ymubs(k,l,i),l_potato_nc
2210  ! stop
2211  ! endif
2212 
2213  ENDDO !Over columns in upper part of matrix
2214 
2215  ENDDO !Over rows
2216 
2217  ENDDO !Over species
2218 
2219  !Fill viscosity matrix using symmetry
2220  DO i=1,ms_nc !Over species
2221 
2222  DO l=1,k_order_nc-1 !Over rows
2223 
2224  DO k=l+1,k_order_nc !Over columns
2225 
2226  ymus_nc(k,l,i)=ymus_nc(l,k,i)
2227 
2228  ENDDO !Over columns
2229 
2230  ENDDO !Over rows
2231 
2232  ENDDO !Over species
2233 
2234 ENDIF
2235 
2236 END SUBROUTINE nclass_mu
2237 
2238 SUBROUTINE nclass_nu(p_ngrth,x,temp_i,ynud_s,ynut_s,ynutis)
2239  !-------------------------------------------------------------------------------
2240  !NCLASS_NU calculates the velocity dependent pitch angle diffusion and
2241  ! anisotropy relaxation rates, nu_D, nu_T, and nu_T*I_Rm
2242  !References:
2243  ! S.P.Hirshman, D.J.Sigmar, Phys Fluids 19 (1976) 1532
2244  ! S.P.Hirshman, D.J.Sigmar, Nucl Fusion 21 (1981) 1079
2245  ! K.C.Shaing, M.Yokoyama, M.Wakatani, C.T.Hsu, Phys Plasmas 3 (1996) 965
2246  ! W.A.Houlberg, K.C.Shaing, S.P.Hirshman, M.C.Zarnstorff, Phys Plasmas 4 (1997)
2247  ! 3230
2248  ! W.A.Houlberg 1/2002
2249  !Input:
2250  ! p_ngrth -<n.grad(Theta)> [/m]
2251  ! x -normalized velocity v/(2kT/m)**0.5 [-]
2252  ! temp_i(i) -temperature [keV]
2253  !Output:
2254  ! ynud_s(s) -pitch angle diffusion rate for s [/s]
2255  ! ynut_s(s) -anisotropy relaxation rate for s [/s]
2256  ! ynutis(:,s) -PS anisotropy relaxation rates for s [-]
2257  !-------------------------------------------------------------------------------
2258 implicit none
2259  !Declaration of input variables
2260 REAL(KIND=rspec), INTENT(IN) :: &
2261  & p_ngrth,x, &
2262  & temp_i(:)
2263 
2264 !Declaration of output variables
2265 REAL(KIND=rspec) :: &
2266  & ynud_s(:),ynut_s(:),ynutis(:,:)
2267 
2268 !Declaration of local variables
2269 INTEGER :: &
2270  & i,im,j,jm,m
2271 
2272 REAL(KIND=rspec) :: &
2273  & c1,c2,c3,g,phi
2274 
2275 !-------------------------------------------------------------------------------
2276 !Initializaton
2277 !-------------------------------------------------------------------------------
2278 ynud_s(:)=0.0_rspec
2279 ynut_s(:)=0.0_rspec
2280 ynutis(:,:)=0.0_rspec
2281 
2282 !-------------------------------------------------------------------------------
2283 !Collision rates
2284 !-------------------------------------------------------------------------------
2285 DO i=1,ms_nc !Over species a
2286 
2287  im=jms_nc(i)
2288 
2289  !nu_D and nu_T
2290  DO j=1,ms_nc !Over species b
2291 
2292  jm=jms_nc(j)
2293  c1=vti_nc(jm)/vti_nc(im)
2294  c2=x/c1
2295  phi=nclass_erf(c2)
2296  g=(phi-c2*(2.0_rspec/sqrt(z_pi))*exp(-c2**2)) &
2297  & /(2.0_rspec*c2**2)
2298  ynud_s(i)=ynud_s(i)+(3.0_rspec*sqrt(z_pi)/4.0_rspec) &
2299  & *(phi-g)/x**3/tauss_nc(i,j)
2300  ynut_s(i)=ynut_s(i)+((3.0_rspec*sqrt(z_pi)/4.0_rspec) &
2301  & *((phi-3.0_rspec*g)/x**3 &
2302  & +4.0_rspec*(temp_i(im)/temp_i(jm) &
2303  & +1.0_rspec/c1**2)*g/x))/tauss_nc(i,j)
2304 
2305  ENDDO !Over species b
2306 
2307  !nu_T*I_m
2308  DO m=1,mf_nc !Over poloidal moments
2309 
2310  IF(abs(p_ngrth) > 0.0_rspec) THEN
2311 
2312  c1=x*vti_nc(im)*m*p_ngrth
2313  c2=(ynut_s(i)/c1)**2
2314 
2315  IF(c2 > 9.0_rspec) THEN
2316 
2317  !Use asymptotic limit for efficiency
2318  ynutis(m,i)=0.4_rspec
2319 
2320  ELSE
2321 
2322  !Use full calculation
2323  c3=ynut_s(i)/c1*atan(c1/ynut_s(i))
2324  ynutis(m,i)=0.5_rspec*c3+c2*(3.0_rspec*(c3-0.5_rspec) &
2325  & +c2*4.5_rspec*(c3-1.0_rspec))
2326 
2327  ENDIF
2328 
2329  ELSE
2330 
2331  ynutis(m,i)=0.4_rspec
2332 
2333  ENDIF
2334 
2335  ENDDO !Over poloidal moments
2336 
2337 ENDDO !Over species a
2338 
2339 END SUBROUTINE nclass_nu
2340 
2341 SUBROUTINE nclass_tau(amu_i,temp_i,den_iz)
2342  !-------------------------------------------------------------------------------
2343  !NCLASS_TAU calculates the collision times for 90 degree scattering and the
2344  ! effective collision rates
2345  !References:
2346  ! S.P.Hirshman, D.J.Sigmar, Nucl Fusion 21 (1981) 1079
2347  ! W.A.Houlberg, K.C.Shaing, S.P.Hirshman, M.C.Zarnstorff, Phys Plasmas 4 (1997)
2348  ! 3230
2349  ! W.A.Houlberg 1/2002
2350  !Input:
2351  ! amu_i(i) -atomic mass number of i [-]
2352  ! temp_i(i) -temperature of i [keV]
2353  ! den_iz(i,z) -density of i,z [/m**3]
2354  !-------------------------------------------------------------------------------
2355 implicit none
2356  !Declaration of input variables
2357 REAL(KIND=rspec), INTENT(IN) :: &
2358  & amu_i(:),temp_i(:), &
2359  & den_iz(:,:)
2360 
2361 !Declaration of local variables
2362 INTEGER :: &
2363  & i,im,iz,iza,j,jm,jz,jza
2364 
2365 REAL(KIND=rspec) :: &
2366  & c1,c2,cln
2367 
2368 !-------------------------------------------------------------------------------
2369 !Initialization
2370 !-------------------------------------------------------------------------------
2371 !Arrays
2372 amntii_nc(:,:)=0.0_rspec
2373 tauss_nc(:,:)=0.0_rspec
2374 
2375 !Use electron Coulomb logarithm for all species
2376 cln=37.8_rspec-log(sqrt(den_iz(imel_nc,1))/temp_i(imel_nc))
2377 
2378 !-------------------------------------------------------------------------------
2379 !Collision times and mass density weighted collision rates
2380 !-------------------------------------------------------------------------------
2381 c1=4.0_rspec/3.0_rspec/sqrt(z_pi)*4.0_rspec*z_pi*cln*(z_coulomb &
2382  & /(4.0_rspec*z_pi*z_epsilon0))**2*(z_coulomb/z_protonmass)**2
2383 ! write(*,*) 'cln=',cln,imel_nc,'den=',den_iz(imel_nc,1)
2384 DO i=1,ms_nc
2385 
2386  im=jms_nc(i)
2387  iz=jzs_nc(i)
2388  iza=iabs(iz)
2389  c2=(vti_nc(im)**3)*amu_i(im)**2/c1
2390  if (im .eq. 12) then
2391  ! write(*,*) 'im=',im,jm,iza,iz
2392  ! write(*,*) 'c1=',c1,amu_i(im)
2393  endif
2394 
2395  DO j=1,ms_nc
2396 
2397  jm=jms_nc(j)
2398  jz=jzs_nc(j)
2399  jza=abs(jz)
2400  tauss_nc(i,j)=c2/iz**2/(den_iz(jm,jza)*jz**2)
2401  amntii_nc(im,jm)=amntii_nc(im,jm)+amu_i(im)*z_protonmass &
2402  & *den_iz(im,iza)/tauss_nc(i,j)
2403 
2404  END DO
2405 
2406 END DO
2407 
2408 END SUBROUTINE nclass_tau
2409 
2410 SUBROUTINE nclass_backsub(a,n,indx,b,ipr)
2411  !-------------------------------------------------------------------------------
2412  !NCLASS_BACKSUB solves the matrix equation a x = b, where a is in LU form as
2413  ! generated by a prior call to NCLASS_DECOMP
2414  !References:
2415  ! Flannery, Teukolsky, Vetterling, Numerical Recipes
2416  ! W.A.Houlberg 1/2002
2417  !Input:
2418  ! a(n,n) -coefficient matrix in lu decomposed form [-]
2419  ! n -number of equations to be solved [-]
2420  ! indx(n) -row permutations due to partial pivoting [-]
2421  !Input/output:
2422  ! b -(input) right hand side of equation [-]
2423  ! -(output) solution vector [-]
2424  !Comments:
2425  ! The complete procedure is, given the equation a*x = b,
2426  ! CALL NCLASS_DECOMP( )
2427  ! CALL NCLASS_BACKSUB( )
2428  ! and b is now the solution vector
2429  ! To solve with a different right hand side, just reload b as desired
2430  ! and use the same LU decomposition
2431  !-------------------------------------------------------------------------------
2432 implicit none
2433  !Declaration of input variables
2434 INTEGER, INTENT(IN) :: &
2435  & n, &
2436  & indx(:)
2437 
2438 REAL(KIND=rspec), INTENT(IN) :: &
2439  & a(:,:)
2440 
2441 !Declaration of input/output variables
2442 REAL(KIND=rspec), INTENT(INOUT) :: &
2443  & b(:)
2444 
2445 !Declaration of local variables
2446 INTEGER :: &
2447  & i,ii,j,k,ipr
2448 
2449 REAL(KIND=rspec) :: &
2450  & sumhoul
2451 
2452 !-------------------------------------------------------------------------------
2453 !Find the index of the first nonzero element of b
2454 !-------------------------------------------------------------------------------
2455 ii = 0
2456 sumhoul = 0.0_rspec
2457 ! if (n .eq. 18) then
2458 ! write(*,*) 'n=',n,indx
2459 ! endif
2460 DO i=1,n
2461 
2462  k = indx(i)
2463  sumhoul = b(k)
2464  b(k) = b(i)
2465  ! if (n .eq. 3 .and. ipr .eq. 1) then
2466  ! write(*,*) 'dans backsub, k=',k,' b(k)=',b(k),' i=',i
2467  ! write(*,*) 'ii=',ii,' sumhoul=',sumhoul
2468  ! endif
2469  IF(ii /= 0) THEN
2470 
2471  if (i .gt. 1) then
2472  DO j=ii,i-1
2473  sumhoul=sumhoul-a(i,j)*b(j)
2474  END DO
2475  endif
2476 
2477  ! ELSEIF (sumhoul /= 0.0_rspec) THEN
2478  ELSEIF (sumhoul /= 0.0) THEN
2479 
2480  ii=i
2481 
2482  ENDIF
2483 
2484  b(i)=sumhoul
2485  ! write(*,*) 'i=',i,'b(i)=',b(i)
2486  ! write(*,*) 'a(i,i)=',a(i,i)
2487  ! endif
2488 
2489 END DO
2490 ! if (iecrit .eq. 1 .and. iindice .eq. 2) then
2491 ! write(*,*) 'rspec=',rspec,'sumhoul=',sumhoul
2492 ! write(*,*) 'iecrit=',iecrit
2493 ! stop
2494 ! endif
2495 
2496 !-------------------------------------------------------------------------------
2497 !Back substitution
2498 !-------------------------------------------------------------------------------
2499 DO i=n,1,-1
2500 
2501  sumhoul=b(i)
2502 
2503  IF(i < n) THEN
2504 
2505  DO j=i+1,n
2506 
2507  sumhoul=sumhoul-a(i,j)*b(j)
2508 
2509  END DO
2510 
2511  ENDIF
2512 
2513  b(i)=sumhoul/a(i,i)
2514  ! if (n .eq. 3 .and. ipr .eq. 1) then
2515  ! write(*,*) 'i=',i,'b(i)=',b(i)
2516  ! write(*,*) 'a(i,i)=',a(i,i)
2517  ! endif
2518 
2519 END DO
2520 ! write(*,*) 'rspec=',rspec,'sumhoul=',sumhoul
2521 ! stop
2522 END SUBROUTINE nclass_backsub
2523 
2524 SUBROUTINE nclass_decomp(a,n,indx,d,iflag,ipr)
2525  !-------------------------------------------------------------------------------
2526  !NCLASS_DECOMP performs an LU decomposition of the matrix a and is called
2527  ! prior to NCLASS_BACKSUB to solve linear equations or to invert a matrix
2528  !References:
2529  ! Flannery, Teukolsky, Vetterling, Numerical Recipes
2530  ! W.A.Houlberg 1/2002
2531  !Input:
2532  ! n -number of equations to be solved [-]
2533  !Input/output:
2534  ! a(n,n) -coefficient matrix, overwritten on return [-]
2535  !Output:
2536  ! indx(n) -row permutations due to partial pivoting [-]
2537  ! d -flag for number of row exchanges [-]
2538  ! =1.0 even number
2539  ! =-1.0 odd number
2540  ! iflag -error and warning flag [-]
2541  ! =-1 warning
2542  ! =0 no warnings or errors
2543  ! =1 error
2544  !-------------------------------------------------------------------------------
2545 implicit none
2546  !Declaration of input variables
2547  INTEGER, INTENT(IN) :: &
2548  & n
2549 
2550  !Declaration of input/output variables
2551  REAL(KIND=rspec), INTENT(INOUT) :: &
2552  & a(:,:)
2553 
2554  !Declaration of output variables
2555  INTEGER, INTENT(OUT) :: &
2556  & iflag, &
2557  & indx(:)
2558 
2559  REAL(KIND=rspec) :: alp(n,n)
2560  REAL(KIND=rspec), INTENT(OUT) :: &
2561  & d
2562 
2563  !Declaration of local variables
2564  INTEGER :: &
2565  & i,imax,j,k,ik,ipr
2566 
2567  REAL(KIND=rspec) :: &
2568  & aamax,dum,sumhoul, &
2569  & vv(n),varia,sumhoult,precis,precis2
2570 
2571 
2572  !-------------------------------------------------------------------------------
2573  !Initialization
2574  !-------------------------------------------------------------------------------
2575  d = 1.0_rspec
2576  sumhoul = 0.0_rspec
2577  precis = 1e-19
2578  precis2 = 1e-19
2579 
2580  !-------------------------------------------------------------------------------
2581  !Loop over rows to get the implicit scaling information
2582  !-------------------------------------------------------------------------------
2583  DO i=1,n
2584 
2585  aamax=0.0_rspec
2586  DO j=1,n
2587  IF(abs(a(i,j)) > aamax) aamax=abs(a(i,j))
2588 
2589  END DO
2590 
2591  IF(aamax == 0.0_rspec) THEN
2592 
2593  iflag=1
2594  write(*,*) 'NCLASS_DECOMP/ERROR:singular matrix(1)',aamax
2595  print*,'erreur nclas_decomp'
2596  goto 9999
2597 
2598  ENDIF
2599 
2600  vv(i)=1.0_rspec/aamax
2601 
2602  ENDDO
2603  ! if (n .eq. 18) stop
2604  !-------------------------------------------------------------------------------
2605  !Use Crout's method for decomposition
2606  !-------------------------------------------------------------------------------
2607  !Loop over columns
2608 
2609  DO j=1,n
2610 
2611  if (j.eq. 1) goto 20
2612  DO i=1,j-1
2613  sumhoul=a(i,j)
2614  ik = i-1
2615  if (ik .ge. 1) then
2616 
2617  DO k=1,ik
2618 
2619  sumhoul=sumhoul-a(i,k)*a(k,j)
2620 
2621  ENDDO
2622 
2623  endif
2624 
2625  a(i,j)=sumhoul
2626 
2627  ENDDO
2628 20 continue
2629 
2630  !Search for largest pivot element using dum as a figure of merit
2631  aamax = 0.0_rspec
2632  dum = 0.0_rspec
2633  DO i=j,n
2634 
2635  sumhoul=a(i,j)
2636  sumhoult = sumhoul
2637  varia = 0.0_rspec
2638  if (j .gt. 1) then
2639  DO k=1,j-1
2640  sumhoul=sumhoul-a(i,k)*a(k,j)
2641  varia = varia + a(i,k)*a(k,j)
2642 
2643  ENDDO
2644  if (abs(sumhoul) .lt. precis) then
2645  sumhoul=sign(precis,sumhoul)
2646  ! sumhoul=0.0_rspec
2647  endif
2648  if (abs(sumhoul) .lt. precis) then
2649  ! write(*,*) 'sumhoul=',sumhoul
2650  ! sumhoul=sign(precis,sumhoul)
2651  ! sumhoul=precis
2652  ! write(*,*) 'sumhoul2=',sumhoul
2653 
2654  endif
2655  ! if (abs(sumhoul) .lt. precis) then
2656  ! write(*,*) 'sumhoul=',sumhoul,i
2657  ! sumhoul=0.0
2658  ! if (sumhoul .ge. 0) then
2659  ! sumhoul=precis
2660  ! else
2661  ! sumhoul=-precis
2662  ! endif
2663 
2664  ! endif
2665  endif
2666  a(i,j)=sumhoul
2667  dum=vv(i)*abs(sumhoul)
2668  ! if (n .eq. 18 .and. j .eq. 1) then
2669  !
2670  ! write(*,*) 'dum=',dum,' sumhoul=',sumhoul
2671  ! write(*,*) 'i=',i,'j=',j,' vv(i)=',vv(i)
2672  ! write(*,*) 'k=',k,' aamax=',aamax,' n=',n
2673  ! write(*,*) 'a(i,1)=',a(i,1),' imax=',imax
2674  !
2675  ! endif
2676  ! if (n .eq. 18) then
2677  ! write(*,*) 'dum=',dum,' aamax=',aamax,' i=',i,imax,j
2678  ! endif
2679  IF((dum >= aamax) .or. (abs(dum-aamax) < precis2)) THEN
2680  if (abs(dum-aamax) < precis2) then
2681  ! write(*,*) 'dum=',(dum-aamax)/precis2,i
2682  ! stop
2683  endif
2684  imax=i
2685  aamax=dum
2686  if (i .eq. 18) then
2687  ! stop
2688 
2689  endif
2690  ENDIF
2691 
2692  ENDDO
2693 
2694  ! write(*,*) '---- imax=',imax,' aamax=',aamax
2695 
2696  IF(j /= imax) THEN
2697 
2698  !Interchange rows
2699  DO k=1,n
2700 
2701  dum=a(imax,k)
2702  a(imax,k)=a(j,k)
2703  a(j,k)=dum
2704 
2705  ENDDO
2706 
2707  d=-d
2708  vv(imax)=vv(j)
2709 
2710  ENDIF
2711 
2712  indx(j)=imax
2713  ! if (n .eq. 18 .and. ipr .eq. 1) then
2714  ! write(*,*) 'j=',j,' imax=',imax
2715  ! write(*,*) 'dum=',dum,' indx(j)=',indx(j)
2716  ! stop
2717  ! endif
2718  IF(a(j,j) == 0.0_rspec) THEN
2719 
2720  iflag=1
2721  write(*,*) 'NCLASS_DECOMP/ERROR:singular matrix(2)'
2722  goto 9999
2723 
2724  ENDIF
2725 
2726  IF(j /= n) THEN
2727 
2728  !Divide by pivot element
2729  dum=1.0_rspec/a(j,j)
2730 
2731  DO i=j+1,n
2732 
2733  a(i,j)=a(i,j)*dum
2734 
2735  ENDDO
2736 
2737  ENDIF
2738 
2739  ENDDO
2740 
2741  ! write(*,*) 'a(2,2)=',a(2,2)
2742  ! stop
2743 
2744 9999 CONTINUE
2745 
2746 END SUBROUTINE nclass_decomp
2747 
2748 FUNCTION nclass_erf(x)
2749  !-------------------------------------------------------------------------------
2750  !NCLASS_ERF evaluates the error function erf(x)
2751  !References:
2752  ! M.Abramowitz, L.A.Stegun, Handbook of Math. Functions, p. 299
2753  ! W.A.Houlberg 1/2002
2754  !Input:
2755  ! x -argument of werror function [-]
2756  !Output:
2757  ! NCLASS_ERF -value of error function [-]
2758  !Comments:
2759  ! The error function can consume a lot of time deep in the multiple species
2760  ! loop so a very efficient calculation is called for
2761  ! This is much more critical than accuracy as suggested by T.Amano
2762  !-------------------------------------------------------------------------------
2763 implicit none
2764 REAL(KIND=rspec), INTENT(IN) :: &
2765  x !argument of error function [-]
2766 
2767 !Declaration of output variables
2768 REAL(KIND=rspec) :: &
2769  nclass_erf !value of error function [-]
2770 
2771 !-------------------------------------------------------------------------------
2772 !Declaration of local variables
2773 REAL(KIND=rspec) :: &
2774  t
2775 
2776 REAL(KIND=rspec), PARAMETER :: &
2777  c = 0.3275911_rspec, &
2778  a1 = 0.254829592_rspec, &
2779  a2 =-0.284496736_rspec, &
2780  a3 = 1.421413741_rspec, &
2781  a4 =-1.453152027_rspec, &
2782  a5 = 1.061405429_rspec
2783 
2784 !-------------------------------------------------------------------------------
2785 !Apply fit
2786 !-------------------------------------------------------------------------------
2787 t = one/(one+c*abs(x))
2788 nclass_erf = sign(one-(a1+(a2+(a3+(a4+a5*t)*t)*t)*t)*t*exp(-x**2),x)
2789 
2790 END FUNCTION nclass_erf
2791 
2792 SUBROUTINE normalisation(RMAT,B,n,VECM)
2793 implicit none
2794 !Declaration of input/output variables
2795 REAL(KIND=rspec), INTENT(INOUT) :: rmat(:,:),b(:)
2796 !Declaration of output variables
2797 REAL(KIND=rspec), INTENT(OUT) :: vecm(:)
2798 !Declaration of input variables
2799 INTEGER, INTENT(in) :: n
2800 
2801 INTEGER :: k,j
2802 
2803 do k=1,n
2804  vecm(k) = -1e10
2805 
2806  do j=1,n
2807  if (vecm(k) .le. rmat(k,j)) vecm(k) = rmat(k,j)
2808  enddo
2809  rmat(k,:) = rmat(k,:)/vecm(k)
2810  b(k) = b(k)/vecm(k)
2811 enddo
2812 
2813 END SUBROUTINE normalisation
2814 
2815 END MODULE nclass_itm
2816 
2817 SUBROUTINE neo(equilibrium, coreprof, neoclassic)
2818 use euitm_schemas
2819 use euitm_routines
2820 use nclass_itm
2821 
2822 implicit none
2823 
2824 type(type_neoclassic),pointer :: neoclassic(:)
2825 type(type_coreprof),pointer :: coreprof(:)
2826 type(type_equilibrium),pointer :: equilibrium(:)
2827 type(type_toroidfield),pointer :: toroidfield(:)
2828 integer :: nrho,nion,neq
2829 
2830 Interface
2831  SUBROUTINE neo_calc(equilibrium, coreprof,neoclassic,nrho,neq,nion)
2832 
2833  use euitm_schemas
2834  use euitm_routines
2835  use nclass_itm
2836  implicit none
2837  type(type_neoclassic),pointer :: neoclassic(:)
2838  type(type_coreprof),pointer :: coreprof(:)
2839  type(type_equilibrium),pointer :: equilibrium(:)
2840  type(type_toroidfield),pointer :: toroidfield(:)
2841  integer :: nrho,nion
2842  END SUBROUTINE neo_calc
2843 end Interface
2844 
2845  nrho=size(coreprof(1)%te%value)
2846  neq=size(equilibrium(1)%profiles_1d%rho_tor)
2847 ! if(nrho.ne.neq) then
2848 ! write(*,*) 'NEO: NRHO must be the same as NEQ ', NRHO, NEQ
2849 ! stop 'ERROR: NRHO != NEQ'
2850 ! endif
2851  nion=size(coreprof(1)%ti%value,2)
2852  call neo_calc(equilibrium, coreprof, neoclassic,nrho,neq,nion)
2853 end SUBROUTINE neo
2854 
2855 
2856 SUBROUTINE neo_calc(equilibrium,coreprof,neoclassic,nrho,neq,nion)
2857 !ITM link
2858 !CPO : neoclassic, coreprof, equilibrium
2859 
2860 use euitm_schemas
2861 use euitm_routines
2862 use cos_util
2863 use nclass_itm
2864 !use itm_toolbox
2865 
2866 implicit none
2867 
2868 type(type_neoclassic),pointer :: neoclassic(:)
2869 type(type_coreprof),pointer :: coreprof(:)
2870 type(type_equilibrium),pointer :: equilibrium(:)
2871 type(type_toroidfield),pointer :: toroidfield(:)
2872 real*8 :: tnp1
2873 
2874 integer :: num,kboot,kpotato,k_order,kpfirsch
2875 integer :: icharge, kbanana, kmmaj, ksmaj
2876 integer :: m_i, m_z, nspec, ipr
2877 integer :: nrho,neq,nion
2878 integer, PARAMETER :: mxmi=100, nmaxspec=10,mxmz=100,mxms=100
2879 REAL*8 p_b2, p_bm2, &
2880  & p_eb, p_fhat, &
2881  & p_fm(3), p_ft, &
2882  & p_grbm2, p_grphi,&
2883  & p_gr2phi, p_ngrth
2884 REAL*8 cden,tabtr(1000,200)
2885 real*8 pcharge(nmaxspec),pmasse(nmaxspec)
2886 REAL*8 amu_i(mxmi) , grti(200,mxmi) , &
2887  & tempi(200,mxmi), ztemp(3) , &
2888  & den(200,mxmi) , grden(200,mxmi) , &
2889  & temp_i(mxmi) , grt_i(mxmi)
2890 REAL*8 rhomax, q0, r0, b0, rkappa0, c_potb, c_pot1, xqs
2891 REAL*8 x(nrho), rho(nrho)
2892 REAL*8 psi(nrho) , psidrho(nrho) , r_in(nrho) , r_out(nrho)
2893 REAL*8 xr0_pr(nrho) , rminx_pr(nrho) , p_ft_pr(nrho) , xqs_pr(nrho)
2894 REAL*8 bt0_pr(nrho) , xvpr_pr(nrho) , gph_pr(nrho) , grho2_pr(nrho) , dpsidrho(nrho)
2895 REAL*8 p_b2_pr(nrho) , p_bm2_pr(nrho) , grho_pr(nrho) , p_eb_pr(nrho)
2896 REAL*8 ergrho(nrho) , inter(nrho) , gr2phi(nrho) , xfs_pr(nrho)
2897 REAL*8 force1(nrho) , force2(nrho) , force3(nrho) , p_grbm2_pr(nrho)
2898 REAL*8 xfs , rap(nrho) , gm1(nrho) , gm9(nrho)
2899 REAL*8 den_iz(mxmi,mxmz), fex_iz(3,mxmi,mxmz), &
2900  & grp_iz(mxmi,mxmz)
2901 REAL*8 denz2(mxmi) , vti(mxmi)
2902 REAL*8 pgrp_iz(mxmi,mxmz) , xi(mxmi,mxmz)
2903 REAL*8 amntii(mxmi,mxmi)
2904 REAL*8 tau_ss(mxms,mxms) , ephi(2)
2905 REAL*8 calm_i(3,3,mxmi)
2906 REAL*8 caln_ii(3,3,mxmi,mxmi),capm_ii(3,3,mxmi,mxmi),&
2907  & capn_ii(3,3,mxmi,mxmi)
2908 REAL*8 bsjbp_s(mxms), bsjbt_s(mxms),&
2909  & DN_S(mxms), GFL_S(5,mxms),&
2910  & QFL_S(5,mxms), SQZ_S(mxms),&
2911  & UPAR_S(3,3,mxms), UTHETA_S(3,3,mxms),&
2912  & VNNT_S(mxms), VNEB_S(mxms),&
2913  & VNEX_S(mxms), VQEB_S(mxms),&
2914  & VQNT_S(mxms), VQEX_S(mxms),&
2915  & CHI_S(mxms), &
2916  & qebs(mxms), XI_S(mxms),&
2917  & YMU_S(3,3,mxms)
2918 REAL*8 chip_ss(mxms,mxms), chit_ss(mxms,mxms),&
2919  & DP_SS(mxms,mxms), DT_SS(mxms,mxms)
2920 REAL*8 ni(200),xr0,xeps
2921 integer :: jm_s(mxms),jz_s(mxms)
2922 logical :: l_banana,l_pfirsch,l_potato
2923 real*8 :: c_potl,p_etap, p_jbbs, p_jbex, p_jboh,c_den,rminx
2924 real*8 :: bt0,dencut,dent,xbeta,xdelp,xvpr,ab,xgph,eps2,b,a,c,xm
2925 real*8 :: dentot, chitp, vconv,etatemp(4)
2926 integer :: m_s,m
2927 integer :: itp1,i,ia,ima,ki,kspec,iza,iflag,j,indsave
2928 
2929 
2930 itp1=1
2931 allocate(neoclassic(1))
2932 neoclassic(itp1)%time=coreprof(itp1)%time
2933 allocate(neoclassic(itp1)%rho_tor(nrho))
2934 allocate(neoclassic(itp1)%rho_tor_norm(nrho))
2935 allocate(neoclassic(itp1)%jboot(nrho))
2936 allocate(neoclassic(itp1)%sigma(nrho))
2937 allocate(neoclassic(itp1)%te_neo%flux(nrho))
2938 allocate(neoclassic(itp1)%ti_neo%flux(nrho,nion))
2939 allocate(neoclassic(itp1)%ne_neo%flux(nrho))
2940 allocate(neoclassic(itp1)%ni_neo%flux(nrho,nion))
2941 allocate(neoclassic(itp1)%te_neo%diff_eff(nrho))
2942 allocate(neoclassic(itp1)%ti_neo%diff_eff(nrho,nion))
2943 allocate(neoclassic(itp1)%ne_neo%diff_eff(nrho))
2944 allocate(neoclassic(itp1)%ni_neo%diff_eff(nrho,nion))
2945 allocate(neoclassic(itp1)%te_neo%vconv_eff(nrho))
2946 allocate(neoclassic(itp1)%ti_neo%vconv_eff(nrho,nion))
2947 allocate(neoclassic(itp1)%ne_neo%vconv_eff(nrho))
2948 allocate(neoclassic(itp1)%ni_neo%vconv_eff(nrho,nion))
2949 
2950 !____________
2951 ! data access
2952 !____________
2953 nspec = nion+1 ! add electron species
2954 pmasse(1) = 1 ! electron
2955 pmasse(2:nspec) = coreprof(itp1)%composition%amn ! species mass
2956 pcharge(1) = -1 ! electron
2957 pcharge(2:nspec) = coreprof(itp1)%composition%zn
2958 !----------------------------------------------------------------------
2959 ! Model options
2960 ! kboot : internal option for metrics
2961 ! = 1 -> Hirshman, Sigmar model: fm(k)=Lck^*,
2962 ! xngrth=<(n.gr(B))^2>/<B^2>.
2963 ! = 2 -> Kessel model: fm(1)=Rq/eps**3/2.
2964 ! xngrth not used.
2965 ! = -> Shaing, et al model: fm(k)=Fk,
2966 ! xngrth=<n.gr(theta)>
2967 !----------------------------------------------------------------------
2968 !kboot = neoclassic(itp1)%codeparam%parameters(1)
2969 kboot = 0
2970 ! by default : kboot = 0nMaxSpec
2971 !----------------------------------------------------
2972 ! kpotato : option to include potato orbits (-)
2973 ! = 0 -> off
2974 ! = -> else on
2975 !-----------------------------------------------------
2976 !kpotato = neoclassic(itp1)%codeparam%parameters(2)
2977 kpotato = 1
2978 ! by default : kpotato = 1
2979 !-------------------------------------------------------
2980 ! k_order : order of v moments to be solved (-)
2981 ! = 2 -> u and q
2982 ! = 3 -> u, q, and u2
2983 ! = -> else error
2984 !-------------------------------------------------------
2985 !k_order = neoclassic(itp1)%codeparam%parameters(3)
2986 k_order = 3
2987 ! by default : k_order = 3
2988 !------------------------------------------------
2989 ! kbanana : option to include banana viscosity (-)
2990 ! = 0 -> off
2991 ! = -> else on
2992 !------------------------------------------------
2993 !kbanana = neoclassic(itp1)%codeparam%parameters(4)
2994 kbanana = 1
2995 ! by default : kbanana = 1
2996 !------------------------------------------------------------------
2997 ! kpfirsch : option to include Pfirsch-Schluter viscosity (-)
2998 ! = 0 -> off
2999 ! = -> else on
3000 !------------------------------------------------------------------
3001 !kpfirsch = neoclassic(itp1)%codeparam%parameters(5)
3002 kpfirsch = 1
3003 m_i = nspec
3004 m_z = 0
3005 DO i=1,nspec
3006  icharge = nint(pcharge(i))
3007  IF (icharge.gt.m_z) m_z = icharge
3008 ENDDO
3009 amu_i(1) = z_electronmass/z_protonmass
3010 DO i=2,nspec
3011  amu_i(i) = pmasse(i)
3012 ENDDO
3013 !------------------------------------------------------------------------
3014 ! c_den : density cutoff below which species is ignored (/m**3)
3015 !------------------------------------------------------------------------
3016 c_den = 1.0e10
3017 !-----------------------------------------------------------------------------
3018 ! p_grphi : radial electric field Phi' (V/rho)
3019 ! p_gr2phi : radial electric field gradient Psi'(Phi'/Psi')' (V/rho**2)
3020 !-----------------------------------------------------------------------------
3021 !!!x(1:nrho) = coreprof(itp1)%rho_tor_norm
3022 !!!rhomax = maxval(coreprof(itp1)%rho_tor(:))
3023 !!!rho = x*rhomax
3024 ! dpc
3025 rho(1:nrho) = coreprof(1)%rho_tor !!! was equilibrium(1)%profiles_1d%rho_tor
3026 rhomax = rho(nrho)
3027 x(1:nrho) = rho/rhomax
3028 ! cpd
3029 !----------------------
3030 ! temprature in keV
3031 !----------------------
3032 tempi(1:nrho,1) = coreprof(itp1)%te%value(1:nrho)/1000.0
3033 do kspec=2,nspec
3034  tempi(1:nrho,kspec) = coreprof(itp1)%ti%value(:,kspec-1)/1000.0
3035 enddo
3036 den(1:nrho,1) = coreprof(itp1)%ne%value
3037 ni = 0
3038 do kspec=2,nspec
3039  den(1:nrho,kspec) = coreprof(itp1)%ni%value(:,kspec-1)
3040  ni(:) = ni(:) + den(:,kspec)
3041 enddo
3042 do kspec = 1,nspec
3043  !print*,"calcul de dTi/dx, Ti= ",tempi(101,kspec)
3044  !print*,"kspec=",kspec
3045  call cos_rpdederive(grti(:,kspec),nrho,x,tempi(:,kspec),0,2,2,1)
3046 enddo
3047 ! den : coreprof(itp1).ne.value
3048 do kspec = 1,nspec
3049  !print*,"calcul de dni/dx, ni= ",den(1:4,kspec)
3050  !print*,"kspec=",kspec
3051  call cos_rpdederive(grden(:,kspec),nrho,x,den(:,kspec),0,2,2,1)
3052 enddo
3053 ! equilibrium data
3054 !!!psi(1:nrho) = coreprof(itp1)%psi%value
3055 ! dpc
3056 !!!psi(1:nrho) = equilibrium(1)%profiles_1d%psi
3057 !psi(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%psi, rho(1:nrho))
3058 CALL l3interp(equilibrium(1)%profiles_1d%psi, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3059  psi, rho, nrho)
3060 ! cpd
3061 call cos_rpdederive(psidrho,nrho,rho,psi,0,2,2,1)
3062 !!!r_in(1:nrho) = equilibrium(1)%profiles_1d%r_inboard
3063 !r_in(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%r_inboard, rho(1:nrho))
3064 CALL l3interp(equilibrium(1)%profiles_1d%r_inboard, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3065  r_in, rho, nrho)
3066 !!!r_out(1:nrho) = equilibrium(1)%profiles_1d%r_outboard
3067 !r_out(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%r_outboard, rho(1:nrho))
3068 CALL l3interp(equilibrium(1)%profiles_1d%r_outboard, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3069  r_out, rho, nrho)
3070 xr0_pr = (r_out+r_in)/2
3071 xr0_pr = r_out !!! ??? DPC
3072 
3073 !rminx_pr(1:nrho) = equilibrium%profiles_1d%rho_rttorfl
3074 rminx_pr(1:nrho) = rho(1:nrho)
3075 !!!p_ft_pr(1:nrho) = equilibrium(1)%profiles_1d%ftrap
3076 !p_ft_pr(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%ftrap, rho(1:nrho))
3077 CALL l3interp(equilibrium(1)%profiles_1d%ftrap, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3078  p_ft_pr, rho, nrho)
3079 !xqs_pr(1:nrho) = coreprof%profiles1d%q%value
3080 !!!xqs_pr(1:nrho) = equilibrium(1)%profiles_1d%q
3081 !xqs_pr(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%q, rho(1:nrho))
3082 CALL l3interp(equilibrium(1)%profiles_1d%q, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3083  xqs_pr, rho, nrho)
3084 q0 = xqs_pr(1)
3085 !!!DPC-EQ-4.08b-problem
3086 !!!xvpr_pr(1:nrho) = equilibrium(1)%profiles_1d%vprime
3087 !dpsidrho(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%dpsidrho_tor, rho(1:nrho))
3088 CALL l3interp(equilibrium(1)%profiles_1d%dpsidrho_tor, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3089  dpsidrho, rho, nrho)
3090 !xvpr_pr(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%vprime, rho(1:nrho))
3091 CALL l3interp(equilibrium(1)%profiles_1d%vprime, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3092  xvpr_pr, rho, nrho)
3093 xvpr_pr = xvpr_pr * dpsidrho * rho(nrho)
3094 
3095 
3096 !!!gm1(1:nrho) = equilibrium(1)%profiles_1d%gm1
3097 !gm1(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%gm1, rho(1:nrho))
3098 CALL l3interp(equilibrium(1)%profiles_1d%gm1, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3099  gm1, rho, nrho)
3100 !!!gm9(1:nrho) = equilibrium(1)%profiles_1d%gm9
3101 !gm9(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%gm9, rho(1:nrho))
3102 CALL l3interp(equilibrium(1)%profiles_1d%gm9, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3103  gm9, rho, nrho)
3104 !!!bt0_pr(1:nrho) = equilibrium(1)%profiles_1d%F_dia * gm1 / gm9
3105 !bt0_pr(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%F_dia, rho(1:nrho)) * gm1 / gm9
3106 CALL l3interp(equilibrium(1)%profiles_1d%F_dia * gm1 / gm9, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3107  bt0_pr, rho, nrho)
3108 r0 = equilibrium(1)%eqgeometry%geom_axis%r
3109 b0 = coreprof(itp1)%toroid_field%b0
3110 !!!gph_pr(1:nrho) = equilibrium(1)%profiles_1d%gm1 * rhomax * xvpr_pr(1:nrho)
3111 gph_pr(1:nrho) = gm1(1:nrho) * rhomax * xvpr_pr
3112 !!!grho2_pr (1:nrho) = equilibrium(1)%profiles_1d%gm3
3113 !grho2_pr(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%gm3, rho(1:nrho))
3114 CALL l3interp(equilibrium(1)%profiles_1d%gm3, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3115  grho2_pr, rho, nrho)
3116 !!!p_b2_pr(1:nrho) = equilibrium(1)%profiles_1d%gm5
3117 !p_b2_pr(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%gm5, rho(1:nrho))
3118 CALL l3interp(equilibrium(1)%profiles_1d%gm5, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3119  p_b2_pr, rho, nrho)
3120 !!!p_bm2_pr(1:nrho) = equilibrium(1)%profiles_1d%gm4
3121 !p_bm2_pr(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%gm4, rho(1:nrho))
3122 CALL l3interp(equilibrium(1)%profiles_1d%gm4, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3123  p_bm2_pr, rho, nrho)
3124 !!!p_grbm2_pr(1:nrho) = equilibrium(1)%profiles_1d%gm6
3125 !p_grbm2_pr(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%gm6, rho(1:nrho))
3126 CALL l3interp(equilibrium(1)%profiles_1d%gm6, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3127  p_grbm2_pr, rho, nrho)
3128 !!!grho_pr(1:nrho) = equilibrium(1)%profiles_1d%gm7
3129 !grho_pr(1:nrho) = interpolate(equilibrium(1)%profiles_1d%rho_tor, equilibrium(1)%profiles_1d%gm7, rho(1:nrho))
3130 CALL l3interp(equilibrium(1)%profiles_1d%gm7, equilibrium(1)%profiles_1d%rho_tor, size(equilibrium(1)%profiles_1d%rho_tor), &
3131  grho_pr, rho, nrho)
3132 rkappa0 = equilibrium(1)%eqgeometry%elongation
3133 write(*,*) 'rkappa0 = ', rkappa0
3134 p_eb_pr(1:nrho) = coreprof(itp1)%profiles1d%eparallel%value
3135 !p_eb_pr = p_eb_pr * bt0_pr
3136 p_eb_pr = p_eb_pr * b0
3137 c_potb = rkappa0*bt0/2/q0/q0
3138 c_potl = q0*xr0
3139 
3140 !attention a prendre de l'equilibre
3141 !ergrho(1:nrho) = neoclassic(itp1)%er / grho_pr
3142 ergrho(1:nrho) = 0
3143 
3144 inter(2:nrho) = -ergrho(2:nrho) / psidrho(2:nrho)
3145 call cos_zconversion(inter,nrho)
3146 !print*,"calcul de dinter/dx, inter= ",inter(1:4)
3147 call cos_rpdederive(gr2phi,nrho,rho,inter,0,2,2,1)
3148 gr2phi = psidrho * gr2phi
3149 ! not in the CPO
3150 ! force1 = 1er moment des forces externes pour l'ion principal (datak.source.totale.q ?)
3151 ! force2 = 2ieme moment des forces externes (?)
3152 ! force3 = 3ime moment des forces externes pour l'ion principal (datak.source.totale.wb ?)
3153 force1 = 0.
3154 force2 = 0.
3155 force3 = 0.
3156 !
3157 ! kmmaj : masse de la particule NBI injectee
3158 ! ksmaj : charge de la particule NBI injectee
3159 !
3160 kmmaj = 2
3161 ksmaj = 2
3162 
3163 !--------------------------------------------------------------------|
3164 ! xfs : courant poloidale externe a une surface de flux (A)|
3165 !--------------------------------------------------------------------|
3166 xfs_pr = 2.0*z_pi*xr0_pr*bt0_pr/z_mu0
3167 !-----------------------
3168 ! Loop over radial nodes
3169 !-----------------------
3170 !DO ipr=1,nrho
3171 DO ipr=2,nrho
3172  !
3173  ! external force set to zero
3174  !
3175  p_grphi = 0.0
3176  p_gr2phi = 0.0
3177  fex_iz(1,kmmaj,ksmaj) = force1(ipr)
3178  fex_iz(2,kmmaj,ksmaj) = force2(ipr)
3179  fex_iz(3,kmmaj,ksmaj) = force3(ipr)
3180  dencut = 1e10
3181  xr0 = xr0_pr(ipr)
3182  rminx = rminx_pr(ipr)
3183  xqs = xqs_pr(ipr)
3184  if (ipr .eq. 4) then
3185 
3186  ! xqs = 0.6772616311415705
3187 
3188  endif
3189  p_ft = p_ft_pr(ipr)
3190  bt0 = bt0_pr(ipr)
3191  xdelp = xvpr_pr(ipr)
3192  xvpr = xdelp
3193  ab = 1.0
3194  p_b2 = p_b2_pr(ipr)
3195  p_bm2 = p_bm2_pr(ipr)
3196  q0 = xqs_pr(1)
3197  xeps = rminx/xr0
3198  xgph = gph_pr(ipr)/4/z_pi/z_pi
3199  xfs = xfs_pr(ipr)
3200  DO ima=1,m_i
3201  DO ia=1,m_z
3202  den_iz(ima,ia) = dencut
3203  ENDDO
3204  ENDDO
3205  DO ima=2,m_i
3206  icharge = nint(pcharge(ima))
3207  den_iz(ima,icharge) = den(ipr,ima)
3208  ENDDO
3209  den_iz(1,1) = 0.0
3210  dent = 0.0
3211  DO ima=2,m_i
3212  denz2(ima) = 0.0
3213  DO ia=1,m_z
3214  den_iz(1,1) = den_iz(1,1)+float(ia)*den_iz(ima,ia)
3215  denz2(ima) = denz2(ima)+den_iz(ima,ia)&
3216  & * float(ia)**2
3217  dent = dent+den_iz(ima,ia)*tempi(ipr,ima)
3218  ENDDO
3219  ENDDO
3220  denz2(1) = den_iz(1,1)
3221  dent = dent+den_iz(1,1)*tempi(ipr,1)
3222  xbeta = dent*z_j7kv/(bt0**2/2.0/z_mu0)
3223  DO ima=1,m_i
3224  DO ia=1,m_z
3225  xi(ima,ia) = den_iz(ima,ia)*float(ia)**2&
3226  & / denz2(ima)
3227  IF(den_iz(ima,ia).eq.0) THEN
3228  grp_iz(ima,ia) = 0
3229  ELSE
3230  grp_iz(ima,ia) = grti(ipr,ima)*den_iz(ima,ia)&
3231  & +tempi(ipr,ima)*grden(ipr,ima)
3232  ENDIF
3233  ENDDO
3234  ENDDO
3235  if (xgph .eq. 0) then
3236  xgph = 0.01
3237  endif
3238  p_fhat = xqs/xgph
3239  p_grbm2 = p_grbm2_pr(ipr)
3240  p_eb = p_eb_pr(ipr)
3241  c_potb = rkappa0*bt0/2/q0/q0
3242  c_potl = q0*xr0
3243  IF(kboot.eq.1) THEN
3244  !-----------------
3245  ! Hirshman
3246  !-----------------
3247  p_ngrth = (xeps/(xqs*xr0))**2/2.0
3248  p_fm(1) = xqs*xr0
3249  p_fm(2) = 0.0
3250  p_fm(3) = 0.0
3251  ELSEIF(kboot.eq.2) then
3252  !---------------
3253  ! Kessel
3254  !---------------
3255  p_ngrth = 0.0
3256  p_fm(1) = xqs*xr0/xeps**1.5
3257  ELSE
3258  !---------------
3259  ! Shaing
3260  !---------------
3261  p_ngrth = 1.0/(xqs*xr0)
3262  DO m=1,3
3263  p_fm(m) = 0
3264  ENDDO
3265  IF (xeps.gt.0.0) THEN
3266  eps2 = xeps**2
3267  b = sqrt(1.0-eps2)
3268  a = (1.0-b)/xeps
3269  c = (b**3.0)*(xqs*xr0)**2
3270  DO m=1,3
3271  xm = float(m)
3272  p_fm(m) = xm*a**(2.0*xm)*(1.0+xm*b)/c
3273  ENDDO
3274  ENDIF
3275  ENDIF
3276  !-------------------------------------------
3277  ! Find significant charge states and mapping
3278  !-------------------------------------------
3279  m_s = 0
3280  DO ima=1,m_i
3281  DO iza=1,m_z
3282  IF(den_iz(ima,iza).gt.c_den) THEN
3283  m_s = m_s+1
3284  !---------------------------------------------------------------
3285  ! Set isotope number and charge state for this species
3286  !---------------------------------------------------------------
3287  jm_s(m_s) = ima
3288  IF(amu_i(ima).lt.0.5) THEN
3289  jz_s(m_s) = -iza
3290  ELSE
3291  jz_s(m_s) = iza
3292  ENDIF
3293  ENDIF
3294  ENDDO
3295  ENDDO
3296  !---------------------------
3297  ! Calculate thermal velocity
3298  !---------------------------
3299  DO ima=1,m_i
3300  vti(ima) = sqrt(2.0*z_j7kv*tempi(ipr,ima)&
3301  & / amu_i(ima)/z_protonmass)
3302  temp_i(ima) = tempi(ipr,ima)
3303  grt_i(ima) = grti(ipr,ima)
3304  ENDDO
3305  if (kbanana .ne. 0) then
3306  l_banana = .true.
3307  else
3308  l_banana = .false.
3309  endif
3310  if (kpotato .ne. 0) then
3311  l_potato = .true.
3312  else
3313  l_potato = .false.
3314  endif
3315  if (kpfirsch .ne. 0) then
3316  l_pfirsch = .true.
3317  else
3318  l_pfirsch = .false.
3319  endif
3320  if (ipr .eq. 4) then
3321 
3322 ! write(*,*) 'p_b2 =',p_b2,' p_bm2=',p_bm2,' p_eb=',p_eb
3323 ! write(*,*) 'p_fhat =',p_fhat,'p_fm=',p_fm
3324 ! write(*,*) 'xqs=',xqs,' xgph=',xgph
3325 ! write(*,*) 'p_ft =',p_ft,' p_grbm2 =', p_grbm2,' p_grphi=',p_grphi
3326 ! write(*,*) 'p_gr2phi=',p_gr2phi,' p_ngrth=',p_ngrth
3327 ! write(*,*) 'grt_i =',grt_i(1),'temp_i =', temp_i(1),' grp_iz=',grp_iz(1,1)
3328 
3329  endif
3330  CALL nclass( m_i , m_z , p_b2 , p_bm2 , p_eb, &
3331  & p_fhat , p_fm , p_ft , p_grbm2 , p_grphi, &
3332  & p_gr2phi , p_ngrth , amu_i , grt_i , temp_i, &
3333  & den_iz , fex_iz , grp_iz , ipr , iflag, &
3334  & l_banana , l_pfirsch , l_potato , k_order , c_den, &
3335  & c_potb , c_potl , p_etap , p_jbbs , p_jbex, &
3336  & p_jboh , m_s , jm_s , jz_s , bsjbp_s, &
3337  & bsjbt_s , gfl_s , dn_s , vnnt_s , vneb_s, &
3338  & vnex_s , dp_ss , dt_ss , upar_s , utheta_s, &
3339  & qfl_s , chi_s , vqnt_s , vqeb_s , vqex_s, &
3340  & chip_ss , chit_ss , calm_i , caln_ii , capm_ii, &
3341  & capn_ii , ymu_s , sqz_s , xi_s , tau_ss)
3342 
3343 
3344  neoclassic(itp1)%rho_tor = rho
3345  neoclassic(itp1)%rho_tor_norm = x
3346  ! mise en forme des sorties
3347  !
3348  !----------------------
3349  ! Variables de sortie
3350  ! <Jbs.B> (A*T/m**2)
3351  !----------------------
3352  ! tabtr(1,ipr) = P_JBBS/zj7kv
3353  tabtr(1,ipr) = p_jbbs
3354  neoclassic(itp1)%jboot(ipr) = p_jbbs / b0
3355  !write(*,*) P_JBBS / b0
3356  !------------------------------------------------------------------
3357  ! <Jbs.B> driven by unit p'/p of species s (A*T*rho/m**3)
3358  ! s : 1 -> electrons, 2 -> espece ionique 1, 3 -> espece ionique 2
3359  !------------------------------------------------------------------
3360  tabtr(2,ipr) = bsjbp_s(1)
3361  tabtr(3,ipr) = bsjbp_s(2)
3362  tabtr(4,ipr) = bsjbp_s(3)
3363  !------------------------------------------------------------------
3364  ! <Jbs.B> driven by unit T'/T of s (A*T*rho/m**3)
3365  ! s : 1 -> electrons, 2 -> espece ionique 1, 3 -> espece ionique 2
3366  !------------------------------------------------------------------
3367  tabtr(5,ipr) = bsjbt_s(1)
3368  tabtr(6,ipr) = bsjbt_s(2)
3369  tabtr(7,ipr) = bsjbt_s(3)
3370  !-----------------------------------------------
3371  ! <Jex.B> current response to fex_iz (A*T/m**2)
3372  !-----------------------------------------------
3373  tabtr(8,ipr) = p_jbex
3374  !-----------------------------------------
3375  ! Parallel electrical resistivity (Ohm*m)
3376  !-----------------------------------------
3377  tabtr(9,ipr) = p_etap
3378  neoclassic(itp1)%sigma(ipr) = 1.0 / p_etap
3379  !write(*,*) 'ipr=',ipr,' sigma=',1.0 / P_ETAP
3380  !------------------------------------------------------------
3381  ! Particle and heat fluxes
3382  ! For each species a and charge state i
3383  ! gfl : particle flux (rho/m**3/s)
3384  ! qfl : conduction heat flux (J*rho/m**3/s)
3385  ! qfl+5.2*T*gfl : total radial heat flux
3386  ! (1,a,i) : p' and T' driven banana-plateau flux
3387  ! (2,a,i) : Pfirsch-Schluter flux
3388  ! (3,a,i) : classical flux
3389  ! (4,a,i) : <E.B> driven flux
3390  ! (5,a,i) : external force driven flux
3391  !------------------------------------------------------------
3392  tabtr(10,ipr) = 0.0d0
3393  tabtr(11,ipr) = 0.0d0
3394  tabtr(12,ipr) = 0.0d0
3395  tabtr(13,ipr) = 0.0d0
3396  !-------------------------------
3397  ! Sum over flux components
3398  !-------------------------------
3399  DO j=1,5
3400  !------------------------------------------
3401  ! Electron conduction heat flow (w)
3402  !------------------------------------------
3403  tabtr(10,ipr) = tabtr(10,ipr)+qfl_s(j,1)
3404  tabtr(12,ipr) = tabtr(12,ipr)+gfl_s(j,1)
3405  DO ima=2,m_i
3406  !-----------------------------------------
3407  ! Ion conduction heat flow (w)
3408  !-----------------------------------------
3409  tabtr(11,ipr) = tabtr(11,ipr)+qfl_s(j,ima)
3410  tabtr(13,ipr) = tabtr(13,ipr)+gfl_s(j,ima)
3411  ENDDO
3412  ENDDO
3413  ! ITM ------------------------->
3414  neoclassic(itp1)%te_neo%flux(ipr) = tabtr(10,ipr) * rhomax / grho2_pr(ipr)
3415 ! neoclassic(itp1)%ti_neo%flux(ipr,1) = tabtr(11,ipr) * rhomax / grho2_pr(ipr)
3416  neoclassic(itp1)%ne_neo%flux(ipr) = tabtr(12,ipr) * rhomax / grho2_pr(ipr)
3417 ! neoclassic(itp1)%ni_neo%flux(ipr,1) = tabtr(13,ipr) * rhomax / grho2_pr(ipr)
3418  do ima=2,m_i
3419  neoclassic(itp1)%ti_neo%flux(ipr,ima-1) = 0.0d0
3420  neoclassic(itp1)%ni_neo%flux(ipr,ima-1) = 0.0d0
3421  do j=1,5
3422  neoclassic(itp1)%ti_neo%flux(ipr,ima-1)=neoclassic(itp1)%ti_neo%flux(ipr,ima-1) + qfl_s(j,ima) * rhomax / grho2_pr(ipr)
3423  neoclassic(itp1)%ni_neo%flux(ipr,ima-1)=neoclassic(itp1)%ni_neo%flux(ipr,ima-1) + gfl_s(j,ima) * rhomax / grho2_pr(ipr)
3424  enddo
3425  enddo
3426  ! ITM ------------------------->
3427  !--------------------------------------------------------------------
3428  !
3429  ! tabtr(14) : diffusion coefficient (diag comp) of s (rho**2/s)
3430  ! tabtr(15) : somme des differentes vitesses electroniques
3431  ! vns -> convection velocity (off diag comps-p', T') of s (rho/s)
3432  ! vebs -> <E.B> particle convection velocity of s (rho/s)
3433  ! qfl(5,1) -> external force driven flux
3434  ! tabtr(16) : heat electronic cond coefficient of s2 on p'/p of s1 (rho**2/s)
3435  ! + heat electronic cond coefficient of s2 on T'/T of s1 (rho**2/s)
3436  ! tabtr(17) : heat electronic convection velocity (rho/s)
3437  ! tabtr(18) : sum [heat ionic cond coefficient of s2 on p'/p of s1 (rho**2/s)
3438  ! + heat ionic cond coefficient of s2 on T'/T of s1 (rho**2/s)]
3439  ! tabtr(19) : sum heat ionic convection velocity (rho/s)
3440  !--------------------------------------------------------------------
3441  tabtr(14,ipr) = dn_s(1)
3442  ! ITM ------------------------->
3443  neoclassic(itp1)%ne_neo%diff_eff(ipr) = tabtr(14,ipr) * rhomax ** 2 / grho2_pr(ipr)
3444  ! ITM ------------------------->
3445  tabtr(15,ipr) = vnnt_s(1)+vneb_s(1)+gfl_s(5,1)/den(ipr,1)&
3446  & +vnex_s(1)
3447  ! ITM ------------------------->
3448  neoclassic(itp1)%ne_neo%vconv_eff(ipr) = - tabtr(15,ipr) * rhomax / grho2_pr(ipr)
3449  ! ITM ------------------------->
3450  tabtr(16,ipr) = chit_ss(1,1) + chip_ss(1,1)
3451  ! ITM ------------------------->
3452  neoclassic(itp1)%te_neo%diff_eff(ipr) = tabtr(16,ipr) * rhomax**2 / grho2_pr(ipr) * den(ipr,1)
3453  !! [DPC:2010-08-21] I think the above multiplication by "den" is wrong
3454  neoclassic(itp1)%te_neo%diff_eff(ipr) = tabtr(16,ipr) * rhomax**2 / grho2_pr(ipr)
3455  ! ITM ------------------------->
3456  tabtr(17,ipr) = tabtr(10,ipr)/&
3457  & (den_iz(1,1)*temp_i(1)*z_j7kv)+&
3458  & tabtr(16,ipr)*grt_i(1)/temp_i(1)
3459  ! ITM ------------------------->
3460  neoclassic(itp1)%te_neo%vconv_eff(ipr) = - tabtr(17,ipr) * rhomax / grho2_pr(ipr)
3461  ! ITM ------------------------->
3462  tabtr(18,ipr) = 0.0
3463  dentot = 0.0
3464  DO ima=2,m_i
3465  chitp = chit_ss(ima,ima)+chip_ss(ima,ima)
3466  tabtr(18,ipr) = tabtr(18,ipr)+&
3467  & chitp*&
3468  & den(ipr,ima)
3469  vconv = tabtr(11,ipr)/(den(ipr,ima)*temp_i(ima)*z_j7kv)+chitp*grt_i(2)/temp_i(2)
3470  ! ITM ------------------------->
3471  neoclassic(itp1)%ti_neo%diff_eff(ipr,ima-1) = chitp*rhomax**2.0/grho2_pr(ipr)
3472  neoclassic(itp1)%ti_neo%vconv_eff(ipr,ima-1) = - vconv * rhomax / grho2_pr(ipr)
3473  ! ITM ------------------------->
3474  dentot = dentot+den(ipr,ima)
3475  ENDDO
3476  tabtr(18,ipr) = tabtr(18,ipr)/dentot
3477  tabtr(19,ipr) = tabtr(11,ipr)/&
3478  & (dentot*temp_i(2)*z_j7kv)+&
3479  & tabtr(18,ipr)*grt_i(2)/temp_i(2)
3480  ! donnees ioniques
3481  tabtr(20,ipr) = dentot
3482  tabtr(21,ipr) = temp_i(2)
3483  tabtr(22,ipr) = grt_i(2)
3484  ! donnees electoniques
3485  tabtr(23,ipr) = temp_i(1)
3486  tabtr(24,ipr) = den_iz(1,1)
3487  tabtr(25,ipr) = grt_i(1)
3488  !-----------------------------------------------
3489  ! <J_OH.B> Ohmic current [A*T/m**2]
3490  !-----------------------------------------------
3491  tabtr(26,ipr) = p_jboh
3492  !-------------------------------------------------------|
3493  ! UPAR_S(3,m,s)-parallel flow of s from force m [T*m/s]|
3494  ! m=1, p', T', Phi' |
3495  ! m=2, <E.B> |
3496  ! m=3, fex_iz |
3497  !-------------------------------------------------------|
3498  indsave = 26
3499  DO m=1,m_i
3500 
3501  tabtr(indsave+m,ipr) = upar_s(1,1,m) + upar_s(1,2,m) &
3502  & + upar_s(1,3,m)
3503 
3504  ENDDO
3505  indsave = indsave + m_i
3506  !---------------------------------------------------------|
3507  ! UTHETA_S(3,m,s)-poloidal flow of s from force m [m/s/T]|
3508  ! m=1, p', T' |
3509  ! m=2, <E.B> |
3510  ! m=3, fex_iz |
3511  !---------------------------------------------------------|
3512 
3513  DO m=1,m_i
3514  tabtr(indsave+m,ipr) = utheta_s(1,1,m) + utheta_s(1,2,m)&
3515  & + utheta_s(1,3,m)
3516  if (ipr .eq. 1) then
3517  ! write(*,*) 'm=',m,UTHETA_S(1,1,m),UTHETA_S(1,2,m)
3518  ! write(*,*) UTHETA_S(1,3,m)
3519  endif
3520  ENDDO
3521  indsave = indsave + m_i
3522  !
3523  ! coefficient de diffusion par espece et etat de charge de la matiere
3524  !
3525  DO m=2,m_i
3526  tabtr(indsave+m,ipr) = dn_s(m)
3527  neoclassic(itp1)%ni_neo%diff_eff(ipr,m-1) = dn_s(m) * rhomax ** 2 / grho2_pr(ipr)
3528  ENDDO
3529  !
3530  ! Vitesse totale de convection par espece et etat de charge
3531  !
3532  indsave = indsave + m_i -1
3533  DO m=2,m_i
3534  tabtr(indsave+m,ipr) = vnnt_s(m)+vneb_s(m)+&
3535  & gfl_s(5,m)/den(ipr,m)+vnex_s(m)
3536  neoclassic(itp1)%ni_neo%vconv_eff(ipr,m-1) = - tabtr(indsave+m,ipr) * rhomax / grho2_pr(ipr)
3537  ENDDO
3538 
3539  ! -- fin boucle sur les profils --
3540  ! write(*,*) 'fin boucle profil ------------------------------'
3541 1000 ENDDO
3542  ! -- lissage au centre --
3543  etatemp(2:4) = 1.0_8/neoclassic(itp1)%sigma(2:4)
3544  etatemp(1) = 0.0d0
3545  call cos_zconversion(etatemp,4)
3546  neoclassic(itp1)%sigma(1) = 1.0_8/etatemp(1)
3547  neoclassic(itp1)%jboot(1) = 0.0_8
3548 
3549 END SUBROUTINE neo_calc
nclass_itm::nclass
subroutine nclass(m_i, m_z, p_b2, p_bm2, p_eb, p_fhat, p_fm, p_ft, p_grbm2, p_grphi, p_gr2phi, p_ngrth, amu_i, grt_i, temp_i, den_iz, fex_iz, grp_iz, ipr, iflag, L_BANANA, L_PFIRSCH, L_POTATO, K_ORDER, C_DEN, C_POTB, C_POTL, P_ETAP, P_JBBS, P_JBEX, P_JBOH, M_S, JM_S, JZ_S, BSJBP_S, BSJBT_S, GFL_S, DN_S, VNNT_S, VNEB_S, VNEX_S, DP_SS, DT_SS, UPAR_S, UTHETA_S, QFL_S, CHI_S, VQNT_S, VQEB_S, VQEX_S, CHIP_SS, CHIT_SS, CALM_I, CALN_II, CAPM_II, CAPN_II, YMU_S, SQZ_S, XI_S, TAU_SS)
Definition: neo.f90:204
neo_calc
subroutine neo_calc(equilibrium, coreprof, neoclassic, nrho, neq, nion)
Definition: neo.f90:2856
nclass_itm
NCLASS-Calculates NeoCLASSical transport properties
Definition: neo.f90:78
l3interp
subroutine l3interp(y_in, x_in, nr_in, y_out, x_out, nr_out)
Definition: l3interp.f90:1
cos_zconversion
subroutine cos_zconversion(e, nbrho)
Definition: cos_zconversion.f90:1
cos_util
cos_util
Definition: cos_util.f90:8
spec_kind_mod
For REAL variables: Use p=6, r=35 for single precision on 32 bit machine Use p=12,r=100 for double precision on 32 bit machine for single precision on 64 bit machine Parameters for SELECTED_REAL_KIND: p -number of digits r -range of exponent from -r to r
Definition: neo.f90:15
neo
subroutine neo(equilibrium, coreprof, neoclassic)
Definition: neo.f90:2817
cos_rpdederive
subroutine cos_rpdederive(dydx, K, x, y, c0, c1, d, o)
Definition: cos_rpdederive_interface.f90:56