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

      SUBROUTINE opa_houdek12(xh,t,ro,kappa,dkapdt,dkapdr,dkapdx)

c routine private du module mod_opacite

c interpolation from Guenter Houdek: bi-rational splines version12
c This routine calls the file SUN_STAR_DATA/HOUDEK/OPINTPATH
c which gives the path to the opacity tables
c The routine also calls three routines:
c maceps (optimization for the compiler)
c opinit
c opints  for the interpolation
c The compilation needs 1 more libraries: opint,
c linker avec :~/SUN_STAR_DATA/HOUDEK12/lib -lopint


c le nom de la table d'opacité OPINTPATH doit être
c indiqué dans la NAMELIST NOM_OPA, du fichier de données mon_modele.don
c Exemple : NOM_OPA='opa_houdek12',
c           F_OPA='opa_yveline.bin','H12/OPINTPATH',6*' ',
c	    NOM_OPA_COND=' '
c l'opacité conductive est incluse dans la valeur délivrée par opa_houdek12
c un flag permettant de l'éliminer est décrit dans la notice de Guenter

c	http://users-phys.au.dk/hg62/opint.html

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

c entrées :
c     xchim(1)=X : comp. chim. en H
c     t : température K
c     ro : densité cgs

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

c     Z est obtenu par 1-X-Y, Y=He4+He3 ou encore Z=Z0

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

      USE mod_donnees, ONLY : f_opa, ihe4, nchim, nom_chemin, z0
      USE mod_kind

      IMPLICIT NONE

      REAL (kind=dp), INTENT(in), DIMENSION(:) :: xh
      REAL (kind=dp), INTENT(in) :: t, ro
      REAL (kind=dp), INTENT(out) ::  kappa, dkapdt, dkapdr, dkapdx

      REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: Lxh
      REAL (kind=dp) :: tlg, rlg, opalg, opr, opt, opx, opz, x, z,
     1 eps, Lt, Lro

      INTEGER, PARAMETER :: iord=4, imode=3
      INTEGER :: iexp, ier

      LOGICAL, SAVE :: init=.TRUE., hxp=.TRUE., hrp=.TRUE., hzp=.TRUE.,
     1 htp=.TRUE., hrm=.TRUE., htm=.TRUE.
      LOGICAL :: hou, ok

      CHARACTER (len=80) :: tabnam

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

2000  FORMAT(8es10.3)

	IF(init)THEN         ! reads the tables at the first call
	 init=.FALSE.

	WRITE(*,1)TRIM(f_opa(2)) ; WRITE(2,1)TRIM(f_opa(2))
1	FORMAT('opacités à Z variable, tables de OPAL 95 :',/,
	1 'interpolation par splines birationnelles de G. Houdek',/,
	2 'version 12 adaptation à CESAM2k de P.Morel, données de ',a)

c lecture et initialisation des tables
       tabnam=TRIM(nom_chemin)//TRIM(f_opa(2))
       INQUIRE(file=TRIM(tabnam),exist=ok)
       IF(.NOT.ok)THEN
        WRITE(*,2)TRIM(f_opa(2)) ; WRITE(2,2)TRIM(f_opa(2))
2	FORMAT('ARRET, le nom de la table d''opacité Houdek12 est inconnu : ',a)
	STOP
       ENDIF

c PRINT*,tabnam ; CALL pause('tabnam')
       CALL maceps(eps)   !precision machine
       CALL opinit(eps,iord,tabnam,imode) !lecture des tables

       x=0.7d0 ; z=.02d0 ; tlg=4.d0 ; rlg=0.d0
       CALL opints(x,z,tlg,rlg,opalg,opr,opt,opx,opz,iexp,ier)
       PRINT* ; PRINT*,'fin des initialisations des opacités Houdek 95'
       PRINT*
       ALLOCATE(Lxh(nchim))
      ENDIF

c transformation of the inputs of OPA  to those of OPINTF:
      x=ABS(xh(1))
      IF(nchim > 1)THEN
       z=ABS(1.d0-xh(1)-xh(ihe4)-xh(ihe4-1))
      ELSE
       z=z0
      ENDIF

