93 lines
2.4 KiB
Fortran
93 lines
2.4 KiB
Fortran
subroutine hedif
|
|
c ================
|
|
c
|
|
c subroutine to calculate the depth dependent abundance profile for
|
|
c a layered H+He atmosphere.
|
|
c
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
common/hediff/ hcmass,radstr
|
|
c real depth(mdepth+1),qs(mdepth+1),
|
|
dimension depth(mdepth+1),qs(mdepth+1),
|
|
* ps(mdepth+1),gams(mdepth+1),
|
|
* abunds(mdepth+1),hms(mdepth+1)
|
|
c
|
|
data smas,srad /1.9891e33,6.9599e10/
|
|
data z1,z2,a1,a2 /1.,2.,1.,4./
|
|
data bigg,pi / 6.6732e-8,3.141592654/
|
|
c
|
|
c Set up starting values
|
|
c
|
|
do id=1,nd
|
|
depth(id+1)=dm(id)
|
|
end do
|
|
if(radstr.lt.1.e3) radius=radstr*srad
|
|
if(hcmass.gt.1.e-10) hcmass=hcmass*1.e-13
|
|
c
|
|
gam=1.e-30
|
|
gams(1)=1.e-30
|
|
c
|
|
10 continue
|
|
depth(1)=1.e-10
|
|
q1=(depth(1)*4*pi*radius**2/smas)
|
|
p1=(q1*grav**2/(4*pi*bigg))
|
|
ps(1)=p1
|
|
qs(1)=q1
|
|
hms(1)=0
|
|
abunds(1)=0.0
|
|
dpsl=-6
|
|
hm=0.0
|
|
do i=2,nd+1,1
|
|
q2=(depth(i)*4*pi*radius**2/smas)
|
|
p2=(q2*grav**2/(4*pi*bigg))
|
|
dp=p2-p1
|
|
dlp=log(p2)-log(p1)
|
|
gam=gam+raph(gam,z1,z2,a1,a2)*dlp
|
|
abun0=gam
|
|
hm=hm+(q2-q1)/((1+gam*a2/a1))
|
|
p1=p2
|
|
ps(i)=p2
|
|
qs(i)=q2
|
|
gams(i)=gam
|
|
abunds(i)=abun0
|
|
hms(i)=hm
|
|
end do
|
|
c
|
|
dh1=(log10(hcmass)-log10(hms(nd+1)))
|
|
dh=hcmass/hms(nd+1)
|
|
if(dh.ge.0.99) go to 20
|
|
gam=gams(1)*1.1
|
|
gams(1)=gam
|
|
hm=0.0
|
|
go to 10
|
|
c
|
|
c Now work backwards to get the full profiles
|
|
c
|
|
20 continue
|
|
q1=(depth(nd+1)*4*pi*radius**2/smas)
|
|
p1=(q1*grav**2/(4*pi*bigg))
|
|
c
|
|
c store new helium abundance and corresponding new YTOT, MMY, WMM
|
|
c
|
|
write(6,600)
|
|
do id=1,nd
|
|
aheold=0.
|
|
ahenew=abunds(id+1)
|
|
if(iathe.gt.0) then
|
|
aheold=abund(iathe,id)
|
|
abund(iathe,id)=ahenew
|
|
end if
|
|
ytot(id)=ytot(id)-aheold+ahenew
|
|
wmy(id)=wmy(id)+(ahenew-aheold)*4.003
|
|
wmm(id)=wmy(id)*hmass/ytot(id)
|
|
write(6,601) id,aheold,ahenew,ytot(id),wmy(id),wmm(id)
|
|
end do
|
|
600 format(' stratified helium abundance'/
|
|
*' id He(old) He(new) ytot wmm wmy'/)
|
|
601 format(i4,1p5e11.3)
|
|
c
|
|
return
|
|
end
|