accdisk.f File Reference

Go to the source code of this file.

Functions/Subroutines
program	at
subroutine	rdprm
subroutine	intmdl
subroutine	sadmk (R, ETA, DL)
subroutine	sadmk2 (R, ETA, DL, DE, VX, VP, TE)
subroutine	sadmn (R, ETA, DL)
subroutine	sadmn2 (R, ETA, DL, DE, VX, VP, TE)
subroutine	intgrt (R, ETA, DL, IC)
subroutine	fun (XI, YI, YP)
subroutine	jac (xi, yi, yp, dfdr, pdj)
double precision function	temp (DE, PR)
double precision function	ondo (DE, PR, taup, tau)
double precision function	fnin (NPL)
double precision function	fnjn (NPL)
subroutine	cnvrbs
subroutine	lw1d
subroutine	cfl
subroutine	prtrb1
subroutine	prtrb2
subroutine	prtrb3
subroutine	prtrb4
subroutine	print
subroutine	file
subroutine	rdfl7
subroutine	stifbs (y, dydx, nv, x, htry, eps, yscal, hdid, hnext, derivs, icon)
subroutine	simpr (y, dydx, dfdx, dfdy, nmax, n, xs, htot, nstep, yout, derivs)
subroutine	pzextr (iest, xest, yest, yz, dy, nv)
subroutine	ludcmp (a, n, np, indx, d)
subroutine	lubksb (a, n, np, indx, b)

Function/Subroutine Documentation

program at

(

)

Definition at line 25 of file accdisk.f.

Here is the call graph for this function:

subroutine cfl

(

)

Definition at line 2406 of file accdisk.f.

C**************************
C
      include 'para.h'
C
      dimension drp(jmax),drm(jmax),dr(jmax)
      dimension dtj(jmax)
      dimension u(jmax,4),f(jmax,4),h(jmax,4),y(jmax,4)
C
c      print *,"CFL"
C
      gm1 = gm-1
      cc2 = cc*cc
     
      DO 5 j=2,jx-1
         drp(j) = x1(j+1)-x1(j)
         drm(j) = x1(j)-x1(j-1)
         dr(j)  = 0.5*(drp(j)+drm(j))
 5    CONTINUE
      drp(1) = x1(2) -x1(1)
      drm(jx)= x1(jx)-x1(jx-1)
      dr(1)  = drp(1)
      dr(jx) = drm(jx)
      taueffjx=sqrt(tauabso(jx)*(tauo(jx)+2./sqrt(3.)))
     
      ji=1
      DO 10 j=ji,jx
         r    = x1(j)
         rm1  = r-1
         r2   = r*r
         r3   = r2*r
         ok   = (r/rm1)*sqrt(0.5/r3)
         de   = w0(j,1)
         vx   = w0(j,2)
         vp   = w0(j,3)
         te   = w0(j,4)
         if (j.eq.jx) then
            pr   = (aa/3.)*(te**4)+rmu*de*te
         else
            if (kcool.eq.1) then 
               fac  =  1./((1.+1./tauabso(j)/(1.5*tauo(j)+sqrt(3.))))
            else 
               fac  = 1.0
            endif
            pr   = (aa/3.)*(te**4)*fac+rmu*de*te
         endif
         beta = rmu*de*te/pr
         hh   = sqrt(b5*pr/de)/ok
         ss   = 2*ain*de*hh*rg
         ww   = 2*ainp1*pr*hh*rg
         wc   = ww/cc2
         ws   = wc/ss
         a1   = 4.-3.*beta
         a2   = beta+12.*gm1*(1.-beta)
         gg1  = beta+(a1*a1*gm1)/a2
         gg3  = 1.+(gg1-beta)/a1
         cs2  = ((3.*gg1-1.)+2.*(gg3-1.)*alpha*alpha)*ws/(gg1+1.)
         cs   = sqrt(cs2)
         dtj(j) = dr(j)/(cs+abs(vx))
         if (nl(j).ge.2) dtj(j)=dtmax
c     
         demn   = ain*de
         temn   = 2*te/3.0
         xkff   = xkkr*demn*(temn**(-3.5))
         tauff  = xkff*de*hh*rg
         taues  = xkes*de*hh*rg
         tau    = taues+tauff

         if (kcool .eq. 1) then 
            taueff = sqrt(tauabso(j)*(tauo(j)+2./sqrt(3.)))
         else 
            taueff = sqrt(3*tau*tauabso(j))
         end if
     
         if (taueff.gt.taueffjx) dtj(j)=dtmax
 10   CONTINUE

      dt = dtmax
      DO 20 j = ji,jx
c     DT = MIN(DT,.04*DTJ(J))
         dt = min(dt,cfc*dtj(j))
 20   CONTINUE
c     
      dt = max(dt,dtmin)
c     
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine cnvrbs

(

)

Definition at line 1495 of file accdisk.f.

C**************************
C
C -----------------------------------------------------------
C     TRANSFORMATION FROM (Y0:ETA,DL) TO (W0:DE,VR,DL,TE)
C     SET OUTER BOUNDARY AT ROUT2
C -----------------------------------------------------------
C
      include 'para.h'
C
c      write(6,*) "sub-CNVRBS"
      j=0
C
 200                         continue
      j  = j+1
      jj = jx-j+1
      r  = x0(jj)
                             if (r.le.rin2) goto 200
100   CONTINUE
      j3 = jj
      j  = 0
C
 2000                        continue
C
      j   = j+1
      jj  = j3-j+1
      r   = x0(jj)
      eta = y0(1,jj)
      IF (eta.LT.0.) THEN
         kpr=1
         return
      ENDIF
      dl  = y0(2,jj)
      rm1 = r-1.
      r2  = r*r
      r3  = r2*r
      ok  = (r/rm1)*sqrt(.5/r3)
      xlk = r2*ok
      xl  = xlk+dl
      xll = xl-xlin
      xll2= xll*xll
      al2 = alpha*alpha
      ws  = xll2*eta/(al2*r2)
      pd  = (ain/ainp1)*ws
      hh  = (sqrt(2.*xnpl1*pd))/ok
      pr  = (xmdot*xll/(4.*pi*r2*alpha*ainp1*hh))*cc/(rg*rg)
      if (pr.le.0.) then
          kpr=1
          return
      endif
      de  = pr/(pd*cc*cc)
      tauff  = xkff*de*hh*rg
      taues  = xkes*de*hh*rg
      tau    = taues+tauff
      te = temp(de,pr)

      qbr   = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
      taup  = qbr/(8.*aa*cc*(te**4))

c      write(6,*) r,tau,taup
c      pause

      if ((kr.eq.1).or.(r.eq.rout1)) then
c          te = temp(de,pr)
c         if (kcool.eq.0) TE  = TEMP(DE,PR)
         if (kcool.eq.1) te  = ondo(de,pr,taup,tau)
      else
         te  = pr/de/rmu
      endif
      s0  = (4.*aa/(3.*de))*(te**3)+rmu*((1./gm1)*log(te)-log(de))
      beta= rmu*de*te/pr
      a   = 3*(1-beta)+beta/gm1
      ss  = 2*ain*de*hh*rg
      ww  = 2*ainp1*hh*pr*rg
      wc  = ww/(cc*cc)
      ws  = wc/ss
      vx  =-xmdot/(2.*pi*r*ss*cc*rg)
      vp  = dl/r
      ee  = (a+0.5)*ws+0.5*(vx*vx+vp*vp)
C
      w0(j,1) = de
      w0(j,2) = vx
      w0(j,3) = vp
      w0(j,4) = te
      x1(j)   = r
      rs      = r*ss
      u0(j,1) = rs
      u0(j,2) = rs*vx
      u0(j,3) = rs*vp
      u0(j,4) = rs*ee
C
      demn   = ain*de
      temn   = 2*te/3.0
      xkff   = xkkr*demn*(temn**(-3.5))
c      TAUES  = 2.*XKES*DE*HH*RG
      taues  = xkes*de*hh*rg
c      TAUFF  = 2.*XKFF*DE*HH*RG
      tauff  = xkff*de*hh*rg
      tau    = taues+tauff
      qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
      taup   = qbr/(8.*aa*cc*(te**4))
      tauabso(j)= taup
      tauo(j)=tau
C
 1    format ("w0=",4e11.3)
 2    format ("x1=",f7.3," rs=",e11.3)
 3    format ("u0=",4e11.3)
C 
                            if (r.lt.rout2) goto 2000
1000  CONTINUE
C
      jx = j
      if (imore.eq.2) 
     &   WRITE (6,*)  '    ### MESH POINTS = ',jx,' ###'
C
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine file

(

)

Definition at line 2833 of file accdisk.f.

C**************************
C
      include 'para.h'
C
      if (kf.eq.0) then
         write (8,902) jx
         write (8,903)
         do 10 j=1,jx
            write (8,904) j,x1(j)
 10      continue   
         kf=1
      ENDIF
C     
      IF ((mod(ns,ksfl).NE.0).AND.(ns.NE.nsx)) RETURN
C     IF ((ns.eq.lns).and.(ns.ne.0)) return
C     
c     WRITE (2) NS,TIME,DT,((W0(J,K),K=1,4),J=1,JX)
      write (8,*)
      write (8,905) ns,time
      write (8,906)
      do 20 j=1,jx
         write (8,907) j,(w0(j,k),k=1,4),tauabso(j),tauo(j),x1(j)
c     write (8,907) j,(w0(j,k),k=1,4),TAUABSO(J),x1(J)
 20   continue
      write (8,908) dt
      write (8,909) xmd
      write (8,901) xlin
C
      RETURN
c     900  format (80x)
 901  format ("        XLIN=",f18.14)
 902  format ("  JX=",i4)
 903  format (2x,"J",10x,"R")
 904  format (i4,1pe15.7)
 905  format (" NS=",i12,"   TIME=",f12.3)
 906  format (5x,"J",11x,"W0(J,1)",15x,"W0(J,2)",15x,"W0(J,3)",15x,
     &     "W0(J,4)",13x,"TAUABSO",15x,"tau",15x,"R")
c     907  format (i7,6(1pe22.14),I3)
 907  format (i7,6(1pe22.14),(1pe15.7))
 908  format (" THE NEXT DT=",e11.3)
 909  format ("       XMDOT=",e11.3)

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function fnin

(

NPL

)

Definition at line 1462 of file accdisk.f.

C***************************
C
      IMPLICIT REAL*8 (a-h,o-z)
      fnin = 1
C
      DO 10 i=1,npl
      odd = 2.*i+1
      evn = 2.*i
      fnin= fnin*evn/odd
10    CONTINUE
C
      RETURN

Here is the caller graph for this function:

double precision function fnjn

(

NPL

)

Definition at line 1478 of file accdisk.f.

C***************************
C
      IMPLICIT REAL*8 (a-h,o-z)
      fnjn = 3.1415926/2.0d0
C
      DO 10 i=1,npl+1
      odd = 2.*i-1
      evn = 2.*i
      fnjn= fnjn*odd/evn
10    CONTINUE
C
      RETURN

Here is the caller graph for this function:

subroutine fun	(		XI,
		dimension(2)	YI,
		dimension(2)	YP
	)

Definition at line 984 of file accdisk.f.

c     input  : XI=r,YI(1)=eta,YI(2)=dl
c     output : YP(1)=d(eta)/dr, YP(2)=d(l)/dr
C*******************************************************
C
      include 'para.h'
C
      dimension yi(2),yp(2)
c     
      cc2=cc*cc
c     ----- rg/c^2 -----
      unit   = rg/cc2
C     
      eta   = yi(1)
c     ---- ETA must be > 0 ----
      IF (eta .LT. 0.) THEN
         kpr=1
         return
      ENDIF
c     -------------------------
      dl    = yi(2)
      r     = xi
      r2    = r*r
      r3    = r2*r
      r4    = r3*r
      rm1   = r-1.
      al2   = alpha**2
      ok    = (r/rm1)*sqrt(.5/r3)
      xlk   = r2*ok
      dlnok =-.5/r-1./rm1
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c             -d     Ok               dlnLk
c     DLNRH = -- ln (---)   ; DLKDR = -----
c             dr      r                dr
c
c     XL    : dL-Lk = L (angular momentum at this radius)     
c     XLL   : L-Lin 
c     XLLK2 : 2*DL*LK+(DL)^2 => (LK+DL)^2 - LK^2 = L^2 - LK^2
c     WS    : Height integrated Pressure/Surface Density 
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++     
      dlnrh = 1./r-dlnok
      dlkdr = xlk*(1.5/r-1./rm1)
      xl    = dl+xlk
      xll   = xl-xlin
      xll2  = xll*xll
      xll3  = xll2*xll
      xll7  = xll**7
      eta2  = eta*eta
      eta5  = eta**5
      xllk2 = 2.*dl*xlk+dl*dl
      ws    = xll2*eta/(al2*r2)

c     +++++++++++++++++++++++++++++++
c     PD    : Pressure(P0)/Density(rho0) ratio at equatorial plane(z=0).
c     HH    : Scale height at z=0.
c     PR    : P0 (Pressure at z=0)
c     DE    : rho (Density at z=0)
c     +++++++++++++++++++++++++++++++
      
      pd    = (ain/ainp1)*ws
      hh    = sqrt(2.*xnpl1*pd)/ok
      pr    = (xmdot*xll/(4.*pi*r2*alpha*ainp1*hh))*(cc/(rg*rg))
      
c     +++ For Error Check +++ 
      if (pr.le.0.) then
         kpr=1
         return
      end if
c     +++++++++++++++++++++++
      de    = pr/(pd*cc*cc)
c     -----------------------------------------------------------------
c     If optically thick or r=rout1 ==> TE derived from P = Pgas+Prad  
c     If optically thin             ==> TE derived from P = Pgas
c     -----------------------------------------------------------------
      if ((kr.eq.1).or.(r.eq.rout1)) then
         te = temp(de,pr)
c     TE = temp2(DE,PR)
      else
         te = pr/de/rmu
      endif
      beta  = rmu*de*te/pr
      beta4 = beta**4
      demn  = ain*de
      temn  = 2.*te/3.
      xkff  = xkkr*demn*(temn**(-3.5))
      xk    = xkes+xkff
      xkxe  = xk/xkes
      xkxf  = xkff/xk
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c     QBR   : Brems cooling [erg/s/cm^2] (BHAD P248 eq.(8.64))
c     taup  : Qbrems/Qradiation(?)
c     B10   : 1-beta
c     B11   : 4-3*beta
c     B12   : beta+12*(gamma-1)*(1-beta)
c     !GG1 and GG3 are the specific heats ratio
c     GG1   : Gamma1
c     GG3   : Gamma3
c     
c     A1    : 2*Gamma1
c     A2    : Gamma1-1
c     A3    : Gamma3-1
c     A4    : 3*Gamma1-1 
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      qbr   = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
      taup  = qbr/(8.*aa*cc*(te**4))
      b10   = 1.-beta
      b11   = 4.-3.*beta
      b12   = beta+12.*gm1*b10
      gg1   = beta+(b11*b11)*gm1/b12
      gg3   = 1.+(gg1-beta)/b11
      a1    = 2.*gg1
      a2    = gg1-1.
      a3    = gg3-1.
      a4    = 3.*gg1-1.
c     ------------------------------------
c     QP    : Q^+ (Energy generation)  
c     QM    : Q^- 
c     FF    : Q^- (Brems cooling)
c     C1LL  : (L^2-Lk^2)/R^3 * SIGMA/PI
c     ------------------------------------
      qp    = 2.*al2*xl/xll
      if ((kr.eq.1).or.(r.eq.rout1)) then
         qm    = ff1*beta4*xll7*eta5/(xkxe*r4)
      else
         ff    = qbr*rg/cc/cc/cc
         qm    = 2.*pi*r2*eta/(xmdot*ws)*ff
      endif
      c1ll  = xllk2/(r3*ws)
C     
C     ***** MATRIX ELEMENTS *****
C
c | D11  D12 | = |      1                    (eta+1)/(L-Lin)           |
c | D21  D22 | = | (3*Gam3-1)/2 (2*Gam1*eta-(Gam3-1)*alpha**2)/(L-Lin) |
c
c     | c1 | = (eta+1)/R + DLNRH + Sigma/Pi * (L^2-Lk^2)/R^3
c     | c2 | =  eta*(2Gam1/R - (Gam1-1)DLNRH) - (Gam3-1)*(Q^+ - Q^-)/R)
c
      d11 = 1.
      d12 = (eta+1.)/xll
      d21 = a4/2.
      d22 = (a1*eta-a3*al2)/xll
      c1  = (eta+1.)/r+dlnrh+c1ll
      c2  = eta*(a1/r-a2*dlnrh)-a3*(qp-qm)/r
      c1d = c1-d12*dlkdr
      c2d = c2-d22*dlkdr
C
c++++++++++++++++++++++++++++++++++    
c    DD   : D (Matsumoto et al. 1984 eq(3.13))
c    XN1D : N_1 eq(3.14)
c    XN2D : N_2 eq(3.14)
c++++++++++++++++++++++++++++++++++
      dd   = d11*d22-d12*d21
      xn1d = c1d*d22-c2d*d12
      xn2d = d11*c2d-d21*c1d
C
C ***** F1 AND F2 *****
C
      yp(1) = xn1d/dd
      yp(2) = xn2d/dd
      kpr=0
C
C
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine intgrt	(		R,
			ETA,
			DL,
			IC
	)

Definition at line 803 of file accdisk.f.

C****************************************
C
C ------------------------------------------------------
C     NUMERICAL INTEGRATION OF STEADY DISK EQUATIONS
C     BY GEAR'S METHOD
C ------------------------------------------------------
C
      include 'para.h'
C
      dimension yy(2)
      dimension f(2),xx(jmax)
      dimension yscal(2)
      dimension oldyy(2)
      EXTERNAL fun,jac
C
c      print *,"INTGRT"
      IF ((iprint(1) .NE. 0).and.(imore.eq.2)) WRITE (6,300)
C
c ----------------------------------------------------------      
c     htry=hini : first dr for steady solution (t=0)
c     ya1,ya2   : eta, dl       
c ----------------------------------------------------------
c
      eps     = 1.e-8
      x0(1)   = r
      y0(1,1) = eta
      y0(2,1) = dl
      yy(1)   = eta
      yy(2)   = dl
      nv=2
      htry=hini
      ya1=abs(yy(1))
      ya2=abs(yy(2))
      if (ya1.ge.eps) yscal(1)=ya1
      if (ya1.lt.eps) yscal(1)=eps
      if (ya2.ge.eps) yscal(2)=ya2
      if (ya2.lt.eps) yscal(2)=eps
c      yscal(2)=abserr
C ---- MESH INTERVAL ----
c      aldr = (log10(ROUT1)-log10(RIN))/(mesh-50.)
      aldr = (log10(rout1)-log10(rin))/mesh
      xx(1)= rout1
      xold = rout1
      jj=2
10    continue
c ++++ Mesh interval 2 ++++
c      XX(JJ)=XOLD/(10.**aldr)
c ---- For r < 1000rg case -----
c      IF (XOLD.GT.1000.) XX(JJ)=XOLD/(10.**(aldr/3.))
      IF (xold.GT.1000.) xx(jj)=xold/(10.**(aldr/1.))
      IF (xold.LE.1000.) xx(jj)=xold/(10.**aldr/1.)
c ---- For r < 3.2rg  case -----
      dxx = xx(jj-1)-xx(jj)
c      write(6,*) "before",jj,xx(jj),dxx
c      pause
c      IF (dxx.LE.0.05) XX(JJ)=XOLD-0.01
c      IF (dxx.lt.0.5) XX(JJ)=XOLD-0.01
c      IF (XOLD.LE.3.0) XX(JJ)=XOLD-0.005      
      IF (dxx.LE.0.3) xx(jj)=xold-0.3*dxx   
cc      IF (XOLD.LE.4.0) xx(jj)=xold-0.3*dxx   
c      IF (XOLD.LE.3.5) xx(jj)=xold-0.01
c      IF (dxx.LE.0.1) XX(JJ)=XOLD-0.01
      if (xold .le. 3.0) xx(jj)=xold-0.01
c      dxx = xx(jj-1)-xx(jj)

