c+++++7++1+++++++++2+++++++++3+++++++++4+++++++++5+++++++++6+++++++++7++
C     < AD1DLW.FORT77 >
C     NUMERICAL HYDRODYNAMICS OF THIN ACCRETION DISKS
C
C               ONE DIMENSION (R)
C               PSEUDO NEWTONIAN POTENTIAL
C               EXPLICIT TIME INTEGRATION
C               BY LAX-WENDROFF METHOD
C
C               CODED BY R.MATUMOTO       83-10-20  V01L01
C                                         83-10-28  V01L02
C                                         83-11-09  V01L03
C                                         83-11-11  V01L04
C                                         88-05-13  VO2L01
c     "s-shaped curve" subroutine  01-10-18  
c     check 03-03-03
c     check 03-04-14
c     include 'para.h'
c 2003-6-11  このプログラムは alpha=0.3の計算専用とする！
C===========================================================
C**********************
C     MAIN ROUTINE
C**********************
C
      program AT
      include 'para.h'

      NS   = 0
      TIME = 0.
      TPTB = 1.
C----- Read Parameter ------
      CALL RDPRM             ! Read Parameter 
c---------------------------
c      KI=0 -> steady calculation 
c      KI=1 -> time-dependent calculation
c----- Begin Steady State Calculation -----
      if (ki.eq.0) then
          if (imore.eq.2) open (8,file=file8,status="unknown")
          lns=0
          ns=0
          nsx=0
          CALL INTMDL
      else

c----- Begin Time-Dependent Calculation -----
          open (7,file=file7,status="old")
          open (8,file=file8,status="unknown")
          open (9,file=file9,status="unknown")
          open (10,file=file10,status="unknown")
c
          kf=0
          call rdfl7
      endif
      CALL PRINT
      IF (KI.EQ.0) CALL FILE
C
c      print *,"end steady"
c      pause
c##### Time-dependent loop ##########
      do 100 ns=lns+1,nsx            ! ns: Time step 
         TIME = TIME+DT
c         print *,"time-depend",ns
c         pause
C     if (time.le.TPTB)
         if (ns.eq.lns+1) then
c ----- Initial perturbation! -----
c         if (ns.eq.0) then
c            print *,"prtrb"
c            pause
c            call PRTRB1         ! If you add a perturbation, 
c     call PRTRB2      ! you can select the form of perturb. 
c             call PRTRB3
             call prtrb4
         end if
c ---------------------------------
c         print *,"LW1D"
         CALL LW1D                 ! Lax-Wendroff subrooutine
c         print *,"PRINT"
         CALL PRINT             ! Print 
c         print *,"CFL"
         CALL CFL               ! CFL condition
c         print *,"FILE"
         CALL FILE              ! Open file
c ---------------------------------
 100  CONTINUE
c      STOP
c##### End of Time-dependent loop #####
      END

C*************************
      SUBROUTINE RDPRM
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
      END

C
C**************************
      SUBROUTINE INTMDL
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
      END
C
C
C******************************************
c      SUBROUTINE SADMK (R,ETA,DL,DE,VX,VP,TE)
      SUBROUTINE SADMK (R,ETA,DL)
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
      END
C
C******************************************
      SUBROUTINE SADMK2 (R,ETA,DL,DE,VX,VP,TE)
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
      END
C
C******************************************
      SUBROUTINE SADMN (R,ETA,DL)
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
      END
C******************************************
      SUBROUTINE SADMN2 (R,ETA,DL,DE,VX,VP,TE)
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
      END
C****************************************
      SUBROUTINE INTGRT (R,ETA,DL,IC)
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
      END
C     
C     
C*******************************************************
      SUBROUTINE FUN (XI,YI,YP)
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
      END
C
C
C**************************************
      subroutine jac(xi,yi,yp,dfdr,pdj)
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
      END
C
C
C*********************************************
      DOUBLE PRECISION FUNCTION TEMP(DE,PR)
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
      END
C
C
C******************************
      DOUBLE PRECISION FUNCTION ondo(DE,PR,taup,tau)
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
      END
C
C***************************
      DOUBLE PRECISION FUNCTION FNIN(NPL)
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
      END
C
C***************************
      DOUBLE PRECISION FUNCTION FNJN(NPL)
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
      END
C
C
C**************************
      SUBROUTINE CNVRBS
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
      END
C
C
C**************************
      SUBROUTINE LW1D
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の定義（中心差分）
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,i10,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)
      END
C
C
C**************************
      SUBROUTINE CFL
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
      END
C
C**************************
      SUBROUTINE PRTRB1
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
      end
C**************************
      SUBROUTINE PRTRB2
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
      END
C**************************
      SUBROUTINE PRTRB3
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=30.0
c      xlamda=pxlin
      xlamda=10.0
      rp=1.0
      do 10 j=1,jx-1
         do 20 k=1,1
            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
      write(6,*) "PERTRB 3"
      return
      end
C**************************
      SUBROUTINE PRTRB4
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

         do 20 k=1,1
            dss = 0.01*exp(-(x1(j)-rc)**2/0.1**2)*sin(xkk*x1(j))*w0(j,k)
            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
      end
C**************************
      SUBROUTINE PRINT
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に追加
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     
      END
C     
C
C**************************
      SUBROUTINE FILE
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)
      END
************************************************************************
      subroutine rdfl7
************************************************************************
      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
      end
***********************************************************************
***** stifbs  from Numerical Recipes (p737)
***********************************************************************
      subroutine stifbs (y,dydx,nv,x,htry,eps,yscal,hdid,hnext,
     &                   derivs,icon)
      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
      end
*******************************************************************
***** simpr (Semi-implicit Extrapolation Method by NR )
*******************************************************************
      subroutine simpr (y,dydx,dfdx,dfdy,nmax,n,xs,htot,nstep,yout,
     &                  derivs)
      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
      end
***************************************************************
***** pzextr (The polynomial extrapolation routine by NR)
***************************************************************
      subroutine pzextr (iest,xest,yest,yz,dy,nv)
      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
      end         
********************************************************************
***** ludcmp
********************************************************************
      subroutine ludcmp (a,n,np,indx,d)
      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
      end
***********************************************************************
***** lubksb
***********************************************************************
      subroutine lubksb (a,n,np,indx,b)
      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
      end     
