************************************************************************
*
      subroutine TtauIntp(Teff,grv,alpha)
c
c  For given Teff, and (linear) surface gravity, grv, interpolates in the table,
c  TtauFeH0.dat, of mixing-length parameters \alpha (Trampedach et al., 2013b)
c  and radiative Hopf functions, q(\tau_Ross), for computing T(\tau_Ross).
c  See Trampedach et al. (2013a) for derivation and details on implementing the
c  T(\tau_Ross)-relations in stellar structure models.
c    These results have been derived for a X=73.6945% and Z=1.8055% mixture
c  with an Anders & Grevesse (1989) metal mixture, with He and Fe adjusted to
c  the modern values of 10.92 and 7.50. This set of abundances is close to
c  that of Grevesse & Noels (1993).
c    The effects of using, e.g., a scaled solar T(\tau)-relation can be explored
c  by removing all but the Solar entry from the table, 'TtauFeH0.dat', changing
c  the dimensions in line 1 of the table and changing the parameters Mqtau and
c  Msims below.
c
c  Arguments:
c  ----------
c  Teff:      Effective temperature/[K] <input>
c  grv:       Acceleration of gravity at the surface/[cm/s^2] <input>
c  alpha:     The primary mixing-length parameter <output>
c
c  Return Values:
c  --------------
c  Puts the radiative Hopf function in qhpf of common-block chopf, and
c  the related Roseeland \tau-scale in qltau of that common-block.
c  The calibrated mixing-length parameter, alpha, is returned in an
c  argument to the routine.
c
c  Side Effects:
c  -------------
c  Changes the values of qltau and qhpf of common-block chopf.
c
c  Coding history:
c  ---------------
c  04.08.2013: Coded and adopted from tg_intrp_hopf.f/R.Trampedach
c              (trampeda@lcd.colorado.edu)
c
c  Known bugs: Nqtau and Nsims of the table, 'TtauFeH0.dat', needs to
c              be identical to Mqtau and Msims in the parameter statement.
c                Undefined behaviour for Nsims=2; Not a problem, unless
c              the user reprograms to Msims=2, which is a less useful
c              case, anyway.
c
c***********************************************************************
c
c  Double precision version.
c  ++++++++++++++++++++++++
c
      implicit double precision (a-h,o-z)
      parameter(Mqtau=82, Msims=37)
      integer     Nqtau, Nsims
      character*2 sNsims, sNsims1
      common/chopf/ Nqtau, qltau(Mqtau), qhpf(Mqtau), dqhpf(Mqtau)
c
      real*8 Teffsl(Msims), gravsl(Msims), FeH(Msims), alph(Msims)
      real*8 sigalph(Msims), LTeff3(Msims), Logg3(Msims)
      real*8 aa(Msims,Mqtau), aa1(Msims+1,Mqtau)
      integer iadj(6*Msims), iend(Msims), isrt(Msims)
c
      data Nsims /0/
      save Nsims, Teffsl,gravsl,FeH,alph,sigalph,aa,iadj,iend,ier
c
      Teffl  = dlog10(Teff)
      gravl  = dlog10(grv)
c
c***********************************************************************
c  This block is only called once and bypassed by subsequent calls
c  to TtauIntp. This block reads the table and performs sanity checks.
c
      if (Nsims .ne. Msims) then
c
c  Open and read the file
        open (37,file="TtauFeH0.dat",status='old')
        read (37,'(12x,2I12)') Nqtau, Nsims
        if (Nsims.ne.Msims .or. Nqtau.ne.Mqtau) then
           write(6,*) 'Msims,Mqtau: ', Msims, Mqtau
           write(6,*) 'Nsims,Nqtau: ', Nsims, Nqtau
           stop
        endif