c      if (xx(jj).lt.3.5) then 
c         write(6,*) "after ",jj,xx(jj),dxx
c         pause
c      end if

c +++++++++++++++++++++++++
      IF (xx(jj).LT.rin) GOTO 11
      xold = xx(jj)
      jj = jj+1
      GOTO 10
11    CONTINUE
c ---- Mesh (in real) ----
      jjx = jj
C ---- Start Integration at each radius ---- 
c
      do 2000 j=1,jjx-1
C
      r  = x0(j)
      xend = xx(j+1)
c ---- IOD : ITR count at each grid. ----
      iod=0
C
 500  continue
C      print *,j,xx(j)
      iod=iod+1
c 600  continue
      oldr=r
      oldhtry=htry
      oldyy(1)=yy(1)
      oldyy(2)=yy(2)
C ----------------------------------------------------------------------
      call fun (r,yy,f)
c              input  => r, eta, dl == r, yy(1),yy(2)
c              output => dy dx        == f 
c      for 
C ----------------------------------------------------------------------
C ---------------------------------------------------------------------
      call stifbs (yy,f,nv,r,htry,relerr,yscal,hdid,hnext,fun,icon)
c     : To solve for stiff problem.     
C ----------------------------------------------------------------------
C
c ----------- CHECK -----------------------------------------------
c      call fun (r,yy,f)
c      if (kpr.eq.1) then
c          r=oldr
c          htry=oldhtry/2.
c          yy(1)=oldyy(1)
c          yy(2)=oldyy(2)
c          goto 600
c      endif
c --------------------------------------------------------------------- 
      IF ((icon .EQ. 10000) .OR. (icon .EQ. 15000)) THEN
         if (imore.eq.2) WRITE (6,100) r,eps
      ELSE IF (icon .EQ. 16000) THEN
         if (imore.eq.2) WRITE (6,120)
         ic = -1
         RETURN
      ELSE IF (icon .EQ. 30000) THEN
         if (imore.eq.2) WRITE (6,140)
         stop
      ENDIF
 100  FORMAT (1h ,'TOLERANCE RESET.     R = ',1pe15.7,5x,'EPSR = ',1pe15
     *     .7)
 120  FORMAT (1h ,'NO CONVERGENCE IN CORRECTOR ITERATION')
 140  FORMAT (1h ,'INVALID INPUT TO ODGE')
      IF ((mod(iod,1000).EQ.0).and.(imore.eq.2))
     &     WRITE (6,*) r,yy(1),yy(2)
C     
C     -- OBTAIN RSTEP --
      if (r.eq.xend) then           ! xend=xx(j+1)
         rstep = r-xx(j+2)
      else
         rstep = r-xx(j+1)
      endif
      htry = -min(abs(hnext),rstep)
      ya1  = abs(yy(1))
      ya2  = abs(yy(2))
      if (ya1.ge.eps) yscal(1) = ya1
      if (ya1.lt.eps) yscal(1) = eps
c     yscal(1)=max(abs(yy(1)),1.e-8)
c     yscal(2)=max(abs(yy(2)),1.e-8)
      if (ya2.ge.eps) yscal(2) = ya2
      if (ya2.lt.eps) yscal(2) = eps
c     yscal(2)=abserr
      if (r.gt.xend) goto 500
C     
      x0(j+1) = r
      y0(1,j+1) = yy(1)
      y0(2,j+1) = yy(2)
C     
      ip1=iprint(1)
      IF (ip1 .EQ. 0) GOTO 1000
      IF (((mod((j+1),ip1) .EQ. 1) .OR. (ip1 .EQ. 1)).and.(imore.eq.2))
     &     WRITE (6,220) r,yy(1),yy(2),iod
 220  FORMAT (1h ,3(1pe15.7),i7)
C     
 2000 continue
 1000 CONTINUE
C     
      eta = yy(1)
      IF (eta .GT. 1.) THEN
         ic = 0
      ELSE
         ic = 1
      ENDIF
C     
c     JX = J+1
      jx = j
C     
 300  format (7x,"R",13x,"ETA",13x,"DL")
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine intmdl

(

)

Definition at line 199 of file accdisk.f.

C**************************
C
C  -------------------------------------------------------------
C     INITIAL CONDITION
C     STEADY MODEL OF THIN ACCRETION DISK AROUND A BLACK HOLE
C     TRANSONIC SOLUTION IS FOUND BY ADJUSTING LIN.
C  -------------------------------------------------------------
C
      include 'para.h'
C
      il = 0
C
C -------------------------
c      DO 100 UNTIL (IC .EQ. 0)
 200                         continue
C -------------------------
c
c ----------------------------------------------------------------------
c     XLIN  : L_in, angular momentum at marginally stable orbit.
c     XLINP : Lower limit of the inner angular momentum
c     XLINM : Upper limit of the inner angular momentum
c --------------------------------------------------------------------- 

      xlin = 0.5*(xlinp+xlinm)
      r = rout1

      if (kr .eq. 1) then 
         CALL sadmk(r,eta,dl)
      else if (kr .eq. 0) then 
         call sadmn(r,eta,dl)
      end if 
cc                do 2 i=1,2
cc                if (i.eq.1) r=rout1
cc                if (i.eq.2) r=rout2
c                call sadmn(r,eta,dl)
cc2               continue
c -------------------------------------------- 
      if (imore.eq.2) print *,"   IL=",il+1 ! IL -> Iteration count 
      CALL intgrt(r,eta,dl,ic)
      IF (ic .EQ. 1 ) xlinm = xlin
      IF (ic .EQ. -1) xlinp = xlin
      il = il+1
      
      IF (il .GT. 50) THEN      ! Set Max Iteration
         if (imore.eq.2) WRITE (6,81) il
         if (imore.eq.2) write (6,82) 
         if (imore.eq.2) write (6,*) "xlin=", xlin, "IC=",ic
         stop
      ENDIF
C
      if (imore.eq.2) WRITE (6,80) xlin,r,ic
 80   FORMAT (1h ,'LIN = ',1pe22.14,5x,'R = ',1pe15.7,5x,'IC = ',i5)
 81   format ("LOOP COUNT (INTMDL) IS OVER MAXIMUM  IL=",i3)
 82   format ("  IL  IC     LIN")
C
      if (ic.ne.0) goto 200   
100   CONTINUE
C
      CALL cnvrbs
C
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine jac	(		xi,
		dimension(2)	yi,
		dimension(2)	yp,
		dimension(2)	dfdr,
		dimension(2,2)	pdj
	)

Definition at line 1151 of file accdisk.f.

C**************************************
C
      include 'para.h'
C
      dimension yi(2),pdj(2,2),dj(2,2,3),cj(2,3),cdj(2,3),dnj(2,3),
     &          ddj(3),dfdr(2),yp(2)
      cc2=cc*cc
      unit   = rg/cc2
C
      eta   = yi(1)
      IF (eta.LT.0.) THEN
      kpr=1
      return
      ENDIF
      dl    = yi(2)
      r     = xi
      r2    = r*r
      r3    = r2*r
      r4    = r2*r2
      r7    = r3*r3*r
      r8    = r7*r
      rm1   = r-1.
      al2   = alpha*alpha
C
      ok    = (r/rm1)*sqrt(.5/r3)
      xlk   = r2*ok
      xl    = xlk+dl
      xll   = xl-xlin
      xll2  = xll*xll
      xll3  = xll**3
      xll7  = xll**7
      eta2  = eta*eta
      eta5  = eta**5
      xl2   = xl*xl
      ws    = xll2*eta/(al2*r2)
      pd    = (ain/ainp1)*ws
      xlk2  = xlk*xlk
      xllk2 = 2.*dl*xlk+dl*dl
      dlnok =-0.5/r-1./rm1
      dlnrh = 1./r-dlnok
      d2lnrh=-1.5/r2-1./(rm1**2)
      dlkdr = xlk*(1.5/r-1./rm1)
      hh    = sqrt(2.*xnpl1*pd)/ok
      pr    = (xmdot*xll/(4.*pi*r2*alpha*ainp1*hh))*(cc/(rg**2))
      if (pr.le.0.) then
          kpr=1
          return
      endif
      de    = pr/(pd*cc*cc)
      if ((kr.eq.1).or.(r.eq.rout1)) then
         te = temp(de,pr)
      else
         te = pr/de/rmu
      endif
      beta  = rmu*de*te/pr
      beta4 = beta**4
      demn  = ain*de
      temn  = 2.*te/3.
      xkff  = xkkr*demn*(temn**(-3.5))
      xk    = xkes+xkff
      xkxe  = xk/xkes
      xkxf  = xkff/xk
      qbr   = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
      taup  = qbr/(8.*aa*cc*(te**4))
      b10   = 1.-beta
      b11   = 4.-3.*beta
      b12   = beta+12.*gm1*b10
      gg1   = beta+(b11**2)*gm1/b12
      gg3   = 1.+(gg1-beta)/b11
      a1    = 2.*gg1
      a2    = gg1-1.
      a3    = gg3-1.
      a4    = 3.*gg1-1.
      qp    = 2.*al2*xl/xll
      if ((kr.eq.1).or.(r.eq.rout1)) then
         qm    = ff1*beta4*xll7*eta5/(xkxe*r4)
      else
         ff    = qbr*rg/cc/cc/cc
         qm    = 2.*pi*r2*eta/(xmdot*ws)*ff
      endif
      c1ll  = xllk2/(r3*ws)
C
C ***** MATRIX ELEMENTS *****
C
      d11 = 1.
      d12 = (eta+1.)/xll
      d21 = a4/2.
      d22 = (a1*eta-a3*al2)/xll
      c1  = (eta+1.)/r+dlnrh+c1ll
      c2  = eta*(a1/r-a2*dlnrh)-a3*(qp-qm)/r
      c1d = c1-d12*dlkdr
      c2d = c2-d22*dlkdr
      dd  = d11*d22-d12*d21
      xn1d= c1d*d22-c2d*d12
      xn2d= d11*c2d-d21*c1d
C
C ***** CAL. JACOBIAN *****
C
      b13   = beta*b10/b11
      dbdl  =-b13*(8./xll)
      dbde  =-b13*(4.5/eta)
      dbdr  = b13*(7./r)
      if ((kr.eq.1).or.(r.eq.rout1)) then
         dg1db = 1.-b11*gm1*(3.*beta+4.+12.*gm1*(2.-3.*beta))/(b12*b12)
         dg3db = (dg1db-1.)/b11+3.*(gg1-beta)/(b11*b11)
      else
         dg1db = 0.
         dg3db = 0.
      endif
      dpdl  = 0.
      dpde  =-0.5/eta
      dpdr  =-1./r
      dddl  =-2./xll
      ddde  =-1.5/eta
      dddr  = 1./r
      dwsde = dpde-ddde
      dwsdl = dpdl-dddl
      dwsdr = dpdr-dddr
      dtdl  = (dpdl-beta*dddl)/b11
      dtde  = (dpde-beta*ddde)/b11
      dtdr  = (dpdr-brta*dddr)/b11
      dtdr  =-2./r 
      dkdl  = xkxf*(dddl-3.5*dtdl)
      dkde  = xkxf*(ddde-3.5*dtde)
      dkdr  = xkxf*(dddr-3.5*dtdr) 
      if ((kr.eq.1).or.(r.eq.rout1)) then
         dqmde = qm*(4.*dbde/beta+5./eta-dkde)
         dqmdl = qm*(4.*dbdl/beta+7./xll-dkdl)
         dqmdr = qm*(4.*dbdr/beta-4./r  -dkdr)
      else
         dqmde = qm*(-1./eta)
         dqmdl = qm*(-2./xll)
         dqmdr = qm*(0.5*(1.-3./r)/rm1)
      endif
      dqpdl = qp*(1./xl-1./xll)
      dqpdr = 0.
      d21db = 1.5*dg1db
      d22db = (2.*dg1db*eta-dg3db*al2)/xll
      dc2db = eta*(2.*dg1db/r-dg1db*dlnrh)-dg3db*(qp-qm)/r
C
C
      dj(1,1,1) = 0.
      dj(1,1,2) = 0.
      dj(1,1,3) = 0.
      dj(1,2,1) = 1./xll
      dj(1,2,2) =-d12/xll
      dj(1,2,3) = 0.
      dj(2,1,1) = d21db*dbde
      dj(2,1,2) = d21db*dbdl
      dj(2,1,3) = d21db*dbdr
      dj(2,2,1) = d22db*dbde+2.*gg1/xll
      dj(2,2,2) = d22db*dbdl-d22/xll
      dj(2,2,3) = d22db*dbdr
C
      cj(1,1)   = 1./r-c1ll*dwsde
      cj(1,2)   = 2.*xl/(r3*ws)-c1ll*dwsdl
      cj(1,3)   =-(eta+2.5)/r2-(3./r)*c1ll-1./(rm1*rm1)
      cj(2,1)   = a1/r-a2*dlnrh+a3*dqmde/r+dc2db*dbde
      cj(2,2)   =-a3*(dqpdl-dqmdl)/r+dc2db*dbdl
      cj(2,3)   =-eta*(a1-1.5*a2)/r2+eta*a2/sqrt(rm1)
     &           +a3*(dqmdr*r+qp-qm)/r2+dc2db*dbdr 
C
      cdj(1,1)  = cj(1,1)-dj(1,2,1)*dlkdr
      cdj(1,2)  = cj(1,2)-dj(1,2,2)*dlkdr
      cdj(1,3)  = cj(1,3)-dj(1,2,3)*dlkdr
     &           -d12*(dlkdr*(1.5/r-1./rm1)+xlk*(-1.5/r2+1./sqrt(rm1)))
      cdj(2,1)  = cj(2,1)-dj(2,2,1)*dlkdr
      cdj(2,2)  = cj(2,2)-dj(2,2,2)*dlkdr
      cdj(2,3)  = cj(2,3)-dj(2,2,3)*dlkdr
     &           -d22*(dlkdr*(1.5/r-1./rm1)+xlk*(-1.5/r2+1./sqrt(rm1)))
C
      dnj(1,1) = cdj(1,1)*d22+c1d*dj(2,2,1)-cdj(2,1)*d12-c2d*dj(1,2,1)
      dnj(1,2) = cdj(1,2)*d22+c1d*dj(2,2,2)-cdj(2,2)*d12-c2d*dj(1,2,2)
      dnj(1,3) = cdj(1,3)*d22+c1d*dj(2,2,3)-cdj(2,3)*d12-c2d*dj(1,2,3)
      dnj(2,1) = dj(1,1,1)*c2d+d11*cdj(2,1)-dj(2,1,1)*c1d-d21*cdj(1,1)
      dnj(2,2) = dj(1,1,2)*c2d+d11*cdj(2,2)-dj(2,1,2)*c1d-d21*cdj(1,2)
      dnj(2,3) = dj(1,1,3)*c2d+d11*cdj(2,3)-dj(2,1,3)*c1d-d21*cdj(1,3)
C
      ddj(1)  = dj(1,1,1)*d22+d11*dj(2,2,1)-dj(1,2,1)*d21-d12*dj(2,1,1)
      ddj(2)  = dj(1,1,2)*d22+d11*dj(2,2,2)-dj(1,2,2)*d21-d12*dj(2,1,2)
      ddj(3)  = dj(1,1,3)*d22+d11*dj(2,2,3)-dj(1,2,3)*d21-d12*dj(2,1,3)
C
      dfdr(1)   = (dnj(1,3)*dd-xn1d*ddj(3))/(dd*dd)
      dfdr(2)   = (dnj(2,3)*dd-xn2d*ddj(3))/(dd*dd)
      pdj(1,1)  = (dnj(1,1)*dd-xn1d*ddj(1))/(dd*dd)
      pdj(1,2)  = (dnj(1,2)*dd-xn1d*ddj(2))/(dd*dd)
      pdj(2,1)  = (dnj(2,1)*dd-xn2d*ddj(1))/(dd*dd)
      pdj(2,2)  = (dnj(2,2)*dd-xn2d*ddj(2))/(dd*dd)
      kpr=0
C
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine lubksb	(	dimension(np,np)	a,
		integer	n,
		integer	np,
		integer, dimension(n)	indx,
		dimension(n)	b
	)

Definition at line 3245 of file accdisk.f.

      implicit real*8 (a-h,o-z)
      common /lscal/ kpr
      integer n,np,indx(n),kpr
c      real a(np,np),b(n)
      dimension a(np,np),b(n)
      integer i,ii,j,ll
      ii=0
      do 12 i=1,n
         ll=indx(i)
         sum=b(ll)
         b(ll)=b(i)
         if (ii.ne.0) then
             do 11 j=ii,i-1
                sum=sum-a(i,j)*b(j)
 11          continue
         else if (sum.ne.0.) then
             ii=i
         endif
         b(i)=sum
 12   continue
      do 14 i=n,1,-1
         sum=b(i)
         do 13 j=i+1,n
            sum=sum-a(i,j)*b(j)
 13      continue
         b(i)=sum/a(i,i)
 14   continue
      return

Here is the caller graph for this function:

subroutine ludcmp	(	dimension(np,np)	a,
		integer	n,
		integer	np,
		integer, dimension(n)	indx,
			d
	)

Definition at line 3181 of file accdisk.f.

      implicit real*8 (a-h,o-z)
      common /lscal/ kpr
      integer n,np,indx(n),nmax,kpr
c      real d,a(np,np),tiny
      dimension a(np,np)
      parameter(nmax=500,tiny=1.0e-20)
      integer i,imax,j,k
c      real aamax,dum,sum,vv(nmax)
      dimension vv(nmax)
      d=1.
      do 12 i=1,n
         aamax=0.
         do 11 j=1,n
            if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
 11      continue
c     if (aamax.eq.0.) pause "singular matrix in ludcmp"
            if (aamax.eq.0.) stop "singular matrix in ludcmp"
         vv(i)=1./aamax
 12   continue
      do 19 j=1,n
         do 14 i=1,j-1
            sum=a(i,j)
            do 13 k=1,i-1
               sum=sum-a(i,k)*a(k,j)
 13         continue
            a(i,j)=sum
 14      continue
         aamax=0.
         do 16 i=j,n
            sum=a(i,j)
            do 15 k=1,j-1
               sum=sum-a(i,k)*a(k,j)
 15         continue
            a(i,j)=sum
            dum=vv(i)*abs(sum)
            if (dum.ge.aamax) then
                imax=i
                aamax=dum
            endif
 16      continue
         if (j.ne.imax) then
             do 17 k=1,n
                dum=a(imax,k)
                a(imax,k)=a(j,k)
                a(j,k)=dum
 17          continue
             d=-d
             vv(imax)=vv(j)
         endif
         indx(j)=imax
         if (a(j,j).eq.0.) a(j,j)=tiny
         if (j.ne.n) then
             dum=1./a(j,j)
             do 18 i=j+1,n
                a(i,j)=a(i,j)*dum
 18          continue
         endif
 19   continue
      return

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine lw1d

(

)

Definition at line 1615 of file accdisk.f.

C**************************
C
C ----------------------------------------------------------------------
C     TIME INTEGRATION OF HYDRODYNAMIC EQUATIONS OF THIN ACCRETION
C     DISK BY LAX-WENDROFF METHOD.
C ----------------------------------------------------------------------
C
      include 'para.h'
C
      dimension u(jmax,4),f(jmax,4),h(jmax,4),y(jmax,4),q1(jmax),
     *q2(jmax),d(jmax,4)
      dimension wk(jmax,4)
      dimension drp(jmax),drm(jmax),dr(jmax)
      dimension betaj(jmax),tej(jmax)
      dimension wt(jmax,4),tauabst(jmax),taut(jmax)
c      dimension dexise(jmax),uo(jmax)
      dimension okc(jmax),hhc(jmax),prc(jmax),tec(jmax),betac(jmax),
     &          ssc(jmax),dec(jmax),betao(jmax)
      dimension des(jmax),oes1(jmax),ou(jmax,4)
C
C
      gm1    = gm-1
      cc2    = cc*cc
      unit   = rg/cc2
