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

	SUBROUTINE calc_grad(qi,grad,gradad,gradldx,gradrad,grad_mu)

c routine public du module mod_static
c calcul des gradients

c Auteur: P.Morel Laboratoire Lagrange O.C.A., CESAM2k

c entrée :
c	qi: indice REAL du point de calcul

c sorties:
c	 grad,gradad,gradldx,gradrad

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

	USE mod_donnees, ONLY : aradia, clight, cpturb, en_m23, g, ihe4, Krot, ledoux,
	1 lsol, msol, m_ch, nchim, ne, nrot, ord_qs, ord_rot, pi, psi0, rsol, t_inf
	USE mod_etat, ONLY : etat, mu_mol
	USE mod_kind
	USE mod_opa, ONLY : opa, dehors
	USE mod_nuc, ONLY : nuc
	USE mod_numerique, ONLY : bsp1dn, no_croiss
	USE mod_variables, ONLY : bp, chim, chim_gram, dim_qs, dv_tr,
	1 id_conv, if_conv, inter, jlim, lconv, knot, knotc, knotr,
	2 mc, mct, mrot, mrott, mstar, m_zc, m_zc23, m_stat,
	3 nc_fixe, n_ch, n_qs, n_rot, omega_jpz, q, qt, rota,
	4 rstar, r_stat, r_ov, r_zc, r_zc_conv, wrot	!, sortie

	IMPLICIT NONE

	REAL (kind=dp), INTENT(in) :: qi
	REAL (kind=dp), INTENT(out) :: grad, gradad, gradldx, gradrad, grad_mu

	REAL (kind=dp), DIMENSION(0,0) :: jac
	REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: dfrot, frot
	REAL (kind=dp), DIMENSION(nchim) :: depsx, dlnmu_x, dgrad_mux,
	1 dmu_x, dxchim, dxchimm, xchim, xchimm
	REAL (kind=dp), DIMENSION(ne) :: df, f
	REAL (kind=dp), SAVE, DIMENSION(5) :: epsilon
	REAL (kind=dp), SAVE :: cte1, cte2, cte8, cte9, cte13
	REAL (kind=dp) :: alfa, beta, be7, b8, cp, dcpp, dcpt, dcpx, delta,
	1 deltap, deltat, deltax, depst, depsro, dgradadp, dgradadt,
	3 dgradadx, dkapro, dkapt, dkapx, dlnmu, drop, drot, drox,
	2 dup, dut, dux, f17, gamma1, gravite, hh, hp, kap,
	4 krad, l, m, mu, mue, nu, n13, o15, p,t, r, ro, u, w, Zbar

	INTEGER, SAVE :: lq=1, lc=1

	LOGICAL, SAVE :: init=.TRUE.

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

2000	FORMAT(8es10.3)

	IF(init)THEN	!initialisations
	 init=.FALSE.
	 cte1=4.d0/3.d0*aradia*clight
	 cte13=g*msol/rsol/rsol
	 cte2=2.d0/3.d0*rsol
	 cte8=lsol/4.d0/pi/rsol/rsol	!de 5.9
	 cte9=3.d0/16.d0/pi/aradia/clight/g
	ENDIF		!initialisation

c les variables en qi
 	CALL bsp1dn(ne,bp,q,qt,n_qs,ord_qs,knot,.TRUE.,qi,lq,f,df)
	IF(no_croiss)PRINT*,'Problème localisé en 1 dans calc_grad'
c	WRITE(*,2000)f

c le gradient
	grad=df(2)/df(1)
	p=EXP(f(1))
	t=EXP(f(2))
	IF(en_m23)THEN
	 r=SQRT(ABS(f(3))) ; l=SQRT(ABS(f(4)))**3 ; m=SQRT(ABS(f(5)))**3
	ELSE
	 r=ABS(f(3)) ; l=f(4) ; m=ABS(f(5))**3 ; f(5)=f(5)**2
	ENDIF
c	PRINT*,'qi f(5)',qi,f(5)

c composition chimique
	CALL bsp1dn(nchim,chim,mc,mct,n_ch,m_ch,
	1 knotc,.TRUE.,MAX(mc(1),MIN(f(5),mc(n_ch))),lc,xchim,dxchim)
	IF(no_croiss)PRINT*,'Problème localisé en 2 dans calc_grad'

c vitesse angulaire
	SELECT CASE(Krot)
	CASE(0,1,2)
	 w=wrot
	CASE(3,4,5)
	 CALL bsp1dn(nrot,rota,mrot,mrott,n_rot,ord_rot,knotr,.TRUE.,
	1 MAX(mrot(1),MIN(f(5),mrot(n_rot))),lc,frot,dfrot)
	 IF(no_croiss)PRINT*,'Problème localisé en 6 dans lim_zc'
	 w=frot(1)
	CASE(6)
	 CALL omega_jpz(r,w)
	END SELECT

c composition chimique /gr
	xchimm=ABS(xchim) ; dxchimm=dxchim
	CALL chim_gram(xchimm,dxchimm)
	CALL etat(p,t,xchimm,.FALSE.,
	1 ro,drop,drot,drox,u,dup,dut,dux,
	2 delta,deltap,deltat,deltax,cp,dcpp,dcpt,dcpx,
	3 gradad,dgradadp,dgradadt,dgradadx,alfa,beta,gamma1)

c opacité avec test sur sortie de table
	CALL opa(xchimm,t,ro,kap,dkapt,dkapro,dkapx)
	IF(dehors)THEN
	 PRINT*,'appel de calc_grad, / P, T, Ro / xchimm'
	 WRITE(*,2000)p,t,ro ; WRITE(*,2000)xchimm
	 RETURN
	ENDIF

	krad=cte1/kap/ro*t**3		!5.1 conductivite radiative
	IF(m*l*r /= 0.d0)THEN		!gradient radiatif
	 gravite=cte13*m/r**2-cte2*w**2*r !gravité effective avec rotation
	 gradrad=cte8*l*p/gravite/ro/r**2/krad/t	!5.9
	 hp=p/gravite/ro		!échelle de hauteur de pression
	ELSE		!au centre
	 IF(t > t_inf)THEN
	  CALL nuc(t,ro,xchim,dxchim,jac,.FALSE.,3,
	1 epsilon,depst,depsro,depsx,hh,be7,b8,n13,o15,f17)
	 ELSE
	  epsilon(1)=0.d0	!total
	 ENDIF
	 gradrad=cte9*kap*epsilon(1)*p/t**4	!au centre l/m ~ epsilon
	 hp=1.d38
	ENDIF

c gradient de Ledoux pour la convection
	nu=m**(2.d0/3.d0)
	CALL mu_mol(dxchim,hp,nu,r,ro,t,xchim,dlnmu,dlnmu_x,grad_mu,
	1 dgrad_mux,mu,dmu_x,mue,Zbar)
	gradldx=gradad+beta/(4.d0-3.d0*beta)*grad_mu
c	PRINT*,'calc_grad',qi,lq,lc,mct(lc),mct(lc+1)
c	WRITE(*,2000)grad,gradad,gradldx,gradrad,grad_mu

	RETURN

	END SUBROUTINE calc_grad