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

	subroutine gr_popion3(t,nel,fparti,ff,popi,j)
	
c routine du module mod_bp_for_alecian_new

! Determination des populations de tous les ions de j,
! a une temperature donnee.
! Correction en fonction du rayon de Debye.

!  Cette subroutine fait partie du module grcesam2p1.
!  Version 3.0 (2/03/2007)
!  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---------------------------------------------------------------------

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

! input du subroutine
	real(DP), Intent(In), dimension(0:pzi) :: fparti
	real(DP), Intent(In), dimension(30,0:pzi)  :: ff
	
! output du subroutine
	real(DP), Intent(out), dimension(0:pzi,pnchim) :: popi

! variables locales
	integer :: k,m,na0,na,nb
	real(DP), parameter :: secu=100.
	real(DP)  :: tke,r,d,an,c
	real(DP), dimension(0:pzi) :: ax,p

!=========== initialisations
	r   = 1.d0
	d   = 1.d0
	na0 = 0
	na  = nint(zi(j))
	an  = 1.d0          ! car on se contente des populations relatives
!===========
	
	tke = t * 8.6173431838d-05     ! temperature en eV
c	c = (1.30364e+26 * 1.38e-16 * t)**(-1.5) 	
	c = SQRT(1.30364d+26 * 1.38d-16 * t)**3 ; c=1.d0/c	!**-3/2
		
! 1.30364e+26 = (2pi*me/h**2)**3/2

	do k=0,nint(zi(j))  ! on boucle sur les ions de j
	
		if(niv_nb(k,j) >= 1) then		
			nb=k
		else
			na0=k+1    ! premier ion de la base de donnees atomiques
		endif
	end do

	do k=nb,na0+1,-1  ! on boucle sur les ions existant dans la base
		if(el_pot(k-1,j)/tke > secu) then
			ax(k) =1.d+35
			nb = k-1
		else
			ax(k) = ff(1,k-1)*(nel/2.)*
     +		      (fparti(k-1)/fparti(k))*c*exp(el_pot(k-1,j)/tke)
		end if
	end do
	
! determination de na
	p=0.
	do k=nb,na0+1,-1
		do m=k,na0+1,-1
			d = d*ax(m)
			r = d + r
		end do
		p(k)=1./r
		if(p(k) < 1.e-30) then
			p(k)=0.
			r=1.
			d=1.
		else
			na=k
			exit
		end if
	end do

	do k=na-1,na0,-1
		p(k)=p(k+1)*ax(k+1)
	end do

	do k=nint(zi(j)),0,-1
		popi(k,j)=p(k)*an
	end do

	RETURN
	
	end subroutine gr_popion3