C
      DO 5 j=2,jx-1                   
         drp(j) = x1(j+1)-x1(j)         ! Á°¿Êº¹Ê¬ 
         drm(j) = x1(j)-x1(j-1)         ! ¸åÂູʬ
         dr(j)  = 0.5*(drp(j)+drm(j))   ! dr¤ÎÄêµÁ¡ÊÃæ¿´º¹Ê¬¡Ë         DRP(1) = X1(2) -X1(1)
         DRM(JX)= X1(JX)-X1(JX-1)
         DR(1)  = DRP(1)
         DR(JX) = DRM(JX)

      do 152 j=1,jx-1
         nl(j)=0
152   continue
c
151   continue
      XLUMI=0.  ! Total luminosity
      xltn=0.   ! optically thin luminosity
      xltk=0.   ! optically thick luminosity

c      print *,W0(1,1)
c      print *,"L1"

      DO 10 J=1,JX
         R    = X1(J)
         RM1  = R-1
         R2   = R*R
         R3   = R2*R
         OK   = (R/RM1)*SQRT(0.5/R3)
         DE   = W0(J,1)
         VX   = W0(J,2)
         VP   = W0(J,3)
         TE   = W0(J,4)
         wt(j,1) = W0(J,1)
         wt(j,2) = W0(J,2)
         wt(j,3) = W0(J,3)
         wt(j,4) = W0(J,4)
         tauabst(j) = tauabso(j)
         taut(j)    = tauo(j)
         tau        = tauo(j)
         taup       = tauabso(j)
c
      if (j.eq.jx) then
         PR   = (AA/3.)*(TE**4)+RMU*DE*TE
      else
         if (kcool.eq.1) fac = 1./((1.+1./taup/(1.5*tau+sqrt(3.))))
         if (kcool.eq.0) fac = 1.
c
         PR   = (AA/3.)*(TE**4)*FAC+RMU*DE*TE
      endif
c
         BETA = RMU*DE*TE/PR
         HH   = SQRT(B5*PR/DE)/OK
         SS   = 2*AIN*DE*HH*RG
         WW   = 2*AINP1*PR*HH*RG
         WC   = WW/CC2
         WS   = WC/SS
         A    = 3*(1-BETA)+BETA/GM1
         EE   = (A+0.5)*WS+0.5*(VX*VX+VP*VP)       ! total energy
         RS   = R*SS
         BETAJ(J)= BETA

         U(J,1) = RS                              ! d(r*SS)/dt 
         U(J,2) = RS*VX                           ! d(r*SS*vr)/dt 
         U(J,3) = RS*VP                           ! d(r*SS*vp)/dt
         U(J,4) = RS*EE                           ! d(r*SS*ee)/dt

         F(J,1) = RS*VX                         ! d(r*SS*vr)/dr
         F(J,2) = RS*(VX*VX+WS)                 ! d(r*SS*vr^2+r*WW)/dr
         F(J,3) = RS*(VX*VP+ALPHA*WS)           ! d(r^2*SS*vr*vp-Trp)/dr
         F(J,4) = RS*(VX*(EE+WS)+ALPHA*WS*VP)  ! d()/dr

         VK     = R*OK
         XLK    = R*VK
         DLNOR  =-(0.5+R/RM1)
         DLNLR  = 2+DLNOR
         DOKDR  = OK*DLNOR/R
         DLKDR  = XLK*DLNLR/R
         DEMN   = AIN*DE
         TEMN   = 2*TE/3.0
         XKFF   = XKKR*DEMN*(TEMN**(-3.5))
         TAUFF  = XKFF*DE*HH*RG
         TAUES  = XKES*DE*HH*RG
         TAU    = TAUES+TAUFF
         qbr    = 2.0*XEFF*AJN*DE*DE*SQRT(TE)*HH*RG
         taup   = qbr/(8.*aa*cc*(te**4))
c
         if (kcool .eq. 1) then 
            taueff = sqrt(taup*(tau+2./sqrt(3.)))
         else if (kcool .eq. 0) then 
            taueff = sqrt(3*tau*tauff)
         end if
c
         if (ns.eq.lns+1) then
c --------- oes1(j): total energy = E - 0.5*(vr^2*vp^2) 
            oes1(j)=u(j,4)/u(j,1)-((u(j,2)/u(j,1))*(u(j,2)/u(j,1))
     &           +(u(j,3)/u(j,1))*(u(j,3)/u(j,1)))/2.
            ou(j,1)=u(j,1)
            ou(j,2)=u(j,2)
            ou(j,3)=u(j,3)
            ou(j,4)=u(j,4)
            tauo(j)=tau
            tauabso(j)=taup
            betao(j)=beta
         endif

         XK     = XKES+XKFF
         FD     = 16*AA*(TE**4)/(3*XK*DE*HH*CC2)
         FF     = 2.0*XEFF*AJN*DE*DE*SQRT(TE)*HH*RG*RG/CC2/CC
c--- 2002-04-11 ---
         if (kcool .eq.1) then 
           FM     = 8.*AA*(TE**4)/(1.5*TAU+SQRT(3.)+1./TAUP)*UNIT
         else
           FM     = FD
         end if
c
         H(J,1) = 0
         H(J,2) = WC*(1-DLNOR)+SS*(2*VP*VK+VP*VP)    ! (11.30)¤Î±¦ÊÕ
         H(J,3) =-SS*(VX*VP+VX*DLKDR+ALPHA*WS)       
         H(J,4) =-R*RS*(VX*VP+ALPHA*WS)*DOKDR-R*FM   ! (11.32)¤Î±¦ÊÕ
         IF (NS.EQ.LNS+1) THEN
cc             FFTR    = FM/UNIT/XMCRIT/CC*16.
             FFTR    = FM/UNIT/XMCRIT/CC
             XLumi   = XLumi+2.*PI*(R*RG)*(DR(J)*RG)*FFTR
          if (j.eq.jx) xljx=2.*PI*(R*RG)*(DR(J)*RG)*FFTR
          if (taueff.le.1.) xltn=xltn+2.*PI*(R*RG)*(DR(J)*RG)*FFTR
          if (taueff.gt.1.) xltk=xltk+2.*PI*(R*RG)*(DR(J)*RG)*FFTR
         ENDIF
         betac(j) = beta
         okc(j)   = ok
         hhc(j)   = hh
         prc(j)   = pr
         tec(j)   = te
         ssc(j)   = ss
         dec(j)   = de
10    CONTINUE
c      IF ((NS.EQ.LNS+1).and.(kcau1+kcau2.eq.0)) 
c     &    WRITE (9,900) TIME-DT,XLumi,xltn,xltk,dt
c      print *,"L1.5"

c###### Thermal Conduction #################################### 
c    xkei: thermal conductivity => K=alpha*xkappa*C_s*Sigma*C_v*H
c
      if (ktcond.eq.1) then
c         do 12 j=2,jx-1
         do 12 j=1,jx-1
c            ef=3.*xkappa*(8.-7.*betac(j))
CCCc            cv=3.*(8.-7.*betac(j))/(8.-3.*betac(j)-3.*betac(j))
            cv=3./2.
cc          xnu=alpha*okc(j)*(hhc(j)**2.)
cc          xkei=xkappa*ssc(j)*cv*xnu
c          gm31=gm1*(4.-3.*betac(j))
c     &         /(betac(j)+12.*gm1*(1.-betac(j)))
c         ef=xkappa*gm31
          ef=xkappa
          xkei=alpha*ef*okc(j)*(hhc(j)**3)*prc(j)/tec(j)
          cn=cv*xkei*(tec(j+1)-tec(j))/drp(j)
     &       *sqrt(2./xnpl1)*ainp1*unit
c         if (nl(j).ge.1) cn=0.
         f(j,4)=f(j,4)-x1(j)*cn
         h(j,4)=h(j,4)-cn
12     continue
      endif
C########################################## 
C
c      print *,"LW2"
      DO 40 J=1,JX-1
         AVE1 = 0.5*(U(J+1,1)+U(J,1))
         AVE2 = 0.5*(U(J+1,2)+U(J,2))
         AVE3 = 0.5*(U(J+1,3)+U(J,3))
         AVE4 = 0.5*(U(J+1,4)+U(J,4))
         FSA1 = F(J+1,1)-F(J,1)
         FSA2 = F(J+1,2)-F(J,2)
         FSA3 = F(J+1,3)-F(J,3)
         FSA4 = F(J+1,4)-F(J,4)
         SWA1 = 0.5*(H(J+1,1)+H(J,1))
         SWA2 = 0.5*(H(J+1,2)+H(J,2))
         SWA3 = 0.5*(H(J+1,3)+H(J,3))
         SWA4 = 0.5*(H(J+1,4)+H(J,4))
         DTDX = DT/DRP(J)
         Y(J,1) = AVE1-0.5*DTDX*FSA1+0.5*DT*SWA1
         Y(J,2) = AVE2-0.5*DTDX*FSA2+0.5*DT*SWA2
         Y(J,3) = AVE3-0.5*DTDX*FSA3+0.5*DT*SWA3
         Y(J,4) = AVE4-0.5*DTDX*FSA4+0.5*DT*SWA4
         es1(j)=y(j,4)/y(j,1)-((y(j,2)/y(j,1))*(y(j,2)/y(j,1))
     &         +(y(j,3)/y(j,1))*(y(j,3)/y(j,1)))/2.
40    CONTINUE
         dtdx    = dt/dr(1)
         fsa1    = f(2,1)-f(1,1)
         fsa2    = f(2,2)-f(1,2)
         fsa3    = f(2,3)-f(1,3)
         fsa4    = f(2,4)-f(1,4)
         swa1    = 0.5*(h(1,1)+h(2,1))
         swa2    = 0.5*(h(1,2)+h(2,2))
         swa3    = 0.5*(h(1,3)+h(2,3))
         swa4    = 0.5*(h(1,4)+h(2,4))
         wk(1,1) = u(1,1)-dtdx*fsa1+dt*swa1
         wk(1,2) = u(1,2)-dtdx*fsa2+dt*swa2
         wk(1,3) = u(1,3)-dtdx*fsa3+dt*swa3
         wk(1,4) = u(1,4)-dtdx*fsa4+dt*swa4
C
c##### Artificial Viscosity #####
c What is kcau1 & kcau2 ? => Ans.: kcau1,2 = 0,or 1
c      if kcau1&2 = 0 -> no problem! 
c      if kcau1&2 = 1 -> problem! -> modify!

      if (kcau1+kcau2.eq.0) then
         QU1=QU
      else
         QU1=16. 
      endif

      DO 330 J = 1,JX-1
         VXP = W0(J+1,2)
         VXM = W0(J,2)
c        qu1=qu
c        if (nl(j).ge.1) qu1=qu*5.
         Q2(J)= QU1*ABS(VXP-VXM)
330   CONTINUE
c
      do 340 J=2,JX-1
         dtdx=dt/dr(j)
         DES(J) = DTDX*(Q2(J)*(es1(J+1)-es1(J))
     &           -Q2(J-1)*(es1(J)-es1(J-1)))
340   CONTINUE
         DTDX = DT/DR(1)
         DES(1)=DTDX*Q2(1)*(es1(2)-es1(1))
c
      DO 350 J=1,JX-1
         es1(J) = es1(J)+DES(J)
c         if (kcau1+kcau2.ne.0) es1(1)=oes1(1)
c         if (es1(j).le.0.) es1(j)=oes1(j)
350   CONTINUE

      do 360 j=1,jx-1
c      if ((y(1,1).le.0.).or.(es1(1).le.0.)) then
c      if ((y(j,1).le.0.).or.(es1(j).le.0.)) then
      if (nl(j).ge.2) then
         es1(j)=oes1(j)
         y(j,1)=ou(j,1)
         y(j,2)=ou(j,2)
         y(j,3)=ou(j,3)
         y(j,4)=ou(j,4)
      endif
360   continue
C
      kcau1=0
      do 553 j=1,jx-1
         if (y(j,1).le.0.e0) then  ! if de < 0, then kcau1=1
            nl(j)=nl(j)+1
            kcau1=1
            nlmax=max(nlmax,nl(j))
            goto 553
         endif
c         es=y(j,4)/y(j,1)-((y(j,2)/y(j,1))*(y(j,2)/y(j,1))
c     &      +(y(j,3)/y(j,1))*(y(j,3)/y(j,1)))/2.
c         if (es1(j).le.0.e14) then
         if (es1(j).le.0.e0) then
            nl(j)=nl(j)+1
            kcau1=1
c            print *,"!",ns
            nlmax=max(nlmax,nl(j))
         endif
553   continue
c------------------------------------------
                       if (kcau1.eq.1) then
c                          write(6,*) "kcau1=",kcau1
c                          pause
c------------------------------------------
      if (nlmax.le.10) then
         do 570 j=1,jx-1
            W0(J,1) = wt(j,1)
            W0(J,2) = wt(j,2)
            W0(J,3) = wt(j,3)
            W0(J,4) = wt(j,4)
            tauabso(j)=tauabst(j)
            tauo(j)=taut(j)
570      continue
         time=time-dt
         dt=dt*0.1
         time=time+dt
c
c         print *,"Check after 570 loop (nlmax < 10)"
c         print *,"time=",time,"dt=",dt
c         pause
c
         kback=1
         goto 151
      endif
         if (y(j,1).le.0.) then
            krs(j)=1
         endif
         if (es1(j).le.0.) then
            kes(j)=1
         endif
         call print
         stop
c---------------------------
                        else
c---------------------------
c      print *,"L3"
      DO 100 J=1,JX-1
         R   = 0.5*(X1(J)+X1(J+1))
         R2  = R*R
         R3  = R2*R
         RM1 = R-1
         OK  = R/RM1*SQRT(0.5/R3)
         RS  = Y(J,1)
         VX  = Y(J,2)/RS
         VP  = Y(J,3)/RS
         EE  = Y(J,4)/RS
         SS  = RS/R
         ES  = EE-0.5*(VX*VX+VP*VP)
         ec2 = es1(j)*cc2
         HHS = SQRT(B5*(AIN/AINP1)/SS)/OK
         beta= betao(j)
         te  = w0(j,4)

c ***** start 98 iteration! *****
c@@@@@ iss > 10 for converge (best value=10) @@@@@
      do 98 iss=1,20
c      if (j.eq.1) print *,beta,te
         obeta  = beta
c  
         A      = 3*(1-BETA)+BETA/GM1
         WW     = EC2*SS/(A+0.5)
         WC     = WW/CC2
         HH     = HHS*SQRT(WW)
         DE     = SS/(2*AIN*HH*RG)
         DEMN   = AIN*DE
         TEMN   = 2*TE/3.0
         XKFF   = XKKR*DEMN*(TEMN**(-3.5))
         TAUFF  = XKFF*DE*HH*RG
c         TAUFF  = 2.*XKFF*DE*HH*RG
c         TAUES  = 2.*XKES*DE*HH*RG
         TAUES  = XKES*DE*HH*RG
         TAU    = TAUES+TAUFF
         qbr    = 2.0*XEFF*AJN*DE*DE*SQRT(TE)*HH*RG
         taup   = qbr/(8.*aa*cc*(te**4))
         PR     = WW/(2*AINP1*HH*RG)
      if (kcool.eq.1) then 
         te = ondo(de,pr,taup,tau)
      else if (kcool.eq.0) then 
         te = temp(de,pr)
      end if
         BETA   = RMU*DE*TE/PR

c***** 2002-7-30 *****
         if (abs((beta-obeta)/beta) .lt. 1.e-6) go to 730
 98   continue
 730     continue
c ***** end 98 iteration! *****
c      taueff = sqrt(taup*(tau+2./sqrt(3.)))
c      taueffjx=sqrt(tauabso(jx)*(tauo(jx)+2./sqrt(3.)))
c
         if (kcool .eq. 1) then 
            taueff = sqrt(taup*(tau+2./sqrt(3.)))
         else 
            taueff = sqrt(3*tau*tauff)
         end if
c
      if (taueff.gt.1.4618e30) then
       beta  = betao(j)
c       tau  = tauo(j)
c       taup = tauabso(j)
c       te   = wt(j,4)
         RS  = ou(J,1)
         VX  = ou(J,2)/RS
         VP  = ou(J,3)/RS
         EE  = ou(J,4)/RS
         SS  = RS/R
         ec2 = oes1(j)*cc2
         HHS = SQRT(B5*(AIN/AINP1)/SS)/OK
         A   = 3*(1-BETA)+BETA/GM1
         WW  = EC2*SS/(A+0.5)
         WC  = WW/CC2
         HH  = HHS*SQRT(WW)
         DE  = SS/(2*AIN*HH*RG)
         DEMN   = AIN*DE
         TEMN   = 2*TE/3.0
         XKFF   = XKKR*DEMN*(TEMN**(-3.5))
         TAUFF  = XKFF*DE*HH*RG
         TAUES  = XKES*DE*HH*RG
         TAU    = TAUES+TAUFF
         qbr    = 2.0*XEFF*AJN*DE*DE*SQRT(TE)*HH*RG
         taup   = qbr/(8.*aa*cc*(te**4))
         PR     = WW/(2*AINP1*HH*RG)
         te     = ondo(de,pr,taup,tau)
      end if

      WS     = WC/SS
      F(J,1) = RS*VX
      F(J,2) = RS*(VX*VX+WS)
      F(J,3) = RS*(VX*VP+ALPHA*WS)
      F(J,4) = RS*(VX*(EE+WS)+ALPHA*WS*VP)
      VK     = R*OK
      XLK    = R*VK
      DLNOR  = -(0.5+R/RM1)
      DLNLR  = 2+DLNOR
      DOKDR  = OK*DLNOR/R
      DLKDR  = XLK*DLNLR/R

      XK     = XKES+XKFF
      fd     = 16*AA*(TE**4)/(3*XK*DE*HH*CC2)
c--- 2002-04-11 ---
      if (kcool .eq.1) then 
        FM     = 8.*AA*(TE**4)/(1.5*TAU+SQRT(3.)+1./TAUP)*UNIT
      else
        FM     = FD
      end if

      H(J,1) = 0
      H(J,2) = WC*(1-DLNOR)+SS*(2*VP*VK+VP*VP)
      H(J,3) = -SS*(VX*VP+VX*DLKDR+ALPHA*WS)
      H(J,4) = -R*RS*(VX*VP+ALPHA*WS)*DOKDR-R*FM
      betaj(j)= beta
      betac(j)= beta
      okc(j) = ok
      hhc(j) = hh
      prc(j) = pr
      tej(j) = te
      tec(j) = te
      ssc(j) = ss
      dec(j) = de
100   CONTINUE
c ----- End of L3 loop -----

c###### Thermal Conduction #################################### 
      if (ktcond.eq.1) then
c      do 103 j=2,jx-1
      do 103 j=1,jx-1
c         ef=3.*xkappa*(8.-7.*betac(j))
CCCc         cv=3.*(8.-7.*betac(j))/(8.-3.*betac(j)-3.*betac(j))
         cv=3./2.
cc         xnu=alpha*okc(j)*(hhc(j)**2.)
cc         xkei=xkappa*ssc(j)*cv*xnu
c         gm31=gm1*(4.-3.*betac(j))
c     &        /(betac(j)+12.*gm1*(1.-betac(j)))
c         ef = xkappa*gm31
         ef = xkappa
         xkei = alpha*ef*okc(j)*(hhc(j)**3)*prc(j)/tec(j)
         cn = cv*xkei*(tec(j+1)-tec(j))/drp(j)
     &      *sqrt(2./xnpl1)*ainp1*unit
c         if (nl(j).ge.1) cn=0.
         f(j,4)=f(j,4)-x1(j)*cn
         h(j,4)=h(j,4)-cn
103   continue
      end if 
c##########################################
c---------------------------------- 
                             end if
C----------------------------------
c      print *,"L4"
      DO 120 J=2,JX-1
         DTDX   = DT/DR(J)
         FSA1   = F(J,1)-F(J-1,1)
         FSA2   = F(J,2)-F(J-1,2)
         FSA3   = F(J,3)-F(J-1,3)
         FSA4   = F(J,4)-F(J-1,4)
         SWA1   = 0.5*(H(J,1)+H(J-1,1))
         SWA2   = 0.5*(H(J,2)+H(J-1,2))
         SWA3   = 0.5*(H(J,3)+H(J-1,3))
         SWA4   = 0.5*(H(J,4)+H(J-1,4))
         WK(J,1) = U(J,1)-DTDX*FSA1+DT*SWA1
         WK(J,2) = U(J,2)-DTDX*FSA2+DT*SWA2
         WK(J,3) = U(J,3)-DTDX*FSA3+DT*SWA3
         WK(J,4) = U(J,4)-DTDX*FSA4+DT*SWA4
