subroutine ccomp * * this routine returns the h, he, and c mass fractions given as * input: xmass (position in the envelope), amhyhe, amheca. * * this is the first diffusion version. the profiles returned are * derived under several assumptions (some questionable). in * short, we can do better, but this is going to be it for now. * we're doing a pulsation study, so it's mostly important that * the profiles be smooth. this version will not mix h into he * by convection * * Modified by PAB to remove discontinuity in He abundance when * helium layer is only 100x thicker than hydrogen layer. 09 05 91 * Here, I let both H and C be trace elements between zone(1,2) * and zone(2,1). 09 09 91 (This one works in giving a smooth profile.) * implicit double precision(a-h,o-z) common/je/je common/comp/amhyhe,amheca,x(4) common/c/az(3),z(3),taux,ts,ps,sm,rs,err,te,gs,rho, 1 ol1x,rl1x,el1x,atgx,rtgx,tl,pl,xmass,acu,am,bk,an0,pi common/alpha/alph(4),qmix1,qmix2,zone(2,2) common/idiffus/idiffus ame = 1.e0-(xmass/sm) * * determine surface composition * if(amhyhe.lt.1.0e-15.and.amheca.lt.1.0e-15) je=3 if(amhyhe.lt.1.0e-15.and.amheca.gt.1.0e-15) je=2 if(amhyhe.gt.1.0e-15) je=1 * * if it's a diffusion equilibrium, then proceed, otherwise, jump * down to the discontinuity code * * determine local composition and load comp vector x accordingly *** hydrogen surface layer *** if ( idiffus .eq. 1 ) then if ( je .eq. 1 ) then if(zone(1,2) .lt. zone(2,1))then * * Normal case. We compute profiles without a care in the world. * pure h * if ( ame .le. zone(1,1) ) then x(1) = 1.0 x(2) = 0.0 x(3) = 0.0 x(4) = 0.0 * * lower h layer (includes he tail) * elseif ( ame .gt. zone(1,1) .and. ame .le. qmix1 ) then xx = ame/qmix1 yy = func1(xx) x(1) = yy x(2) = 1.0 - yy x(3) = 0.0 x(4) = 0.0 * * upper he layer * elseif ( ame .gt. qmix1 .and. ame .le. zone(1,2) ) then xx = ame/qmix1 yy = func2(xx) x(1) = yy x(2) = 1.0 - yy x(3) = 0.0 x(4) = 0.0 * * pure he * elseif ( ame .gt. zone(1,2) .and. ame .le. zone(2,1) ) then x(1) = 0. x(2) = 1.0 x(3) = 0.0 x(4) = 0.0 * * lower he layer * elseif ( ame .gt. zone(2,1) .and. ame .le. qmix2) then xx = ame/qmix2 yy = func3(xx) x(1) = 0. x(2) = yy x(3) = 1.0 - yy x(4) = 0.0 * * upper c layer * elseif ( ame .gt. qmix2 .and. ame .le. zone(2,2) ) then xx = ame/qmix2 yy = func4(xx) x(1) = 0. x(2) = yy x(3) = 1.0 - yy x(4) = 0.0 * * pure c * elseif ( ame .gt. zone(2,2) ) then x(1) = 0.0 x(2) = 0.0 x(3) = 1.0 x(4) = 0.0 endif else * * Special case where M_H is only 100x less than M_He. * We compute profiles very carefully to avoid discontinuous slopes. * * compute changeover point from H as trace to C as trace * zond1 = dlog10(zone(1,2)) zond2 = dlog10(zone(2,1)) zond3 = (zond1 + zond2) / 2. change = 10.**(zond3) * * pure h (same as before) * if ( ame .le. zone(1,1) ) then x(1) = 1.0 x(2) = 0.0 x(3) = 0.0 x(4) = 0.0 * * lower h layer (includes he tail) (again, same as before) * elseif ( ame .gt. zone(1,1) .and. ame .le. qmix1 ) then xx = ame/qmix1 yy = func1(xx) x(1) = yy x(2) = 1.0 - yy x(3) = 0.0 x(4) = 0.0 * * upper he layer (up to point where C should become a trace element) * elseif ( ame .gt. qmix1 .and. ame .le. zone(2,1) ) then xx = ame/qmix1 yy = func2(xx) x(1) = yy x(2) = 1.0 - yy x(3) = 0.0 x(4) = 0.0 * * Middle of He layer (between zone(1,2) and Zone (2,1)) * elseif ( ame .gt. zone(2,1) .and. ame .le. zone(1,2) ) then xx = ame/qmix1 yy = func2(xx) x(1) = yy yhe = 1.0 - yy zz = ame/qmix2 ww = func3(zz) if (zz.le.1.0) then ww = func3(zz) else ww = func4(zz) endif * * changing following lines so that there is no negative helium abundance * x(2) = ww - yy x(3) = 1.0 - ww x(4) = 0.0 * * lower He layer. (Below where H is a trace element) * elseif ( ame .gt. zone(1,2) .and. ame .le. qmix2) then xx = ame/qmix2 yy = func3(xx) x(1) = 0. x(2) = yy x(3) = 1.0 - yy * * upper c layer (same as before) * elseif ( ame .gt. qmix2 .and. ame .le. zone(2,2) ) then xx = ame/qmix2 yy = func4(xx) x(1) = 0. x(2) = yy x(3) = 1.0 - yy x(4) = 0.0 * * pure c (same as before): * elseif ( ame .gt. zone(2,2) ) then x(1) = 0.0 x(2) = 0.0 x(3) = 1.0 x(4) = 0.0 endif endif *** helium surface layer *** elseif ( je .eq. 2 ) then * * pure he * if ( ame .le. zone(2,1) ) then x(1) = 0. x(2) = 1.0 x(3) = 0.0 x(4) = 0.0 * * lower he layer * elseif ( ame .gt. zone(2,1) .and. ame .le. qmix2) then xx = ame/qmix2 yy = func3(xx) x(1) = 0. x(2) = yy x(3) = 1.0 - yy x(4) = 0.0 * * upper c layer * elseif ( ame .gt. qmix2 .and. ame .le. zone(2,2) ) then xx = ame/qmix2 yy = func4(xx) x(1) = 0. x(2) = yy x(3) = 1.0 - yy x(4) = 0.0 * * pure c * elseif ( ame .gt. zone(2,2) ) then x(1) = 0.0 x(2) = 0.0 x(3) = 1.0 x(4) = 0.0 endif *** carbon surface layer *** elseif ( je .eq. 3 ) then x(1) = 0.0 x(2) = 0.0 x(3) = 1.0 x(4) = 0.0 endif else * * come here for discontinuity profiles * * determine local composition * if(je.ne.1) goto 100 if(ame.gt.amhyhe.and.amheca.gt.1.0e-15) je=2 if(ame.gt.amhyhe.and.amheca.lt.1.0e-15) je=3 100 continue if(je.ne.2) goto 101 if(ame.gt.amheca) je=3 101 continue * * load comp vector X accordingly * if ( je .eq. 1 ) then x(1) = 1.0 x(2) = 0.0 x(3) = 0.0 x(4) = 0.0 else if ( je .eq. 2 ) then x(1) = 0.0 x(2) = 1.0 x(3) = 0.0 x(4) = 0.0 else if ( je .eq. 3 ) then x(1) = 0.0 x(2) = 0.0 x(3) = 1.0 x(4) = 0.0 endif endif return end ************************************************************************