The solar diffusion routine by Thoul et al. is available by 
anonymous ftp at the following address:
         ftp ftp.sns.ias.edu
         cd pub
	 cd thoul
         get routine.f
The routine is based on the calculations developed in the following paper:
    Thoul, A. A., Bahcall, J. N., and Loeb, A. 1994, Ap.J. 421, 828.

Here are a few additional explanatory notes about the diffusion routine.

Questions and/or comments are welcome. Please send them to
           anne.thoul@ulg.ac.be

-------------------------------------------------------------------------
-------------------------------------------------------------------------

Given M elements of atomic mass numbers A(M), charges Z(M), 
and mass fractions X(M), and given the Coulomb logarithms
CL(M,M), the routine returns the diffusion coefficients AP(M), 
AT(M), and AX(M,M).

(NOTE: in the routine, the maximum number of elements MMAX has been 
set to 20. If the number of elements is larger, this number
should be increased. And NMAX=2*MMAX+2 should also be changed.)

The diffusion function $\xi$, as defined in TBL (Thoul \& al, Ap.J. 
421, p.828, 1994), can then be obtained by multiplying these coefficients 
by the appropriate gradients:
\xi(M)=AP(M) {d lnp/dr} + AT(M) {d lnT/dr} + \sum_N AX(M,N) {d lnC(N)/dr}
where N is different from 2 (Helium) and M (electrons).

Note that these gradients are dimensionless. The radii are in units of 
R_\sun, as defined in TBL.

-------------------------------------------------------------------------
If mass fraction gradients are available instead of concentration fraction
gradients, it is necessary to convert the latter.
To get (dimensionless) concentration gradients from (dimensionless) mass 
fraction gradients: 
(here, GRADX contains the mass fraction gradients dlnX(N)/dr, and GRADC 
contains the concentration fraction gradients dlnC(N)/dr)
      TEMP=0.
      DO I=1,M-1
         TEMP=TEMP+Z(I)*X(I)/A(I)*GRADX(I)
      ENDDO
      DO J=1,M-1
         GRADC(J)=TEMP
      ENDDO
      TEMP=0.
      DO I=1,M-1
         TEMP=TEMP+Z(I)*X(I)/A(I)
      ENDDO
      DO J=1,M-1
         GRADC(J)=GRADX(J)-GRADC(J)/TEMP
      ENDDO

Alternatively, using the formula given in footnote #2 in TBL, it is possible to
express the diffusion function $\xi$ in terms of the mass fraction gradients
rather than the concentration gradients (NOTE: again, these are dimensionless 
gradients):
\xi(M)=AP(M) {d lnp/dr} + AT(M) {d lnT/dr} + \sum_N AX(M,N) {d lnX(N)/dr}
     -\sum_N AX(M,N))*(\sum_J Z(J)X(J)/A(J) {d lnX(J)/dr})/(\sum_I Z(I)X(I)/A(I))
where N is different from 2 (Helium) and M (electrons), and the sums over I and J
are over all ionic species. 
      TEMP1=0.
      TEMP2=0.
      DO I=1,M-1
	 TEMP1=TEMP1+Z(I)*X(I)/A(I)
         TEMP2=TEMP2+Z(I)*X(I)/A(I)*GRADX(I)
      ENDDO
      DO I=1,M
         SUM1=0.
         SUM2=0.
         DO J=1,M-1
            IF (J.NE.2) THEN
               SUM1=SUM1+AX(I,J)*GRADX(J)
	       SUM2=SUM2+AX(I,J)
            ENDIF
         ENDDO
         XI(I)=AP(I)*GRADP + AT(I)*GRADT + SUM1 - SUM2*TEMP2/TEMP1
      ENDDO
            


To obtain the dimensionless diffusion velocity, simply multiply \xi by 
T^{5/2}/\rho where T is in units of 10^7 K, and \rho is in units of 100 g/cm^3.
To obtain the velocity in cgs units, simply multiply the dimensionless
velocity by R_\sun/\tau_0 (see TBL).

-------------------------------------------------------------------------
To obtain the values of the Coulomb logarithms used in TBL,
add the following lines just before calling the diffusion subroutine 
in your main program. 
NOTES: (1) in the following, RHO and T are the mass density and the 
           temperature in CGS UNITS!
       (2) there is a typo in TBL paper: in eq. 9, defining the 
           Coulomb logarithms, the quantity in the parenthesis 
           (4kT\lambda/Z_s Z_t e^2) should be raised to the power 1.2.

c additional variables to be declared:
c     intermediate variables:
      REAL ZXA,AC,NI,CZ,XIJ,NE,AO,LAMBDAD,LAMBDA,C(M)

c calculate concentrations from mass fractions:
      ZXA=0.
      DO I=1,M-1
	 ZXA=ZXA+Z(I)*X(I)/A(I)
      ENDDO
      DO I=1,M-1
         C(I)=X(I)/(A(I)*ZXA)
      ENDDO
      C(M)=1.
c calculate density of electrons (NE) from mass density (RHO):
      AC=0.
      DO I=1,M
	 AC=AC+A(I)*C(I)
      ENDDO	
      NE=RHO/(1.6726E-24*AC)	
c calculate interionic distance (AO): 
      NI=0.
      DO I=1,M-1
         NI=NI+C(I)*NE
      ENDDO
      AO=(0.23873/NI)**(1./3.)	
c calculate Debye length (LAMBDAD):
      CZ=0.
      DO I=1,M
	 CZ=CZ+C(I)*Z(I)**2
      ENDDO
      LAMBDAD=6.9010*SQRT(T/(NE*CZ))
c calculate LAMBDA to use in Coulomb logarithm:
      LAMBDA=MAX(LAMBDAD,AO)
c calculate Coulomb logarithms:
      DO I=1,M
	 DO J=1,M
	    XIJ=2.3939E3*T*LAMBDA/ABS(Z(I)*Z(J))
	    CL(I,J)=0.81245*LOG(1.+0.18769*XIJ**1.2)
	 ENDDO
      ENDDO
 
---------------------------------------------------------------------