c      endif
120   CONTINUE
C
C****** ARTIFICIAL VISCOSITY *****
C
c      print *,"L6"

      if (kcau1+kcau2.eq.0) then
         QB1=QB
      else
         QB1=16.
      endif

      DO 130 J = 1,JX-1
         VXP  = W0(J+1,2)
         VXM  = W0(J,2)
         Q1(J)= QB1*ABS(VXP-VXM)
130   CONTINUE
C
      DO 140 J=2,JX-1
      DTDX   = DT/DR(J)
      D(J,1) = DTDX*(Q1(J)*(U(J+1,1)-U(J,1))-Q1(J-1)*(U(J,1)-U(J-1,1)))
      D(J,2) = DTDX*(Q1(J)*(U(J+1,2)-U(J,2))-Q1(J-1)*(U(J,2)-U(J-1,2)))
      D(J,3) = DTDX*(Q1(J)*(U(J+1,3)-U(J,3))-Q1(J-1)*(U(J,3)-U(J-1,3)))
      D(J,4) = DTDX*(Q1(J)*(U(J+1,4)-U(J,4))-Q1(J-1)*(U(J,4)-U(J-1,4)))
140   CONTINUE
C
      DTDX   = DT/DR(1)
      D(1,1) = DTDX*Q1(1)*(U(2,1)-U(1,1))
      D(1,2) = DTDX*Q1(1)*(U(2,2)-U(1,2))
      D(1,3) = DTDX*Q1(1)*(U(2,3)-U(1,3))
      D(1,4) = DTDX*Q1(1)*(U(2,4)-U(1,4))
C
      DO 150 J=1,JX-1
      U(J,1) = WK(J,1)+D(J,1)
      U(J,2) = WK(J,2)+D(J,2)
      U(J,3) = WK(J,3)+D(J,3)
      U(J,4) = WK(J,4)+D(J,4)
      es1(j) = u(j,4)/u(j,1)-((u(j,2)/u(j,1))*(u(j,2)/u(j,1))
     &      +(u(j,3)/u(j,1))*(u(j,3)/u(j,1)))/2.
150   CONTINUE

      do 240 J=2,JX-1
         dtdx = dt/dr(j)
         DES(J) = DTDX*(Q2(J)*(es1(J+1)-es1(J))
     &           -Q2(J-1)*(es1(J)-es1(J-1)))
240   CONTINUE
      DTDX  = DT/DR(1)
      DES(1)= DTDX*Q2(1)*(es1(2)-es1(1))

      DO 250 J=1,JX-1
         es1(J) = es1(J)+DES(J)
250   CONTINUE

      do 260 j=1,jx-1
      if (nl(j).ge.2) then
         es1(j) = oes1(j)
c
         u(j,1) = ou(j,1)
         u(j,2) = ou(j,2)
         u(j,3) = ou(j,3)
         u(j,4) = ou(j,4)
      endif
260   continue
C
      kcau2=0
      do 153 j=1,jx-1
         if (u(j,1).le.0.e14) then
            nl(j)=nl(j)+1
            kcau2=1
            nlmax=max(nlmax,nl(j))
            goto 153
         endif
c         es=u(j,4)/u(j,1)-((u(j,2)/u(j,1))*(u(j,2)/u(j,1))
c     &      +(u(j,3)/u(j,1))*(u(j,3)/u(j,1)))/2.
c         if (es.le.0.e14) then
         if (es1(j).le.0.e14) then
            nl(j)=nl(j)+1
            kcau2=1
c            print *,"@",ns
            nlmax=max(nlmax,nl(j))
         endif
153   continue
c------------------------------------------
                       if (kcau2.eq.1) then
c                          write(6,*) "kcau2=",kcau2
c                          pause
c------------------------------------------
      if (nlmax.le.10) then
         do 170 j=1,jx-1
            W0(J,1) = wt(j,1)
            W0(J,2) = wt(j,2)
            W0(J,3) = wt(j,3)
            W0(J,4) = wt(j,4)
            tauabso(j)=tauabst(j)
            tauo(j)=taut(j)
170      continue
         time=time-dt
         dt=dt*0.1
         time=time+dt
c
c         print *,"Check after 170 loop (nlmax>=10)"
c         print *,"time=",time,"dt=",dt
c         pause
c
         kback=1
         goto 151
      endif
         if (u(j,1).le.0.) then
            krs(j)=1
         endif
         if (es1(j).le.0.) then
            kes(j)=1
         endif
         call print
         stop
c---------------------------
                        else
c---------------------------
C
C
c      print *,"LW7"
      DO 160 J=1,JX-1
         R   = X1(J)
         R2  = R*R
         R3  = R2*R
         RM1 = R-1
         OK  = R/RM1*SQRT(0.5/R3)
         RS  = U(J,1)
         VX  = U(J,2)/RS
         VP  = U(J,3)/RS
         EE  = U(J,4)/RS
         SS  = RS/R
         ES  = EE-0.5*(VX*VX+VP*VP)
         ec2=es1(j)*cc2
         HHS = SQRT(B5*(AIN/AINP1)/SS)/OK
         beta=betao(j)
         te=tej(j)

c ----- start 158 iteration! -----
c@@@@@ iss > 10 for converge @@@@@      
      do 158 iss=1,10
c      if (j.eq.1) print *,beta,te
         obeta = beta
c
         A     = 3*(1-BETA)+BETA/GM1
         WW    = EC2*SS/(A+0.5)
         WC    = WW/CC2
         HH    = HHS*SQRT(WW)
         DE    = SS/(2*AIN*HH*RG)
         DEMN  = AIN*DE
         TEMN  = 2*TE/3.0
         XKFF  = XKKR*DEMN*(TEMN**(-3.5))
c         TAUFF = 2.*XKFF*DE*HH*RG
         TAUFF = XKFF*DE*HH*RG
c      if (tauff.le.0.) print *,"",tauff,de,hh
      if (tauff.le.0.) print *,"TAUABS=",tauff,"LW1D2"
c      TAUES  = 2.*XKES*DE*HH*RG
         TAUES = XKES*DE*HH*RG
         TAU   = TAUES+TAUFF
         qbr   = 2.0*XEFF*AJN*DE*DE*SQRT(TE)*HH*RG
         taup  = qbr/(8.*aa*cc*(te**4))
         PR    = WW/(2*AINP1*HH*RG)
      if (kcool.eq.1) then 
         te=ondo(de,pr,taup,tau)
      else if (kcool.eq.0) then 
         te=temp(de,pr) 
      end if
cc         te=temp(de,pr)
         BETA = RMU*DE*TE/PR
c***** 2002-7-30 *****
      if (abs((beta-obeta)/beta) .lt. 1.e-6) go to 731
158   continue
 731  continue
c ***** end of 158 iteration !*****
c      taueff=sqrt(taup*(tau+2./sqrt(3.)))
c      taueffjx=sqrt(tauabso(jx)*(tauo(jx)+2./sqrt(3.)))
c
         if (kcool .eq. 1) then 
            taueff = sqrt(taup*(tau+2./sqrt(3.)))
         else 
            taueff = sqrt(3*tau*tauff)
         end if
c
      if (taueff.gt.1.4618e30) then
         beta=betao(j)
c         tau=tauo(j)
c         taup=tauabso(j)
c         te=wt(j,4)
         RS  = ou(J,1)
         VX  = ou(J,2)/RS
         VP  = ou(J,3)/RS
         EE  = ou(J,4)/RS
         SS  = RS/R
         ec2 = oes1(j)*cc2
         A      = 3*(1-BETA)+BETA/GM1
         WW     = EC2*SS/(A+0.5)
         WC     = WW/CC2
         HH     = HHS*SQRT(WW)
         DE     = SS/(2*AIN*HH*RG)
         DEMN   = AIN*DE
         TEMN   = 2*TE/3.0
         XKFF   = XKKR*DEMN*(TEMN**(-3.5))
         TAUFF  = XKFF*DE*HH*RG
         TAUES  = XKES*DE*HH*RG
         TAU    = TAUES+TAUFF
         qbr    = 2.0*XEFF*AJN*DE*DE*SQRT(TE)*HH*RG
         taup   = qbr/(8.*aa*cc*(te**4))
         PR     = WW/(2*AINP1*HH*RG)
         te     = ondo(de,pr,taup,tau)
      end if
      WS    = WC/SS
c
      XK      = XKES+XKFF
        fd = 16*AA*(TE**4)/(3.*XK*DE*HH*CC2)

c--- 2002-04-11 ---
      if (kcool .eq.1) then 
        FM     = 8.*AA*(TE**4)/(1.5*TAU+SQRT(3.)+1./TAUP)*UNIT
      else
        FM     = FD
      end if
c
      W0(J,1) = DE
      W0(J,2) = VX
      W0(J,3) = VP
      W0(J,4) = TE
      tauo(j) = tau
      tauabso(j) = taup
      betao(j)   = beta
c      taueff     = sqrt(taup*(tau+2./sqrt(3.)))
c
         if (kcool .eq. 1) then 
            taueff = sqrt(taup*(tau+2./sqrt(3.)))
         else 
            taueff = sqrt(3*tau*tauff)
         end if
c
      oes1(j)    = es1(j)
      ou(j,1) = u(j,1)
      ou(j,2) = u(j,2)
      ou(j,3) = u(j,3)
      ou(j,4) = u(j,4)
cc      FFTR    = FM/UNIT/XMCRIT/CC*16.
      FFTR    = FM/UNIT/XMCRIT/CC
      XLumi   = XLumi+2.*PI*(R*RG)*(DR(J)*RG)*FFTR
      if (taueff.le.1.) xltn=xltn+2.*PI*(R*RG)*(DR(J)*RG)*FFTR
      if (taueff.gt.1.) xltk=xltk+2.*PI*(R*RG)*(DR(J)*RG)*FFTR
c
c      xmdeff = 2.*pi*r2*alpha*ww/(xl-xlin)/xmcrit*rg/cc
      xmdeff = -2.*pi*r*rg*vx*cc*ss/xmcrit
C      xmdeff = 2.*pi*r2*alpha*ww/xlk/xmcrit*rg/cc
c --- Data for S-shaped curve ---
      IF (MOD(NS,KLC).EQ.0) then 
         frag1=0
      else if (frag1 .eq. 1) then 
            go to 160
      else if (r .ge. 5.) then 
            frag1=1
            go to 901
c
 901  WRITE (10,910) TIME,r,ss,te,xmdeff
 910  FORMAT (1PE10.3,1X,1PE10.3,1X,1PE10.3,1X,1PE10.3,1x,1pe10.3,
     &        1x,1pe10.3)
      end if 
c ---
c
160   CONTINUE
      xlumi = xlumi+xljx
      xltk  = xltk+xljx
C
c##### Data for yy.lc (Light Curve) ##### 
c     time : Time 
c     xlumi: Total Disk luminosity 
c     xltn : Disk luminosity for [tau < 1] 
c     xttk : Disk luminosity for [tau > 1]
c     dt   : Delta t at that time
c     NS   : Time step
c     KLC  : Light curve => file9 every KLC time steps!
c########################################
c
      if (imore .eq. 2) then 
         if (mod(ns,200) .eq. 0) then 
            write(6,88) "ns=",ns," t=",time," dt=",dt, " L=",xlumi
 88         format (a3,i12,2x,a3,1pe9.3,2x,a3,1pe9.3,2x,a3,1pe9.3)
         end if
      end if
c
      IF (MOD(NS,KLC).EQ.0) WRITE (9,900) TIME,xlumi,xltn,xltk,dt,ns
      RETURN
c-----------------------
                   endif
c-----------------------
cc900   FORMAT (1PE16.8,1X,1PE15.7,1X,1PE15.7,1X,1PE15.7,1x,1pe15.7)
900   FORMAT (1PE10.4,2X,1PE10.4,2X,1PE10.4,2X,1PE10.4,2x,1pe10.4,2x,i9)

5     CONTINUE
         drp(1) = x1(2) -x1(1)
         drm(jx)= x1(jx)-x1(jx-1)
         dr(1)  = drp(1)
         dr(jx) = drm(jx)

      do 152 j=1,jx-1
         nl(j)=0
152   continue
c
151   continue
      xlumi=0.  ! Total luminosity
      xltn=0.   ! optically thin luminosity
      xltk=0.   ! optically thick luminosity

c      print *,W0(1,1)
c      print *,"L1"

      DO 10 j=1,jx
         r    = x1(j)
         rm1  = r-1
         r2   = r*r
         r3   = r2*r
         ok   = (r/rm1)*sqrt(0.5/r3)
         de   = w0(j,1)
         vx   = w0(j,2)
         vp   = w0(j,3)
         te   = w0(j,4)
         wt(j,1) = w0(j,1)
         wt(j,2) = w0(j,2)
         wt(j,3) = w0(j,3)
         wt(j,4) = w0(j,4)
         tauabst(j) = tauabso(j)
         taut(j)    = tauo(j)
         tau        = tauo(j)
         taup       = tauabso(j)
c
      if (j.eq.jx) then
         pr   = (aa/3.)*(te**4)+rmu*de*te
      else
         if (kcool.eq.1) fac = 1./((1.+1./taup/(1.5*tau+sqrt(3.))))
         if (kcool.eq.0) fac = 1.
c
         pr   = (aa/3.)*(te**4)*fac+rmu*de*te
      endif
c
         beta = rmu*de*te/pr
         hh   = sqrt(b5*pr/de)/ok
         ss   = 2*ain*de*hh*rg
         ww   = 2*ainp1*pr*hh*rg
         wc   = ww/cc2
         ws   = wc/ss
         a    = 3*(1-beta)+beta/gm1
         ee   = (a+0.5)*ws+0.5*(vx*vx+vp*vp)       ! total energy
         rs   = r*ss
         betaj(j)= beta

         u(j,1) = rs                              ! d(r*SS)/dt 
         u(j,2) = rs*vx                           ! d(r*SS*vr)/dt 
         u(j,3) = rs*vp                           ! d(r*SS*vp)/dt
         u(j,4) = rs*ee                           ! d(r*SS*ee)/dt

         f(j,1) = rs*vx                         ! d(r*SS*vr)/dr
         f(j,2) = rs*(vx*vx+ws)                 ! d(r*SS*vr^2+r*WW)/dr
         f(j,3) = rs*(vx*vp+alpha*ws)           ! d(r^2*SS*vr*vp-Trp)/dr
         f(j,4) = rs*(vx*(ee+ws)+alpha*ws*vp)  ! d()/dr

         vk     = r*ok
         xlk    = r*vk
         dlnor  =-(0.5+r/rm1)
         dlnlr  = 2+dlnor
         dokdr  = ok*dlnor/r
         dlkdr  = xlk*dlnlr/r
         demn   = ain*de
         temn   = 2*te/3.0
         xkff   = xkkr*demn*(temn**(-3.5))
         tauff  = xkff*de*hh*rg
         taues  = xkes*de*hh*rg
         tau    = taues+tauff
         qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
         taup   = qbr/(8.*aa*cc*(te**4))
c
         if (kcool .eq. 1) then 
            taueff = sqrt(taup*(tau+2./sqrt(3.)))
         else if (kcool .eq. 0) then 
            taueff = sqrt(3*tau*tauff)
         end if
c
         if (ns.eq.lns+1) then
c --------- oes1(j): total energy = E - 0.5*(vr^2*vp^2) 
            oes1(j)=u(j,4)/u(j,1)-((u(j,2)/u(j,1))*(u(j,2)/u(j,1))
     &           +(u(j,3)/u(j,1))*(u(j,3)/u(j,1)))/2.
            ou(j,1)=u(j,1)
            ou(j,2)=u(j,2)
            ou(j,3)=u(j,3)
            ou(j,4)=u(j,4)
            tauo(j)=tau
            tauabso(j)=taup
            betao(j)=beta
         endif

         xk     = xkes+xkff
         fd     = 16*aa*(te**4)/(3*xk*de*hh*cc2)
         ff     = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg*rg/cc2/cc
c--- 2002-04-11 ---
         if (kcool .eq.1) then 
           fm     = 8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)*unit
         else
           fm     = fd
         end if
c
         h(j,1) = 0
         h(j,2) = wc*(1-dlnor)+ss*(2*vp*vk+vp*vp)    ! (11.30)¤Î±¦ÊÕ
         h(j,3) =-ss*(vx*vp+vx*dlkdr+alpha*ws)       
         h(j,4) =-r*rs*(vx*vp+alpha*ws)*dokdr-r*fm   ! (11.32)¤Î±¦ÊÕ
         IF (ns.EQ.lns+1) THEN
cc             FFTR    = FM/UNIT/XMCRIT/CC*16.
             fftr    = fm/unit/xmcrit/cc
             xlumi   = xlumi+2.*pi*(r*rg)*(dr(j)*rg)*fftr
          if (j.eq.jx) xljx=2.*pi*(r*rg)*(dr(j)*rg)*fftr
          if (taueff.le.1.) xltn=xltn+2.*pi*(r*rg)*(dr(j)*rg)*fftr
          if (taueff.gt.1.) xltk=xltk+2.*pi*(r*rg)*(dr(j)*rg)*fftr
         ENDIF
         betac(j) = beta
         okc(j)   = ok
         hhc(j)   = hh
         prc(j)   = pr
         tec(j)   = te
         ssc(j)   = ss
         dec(j)   = de
10    CONTINUE
c      IF ((NS.EQ.LNS+1).and.(kcau1+kcau2.eq.0)) 
c     &    WRITE (9,900) TIME-DT,XLumi,xltn,xltk,dt
c      print *,"L1.5"

c###### Thermal Conduction #################################### 
c    xkei: thermal conductivity => K=alpha*xkappa*C_s*Sigma*C_v*H
c
      if (ktcond.eq.1) then
c         do 12 j=2,jx-1
         do 12 j=1,jx-1
c            ef=3.*xkappa*(8.-7.*betac(j))
CCCc            cv=3.*(8.-7.*betac(j))/(8.-3.*betac(j)-3.*betac(j))
            cv=3./2.
cc          xnu=alpha*okc(j)*(hhc(j)**2.)
cc          xkei=xkappa*ssc(j)*cv*xnu
c          gm31=gm1*(4.-3.*betac(j))
c     &         /(betac(j)+12.*gm1*(1.-betac(j)))
c         ef=xkappa*gm31
          ef=xkappa
          xkei=alpha*ef*okc(j)*(hhc(j)**3)*prc(j)/tec(j)
          cn=cv*xkei*(tec(j+1)-tec(j))/drp(j)
     &       *sqrt(2./xnpl1)*ainp1*unit
c         if (nl(j).ge.1) cn=0.
         f(j,4)=f(j,4)-x1(j)*cn
         h(j,4)=h(j,4)-cn
12     continue
      endif
C########################################## 
C
c      print *,"LW2"
      DO 40 j=1,jx-1
         ave1 = 0.5*(u(j+1,1)+u(j,1))
         ave2 = 0.5*(u(j+1,2)+u(j,2))
         ave3 = 0.5*(u(j+1,3)+u(j,3))
         ave4 = 0.5*(u(j+1,4)+u(j,4))
         fsa1 = f(j+1,1)-f(j,1)
         fsa2 = f(j+1,2)-f(j,2)
         fsa3 = f(j+1,3)-f(j,3)
         fsa4 = f(j+1,4)-f(j,4)
         swa1 = 0.5*(h(j+1,1)+h(j,1))
         swa2 = 0.5*(h(j+1,2)+h(j,2))
         swa3 = 0.5*(h(j+1,3)+h(j,3))
         swa4 = 0.5*(h(j+1,4)+h(j,4))
         dtdx = dt/drp(j)
         y(j,1) = ave1-0.5*dtdx*fsa1+0.5*dt*swa1
         y(j,2) = ave2-0.5*dtdx*fsa2+0.5*dt*swa2
         y(j,3) = ave3-0.5*dtdx*fsa3+0.5*dt*swa3
         y(j,4) = ave4-0.5*dtdx*fsa4+0.5*dt*swa4
         es1(j)=y(j,4)/y(j,1)-((y(j,2)/y(j,1))*(y(j,2)/y(j,1))
     &         +(y(j,3)/y(j,1))*(y(j,3)/y(j,1)))/2.
