190 lines
6.2 KiB
Fortran
190 lines
6.2 KiB
Fortran
SUBROUTINE INPDIS
|
|
C =================
|
|
C
|
|
C driver for input specific for disks
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
INCLUDE 'ODFPAR.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
PARAMETER (VELC=2.997925E10,
|
|
* pi4=12.5663706d0)
|
|
PARAMETER (GRCON = 6.668D-8)
|
|
common/relcor/arh,brh,crh,drh
|
|
C
|
|
C ----------------------
|
|
C Basic input parameters
|
|
C ----------------------
|
|
C
|
|
C The user may choose one of the two
|
|
C following possibilities to input the basic physical parameters:
|
|
C
|
|
C XMSTAR - M(star), either in M(Sun), or in grams;
|
|
C XMDOT - M(dot), either in M(Sun)/year; or in g/s
|
|
C RSTAR - R(star), either in R(Sun), or in cm
|
|
C RELDST - R/R(star)
|
|
C
|
|
C or to directly give parameters with which the program works, ie
|
|
C
|
|
C TEFF - effective temperature
|
|
C QGRAV - coefficient in the hydrostatic equilibrium equation
|
|
C QGRAV=G*M(star)/R**3
|
|
C DMTOT - total column mass at the midplane
|
|
C
|
|
WRITE(6,660)
|
|
IF(FRACTV.LT.0.) THEN
|
|
AMUV0=DMVISC**(ZETA0+UN)
|
|
FRACTV=UN/(UN+(ZETA0+UN)/(ZETA1+UN)*AMUV0/(UN-AMUV0))
|
|
END IF
|
|
IF(DMVISC.LT.0.) DMVISC=(UN/(UN+(ZETA0+UN)/(ZETA1+UN)*
|
|
* FRACTV/(UN-FRACTV)))**(UN/(ZETA0+UN))
|
|
alpha0=alphav
|
|
WRITE(6,600) XMSTAR,XMDOT,RSTAR,RELDST,ALPHAV
|
|
C
|
|
c if XMSTAR<0, turn on general relativistic corrections
|
|
c RSTAR now has the input meaning of dimensionless angular
|
|
c momentum of the Kerr black hole, but we will call this
|
|
c value AA and RSTAR will take on the meaning of 1 radius-
|
|
c equivalent of mass of the black hole
|
|
c RELDST now is expressed in multiples of 1*G*XMSTAR/c^2
|
|
c (note that for Schwarzschild black hole, the
|
|
c horizon is at RELDST=2 and the smallest radius of
|
|
c stable circular orbit is at RELDST=6)
|
|
C
|
|
IF(XMSTAR.NE.0.) THEN
|
|
IF(XMSTAR.LT.0) THEN
|
|
AA=RSTAR
|
|
RSTAR=-XMSTAR*1.989D33*GRCON/VELC/VELC
|
|
ELSE
|
|
AA=0.
|
|
END IF
|
|
IF(abs(XMSTAR).GT.1.D16) XMSTAR=XMSTAR/1.989D33
|
|
IF(XMDOT.GT.1.D3) XMDOT=XMDOT/6.3029D25
|
|
IF(RSTAR.GT.1.D3) RSTAR=RSTAR/6.9598D10
|
|
R=RSTAR*ABS(RELDST)
|
|
QGRAV=5.9D0*GRCON*abs(XMSTAR)/R**3
|
|
OMEG32=SQRT(QGRAV)*1.5
|
|
c
|
|
c apply general relativistic corrections to
|
|
c QGRAV and TEFF; keep MSTAR<0 for future use
|
|
c
|
|
RR0=RELDST
|
|
CALL GRCOR(AA,RR0,XMSTAR,QCOR,TCOR,ARH,BRH,CRH,DRH)
|
|
TEFF0=(1.79049311D-1*QGRAV*3.34379D24*XMDOT/SIG4P)**0.25
|
|
TEFF=TEFF0*TCOR
|
|
QGRAV=QGRAV*QCOR
|
|
OMEG32=OMEG32*ARH/BRH
|
|
rr0=abs(rr0)
|
|
xmas9=abs(Xmstar)*1.d-9
|
|
xmdt=Xmdot/xmas9/2.22
|
|
XMD=XMDOT*6.3029D25
|
|
alpav=abs(alphav)
|
|
c
|
|
c following is the old procedure
|
|
c
|
|
if(alphav.le.0.) then
|
|
C -------------------------
|
|
chih=0.39
|
|
if(reynum.le.0.) then
|
|
reynum=(rr0/xmdt)**2/alpav*arh*crh/drh/drh
|
|
c REYNUM=(R/XMDOT*1.10422E-15*12.5663*VELC/CHIH)**2*
|
|
c * 2./ALPAV*ARH*CRH/DRH/DRH
|
|
else
|
|
alpav=(rr0/xmdt)**2/reynum*arh*crh/drh/drh
|
|
c ALPAV=(R/XMDOT*1.10422E-15*12.5663*VELC/CHIH)**2*
|
|
c * 2./REYNUM*ARH*CRH/DRH/DRH
|
|
endif
|
|
VISC=1.176565D22*SQRT(GRCON*abs(XMSTAR)*R)/REYNUM
|
|
DMTOT=3.34379D24*XMDOT/VISC*BRH*DRH/ARH/ARH
|
|
C
|
|
if(alphav.lt.0.) then
|
|
c
|
|
C ****************************************************************
|
|
C Compute the Keplerian rotation frequency (omega=c/r_g/x^1.5):
|
|
C
|
|
RE=ABS(RR0)
|
|
OMEGA=VELC/RSTAR/6.9698D10/RE**1.5D0
|
|
C
|
|
C Compute Relativistic factors, using Krolik (1998) notation:
|
|
C
|
|
RELT=DRH/ARH
|
|
RELR=DRH/BRH
|
|
RL2=RE*(1.D0-2.D0*AA/RE**1.5D0+AA*AA/RE/RE)**2/BRH
|
|
EINF=(1.D0-2.D0/RE+AA/RE**1.5D0)/SQRT(BRH)
|
|
RELZ=(RL2-AA*AA*(EINF-1.D0))/RE
|
|
C
|
|
C Compute the surface mass density (assuming pure electron scattering,
|
|
C pure hydrogen composition, and that T_rphi = \alpha P_tot):
|
|
C
|
|
DMTOT=0.5D0*SIGMAR(ALPAV,XMD,TEFF,OMEGA,RELR,RELT,RELZ)
|
|
end if
|
|
c ----------------------------
|
|
c
|
|
c new procedure
|
|
c
|
|
else
|
|
call column
|
|
end if
|
|
c
|
|
EDISC=SIG4P*TEFF**4/DMTOT
|
|
WBARM=XMDOT*6.3029D25/6./3.1415926*BRH*DRH/(ARH*ARH)
|
|
reynum=dmtot/wbarm*sqrt(xmstar*r)*3.03818e18
|
|
WRITE(6,601) TEFF,QGRAV,DMTOT,
|
|
* ZETA0,ZETA1,FRACTV,DMVISC,TSTAR
|
|
write(6,321) tcor,qcor,arh,brh,crh,drh,
|
|
* reynum,alpav,
|
|
* dmtot,edisc,wbarm
|
|
321 FORMAT(
|
|
* ' tcor =',1PD10.3/
|
|
* ' qcor =',D10.3/
|
|
* ' A(RH) =',D10.3/
|
|
* ' B(RB) =',D10.3/
|
|
* ' C(RH) =',D10.3/
|
|
* ' D(RH) =',D10.3/
|
|
* ' Re =',D10.3/
|
|
* ' alpha =',D10.3//
|
|
* ' DMTOT =',D10.3/
|
|
* ' EDISC =',D10.3/
|
|
* ' WBARM =',D10.3//)
|
|
C
|
|
ELSE
|
|
TEFF=XMDOT
|
|
QGRAV=RSTAR
|
|
DMTOT=RELDST
|
|
EDISC=SIG4P*TEFF**4/DMTOT
|
|
OMEG32=SQRT(QGRAV)*1.5
|
|
WRITE(6,601) TEFF,QGRAV,DMTOT,
|
|
* ZETA0,ZETA1,FRACTV,DMVISC,TSTAR
|
|
END IF
|
|
C
|
|
C set up the maximum frequency
|
|
C
|
|
c if(idgrey.le.2) then
|
|
c IF(FRCMAX.EQ.0.) FRCMAX=2.83e11*(dmtot*0.39)**0.25*teff
|
|
c end if
|
|
IF(FRLMAX.EQ.0.) FRLMAX=1.D11*CNU1*TEFF
|
|
C
|
|
660 FORMAT(1H1,'***************************************'//
|
|
* ' M O D E L O F A D I S K R I N G'//
|
|
* ' ***************************************'//)
|
|
600 FORMAT(
|
|
* ' M(STAR) =',1PD10.3/
|
|
* ' M(DOT) =',D10.3/
|
|
* ' R(STAR) =',D10.3/
|
|
* ' R/R(STAR) =',D10.3/
|
|
* ' ALPHA =',D10.3//)
|
|
601 FORMAT(
|
|
* ' TEFF =',F10.0/
|
|
* ' QGRAV =',1PD10.3/
|
|
* ' DMTOT =',D10.3/
|
|
* ' ZETA0 =',D10.3/
|
|
* ' ZETA1 =',D10.3/
|
|
* ' FRACTV =',D10.3/
|
|
* ' DMVISC =',D10.3/
|
|
* ' TSTAR =',0PF10.0//)
|
|
RETURN
|
|
END
|