c********************************************************************

	SUBROUTINE cond_mestel(xchim1,t,ro,kappa_c)

c subroutine PRIVATE du module mod_opa

c opacité conductive selon Mestel rapporté par Cox & Guili §16.7, p.384

c entrées:
c	xchim1 : comp. chim. en fraction de masse
c	t : température K
c	ro : densité cgs

c sortie:
c	kapa : opacité gr / cm2)
c	dkapdt : kappa / d t
c	dkapdr : kappa / d densité
c       dkapdx : kappa / d xchim(1)

c Auteur: P. Morel, Département Cassiopée, O.C.A. CESAM2k

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

	USE mod_donnees, ONLY : amu, hpl, kbol, langue, me, mu_saha, nchim,
	1 nucleo, pi, zi
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn, ferdir => fermi_dirac,
	1 inside, no_croiss

	REAL (kind=dp), INTENT(in), DIMENSION(nchim) :: xchim1
	REAL (kind=dp), INTENT(in) :: t, ro
	REAL (kind=dp), INTENT(out) :: kappa_c

	INTEGER, PARAMETER :: m=2, n=16, nft=2, pf12=101
	REAL (kind=dp), SAVE, DIMENSION(nft,n) :: ft
	REAL (kind=dp), SAVE, DIMENSION(1,pf12) :: dgce
	REAL (kind=dp), SAVE, DIMENSION(pf12+m) :: f12t
	REAL (kind=dp), SAVE, DIMENSION(pf12) :: f12
	REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: zi13, zi2
	REAL (kind=dp), SAVE, DIMENSION(n+m) :: eta_t
	REAL (kind=dp), SAVE, DIMENSION(n) :: eta
	REAL (kind=dp), DIMENSION(nchim) :: xchim, teta_i
	REAL (kind=dp), DIMENSION(nft) :: dfxdx, fx
	REAL (kind=dp), DIMENSION(nft*n) :: bid
	REAL (kind=dp), DIMENSION(5) :: fd

	REAL (kind=dp), SAVE :: cte1, cte2
	REAL (kind=dp) :: degene, nel, t7


	INTEGER, SAVE :: dep, fin, knot, kno, l=m

	LOGICAL, SAVE :: init=.TRUE.

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

2000	FORMAT(8es10.3)

	IF(init)THEN	!initialisations
	 init=.FALSE.

c avec les opacités conductives de Mestel mu_saha=.TRUE. pour la dégérescence
c	 mu_saha=.TRUE.

c allocations
	 ALLOCATE(zi13(nchim),zi2(nchim))

c tabulation de ft
	 eta=(/ -4.0d0, -3.0d0, -2.0d0, -1.0d0, -0.2d0, +0.4d0, +1.0d0, +1.6d0,
	1       +2.2d0, +2.8d0, +3.4d0, +4.0d0, +5.0d0, +6.0d0, +7.0d0, +8.0d0 /)
	 bid=(/ 0.15479d0, 23.85d0, 0.21466d0, 23.81d0, 0.29766d0, 23.53d0,
	1       0.39518d0, 22.79d0, 0.48637d0, 21.58d0, 0.54945d0, 20.17d0,
	2       0.60840d0, 18.30d0, 0.65825d0, 16.04d0, 0.69820d0, 13.57d0,
	3       0.72929d0, 11.09d0, 0.75305d0, 8.770d0, 0.77132d0, 6.738d0,
	4       0.79263d0, 4.100d0, 0.80670d0, 2.349d0, 0.81615d0, 1.282d0,
	5       0.82284d0, 0.6713d0 /)
	 ft=RESHAPE(bid,SHAPE(ft))
	 CALL bsp1dn(nft,ft,eta,eta_t,n,m,knot,.FALSE.,eta(2),l,fx,dfxdx)
	 IF(no_croiss)PRINT*,'Pb. at 1 in cond_mestel'