c       write(6,*) 'Nsims,Nqtau: ', Nsims, Nqtau
        write(sNsims, '(I2)') Nsims
        write(sNsims1,'(I2)') Nsims+1
        read (37,'(12x,'//sNsims//'f12.3)') Teffsl
c       write(6,*) 'Teff: ', Teffsl
        read (37,'(12x,'//sNsims//'f12.3)') gravsl
c       write(6,*) 'grav: ', gravsl
        read (37,'(12x,'//sNsims//'f12.3)') FeH
        read (37,'(12x,'//sNsims//'f12.3)') alph
        read (37,'(12x,'//sNsims//'f12.3)') sigalph
        read (37,'('//sNsims1//'f12.8)') aa1
        close(37)
c
c  Split off the log10(tau) array from aa1 and store in qltau
        do j=1, Nqtau
          qltau(j) = aa1(1,j)
          do i=1, Nsims
            aa(i,j) = aa1(i+1,j)
          enddo
        enddo
        do i=1, Nsims
          Teffsl(i) = dlog10(Teffsl(i))
        enddo
        if (Nsims .ge. 3) then
c
c  More than 3 simulations
c
c  Use LTeff3 as basis for triangulation since Teff varies a lot
c  less than g - this ensures less acute angles.
          do i=1, Nsims
            LTeff3(i) = (Teffsl(i)-3d0)*8d0
            Logg3 (i) =  gravsl(i)
          enddo
c  Remove these 11 lines if used with tables other than TtauFeH0.dat
c  from Trampedach et a. (2014).
          LTeff3(19) = LTeff3(19) - 0.15d0
          Logg3 (29) = Logg3 (29) + 0.08d0
          Logg3 ( 9) = Logg3 ( 9) + 0.3d0
          Logg3 ( 7) = Logg3 ( 7) - 0.2d0
          Logg3 (11) = Logg3 (11) + 0.2d0
          Logg3 (33) = Logg3 (33) - 0.07d0
          Logg3 (28) = Logg3 (28) - 0.4d0
          LTeff3(28) = LTeff3(28) - 0.02d0
          LTeff3(32) = LTeff3(32) - 0.1d0
          LTeff3(37) = LTeff3(37) + 0.03d0
          Logg3 (36) = Logg3 (36) - 0.15d0
c
c  Sort Teffsl, gravsl and alpha [=aa(*,1)] w.r.t. Teffsl and store
c  indexes in isrt
          call reordr (Nsims,3,Teffsl,gravsl,alph,isrt)
          call permut (Nsims,isrt,LTeff3)
          call permut (Nsims,isrt,Logg3)
c         print *,'After sorting by Teff:'
c         print *,Teffsl
c         print *,gravsl

c
c  - Set-up triangulation in the reordered log10(Teff) and log10(grav)
          call trmesh (Nsims,LTeff3,Logg3, iadj,iend,ier)
c         print *,'iadj: ', iadj
c         print *,'iend: ', iend
c         print *,'ier:  ', ier
c  and reorder the data array (as was done for the Teffsl, gravsl and alph
c  arrays in the above call to reordr).
          do j=1, Nqtau
            call permut (Nsims,isrt,aa(1,j))
          enddo
c         print *,(aa(i,1),i=1,Nsims)
        endif
        print *," Read the table, 'TtauFeH0.dat', of calibrated
     & mixing-length parameters"
        print *," and radiative Hopf functions (for T(tau)-relations)."
      endif
c
c  End of I/O-block.
c***********************************************************************
c  Start of interpolation.
c
c  Only results of one simulation in the table - keep it constant...
      if (Nsims .eq. 1) then
        alpha     = alph(1)
        do j=1, Nqtau
          qhpf(j) = aa(1,j)
        enddo
      endif
      if (Nsims .ge. 3) then
        ist = 1
c
c  Perform linear interpolation of alpha, in log10(Teff) and log10(grav)
        call   intrc0 (Nsims,Teffl,gravl,Teffsl,gravsl,alph,
     &                            iadj,iend,ist, alpha,   ier)
c
c  Perform linear interpolation of Hopf function, in log10(Teff) and log10(grav)
        do j=1, Nqtau
          call intrc0 (Nsims,Teffl,gravl,Teffsl,gravsl,aa(1,j),
     &                            iadj,iend,ist, qhpf(j), ier)
        enddo
      endif
c     print '("Hopf function: ")'
c     print *,qltau,qhpf
c
c  Flag to indicate that the derivative has not been evaluated yet.
      dqhpf(1) = -42.
c
      return
      end