40    CONTINUE
         dtdx    = dt/dr(1)
         fsa1    = f(2,1)-f(1,1)
         fsa2    = f(2,2)-f(1,2)
         fsa3    = f(2,3)-f(1,3)
         fsa4    = f(2,4)-f(1,4)
         swa1    = 0.5*(h(1,1)+h(2,1))
         swa2    = 0.5*(h(1,2)+h(2,2))
         swa3    = 0.5*(h(1,3)+h(2,3))
         swa4    = 0.5*(h(1,4)+h(2,4))
         wk(1,1) = u(1,1)-dtdx*fsa1+dt*swa1
         wk(1,2) = u(1,2)-dtdx*fsa2+dt*swa2
         wk(1,3) = u(1,3)-dtdx*fsa3+dt*swa3
         wk(1,4) = u(1,4)-dtdx*fsa4+dt*swa4
C
c##### Artificial Viscosity #####
c What is kcau1 & kcau2 ? => Ans.: kcau1,2 = 0,or 1
c      if kcau1&2 = 0 -> no problem! 
c      if kcau1&2 = 1 -> problem! -> modify!

      if (kcau1+kcau2.eq.0) then
         qu1=qu
      else
         qu1=16. 
      endif

      DO 330 j = 1,jx-1
         vxp = w0(j+1,2)
         vxm = w0(j,2)
c        qu1=qu
c        if (nl(j).ge.1) qu1=qu*5.
         q2(j)= qu1*abs(vxp-vxm)
330   CONTINUE
c
      do 340 j=2,jx-1
         dtdx=dt/dr(j)
         des(j) = dtdx*(q2(j)*(es1(j+1)-es1(j))
     &           -q2(j-1)*(es1(j)-es1(j-1)))
340   CONTINUE
         dtdx = dt/dr(1)
         des(1)=dtdx*q2(1)*(es1(2)-es1(1))
c
      DO 350 j=1,jx-1
         es1(j) = es1(j)+des(j)
c         if (kcau1+kcau2.ne.0) es1(1)=oes1(1)
c         if (es1(j).le.0.) es1(j)=oes1(j)
350   CONTINUE

      do 360 j=1,jx-1
c      if ((y(1,1).le.0.).or.(es1(1).le.0.)) then
c      if ((y(j,1).le.0.).or.(es1(j).le.0.)) then
      if (nl(j).ge.2) then
         es1(j)=oes1(j)
         y(j,1)=ou(j,1)
         y(j,2)=ou(j,2)
         y(j,3)=ou(j,3)
         y(j,4)=ou(j,4)
      endif
360   continue
C
      kcau1=0
      do 553 j=1,jx-1
         if (y(j,1).le.0.e0) then  ! if de < 0, then kcau1=1
            nl(j)=nl(j)+1
            kcau1=1
            nlmax=max(nlmax,nl(j))
            goto 553
         endif
c         es=y(j,4)/y(j,1)-((y(j,2)/y(j,1))*(y(j,2)/y(j,1))
c     &      +(y(j,3)/y(j,1))*(y(j,3)/y(j,1)))/2.
c         if (es1(j).le.0.e14) then
         if (es1(j).le.0.e0) then
            nl(j)=nl(j)+1
            kcau1=1
c            print *,"!",ns
            nlmax=max(nlmax,nl(j))
         endif
553   continue
c------------------------------------------
                       if (kcau1.eq.1) then
c                          write(6,*) "kcau1=",kcau1
c                          pause
c------------------------------------------
      if (nlmax.le.10) then
         do 570 j=1,jx-1
            w0(j,1) = wt(j,1)
            w0(j,2) = wt(j,2)
            w0(j,3) = wt(j,3)
            w0(j,4) = wt(j,4)
            tauabso(j)=tauabst(j)
            tauo(j)=taut(j)
570      continue
         time=time-dt
         dt=dt*0.1
         time=time+dt
c
c         print *,"Check after 570 loop (nlmax < 10)"
c         print *,"time=",time,"dt=",dt
c         pause
c
         kback=1
         goto 151
      endif
         if (y(j,1).le.0.) then
            krs(j)=1
         endif
         if (es1(j).le.0.) then
            kes(j)=1
         endif
         call print
         stop
c---------------------------
                        else
c---------------------------
c      print *,"L3"
      DO 100 j=1,jx-1
         r   = 0.5*(x1(j)+x1(j+1))
         r2  = r*r
         r3  = r2*r
         rm1 = r-1
         ok  = r/rm1*sqrt(0.5/r3)
         rs  = y(j,1)
         vx  = y(j,2)/rs
         vp  = y(j,3)/rs
         ee  = y(j,4)/rs
         ss  = rs/r
         es  = ee-0.5*(vx*vx+vp*vp)
         ec2 = es1(j)*cc2
         hhs = sqrt(b5*(ain/ainp1)/ss)/ok
         beta= betao(j)
         te  = w0(j,4)

c ***** start 98 iteration! *****
c@@@@@ iss > 10 for converge (best value=10) @@@@@
      do 98 iss=1,20
c      if (j.eq.1) print *,beta,te
         obeta  = beta
c  
         a      = 3*(1-beta)+beta/gm1
         ww     = ec2*ss/(a+0.5)
         wc     = ww/cc2
         hh     = hhs*sqrt(ww)
         de     = ss/(2*ain*hh*rg)
         demn   = ain*de
         temn   = 2*te/3.0
         xkff   = xkkr*demn*(temn**(-3.5))
         tauff  = xkff*de*hh*rg
c         TAUFF  = 2.*XKFF*DE*HH*RG
c         TAUES  = 2.*XKES*DE*HH*RG
         taues  = xkes*de*hh*rg
         tau    = taues+tauff
         qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
         taup   = qbr/(8.*aa*cc*(te**4))
         pr     = ww/(2*ainp1*hh*rg)
      if (kcool.eq.1) then 
         te = ondo(de,pr,taup,tau)
      else if (kcool.eq.0) then 
         te = temp(de,pr)
      end if
         beta   = rmu*de*te/pr

c***** 2002-7-30 *****
         if (abs((beta-obeta)/beta) .lt. 1.e-6) go to 730
 98   continue
 730     continue
c ***** end 98 iteration! *****
c      taueff = sqrt(taup*(tau+2./sqrt(3.)))
c      taueffjx=sqrt(tauabso(jx)*(tauo(jx)+2./sqrt(3.)))
c
         if (kcool .eq. 1) then 
            taueff = sqrt(taup*(tau+2./sqrt(3.)))
         else 
            taueff = sqrt(3*tau*tauff)
         end if
c
      if (taueff.gt.1.4618e30) then
       beta  = betao(j)
c       tau  = tauo(j)
c       taup = tauabso(j)
c       te   = wt(j,4)
         rs  = ou(j,1)
         vx  = ou(j,2)/rs
         vp  = ou(j,3)/rs
         ee  = ou(j,4)/rs
         ss  = rs/r
         ec2 = oes1(j)*cc2
         hhs = sqrt(b5*(ain/ainp1)/ss)/ok
         a   = 3*(1-beta)+beta/gm1
         ww  = ec2*ss/(a+0.5)
         wc  = ww/cc2
         hh  = hhs*sqrt(ww)
         de  = ss/(2*ain*hh*rg)
         demn   = ain*de
         temn   = 2*te/3.0
         xkff   = xkkr*demn*(temn**(-3.5))
         tauff  = xkff*de*hh*rg
         taues  = xkes*de*hh*rg
         tau    = taues+tauff
         qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
         taup   = qbr/(8.*aa*cc*(te**4))
         pr     = ww/(2*ainp1*hh*rg)
         te     = ondo(de,pr,taup,tau)
      end if

      ws     = wc/ss
      f(j,1) = rs*vx
      f(j,2) = rs*(vx*vx+ws)
      f(j,3) = rs*(vx*vp+alpha*ws)
      f(j,4) = rs*(vx*(ee+ws)+alpha*ws*vp)
      vk     = r*ok
      xlk    = r*vk
      dlnor  = -(0.5+r/rm1)
      dlnlr  = 2+dlnor
      dokdr  = ok*dlnor/r
      dlkdr  = xlk*dlnlr/r

      xk     = xkes+xkff
      fd     = 16*aa*(te**4)/(3*xk*de*hh*cc2)
c--- 2002-04-11 ---
      if (kcool .eq.1) then 
        fm     = 8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)*unit
      else
        fm     = fd
      end if

      h(j,1) = 0
      h(j,2) = wc*(1-dlnor)+ss*(2*vp*vk+vp*vp)
      h(j,3) = -ss*(vx*vp+vx*dlkdr+alpha*ws)
      h(j,4) = -r*rs*(vx*vp+alpha*ws)*dokdr-r*fm
      betaj(j)= beta
      betac(j)= beta
      okc(j) = ok
      hhc(j) = hh
      prc(j) = pr
      tej(j) = te
      tec(j) = te
      ssc(j) = ss
      dec(j) = de
100   CONTINUE
c ----- End of L3 loop -----

c###### Thermal Conduction #################################### 
      if (ktcond.eq.1) then
c      do 103 j=2,jx-1
      do 103 j=1,jx-1
c         ef=3.*xkappa*(8.-7.*betac(j))
CCCc         cv=3.*(8.-7.*betac(j))/(8.-3.*betac(j)-3.*betac(j))
         cv=3./2.
cc         xnu=alpha*okc(j)*(hhc(j)**2.)
cc         xkei=xkappa*ssc(j)*cv*xnu
c         gm31=gm1*(4.-3.*betac(j))
c     &        /(betac(j)+12.*gm1*(1.-betac(j)))
c         ef = xkappa*gm31
         ef = xkappa
         xkei = alpha*ef*okc(j)*(hhc(j)**3)*prc(j)/tec(j)
         cn = cv*xkei*(tec(j+1)-tec(j))/drp(j)
     &      *sqrt(2./xnpl1)*ainp1*unit
c         if (nl(j).ge.1) cn=0.
         f(j,4)=f(j,4)-x1(j)*cn
         h(j,4)=h(j,4)-cn
103   continue
      end if 
c##########################################
c---------------------------------- 
                             end if
C----------------------------------
c      print *,"L4"
      DO 120 j=2,jx-1
         dtdx   = dt/dr(j)
         fsa1   = f(j,1)-f(j-1,1)
         fsa2   = f(j,2)-f(j-1,2)
         fsa3   = f(j,3)-f(j-1,3)
         fsa4   = f(j,4)-f(j-1,4)
         swa1   = 0.5*(h(j,1)+h(j-1,1))
         swa2   = 0.5*(h(j,2)+h(j-1,2))
         swa3   = 0.5*(h(j,3)+h(j-1,3))
         swa4   = 0.5*(h(j,4)+h(j-1,4))
         wk(j,1) = u(j,1)-dtdx*fsa1+dt*swa1
         wk(j,2) = u(j,2)-dtdx*fsa2+dt*swa2
         wk(j,3) = u(j,3)-dtdx*fsa3+dt*swa3
         wk(j,4) = u(j,4)-dtdx*fsa4+dt*swa4
c      endif
120   CONTINUE
C
C****** ARTIFICIAL VISCOSITY *****
C
c      print *,"L6"

      if (kcau1+kcau2.eq.0) then
         qb1=qb
      else
         qb1=16.
      endif

      DO 130 j = 1,jx-1
         vxp  = w0(j+1,2)
         vxm  = w0(j,2)
         q1(j)= qb1*abs(vxp-vxm)
130   CONTINUE
C
      DO 140 j=2,jx-1
      dtdx   = dt/dr(j)
      d(j,1) = dtdx*(q1(j)*(u(j+1,1)-u(j,1))-q1(j-1)*(u(j,1)-u(j-1,1)))
      d(j,2) = dtdx*(q1(j)*(u(j+1,2)-u(j,2))-q1(j-1)*(u(j,2)-u(j-1,2)))
      d(j,3) = dtdx*(q1(j)*(u(j+1,3)-u(j,3))-q1(j-1)*(u(j,3)-u(j-1,3)))
      d(j,4) = dtdx*(q1(j)*(u(j+1,4)-u(j,4))-q1(j-1)*(u(j,4)-u(j-1,4)))
140   CONTINUE
C
      dtdx   = dt/dr(1)
      d(1,1) = dtdx*q1(1)*(u(2,1)-u(1,1))
      d(1,2) = dtdx*q1(1)*(u(2,2)-u(1,2))
      d(1,3) = dtdx*q1(1)*(u(2,3)-u(1,3))
      d(1,4) = dtdx*q1(1)*(u(2,4)-u(1,4))
C
      DO 150 j=1,jx-1
      u(j,1) = wk(j,1)+d(j,1)
      u(j,2) = wk(j,2)+d(j,2)
      u(j,3) = wk(j,3)+d(j,3)
      u(j,4) = wk(j,4)+d(j,4)
      es1(j) = u(j,4)/u(j,1)-((u(j,2)/u(j,1))*(u(j,2)/u(j,1))
     &      +(u(j,3)/u(j,1))*(u(j,3)/u(j,1)))/2.
150   CONTINUE

      do 240 j=2,jx-1
         dtdx = dt/dr(j)
         des(j) = dtdx*(q2(j)*(es1(j+1)-es1(j))
     &           -q2(j-1)*(es1(j)-es1(j-1)))
240   CONTINUE
      dtdx  = dt/dr(1)
      des(1)= dtdx*q2(1)*(es1(2)-es1(1))

      DO 250 j=1,jx-1
         es1(j) = es1(j)+des(j)
250   CONTINUE

      do 260 j=1,jx-1
      if (nl(j).ge.2) then
         es1(j) = oes1(j)
c
         u(j,1) = ou(j,1)
         u(j,2) = ou(j,2)
         u(j,3) = ou(j,3)
         u(j,4) = ou(j,4)
      endif
260   continue
C
      kcau2=0
      do 153 j=1,jx-1
         if (u(j,1).le.0.e14) then
            nl(j)=nl(j)+1
            kcau2=1
            nlmax=max(nlmax,nl(j))
            goto 153
         endif
c         es=u(j,4)/u(j,1)-((u(j,2)/u(j,1))*(u(j,2)/u(j,1))
c     &      +(u(j,3)/u(j,1))*(u(j,3)/u(j,1)))/2.
c         if (es.le.0.e14) then
         if (es1(j).le.0.e14) then
            nl(j)=nl(j)+1
            kcau2=1
c            print *,"@",ns
            nlmax=max(nlmax,nl(j))
         endif
153   continue
c------------------------------------------
                       if (kcau2.eq.1) then
c                          write(6,*) "kcau2=",kcau2
c                          pause
c------------------------------------------
      if (nlmax.le.10) then
         do 170 j=1,jx-1
            w0(j,1) = wt(j,1)
            w0(j,2) = wt(j,2)
            w0(j,3) = wt(j,3)
            w0(j,4) = wt(j,4)
            tauabso(j)=tauabst(j)
            tauo(j)=taut(j)
170      continue
         time=time-dt
         dt=dt*0.1
         time=time+dt
c
c         print *,"Check after 170 loop (nlmax>=10)"
c         print *,"time=",time,"dt=",dt
c         pause
c
         kback=1
         goto 151
      endif
         if (u(j,1).le.0.) then
            krs(j)=1
         endif
         if (es1(j).le.0.) then
            kes(j)=1
         endif
         call print
         stop
c---------------------------
                        else
c---------------------------
C
C
c      print *,"LW7"
      DO 160 j=1,jx-1
         r   = x1(j)
         r2  = r*r
         r3  = r2*r
         rm1 = r-1
         ok  = r/rm1*sqrt(0.5/r3)
         rs  = u(j,1)
         vx  = u(j,2)/rs
         vp  = u(j,3)/rs
         ee  = u(j,4)/rs
         ss  = rs/r
         es  = ee-0.5*(vx*vx+vp*vp)
         ec2=es1(j)*cc2
         hhs = sqrt(b5*(ain/ainp1)/ss)/ok
         beta=betao(j)
         te=tej(j)

c ----- start 158 iteration! -----
c@@@@@ iss > 10 for converge @@@@@      
      do 158 iss=1,10
c      if (j.eq.1) print *,beta,te
         obeta = beta
c
         a     = 3*(1-beta)+beta/gm1
         ww    = ec2*ss/(a+0.5)
         wc    = ww/cc2
         hh    = hhs*sqrt(ww)
         de    = ss/(2*ain*hh*rg)
         demn  = ain*de
         temn  = 2*te/3.0
         xkff  = xkkr*demn*(temn**(-3.5))
c         TAUFF = 2.*XKFF*DE*HH*RG
         tauff = xkff*de*hh*rg
c      if (tauff.le.0.) print *,"",tauff,de,hh
      if (tauff.le.0.) print *,"TAUABS=",tauff,"LW1D2"
c      TAUES  = 2.*XKES*DE*HH*RG
         taues = xkes*de*hh*rg
         tau   = taues+tauff
         qbr   = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
         taup  = qbr/(8.*aa*cc*(te**4))
         pr    = ww/(2*ainp1*hh*rg)
      if (kcool.eq.1) then 
         te=ondo(de,pr,taup,tau)
      else if (kcool.eq.0) then 
         te=temp(de,pr) 
      end if
cc         te=temp(de,pr)
         beta = rmu*de*te/pr
c***** 2002-7-30 *****
      if (abs((beta-obeta)/beta) .lt. 1.e-6) go to 731
158   continue
 731  continue
c ***** end of 158 iteration !*****
c      taueff=sqrt(taup*(tau+2./sqrt(3.)))
c      taueffjx=sqrt(tauabso(jx)*(tauo(jx)+2./sqrt(3.)))
c
         if (kcool .eq. 1) then 
            taueff = sqrt(taup*(tau+2./sqrt(3.)))
         else 
            taueff = sqrt(3*tau*tauff)
         end if
c
      if (taueff.gt.1.4618e30) then
         beta=betao(j)
c         tau=tauo(j)
c         taup=tauabso(j)
c         te=wt(j,4)
         rs  = ou(j,1)
         vx  = ou(j,2)/rs
         vp  = ou(j,3)/rs
         ee  = ou(j,4)/rs
         ss  = rs/r
         ec2 = oes1(j)*cc2
         a      = 3*(1-beta)+beta/gm1
         ww     = ec2*ss/(a+0.5)
         wc     = ww/cc2
         hh     = hhs*sqrt(ww)
         de     = ss/(2*ain*hh*rg)
         demn   = ain*de
         temn   = 2*te/3.0
         xkff   = xkkr*demn*(temn**(-3.5))
         tauff  = xkff*de*hh*rg
         taues  = xkes*de*hh*rg
         tau    = taues+tauff
         qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
         taup   = qbr/(8.*aa*cc*(te**4))
         pr     = ww/(2*ainp1*hh*rg)
         te     = ondo(de,pr,taup,tau)
      end if
      ws    = wc/ss
c
      xk      = xkes+xkff
        fd = 16*aa*(te**4)/(3.*xk*de*hh*cc2)

c--- 2002-04-11 ---
      if (kcool .eq.1) then 
        fm     = 8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)*unit
      else
        fm     = fd
      end if
c
      w0(j,1) = de
      w0(j,2) = vx
      w0(j,3) = vp
      w0(j,4) = te
      tauo(j) = tau
      tauabso(j) = taup
      betao(j)   = beta
c      taueff     = sqrt(taup*(tau+2./sqrt(3.)))
c
         if (kcool .eq. 1) then 
            taueff = sqrt(taup*(tau+2./sqrt(3.)))
         else 
            taueff = sqrt(3*tau*tauff)
         end if
c
      oes1(j)    = es1(j)
      ou(j,1) = u(j,1)
      ou(j,2) = u(j,2)
      ou(j,3) = u(j,3)
      ou(j,4) = u(j,4)
