************************************************************************
*
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]
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