c tabulation des fonctions F12 de Fermi-Dirac de dep à fin
c signe opposé de celui de Clayton
	 dep=-30.d0 ; fin=30.d0 ; t7=(fin-dep)/REAL(pf12-1,dp)	!vT
	 DO i=1,pf12
	  dgce(1,i)=dep+t7*(i-1) ; CALL ferdir(dgce(1,i),fd) ; f12(i)=fd(2)
	 ENDDO
	 CALL bsp1dn(1,dgce,f12,f12t,pf12,m,kno,.FALSE.,f12(2),l,fx,dfxdx)
	 IF(no_croiss)PRINT*,'Pb. at 2 in cond_mestel'

c constantes
	 cte1=pi**2/3.d0 ; cte2=SQRT(hpl/2.d0/me/kbol*hpl)**3/4.d0/pi
	 zi2=zi**2 ; zi13=ABS(zi)**(1.d0/3.d0)

	 SELECT CASE(langue)
	 CASE('english')
	  WRITE(*,1001) ; WRITE(2,1001)
1001	  FORMAT('Conductives Opacities of Mestel, Cox & Guili §16.7, p.384',/)
	 CASE DEFAULT
	  WRITE(*,1) ; WRITE(2,1)
1	  FORMAT('Opacités conductives de Mestel (Cox & Guili §16.7, p.384',/)
	 END SELECT
	ENDIF

c-----------------------initialisations fin-------------------------

c calcul direct de la dégénérescence, cas totalement ionisé
c t7 vT, SQRT(hpl/2.d0/me/kbol*hpl)**3/4.d0/pi*nel/SQRT(t)**3
	xchim=xchim1/nucleo ; nel=DOT_PRODUCT(zi,xchim)*ro/amu
	t7=inside(f12(1),f12(pf12),cte2*nel/SQRT(t)**3)
	CALL bsp1dn(1,dgce,f12,f12t,pf12,m,kno,.TRUE.,t7,l,fx,dfxdx)
	IF(no_croiss)PRINT*,'Pb. at 3 in cond_mestel'
	degene=fx(1)

c opacité conductive
	t7=t*1.d-7
	IF(degene <= -4.d0)THEN
	 teta_i=0.589d0/zi13*EXP(degene/3.d0)		!eq 16.130
	 teta_i=LOG(SQRT(2.d0/(1.d0-COS(teta_i))))	!16.127
	 kappa_c=48.2d0/t7*SUM(zi2*xchim*teta_i)/EXP(degene)	!16.129
	ELSEIF(degene >= 8.d0)THEN
	 teta_i=0.848d0/zi13*(1.d0-2.06d0/degene**2)	!eq. 16.137
	 teta_i=LOG(SQRT(2.d0/(1.d0-COS(teta_i))))	!16.127
	 fx(2)=cte1*(degene**2*(degene**2*(degene+21.71d0)+136.4d0)-628.1d0)/
	1 (degene*(degene**2+9.870d0))		!eq. 16.134
	 kappa_c=5.12d-3*SUM(zi2*xchim*teta_i)*(t7/(ro*1.d-5)/(1.d0+xchim(1)))**2
	ELSE
	 CALL bsp1dn(nft,ft,eta,eta_t,n,m,knot,.TRUE.,degene,l,fx,dfxdx)
	 IF(no_croiss)PRINT*,'Pb. at 4 in cond_mestel'
	 fx(2)=fx(2)*EXP(degene)	!f(eta) tabulation16.3
	 teta_i=fx(1)/zi13		!tabulation16.3
	 teta_i=LOG(SQRT(2.d0/(1.d0-COS(teta_i))))	!16.127
	 kappa_c=1.158d3/t7/fx(2)*SUM(zi2*xchim*teta_i)		!16.126
	ENDIF
c	PRINT*,'t,nel,degene,kappa_c' ; WRITE(*,2000)t,nel,degene,kappa_c

	RETURN

	END SUBROUTINE cond_mestel