cc      FFTR    = FM/UNIT/XMCRIT/CC*16.
      fftr    = fm/unit/xmcrit/cc
      xlumi   = xlumi+2.*pi*(r*rg)*(dr(j)*rg)*fftr
      if (taueff.le.1.) xltn=xltn+2.*pi*(r*rg)*(dr(j)*rg)*fftr
      if (taueff.gt.1.) xltk=xltk+2.*pi*(r*rg)*(dr(j)*rg)*fftr
c
c      xmdeff = 2.*pi*r2*alpha*ww/(xl-xlin)/xmcrit*rg/cc
      xmdeff = -2.*pi*r*rg*vx*cc*ss/xmcrit
C      xmdeff = 2.*pi*r2*alpha*ww/xlk/xmcrit*rg/cc
c --- Data for S-shaped curve ---
      IF (mod(ns,klc).EQ.0) then 
         frag1=0
      else if (frag1 .eq. 1) then 
            go to 160
      else if (r .ge. 5.) then 
            frag1=1
            go to 901
c
 901  WRITE (10,910) time,r,ss,te,xmdeff
 910  FORMAT (1pe10.3,1x,1pe10.3,1x,1pe10.3,1x,1pe10.3,1x,1pe10.3,
     &        1x,1pe10.3)
      end if 
c ---
c
160   CONTINUE
      xlumi = xlumi+xljx
      xltk  = xltk+xljx
C
c##### Data for yy.lc (Light Curve) ##### 
c     time : Time 
c     xlumi: Total Disk luminosity 
c     xltn : Disk luminosity for [tau < 1] 
c     xttk : Disk luminosity for [tau > 1]
c     dt   : Delta t at that time
c     NS   : Time step
c     KLC  : Light curve => file9 every KLC time steps!
c########################################
c
      if (imore .eq. 2) then 
         if (mod(ns,200) .eq. 0) then 
            write(6,88) "ns=",ns," t=",time," dt=",dt, " L=",xlumi
 88         format (a3,i12,2x,a3,1pe9.3,2x,a3,1pe9.3,2x,a3,1pe9.3)
         end if
      end if
c
      IF (mod(ns,klc).EQ.0) WRITE (9,900) time,xlumi,xltn,xltk,dt,ns
      RETURN
c-----------------------
                   endif
c-----------------------
cc900   FORMAT (1PE16.8,1X,1PE15.7,1X,1PE15.7,1X,1PE15.7,1x,1pe15.7)
900   FORMAT (1pe10.4,2x,1pe10.4,2x,1pe10.4,2x,1pe10.4,2x,1pe10.4,2x,i9)

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function ondo	(		DE,
			PR,
			taup,
			tau
	)

Definition at line 1420 of file accdisk.f.

C******************************
C
      include 'para.h'
C
      a1  = 3.*rmu*de/aa*(1.+1./taup/(1.5*tau+sqrt(3.)))
      a2  =-3.*pr/aa*(1.+1./taup/(1.5*tau+sqrt(3.)))
      eps = 1.0e-12
      tmax= (-a2)**0.25
      tmin= 0.
      fx  = tmax**4+a1*tmax+a2
      fn  = tmin**4+a1*tmin+a2
      DO 40 itr=1,10
         th = 0.5*(tmax+tmin)
         fh = th**4+a1*th+a2
         IF (fh*fn .LT. 0.) tmax=th
         IF (fh*fn .LT. 0.) fx=fh
         IF (fh*fx .LT. 0.) tmin=th
         IF (fh*fx .LT. 0.) fn=fh
 40   CONTINUE
C
      xxx = th
C
      DO 100 itr=1,20
         xx1 = xxx-(xxx**4+a1*xxx+a2)/(4.*xxx**3+a1)
         IF (abs((xx1-xxx)/xx1) .LT. eps) GO TO 120
         xxx = xx1
 100  CONTINUE
C
      WRITE (6,110)
      print *,de,pr,tmax,tmin,fx,fn,th,fh,xxx,xx1
 110  FORMAT (1h ,'LOOP COUNT (TEMP) IS OVER MAXIMUM')
      stop
C
 120  CONTINUE
C
      ondo=xx1
C
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine print

(

)

Definition at line 2615 of file accdisk.f.

C**************************
C
      include 'para.h'
C     
      dimension drp(jmax),drm(jmax),dr(jmax)
C     
      if ((krs(j).eq.1).or.(kes(j).eq.1)) goto 100
      IF (mod(ns,iprint(4)).NE.0) RETURN
      js = iprint(2)
      IF (iprint(3).NE.0) THEN
         IF (mod(ns,iprint(3)).EQ.0) js = 1
      END IF
 100  continue
C     
c      print *,"PRINT"
      gm1 = gm-1
      cc2 = cc*cc
C     
      DO 5 j=2,jx-1
         drp(j) = x1(j+1)-x1(j)
         drm(j) = x1(j)-x1(j-1)
         dr(j)  = 0.5*(drp(j)+drm(j))
 5    CONTINUE
      drp(1) = x1(2) -x1(1)
      drm(jx)= x1(jx)-x1(jx-1)
      dr(1)  = drp(1)
      dr(jx) = drm(jx)
      
      WRITE(6,*)
      WRITE (6,1000) ns,time,dt,nlmax

      IF (imore.EQ.1) THEN
         WRITE (6,1002)
      ELSE IF (imore.EQ.2) THEN
         WRITE (6,1004)
      ELSE
         WRITE (6,1005)
      ENDIF

      DO 10 j=1,jx,js
         r    = x1(j)
         rm1  = r-1
         r2   = r*r
         r3   = r2*r
         ok   = (r/rm1)*sqrt(0.5/r3)
         de   = w0(j,1)
         vx   = w0(j,2)
         vp   = w0(j,3)
         te   = w0(j,4)
         if (ki.eq.1) then
            if (j.eq.jx) then
               pr   = (aa/3.)*(te**4)+rmu*de*te
            else            
               if (kcool .eq.1) then 
                  fac =  1./((1.+1./tauabso(j)/(1.5*tauo(j)+sqrt(3.))))
               else 
                  fac = 1. 
               end if
               pr   = (aa/3.)*(te**4)*fac+rmu*de*te
            endif
         else if (kr.eq.1) then
            pr   = (aa/3.)*(te**4)+rmu*de*te
         else
            if (j.eq.jx) then
               pr   = (aa/3.)*(te**4)+rmu*de*te
            else
               pr   = rmu*de*te
            endif
         endif
c--------------------------------------------------------
c     BETA : (gas pressure)/(total pressure)
c     HH   : Disk Half Thickness 
c     SS   : Surface Density 
c     WW   : Height Integrated Pressure
c     WC   : WW/c^2 
c     WS   : WC/SS 
c     
         beta = rmu*de*te/pr
         hh   = sqrt(b5*pr/de)/ok
         ss   = 2*ain*de*hh*rg
         ww   = 2*ainp1*pr*hh*rg
         wc   = ww/cc2
         ws   = wc/ss
c------------------------------------------ 
      a1   = 4.-3.*beta
      a2   = beta+12.*gm1*(1.-beta)
      gg1  = beta+(a1*a1*gm1)/a2
      gg3  = 1.+(gg1-beta)/a1

      cs2  = ((3.*gg1-1.)+2.*(gg3-1.)*alpha*alpha)*ws/(gg1+1.)
      cs   = sqrt(cs2)
      cfln = dt/dr(j)*(cs+abs(vx))
      vk   = ok*r
      xlk  = ok*r2
      xl   = r*(vp+vk)
c      xmdeff = 2.*pi*r2*alpha*ww/(xl-xlin)/xmcrit*rg/cc
c      xmdeff = 2.*pi*r2*alpha*ww/xlk/xmcrit*rg/cc
      xmdeff = -2.*pi*r*rg*vx*cc*ss/xmcrit
      xmd    = -2.*pi*r*ss*vx*rg*cc/xmcrit
c
      demn   = ain*de
      temn   = 2*te/3.0
      xkff   = xkkr*demn*(temn**(-3.5))
      tauff  = xkff*de*hh*rg
      taues  = xkes*de*hh*rg
      tau    = taues+tauff

      qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
      taup   = qbr/(8.*aa*cc*(te**4))
c
      taueff = sqrt(tauabso(j)*(tauo(j)+2./sqrt(3.)))
c
         if (kcool .eq. 1) then 
            taueff = sqrt(tauabso(j)*(tauo(j)+2./sqrt(3.)))
         else 
            taueff = sqrt(3*tau*tauff)
         end if
c
      t1 = tauo(j)+2./sqrt(3.)
      t2 = 1./tauabso(j)
c      if (t1.ge.t2) then
c         taueff=1.
c      else
c         taueff=2.
c      endif 
      a    = 3*(1-beta)+beta/gm1
      ee   = (a+0.5)*ws+0.5*(vx*vx+vp*vp)
c      FMC  = 8.*AA*CC*(TE**4.)/(3.0*tauo(j))
      fmc   = 8.*aa*cc*(te**4)/(1.5*tauo(j)+sqrt(3.)+1./taup)

      teff = (4.*fmc/aa/cc)**0.25 
      vphi = vp+vk

      xll = xlk-xlin
      qp    = 2.*alpha*alpha*xl/xll
c      if ((kr.eq.1).or.(r.eq.rout1)) then
c         QM    = FF1*BETA4*XLL7*ETA5/(XKXE*R4)
c      else
      eta   = vx*vx/(b2*pr/de)
      ff    = qbr*rg/cc/cc/cc
      qm    = 2.*pi*r*r*eta/(xmdot*ws)*ff
c      endif
      qadv  = qp-qm
      fene  = qadv/qp
c      write(6,*) "check!",j,r,qp,qm,fene,b2,xmdot,ws
c      pause
c
c     ----- For print routine -----   
      IF (imore.EQ.1) THEN
         WRITE(6,1010) j,r,dr(j),de,vx,vp,te,cs,cfln,ww,ss,xmd,hh,es1(j)
     *        ,ws,beta,tauo(j)
      ELSE IF (imore.EQ.2) THEN
         IF (ns.EQ.lns) THEN
            WRITE(6,1010) j,r,beta,taueff,te,ss,xmdeff,xmd,t1,1/t2
         ELSE
            WRITE(6,1010) j,r,beta,taueff,te,ss,xmdeff,xmd,t1,1/t2
         ENDIF
      else if (imore .eq. 3) then 
         write(6,1010) j,r,dr(j),de,te,teff,ss,vx,vp,cs,hh,tauo(j)
     *        ,taueff,ww,xmd      
      ELSE if (imore .eq. 0) then 
         IF (ns.EQ.lns) THEN
            WRITE(6,1010) j,r,beta,teff,te,ss,-vx,vphi,cs,hh,tauo(j) 
     *           ,taueff,ww,fene,xmd
         ELSE
            WRITE(6,1010) j,r,beta,teff,te,ss,-vx,vphi,cs,hh,tauo(j)
     *           ,taueff,ww,fene,xmd
         ENDIF
      ENDIF

 10   CONTINUE
      do 20 j=1,jx-1
         if (krs(j).ne.0) then
            write (6,1030) j,ns+1,time
        krs(j)=0
         end if
 20   continue
      do 21 j=1,jx-1
         if (kes(j).ne.0) then
            write (6,1031) j,ns+1,time
            kes(j)=0
         end if
21    continue

      nlmax =0
C
      RETURN