c variables de OPAL en LOG10
      tlg=LOG10(t) ; rlg=LOG10(ro/(t*1.d-06)**3) ; hou=.TRUE.

c encadrement des inputs
      IF(x > 0.8d0)THEN
       IF(x > 0.81d0)THEN
        IF(hxp)THEN
         WRITE(2,11) ; WRITE(*,11)
11       FORMAT(/,'opa_houdek12, on rencontre au moins une fois:',/,
     1   'X > 0.81, appel à opa_yveline')
         hxp=.FALSE.
        ENDIF
        hou=.FALSE.
       ENDIF
       x=0.8d0
      ENDIF

      IF(z > 0.1d0)THEN
       IF(hzp)THEN
        WRITE(2,12) ; WRITE(*,12)
12      FORMAT(/,'opa_houdek12, on rencontre au moins une fois:',/,
     1  'Z > 0.10, appel à opa_yveline')
        hzp=.FALSE.
       ENDIF
       hou=.FALSE.
      ENDIF

      IF(tlg < 3.d0)THEN
       IF(tlg < 2.9d0)THEN
        IF(htm)THEN
         WRITE(2,13) ; WRITE(*,13)
13       FORMAT(/,'opa_houdek12, on rencontre au moins une fois:',/,
     1   'T < 1000, appel à opa_yveline')
         htm=.FALSE.
        ENDIF
        hou=.FALSE.
       ENDIF
       tlg=3.d0
      ELSEIF(tlg > 8.7d0)THEN
       IF(tlg > 8.8d0)THEN
        IF(htp)THEN
         WRITE(2,14) ; WRITE(*,14)
14       FORMAT(/,'opa_houdek12, on rencontre au moins une fois:',/,
     1   'T > 6.4d8, appel à opa_yveline')
         htp=.FALSE.
        ENDIF
        hou=.FALSE.
       ENDIF
       tlg=8.7d0
      ENDIF

      IF(rlg < -8.d0)THEN
       IF(rlg < -8.05d0)THEN
        IF(hrp)THEN
         WRITE(2,15) ; WRITE(*,15)
15       FORMAT(/,'opa_houdek12, on rencontre au moins une fois:',/,
     1   'rlg <  -8, appel à opa_yveline')
         hrp=.FALSE.
        ENDIF
        hou=.FALSE.
       ENDIF
       rlg=-8.d0
      ELSEIF(rlg > 1.d0)THEN
       IF(rlg > 1.5d0)THEN
        IF(hrm)THEN
         WRITE(2,16) ; WRITE(*,16)
16       FORMAT(/,'opa_houdek12, on rencontre au moins une fois:',/,
     1   'rlg > 1, appel à opa_yveline')
         hrm=.FALSE.
        ENDIF
        hou=.FALSE.
       ENDIF
       rlg=1.d0
      ENDIF

c CALL OF THE SUBROUTINE OF G. Houdek:
      IF(hou)THEN
       CALL opints(x,z,tlg,rlg,opalg,opr,opt,opx,opz,iexp,ier)
       IF(ier == 1)THEN
        WRITE(*,*)'opa_houdek12, sortie de table et extrapolation'
       ENDIF

       kappa =10.d0**opalg            ! opacity
       dkapdt =(kappa /t )*opt         ! d kappa / d T
       dkapdr =(kappa /ro )*opr        ! d kappa / d ro
       IF(x > 0.d0)THEN
        dkapdx =(kappa /x )*opx         ! d kappa / d X
       ELSE
        dkapdx =0.d0
       ENDIF

c      WRITE(*,*)'kappa,dkapdr,dkapdt,dkapdx'
c      WRITE(*,2000)kappa,dkapdr,dkapdt,dkapdx
c      CALL pause5('les kappa'

      ELSE
       Lxh(1:nchim)=xh(1:nchim) ; Lt=t ; Lro=ro
       CALL opa_yveline(Lxh,Lt,Lro,kappa,dkapdt,dkapdr,dkapdx)
      ENDIF

      RETURN

      END SUBROUTINE opa_houdek12