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

	subroutine gr_theta(t,fparti,ff,theta,j)
	
c routine du module mod_bp_for_alecian_new

! Version 1.1 (05/03/2007)
! Determination de la fonction Theta pour le g_cont de tous les ions de j,
! a une temperature donnee

!  Cette subroutine fait partie du module grcesam2.
!  Version 2.0 (7/01/2004)
!  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
	real(DP), dimension(0:pzi)        :: fparti
	real(DP), dimension(30,0:pzi)     :: ff

! output du subroutine
	real(DP), dimension(0:pzi) :: theta

! variables locales
	integer                             :: k,m,ks
c	real, parameter                     :: secu=100.
	real(DP)                :: tke,uk,Qk

!=========== initialisation tableaux
	theta=0.
!===========
	

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

	do k=1,nint(zi(j)) ! on boucle sur les ions de j (sauf le premier)

		ks = niv_nb(k,j)
		if(ks.ge.1) then
			do m=1,ks   ! on boucle sur les niveaux de chaque ion k
				uk  = (el_pot(k-1,j)-niv_e(m,k-1,j))/tke
				if(uk.lt.1.e-2)  CYCLE
				if(uk.gt.20.)    CYCLE
				
				Qk=(uk**3)*(uk/(1.-exp(-uk)) - exp(uk)*log(1.-exp(-uk)))
				Qk= ff(m,k-1) * Qk

				theta(k) = theta(k) + 
     +	              niv_q(m,k-1,j)*(niv_g(m,k-1,j)/niv_g(1,k-1,j)) * Qk
	
			end do
			theta(k)=niv_g(1,k-1,j)/fparti(k)*theta(k)
		end if
	end do

	end subroutine gr_theta