108 lines
2.9 KiB
Fortran
108 lines
2.9 KiB
Fortran
SUBROUTINE ZMRHO(R,HG)
|
|
C ======================
|
|
C
|
|
C Initial estimate of DM, DENS, and ZD
|
|
C by an approximate solution of the hydrostatic equilibrium
|
|
C equation
|
|
C Both gas and radiation pressure contribute;
|
|
C
|
|
C Input: R - ration of radiation to gas pressure scale heights
|
|
C HG - gas pressure scale height
|
|
C DM1 - mass at the first depth point
|
|
C DMTOT - mass at the last depth point (central plane)
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
PARAMETER (PISQ=1.77245385090551D0,pisq2=pisq*half)
|
|
C
|
|
C Mass-depth grid - logarithmicaly equidistant
|
|
C
|
|
if(nd.gt.mdepth) nd=mdepth
|
|
C
|
|
if(dm1.gt.0) then
|
|
DM(ND)=DMTOT
|
|
DM(ND-1)=0.99D0*DMTOT
|
|
DML=LOG(DM(ND-1)/DM1)/(ND-2)
|
|
DML1=LOG(DM1)
|
|
DO ID=1,ND-1
|
|
DM(ID)=EXP(DML1+(ID-1)*DML)
|
|
END DO
|
|
else if(dm1.lt.-1.e-20.and.dm1.gt.-1.e-10) then
|
|
DM(ND)=DMTOT
|
|
DM1=ABS(DM1*1.e20)
|
|
DML=LOG(DM(ND)/DM1)/(ND-1)
|
|
DML1=LOG(DM1)
|
|
DO ID=1,ND-1
|
|
DM(ID)=EXP(DML1+(ID-1)*DML)
|
|
end do
|
|
else if(dm1.gt.-1.e-10) then
|
|
if(mod(nd,2).eq.0) nd=nd-1
|
|
dmha=dmtot*half
|
|
dm1=abs(dm1*1.e10)
|
|
ndha=nd/2
|
|
dml=log(dmha/dm1)/ndha
|
|
dml1=log(dm1)
|
|
dm(nd)=dmtot
|
|
do id=1,ndha
|
|
dm(id)=exp(dml1+(id-1)*dml)
|
|
dm(nd-id)=dm(nd)-exp(dml1+id*dml)
|
|
end do
|
|
C
|
|
C Determination of the total pressure scale height - function BETAH
|
|
C HH - total pressure scale height
|
|
C
|
|
end if
|
|
HH=BETAH(R)*R
|
|
DMH=PISQ/2.D0*EXP(-R*(HH-R))*ERFCX(HH-R)/HH
|
|
RHO0=DM(ND)/HH/HG
|
|
C
|
|
C Approximate solution of the hydrostatic equilibrium
|
|
C
|
|
DO ID=1,ND
|
|
DMREL=DM(ID)/DM(ND)
|
|
IF(DMREL.LE.DMH) THEN
|
|
X=R+ERFCIN(DMREL*2.D0/PISQ*HH*EXP(R*(HH-R)))
|
|
RHO=EXP(-(X-R)*(X-R)-(HH-R)*R)
|
|
ELSE IF(DMREL.LT.UN) THEN
|
|
HSQ=SQRT(HH*(HH-R))
|
|
X=ERFCIN(2.D0/PISQ*HSQ*(DMREL-DMH)+ERFCX(HSQ))*HH/HSQ
|
|
RHO=EXP(-X*X*(UN-R/HH))
|
|
ELSE
|
|
X=0.
|
|
RHO=UN
|
|
END IF
|
|
DENS(ID)=RHO0*RHO
|
|
ZD(ID)=X*HG
|
|
END DO
|
|
c
|
|
if(dm1.lt.-1.e-10) then
|
|
DM(ND)=DMTOT
|
|
DM(ND-1)=0.99D0*DMTOT
|
|
dm1=abs(dm1)
|
|
DML=LOG(DM(ND-1)/DM1)/(ND-2)
|
|
DML1=LOG(DM1)
|
|
DO ID=1,ND-1
|
|
DM(ID)=EXP(DML1+(ID-1)*DML)
|
|
end do
|
|
c
|
|
hr=r*hg
|
|
hg2=hg*pisq2
|
|
hrg=hr+hg2
|
|
dmh=hg2/hrg
|
|
rho0=dmtot/hrg
|
|
do id=1,nd
|
|
dmrel=dm(id)/dm(nd)
|
|
if(dmrel.ge.dmh) then
|
|
zd(id)=hrg*(un-dmrel)
|
|
dens(id)=rho0
|
|
else
|
|
zd(id)=hr+hg*erfcin(dmrel*hrg/hg2)
|
|
x=(zd(id)-hr)/hg
|
|
dens(id)=rho0*exp(-x*x)
|
|
end if
|
|
end do
|
|
end if
|
|
RETURN
|
|
END
|