C
 1000 FORMAT ('NS=',i12,4x,'TIME=',f12.3,4x,'DT=',e10.3,4x,
     &     "NLMAX=",i5//)
 1002 format (4x,"J",5x,"R",11x,"DR",10x,"DE",10x,"VX",
     &     10x,"VP",10x,"TE",10x,"CS",9x,"CFLN",9x,"WW",10x,"SS",
     &     9x,"MDOT",9x,"HH",10x,"ES1",9x,"WS",9x,"BETA",9x,"TAU")
c     1002   format (4x,"J",5x,"R",11x,"DR",10x,"DE",10x,"VX",
c     &        10x,"VP",10x,"TE",10x,"CS",9x,"WW",10x,"SS",
c     &        9x,"MDOT",9x,"HH",10x,"BETA",9x,"TAU")
 1004 format (4x,"J",5x,"R",10x,"BETA",8x,"TAUEFF",7x,"TE",
     &     10x,"SS",8x,"MDOTEFF",7x,"MDOT")
c     
c     ¤³¤Î¹Ô¤Ï2000/08/07¤ËÄɲà 1005 format (4x,"J",5x,"R",10x,"BETA",8x,"TEFF",8x,"TE",
     &     10x,"SS",10x,"VX",10x,"VP",10x,"CS",10x,"HH",
     &     10x,"tauo",8x,"taueff",6x,"WW",10x,"WS",10x,"MDOT")
c     1005  format (4x,"J",5x,"R",10x,"BETA",9x,"TEFF",8x,"TE",
c     &        10x,"SS",9x,"MDOT")
 1010 FORMAT (I5,20(1PE12.4))
 1020 FORMAT (I5,5(1PE12.4),I8)
 1030 format ("CAUTION   RS < 0  AT  J =",I3,"  NS =",I8,
     &     "  TIME =",E10.3)
 1031 format ("CAUTION   ES < 0  AT  J =",I3,"  NS =",I8,
     &     "  TIME =",E10.3)
C     

c     imore=0¤Ç¤Îoutput data!
 1005 format (4x,"J",5x,"R",10x,"BETA",8x,"TEFF",8x,"TE",
     &     10x,"SS",10x,"VX",10x,"VP",10x,"CS",10x,"HH",
     &     10x,"tauo",8x,"taueff",6x,"WW",10x,"WS",10x,"MDOT")
c     1005  format (4x,"J",5x,"R",10x,"BETA",9x,"TEFF",8x,"TE",
c     &        10x,"SS",9x,"MDOT")
 1010 FORMAT (i5,20(1pe12.4))
 1020 FORMAT (i5,5(1pe12.4),i8)
 1030 format ("CAUTION   RS < 0  AT  J =",i3,"  NS =",i8,
     &     "  TIME =",e10.3)
 1031 format ("CAUTION   ES < 0  AT  J =",i3,"  NS =",i8,
     &     "  TIME =",e10.3)
C     

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine prtrb1

(

)

Definition at line 2495 of file accdisk.f.

C**************************
      include 'para.h'

      xmdot=pxmdot
      xmd=pxmd
      xlin=pxlin
      do 10 j=1,jx-1
         de   = w0(j,1)
         vx   = w0(j,2)   
         te   = w0(j,4)
         r    = x1(j)
         rm1  = r-1
         r2   = r*r
         r3   = r2*r
         ok   = (r/rm1)*sqrt(0.5/r3)
         pr   = (aa/3.)*(te**4)+rmu*de*te
         hh   = sqrt(b5*pr/de)/ok
         ss   = 2*ain*de*hh*rg
         vx   = -xmdot/(2.*pi*r*ss)/cc
         w0(j,2)=vx
 10   continue
      return

Here is the caller graph for this function:

subroutine prtrb2

(

)

Definition at line 2520 of file accdisk.f.

C**************************
      include 'para.h'
      xmdot=pxmdot
      xmd=pxmd
c      xlin=pxlin
c      rc = x1(jx)
c      rp = 10.0

      do 10 j=jx,jx
c         CALL SADMK(X1(J),ETA,DL,W0(J,1),W0(J,2),W0(J,3),W0(J,4))
c         CALL SADMK2(X1(J),ETA,DL,W0(J,1),W0(J,2),W0(J,3),W0(J,4))
         CALL sadmn2(x1(j),eta,dl,w0(j,1),w0(j,2),w0(j,3),w0(j,4))
c         do 20 k=1,1
c            w0(j,k)=w0(j,k)
c     &      +rp*w0(j,k)*exp(-(x1(j)-rc)*(x1(j)-rc))
c 20      continue

 10   continue
      RETURN

Here is the call graph for this function:

subroutine prtrb3

(

)

Definition at line 2542 of file accdisk.f.

C**************************
      include 'para.h'
      dimension ww(4)
c      xmdot=pxmdot
c      xmd=pxmd
c      xlin=pxlin
c      CALL SADMK(X1(JX),ETA,DL,ww(1),ww(2),ww(3),ww(4))
c-----------------------
c      rc     : Radius 
c      xlamda : wave length
c      rp     : amplitude of the wave 
c-----------------------
c     rc=x1(JX)
      rc=50.0
c      xlamda=pxlin
      xlamda=1.0
      rp=1.0
      do 10 j=1,jx-1
         do 20 k=1,4
            w0(j,k)=w0(j,k)
     &      +rp*w0(j,k)*exp(-(x1(j)-rc)*(x1(j)-rc)/xlamda/xlamda)
c     &           -rp*w0(j,k)*exp(-(x1(j)-rc)*(x1(j)-rc)/xlamda/xlamda)
 20      continue
 10   continue
c      write(6,*) "PERTRB 3"
      return

subroutine prtrb4

(

)

Definition at line 2572 of file accdisk.f.

C**************************
      include 'para.h'
      dimension ww(4)
c      xmdot=pxmdot
c      xmd=pxmd
c-----------------------
c      rc     : Radius 
c      xlamda : wave length
c      rp     : amplitude of the wave 
c-----------------------
      rc=30.0
c     xlamda=pxlin
      xkk=30
      rp=2.0

      do 10 j=1,jx-1
         de   = w0(j,1)
         vx   = w0(j,2)   
         te   = w0(j,4)
         r    = x1(j)
         rm1  = r-1
         r2   = r*r
         r3   = r2*r
         ok   = (r/rm1)*sqrt(0.5/r3)
         pr   = (aa/3.)*(te**4)+rmu*de*te
         hh   = sqrt(b5*pr/de)/ok
         ss   = 2*ain*de*hh*rg
         vx   = -xmdot/(2.*pi*r*ss)/cc
         w0(j,2)=vx

c         do 20 k=1,1
c            dss = 0.01*exp(-(x1(j)-rc)**2/0.1**2)*sin(xkk*x1(j))*w0(j,k)
c            w0(j,k)=w0(j,k)+dss
 20      continue
         dww = dss/xkk/x1(j)*w0(j,1)
         dvr = dss/sqrt(xkk*x1(j))*w0(j,2)
         dvp = dss/sqrt(xkk*x1(j))*w0(j,3)
 10   continue
      write(6,*) "PERTRB 4"
      return

subroutine pzextr	(	integer	iest,
			xest,
		dimension(nv)	yest,
		dimension(nv)	yz,
		dimension(nv)	dy,
		integer	nv
	)

Definition at line 3132 of file accdisk.f.

      implicit real*8 (a-h,o-z)
      common /lscal/ kpr
      integer iest,nv,imax,nmax,kpr
c      real xest,dy(nv),yest(nv),yz(nv)
      dimension dy(nv),yest(nv),yz(nv)
      parameter(imax=13,nmax=50)
      integer j,k1
c      real delta,f1,f2,q,d(nmax),qcol(nmax,imax),x(imax)
      dimension d(nmax),qcol(nmax,imax),x(imax)
      save qcol,x
c --- Save current independent variable.
      x(iest)=xest
      do 11 j=1,nv
         dy(j)=yest(j)
         yz(j)=yest(j)
 11       continue
c --- Store first estimate in first column.
      if (iest.eq.1) then
          do 12 j=1,nv
             qcol(j,1)=yest(j)
 12       continue
      else
          do 13 j=1,nv
             d(j)=yest(j)
 13       continue
          do 15 k1=1,iest-1
             delta=1./(x(iest-k1)-xest)
             f1=xest*delta
             f2=x(iest-k1)*delta   
c --- propagate tableau 1 diagonal more.
             do 14 j=1,nv
                q=qcol(j,k1)
                qcol(j,k1)=dy(j)
                delta=d(j)-q
                dy(j)=f1*delta
                d(j)=f2*delta
                yz(j)=yz(j)+dy(j)
 14          continue
 15       continue
          do 16 j=1,nv
             qcol(j,iest)=dy(j)
 16       continue
      endif
      return

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine rdfl7

(

)

Definition at line 2877 of file accdisk.f.

************************************************************************
      include 'para.h'
      
      read (7,802) jx
      read (7,800)
      do 10 j=1,jx
         read (7,803) x1(j)
 10   continue
 100  continue
C     
      read (7,800)
      read (7,804) ns,time
      read (7,800)
      do 30 j=1,jx
         read (7,805) (w0(j,i),i=1,4),tauabso(j),tauo(j)
         if (ns.eq.0) then
            deini(j)=w0(j,1)
            vxini(j)=w0(j,2)
            vpini(j)=w0(j,3)
            teini(j)=w0(j,4)
         endif
 30   continue
      read (7,806) dt
      read (7,806) xmd
      read (7,801) xlin
      if (ns.eq.0) xmdi=xmd
C     
      call file
c     call print
      if (ns.ne.lns) goto 100
      xmdot=xmd*xmcrit
 800  format (80x)
 801  format (13x,f18.14)
 802  format (5x,i4)
 803  format (4x,1pe15.7)
 804  format (4x,i12,8x,f12.3)
 805  format (7x,6(1pe22.14),i3)
 806  format (13x,e11.3)
      return

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine rdprm

(

)

Definition at line 91 of file accdisk.f.

C*************************

      include 'para.h'

C***** CONSTANTS *****
      aa    = 7.56464e-15       ! Radiation constant (erg/cm^3/K^4)
      cc    = 2.997925e10       ! Speed of light
      gg    = 6.672e-8          ! Gravitational constant
      xmu   = 0.5               ! 1/2
      xmh   = 1.660531e-24      ! Atomic mass (g)
      xkb   = 1.38062e-16       ! Boltzman constant (erg/K)
      rmu   = xkb/(xmu*xmh)     ! 2 xkb/xmh (coefficient of gas pres) 
      xmsun = 1.989e33          ! Solar mass (g)
      xkes  = 0.34           ! Electron Scattering Opacity:
c     XKKR  = 0.34*6E24      ! Free-free opacity coefficient
c      xkes  = 0.4
      xkkr  = 0.64e23
      xeff  = 6.22e20
c
c----- Ho-shi's parameter(1977) ----------------------------------------
      npl   = 3                 ! Polytropic index (gamma=1+1/N)
      xnpl1 = npl+1             ! XNPL1:N+1
      ain   = fnin(npl)         ! I_N -> related to heght integration. 
      ainp1 = fnin(npl+1)       ! AIN   : I_N
      ainp1 = fnin(npl+1)       ! AINP1 : I_N+1
      ajn   = fnjn(2*npl)       ! AJN   : I_2N 
c----- Specific heat ? -----
      gm    = 5./3.
      gm1   = gm-1.
      pi    = 3.14159263
C     
C     ***** READ PARAMETERS *****
C     
      open (9,file="at.job",status="old") ! Data input 
      read (9,*)
      read (9,*)
      read (9,901) date
      read (9,*)
      read (9,*) bhm,alpha,xmd,xlinm,xlinp
      read (9,*)
      read (9,*) mesh,rin,rin2,rout1,rout2
      read (9,*)
      read (9,*) dt,dtmax,dtmin,nsx,lns
      read (9,*)
      read (9,*) (iprint(i),i=1,4),ksfl,irstrt
      read (9,*)
      read (9,*) hini,relerr,abserr,qb,qu
      read (9,*)
      read (9,902) ki,file7,file8,imore,xkappa
      read (9,*)
      read (9,903) klc,file9,file10
      read (9,*)
      read (9,*) kr,pxmd,pxlin,kcool
      read (9,*) 
      read (9,*) ktcond
      read (9,*)
      read (9,*) cfc 

 901  format (4x,a10)
 902  format (4x,i2,7x,a6,4x,a6,4x,i2,8x,e12.5)
 903  format (4x,i6,3x,a6,4x,a8)
C
C ***** WRITE INPUT PARAMETERS *****
C
      WRITE (6,*) '***** INPUT PARAMETERS *****', date
      WRITE (6,100) bhm,alpha,xmd,xlinm,xlinp
      WRITE (6,200) mesh,rin,rout1,rout2
      WRITE (6,300) dt,dtmax,dtmin,nsx
c     WRITE (6,400) (IPRINT(I),I=1,4),KSFL
      WRITE (6,400) qb,qu,xkappa
      write (6,500) hini,relerr,abserr
C
 100  FORMAT (/1h ,'BHM  =',1pe15.7,5x,'ALPHA=',1pe15.7,5x,'XMDOT=',
     *     1pe15.7/1h ,'XLINM=',1pe15.7,5x,'XLINP=',1pe15.7)
 200  FORMAT (1h ,'MESH =',i15,5x,'RIN  ='
     *     ,1pe15.7/1h ,'ROUT1=',1pe15.7,5x,'ROUT2=',1pe15.7)
 300  FORMAT (1h ,'DT   =',1pe15.7,5x,'DTMAX=',1pe15.7,5x,'DTMIN=',
     *     1pe15.7/1h ,'NSX  =',i15)
c     400   FORMAT (1H ,'IPRINT=',4(I6),5X,'KSFL =',I8//)
 400  FORMAT (1h ,'QB=',1pe15.7,5x,"QU=",1pe15.7,5x,'XKAPPA=',
     *     1pe15.7,//)
 500  format (1h ,"HINI=",1pe9.1,5x,"RELERR=",1pe9.1,5x
     &     ,"ABSERR=",1pe9.1//) 
C     
c--- B1,B2,B3,B4,B5 -> vertical integration factor (for normalize). ---
c
      rg     = 2.*gg*bhm*xmsun/(cc*cc) !  Schwarzschild radius (cm) 
c      XMCRIT = 32.*PI*CC*RG/XKES
      xmcrit = 2.*pi*cc*rg/xkes ! Critical accretion rate (L_E/eta/c^2)
      xmdi   = xmd              ! normalized accretion rate
      xmdot  = xmd*xmcrit              
      pxmdot = pxmd*xmcrit
      al7    = alpha**7
      ai5    = (ain**5)/(ainp1**4)
      ff1    = (128.*pi*pi*aa*cc/(3.*xkes*xmdot*xmdot*al7))*ai5
     *     *((cc/rmu)**4)*((cc*rg)**3)
      b1 = ain*cc*rg*rg
      b2 = (ainp1/ain)/(cc*cc)
      b3 = ainp1*rg*rg/cc
      b4 = rg/(cc*cc)
      b5 = 2.*xnpl1/(cc*cc)

      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sadmk	(		R,
			ETA,
			DL
	)

Definition at line 266 of file accdisk.f.

C******************************************
C
C -------------------------------------------------------------------
C     STANDARD ACCRETION DISK MODEL
C     OPTICALLY THICK, GEOMETRICALLY THIN, KEPLERIAN, ALPHA-MODEL
C     PSEUDO NEWTONIAN POTENTIAL
C
C     ETA = (VR**2)/(WW/SS)     DL = XL-XLK
C --------------------------------------------------------------------
C
      include 'para.h'
C
      dimension ph(4)
C
c      print *,"SADMK"
      xkes1 = xkes
      xkkr1 = xkkr

      rm1   = r-1
      r2    = r*r
      r3    = r2*r
      eps   = 1.0e-6

      ok    = (r/rm1)*sqrt(.5/r3) ! OK : Kepler angular velocity 
      vk    = r*ok              ! VK : Kepler velocity
      xlk   = r*vk              ! XLK : Kepler angular momentum  
      dlkdr = xlk*(1.5/r - 1./rm1) !DLKDR : dlnLk/dr 
      dlnok = -1./rm1 - .5/r    !DLNOK : d(lnOK)/dlnR
      ww    = xmdot*(xlk-xlin)/(2.*pi*b3*r*r*alpha)
      ff2   = (-3.*xkes1*b3*r*alpha*ww*ok*dlnok/(32.*b4*aa*cc))*ww
      ff3   = ok*ok*ww*ww/(4.*b5)
      ff4   = (aa/3.)*(ff2**0.9)/(sqrt(ff3)*(rmu**0.4))
      betan = 0.01
      betax = 1.00
      DO 10 i=1,10
         beta  = 0.5*(betan+betax)
         fff   = 1.-beta-ff4*(beta**0.4)
         IF (fff .GT. 0.) betan=beta
         IF (fff .LE. 0.) betax=beta
10    CONTINUE
C
c      print *,"S1"
      pr = (3./aa)*(ff3/ff2)*(1.-beta)! Pressure (at equatrial plane)
      de = (pr**3)/ff3          !DE : density (g/cm^3)
      te = (de*ff2/pr)**0.25    !TE : temperature (K)
      hh = (sqrt(b5*pr/de))/ok  !HH : Scale height    
      vr = xmdot/(4.*pi*b1*r*de*hh) !VR : radial velocity 
c      write(6,*) pr,de,te
c      write(6,*) hh,vr
C
C ***** CORRECTION FOR OPACITY *****
C
C    --- ITERATION ---
c ----------------------------------------------------------------------
c     PR   : a/3 * TE^4 (radiation pressure)+2kb/mh*de*TE (gas pressure)
c     BETA : P_gas (Gas pressure)/P_total (Total pressure)
c     DEMN : I_N DE (mean density)
c     TEMN : 2/3*TE (mean temperature)
c     XKFF : Free-Free opacity
c     XK   : Mean Opacity
c     XKXF : free-free opacity/mean opacity
c     HH   : Scale Height 
c     FM   : Radiative flux
c     FF3  :
c     FF4  : 
c ----------------------------------------------------------------------
c      print *,"S2"
      DO 200 itr=1,20
         pr   = (aa/3.)*(te**4)+rmu*de*te
c     if (pr.le.0.) print *,"PR=",pr,itr
         beta = rmu*de*te/pr
         demn = ain*de
         temn = 2.*te/3.
         xkff = xkkr1*demn*(temn**(-3.5))
         xk   = xkes1+xkff
c     if (xk.le.0.) print *,"XK=",xk
         xkxf = xkff/xk
         hh   = (sqrt(b5*pr/de))/ok
c     
c      write(6,*) itr,pr,beta
c
         if (de.le.0.) print *,"DE1=",de
         if (hh.le.0.) print *,"HH=",hh

         fm   = 2.*pi*b4*r*16.*aa*cc*(te**4)/(3.*xk*de*hh)
         ff3  = 4.*pi*b3*r*r*alpha*pr*hh
         ff4  = -xmdot*(xlk-xlin)*ok*dlnok
C     
         ph(3)= xmdot*(xlk-xlin)-ff3
         ph(4)= ff4-fm
C     
         a1 = 1.-beta
         a4 = 4.-3.*beta
C
c     if (te.le.0.) print *,"TE=",te
         dlnhdd = -0.5*a1/de
         dlnhdt = 0.5*a4/te
         dlnpdd = beta/de
         dlnpdt = a4/te
         dlnkdd = xkxf/de
         dlnkdt = -3.5*xkxf/te
         dlnfdd = -dlnkdd-dlnhdd-1./de
         dlnfdt = 4./te-dlnkdt-dlnhdt
C
c ---------------------------------------------
c         | D3DD D3DT |
c    det. |           | == D3DD*D4DT-D3DT*D4DD
c         | D4DD D4DT |
c ---------------------------------------------
C   ** JACOBIANS **
C
         d3dd = -ff3*(dlnpdd+dlnhdd)
         d3dt = -ff3*(dlnpdt+dlnhdt)
         d4dd = -fm*dlnfdd
         d4dt = -fm*dlnfdt
         det  = d3dd*d4dt-d3dt*d4dd
C
c --------------------------------------------
c    |D3DD D3DT| |DELD| = |PH(3)|
c    |D4DD D4DT| |DELT|   |PH(4)|
c
c    |DELD|    1  |D4DT -D3DT| |PH(3)|  
c    |    | = --- |          |*|     |
c    |DELT|   det.|-D4DD D3DD| |PH(4)|
c --------------------------------------------
C   ** CORRECTIONS **
C
c      if (det.le.0.) print *,"DET=",det
         deld = (ph(3)*d4dt-ph(4)*d3dt)/det
         delt = (d3dd*ph(4)-d4dd*ph(3))/det
C
C   ** APPLY CORRECTIONS **
C
         de = de-deld
         te = te-delt
c     WRITE (6,40) ITR,DE,TE
 40      FORMAT (1h ,i7,5x,'DE =',1pe12.4,5x,'TE =',1pe12.4)
C
C   ** CHECK CONVERGENCE **
C
c       if (de.le.0.) print *,"DE2=",de
c       if (te.le.0.) print *,"TE=",te
         IF (abs(deld/de) .GE. eps) GO TO 100
         IF (abs(delt/te) .GE. eps) GO TO 100
C     
         GO TO 300
C     
 100     CONTINUE
 200  CONTINUE
C
      if (imore.eq.2) WRITE (6,210)
 210  FORMAT (1h ,'LOOP COUNT (SADM) IS OVER MAXIMUM !')
      stop
C
 300  CONTINUE
C
c----- Physical Quantities at ROUT (in SSD). -----
c     
      pr    = (aa/3.)*(te**4)+rmu*de*te
      hh    = (sqrt(b5*pr/de))/ok
c      if (hh.le.0.) print *,"HH=",hh
      vr    = xmdot/(4.*pi*b1*r*hh*de)
      eta   = vr*vr/(b2*pr/de)
      ws    = ((xlk-xlin)**2)*eta/(r2*alpha*alpha)
c      if (xlk-xlin.eq.0.) print *,"XLK-XLIN",xlk-xlin
      dlnwr = -2./r+dlkdr/(xlk-xlin)
      dl    = r3*ws*(dlnwr+dlnok)/(2.*xlk)
C
      vx =-vr                   ! radial velocity
      vp = dl/r                 ! rotational velocity
      demn  = ain*de            ! Mean density
      temn  = 2*te/3.0          ! Mean temperature
c      TAUES  = 2.*XKES*DE*HH*RG
      taues = xkes*de*hh*rg     ! Optical depth (electron scattering)
      xkff  = xkkr*demn*(temn**(-3.5)) ! Mean Free-Free opacity
c      TAUFF = 2.*XKFF*DE*HH*RG
      tauff = xkff*de*hh*rg     ! Optical depth (free-free absorption)
c      if (tauff.le.0) print *,"TAUABS=",tauff,"LW1D4"
      qbr   = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg ! Brems cooling 
      taup  = qbr/(8.*aa*cc*(te**4)) ! Optical depth (planck)
      tau   = taues+tauff       ! Total optical depth 
      tauabso(jx)=taup 
      tauo(jx)=tau
      ss    = 2*ain*de*hh*rg    ! Surface density
c      WRITE(6,1010) TIME,SS,TE,XMD
C
      IF (iprint(1) .EQ. 0) RETURN
C
c
      if (imore.eq.2)  WRITE (6,400) de,vr,dl,te,eta
400   FORMAT (/1h ,'***** STANDARD MODEL *****'//1h ,'DE   =',1pe15.7,
     *            5x,'VR   =',1pe15.7,5x,'DL   =',1pe15.7/1h ,'TE   =',
     *               1pe15.7,5x,'ETA  =',1pe15.7//)
1010  FORMAT (20(1pe12.4))
C
c      print *,"S-END"
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine sadmk2	(		R,
			ETA,
			DL,
			DE,
			VX,
			VP,
			TE
	)

Definition at line 467 of file accdisk.f.

C******************************************
C
      include 'para.h'
C
      dimension ph(4)
C
c      print *,"SADMK"
      xkes1 = xkes
      xkkr1 = xkkr

      rm1   = r-1
      r2    = r*r
      r3    = r2*r
      eps   = 1.0e-6

      ok    = (r/rm1)*sqrt(.5/r3) ! OK : Kepler angular velocity 
      vk    = r*ok              ! VK : Kepler velocity
      xlk   = r*vk              ! XLK : Kepler angular momentum  
      dlkdr = xlk*(1.5/r - 1./rm1) !DLKDR : dlnLk/dr 
      dlnok = -1./rm1 - .5/r    !DLNOK : d(lnOK)/dlnR
      ww    = xmdot*(xlk-xlin)/(2.*pi*b3*r*r*alpha)
      ff2   = (-3.*xkes1*b3*r*alpha*ww*ok*dlnok/(32.*b4*aa*cc))*ww
      ff3   = ok*ok*ww*ww/(4.*b5)
      ff4   = (aa/3.)*(ff2**0.9)/(sqrt(ff3)*(rmu**0.4))
      betan = 0.01
      betax = 1.00
      DO 10 i=1,10
         beta  = 0.5*(betan+betax)
         fff   = 1.-beta-ff4*(beta**0.4)
         IF (fff .GT. 0.) betan=beta
         IF (fff .LE. 0.) betax=beta
10    CONTINUE
C
c      print *,"S1"
      pr = (3./aa)*(ff3/ff2)*(1.-beta)! Pressure (at equatrial plane)
      de = (pr**3)/ff3          !DE : density (g/cm^3)
      te = (de*ff2/pr)**0.25    !TE : temperature (K)
      hh = (sqrt(b5*pr/de))/ok  !HH : Scale height    
      vr = xmdot/(4.*pi*b1*r*de*hh) !VR : radial velocity 
c      print *,"S2"
      DO 200 itr=1,20
         pr   = (aa/3.)*(te**4)+rmu*de*te
c     if (pr.le.0.) print *,"PR=",pr,itr
         beta = rmu*de*te/pr
         demn = ain*de
         temn = 2.*te/3.
         xkff = xkkr1*demn*(temn**(-3.5))
         xk   = xkes1+xkff
c     if (xk.le.0.) print *,"XK=",xk
         xkxf = xkff/xk
         hh   = (sqrt(b5*pr/de))/ok
c     
c      write(6,*) itr,pr,beta
c
         if (de.le.0.) print *,"DE1=",de
         if (hh.le.0.) print *,"HH=",hh

         fm   = 2.*pi*b4*r*16.*aa*cc*(te**4)/(3.*xk*de*hh)
         ff3  = 4.*pi*b3*r*r*alpha*pr*hh
         ff4  = -xmdot*(xlk-xlin)*ok*dlnok
C     
         ph(3)= xmdot*(xlk-xlin)-ff3
         ph(4)= ff4-fm
C     
         a1 = 1.-beta
         a4 = 4.-3.*beta
C
c     if (te.le.0.) print *,"TE=",te
         dlnhdd = -0.5*a1/de
         dlnhdt = 0.5*a4/te
         dlnpdd = beta/de
         dlnpdt = a4/te
         dlnkdd = xkxf/de
         dlnkdt = -3.5*xkxf/te
         dlnfdd = -dlnkdd-dlnhdd-1./de
         dlnfdt = 4./te-dlnkdt-dlnhdt
C
c ---------------------------------------------
c         | D3DD D3DT |
c    det. |           | == D3DD*D4DT-D3DT*D4DD
c         | D4DD D4DT |
c ---------------------------------------------
C   ** JACOBIANS **
C
         d3dd = -ff3*(dlnpdd+dlnhdd)
         d3dt = -ff3*(dlnpdt+dlnhdt)
         d4dd = -fm*dlnfdd
         d4dt = -fm*dlnfdt
         det  = d3dd*d4dt-d3dt*d4dd
C
C   ** CORRECTIONS **
C
c      if (det.le.0.) print *,"DET=",det
         deld = (ph(3)*d4dt-ph(4)*d3dt)/det
         delt = (d3dd*ph(4)-d4dd*ph(3))/det
C
C   ** APPLY CORRECTIONS **
C
         de = de-deld
         te = te-delt
c     WRITE (6,40) ITR,DE,TE
 40      FORMAT (1h ,i7,5x,'DE =',1pe12.4,5x,'TE =',1pe12.4)
C
C   ** CHECK CONVERGENCE **
C
c       if (de.le.0.) print *,"DE2=",de
c       if (te.le.0.) print *,"TE=",te
         IF (abs(deld/de) .GE. eps) GO TO 100
         IF (abs(delt/te) .GE. eps) GO TO 100
C     
         GO TO 300
C     
 100     CONTINUE
 200  CONTINUE
C
      if (imore.eq.2) WRITE (6,210)
 210  FORMAT (1h ,'LOOP COUNT (SADM) IS OVER MAXIMUM !')
      stop
C
 300  CONTINUE
C
c----- Physical Quantities at ROUT (in SSD). -----
c     
      pr    = (aa/3.)*(te**4)+rmu*de*te
      hh    = (sqrt(b5*pr/de))/ok
c      if (hh.le.0.) print *,"HH=",hh
      vr    = xmdot/(4.*pi*b1*r*hh*de)
      eta   = vr*vr/(b2*pr/de)
      ws    = ((xlk-xlin)**2)*eta/(r2*alpha*alpha)
c      if (xlk-xlin.eq.0.) print *,"XLK-XLIN",xlk-xlin
      dlnwr = -2./r+dlkdr/(xlk-xlin)
      dl    = r3*ws*(dlnwr+dlnok)/(2.*xlk)
C
      vx =-vr                   ! radial velocity
      vp = dl/r                 ! rotational velocity
      demn  = ain*de            ! Mean density
      temn  = 2*te/3.0          ! Mean temperature
c      TAUES  = 2.*XKES*DE*HH*RG
      taues = xkes*de*hh*rg     ! Optical depth (electron scattering)
      xkff  = xkkr*demn*(temn**(-3.5)) ! Mean Free-Free opacity
c      TAUFF = 2.*XKFF*DE*HH*RG
      tauff = xkff*de*hh*rg     ! Optical depth (free-free absorption)
c      if (tauff.le.0) print *,"TAUABS=",tauff,"LW1D4"
      qbr   = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg ! Brems cooling 
      taup  = qbr/(8.*aa*cc*(te**4)) ! Optical depth (planck)
      tau   = taues+tauff       ! Total optical depth 
      tauabso(jx)=taup 
      tauo(jx)=tau
      ss    = 2*ain*de*hh*rg    ! Surface density
c      WRITE(6,1010) TIME,SS,TE,XMD
C
      IF (iprint(1) .EQ. 0) RETURN
C
c      print *,"S-END"
      RETURN

Here is the call graph for this function:

subroutine sadmn	(		R,
			ETA,
			DL
	)

Definition at line 626 of file accdisk.f.

C******************************************
c 2001/02/07 comments 
C     This subroutine is not used ! (except subroutine INTMDL)
C -------------------------------------------------------------------
C     STANDARD ACCRETION DISK MODEL
C     OPTICALLY THIN, GEOMETRICALLY THIN, KEPLERIAN, ALPHA-MODEL
C     PSEUDO NEWTONIAN POTENTIAL
C
C     ETA = (VR**2)/(WW/SS)     DL = XL-XLK
C
C --------------------------------------------------------------------
C
      include 'para.h'
C
      rm1 = r-1.
      r2  = r*r
      r3  = r2*r
      ok  =(r/rm1)*sqrt(.5/r3)
      vk  = r*ok
      xlk = r*vk
      xll = xlk-xlin
      dlnok =-1./rm1-.5/r
      dokdr = ok*dlnok
      cc2 = cc*cc
      rg2 = rg*rg
      ww  =(xmdot*xll)/(2.*pi*r2*alpha)*cc/rg
      wwtrue = ww
      t1=1.e4
      t2=1.e30
c100   continue
c ---------------------------------------
c      te=1.e6
      do 1 iss=1,10
      te = sqrt(t1*t2)
      wweps = 1.e-6
cc      print *,"TE=",te
      pd  =  rmu*te
      hh  =  sqrt(2.*xnpl1*pd)/ok/cc
      pr  =  ww/(2.*ainp1*hh)/rg
      de  =  pr/pd
      vr  =  xmdot/(4.*pi*r*ain*de*hh)/(cc*rg2)
      demn   = ain*de
      temn   = 2*te/3.0
      xkff   = xkkr*demn*(temn**(-3.5))
      tauff  = xkff*de*hh*rg
      taues  = xkes*de*hh*rg
      tau    = taues+tauff
      qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
      taup   = qbr/(8.*aa*cc*(te**4))
      wwfalse= -8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)
     *        /(alpha*r*dokdr)*(rg/cc)
      te4=-ww*(1.5*tau+sqrt(3.)+1./taup)*(alpha*r*dokdr)/(8.*aa)
     *                                              /(rg/cc)
cc      wwfalse=-8.*AA*(TE**4)*TAUP/(ALPHA*R*DOKDR)*(RG/CC)
cc      te4=-ww/TAUP*(ALPHA*R*DOKDR)/(8.*AA)/(RG/CC)
      te=te4**0.25
      if (wwfalse.gt.wwtrue) t1=te
      if (wwfalse.lt.wwtrue) t2=te
      if (abs(wwfalse-wwtrue)/wwtrue.lt.wweps) goto 100
c      if (abs(wwfalse-wwtrue).ge.1.e4) goto 100
1     continue
C --------------------------------------------
100   continue
      eta = ain*de*vr*vr/(ainp1*pr)*cc*cc
      dl=0.
      ss  = 2*ain*de*hh*rg
C -----------------------------
c      if (imore.eq.2) print *,"PR=",pr," HH=",hh," WW=",ww
c      if (imore.eq.2) WRITE (6,400) DE,VR,DL,TE,SS
      demn   = ain*de
      temn   = 2*te/3.0
      xkff   = xkkr*demn*(temn**(-3.5))
      tauff  = xkff*de*hh*rg
      taues  = xkes*de*hh*rg
      tau    = taues+tauff
      qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
      taup   = qbr/(8.*aa*cc*(te**4))
      taueff = sqrt(taup*(tau+2./sqrt(3.)))
      beta=1.
c      t1=tau+2./sqrt(3.)
c      t2=1./taup
c      if (t1.ge.t2) then
c         taueff=1.
c      else
c         taueff=2.
c      endif
c      WRITE(6,1000) TIME,SS,TE,XMD,taueff,tau,taup
cc      WRITE(6,1010) JX,R,BETA,TAUEFF,TE,SS,XMD
400   FORMAT (/1h ,'***** STANDARD MODEL *****'//1h ,'DE   =',1pe15.7,
     *            5x,'VR   =',1pe15.7,5x,'DL   =',1pe15.7/1h ,'TE   =',
     *               1pe15.7,5x,'SS  =',1pe15.7//)
1000  FORMAT (20(1pe12.4))
1010  FORMAT (i5,20(1pe12.4))
C
C
      RETURN

Here is the caller graph for this function:

subroutine sadmn2	(		R,
			ETA,
			DL,
			DE,
			VX,
			VP,
			TE
	)

Definition at line 725 of file accdisk.f.

C******************************************
c     for optically thin accretion disk
C
      include 'para.h'
C
      rm1 = r-1.
      r2  = r*r
      r3  = r2*r
      ok  =(r/rm1)*sqrt(.5/r3)
      vk  = r*ok
      xlk = r*vk
      xll = xlk-xlin
      dlnok =-1./rm1-.5/r
      dokdr = ok*dlnok
      cc2 = cc*cc
      rg2 = rg*rg
      ww  =(xmdot*xll)/(2.*pi*r2*alpha)*cc/rg
      wwtrue = ww
      t1=1.e4
      t2=1.e30
c100   continue
c ---------------------------------------
c      te=1.e6
      do 1 iss=1,10
      te = sqrt(t1*t2)
      wweps = 1.e-6
cc      print *,"TE=",te
      pd  =  rmu*te
      hh  =  sqrt(2.*xnpl1*pd)/ok/cc
      pr  =  ww/(2.*ainp1*hh)/rg
      de  =  pr/pd
      vr  =  xmdot/(4.*pi*r*ain*de*hh)/(cc*rg2)
      demn   = ain*de
      temn   = 2*te/3.0
      xkff   = xkkr*demn*(temn**(-3.5))
      tauff  = xkff*de*hh*rg
      taues  = xkes*de*hh*rg
      tau    = taues+tauff
      qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
      taup   = qbr/(8.*aa*cc*(te**4))
      wwfalse= -8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)
     *        /(alpha*r*dokdr)*(rg/cc)
      te4=-ww*(1.5*tau+sqrt(3.)+1./taup)*(alpha*r*dokdr)/(8.*aa)
     *                                              /(rg/cc)
cc      wwfalse=-8.*AA*(TE**4)*TAUP/(ALPHA*R*DOKDR)*(RG/CC)
cc      te4=-ww/TAUP*(ALPHA*R*DOKDR)/(8.*AA)/(RG/CC)
      te=te4**0.25
      if (wwfalse.gt.wwtrue) t1=te
      if (wwfalse.lt.wwtrue) t2=te
      if (abs(wwfalse-wwtrue)/wwtrue.lt.wweps) goto 100
c      if (abs(wwfalse-wwtrue).ge.1.e4) goto 100
1     continue
C --------------------------------------------
100   continue
      eta = ain*de*vr*vr/(ainp1*pr)*cc*cc
      dl  = 0.
      ss  = 2*ain*de*hh*rg
C -----------------------------
c      if (imore.eq.2) print *,"PR=",pr," HH=",hh," WW=",ww
c      if (imore.eq.2) WRITE (6,400) DE,VR,DL,TE,SS
      demn   = ain*de
      temn   = 2*te/3.0
      xkff   = xkkr*demn*(temn**(-3.5))
      tauff  = xkff*de*hh*rg
      taues  = xkes*de*hh*rg
      tau    = taues+tauff
      qbr    = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
      taup   = qbr/(8.*aa*cc*(te**4))
      taueff = sqrt(taup*(tau+2./sqrt(3.)))
      beta=1.

      vx = -vr
      vp = dl/r
C
      RETURN

Here is the caller graph for this function:

subroutine simpr	(	dimension(n)	y,
		dimension(n)	dydx,
		dimension(n)	dfdx,
		dimension(nmax,nmax)	dfdy,
		integer	nmax,
		integer	n,
			xs,
			htot,
		integer	nstep,
		dimension(n)	yout,
		external	derivs
	)

Definition at line 3064 of file accdisk.f.

      implicit real*8 (a-h,o-z)
      common /lscal/ kpr
      integer n,nmax,nstep,nmaxx,kpr
c      real htot,xs,dfdx(n),dfdy(nmax,nmax),dydx(n),y(n),yout(n)
      dimension dfdx(n),dfdy(nmax,nmax),dydx(n),y(n),yout(n)
      external derivs
c --- Maximum expected value of n.
      parameter(nmaxx=50)
      integer i,j,nn,indx(nmaxx)
c      real d,h,x,a(nmaxx,nmaxx),del(nmaxx),ytemp(nmaxx)
      dimension a(nmaxx,nmaxx),del(nmaxx),ytemp(nmaxx)
c --- Stepsize 
      h=htot/nstep
c --- Set up the matrix 1-hf'.
      do 12 i=1,n
         do 11 j=1,n
            a(i,j)=-h*dfdy(i,j)
 11      continue   
         a(i,i)=a(i,i)+1.
 12   continue
c --- LU decommposition of the matrix.
      call ludcmp (a,n,nmaxx,indx,d)
c --- Set up right hand side for first step. 
c --- Use yout for temporary storage.
      do 13 i=1,n
         yout(i)=h*(dydx(i)+h*dfdx(i))
 13   continue
      call lubksb (a,n,nmaxx,indx,yout)
c --- First step.
      do 14 i=1,n
         del(i)=yout(i)
         ytemp(i)=y(i)+del(i)
 14   continue
      x=xs+h
c --- Use yout for temporary storage of derivatives.
      call derivs(x,ytemp,yout)
      if (kpr.eq.1) return
c --- General step.
      do 17 nn=2,nstep
c --- set up right hand side for general step.
         do 15 i=1,n
            yout(i)=h*yout(i)-del(i)
 15      continue
         call lubksb (a,n,nmaxx,indx,yout)
         do 16 i=1,n
            del(i)=del(i)+2.*yout(i)
            ytemp(i)=ytemp(i)+del(i)
 16      continue
         x=x+h
         call derivs (x,ytemp,yout)
         if (kpr.eq.1) return
 17   continue
c --- Set up right hand side for last step.
      do 18 i=1,n
         yout(i)=h*yout(i)-del(i)
 18   continue
      call lubksb(a,n,nmaxx,indx,yout)
c --- take last step.
      do 19 i=1,n
         yout(i)=ytemp(i)+yout(i)
 19   continue
      kpr=0
      return

Here is the call graph for this function:

Here is the caller graph for this function:

subroutine stifbs	(	dimension(nv)	y,
		dimension(nv)	dydx,
		integer	nv,
			x,
			htry,
			eps,
		dimension(nv)	yscal,
			hdid,
			hnext,
		external	derivs,
			icon
	)

Definition at line 2922 of file accdisk.f.

      implicit real*8 (a-h,o-z)
      common /lscal/ kpr
      integer nv,nmax,kmaxx,imax,kpr
c      real eps,hdid,hnext,htry,x,dydx(nv),y(nv),yscal(nv),safe1,
c     &     safe2,redmax,redmin,tiny,scalmx
      dimension dydx(nv),y(nv),yscal(nv)
      external derivs
      parameter(nmax=50,kmaxx=7,imax=kmaxx+1,safe1=.25,safe2=.7,
     &           redmax=1.e-5,redmin=.7,tiny=1.e-30,scalmx=.1)
      integer i,iq,k,kk,km,kmax,kopt,nvold,nseq(imax)
c      real eps1,epsold,errmax,fact,h,red,scale,work,wrkmin,xest,xnew,    
c     &     a(imax),alf(kmaxx,kmaxx),dfdx(nmax),dfdy(nmax,nmax),
c     &     err(kmaxx),yerr(nmax),ysav(nmax),yseq(nmax)
      dimension a(imax),alf(kmaxx,kmaxx),dfdx(nmax),dfdy(nmax,nmax),
     &     err(kmaxx),yerr(nmax),ysav(nmax),yseq(nmax)
      logical first,reduct
      save a,alf,epsold,first,kmax,kopt,nseq,nvold,xnew
      data first/.true./,epsold/-1./,nvold/-1/
c----- Sequence is different from bsstep.      
      data nseq /2,6,10,14,22,34,50,70/
c----- Reinitialize also if nv has changed.   
      if (eps.ne.epsold.or.nv.ne.nvold) then
          hnext=-1.e29
          xnew=-1.e29
          eps1=safe1*eps
c----- Compute work coefficients A_k. (NR eq. 16.4.6)     
          a(1)=nseq(1)+1
          do 11 k=1,kmaxx
                a(k+1)=a(k)+nseq(k+1)
 11       continue
c----- Compute alpha(k,q). (NR eq. 16.4.11)
          do 13 iq=2,kmaxx
                do 12 k=1,iq-1
                      alf(k,iq)=eps1**((a(k+1)-a(iq+1))/
     &                          ((a(iq+1)-a(1)+1.)*(2*k+1)))
 12             continue
 13       continue
c-----
          epsold=eps
          nvold=nv
          a(1)=nv+a(1)
          do 14 k=1,kmaxx
                a(k+1)=a(k)+nseq(k+1)
 14       continue
c----- Determine optimal row number for convergence. 
          do 15 kopt=2,kmaxx-1
                if (a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt)) goto 1
 15       continue
 1        kmax=kopt
      endif
      h=htry
c----- Save the starting values. 
      do 16 i=1,nv
            ysav(i)=y(i)
 16   continue
c----- Evaluate Jacobian. 
      call jac(x,y,dydx,dfdx,dfdy)
      if (kpr.eq.1) return
      if (h.ne.hnext.or.x.ne.xnew) then
          first=.true.
          kopt=kmax
      endif
      reduct=.false.
 2    do 18 k=1,kmax
            xnew=x+h
            if (xnew.eq.x) then
                icon=16000
                return
            endif
            call simpr (ysav,dydx,dfdx,dfdy,nmax,nv,x,h,nseq(k),yseq,
     &                  derivs)
            xest=(h/nseq(k))**2
            call pzextr (k,xest,yseq,y,yerr,nv)
            if (k.ne.1) then
                errmax=tiny
                do 17 i=1,nv
                      errmax=max(errmax,abs(yerr(i)/yscal(i)))
 17             continue
                errmax=errmax/eps
                km=k-1
                err(km)=(errmax/safe1)**(1./(2*km+1))
            endif
            if (k.ne.1.and.(k.ge.kopt-1.or.first)) then
                if (errmax.lt.1.) goto 4
                if (k.eq.kmax.or.k.eq.kopt+1) then
                    red=safe2/err(km)
                    goto 3
                else if (k.eq.kopt) then
                    if (alf(kopt-1,kopt).lt.err(km)) then
                        red=1./err(km)
                        goto 3
                    endif
                else if (kopt.eq.kmax) then
                    if (alf(km,kmax-1).lt.err(km)) then
                        red=alf(km,kmax-1)*safe2/err(km)
                        goto 3
                    endif
                else if (alf(km,kopt).lt.err(km)) then
                    red=alf(km,kopt-1)/err(km)
                    goto 3
                endif
            endif
 18   continue
c----- Reduce stepsize by at least REDMIN and at most REDMAX.
 3    red=min(red,redmin)
      red=max(red,redmax)
      h=h*red
      reduct=.true.
c----- Try again! 
      goto 2
c----- Successful step taken. 
 4    x=xnew
      hdid=h
      first=.false.
c----- Compute optimalrow for convergence and corresponding stepsize.  
      wrkmin=1.e35
      do 19 kk=1,km
         fact=max(err(kk),scalmx)
         work=fact*a(kk+1)
         if (work.lt.wrkmin) then
             scale=fact
             wrkmin=work
             kopt=kk+1
         endif
 19   continue
      hnext=h/scale
c --- Check for possible order increase, but not if stepsize was just reduced. 
      if (kopt.ge.k.and.kopt.ne.kmax.and..not.reduct) then
          fact=max(scale/alf(kopt-1,kopt),scalmx)
          if (a(kopt+1)*fact.le.wrkmin) then
              hnext=h/fact
              kopt=kopt+1
          endif
      endif
      icon=10
      return

Here is the call graph for this function:

Here is the caller graph for this function:

double precision function temp	(		DE,
			PR
	)

Definition at line 1347 of file accdisk.f.

C*********************************************
C
      include 'para.h'
C
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c     A1    : 3/a*(RMU*rho0)
c     A2    : -3/a*P0 (=T0^4)
c     TMAX  : Temperature for Ptotal=Prad (P0=a/3*T0^4)  
c     TMIN  : 0.
c
c     FX    : TMAX^4 + A1*TMAX + A2 == 3/a*(Prad + Pgas - Ptotal) > 0
c     FN    : TMIN^4 + A1*TMIN + A2 ==                            < 0
c
c     TH    : Mean Temperature
c     FH    : TH^4 + A1*TH + A2 ==> Basically, this quantity = 0 !
c
c           if T > Treal -> F(T) > 0 
c           if T < Treal -> F(T) < 0
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      a1  = 3.*rmu*de/aa
      a2  =-3.*pr/aa
      eps = 1.0e-12
c      eps = 1.0e-10
      tmax= (-a2)**0.25
      tmin= 0.
      fx  = tmax**4+a1*tmax+a2
      fn  = tmin**4+a1*tmin+a2
      DO 40 itr=1,10
         th  = 0.5*(tmax+tmin)
         fh  = th**4+a1*th+a2
c
         IF (fh*fn .LT. 0.) tmax=th
         IF (fh*fn .LT. 0.) fx=fh
c     
         IF (fh*fx .LT. 0.) tmin=th
         IF (fh*fx .LT. 0.) fn=fh
 40   CONTINUE
C
      xxx = th
C
c +++++ For obtain the solutiono of TEMPERATURE using Newton Method ++++
c
c     xxx : Rough solution of TE
c     xx1 : Better solution of TE
c                 d(Prad+Pgas)
c     xx1 = xxx - ------------    
c                     dT
c     TEMP: Final solution of TE
c     
      DO 100 itr=1,20
         xx1 = xxx-(xxx**4+a1*xxx+a2)/(4.*xxx**3+a1)
         IF (abs((xx1-xxx)/xx1) .LT. eps) GO TO 120
         xxx = xx1
 100  CONTINUE
c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c +++++ For not converge ! +++++
      WRITE (6,110)
      print *,de,pr,tmax,tmin,fx,fn,th,fh,xxx,xx1
110   FORMAT (1h ,'LOOP COUNT (TEMP) IS OVER MAXIMUM')
      stop
c ++++++++++++++++++++++++++++++
c
120   CONTINUE
c
      temp=xx1
C
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function: