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

	SUBROUTINE inter(m_ou_r,bp,q,qt,n,knot,x,f,dfdx,r_stat,m_stat)

c routine public du module mod_variables

c en (x=m^2/3 ou x=r^2) ou bien (x=m^1/3 ou x=r)
c on fait une interpolation inverse x-->q-->f
c avec m_ou_r=m la localisation de x en indice q0 est faite dans
c m_stat=m^2/3 ou m_stat=m^1/3
c avec m_ou_r=r la localisation de x en indice q0 est faite dans
c r_stat=r^2 ou r_stat=r

c adaptation de l'algorithme de Brent Num. Receip. p.354

c Auteur: P.Morel, Département J.D. Cassini, O.C.A., CESAM2k

c ATTENTION : il faut initialiser et utiliser les bons bp, m_stat et r_stat !!!

c entrées
c	m_ou_r = 'm' : en m^2/3 ou m^1/3, m_ou_r = 'r' : en r^2 ou r
c	ne, bp,  : nb de fonctions (6)
c	bp,q,qt,n,m,knot : spline a interpoler
c	x : point d'interpolation
c	r_stat, m_stat : abscisses

c sorties
c	f, dfdx : valeurs et dérivées /m^2/3 ou m^1/3 si 'm' ou /r^2
c	ou /r si 'r'
c	ATTENTION f(6)=q_int : point d'interpolation

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

	USE mod_donnees, ONLY: ne, ord_qs
	USE mod_kind
	USE mod_numerique, ONLY : bsp1dn, linf

	IMPLICIT NONE

	REAL (kind=dp), INTENT(in) :: x
	INTEGER, INTENT(in) :: n
	CHARACTER (len=1), INTENT(in) :: m_ou_r

	REAL (kind=dp), INTENT(inout), DIMENSION(:,:) :: bp
	REAL (kind=dp), INTENT(inout), DIMENSION(:) :: m_stat, q, qt, r_stat
	INTEGER, INTENT(inout) :: knot
	REAL (kind=dp), INTENT(out), DIMENSION(:) :: f, dfdx

        REAL (kind=dp), PARAMETER :: epsi=1.d-6
	REAL (kind=dp) :: corr, qi, qs, q_int
	INTEGER, PARAMETER :: ntour_max=40

	INTEGER :: i, ntour, inc, l=1

	LOGICAL :: iterate

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

2000	FORMAT(8es10.3)

	SELECT CASE(m_ou_r)
	CASE('m')		!interpolation inverse en m^2/3 ou m^1/3
	 inc=5	!indice de l'inconnue pour l'interpolation inverse
	 IF(x <= m_stat(1))THEN
	  q_int=1.d0 ; iterate=.FALSE.
	 ELSEIF(x >= m_stat(n))THEN
	  q_int=n ; iterate=.FALSE.
	 ELSE
	  CALL linf(x,m_stat,n,l)		!entre q_int et q_int+1
	  q_int=l
	  IF(x == m_stat(l))THEN
	   iterate=.FALSE.
	  ELSEIF(x == m_stat(l+1))THEN
	   q_int=l+1 ; iterate=.FALSE.
	  ELSE
	   iterate=.TRUE.
	  ENDIF
	 ENDIF

c on précise la localisation à mieux que 0.1% de l'écart
c	 IF(iterate)epsi=(1.d0-m_stat(q_int)/m_stat(q_int+1))*1.d-3

	CASE('r')			!interpolation inverse en r^2 ou r
	 inc=3
	 IF(x <= r_stat(1))THEN
	  q_int=1.d0 ; iterate=.FALSE.
	 ELSEIF(x >= r_stat(n))THEN
	  q_int=n ; iterate=.FALSE.
	 ELSE
	  CALL linf(x,r_stat,n,l)		!entre q_int et q_int+1
	  q_int=l
	  IF(x == r_stat(l))THEN
	   iterate=.FALSE.
	  ELSEIF(x== r_stat(l+1))THEN
	   q_int=l+1 ; iterate=.FALSE.
	  ELSE
	   iterate=.TRUE.
	  ENDIF
	 ENDIF

c on précise la localisation à mieux que 0.1% de l'écart
c	 IF(iterate)epsi=(1.d0-r_stat(q_int)/r_stat(q_int+1))*1.d-3

	CASE DEFAULT
	 WRITE(*,1)m_ou_r ; WRITE(2,1)m_ou_r ; CALL sortie('inter')
1	 FORMAT('ARRET appel a inter avec m_stat_ou_r=',a3)
	END SELECT

c hors des extémités
	IF(iterate)THEN		!dichotomie
	 l=q_int ; corr=1 ; ntour=0
	 qi=q_int	!q_int inferieur
	 qs=q_int+1	!q_int superieur

c détermination de l'indice par dichotomie pour interpolation inverse
	 B1: DO
	  q_int=(qs+qi)/2.d0
	  CALL bsp1dn(ne,bp,q,qt,n,ord_qs,knot,.TRUE.,q_int,l,f,dfdx)
	  IF(f(inc) == x)THEN
	   qi=q_int
	   qs=q_int
	  ELSEIF(f(inc) < x)THEN
	   qi=q_int
	  ELSE
	   qs=q_int
	  ENDIF
	  corr=qs-qi	!qs >= qi ==> corr >= 0
  	  IF(corr <= epsi)EXIT B1
	  ntour=ntour+1
	  IF(ntour >= ntour_max)THEN
	   WRITE(*,2)m_ou_r,q_int,corr,epsi,x,x-f(inc)
2	   FORMAT('pas de conv. dicho. dans inter en',a3,/,
	1  'q_int=',es10.3,', corr=',es10.3,', epsi=',es10.3,', x=',es10.3,
	2  ', x-f(inc)=',es10.3)
	   EXIT B1
	  ENDIF
c	  PRINT*,'corr,x-f(inc),qi,qs,x,f(inc)'
c	  WRITE(*,2000)corr,x-f(inc),qi,qs,x,f(inc) !; PAUSE'inter1'
	 ENDDO B1
	ENDIF	!iterate

c interpolation
	CALL bsp1dn(ne,bp,q,qt,n,ord_qs,knot,.TRUE.,q_int,l,f,dfdx)

c dérivées par rapport à la variable d'interpolation inverse
c inc=3 R^2 ou R, inc=5 M^2/3 ou M^1/3
	DO i=1,ne
	 IF(i /= inc)dfdx(i)=dfdx(i)/dfdx(inc)
	ENDDO

c astuce pour sortir l'indice: en sortie f(6)=q et non la fonction
c de répartition
	f(6)=q_int

	RETURN

	END SUBROUTINE inter