subroutine ledoux * * compute modified Ledoux B term to include composition gradient * effects in the Brunt-Vaisala frequency * Pressure is the dependent variable throught and I sum partial pressures. * In this version I change only the He abundance at the H/He interface. * while I keep the sum of the abundances unity at the He/C interface. * This duplicates the numerical differencing results, but the physical * justification eludes me at present. * PAB Sept 10, 1991 * implicit double precision(a-h,o-z) * * EOS data * common/tenv/rhok(3,60,35),eta(3,60,35),pp(3,60,35),datg(3,60,35), 1 od(3,18,35),tk(3,60),uu(3,60,35),xtt(3,60,35),xrt(3,60,35) * * composition info * common/je/je common/comp/amhyhe,amheca,x(4) * * need for lmax * common/d/ ml(3,60),nmax,kind,lmax,nsuff,nlim,nlim1,j,iter,lmam, 1 nlum,ko,nhp * * has log T=tl, log Rho=rho and log P=pl * also has stellar mass=sm and current mass=xmass * common/c/az(3),z(3),taux,ts,ps,sm,rs,err,te,gs,rho,ol1x,rl1x,el1x 1 ,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi * * diffusion profile info * common/alpha/alph(4),qmix1,qmix2,zone(2,2) * * ledoux term value * common/bldx/bled real*8 rhot1(35),rhot2(35),rhot3(35),rhot4(35) real*8 ptmp1(35),ptmp2(35),ptmp3(35),ptmp4(35) real*8 y21(35),y22(35),y23(35),y24(35) real*8 temp(60) * * initialize summation variables * bled = 0.0 * * step size increments in composition * del = 0.01 * * This routine is complicated because we have to consider the * H/He interface and He/C interfaces. * We also have to make sure to avoid non-physical compositions * such as X<0 or X>1. * * H/He interface * if(x(1).gt.0.0 .and. x(2).gt.0.0)then if(x(3) .gt. x(1))then * * Very thin He buffer layer. Call H/He or He/C transition zone * depending on case involved. * go to 45 endif * * first, locate T point on grid for Helium * do 992 i=1,lmax temp(i) = tk(2,i) 992 continue call hunt(temp,lmax,tl,ktlo) kthi = ktlo + 1 * * set size of temporary arrays * kmax1 = ml(2,ktlo) kmax2 = ml(2,kthi) * * fill up temporary 1-D arrays for rho and T * He arrays * do 10 k=1,kmax1 rhot1(k) = rhok(2,ktlo,k) ptmp1(k) = pp(2,ktlo,k) 10 continue do 20 k=1,kmax2 rhot2(k) = rhok(2,kthi,k) ptmp2(k) = pp(2,kthi,k) 20 continue * * first, locate T point on grid for Hydrogen * do 993 i=1,lmax temp(i) = tk(1,i) 993 continue call hunt(temp,lmax,tl,ktlo) kthi = ktlo + 1 * * set size of temporary arrays * kmax3 = ml(1,ktlo) kmax4 = ml(1,kthi) * * H arrays * do 30 k=1,kmax3 rhot3(k) = rhok(1,ktlo,k) ptmp3(k) = pp(1,ktlo,k) 30 continue do 40 k=1,kmax4 rhot4(k) = rhok(1,kthi,k) ptmp4(k) = pp(1,kthi,k) 40 continue * * set BC's so second der is zero there * yp1 = 1.0e+40 ypn = 1.0e+40 * * call for spline coefficients * call spline(rhot1,ptmp1,kmax1,yp1,ypn,y21) call spline(rhot2,ptmp2,kmax2,yp1,ypn,y22) call spline(rhot3,ptmp3,kmax3,yp1,ypn,y23) call spline(rhot4,ptmp4,kmax4,yp1,ypn,y24) * * now compute P for each isotherm at a given rho,comp point * call splint(rhot1,ptmp1,y21,kmax1,rho,p11) call splint(rhot2,ptmp2,y22,kmax2,rho,p12) call splint(rhot3,ptmp3,y23,kmax3,rho,p13) call splint(rhot4,ptmp4,y24,kmax4,rho,p14) * * linear interpolation for P at the actual temp. * frac=(tl-tk(je,ktlo))/(tk(je,kthi)-tk(je,ktlo)) pone=p11+frac*(p12-p11) ptwo=p13+frac*(p14-p13) pone=10.**(pone) ptwo=10.**(ptwo) * * now to consider specific composition mixtures! * CASE 1: x(1) or x(2) is less than del (use a 3 point formula) * CASE 1.1: Nearly pure helium: Treat first point as pure He. * treat second point as H / He mixture * if(x(1).lt.del)then diff = 1.0-x(2) xone = x(2) - diff xtwo = 1.0 * * compute P for mixture by partial pressures. * pmix1 is for x(2)-diff; pmix2 is for x(2)+diff * pmix1 = xone*pone + (diff)*ptwo pmix1 = dlog(pmix1) pmix2 = pone + (diff)*ptwo pmix2 = dlog(pmix2) * * compute d ln P/d ln Y and return * bled=(x(2)/(2.*diff))*(pmix2-pmix1) return * * CASE 1.2: Nearly pure hydrogen: Treat first point as pure H. * treat second point as H / He mixture * else if(x(2).lt.del)then diff=x(2) xone = 0.0 xtwo = x(2) + diff * * compute P for mixture by summing partial pressures. * pmix1 is for x(2)-diff; pmix2 is for x(2)+diff * pmix1 = x(1)*ptwo pmix1 = dlog(pmix1) pmix2 = xtwo*pone + x(1)*ptwo pmix2 = dlog(pmix2) bled=(x(2)/(2.*diff))*(pmix2-pmix1) return endif * * CASE 2: We are not near 0 or 1 in composition * we can use 3-point differencing formula with impunity. * xone = x(2) - del xtwo = x(2) + del pmix1 = xone*pone + x(1)*ptwo pmix2 = xtwo*pone + x(1)*ptwo pmix1 = dlog(pmix1) pmix2 = dlog(pmix2) * * compute d ln P/d ln Y and return * bled=(x(2)/(2.*del))*(pmix2-pmix1) return endif * * finished with H/He mixtures * * He/C interface * if(x(2).gt.0.0 .and. x(3).gt.0.0)then * * first, locate T point on grid for Helium * 45 do 994 i=1,lmax temp(i) = tk(2,i) 994 continue call hunt(temp,lmax,tl,ktlo) kthi = ktlo + 1 * * set size of temporary arrays * kmax1 = ml(2,ktlo) kmax2 = ml(2,kthi) * * fill up temporary 1-D arrays for rho and T * He arrays * do 50 k=1,kmax1 rhot1(k) = rhok(2,ktlo,k) ptmp1(k) = pp(2,ktlo,k) 50 continue do 60 k=1,kmax2 rhot2(k) = rhok(2,kthi,k) ptmp2(k) = pp(2,kthi,k) 60 continue * * first, locate T point on grid for Carbon * do 995 i=1,lmax temp(i) = tk(3,i) 995 continue call hunt(temp,lmax,tl,ktlo) kthi = ktlo + 1 * * set size of temporary arrays * kmax3 = ml(3,ktlo) kmax4 = ml(3,kthi) * * C arrays * do 70 k=1,kmax3 rhot3(k) = rhok(3,ktlo,k) ptmp3(k) = pp(3,ktlo,k) 70 continue do 80 k=1,kmax4 rhot4(k) = rhok(3,kthi,k) ptmp4(k) = pp(3,kthi,k) 80 continue * * set BC's so second der is zero there * yp1 = 1.0e+40 ypn = 1.0e+40 * * call for spline coefficients * call spline(rhot1,ptmp1,kmax1,yp1,ypn,y21) call spline(rhot2,ptmp2,kmax2,yp1,ypn,y22) call spline(rhot3,ptmp3,kmax3,yp1,ypn,y23) call spline(rhot4,ptmp4,kmax4,yp1,ypn,y24) * * now compute P for each isotherm at a given rho,comp point * call splint(rhot1,ptmp1,y21,kmax1,rho,p11) call splint(rhot2,ptmp2,y22,kmax2,rho,p12) call splint(rhot3,ptmp3,y23,kmax3,rho,p13) call splint(rhot4,ptmp4,y24,kmax4,rho,p14) * * linear interpolation for P at the actual temp. * frac=(tl-tk(je,ktlo))/(tk(je,kthi)-tk(je,ktlo)) pone=p11+frac*(p12-p11) ptwo=p13+frac*(p14-p13) pone=10.**(pone) ptwo=10.**(ptwo) * * now to consider specific composition mixtures! * CASE 3: x(3) or x(2) is less than del (use a 3 point formula) * CASE 3.1: Nearly pure helium: Treat first point as pure He. * treat second point as He / C mixture * if(x(3).lt.del)then diff = 1.0-x(2) xone = x(2) - diff xtwo = 1.0 * * compute P for mixture by summing partial pressures. * pmix1 is for x(2)-diff; pmix2 is for x(2)+diff * pmix1 = xone*pone + (1.0-xone)*ptwo pmix2 = xtwo*pone + (1.0-xtwo)*ptwo pmix1 = dlog(pmix1) pmix2 = dlog(pmix2) bled=(x(2)/(2.*diff))*(pmix1-pmix2) return * * CASE 3.2: Nearly pure Carbon: Treat first point as pure C. * treat second point as He / C mixture * else if(x(2).lt.del)then diff=x(2) xone = 0.0 xtwo = x(2) + diff * * compute P for mixture by additive volume method. * pmix1 is for x(2)-diff; pmix2 is for x(2)+diff * pmix1 = xone*pone + (1.0-xone)*ptwo pmix2 = xtwo*pone + (1.0-xtwo)*ptwo pmix1 = dlog(pmix1) pmix2 = dlog(pmix2) * * compute d ln P/d ln Y and return. * bled=(x(2)/(2.*diff))*(pmix1-pmix2) return endif * * CASE 2: We are not near 0 or 1 in composition * we can use 3-point differencing formula with impunity. * xone = x(2) - del xtwo = x(2) + del pmix1 = xone*pone + (1.0-xone)*ptwo pmix2 = xtwo*pone + (1.0-xtwo)*ptwo pmix1 = dlog(pmix1) pmix2 = dlog(pmix2) * * compute d ln P/d ln Y and return. * bled=(x(2)/(2.*del))*(pmix1-pmix2) return endif * * finished with He/C mixtures. * * If no transition zone, then set bled to 0.0 * bled = 0.0 return end ************************************************************************