c-----------------------------------------------------------------------

	subroutine gr_parti3(t,nel,fparti,ff,j)
	
c routine PRIVATE du module mod_bp_for_alecian_new	

! Determination des fonctions de partition de tous les ions de j,
! a une temperature donnee. Methode simple, suffisante pour le calcul des
! accelerations, a n'utiliser qu'avec les niveaux d'energies lus avec niv.
! Troncage des niveaux en fonction du rayon de Debye.

!  Auteur:
!     Georges ALECIAN
!     LUTH - UMR 8102, Observatoire de Meudon
!     F-92195 MEUDON CEDEX, FRANCE
!     Tel: 01 45 07 74 20, + 33 1 45 07 74 20

c Adaptation à CESAM2k B.Pichon, P.Morel

c----------------------------------------------------------------------

	implicit none
	
	integer :: j
	real(DP) :: t, nel

! output du subroutine
	real(DP), dimension(0:pzi)     :: fparti
	real(DP), dimension(30,0:pzi)  :: ff

! variables locales
	integer                             :: k,m,ks
	real, parameter                     :: Ze=4.8032044E-10,secu=100.
	real, parameter                     :: ajust=5.
	real(DP)                :: tke,Rd,x_C,x_cesam
c	real(DP), dimension(30) :: f,xx

!=========== initialisation tableaux
	fparti=0.
!===========
	
	Rd = SQRT(23.8112260522911 * t / nel)
! Rayon Debye (suppose ionisation complete, et nb de charges = 2*ne)
! 23.8112260522911=kB/(8*pi*(e^2)), approximation "livre, Eq.4.10 "eq:coefRD"

	tke=t * 8.6173431838E-05     ! temperature en eV
	


	ff = 1.  ! initialisation tableau
	do k=0,nint(zi(j)) ! on boucle sur les ions (charge k) de j

		ks = niv_nb(k,j)
c		xx= 0.d0
		if(k.eq.0) then         ! test sur la charge de l'ion
			x_C = 0.001*1.6E-12  ! pour imposer x_cesam>4.
		else
			x_C = real(k)*(Ze**2)/Rd
		endif
		
		do m=1,ks   ! on boucle sur les niveaux de chaque ion k

			if(niv_e(m,k,j)/tke.gt.secu) CYCLE   ! m suivant

			x_cesam = ( ( (el_pot(k,j)-niv_e(m,k,j)) *
     +		          1.6E-12/x_C)/ajust - 1.)

!	xx(m)=x_cesam
			if(x_cesam.ge.4.) then
				ff(m,k) = 1.
			else if(x_cesam.lt.0.) then
				ff(m,k) = 0.
			else
				ff(m,k) = ((x_cesam/4.)**2)*(-2*(x_cesam/4.)+3.)
			end if

			if(m.eq.1) then
				fparti(k) = fparti(k) + 
     +			            niv_g(m,k,j) * exp(-niv_e(m,k,j)/tke)
			else
				fparti(k) = fparti(k) + ff(m,k) *
     +			            niv_g(m,k,j) * exp(-niv_e(m,k,j)/tke)
			end if



		end do


	end do
	
	end subroutine gr_parti3