black hole accretion calculation
accdisk.f
Go to the documentation of this file.
1 c+++++7++1+++++++++2+++++++++3+++++++++4+++++++++5+++++++++6+++++++++7++
2 C < AD1DLW.FORT77 >
3 C NUMERICAL HYDRODYNAMICS OF THIN ACCRETION DISKS
4 C
5 C ONE DIMENSION (R)
6 C PSEUDO NEWTONIAN POTENTIAL
7 C EXPLICIT TIME INTEGRATION
8 C BY LAX-WENDROFF METHOD
9 C
10 C CODED BY R.MATUMOTO 83-10-20 V01L01
11 C 83-10-28 V01L02
12 C 83-11-09 V01L03
13 C 83-11-11 V01L04
14 C 88-05-13 VO2L01
15 c "s-shaped curve" subroutine 01-10-18
16 c check 03-03-03
17 c check 03-04-14
18 c include 'para.h'
19 c 2003-6-11 このプログラムは alpha=0.3の計算専用とする！
20 C===========================================================
21 C**********************
22 C MAIN ROUTINE
23 C**********************
24 C
25  program at
26  include 'para.h'
27 
28  ns = 0
29  time = 0.
30  tptb = 1.
31 C----- Read Parameter ------
32  CALL rdprm ! Read Parameter
33 c---------------------------
34 c KI=0 -> steady calculation
35 c KI=1 -> time-dependent calculation
36 c----- Begin Steady State Calculation -----
37  if (ki.eq.0) then
38  if (imore.eq.2) open (8,file=file8,status="unknown")
39  lns=0
40  ns=0
41  nsx=0
42  CALL intmdl
43  else
44 
45 c----- Begin Time-Dependent Calculation -----
46  open (7,file=file7,status="old")
47  open (8,file=file8,status="unknown")
48  open (9,file=file9,status="unknown")
49  open (10,file=file10,status="unknown")
50 c
51  kf=0
52  call rdfl7
53  endif
54  CALL print
55  IF (ki.EQ.0) CALL file
56 C
57  print *,"end steady"
58  pause
59 c##### Time-dependent loop ##########
60  do 100 ns=lns+1,nsx ! ns: Time step
61  time = time+dt
62  print *,"time-depend",ns
63  pause
64 C if (time.le.TPTB)
65 c if (ns.eq.lns+1) then
66 c ----- Initial perturbation! -----
67  if (ns.eq.0) then
68  print *,"prtrb"
69  pause
70  call prtrb1 ! If you add a perturbation,
71 c call PRTRB2 ! you can select the form of perturb.
72 c call PRTRB3
73  end if
74 c ---------------------------------
75  print *,"LW1D"
76  CALL lw1d ! Lax-Wendroff subrooutine
77  print *,"PRINT"
78  CALL print ! Print
79  print *,"CFL"
80  CALL cfl ! CFL condition
81  print *,"FILE"
82  CALL file ! Open file
83 c ---------------------------------
84  100 CONTINUE
85 c STOP
86 c##### End of Time-dependent loop #####
87  END
88 
89 C*************************
90  SUBROUTINE rdprm
91 C*************************
92 
93  include 'para.h'
94 
95 C***** CONSTANTS *****
96  aa = 7.56464e-15 ! Radiation constant (erg/cm^3/K^4)
97  cc = 2.997925e10 ! Speed of light
98  gg = 6.672e-8 ! Gravitational constant
99  xmu = 0.5 ! 1/2
100  xmh = 1.660531e-24 ! Atomic mass (g)
101  xkb = 1.38062e-16 ! Boltzman constant (erg/K)
102  rmu = xkb/(xmu*xmh) ! 2 xkb/xmh (coefficient of gas pres)
103  xmsun = 1.989e33 ! Solar mass (g)
104  xkes = 0.34 ! Electron Scattering Opacity:
105 c XKKR = 0.34*6E24 ! Free-free opacity coefficient
106 c xkes = 0.4
107  xkkr = 0.64e23
108  xeff = 6.22e20
109 c
110 c----- Ho-shi's parameter(1977) ----------------------------------------
111  npl = 3 ! Polytropic index (gamma=1+1/N)
112  xnpl1 = npl+1 ! XNPL1:N+1
113  ain = fnin(npl) ! I_N -> related to heght integration.
114  ainp1 = fnin(npl+1) ! AIN : I_N
115  ainp1 = fnin(npl+1) ! AINP1 : I_N+1
116  ajn = fnjn(2*npl) ! AJN : I_2N
117 c----- Specific heat ? -----
118  gm = 5./3.
119  gm1 = gm-1.
120  pi = 3.14159263
121 C
122 C ***** READ PARAMETERS *****
123 C
124  open (9,file="at.job",status="old") ! Data input
125  read (9,*)
126  read (9,*)
127  read (9,901) date
128  read (9,*)
129  read (9,*) bhm,alpha,xmd,xlinm,xlinp
130  read (9,*)
131  read (9,*) mesh,rin,rin2,rout1,rout2
132  read (9,*)
133  read (9,*) dt,dtmax,dtmin,nsx,lns
134  read (9,*)
135  read (9,*) (iprint(i),i=1,4),ksfl,irstrt
136  read (9,*)
137  read (9,*) hini,relerr,abserr,qb,qu
138  read (9,*)
139  read (9,902) ki,file7,file8,imore,xkappa
140  read (9,*)
141  read (9,903) klc,file9,file10
142  read (9,*)
143  read (9,*) kr,pxmd,pxlin,kcool
144  read (9,*)
145  read (9,*) ktcond
146  read (9,*)
147  read (9,*) cfc
148 
149  901 format (4x,a10)
150  902 format (4x,i2,7x,a6,4x,a6,4x,i2,8x,e12.5)
151  903 format (4x,i6,3x,a6,4x,a8)
152 C
153 C ***** WRITE INPUT PARAMETERS *****
154 C
155  WRITE (6,*) '***** INPUT PARAMETERS *****', date
156  WRITE (6,100) bhm,alpha,xmd,xlinm,xlinp
157  WRITE (6,200) mesh,rin,rout1,rout2
158  WRITE (6,300) dt,dtmax,dtmin,nsx
159 c WRITE (6,400) (IPRINT(I),I=1,4),KSFL
160  WRITE (6,400) qb,qu,xkappa
161  write (6,500) hini,relerr,abserr
162 C
163  100 FORMAT (/1h ,'BHM =',1pe15.7,5x,'ALPHA=',1pe15.7,5x,'XMDOT=',
164  * 1pe15.7/1h ,'XLINM=',1pe15.7,5x,'XLINP=',1pe15.7)
165  200 FORMAT (1h ,'MESH =',i15,5x,'RIN ='
166  * ,1pe15.7/1h ,'ROUT1=',1pe15.7,5x,'ROUT2=',1pe15.7)
167  300 FORMAT (1h ,'DT =',1pe15.7,5x,'DTMAX=',1pe15.7,5x,'DTMIN=',
168  * 1pe15.7/1h ,'NSX =',i15)
169 c 400 FORMAT (1H ,'IPRINT=',4(I6),5X,'KSFL =',I8//)
170  400 FORMAT (1h ,'QB=',1pe15.7,5x,"QU=",1pe15.7,5x,'XKAPPA=',
171  * 1pe15.7,//)
172  500 format (1h ,"HINI=",1pe9.1,5x,"RELERR=",1pe9.1,5x
173  & ,"ABSERR=",1pe9.1//)
174 C
175 c--- B1,B2,B3,B4,B5 -> vertical integration factor (for normalize). ---
176 c
177  rg = 2.*gg*bhm*xmsun/(cc*cc) ! Schwarzschild radius (cm)
178 c XMCRIT = 32.*PI*CC*RG/XKES
179  xmcrit = 2.*pi*cc*rg/xkes ! Critical accretion rate (L_E/eta/c^2)
180  xmdi = xmd ! normalized accretion rate
181  xmdot = xmd*xmcrit
182  pxmdot = pxmd*xmcrit
183  al7 = alpha**7
184  ai5 = (ain**5)/(ainp1**4)
185  ff1 = (128.*pi*pi*aa*cc/(3.*xkes*xmdot*xmdot*al7))*ai5
186  * *((cc/rmu)**4)*((cc*rg)**3)
187  b1 = ain*cc*rg*rg
188  b2 = (ainp1/ain)/(cc*cc)
189  b3 = ainp1*rg*rg/cc
190  b4 = rg/(cc*cc)
191  b5 = 2.*xnpl1/(cc*cc)
192 
193  RETURN
194  END
195 
196 C
197 C**************************
198  SUBROUTINE intmdl
199 C**************************
200 C
201 C -------------------------------------------------------------
202 C INITIAL CONDITION
203 C STEADY MODEL OF THIN ACCRETION DISK AROUND A BLACK HOLE
204 C TRANSONIC SOLUTION IS FOUND BY ADJUSTING LIN.
205 C -------------------------------------------------------------
206 C
207  include 'para.h'
208 C
209  il = 0
210 C
211 C -------------------------
212 c DO 100 UNTIL (IC .EQ. 0)
213  200 continue
214 C -------------------------
215 c
216 c ----------------------------------------------------------------------
217 c XLIN : L_in, angular momentum at marginally stable orbit.
218 c XLINP : Lower limit of the inner angular momentum
219 c XLINM : Upper limit of the inner angular momentum
220 c ---------------------------------------------------------------------
221 
222  xlin = 0.5*(xlinp+xlinm)
223  r = rout1
224 
225  if (kr .eq. 1) then
226  CALL sadmk(r,eta,dl)
227  else if (kr .eq. 0) then
228  call sadmn(r,eta,dl)
229  end if
230 cc do 2 i=1,2
231 cc if (i.eq.1) r=rout1
232 cc if (i.eq.2) r=rout2
233 c call sadmn(r,eta,dl)
234 cc2 continue
235 c --------------------------------------------
236  if (imore.eq.2) print *," IL=",il+1 ! IL -> Iteration count
237  CALL intgrt(r,eta,dl,ic)
238  IF (ic .EQ. 1 ) xlinm = xlin
239  IF (ic .EQ. -1) xlinp = xlin
240  il = il+1
241 
242  IF (il .GT. 50) THEN ! Set Max Iteration
243  if (imore.eq.2) WRITE (6,81) il
244  if (imore.eq.2) write (6,82)
245  if (imore.eq.2) write (6,*) "xlin=", xlin, "IC=",ic
246  stop
247  ENDIF
248 C
249  if (imore.eq.2) WRITE (6,80) xlin,r,ic
250  80 FORMAT (1h ,'LIN = ',1pe22.14,5x,'R = ',1pe15.7,5x,'IC = ',i5)
251  81 format ("LOOP COUNT (INTMDL) IS OVER MAXIMUM IL=",i3)
252  82 format (" IL IC LIN")
253 C
254  if (ic.ne.0) goto 200
255 100 CONTINUE
256 C
257  CALL cnvrbs
258 C
259  RETURN
260  END
261 C
262 C
263 C******************************************
264 c SUBROUTINE SADMK (R,ETA,DL,DE,VX,VP,TE)
265  SUBROUTINE sadmk (R,ETA,DL)
266 C******************************************
267 C
268 C -------------------------------------------------------------------
269 C STANDARD ACCRETION DISK MODEL
270 C OPTICALLY THICK, GEOMETRICALLY THIN, KEPLERIAN, ALPHA-MODEL
271 C PSEUDO NEWTONIAN POTENTIAL
272 C
273 C ETA = (VR**2)/(WW/SS) DL = XL-XLK
274 C --------------------------------------------------------------------
275 C
276  include 'para.h'
277 C
278  dimension ph(4)
279 C
280 c print *,"SADMK"
281  xkes1 = xkes
282  xkkr1 = xkkr
283 
284  rm1 = r-1
285  r2 = r*r
286  r3 = r2*r
287  eps = 1.0e-6
288 
289  ok = (r/rm1)*sqrt(.5/r3) ! OK : Kepler angular velocity
290  vk = r*ok ! VK : Kepler velocity
291  xlk = r*vk ! XLK : Kepler angular momentum
292  dlkdr = xlk*(1.5/r - 1./rm1) !DLKDR : dlnLk/dr
293  dlnok = -1./rm1 - .5/r !DLNOK : d(lnOK)/dlnR
294  ww = xmdot*(xlk-xlin)/(2.*pi*b3*r*r*alpha)
295  ff2 = (-3.*xkes1*b3*r*alpha*ww*ok*dlnok/(32.*b4*aa*cc))*ww
296  ff3 = ok*ok*ww*ww/(4.*b5)
297  ff4 = (aa/3.)*(ff2**0.9)/(sqrt(ff3)*(rmu**0.4))
298  betan = 0.01
299  betax = 1.00
300  DO 10 i=1,10
301  beta = 0.5*(betan+betax)
302  fff = 1.-beta-ff4*(beta**0.4)
303  IF (fff .GT. 0.) betan=beta
304  IF (fff .LE. 0.) betax=beta
305 10 CONTINUE
306 C
307 c print *,"S1"
308  pr = (3./aa)*(ff3/ff2)*(1.-beta)! Pressure (at equatrial plane)
309  de = (pr**3)/ff3 !DE : density (g/cm^3)
310  te = (de*ff2/pr)**0.25 !TE : temperature (K)
311  hh = (sqrt(b5*pr/de))/ok !HH : Scale height
312  vr = xmdot/(4.*pi*b1*r*de*hh) !VR : radial velocity
313 c write(6,*) pr,de,te
314 c write(6,*) hh,vr
315 C
316 C ***** CORRECTION FOR OPACITY *****
317 C
318 C --- ITERATION ---
319 c ----------------------------------------------------------------------
320 c PR : a/3 * TE^4 (radiation pressure)+2kb/mh*de*TE (gas pressure)
321 c BETA : P_gas (Gas pressure)/P_total (Total pressure)
322 c DEMN : I_N DE (mean density)
323 c TEMN : 2/3*TE (mean temperature)
324 c XKFF : Free-Free opacity
325 c XK : Mean Opacity
326 c XKXF : free-free opacity/mean opacity
327 c HH : Scale Height
328 c FM : Radiative flux
329 c FF3 :
330 c FF4 :
331 c ----------------------------------------------------------------------
332 c print *,"S2"
333  DO 200 itr=1,20
334  pr = (aa/3.)*(te**4)+rmu*de*te
335 c if (pr.le.0.) print *,"PR=",pr,itr
336  beta = rmu*de*te/pr
337  demn = ain*de
338  temn = 2.*te/3.
339  xkff = xkkr1*demn*(temn**(-3.5))
340  xk = xkes1+xkff
341 c if (xk.le.0.) print *,"XK=",xk
342  xkxf = xkff/xk
343  hh = (sqrt(b5*pr/de))/ok
344 c
345 c write(6,*) itr,pr,beta
346 c
347  if (de.le.0.) print *,"DE1=",de
348  if (hh.le.0.) print *,"HH=",hh
349 
350  fm = 2.*pi*b4*r*16.*aa*cc*(te**4)/(3.*xk*de*hh)
351  ff3 = 4.*pi*b3*r*r*alpha*pr*hh
352  ff4 = -xmdot*(xlk-xlin)*ok*dlnok
353 C
354  ph(3)= xmdot*(xlk-xlin)-ff3
355  ph(4)= ff4-fm
356 C
357  a1 = 1.-beta
358  a4 = 4.-3.*beta
359 C
360 c if (te.le.0.) print *,"TE=",te
361  dlnhdd = -0.5*a1/de
362  dlnhdt = 0.5*a4/te
363  dlnpdd = beta/de
364  dlnpdt = a4/te
365  dlnkdd = xkxf/de
366  dlnkdt = -3.5*xkxf/te
367  dlnfdd = -dlnkdd-dlnhdd-1./de
368  dlnfdt = 4./te-dlnkdt-dlnhdt
369 C
370 c ---------------------------------------------
371 c | D3DD D3DT |
372 c det. | | == D3DD*D4DT-D3DT*D4DD
373 c | D4DD D4DT |
374 c ---------------------------------------------
375 C ** JACOBIANS **
376 C
377  d3dd = -ff3*(dlnpdd+dlnhdd)
378  d3dt = -ff3*(dlnpdt+dlnhdt)
379  d4dd = -fm*dlnfdd
380  d4dt = -fm*dlnfdt
381  det = d3dd*d4dt-d3dt*d4dd
382 C
383 c --------------------------------------------
384 c |D3DD D3DT| |DELD| = |PH(3)|
385 c |D4DD D4DT| |DELT| |PH(4)|
386 c
387 c |DELD| 1 |D4DT -D3DT| |PH(3)|
388 c | | = --- | |*| |
389 c |DELT| det.|-D4DD D3DD| |PH(4)|
390 c --------------------------------------------
391 C ** CORRECTIONS **
392 C
393 c if (det.le.0.) print *,"DET=",det
394  deld = (ph(3)*d4dt-ph(4)*d3dt)/det
395  delt = (d3dd*ph(4)-d4dd*ph(3))/det
396 C
397 C ** APPLY CORRECTIONS **
398 C
399  de = de-deld
400  te = te-delt
401 c WRITE (6,40) ITR,DE,TE
402  40 FORMAT (1h ,i7,5x,'DE =',1pe12.4,5x,'TE =',1pe12.4)
403 C
404 C ** CHECK CONVERGENCE **
405 C
406 c if (de.le.0.) print *,"DE2=",de
407 c if (te.le.0.) print *,"TE=",te
408  IF (abs(deld/de) .GE. eps) GO TO 100
409  IF (abs(delt/te) .GE. eps) GO TO 100
410 C
411  GO TO 300
412 C
413  100 CONTINUE
414  200 CONTINUE
415 C
416  if (imore.eq.2) WRITE (6,210)
417  210 FORMAT (1h ,'LOOP COUNT (SADM) IS OVER MAXIMUM !')
418  stop
419 C
420  300 CONTINUE
421 C
422 c----- Physical Quantities at ROUT (in SSD). -----
423 c
424  pr = (aa/3.)*(te**4)+rmu*de*te
425  hh = (sqrt(b5*pr/de))/ok
426 c if (hh.le.0.) print *,"HH=",hh
427  vr = xmdot/(4.*pi*b1*r*hh*de)
428  eta = vr*vr/(b2*pr/de)
429  ws = ((xlk-xlin)**2)*eta/(r2*alpha*alpha)
430 c if (xlk-xlin.eq.0.) print *,"XLK-XLIN",xlk-xlin
431  dlnwr = -2./r+dlkdr/(xlk-xlin)
432  dl = r3*ws*(dlnwr+dlnok)/(2.*xlk)
433 C
434  vx =-vr ! radial velocity
435  vp = dl/r ! rotational velocity
436  demn = ain*de ! Mean density
437  temn = 2*te/3.0 ! Mean temperature
438 c TAUES = 2.*XKES*DE*HH*RG
439  taues = xkes*de*hh*rg ! Optical depth (electron scattering)
440  xkff = xkkr*demn*(temn**(-3.5)) ! Mean Free-Free opacity
441 c TAUFF = 2.*XKFF*DE*HH*RG
442  tauff = xkff*de*hh*rg ! Optical depth (free-free absorption)
443 c if (tauff.le.0) print *,"TAUABS=",tauff,"LW1D4"
444  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg ! Brems cooling
445  taup = qbr/(8.*aa*cc*(te**4)) ! Optical depth (planck)
446  tau = taues+tauff ! Total optical depth
447  tauabso(jx)=taup
448  tauo(jx)=tau
449  ss = 2*ain*de*hh*rg ! Surface density
450 c WRITE(6,1010) TIME,SS,TE,XMD
451 C
452  IF (iprint(1) .EQ. 0) RETURN
453 C
454 c
455  if (imore.eq.2) WRITE (6,400) de,vr,dl,te,eta
456 400 FORMAT (/1h ,'***** STANDARD MODEL *****'//1h ,'DE =',1pe15.7,
457  * 5x,'VR =',1pe15.7,5x,'DL =',1pe15.7/1h ,'TE =',
458  * 1pe15.7,5x,'ETA =',1pe15.7//)
459 1010 FORMAT (20(1pe12.4))
460 C
461 c print *,"S-END"
462  RETURN
463  END
464 C
465 C******************************************
466  SUBROUTINE sadmk2 (R,ETA,DL,DE,VX,VP,TE)
467 C******************************************
468 C
469  include 'para.h'
470 C
471  dimension ph(4)
472 C
473 c print *,"SADMK"
474  xkes1 = xkes
475  xkkr1 = xkkr
476 
477  rm1 = r-1
478  r2 = r*r
479  r3 = r2*r
480  eps = 1.0e-6
481 
482  ok = (r/rm1)*sqrt(.5/r3) ! OK : Kepler angular velocity
483  vk = r*ok ! VK : Kepler velocity
484  xlk = r*vk ! XLK : Kepler angular momentum
485  dlkdr = xlk*(1.5/r - 1./rm1) !DLKDR : dlnLk/dr
486  dlnok = -1./rm1 - .5/r !DLNOK : d(lnOK)/dlnR
487  ww = xmdot*(xlk-xlin)/(2.*pi*b3*r*r*alpha)
488  ff2 = (-3.*xkes1*b3*r*alpha*ww*ok*dlnok/(32.*b4*aa*cc))*ww
489  ff3 = ok*ok*ww*ww/(4.*b5)
490  ff4 = (aa/3.)*(ff2**0.9)/(sqrt(ff3)*(rmu**0.4))
491  betan = 0.01
492  betax = 1.00
493  DO 10 i=1,10
494  beta = 0.5*(betan+betax)
495  fff = 1.-beta-ff4*(beta**0.4)
496  IF (fff .GT. 0.) betan=beta
497  IF (fff .LE. 0.) betax=beta
498 10 CONTINUE
499 C
500 c print *,"S1"
501  pr = (3./aa)*(ff3/ff2)*(1.-beta)! Pressure (at equatrial plane)
502  de = (pr**3)/ff3 !DE : density (g/cm^3)
503  te = (de*ff2/pr)**0.25 !TE : temperature (K)
504  hh = (sqrt(b5*pr/de))/ok !HH : Scale height
505  vr = xmdot/(4.*pi*b1*r*de*hh) !VR : radial velocity
506 c print *,"S2"
507  DO 200 itr=1,20
508  pr = (aa/3.)*(te**4)+rmu*de*te
509 c if (pr.le.0.) print *,"PR=",pr,itr
510  beta = rmu*de*te/pr
511  demn = ain*de
512  temn = 2.*te/3.
513  xkff = xkkr1*demn*(temn**(-3.5))
514  xk = xkes1+xkff
515 c if (xk.le.0.) print *,"XK=",xk
516  xkxf = xkff/xk
517  hh = (sqrt(b5*pr/de))/ok
518 c
519 c write(6,*) itr,pr,beta
520 c
521  if (de.le.0.) print *,"DE1=",de
522  if (hh.le.0.) print *,"HH=",hh
523 
524  fm = 2.*pi*b4*r*16.*aa*cc*(te**4)/(3.*xk*de*hh)
525  ff3 = 4.*pi*b3*r*r*alpha*pr*hh
526  ff4 = -xmdot*(xlk-xlin)*ok*dlnok
527 C
528  ph(3)= xmdot*(xlk-xlin)-ff3
529  ph(4)= ff4-fm
530 C
531  a1 = 1.-beta
532  a4 = 4.-3.*beta
533 C
534 c if (te.le.0.) print *,"TE=",te
535  dlnhdd = -0.5*a1/de
536  dlnhdt = 0.5*a4/te
537  dlnpdd = beta/de
538  dlnpdt = a4/te
539  dlnkdd = xkxf/de
540  dlnkdt = -3.5*xkxf/te
541  dlnfdd = -dlnkdd-dlnhdd-1./de
542  dlnfdt = 4./te-dlnkdt-dlnhdt
543 C
544 c ---------------------------------------------
545 c | D3DD D3DT |
546 c det. | | == D3DD*D4DT-D3DT*D4DD
547 c | D4DD D4DT |
548 c ---------------------------------------------
549 C ** JACOBIANS **
550 C
551  d3dd = -ff3*(dlnpdd+dlnhdd)
552  d3dt = -ff3*(dlnpdt+dlnhdt)
553  d4dd = -fm*dlnfdd
554  d4dt = -fm*dlnfdt
555  det = d3dd*d4dt-d3dt*d4dd
556 C
557 C ** CORRECTIONS **
558 C
559 c if (det.le.0.) print *,"DET=",det
560  deld = (ph(3)*d4dt-ph(4)*d3dt)/det
561  delt = (d3dd*ph(4)-d4dd*ph(3))/det
562 C
563 C ** APPLY CORRECTIONS **
564 C
565  de = de-deld
566  te = te-delt
567 c WRITE (6,40) ITR,DE,TE
568  40 FORMAT (1h ,i7,5x,'DE =',1pe12.4,5x,'TE =',1pe12.4)
569 C
570 C ** CHECK CONVERGENCE **
571 C
572 c if (de.le.0.) print *,"DE2=",de
573 c if (te.le.0.) print *,"TE=",te
574  IF (abs(deld/de) .GE. eps) GO TO 100
575  IF (abs(delt/te) .GE. eps) GO TO 100
576 C
577  GO TO 300
578 C
579  100 CONTINUE
580  200 CONTINUE
581 C
582  if (imore.eq.2) WRITE (6,210)
583  210 FORMAT (1h ,'LOOP COUNT (SADM) IS OVER MAXIMUM !')
584  stop
585 C
586  300 CONTINUE
587 C
588 c----- Physical Quantities at ROUT (in SSD). -----
589 c
590  pr = (aa/3.)*(te**4)+rmu*de*te
591  hh = (sqrt(b5*pr/de))/ok
592 c if (hh.le.0.) print *,"HH=",hh
593  vr = xmdot/(4.*pi*b1*r*hh*de)
594  eta = vr*vr/(b2*pr/de)
595  ws = ((xlk-xlin)**2)*eta/(r2*alpha*alpha)
596 c if (xlk-xlin.eq.0.) print *,"XLK-XLIN",xlk-xlin
597  dlnwr = -2./r+dlkdr/(xlk-xlin)
598  dl = r3*ws*(dlnwr+dlnok)/(2.*xlk)
599 C
600  vx =-vr ! radial velocity
601  vp = dl/r ! rotational velocity
602  demn = ain*de ! Mean density
603  temn = 2*te/3.0 ! Mean temperature
604 c TAUES = 2.*XKES*DE*HH*RG
605  taues = xkes*de*hh*rg ! Optical depth (electron scattering)
606  xkff = xkkr*demn*(temn**(-3.5)) ! Mean Free-Free opacity
607 c TAUFF = 2.*XKFF*DE*HH*RG
608  tauff = xkff*de*hh*rg ! Optical depth (free-free absorption)
609 c if (tauff.le.0) print *,"TAUABS=",tauff,"LW1D4"
610  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg ! Brems cooling
611  taup = qbr/(8.*aa*cc*(te**4)) ! Optical depth (planck)
612  tau = taues+tauff ! Total optical depth
613  tauabso(jx)=taup
614  tauo(jx)=tau
615  ss = 2*ain*de*hh*rg ! Surface density
616 c WRITE(6,1010) TIME,SS,TE,XMD
617 C
618  IF (iprint(1) .EQ. 0) RETURN
619 C
620 c print *,"S-END"
621  RETURN
622  END
623 C
624 C******************************************
625  SUBROUTINE sadmn (R,ETA,DL)
626 C******************************************
627 c 2001/02/07 comments
628 C This subroutine is not used ! (except subroutine INTMDL)
629 C -------------------------------------------------------------------
630 C STANDARD ACCRETION DISK MODEL
631 C OPTICALLY THIN, GEOMETRICALLY THIN, KEPLERIAN, ALPHA-MODEL
632 C PSEUDO NEWTONIAN POTENTIAL
633 C
634 C ETA = (VR**2)/(WW/SS) DL = XL-XLK
635 C
636 C --------------------------------------------------------------------
637 C
638  include 'para.h'
639 C
640  rm1 = r-1.
641  r2 = r*r
642  r3 = r2*r
643  ok =(r/rm1)*sqrt(.5/r3)
644  vk = r*ok
645  xlk = r*vk
646  xll = xlk-xlin
647  dlnok =-1./rm1-.5/r
648  dokdr = ok*dlnok
649  cc2 = cc*cc
650  rg2 = rg*rg
651  ww =(xmdot*xll)/(2.*pi*r2*alpha)*cc/rg
652  wwtrue = ww
653  t1=1.e4
654  t2=1.e30
655 c100 continue
656 c ---------------------------------------
657 c te=1.e6
658  do 1 iss=1,10
659  te = sqrt(t1*t2)
660  wweps = 1.e-6
661 cc print *,"TE=",te
662  pd = rmu*te
663  hh = sqrt(2.*xnpl1*pd)/ok/cc
664  pr = ww/(2.*ainp1*hh)/rg
665  de = pr/pd
666  vr = xmdot/(4.*pi*r*ain*de*hh)/(cc*rg2)
667  demn = ain*de
668  temn = 2*te/3.0
669  xkff = xkkr*demn*(temn**(-3.5))
670  tauff = xkff*de*hh*rg
671  taues = xkes*de*hh*rg
672  tau = taues+tauff
673  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
674  taup = qbr/(8.*aa*cc*(te**4))
675  wwfalse= -8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)
676  * /(alpha*r*dokdr)*(rg/cc)
677  te4=-ww*(1.5*tau+sqrt(3.)+1./taup)*(alpha*r*dokdr)/(8.*aa)
678  * /(rg/cc)
679 cc wwfalse=-8.*AA*(TE**4)*TAUP/(ALPHA*R*DOKDR)*(RG/CC)
680 cc te4=-ww/TAUP*(ALPHA*R*DOKDR)/(8.*AA)/(RG/CC)
681  te=te4**0.25
682  if (wwfalse.gt.wwtrue) t1=te
683  if (wwfalse.lt.wwtrue) t2=te
684  if (abs(wwfalse-wwtrue)/wwtrue.lt.wweps) goto 100
685 c if (abs(wwfalse-wwtrue).ge.1.e4) goto 100
686 1 continue
687 C --------------------------------------------
688 100 continue
689  eta = ain*de*vr*vr/(ainp1*pr)*cc*cc
690  dl=0.
691  ss = 2*ain*de*hh*rg
692 C -----------------------------
693 c if (imore.eq.2) print *,"PR=",pr," HH=",hh," WW=",ww
694 c if (imore.eq.2) WRITE (6,400) DE,VR,DL,TE,SS
695  demn = ain*de
696  temn = 2*te/3.0
697  xkff = xkkr*demn*(temn**(-3.5))
698  tauff = xkff*de*hh*rg
699  taues = xkes*de*hh*rg
700  tau = taues+tauff
701  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
702  taup = qbr/(8.*aa*cc*(te**4))
703  taueff = sqrt(taup*(tau+2./sqrt(3.)))
704  beta=1.
705 c t1=tau+2./sqrt(3.)
706 c t2=1./taup
707 c if (t1.ge.t2) then
708 c taueff=1.
709 c else
710 c taueff=2.
711 c endif
712 c WRITE(6,1000) TIME,SS,TE,XMD,taueff,tau,taup
713 cc WRITE(6,1010) JX,R,BETA,TAUEFF,TE,SS,XMD
714 400 FORMAT (/1h ,'***** STANDARD MODEL *****'//1h ,'DE =',1pe15.7,
715  * 5x,'VR =',1pe15.7,5x,'DL =',1pe15.7/1h ,'TE =',
716  * 1pe15.7,5x,'SS =',1pe15.7//)
717 1000 FORMAT (20(1pe12.4))
718 1010 FORMAT (i5,20(1pe12.4))
719 C
720 C
721  RETURN
722  END
723 C******************************************
724  SUBROUTINE sadmn2 (R,ETA,DL,DE,VX,VP,TE)
725 C******************************************
726 c for optically thin accretion disk
727 C
728  include 'para.h'
729 C
730  rm1 = r-1.
731  r2 = r*r
732  r3 = r2*r
733  ok =(r/rm1)*sqrt(.5/r3)
734  vk = r*ok
735  xlk = r*vk
736  xll = xlk-xlin
737  dlnok =-1./rm1-.5/r
738  dokdr = ok*dlnok
739  cc2 = cc*cc
740  rg2 = rg*rg
741  ww =(xmdot*xll)/(2.*pi*r2*alpha)*cc/rg
742  wwtrue = ww
743  t1=1.e4
744  t2=1.e30
745 c100 continue
746 c ---------------------------------------
747 c te=1.e6
748  do 1 iss=1,10
749  te = sqrt(t1*t2)
750  wweps = 1.e-6
751 cc print *,"TE=",te
752  pd = rmu*te
753  hh = sqrt(2.*xnpl1*pd)/ok/cc
754  pr = ww/(2.*ainp1*hh)/rg
755  de = pr/pd
756  vr = xmdot/(4.*pi*r*ain*de*hh)/(cc*rg2)
757  demn = ain*de
758  temn = 2*te/3.0
759  xkff = xkkr*demn*(temn**(-3.5))
760  tauff = xkff*de*hh*rg
761  taues = xkes*de*hh*rg
762  tau = taues+tauff
763  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
764  taup = qbr/(8.*aa*cc*(te**4))
765  wwfalse= -8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)
766  * /(alpha*r*dokdr)*(rg/cc)
767  te4=-ww*(1.5*tau+sqrt(3.)+1./taup)*(alpha*r*dokdr)/(8.*aa)
768  * /(rg/cc)
769 cc wwfalse=-8.*AA*(TE**4)*TAUP/(ALPHA*R*DOKDR)*(RG/CC)
770 cc te4=-ww/TAUP*(ALPHA*R*DOKDR)/(8.*AA)/(RG/CC)
771  te=te4**0.25
772  if (wwfalse.gt.wwtrue) t1=te
773  if (wwfalse.lt.wwtrue) t2=te
774  if (abs(wwfalse-wwtrue)/wwtrue.lt.wweps) goto 100
775 c if (abs(wwfalse-wwtrue).ge.1.e4) goto 100
776 1 continue
777 C --------------------------------------------
778 100 continue
779  eta = ain*de*vr*vr/(ainp1*pr)*cc*cc
780  dl = 0.
781  ss = 2*ain*de*hh*rg
782 C -----------------------------
783 c if (imore.eq.2) print *,"PR=",pr," HH=",hh," WW=",ww
784 c if (imore.eq.2) WRITE (6,400) DE,VR,DL,TE,SS
785  demn = ain*de
786  temn = 2*te/3.0
787  xkff = xkkr*demn*(temn**(-3.5))
788  tauff = xkff*de*hh*rg
789  taues = xkes*de*hh*rg
790  tau = taues+tauff
791  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
792  taup = qbr/(8.*aa*cc*(te**4))
793  taueff = sqrt(taup*(tau+2./sqrt(3.)))
794  beta=1.
795 
796  vx = -vr
797  vp = dl/r
798 C
799  RETURN
800  END
801 C****************************************
802  SUBROUTINE intgrt (R,ETA,DL,IC)
803 C****************************************
804 C
805 C ------------------------------------------------------
806 C NUMERICAL INTEGRATION OF STEADY DISK EQUATIONS
807 C BY GEAR'S METHOD
808 C ------------------------------------------------------
809 C
810  include 'para.h'
811 C
812  dimension yy(2)
813  dimension f(2),xx(jmax)
814  dimension yscal(2)
815  dimension oldyy(2)
816  EXTERNAL fun,jac
817 C
818 c print *,"INTGRT"
819  IF ((iprint(1) .NE. 0).and.(imore.eq.2)) WRITE (6,300)
820 C
821 c ----------------------------------------------------------
822 c htry=hini : first dr for steady solution (t=0)
823 c ya1,ya2 : eta, dl
824 c ----------------------------------------------------------
825 c
826  eps = 1.e-8
827  x0(1) = r
828  y0(1,1) = eta
829  y0(2,1) = dl
830  yy(1) = eta
831  yy(2) = dl
832  nv=2
833  htry=hini
834  ya1=abs(yy(1))
835  ya2=abs(yy(2))
836  if (ya1.ge.eps) yscal(1)=ya1
837  if (ya1.lt.eps) yscal(1)=eps
838  if (ya2.ge.eps) yscal(2)=ya2
839  if (ya2.lt.eps) yscal(2)=eps
840 c yscal(2)=abserr
841 C ---- MESH INTERVAL ----
842 c aldr = (log10(ROUT1)-log10(RIN))/(mesh-50.)
843  aldr = (log10(rout1)-log10(rin))/mesh
844  xx(1)= rout1
845  xold = rout1
846  jj=2
847 10 continue
848 c ++++ Mesh interval 2 ++++
849 c XX(JJ)=XOLD/(10.**aldr)
850 c ---- For r < 1000rg case -----
851 c IF (XOLD.GT.1000.) XX(JJ)=XOLD/(10.**(aldr/3.))
852  IF (xold.GT.1000.) xx(jj)=xold/(10.**(aldr/1.))
853  IF (xold.LE.1000.) xx(jj)=xold/(10.**aldr/1.)
854 c ---- For r < 3.2rg case -----
855  dxx = xx(jj-1)-xx(jj)
856 c write(6,*) "before",jj,xx(jj),dxx
857 c pause
858 c IF (dxx.LE.0.05) XX(JJ)=XOLD-0.01
859 c IF (dxx.lt.0.5) XX(JJ)=XOLD-0.01
860 c IF (XOLD.LE.3.0) XX(JJ)=XOLD-0.005
861  IF (dxx.LE.0.3) xx(jj)=xold-0.3*dxx
862 cc IF (XOLD.LE.4.0) xx(jj)=xold-0.3*dxx
863 c IF (XOLD.LE.3.5) xx(jj)=xold-0.01
864 c IF (dxx.LE.0.1) XX(JJ)=XOLD-0.01
865  if (xold .le. 3.0) xx(jj)=xold-0.01
866 c dxx = xx(jj-1)-xx(jj)
867 
868 c if (xx(jj).lt.3.5) then
869 c write(6,*) "after ",jj,xx(jj),dxx
870 c pause
871 c end if
872 
873 c +++++++++++++++++++++++++
874  IF (xx(jj).LT.rin) GOTO 11
875  xold = xx(jj)
876  jj = jj+1
877  GOTO 10
878 11 CONTINUE
879 c ---- Mesh (in real) ----
880  jjx = jj
881 C ---- Start Integration at each radius ----
882 c
883  do 2000 j=1,jjx-1
884 C
885  r = x0(j)
886  xend = xx(j+1)
887 c ---- IOD : ITR count at each grid. ----
888  iod=0
889 C
890  500 continue
891 C print *,j,xx(j)
892  iod=iod+1
893 c 600 continue
894  oldr=r
895  oldhtry=htry
896  oldyy(1)=yy(1)
897  oldyy(2)=yy(2)
898 C ----------------------------------------------------------------------
899  call fun (r,yy,f)
900 c input => r, eta, dl == r, yy(1),yy(2)
901 c output => dy dx == f
902 c for
903 C ----------------------------------------------------------------------
904 C ---------------------------------------------------------------------
905  call stifbs (yy,f,nv,r,htry,relerr,yscal,hdid,hnext,fun,icon)
906 c : To solve for stiff problem.
907 C ----------------------------------------------------------------------
908 C
909 c ----------- CHECK -----------------------------------------------
910 c call fun (r,yy,f)
911 c if (kpr.eq.1) then
912 c r=oldr
913 c htry=oldhtry/2.
914 c yy(1)=oldyy(1)
915 c yy(2)=oldyy(2)
916 c goto 600
917 c endif
918 c ---------------------------------------------------------------------
919  IF ((icon .EQ. 10000) .OR. (icon .EQ. 15000)) THEN
920  if (imore.eq.2) WRITE (6,100) r,eps
921  ELSE IF (icon .EQ. 16000) THEN
922  if (imore.eq.2) WRITE (6,120)
923  ic = -1
924  RETURN
925  ELSE IF (icon .EQ. 30000) THEN
926  if (imore.eq.2) WRITE (6,140)
927  stop
928  ENDIF
929  100 FORMAT (1h ,'TOLERANCE RESET. R = ',1pe15.7,5x,'EPSR = ',1pe15
930  * .7)
931  120 FORMAT (1h ,'NO CONVERGENCE IN CORRECTOR ITERATION')
932  140 FORMAT (1h ,'INVALID INPUT TO ODGE')
933  IF ((mod(iod,1000).EQ.0).and.(imore.eq.2))
934  & WRITE (6,*) r,yy(1),yy(2)
935 C
936 C -- OBTAIN RSTEP --
937  if (r.eq.xend) then ! xend=xx(j+1)
938  rstep = r-xx(j+2)
939  else
940  rstep = r-xx(j+1)
941  endif
942  htry = -min(abs(hnext),rstep)
943  ya1 = abs(yy(1))
944  ya2 = abs(yy(2))
945  if (ya1.ge.eps) yscal(1) = ya1
946  if (ya1.lt.eps) yscal(1) = eps
947 c yscal(1)=max(abs(yy(1)),1.e-8)
948 c yscal(2)=max(abs(yy(2)),1.e-8)
949  if (ya2.ge.eps) yscal(2) = ya2
950  if (ya2.lt.eps) yscal(2) = eps
951 c yscal(2)=abserr
952  if (r.gt.xend) goto 500
953 C
954  x0(j+1) = r
955  y0(1,j+1) = yy(1)
956  y0(2,j+1) = yy(2)
957 C
958  ip1=iprint(1)
959  IF (ip1 .EQ. 0) GOTO 1000
960  IF (((mod((j+1),ip1) .EQ. 1) .OR. (ip1 .EQ. 1)).and.(imore.eq.2))
961  & WRITE (6,220) r,yy(1),yy(2),iod
962  220 FORMAT (1h ,3(1pe15.7),i7)
963 C
964  2000 continue
965  1000 CONTINUE
966 C
967  eta = yy(1)
968  IF (eta .GT. 1.) THEN
969  ic = 0
970  ELSE
971  ic = 1
972  ENDIF
973 C
974 c JX = J+1
975  jx = j
976 C
977  300 format (7x,"R",13x,"ETA",13x,"DL")
978  RETURN
979  END
980 C
981 C
982 C*******************************************************
983  SUBROUTINE fun (XI,YI,YP)
984 c input : XI=r,YI(1)=eta,YI(2)=dl
985 c output : YP(1)=d(eta)/dr, YP(2)=d(l)/dr
986 C*******************************************************
987 C
988  include 'para.h'
989 C
990  dimension yi(2),yp(2)
991 c
992  cc2=cc*cc
993 c ----- rg/c^2 -----
994  unit = rg/cc2
995 C
996  eta = yi(1)
997 c ---- ETA must be > 0 ----
998  IF (eta .LT. 0.) THEN
999  kpr=1
1000  return
1001  ENDIF
1002 c -------------------------
1003  dl = yi(2)
1004  r = xi
1005  r2 = r*r
1006  r3 = r2*r
1007  r4 = r3*r
1008  rm1 = r-1.
1009  al2 = alpha**2
1010  ok = (r/rm1)*sqrt(.5/r3)
1011  xlk = r2*ok
1012  dlnok =-.5/r-1./rm1
1013 c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1014 c -d Ok dlnLk
1015 c DLNRH = -- ln (---) ; DLKDR = -----
1016 c dr r dr
1017 c
1018 c XL : dL-Lk = L (angular momentum at this radius)
1019 c XLL : L-Lin
1020 c XLLK2 : 2*DL*LK+(DL)^2 => (LK+DL)^2 - LK^2 = L^2 - LK^2
1021 c WS : Height integrated Pressure/Surface Density
1022 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1023  dlnrh = 1./r-dlnok
1024  dlkdr = xlk*(1.5/r-1./rm1)
1025  xl = dl+xlk
1026  xll = xl-xlin
1027  xll2 = xll*xll
1028  xll3 = xll2*xll
1029  xll7 = xll**7
1030  eta2 = eta*eta
1031  eta5 = eta**5
1032  xllk2 = 2.*dl*xlk+dl*dl
1033  ws = xll2*eta/(al2*r2)
1034 
1035 c +++++++++++++++++++++++++++++++
1036 c PD : Pressure(P0)/Density(rho0) ratio at equatorial plane(z=0).
1037 c HH : Scale height at z=0.
1038 c PR : P0 (Pressure at z=0)
1039 c DE : rho (Density at z=0)
1040 c +++++++++++++++++++++++++++++++
1041 
1042  pd = (ain/ainp1)*ws
1043  hh = sqrt(2.*xnpl1*pd)/ok
1044  pr = (xmdot*xll/(4.*pi*r2*alpha*ainp1*hh))*(cc/(rg*rg))
1045 
1046 c +++ For Error Check +++
1047  if (pr.le.0.) then
1048  kpr=1
1049  return
1050  end if
1051 c +++++++++++++++++++++++
1052  de = pr/(pd*cc*cc)
1053 c -----------------------------------------------------------------
1054 c If optically thick or r=rout1 ==> TE derived from P = Pgas+Prad
1055 c If optically thin ==> TE derived from P = Pgas
1056 c -----------------------------------------------------------------
1057  if ((kr.eq.1).or.(r.eq.rout1)) then
1058  te = temp(de,pr)
1059 c TE = temp2(DE,PR)
1060  else
1061  te = pr/de/rmu
1062  endif
1063  beta = rmu*de*te/pr
1064  beta4 = beta**4
1065  demn = ain*de
1066  temn = 2.*te/3.
1067  xkff = xkkr*demn*(temn**(-3.5))
1068  xk = xkes+xkff
1069  xkxe = xk/xkes
1070  xkxf = xkff/xk
1071 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1072 c QBR : Brems cooling [erg/s/cm^2] (BHAD P248 eq.(8.64))
1073 c taup : Qbrems/Qradiation(?)
1074 c B10 : 1-beta
1075 c B11 : 4-3*beta
1076 c B12 : beta+12*(gamma-1)*(1-beta)
1077 c !GG1 and GG3 are the specific heats ratio
1078 c GG1 : Gamma1
1079 c GG3 : Gamma3
1080 c
1081 c A1 : 2*Gamma1
1082 c A2 : Gamma1-1
1083 c A3 : Gamma3-1
1084 c A4 : 3*Gamma1-1
1085 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1086  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
1087  taup = qbr/(8.*aa*cc*(te**4))
1088  b10 = 1.-beta
1089  b11 = 4.-3.*beta
1090  b12 = beta+12.*gm1*b10
1091  gg1 = beta+(b11*b11)*gm1/b12
1092  gg3 = 1.+(gg1-beta)/b11
1093  a1 = 2.*gg1
1094  a2 = gg1-1.
1095  a3 = gg3-1.
1096  a4 = 3.*gg1-1.
1097 c ------------------------------------
1098 c QP : Q^+ (Energy generation)
1099 c QM : Q^-
1100 c FF : Q^- (Brems cooling)
1101 c C1LL : (L^2-Lk^2)/R^3 * SIGMA/PI
1102 c ------------------------------------
1103  qp = 2.*al2*xl/xll
1104  if ((kr.eq.1).or.(r.eq.rout1)) then
1105  qm = ff1*beta4*xll7*eta5/(xkxe*r4)
1106  else
1107  ff = qbr*rg/cc/cc/cc
1108  qm = 2.*pi*r2*eta/(xmdot*ws)*ff
1109  endif
1110  c1ll = xllk2/(r3*ws)
1111 C
1112 C ***** MATRIX ELEMENTS *****
1113 C
1114 c | D11 D12 | = | 1 (eta+1)/(L-Lin) |
1115 c | D21 D22 | = | (3*Gam3-1)/2 (2*Gam1*eta-(Gam3-1)*alpha**2)/(L-Lin) |
1116 c
1117 c | c1 | = (eta+1)/R + DLNRH + Sigma/Pi * (L^2-Lk^2)/R^3
1118 c | c2 | = eta*(2Gam1/R - (Gam1-1)DLNRH) - (Gam3-1)*(Q^+ - Q^-)/R)
1119 c
1120  d11 = 1.
1121  d12 = (eta+1.)/xll
1122  d21 = a4/2.
1123  d22 = (a1*eta-a3*al2)/xll
1124  c1 = (eta+1.)/r+dlnrh+c1ll
1125  c2 = eta*(a1/r-a2*dlnrh)-a3*(qp-qm)/r
1126  c1d = c1-d12*dlkdr
1127  c2d = c2-d22*dlkdr
1128 C
1129 c++++++++++++++++++++++++++++++++++
1130 c DD : D (Matsumoto et al. 1984 eq(3.13))
1131 c XN1D : N_1 eq(3.14)
1132 c XN2D : N_2 eq(3.14)
1133 c++++++++++++++++++++++++++++++++++
1134  dd = d11*d22-d12*d21
1135  xn1d = c1d*d22-c2d*d12
1136  xn2d = d11*c2d-d21*c1d
1137 C
1138 C ***** F1 AND F2 *****
1139 C
1140  yp(1) = xn1d/dd
1141  yp(2) = xn2d/dd
1142  kpr=0
1143 C
1144 C
1145  RETURN
1146  END
1147 C
1148 C
1149 C**************************************
1150  subroutine jac(xi,yi,yp,dfdr,pdj)
1151 C**************************************
1152 C
1153  include 'para.h'
1154 C
1155  dimension yi(2),pdj(2,2),dj(2,2,3),cj(2,3),cdj(2,3),dnj(2,3),
1156  & ddj(3),dfdr(2),yp(2)
1157  cc2=cc*cc
1158  unit = rg/cc2
1159 C
1160  eta = yi(1)
1161  IF (eta.LT.0.) THEN
1162  kpr=1
1163  return
1164  ENDIF
1165  dl = yi(2)
1166  r = xi
1167  r2 = r*r
1168  r3 = r2*r
1169  r4 = r2*r2
1170  r7 = r3*r3*r
1171  r8 = r7*r
1172  rm1 = r-1.
1173  al2 = alpha*alpha
1174 C
1175  ok = (r/rm1)*sqrt(.5/r3)
1176  xlk = r2*ok
1177  xl = xlk+dl
1178  xll = xl-xlin
1179  xll2 = xll*xll
1180  xll3 = xll**3
1181  xll7 = xll**7
1182  eta2 = eta*eta
1183  eta5 = eta**5
1184  xl2 = xl*xl
1185  ws = xll2*eta/(al2*r2)
1186  pd = (ain/ainp1)*ws
1187  xlk2 = xlk*xlk
1188  xllk2 = 2.*dl*xlk+dl*dl
1189  dlnok =-0.5/r-1./rm1
1190  dlnrh = 1./r-dlnok
1191  d2lnrh=-1.5/r2-1./(rm1**2)
1192  dlkdr = xlk*(1.5/r-1./rm1)
1193  hh = sqrt(2.*xnpl1*pd)/ok
1194  pr = (xmdot*xll/(4.*pi*r2*alpha*ainp1*hh))*(cc/(rg**2))
1195  if (pr.le.0.) then
1196  kpr=1
1197  return
1198  endif
1199  de = pr/(pd*cc*cc)
1200  if ((kr.eq.1).or.(r.eq.rout1)) then
1201  te = temp(de,pr)
1202  else
1203  te = pr/de/rmu
1204  endif
1205  beta = rmu*de*te/pr
1206  beta4 = beta**4
1207  demn = ain*de
1208  temn = 2.*te/3.
1209  xkff = xkkr*demn*(temn**(-3.5))
1210  xk = xkes+xkff
1211  xkxe = xk/xkes
1212  xkxf = xkff/xk
1213  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
1214  taup = qbr/(8.*aa*cc*(te**4))
1215  b10 = 1.-beta
1216  b11 = 4.-3.*beta
1217  b12 = beta+12.*gm1*b10
1218  gg1 = beta+(b11**2)*gm1/b12
1219  gg3 = 1.+(gg1-beta)/b11
1220  a1 = 2.*gg1
1221  a2 = gg1-1.
1222  a3 = gg3-1.
1223  a4 = 3.*gg1-1.
1224  qp = 2.*al2*xl/xll
1225  if ((kr.eq.1).or.(r.eq.rout1)) then
1226  qm = ff1*beta4*xll7*eta5/(xkxe*r4)
1227  else
1228  ff = qbr*rg/cc/cc/cc
1229  qm = 2.*pi*r2*eta/(xmdot*ws)*ff
1230  endif
1231  c1ll = xllk2/(r3*ws)
1232 C
1233 C ***** MATRIX ELEMENTS *****
1234 C
1235  d11 = 1.
1236  d12 = (eta+1.)/xll
1237  d21 = a4/2.
1238  d22 = (a1*eta-a3*al2)/xll
1239  c1 = (eta+1.)/r+dlnrh+c1ll
1240  c2 = eta*(a1/r-a2*dlnrh)-a3*(qp-qm)/r
1241  c1d = c1-d12*dlkdr
1242  c2d = c2-d22*dlkdr
1243  dd = d11*d22-d12*d21
1244  xn1d= c1d*d22-c2d*d12
1245  xn2d= d11*c2d-d21*c1d
1246 C
1247 C ***** CAL. JACOBIAN *****
1248 C
1249  b13 = beta*b10/b11
1250  dbdl =-b13*(8./xll)
1251  dbde =-b13*(4.5/eta)
1252  dbdr = b13*(7./r)
1253  if ((kr.eq.1).or.(r.eq.rout1)) then
1254  dg1db = 1.-b11*gm1*(3.*beta+4.+12.*gm1*(2.-3.*beta))/(b12*b12)
1255  dg3db = (dg1db-1.)/b11+3.*(gg1-beta)/(b11*b11)
1256  else
1257  dg1db = 0.
1258  dg3db = 0.
1259  endif
1260  dpdl = 0.
1261  dpde =-0.5/eta
1262  dpdr =-1./r
1263  dddl =-2./xll
1264  ddde =-1.5/eta
1265  dddr = 1./r
1266  dwsde = dpde-ddde
1267  dwsdl = dpdl-dddl
1268  dwsdr = dpdr-dddr
1269  dtdl = (dpdl-beta*dddl)/b11
1270  dtde = (dpde-beta*ddde)/b11
1271  dtdr = (dpdr-brta*dddr)/b11
1272  dtdr =-2./r
1273  dkdl = xkxf*(dddl-3.5*dtdl)
1274  dkde = xkxf*(ddde-3.5*dtde)
1275  dkdr = xkxf*(dddr-3.5*dtdr)
1276  if ((kr.eq.1).or.(r.eq.rout1)) then
1277  dqmde = qm*(4.*dbde/beta+5./eta-dkde)
1278  dqmdl = qm*(4.*dbdl/beta+7./xll-dkdl)
1279  dqmdr = qm*(4.*dbdr/beta-4./r -dkdr)
1280  else
1281  dqmde = qm*(-1./eta)
1282  dqmdl = qm*(-2./xll)
1283  dqmdr = qm*(0.5*(1.-3./r)/rm1)
1284  endif
1285  dqpdl = qp*(1./xl-1./xll)
1286  dqpdr = 0.
1287  d21db = 1.5*dg1db
1288  d22db = (2.*dg1db*eta-dg3db*al2)/xll
1289  dc2db = eta*(2.*dg1db/r-dg1db*dlnrh)-dg3db*(qp-qm)/r
1290 C
1291 C
1292  dj(1,1,1) = 0.
1293  dj(1,1,2) = 0.
1294  dj(1,1,3) = 0.
1295  dj(1,2,1) = 1./xll
1296  dj(1,2,2) =-d12/xll
1297  dj(1,2,3) = 0.
1298  dj(2,1,1) = d21db*dbde
1299  dj(2,1,2) = d21db*dbdl
1300  dj(2,1,3) = d21db*dbdr
1301  dj(2,2,1) = d22db*dbde+2.*gg1/xll
1302  dj(2,2,2) = d22db*dbdl-d22/xll
1303  dj(2,2,3) = d22db*dbdr
1304 C
1305  cj(1,1) = 1./r-c1ll*dwsde
1306  cj(1,2) = 2.*xl/(r3*ws)-c1ll*dwsdl
1307  cj(1,3) =-(eta+2.5)/r2-(3./r)*c1ll-1./(rm1*rm1)
1308  cj(2,1) = a1/r-a2*dlnrh+a3*dqmde/r+dc2db*dbde
1309  cj(2,2) =-a3*(dqpdl-dqmdl)/r+dc2db*dbdl
1310  cj(2,3) =-eta*(a1-1.5*a2)/r2+eta*a2/sqrt(rm1)
1311  & +a3*(dqmdr*r+qp-qm)/r2+dc2db*dbdr
1312 C
1313  cdj(1,1) = cj(1,1)-dj(1,2,1)*dlkdr
1314  cdj(1,2) = cj(1,2)-dj(1,2,2)*dlkdr
1315  cdj(1,3) = cj(1,3)-dj(1,2,3)*dlkdr
1316  & -d12*(dlkdr*(1.5/r-1./rm1)+xlk*(-1.5/r2+1./sqrt(rm1)))
1317  cdj(2,1) = cj(2,1)-dj(2,2,1)*dlkdr
1318  cdj(2,2) = cj(2,2)-dj(2,2,2)*dlkdr
1319  cdj(2,3) = cj(2,3)-dj(2,2,3)*dlkdr
1320  & -d22*(dlkdr*(1.5/r-1./rm1)+xlk*(-1.5/r2+1./sqrt(rm1)))
1321 C
1322  dnj(1,1) = cdj(1,1)*d22+c1d*dj(2,2,1)-cdj(2,1)*d12-c2d*dj(1,2,1)
1323  dnj(1,2) = cdj(1,2)*d22+c1d*dj(2,2,2)-cdj(2,2)*d12-c2d*dj(1,2,2)
1324  dnj(1,3) = cdj(1,3)*d22+c1d*dj(2,2,3)-cdj(2,3)*d12-c2d*dj(1,2,3)
1325  dnj(2,1) = dj(1,1,1)*c2d+d11*cdj(2,1)-dj(2,1,1)*c1d-d21*cdj(1,1)
1326  dnj(2,2) = dj(1,1,2)*c2d+d11*cdj(2,2)-dj(2,1,2)*c1d-d21*cdj(1,2)
1327  dnj(2,3) = dj(1,1,3)*c2d+d11*cdj(2,3)-dj(2,1,3)*c1d-d21*cdj(1,3)
1328 C
1329  ddj(1) = dj(1,1,1)*d22+d11*dj(2,2,1)-dj(1,2,1)*d21-d12*dj(2,1,1)
1330  ddj(2) = dj(1,1,2)*d22+d11*dj(2,2,2)-dj(1,2,2)*d21-d12*dj(2,1,2)
1331  ddj(3) = dj(1,1,3)*d22+d11*dj(2,2,3)-dj(1,2,3)*d21-d12*dj(2,1,3)
1332 C
1333  dfdr(1) = (dnj(1,3)*dd-xn1d*ddj(3))/(dd*dd)
1334  dfdr(2) = (dnj(2,3)*dd-xn2d*ddj(3))/(dd*dd)
1335  pdj(1,1) = (dnj(1,1)*dd-xn1d*ddj(1))/(dd*dd)
1336  pdj(1,2) = (dnj(1,2)*dd-xn1d*ddj(2))/(dd*dd)
1337  pdj(2,1) = (dnj(2,1)*dd-xn2d*ddj(1))/(dd*dd)
1338  pdj(2,2) = (dnj(2,2)*dd-xn2d*ddj(2))/(dd*dd)
1339  kpr=0
1340 C
1341  RETURN
1342  END
1343 C
1344 C
1345 C*********************************************
1346  DOUBLE PRECISION FUNCTION temp(DE,PR)
1347 C*********************************************
1348 C
1349  include 'para.h'
1350 C
1351 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1352 c A1 : 3/a*(RMU*rho0)
1353 c A2 : -3/a*P0 (=T0^4)
1354 c TMAX : Temperature for Ptotal=Prad (P0=a/3*T0^4)
1355 c TMIN : 0.
1356 c
1357 c FX : TMAX^4 + A1*TMAX + A2 == 3/a*(Prad + Pgas - Ptotal) > 0
1358 c FN : TMIN^4 + A1*TMIN + A2 == < 0
1359 c
1360 c TH : Mean Temperature
1361 c FH : TH^4 + A1*TH + A2 ==> Basically, this quantity = 0 !
1362 c
1363 c if T > Treal -> F(T) > 0
1364 c if T < Treal -> F(T) < 0
1365 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1366  a1 = 3.*rmu*de/aa
1367  a2 =-3.*pr/aa
1368  eps = 1.0e-12
1369 c eps = 1.0e-10
1370  tmax= (-a2)**0.25
1371  tmin= 0.
1372  fx = tmax**4+a1*tmax+a2
1373  fn = tmin**4+a1*tmin+a2
1374  DO 40 itr=1,10
1375  th = 0.5*(tmax+tmin)
1376  fh = th**4+a1*th+a2
1377 c
1378  IF (fh*fn .LT. 0.) tmax=th
1379  IF (fh*fn .LT. 0.) fx=fh
1380 c
1381  IF (fh*fx .LT. 0.) tmin=th
1382  IF (fh*fx .LT. 0.) fn=fh
1383  40 CONTINUE
1384 C
1385  xxx = th
1386 C
1387 c +++++ For obtain the solutiono of TEMPERATURE using Newton Method ++++
1388 c
1389 c xxx : Rough solution of TE
1390 c xx1 : Better solution of TE
1391 c d(Prad+Pgas)
1392 c xx1 = xxx - ------------
1393 c dT
1394 c TEMP: Final solution of TE
1395 c
1396  DO 100 itr=1,20
1397  xx1 = xxx-(xxx**4+a1*xxx+a2)/(4.*xxx**3+a1)
1398  IF (abs((xx1-xxx)/xx1) .LT. eps) GO TO 120
1399  xxx = xx1
1400  100 CONTINUE
1401 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1402 c
1403 c +++++ For not converge ! +++++
1404  WRITE (6,110)
1405  print *,de,pr,tmax,tmin,fx,fn,th,fh,xxx,xx1
1406 110 FORMAT (1h ,'LOOP COUNT (TEMP) IS OVER MAXIMUM')
1407  stop
1408 c ++++++++++++++++++++++++++++++
1409 c
1410 120 CONTINUE
1411 c
1412  temp=xx1
1413 C
1414  RETURN
1415  END
1416 C
1417 C
1418 C******************************
1419  DOUBLE PRECISION FUNCTION ondo(DE,PR,taup,tau)
1420 C******************************
1421 C
1422  include 'para.h'
1423 C
1424  a1 = 3.*rmu*de/aa*(1.+1./taup/(1.5*tau+sqrt(3.)))
1425  a2 =-3.*pr/aa*(1.+1./taup/(1.5*tau+sqrt(3.)))
1426  eps = 1.0e-12
1427  tmax= (-a2)**0.25
1428  tmin= 0.
1429  fx = tmax**4+a1*tmax+a2
1430  fn = tmin**4+a1*tmin+a2
1431  DO 40 itr=1,10
1432  th = 0.5*(tmax+tmin)
1433  fh = th**4+a1*th+a2
1434  IF (fh*fn .LT. 0.) tmax=th
1435  IF (fh*fn .LT. 0.) fx=fh
1436  IF (fh*fx .LT. 0.) tmin=th
1437  IF (fh*fx .LT. 0.) fn=fh
1438  40 CONTINUE
1439 C
1440  xxx = th
1441 C
1442  DO 100 itr=1,20
1443  xx1 = xxx-(xxx**4+a1*xxx+a2)/(4.*xxx**3+a1)
1444  IF (abs((xx1-xxx)/xx1) .LT. eps) GO TO 120
1445  xxx = xx1
1446  100 CONTINUE
1447 C
1448  WRITE (6,110)
1449  print *,de,pr,tmax,tmin,fx,fn,th,fh,xxx,xx1
1450  110 FORMAT (1h ,'LOOP COUNT (TEMP) IS OVER MAXIMUM')
1451  stop
1452 C
1453  120 CONTINUE
1454 C
1455  ondo=xx1
1456 C
1457  RETURN
1458  END
1459 C
1460 C***************************
1461  DOUBLE PRECISION FUNCTION fnin(NPL)
1462 C***************************
1463 C
1464  IMPLICIT REAL*8 (a-h,o-z)
1465  fnin = 1
1466 C
1467  DO 10 i=1,npl
1468  odd = 2.*i+1
1469  evn = 2.*i
1470  fnin= fnin*evn/odd
1471 10 CONTINUE
1472 C
1473  RETURN
1474  END
1475 C
1476 C***************************
1477  DOUBLE PRECISION FUNCTION fnjn(NPL)
1478 C***************************
1479 C
1480  IMPLICIT REAL*8 (a-h,o-z)
1481  fnjn = 3.1415926/2.0d0
1482 C
1483  DO 10 i=1,npl+1
1484  odd = 2.*i-1
1485  evn = 2.*i
1486  fnjn= fnjn*odd/evn
1487 10 CONTINUE
1488 C
1489  RETURN
1490  END
1491 C
1492 C
1493 C**************************
1494  SUBROUTINE cnvrbs
1495 C**************************
1496 C
1497 C -----------------------------------------------------------
1498 C TRANSFORMATION FROM (Y0:ETA,DL) TO (W0:DE,VR,DL,TE)
1499 C SET OUTER BOUNDARY AT ROUT2
1500 C -----------------------------------------------------------
1501 C
1502  include 'para.h'
1503 C
1504 c write(6,*) "sub-CNVRBS"
1505  j=0
1506 C
1507  200 continue
1508  j = j+1
1509  jj = jx-j+1
1510  r = x0(jj)
1511  if (r.le.rin2) goto 200
1512 100 CONTINUE
1513  j3 = jj
1514  j = 0
1515 C
1516  2000 continue
1517 C
1518  j = j+1
1519  jj = j3-j+1
1520  r = x0(jj)
1521  eta = y0(1,jj)
1522  IF (eta.LT.0.) THEN
1523  kpr=1
1524  return
1525  ENDIF
1526  dl = y0(2,jj)
1527  rm1 = r-1.
1528  r2 = r*r
1529  r3 = r2*r
1530  ok = (r/rm1)*sqrt(.5/r3)
1531  xlk = r2*ok
1532  xl = xlk+dl
1533  xll = xl-xlin
1534  xll2= xll*xll
1535  al2 = alpha*alpha
1536  ws = xll2*eta/(al2*r2)
1537  pd = (ain/ainp1)*ws
1538  hh = (sqrt(2.*xnpl1*pd))/ok
1539  pr = (xmdot*xll/(4.*pi*r2*alpha*ainp1*hh))*cc/(rg*rg)
1540  if (pr.le.0.) then
1541  kpr=1
1542  return
1543  endif
1544  de = pr/(pd*cc*cc)
1545  tauff = xkff*de*hh*rg
1546  taues = xkes*de*hh*rg
1547  tau = taues+tauff
1548  te = temp(de,pr)
1549 
1550  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
1551  taup = qbr/(8.*aa*cc*(te**4))
1552 
1553 c write(6,*) r,tau,taup
1554 c pause
1555 
1556  if ((kr.eq.1).or.(r.eq.rout1)) then
1557 c te = temp(de,pr)
1558 c if (kcool.eq.0) TE = TEMP(DE,PR)
1559  if (kcool.eq.1) te = ondo(de,pr,taup,tau)
1560  else
1561  te = pr/de/rmu
1562  endif
1563  s0 = (4.*aa/(3.*de))*(te**3)+rmu*((1./gm1)*log(te)-log(de))
1564  beta= rmu*de*te/pr
1565  a = 3*(1-beta)+beta/gm1
1566  ss = 2*ain*de*hh*rg
1567  ww = 2*ainp1*hh*pr*rg
1568  wc = ww/(cc*cc)
1569  ws = wc/ss
1570  vx =-xmdot/(2.*pi*r*ss*cc*rg)
1571  vp = dl/r
1572  ee = (a+0.5)*ws+0.5*(vx*vx+vp*vp)
1573 C
1574  w0(j,1) = de
1575  w0(j,2) = vx
1576  w0(j,3) = vp
1577  w0(j,4) = te
1578  x1(j) = r
1579  rs = r*ss
1580  u0(j,1) = rs
1581  u0(j,2) = rs*vx
1582  u0(j,3) = rs*vp
1583  u0(j,4) = rs*ee
1584 C
1585  demn = ain*de
1586  temn = 2*te/3.0
1587  xkff = xkkr*demn*(temn**(-3.5))
1588 c TAUES = 2.*XKES*DE*HH*RG
1589  taues = xkes*de*hh*rg
1590 c TAUFF = 2.*XKFF*DE*HH*RG
1591  tauff = xkff*de*hh*rg
1592  tau = taues+tauff
1593  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
1594  taup = qbr/(8.*aa*cc*(te**4))
1595  tauabso(j)= taup
1596  tauo(j)=tau
1597 C
1598  1 format ("w0=",4e11.3)
1599  2 format ("x1=",f7.3," rs=",e11.3)
1600  3 format ("u0=",4e11.3)
1601 C
1602  if (r.lt.rout2) goto 2000
1603 1000 CONTINUE
1604 C
1605  jx = j
1606  if (imore.eq.2)
1607  & WRITE (6,*) ' ### MESH POINTS = ',jx,' ###'
1608 C
1609  RETURN
1610  END
1611 C
1612 C
1613 C**************************
1614  SUBROUTINE lw1d
1615 C**************************
1616 C
1617 C ----------------------------------------------------------------------
1618 C TIME INTEGRATION OF HYDRODYNAMIC EQUATIONS OF THIN ACCRETION
1619 C DISK BY LAX-WENDROFF METHOD.
1620 C ----------------------------------------------------------------------
1621 C
1622  include 'para.h'
1623 C
1624  dimension u(jmax,4),f(jmax,4),h(jmax,4),y(jmax,4),q1(jmax),
1625  *q2(jmax),d(jmax,4)
1626  dimension wk(jmax,4)
1627  dimension drp(jmax),drm(jmax),dr(jmax)
1628  dimension betaj(jmax),tej(jmax)
1629  dimension wt(jmax,4),tauabst(jmax),taut(jmax)
1630 c dimension dexise(jmax),uo(jmax)
1631  dimension okc(jmax),hhc(jmax),prc(jmax),tec(jmax),betac(jmax),
1632  & ssc(jmax),dec(jmax),betao(jmax)
1633  dimension des(jmax),oes1(jmax),ou(jmax,4)
1634 C
1635 C
1636  gm1 = gm-1
1637  cc2 = cc*cc
1638  unit = rg/cc2
1639 C
1640  DO 5 j=2,jx-1
1641  drp(j) = x1(j+1)-x1(j) ! 前進差分
1642  drm(j) = x1(j)-x1(j-1) ! 後退差分
1643  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,i12,2x,a3,1pe9.3,2x,a3,1pe9.3,2x,a3,1pe9.3) end if end if c IF (MOD(NS,KLC).EQ.0) WRITE (9,900) TIME,xlumi,xltn,xltk,dt,ns RETURN c----------------------- endif c----------------------- cc900 FORMAT (1PE16.8,1X,1PE15.7,1X,1PE15.7,1X,1PE15.7,1x,1pe15.7) 900 FORMAT (1PE10.4,2X,1PE10.4,2X,1PE10.4,2X,1PE10.4,2x,1pe10.4,2x,i9) 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=50.0 c xlamda=pxlin xlamda=1.0 rp=1.0 do 10 j=1,jx-1 do 20 k=1,4 w0(j,k)=w0(j,k) & +rp*w0(j,k)*exp(-(x1(j)-rc)*(x1(j)-rc)/xlamda/xlamda) c & -rp*w0(j,k)*exp(-(x1(j)-rc)*(x1(j)-rc)/xlamda/xlamda) 20 continue 10 continue c write(6,*) "PERTRB 3" return 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 c do 20 k=1,1 c dss = 0.01*exp(-(x1(j)-rc)**2/0.1**2)*sin(xkk*x1(j))*w0(j,k) c w0(j,k)=w0(j,k)+dss 20 continue dww = dss/xkk/x1(j)*w0(j,1) dvr = dss/sqrt(xkk*x1(j))*w0(j,2) dvp = dss/sqrt(xkk*x1(j))*w0(j,3) 10 continue write(6,*) "PERTRB 4" return 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
1644 5 CONTINUE
1645  drp(1) = x1(2) -x1(1)
1646  drm(jx)= x1(jx)-x1(jx-1)
1647  dr(1) = drp(1)
1648  dr(jx) = drm(jx)
1649 
1650  do 152 j=1,jx-1
1651  nl(j)=0
1652 152 continue
1653 c
1654 151 continue
1655  xlumi=0. ! Total luminosity
1656  xltn=0. ! optically thin luminosity
1657  xltk=0. ! optically thick luminosity
1658 
1659 c print *,W0(1,1)
1660 c print *,"L1"
1661 
1662  DO 10 j=1,jx
1663  r = x1(j)
1664  rm1 = r-1
1665  r2 = r*r
1666  r3 = r2*r
1667  ok = (r/rm1)*sqrt(0.5/r3)
1668  de = w0(j,1)
1669  vx = w0(j,2)
1670  vp = w0(j,3)
1671  te = w0(j,4)
1672  wt(j,1) = w0(j,1)
1673  wt(j,2) = w0(j,2)
1674  wt(j,3) = w0(j,3)
1675  wt(j,4) = w0(j,4)
1676  tauabst(j) = tauabso(j)
1677  taut(j) = tauo(j)
1678  tau = tauo(j)
1679  taup = tauabso(j)
1680 c
1681  if (j.eq.jx) then
1682  pr = (aa/3.)*(te**4)+rmu*de*te
1683  else
1684  if (kcool.eq.1) fac = 1./((1.+1./taup/(1.5*tau+sqrt(3.))))
1685  if (kcool.eq.0) fac = 1.
1686 c
1687  pr = (aa/3.)*(te**4)*fac+rmu*de*te
1688  endif
1689 c
1690  beta = rmu*de*te/pr
1691  hh = sqrt(b5*pr/de)/ok
1692  ss = 2*ain*de*hh*rg
1693  ww = 2*ainp1*pr*hh*rg
1694  wc = ww/cc2
1695  ws = wc/ss
1696  a = 3*(1-beta)+beta/gm1
1697  ee = (a+0.5)*ws+0.5*(vx*vx+vp*vp) ! total energy
1698  rs = r*ss
1699  betaj(j)= beta
1700 
1701  u(j,1) = rs ! d(r*SS)/dt
1702  u(j,2) = rs*vx ! d(r*SS*vr)/dt
1703  u(j,3) = rs*vp ! d(r*SS*vp)/dt
1704  u(j,4) = rs*ee ! d(r*SS*ee)/dt
1705 
1706  f(j,1) = rs*vx ! d(r*SS*vr)/dr
1707  f(j,2) = rs*(vx*vx+ws) ! d(r*SS*vr^2+r*WW)/dr
1708  f(j,3) = rs*(vx*vp+alpha*ws) ! d(r^2*SS*vr*vp-Trp)/dr
1709  f(j,4) = rs*(vx*(ee+ws)+alpha*ws*vp) ! d()/dr
1710 
1711  vk = r*ok
1712  xlk = r*vk
1713  dlnor =-(0.5+r/rm1)
1714  dlnlr = 2+dlnor
1715  dokdr = ok*dlnor/r
1716  dlkdr = xlk*dlnlr/r
1717  demn = ain*de
1718  temn = 2*te/3.0
1719  xkff = xkkr*demn*(temn**(-3.5))
1720  tauff = xkff*de*hh*rg
1721  taues = xkes*de*hh*rg
1722  tau = taues+tauff
1723  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
1724  taup = qbr/(8.*aa*cc*(te**4))
1725 c
1726  if (kcool .eq. 1) then
1727  taueff = sqrt(taup*(tau+2./sqrt(3.)))
1728  else if (kcool .eq. 0) then
1729  taueff = sqrt(3*tau*tauff)
1730  end if
1731 c
1732  if (ns.eq.lns+1) then
1733 c --------- oes1(j): total energy = E - 0.5*(vr^2*vp^2)
1734  oes1(j)=u(j,4)/u(j,1)-((u(j,2)/u(j,1))*(u(j,2)/u(j,1))
1735  & +(u(j,3)/u(j,1))*(u(j,3)/u(j,1)))/2.
1736  ou(j,1)=u(j,1)
1737  ou(j,2)=u(j,2)
1738  ou(j,3)=u(j,3)
1739  ou(j,4)=u(j,4)
1740  tauo(j)=tau
1741  tauabso(j)=taup
1742  betao(j)=beta
1743  endif
1744 
1745  xk = xkes+xkff
1746  fd = 16*aa*(te**4)/(3*xk*de*hh*cc2)
1747  ff = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg*rg/cc2/cc
1748 c--- 2002-04-11 ---
1749  if (kcool .eq.1) then
1750  fm = 8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)*unit
1751  else
1752  fm = fd
1753  end if
1754 c
1755  h(j,1) = 0
1756  h(j,2) = wc*(1-dlnor)+ss*(2*vp*vk+vp*vp) ! (11.30)の右辺
1757  h(j,3) =-ss*(vx*vp+vx*dlkdr+alpha*ws)
1758  h(j,4) =-r*rs*(vx*vp+alpha*ws)*dokdr-r*fm ! (11.32)の右辺
1759  IF (ns.EQ.lns+1) THEN
1760 cc FFTR = FM/UNIT/XMCRIT/CC*16.
1761  fftr = fm/unit/xmcrit/cc
1762  xlumi = xlumi+2.*pi*(r*rg)*(dr(j)*rg)*fftr
1763  if (j.eq.jx) xljx=2.*pi*(r*rg)*(dr(j)*rg)*fftr
1764  if (taueff.le.1.) xltn=xltn+2.*pi*(r*rg)*(dr(j)*rg)*fftr
1765  if (taueff.gt.1.) xltk=xltk+2.*pi*(r*rg)*(dr(j)*rg)*fftr
1766  ENDIF
1767  betac(j) = beta
1768  okc(j) = ok
1769  hhc(j) = hh
1770  prc(j) = pr
1771  tec(j) = te
1772  ssc(j) = ss
1773  dec(j) = de
1774 10 CONTINUE
1775 c IF ((NS.EQ.LNS+1).and.(kcau1+kcau2.eq.0))
1776 c & WRITE (9,900) TIME-DT,XLumi,xltn,xltk,dt
1777 c print *,"L1.5"
1778 
1779 c###### Thermal Conduction ####################################
1780 c xkei: thermal conductivity => K=alpha*xkappa*C_s*Sigma*C_v*H
1781 c
1782  if (ktcond.eq.1) then
1783 c do 12 j=2,jx-1
1784  do 12 j=1,jx-1
1785 c ef=3.*xkappa*(8.-7.*betac(j))
1786 CCCc cv=3.*(8.-7.*betac(j))/(8.-3.*betac(j)-3.*betac(j))
1787  cv=3./2.
1788 cc xnu=alpha*okc(j)*(hhc(j)**2.)
1789 cc xkei=xkappa*ssc(j)*cv*xnu
1790 c gm31=gm1*(4.-3.*betac(j))
1791 c & /(betac(j)+12.*gm1*(1.-betac(j)))
1792 c ef=xkappa*gm31
1793  ef=xkappa
1794  xkei=alpha*ef*okc(j)*(hhc(j)**3)*prc(j)/tec(j)
1795  cn=cv*xkei*(tec(j+1)-tec(j))/drp(j)
1796  & *sqrt(2./xnpl1)*ainp1*unit
1797 c if (nl(j).ge.1) cn=0.
1798  f(j,4)=f(j,4)-x1(j)*cn
1799  h(j,4)=h(j,4)-cn
1800 12 continue
1801  endif
1802 C##########################################
1803 C
1804 c print *,"LW2"
1805  DO 40 j=1,jx-1
1806  ave1 = 0.5*(u(j+1,1)+u(j,1))
1807  ave2 = 0.5*(u(j+1,2)+u(j,2))
1808  ave3 = 0.5*(u(j+1,3)+u(j,3))
1809  ave4 = 0.5*(u(j+1,4)+u(j,4))
1810  fsa1 = f(j+1,1)-f(j,1)
1811  fsa2 = f(j+1,2)-f(j,2)
1812  fsa3 = f(j+1,3)-f(j,3)
1813  fsa4 = f(j+1,4)-f(j,4)
1814  swa1 = 0.5*(h(j+1,1)+h(j,1))
1815  swa2 = 0.5*(h(j+1,2)+h(j,2))
1816  swa3 = 0.5*(h(j+1,3)+h(j,3))
1817  swa4 = 0.5*(h(j+1,4)+h(j,4))
1818  dtdx = dt/drp(j)
1819  y(j,1) = ave1-0.5*dtdx*fsa1+0.5*dt*swa1
1820  y(j,2) = ave2-0.5*dtdx*fsa2+0.5*dt*swa2
1821  y(j,3) = ave3-0.5*dtdx*fsa3+0.5*dt*swa3
1822  y(j,4) = ave4-0.5*dtdx*fsa4+0.5*dt*swa4
1823  es1(j)=y(j,4)/y(j,1)-((y(j,2)/y(j,1))*(y(j,2)/y(j,1))
1824  & +(y(j,3)/y(j,1))*(y(j,3)/y(j,1)))/2.
1825 40 CONTINUE
1826  dtdx = dt/dr(1)
1827  fsa1 = f(2,1)-f(1,1)
1828  fsa2 = f(2,2)-f(1,2)
1829  fsa3 = f(2,3)-f(1,3)
1830  fsa4 = f(2,4)-f(1,4)
1831  swa1 = 0.5*(h(1,1)+h(2,1))
1832  swa2 = 0.5*(h(1,2)+h(2,2))
1833  swa3 = 0.5*(h(1,3)+h(2,3))
1834  swa4 = 0.5*(h(1,4)+h(2,4))
1835  wk(1,1) = u(1,1)-dtdx*fsa1+dt*swa1
1836  wk(1,2) = u(1,2)-dtdx*fsa2+dt*swa2
1837  wk(1,3) = u(1,3)-dtdx*fsa3+dt*swa3
1838  wk(1,4) = u(1,4)-dtdx*fsa4+dt*swa4
1839 C
1840 c##### Artificial Viscosity #####
1841 c What is kcau1 & kcau2 ? => Ans.: kcau1,2 = 0,or 1
1842 c if kcau1&2 = 0 -> no problem!
1843 c if kcau1&2 = 1 -> problem! -> modify!
1844 
1845  if (kcau1+kcau2.eq.0) then
1846  qu1=qu
1847  else
1848  qu1=16.
1849  endif
1850 
1851  DO 330 j = 1,jx-1
1852  vxp = w0(j+1,2)
1853  vxm = w0(j,2)
1854 c qu1=qu
1855 c if (nl(j).ge.1) qu1=qu*5.
1856  q2(j)= qu1*abs(vxp-vxm)
1857 330 CONTINUE
1858 c
1859  do 340 j=2,jx-1
1860  dtdx=dt/dr(j)
1861  des(j) = dtdx*(q2(j)*(es1(j+1)-es1(j))
1862  & -q2(j-1)*(es1(j)-es1(j-1)))
1863 340 CONTINUE
1864  dtdx = dt/dr(1)
1865  des(1)=dtdx*q2(1)*(es1(2)-es1(1))
1866 c
1867  DO 350 j=1,jx-1
1868  es1(j) = es1(j)+des(j)
1869 c if (kcau1+kcau2.ne.0) es1(1)=oes1(1)
1870 c if (es1(j).le.0.) es1(j)=oes1(j)
1871 350 CONTINUE
1872 
1873  do 360 j=1,jx-1
1874 c if ((y(1,1).le.0.).or.(es1(1).le.0.)) then
1875 c if ((y(j,1).le.0.).or.(es1(j).le.0.)) then
1876  if (nl(j).ge.2) then
1877  es1(j)=oes1(j)
1878  y(j,1)=ou(j,1)
1879  y(j,2)=ou(j,2)
1880  y(j,3)=ou(j,3)
1881  y(j,4)=ou(j,4)
1882  endif
1883 360 continue
1884 C
1885  kcau1=0
1886  do 553 j=1,jx-1
1887  if (y(j,1).le.0.e0) then ! if de < 0, then kcau1=1
1888  nl(j)=nl(j)+1
1889  kcau1=1
1890  nlmax=max(nlmax,nl(j))
1891  goto 553
1892  endif
1893 c es=y(j,4)/y(j,1)-((y(j,2)/y(j,1))*(y(j,2)/y(j,1))
1894 c & +(y(j,3)/y(j,1))*(y(j,3)/y(j,1)))/2.
1895 c if (es1(j).le.0.e14) then
1896  if (es1(j).le.0.e0) then
1897  nl(j)=nl(j)+1
1898  kcau1=1
1899 c print *,"!",ns
1900  nlmax=max(nlmax,nl(j))
1901  endif
1902 553 continue
1903 c------------------------------------------
1904  if (kcau1.eq.1) then
1905 c write(6,*) "kcau1=",kcau1
1906 c pause
1907 c------------------------------------------
1908  if (nlmax.le.10) then
1909  do 570 j=1,jx-1
1910  w0(j,1) = wt(j,1)
1911  w0(j,2) = wt(j,2)
1912  w0(j,3) = wt(j,3)
1913  w0(j,4) = wt(j,4)
1914  tauabso(j)=tauabst(j)
1915  tauo(j)=taut(j)
1916 570 continue
1917  time=time-dt
1918  dt=dt*0.1
1919  time=time+dt
1920 c
1921 c print *,"Check after 570 loop (nlmax < 10)"
1922 c print *,"time=",time,"dt=",dt
1923 c pause
1924 c
1925  kback=1
1926  goto 151
1927  endif
1928  if (y(j,1).le.0.) then
1929  krs(j)=1
1930  endif
1931  if (es1(j).le.0.) then
1932  kes(j)=1
1933  endif
1934  call print
1935  stop
1936 c---------------------------
1937  else
1938 c---------------------------
1939 c print *,"L3"
1940  DO 100 j=1,jx-1
1941  r = 0.5*(x1(j)+x1(j+1))
1942  r2 = r*r
1943  r3 = r2*r
1944  rm1 = r-1
1945  ok = r/rm1*sqrt(0.5/r3)
1946  rs = y(j,1)
1947  vx = y(j,2)/rs
1948  vp = y(j,3)/rs
1949  ee = y(j,4)/rs
1950  ss = rs/r
1951  es = ee-0.5*(vx*vx+vp*vp)
1952  ec2 = es1(j)*cc2
1953  hhs = sqrt(b5*(ain/ainp1)/ss)/ok
1954  beta= betao(j)
1955  te = w0(j,4)
1956 
1957 c ***** start 98 iteration! *****
1958 c@@@@@ iss > 10 for converge (best value=10) @@@@@
1959  do 98 iss=1,20
1960 c if (j.eq.1) print *,beta,te
1961  obeta = beta
1962 c
1963  a = 3*(1-beta)+beta/gm1
1964  ww = ec2*ss/(a+0.5)
1965  wc = ww/cc2
1966  hh = hhs*sqrt(ww)
1967  de = ss/(2*ain*hh*rg)
1968  demn = ain*de
1969  temn = 2*te/3.0
1970  xkff = xkkr*demn*(temn**(-3.5))
1971  tauff = xkff*de*hh*rg
1972 c TAUFF = 2.*XKFF*DE*HH*RG
1973 c TAUES = 2.*XKES*DE*HH*RG
1974  taues = xkes*de*hh*rg
1975  tau = taues+tauff
1976  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
1977  taup = qbr/(8.*aa*cc*(te**4))
1978  pr = ww/(2*ainp1*hh*rg)
1979  if (kcool.eq.1) then
1980  te = ondo(de,pr,taup,tau)
1981  else if (kcool.eq.0) then
1982  te = temp(de,pr)
1983  end if
1984  beta = rmu*de*te/pr
1985 
1986 c***** 2002-7-30 *****
1987  if (abs((beta-obeta)/beta) .lt. 1.e-6) go to 730
1988  98 continue
1989  730 continue
1990 c ***** end 98 iteration! *****
1991 c taueff = sqrt(taup*(tau+2./sqrt(3.)))
1992 c taueffjx=sqrt(tauabso(jx)*(tauo(jx)+2./sqrt(3.)))
1993 c
1994  if (kcool .eq. 1) then
1995  taueff = sqrt(taup*(tau+2./sqrt(3.)))
1996  else
1997  taueff = sqrt(3*tau*tauff)
1998  end if
1999 c
2000  if (taueff.gt.1.4618e30) then
2001  beta = betao(j)
2002 c tau = tauo(j)
2003 c taup = tauabso(j)
2004 c te = wt(j,4)
2005  rs = ou(j,1)
2006  vx = ou(j,2)/rs
2007  vp = ou(j,3)/rs
2008  ee = ou(j,4)/rs
2009  ss = rs/r
2010  ec2 = oes1(j)*cc2
2011  hhs = sqrt(b5*(ain/ainp1)/ss)/ok
2012  a = 3*(1-beta)+beta/gm1
2013  ww = ec2*ss/(a+0.5)
2014  wc = ww/cc2
2015  hh = hhs*sqrt(ww)
2016  de = ss/(2*ain*hh*rg)
2017  demn = ain*de
2018  temn = 2*te/3.0
2019  xkff = xkkr*demn*(temn**(-3.5))
2020  tauff = xkff*de*hh*rg
2021  taues = xkes*de*hh*rg
2022  tau = taues+tauff
2023  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
2024  taup = qbr/(8.*aa*cc*(te**4))
2025  pr = ww/(2*ainp1*hh*rg)
2026  te = ondo(de,pr,taup,tau)
2027  end if
2028 
2029  ws = wc/ss
2030  f(j,1) = rs*vx
2031  f(j,2) = rs*(vx*vx+ws)
2032  f(j,3) = rs*(vx*vp+alpha*ws)
2033  f(j,4) = rs*(vx*(ee+ws)+alpha*ws*vp)
2034  vk = r*ok
2035  xlk = r*vk
2036  dlnor = -(0.5+r/rm1)
2037  dlnlr = 2+dlnor
2038  dokdr = ok*dlnor/r
2039  dlkdr = xlk*dlnlr/r
2040 
2041  xk = xkes+xkff
2042  fd = 16*aa*(te**4)/(3*xk*de*hh*cc2)
2043 c--- 2002-04-11 ---
2044  if (kcool .eq.1) then
2045  fm = 8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)*unit
2046  else
2047  fm = fd
2048  end if
2049 
2050  h(j,1) = 0
2051  h(j,2) = wc*(1-dlnor)+ss*(2*vp*vk+vp*vp)
2052  h(j,3) = -ss*(vx*vp+vx*dlkdr+alpha*ws)
2053  h(j,4) = -r*rs*(vx*vp+alpha*ws)*dokdr-r*fm
2054  betaj(j)= beta
2055  betac(j)= beta
2056  okc(j) = ok
2057  hhc(j) = hh
2058  prc(j) = pr
2059  tej(j) = te
2060  tec(j) = te
2061  ssc(j) = ss
2062  dec(j) = de
2063 100 CONTINUE
2064 c ----- End of L3 loop -----
2065 
2066 c###### Thermal Conduction ####################################
2067  if (ktcond.eq.1) then
2068 c do 103 j=2,jx-1
2069  do 103 j=1,jx-1
2070 c ef=3.*xkappa*(8.-7.*betac(j))
2071 CCCc cv=3.*(8.-7.*betac(j))/(8.-3.*betac(j)-3.*betac(j))
2072  cv=3./2.
2073 cc xnu=alpha*okc(j)*(hhc(j)**2.)
2074 cc xkei=xkappa*ssc(j)*cv*xnu
2075 c gm31=gm1*(4.-3.*betac(j))
2076 c & /(betac(j)+12.*gm1*(1.-betac(j)))
2077 c ef = xkappa*gm31
2078  ef = xkappa
2079  xkei = alpha*ef*okc(j)*(hhc(j)**3)*prc(j)/tec(j)
2080  cn = cv*xkei*(tec(j+1)-tec(j))/drp(j)
2081  & *sqrt(2./xnpl1)*ainp1*unit
2082 c if (nl(j).ge.1) cn=0.
2083  f(j,4)=f(j,4)-x1(j)*cn
2084  h(j,4)=h(j,4)-cn
2085 103 continue
2086  end if
2087 c##########################################
2088 c----------------------------------
2089  end if
2090 C----------------------------------
2091 c print *,"L4"
2092  DO 120 j=2,jx-1
2093  dtdx = dt/dr(j)
2094  fsa1 = f(j,1)-f(j-1,1)
2095  fsa2 = f(j,2)-f(j-1,2)
2096  fsa3 = f(j,3)-f(j-1,3)
2097  fsa4 = f(j,4)-f(j-1,4)
2098  swa1 = 0.5*(h(j,1)+h(j-1,1))
2099  swa2 = 0.5*(h(j,2)+h(j-1,2))
2100  swa3 = 0.5*(h(j,3)+h(j-1,3))
2101  swa4 = 0.5*(h(j,4)+h(j-1,4))
2102  wk(j,1) = u(j,1)-dtdx*fsa1+dt*swa1
2103  wk(j,2) = u(j,2)-dtdx*fsa2+dt*swa2
2104  wk(j,3) = u(j,3)-dtdx*fsa3+dt*swa3
2105  wk(j,4) = u(j,4)-dtdx*fsa4+dt*swa4
2106 c endif
2107 120 CONTINUE
2108 C
2109 C****** ARTIFICIAL VISCOSITY *****
2110 C
2111 c print *,"L6"
2112 
2113  if (kcau1+kcau2.eq.0) then
2114  qb1=qb
2115  else
2116  qb1=16.
2117  endif
2118 
2119  DO 130 j = 1,jx-1
2120  vxp = w0(j+1,2)
2121  vxm = w0(j,2)
2122  q1(j)= qb1*abs(vxp-vxm)
2123 130 CONTINUE
2124 C
2125  DO 140 j=2,jx-1
2126  dtdx = dt/dr(j)
2127  d(j,1) = dtdx*(q1(j)*(u(j+1,1)-u(j,1))-q1(j-1)*(u(j,1)-u(j-1,1)))
2128  d(j,2) = dtdx*(q1(j)*(u(j+1,2)-u(j,2))-q1(j-1)*(u(j,2)-u(j-1,2)))
2129  d(j,3) = dtdx*(q1(j)*(u(j+1,3)-u(j,3))-q1(j-1)*(u(j,3)-u(j-1,3)))
2130  d(j,4) = dtdx*(q1(j)*(u(j+1,4)-u(j,4))-q1(j-1)*(u(j,4)-u(j-1,4)))
2131 140 CONTINUE
2132 C
2133  dtdx = dt/dr(1)
2134  d(1,1) = dtdx*q1(1)*(u(2,1)-u(1,1))
2135  d(1,2) = dtdx*q1(1)*(u(2,2)-u(1,2))
2136  d(1,3) = dtdx*q1(1)*(u(2,3)-u(1,3))
2137  d(1,4) = dtdx*q1(1)*(u(2,4)-u(1,4))
2138 C
2139  DO 150 j=1,jx-1
2140  u(j,1) = wk(j,1)+d(j,1)
2141  u(j,2) = wk(j,2)+d(j,2)
2142  u(j,3) = wk(j,3)+d(j,3)
2143  u(j,4) = wk(j,4)+d(j,4)
2144  es1(j) = u(j,4)/u(j,1)-((u(j,2)/u(j,1))*(u(j,2)/u(j,1))
2145  & +(u(j,3)/u(j,1))*(u(j,3)/u(j,1)))/2.
2146 150 CONTINUE
2147 
2148  do 240 j=2,jx-1
2149  dtdx = dt/dr(j)
2150  des(j) = dtdx*(q2(j)*(es1(j+1)-es1(j))
2151  & -q2(j-1)*(es1(j)-es1(j-1)))
2152 240 CONTINUE
2153  dtdx = dt/dr(1)
2154  des(1)= dtdx*q2(1)*(es1(2)-es1(1))
2155 
2156  DO 250 j=1,jx-1
2157  es1(j) = es1(j)+des(j)
2158 250 CONTINUE
2159 
2160  do 260 j=1,jx-1
2161  if (nl(j).ge.2) then
2162  es1(j) = oes1(j)
2163 c
2164  u(j,1) = ou(j,1)
2165  u(j,2) = ou(j,2)
2166  u(j,3) = ou(j,3)
2167  u(j,4) = ou(j,4)
2168  endif
2169 260 continue
2170 C
2171  kcau2=0
2172  do 153 j=1,jx-1
2173  if (u(j,1).le.0.e14) then
2174  nl(j)=nl(j)+1
2175  kcau2=1
2176  nlmax=max(nlmax,nl(j))
2177  goto 153
2178  endif
2179 c es=u(j,4)/u(j,1)-((u(j,2)/u(j,1))*(u(j,2)/u(j,1))
2180 c & +(u(j,3)/u(j,1))*(u(j,3)/u(j,1)))/2.
2181 c if (es.le.0.e14) then
2182  if (es1(j).le.0.e14) then
2183  nl(j)=nl(j)+1
2184  kcau2=1
2185 c print *,"@",ns
2186  nlmax=max(nlmax,nl(j))
2187  endif
2188 153 continue
2189 c------------------------------------------
2190  if (kcau2.eq.1) then
2191 c write(6,*) "kcau2=",kcau2
2192 c pause
2193 c------------------------------------------
2194  if (nlmax.le.10) then
2195  do 170 j=1,jx-1
2196  w0(j,1) = wt(j,1)
2197  w0(j,2) = wt(j,2)
2198  w0(j,3) = wt(j,3)
2199  w0(j,4) = wt(j,4)
2200  tauabso(j)=tauabst(j)
2201  tauo(j)=taut(j)
2202 170 continue
2203  time=time-dt
2204  dt=dt*0.1
2205  time=time+dt
2206 c
2207 c print *,"Check after 170 loop (nlmax>=10)"
2208 c print *,"time=",time,"dt=",dt
2209 c pause
2210 c
2211  kback=1
2212  goto 151
2213  endif
2214  if (u(j,1).le.0.) then
2215  krs(j)=1
2216  endif
2217  if (es1(j).le.0.) then
2218  kes(j)=1
2219  endif
2220  call print
2221  stop
2222 c---------------------------
2223  else
2224 c---------------------------
2225 C
2226 C
2227 c print *,"LW7"
2228  DO 160 j=1,jx-1
2229  r = x1(j)
2230  r2 = r*r
2231  r3 = r2*r
2232  rm1 = r-1
2233  ok = r/rm1*sqrt(0.5/r3)
2234  rs = u(j,1)
2235  vx = u(j,2)/rs
2236  vp = u(j,3)/rs
2237  ee = u(j,4)/rs
2238  ss = rs/r
2239  es = ee-0.5*(vx*vx+vp*vp)
2240  ec2=es1(j)*cc2
2241  hhs = sqrt(b5*(ain/ainp1)/ss)/ok
2242  beta=betao(j)
2243  te=tej(j)
2244 
2245 c ----- start 158 iteration! -----
2246 c@@@@@ iss > 10 for converge @@@@@
2247  do 158 iss=1,10
2248 c if (j.eq.1) print *,beta,te
2249  obeta = beta
2250 c
2251  a = 3*(1-beta)+beta/gm1
2252  ww = ec2*ss/(a+0.5)
2253  wc = ww/cc2
2254  hh = hhs*sqrt(ww)
2255  de = ss/(2*ain*hh*rg)
2256  demn = ain*de
2257  temn = 2*te/3.0
2258  xkff = xkkr*demn*(temn**(-3.5))
2259 c TAUFF = 2.*XKFF*DE*HH*RG
2260  tauff = xkff*de*hh*rg
2261 c if (tauff.le.0.) print *,"",tauff,de,hh
2262  if (tauff.le.0.) print *,"TAUABS=",tauff,"LW1D2"
2263 c TAUES = 2.*XKES*DE*HH*RG
2264  taues = xkes*de*hh*rg
2265  tau = taues+tauff
2266  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
2267  taup = qbr/(8.*aa*cc*(te**4))
2268  pr = ww/(2*ainp1*hh*rg)
2269  if (kcool.eq.1) then
2270  te=ondo(de,pr,taup,tau)
2271  else if (kcool.eq.0) then
2272  te=temp(de,pr)
2273  end if
2274 cc te=temp(de,pr)
2275  beta = rmu*de*te/pr
2276 c***** 2002-7-30 *****
2277  if (abs((beta-obeta)/beta) .lt. 1.e-6) go to 731
2278 158 continue
2279  731 continue
2280 c ***** end of 158 iteration !*****
2281 c taueff=sqrt(taup*(tau+2./sqrt(3.)))
2282 c taueffjx=sqrt(tauabso(jx)*(tauo(jx)+2./sqrt(3.)))
2283 c
2284  if (kcool .eq. 1) then
2285  taueff = sqrt(taup*(tau+2./sqrt(3.)))
2286  else
2287  taueff = sqrt(3*tau*tauff)
2288  end if
2289 c
2290  if (taueff.gt.1.4618e30) then
2291  beta=betao(j)
2292 c tau=tauo(j)
2293 c taup=tauabso(j)
2294 c te=wt(j,4)
2295  rs = ou(j,1)
2296  vx = ou(j,2)/rs
2297  vp = ou(j,3)/rs
2298  ee = ou(j,4)/rs
2299  ss = rs/r
2300  ec2 = oes1(j)*cc2
2301  a = 3*(1-beta)+beta/gm1
2302  ww = ec2*ss/(a+0.5)
2303  wc = ww/cc2
2304  hh = hhs*sqrt(ww)
2305  de = ss/(2*ain*hh*rg)
2306  demn = ain*de
2307  temn = 2*te/3.0
2308  xkff = xkkr*demn*(temn**(-3.5))
2309  tauff = xkff*de*hh*rg
2310  taues = xkes*de*hh*rg
2311  tau = taues+tauff
2312  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
2313  taup = qbr/(8.*aa*cc*(te**4))
2314  pr = ww/(2*ainp1*hh*rg)
2315  te = ondo(de,pr,taup,tau)
2316  end if
2317  ws = wc/ss
2318 c
2319  xk = xkes+xkff
2320  fd = 16*aa*(te**4)/(3.*xk*de*hh*cc2)
2321 
2322 c--- 2002-04-11 ---
2323  if (kcool .eq.1) then
2324  fm = 8.*aa*(te**4)/(1.5*tau+sqrt(3.)+1./taup)*unit
2325  else
2326  fm = fd
2327  end if
2328 c
2329  w0(j,1) = de
2330  w0(j,2) = vx
2331  w0(j,3) = vp
2332  w0(j,4) = te
2333  tauo(j) = tau
2334  tauabso(j) = taup
2335  betao(j) = beta
2336 c taueff = sqrt(taup*(tau+2./sqrt(3.)))
2337 c
2338  if (kcool .eq. 1) then
2339  taueff = sqrt(taup*(tau+2./sqrt(3.)))
2340  else
2341  taueff = sqrt(3*tau*tauff)
2342  end if
2343 c
2344  oes1(j) = es1(j)
2345  ou(j,1) = u(j,1)
2346  ou(j,2) = u(j,2)
2347  ou(j,3) = u(j,3)
2348  ou(j,4) = u(j,4)
2349 cc FFTR = FM/UNIT/XMCRIT/CC*16.
2350  fftr = fm/unit/xmcrit/cc
2351  xlumi = xlumi+2.*pi*(r*rg)*(dr(j)*rg)*fftr
2352  if (taueff.le.1.) xltn=xltn+2.*pi*(r*rg)*(dr(j)*rg)*fftr
2353  if (taueff.gt.1.) xltk=xltk+2.*pi*(r*rg)*(dr(j)*rg)*fftr
2354 c
2355 c xmdeff = 2.*pi*r2*alpha*ww/(xl-xlin)/xmcrit*rg/cc
2356  xmdeff = -2.*pi*r*rg*vx*cc*ss/xmcrit
2357 C xmdeff = 2.*pi*r2*alpha*ww/xlk/xmcrit*rg/cc
2358 c --- Data for S-shaped curve ---
2359  IF (mod(ns,klc).EQ.0) then
2360  frag1=0
2361  else if (frag1 .eq. 1) then
2362  go to 160
2363  else if (r .ge. 5.) then
2364  frag1=1
2365  go to 901
2366 c
2367  901 WRITE (10,910) time,r,ss,te,xmdeff
2368  910 FORMAT (1pe10.3,1x,1pe10.3,1x,1pe10.3,1x,1pe10.3,1x,1pe10.3,
2369  & 1x,1pe10.3)
2370  end if
2371 c ---
2372 c
2373 160 CONTINUE
2374  xlumi = xlumi+xljx
2375  xltk = xltk+xljx
2376 C
2377 c##### Data for yy.lc (Light Curve) #####
2378 c time : Time
2379 c xlumi: Total Disk luminosity
2380 c xltn : Disk luminosity for [tau < 1]
2381 c xttk : Disk luminosity for [tau > 1]
2382 c dt : Delta t at that time
2383 c NS : Time step
2384 c KLC : Light curve => file9 every KLC time steps!
2385 c########################################
2386 c
2387  if (imore .eq. 2) then
2388  if (mod(ns,200) .eq. 0) then
2389  write(6,88) "ns=",ns," t=",time," dt=",dt, " L=",xlumi
2390  88 format (a3,i12,2x,a3,1pe9.3,2x,a3,1pe9.3,2x,a3,1pe9.3)
2391  end if
2392  end if
2393 c
2394  IF (mod(ns,klc).EQ.0) WRITE (9,900) time,xlumi,xltn,xltk,dt,ns
2395  RETURN
2396 c-----------------------
2397  endif
2398 c-----------------------
2399 cc900 FORMAT (1PE16.8,1X,1PE15.7,1X,1PE15.7,1X,1PE15.7,1x,1pe15.7)
2400 900 FORMAT (1pe10.4,2x,1pe10.4,2x,1pe10.4,2x,1pe10.4,2x,1pe10.4,2x,i9)
2401  END
2402 C
2403 C
2404 C**************************
2405  SUBROUTINE cfl
2406 C**************************
2407 C
2408  include 'para.h'
2409 C
2410  dimension drp(jmax),drm(jmax),dr(jmax)
2411  dimension dtj(jmax)
2412  dimension u(jmax,4),f(jmax,4),h(jmax,4),y(jmax,4)
2413 C
2414 c print *,"CFL"
2415 C
2416  gm1 = gm-1
2417  cc2 = cc*cc
2418 
2419  DO 5 j=2,jx-1
2420  drp(j) = x1(j+1)-x1(j)
2421  drm(j) = x1(j)-x1(j-1)
2422  dr(j) = 0.5*(drp(j)+drm(j))
2423  5 CONTINUE
2424  drp(1) = x1(2) -x1(1)
2425  drm(jx)= x1(jx)-x1(jx-1)
2426  dr(1) = drp(1)
2427  dr(jx) = drm(jx)
2428  taueffjx=sqrt(tauabso(jx)*(tauo(jx)+2./sqrt(3.)))
2429 
2430  ji=1
2431  DO 10 j=ji,jx
2432  r = x1(j)
2433  rm1 = r-1
2434  r2 = r*r
2435  r3 = r2*r
2436  ok = (r/rm1)*sqrt(0.5/r3)
2437  de = w0(j,1)
2438  vx = w0(j,2)
2439  vp = w0(j,3)
2440  te = w0(j,4)
2441  if (j.eq.jx) then
2442  pr = (aa/3.)*(te**4)+rmu*de*te
2443  else
2444  if (kcool.eq.1) then
2445  fac = 1./((1.+1./tauabso(j)/(1.5*tauo(j)+sqrt(3.))))
2446  else
2447  fac = 1.0
2448  endif
2449  pr = (aa/3.)*(te**4)*fac+rmu*de*te
2450  endif
2451  beta = rmu*de*te/pr
2452  hh = sqrt(b5*pr/de)/ok
2453  ss = 2*ain*de*hh*rg
2454  ww = 2*ainp1*pr*hh*rg
2455  wc = ww/cc2
2456  ws = wc/ss
2457  a1 = 4.-3.*beta
2458  a2 = beta+12.*gm1*(1.-beta)
2459  gg1 = beta+(a1*a1*gm1)/a2
2460  gg3 = 1.+(gg1-beta)/a1
2461  cs2 = ((3.*gg1-1.)+2.*(gg3-1.)*alpha*alpha)*ws/(gg1+1.)
2462  cs = sqrt(cs2)
2463  dtj(j) = dr(j)/(cs+abs(vx))
2464  if (nl(j).ge.2) dtj(j)=dtmax
2465 c
2466  demn = ain*de
2467  temn = 2*te/3.0
2468  xkff = xkkr*demn*(temn**(-3.5))
2469  tauff = xkff*de*hh*rg
2470  taues = xkes*de*hh*rg
2471  tau = taues+tauff
2472 
2473  if (kcool .eq. 1) then
2474  taueff = sqrt(tauabso(j)*(tauo(j)+2./sqrt(3.)))
2475  else
2476  taueff = sqrt(3*tau*tauabso(j))
2477  end if
2478 
2479  if (taueff.gt.taueffjx) dtj(j)=dtmax
2480  10 CONTINUE
2481 
2482  dt = dtmax
2483  DO 20 j = ji,jx
2484 c DT = MIN(DT,.04*DTJ(J))
2485  dt = min(dt,cfc*dtj(j))
2486  20 CONTINUE
2487 c
2488  dt = max(dt,dtmin)
2489 c
2490  RETURN
2491  END
2492 C
2493 C**************************
2494  SUBROUTINE prtrb1
2495 C**************************
2496  include 'para.h'
2497 
2498  xmdot=pxmdot
2499  xmd=pxmd
2500  xlin=pxlin
2501  do 10 j=1,jx-1
2502  de = w0(j,1)
2503  vx = w0(j,2)
2504  te = w0(j,4)
2505  r = x1(j)
2506  rm1 = r-1
2507  r2 = r*r
2508  r3 = r2*r
2509  ok = (r/rm1)*sqrt(0.5/r3)
2510  pr = (aa/3.)*(te**4)+rmu*de*te
2511  hh = sqrt(b5*pr/de)/ok
2512  ss = 2*ain*de*hh*rg
2513  vx = -xmdot/(2.*pi*r*ss)/cc
2514  w0(j,2)=vx
2515  10 continue
2516  return
2517  end
2518 C**************************
2519  SUBROUTINE prtrb2
2520 C**************************
2521  include 'para.h'
2522  xmdot=pxmdot
2523  xmd=pxmd
2524 c xlin=pxlin
2525 c rc = x1(jx)
2526 c rp = 10.0
2527 
2528  do 10 j=jx,jx
2529 c CALL SADMK(X1(J),ETA,DL,W0(J,1),W0(J,2),W0(J,3),W0(J,4))
2530 c CALL SADMK2(X1(J),ETA,DL,W0(J,1),W0(J,2),W0(J,3),W0(J,4))
2531  CALL sadmn2(x1(j),eta,dl,w0(j,1),w0(j,2),w0(j,3),w0(j,4))
2532 c do 20 k=1,1
2533 c w0(j,k)=w0(j,k)
2534 c & +rp*w0(j,k)*exp(-(x1(j)-rc)*(x1(j)-rc))
2535 c 20 continue
2536 
2537  10 continue
2538  RETURN
2539  END
2540 C**************************
2541  SUBROUTINE prtrb3
2542 C**************************
2543  include 'para.h'
2544  dimension ww(4)
2545 c xmdot=pxmdot
2546 c xmd=pxmd
2547 c xlin=pxlin
2548 c CALL SADMK(X1(JX),ETA,DL,ww(1),ww(2),ww(3),ww(4))
2549 c-----------------------
2550 c rc : Radius
2551 c xlamda : wave length
2552 c rp : amplitude of the wave
2553 c-----------------------
2554 c rc=x1(JX)
2555  rc=50.0
2556 c xlamda=pxlin
2557  xlamda=1.0
2558  rp=1.0
2559  do 10 j=1,jx-1
2560  do 20 k=1,4
2561  w0(j,k)=w0(j,k)
2562  & +rp*w0(j,k)*exp(-(x1(j)-rc)*(x1(j)-rc)/xlamda/xlamda)
2563 c & -rp*w0(j,k)*exp(-(x1(j)-rc)*(x1(j)-rc)/xlamda/xlamda)
2564  20 continue
2565  10 continue
2566 c write(6,*) "PERTRB 3"
2567  return
2568  end
2569 
2570 C**************************
2571  SUBROUTINE prtrb4
2572 C**************************
2573  include 'para.h'
2574  dimension ww(4)
2575 c xmdot=pxmdot
2576 c xmd=pxmd
2577 c-----------------------
2578 c rc : Radius
2579 c xlamda : wave length
2580 c rp : amplitude of the wave
2581 c-----------------------
2582  rc=30.0
2583 c xlamda=pxlin
2584  xkk=30
2585  rp=2.0
2586 
2587  do 10 j=1,jx-1
2588  de = w0(j,1)
2589  vx = w0(j,2)
2590  te = w0(j,4)
2591  r = x1(j)
2592  rm1 = r-1
2593  r2 = r*r
2594  r3 = r2*r
2595  ok = (r/rm1)*sqrt(0.5/r3)
2596  pr = (aa/3.)*(te**4)+rmu*de*te
2597  hh = sqrt(b5*pr/de)/ok
2598  ss = 2*ain*de*hh*rg
2599  vx = -xmdot/(2.*pi*r*ss)/cc
2600  w0(j,2)=vx
2601 
2602 c do 20 k=1,1
2603 c dss = 0.01*exp(-(x1(j)-rc)**2/0.1**2)*sin(xkk*x1(j))*w0(j,k)
2604 c w0(j,k)=w0(j,k)+dss
2605  20 continue
2606  dww = dss/xkk/x1(j)*w0(j,1)
2607  dvr = dss/sqrt(xkk*x1(j))*w0(j,2)
2608  dvp = dss/sqrt(xkk*x1(j))*w0(j,3)
2609  10 continue
2610  write(6,*) "PERTRB 4"
2611  return
2612  end
2613 C**************************
2614  SUBROUTINE print
2615 C**************************
2616 C
2617  include 'para.h'
2618 C
2619  dimension drp(jmax),drm(jmax),dr(jmax)
2620 C
2621  if ((krs(j).eq.1).or.(kes(j).eq.1)) goto 100
2622  IF (mod(ns,iprint(4)).NE.0) RETURN
2623  js = iprint(2)
2624  IF (iprint(3).NE.0) THEN
2625  IF (mod(ns,iprint(3)).EQ.0) js = 1
2626  END IF
2627  100 continue
2628 C
2629 c print *,"PRINT"
2630  gm1 = gm-1
2631  cc2 = cc*cc
2632 C
2633  DO 5 j=2,jx-1
2634  drp(j) = x1(j+1)-x1(j)
2635  drm(j) = x1(j)-x1(j-1)
2636  dr(j) = 0.5*(drp(j)+drm(j))
2637  5 CONTINUE
2638  drp(1) = x1(2) -x1(1)
2639  drm(jx)= x1(jx)-x1(jx-1)
2640  dr(1) = drp(1)
2641  dr(jx) = drm(jx)
2642 
2643  WRITE(6,*)
2644  WRITE (6,1000) ns,time,dt,nlmax
2645 
2646  IF (imore.EQ.1) THEN
2647  WRITE (6,1002)
2648  ELSE IF (imore.EQ.2) THEN
2649  WRITE (6,1004)
2650  ELSE
2651  WRITE (6,1005)
2652  ENDIF
2653 
2654  DO 10 j=1,jx,js
2655  r = x1(j)
2656  rm1 = r-1
2657  r2 = r*r
2658  r3 = r2*r
2659  ok = (r/rm1)*sqrt(0.5/r3)
2660  de = w0(j,1)
2661  vx = w0(j,2)
2662  vp = w0(j,3)
2663  te = w0(j,4)
2664  if (ki.eq.1) then
2665  if (j.eq.jx) then
2666  pr = (aa/3.)*(te**4)+rmu*de*te
2667  else
2668  if (kcool .eq.1) then
2669  fac = 1./((1.+1./tauabso(j)/(1.5*tauo(j)+sqrt(3.))))
2670  else
2671  fac = 1.
2672  end if
2673  pr = (aa/3.)*(te**4)*fac+rmu*de*te
2674  endif
2675  else if (kr.eq.1) then
2676  pr = (aa/3.)*(te**4)+rmu*de*te
2677  else
2678  if (j.eq.jx) then
2679  pr = (aa/3.)*(te**4)+rmu*de*te
2680  else
2681  pr = rmu*de*te
2682  endif
2683  endif
2684 c--------------------------------------------------------
2685 c BETA : (gas pressure)/(total pressure)
2686 c HH : Disk Half Thickness
2687 c SS : Surface Density
2688 c WW : Height Integrated Pressure
2689 c WC : WW/c^2
2690 c WS : WC/SS
2691 c
2692  beta = rmu*de*te/pr
2693  hh = sqrt(b5*pr/de)/ok
2694  ss = 2*ain*de*hh*rg
2695  ww = 2*ainp1*pr*hh*rg
2696  wc = ww/cc2
2697  ws = wc/ss
2698 c------------------------------------------
2699  a1 = 4.-3.*beta
2700  a2 = beta+12.*gm1*(1.-beta)
2701  gg1 = beta+(a1*a1*gm1)/a2
2702  gg3 = 1.+(gg1-beta)/a1
2703 
2704  cs2 = ((3.*gg1-1.)+2.*(gg3-1.)*alpha*alpha)*ws/(gg1+1.)
2705  cs = sqrt(cs2)
2706  cfln = dt/dr(j)*(cs+abs(vx))
2707  vk = ok*r
2708  xlk = ok*r2
2709  xl = r*(vp+vk)
2710 c xmdeff = 2.*pi*r2*alpha*ww/(xl-xlin)/xmcrit*rg/cc
2711 c xmdeff = 2.*pi*r2*alpha*ww/xlk/xmcrit*rg/cc
2712  xmdeff = -2.*pi*r*rg*vx*cc*ss/xmcrit
2713  xmd = -2.*pi*r*ss*vx*rg*cc/xmcrit
2714 c
2715  demn = ain*de
2716  temn = 2*te/3.0
2717  xkff = xkkr*demn*(temn**(-3.5))
2718  tauff = xkff*de*hh*rg
2719  taues = xkes*de*hh*rg
2720  tau = taues+tauff
2721 
2722  qbr = 2.0*xeff*ajn*de*de*sqrt(te)*hh*rg
2723  taup = qbr/(8.*aa*cc*(te**4))
2724 c
2725  taueff = sqrt(tauabso(j)*(tauo(j)+2./sqrt(3.)))
2726 c
2727  if (kcool .eq. 1) then
2728  taueff = sqrt(tauabso(j)*(tauo(j)+2./sqrt(3.)))
2729  else
2730  taueff = sqrt(3*tau*tauff)
2731  end if
2732 c
2733  t1 = tauo(j)+2./sqrt(3.)
2734  t2 = 1./tauabso(j)
2735 c if (t1.ge.t2) then
2736 c taueff=1.
2737 c else
2738 c taueff=2.
2739 c endif
2740  a = 3*(1-beta)+beta/gm1
2741  ee = (a+0.5)*ws+0.5*(vx*vx+vp*vp)
2742 c FMC = 8.*AA*CC*(TE**4.)/(3.0*tauo(j))
2743  fmc = 8.*aa*cc*(te**4)/(1.5*tauo(j)+sqrt(3.)+1./taup)
2744 
2745  teff = (4.*fmc/aa/cc)**0.25
2746  vphi = vp+vk
2747 
2748  xll = xlk-xlin
2749  qp = 2.*alpha*alpha*xl/xll
2750 c if ((kr.eq.1).or.(r.eq.rout1)) then
2751 c QM = FF1*BETA4*XLL7*ETA5/(XKXE*R4)
2752 c else
2753  eta = vx*vx/(b2*pr/de)
2754  ff = qbr*rg/cc/cc/cc
2755  qm = 2.*pi*r*r*eta/(xmdot*ws)*ff
2756 c endif
2757  qadv = qp-qm
2758  fene = qadv/qp
2759 c write(6,*) "check!",j,r,qp,qm,fene,b2,xmdot,ws
2760 c pause
2761 c
2762 c ----- For print routine -----
2763  IF (imore.EQ.1) THEN
2764  WRITE(6,1010) j,r,dr(j),de,vx,vp,te,cs,cfln,ww,ss,xmd,hh,es1(j)
2765  * ,ws,beta,tauo(j)
2766  ELSE IF (imore.EQ.2) THEN
2767  IF (ns.EQ.lns) THEN
2768  WRITE(6,1010) j,r,beta,taueff,te,ss,xmdeff,xmd,t1,1/t2
2769  ELSE
2770  WRITE(6,1010) j,r,beta,taueff,te,ss,xmdeff,xmd,t1,1/t2
2771  ENDIF
2772  else if (imore .eq. 3) then
2773  write(6,1010) j,r,dr(j),de,te,teff,ss,vx,vp,cs,hh,tauo(j)
2774  * ,taueff,ww,xmd
2775  ELSE if (imore .eq. 0) then
2776  IF (ns.EQ.lns) THEN
2777  WRITE(6,1010) j,r,beta,teff,te,ss,-vx,vphi,cs,hh,tauo(j)
2778  * ,taueff,ww,fene,xmd
2779  ELSE
2780  WRITE(6,1010) j,r,beta,teff,te,ss,-vx,vphi,cs,hh,tauo(j)
2781  * ,taueff,ww,fene,xmd
2782  ENDIF
2783  ENDIF
2784 
2785  10 CONTINUE
2786  do 20 j=1,jx-1
2787  if (krs(j).ne.0) then
2788  write (6,1030) j,ns+1,time
2789  krs(j)=0
2790  end if
2791  20 continue
2792  do 21 j=1,jx-1
2793  if (kes(j).ne.0) then
2794  write (6,1031) j,ns+1,time
2795  kes(j)=0
2796  end if
2797 21 continue
2798 
2799  nlmax =0
2800 C
2801  RETURN
2802 C
2803  1000 FORMAT ('NS=',i12,4x,'TIME=',f12.3,4x,'DT=',e10.3,4x,
2804  & "NLMAX=",i5//)
2805  1002 format (4x,"J",5x,"R",11x,"DR",10x,"DE",10x,"VX",
2806  & 10x,"VP",10x,"TE",10x,"CS",9x,"CFLN",9x,"WW",10x,"SS",
2807  & 9x,"MDOT",9x,"HH",10x,"ES1",9x,"WS",9x,"BETA",9x,"TAU")
2808 c 1002 format (4x,"J",5x,"R",11x,"DR",10x,"DE",10x,"VX",
2809 c & 10x,"VP",10x,"TE",10x,"CS",9x,"WW",10x,"SS",
2810 c & 9x,"MDOT",9x,"HH",10x,"BETA",9x,"TAU")
2811  1004 format (4x,"J",5x,"R",10x,"BETA",8x,"TAUEFF",7x,"TE",
2812  & 10x,"SS",8x,"MDOTEFF",7x,"MDOT")
2813 c
2814 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
2815 c imore=0でのoutput data!
2816  1005 format (4x,"J",5x,"R",10x,"BETA",8x,"TEFF",8x,"TE",
2817  & 10x,"SS",10x,"VX",10x,"VP",10x,"CS",10x,"HH",
2818  & 10x,"tauo",8x,"taueff",6x,"WW",10x,"WS",10x,"MDOT")
2819 c 1005 format (4x,"J",5x,"R",10x,"BETA",9x,"TEFF",8x,"TE",
2820 c & 10x,"SS",9x,"MDOT")
2821  1010 FORMAT (i5,20(1pe12.4))
2822  1020 FORMAT (i5,5(1pe12.4),i8)
2823  1030 format ("CAUTION RS < 0 AT J =",i3," NS =",i8,
2824  & " TIME =",e10.3)
2825  1031 format ("CAUTION ES < 0 AT J =",i3," NS =",i8,
2826  & " TIME =",e10.3)
2827 C
2828  END
2829 C
2830 C
2831 C**************************
2832  SUBROUTINE file
2833 C**************************
2834 C
2835  include 'para.h'
2836 C
2837  if (kf.eq.0) then
2838  write (8,902) jx
2839  write (8,903)
2840  do 10 j=1,jx
2841  write (8,904) j,x1(j)
2842  10 continue
2843  kf=1
2844  ENDIF
2845 C
2846  IF ((mod(ns,ksfl).NE.0).AND.(ns.NE.nsx)) RETURN
2847 C IF ((ns.eq.lns).and.(ns.ne.0)) return
2848 C
2849 c WRITE (2) NS,TIME,DT,((W0(J,K),K=1,4),J=1,JX)
2850  write (8,*)
2851  write (8,905) ns,time
2852  write (8,906)
2853  do 20 j=1,jx
2854  write (8,907) j,(w0(j,k),k=1,4),tauabso(j),tauo(j),x1(j)
2855 c write (8,907) j,(w0(j,k),k=1,4),TAUABSO(J),x1(J)
2856  20 continue
2857  write (8,908) dt
2858  write (8,909) xmd
2859  write (8,901) xlin
2860 C
2861  RETURN
2862 c 900 format (80x)
2863  901 format (" XLIN=",f18.14)
2864  902 format (" JX=",i4)
2865  903 format (2x,"J",10x,"R")
2866  904 format (i4,1pe15.7)
2867  905 format (" NS=",i12," TIME=",f12.3)
2868  906 format (5x,"J",11x,"W0(J,1)",15x,"W0(J,2)",15x,"W0(J,3)",15x,
2869  & "W0(J,4)",13x,"TAUABSO",15x,"tau",15x,"R")
2870 c 907 format (i7,6(1pe22.14),I3)
2871  907 format (i7,6(1pe22.14),(1pe15.7))
2872  908 format (" THE NEXT DT=",e11.3)
2873  909 format (" XMDOT=",e11.3)
2874  END
2875 ************************************************************************
2876  subroutine rdfl7
2877 ************************************************************************
2878  include 'para.h'
2879 
2880  read (7,802) jx
2881  read (7,800)
2882  do 10 j=1,jx
2883  read (7,803) x1(j)
2884  10 continue
2885  100 continue
2886 C
2887  read (7,800)
2888  read (7,804) ns,time
2889  read (7,800)
2890  do 30 j=1,jx
2891  read (7,805) (w0(j,i),i=1,4),tauabso(j),tauo(j)
2892  if (ns.eq.0) then
2893  deini(j)=w0(j,1)
2894  vxini(j)=w0(j,2)
2895  vpini(j)=w0(j,3)
2896  teini(j)=w0(j,4)
2897  endif
2898  30 continue
2899  read (7,806) dt
2900  read (7,806) xmd
2901  read (7,801) xlin
2902  if (ns.eq.0) xmdi=xmd
2903 C
2904  call file
2905 c call print
2906  if (ns.ne.lns) goto 100
2907  xmdot=xmd*xmcrit
2908  800 format (80x)
2909  801 format (13x,f18.14)
2910  802 format (5x,i4)
2911  803 format (4x,1pe15.7)
2912  804 format (4x,i12,8x,f12.3)
2913  805 format (7x,6(1pe22.14),i3)
2914  806 format (13x,e11.3)
2915  return
2916  end
2917 ***********************************************************************
2918 ***** stifbs from Numerical Recipes (p737)
2919 ***********************************************************************
2920  subroutine stifbs (y,dydx,nv,x,htry,eps,yscal,hdid,hnext,
2921  & derivs,icon)
2922  implicit real*8 (a-h,o-z)
2923  common /lscal/ kpr
2924  integer nv,nmax,kmaxx,imax,kpr
2925 c real eps,hdid,hnext,htry,x,dydx(nv),y(nv),yscal(nv),safe1,
2926 c & safe2,redmax,redmin,tiny,scalmx
2927  dimension dydx(nv),y(nv),yscal(nv)
2928  external derivs
2929  parameter(nmax=50,kmaxx=7,imax=kmaxx+1,safe1=.25,safe2=.7,
2930  & redmax=1.e-5,redmin=.7,tiny=1.e-30,scalmx=.1)
2931  integer i,iq,k,kk,km,kmax,kopt,nvold,nseq(imax)
2932 c real eps1,epsold,errmax,fact,h,red,scale,work,wrkmin,xest,xnew,
2933 c & a(imax),alf(kmaxx,kmaxx),dfdx(nmax),dfdy(nmax,nmax),
2934 c & err(kmaxx),yerr(nmax),ysav(nmax),yseq(nmax)
2935  dimension a(imax),alf(kmaxx,kmaxx),dfdx(nmax),dfdy(nmax,nmax),
2936  & err(kmaxx),yerr(nmax),ysav(nmax),yseq(nmax)
2937  logical first,reduct
2938  save a,alf,epsold,first,kmax,kopt,nseq,nvold,xnew
2939  data first/.true./,epsold/-1./,nvold/-1/
2940 c----- Sequence is different from bsstep.
2941  data nseq /2,6,10,14,22,34,50,70/
2942 c----- Reinitialize also if nv has changed.
2943  if (eps.ne.epsold.or.nv.ne.nvold) then
2944  hnext=-1.e29
2945  xnew=-1.e29
2946  eps1=safe1*eps
2947 c----- Compute work coefficients A_k. (NR eq. 16.4.6)
2948  a(1)=nseq(1)+1
2949  do 11 k=1,kmaxx
2950  a(k+1)=a(k)+nseq(k+1)
2951  11 continue
2952 c----- Compute alpha(k,q). (NR eq. 16.4.11)
2953  do 13 iq=2,kmaxx
2954  do 12 k=1,iq-1
2955  alf(k,iq)=eps1**((a(k+1)-a(iq+1))/
2956  & ((a(iq+1)-a(1)+1.)*(2*k+1)))
2957  12 continue
2958  13 continue
2959 c-----
2960  epsold=eps
2961  nvold=nv
2962  a(1)=nv+a(1)
2963  do 14 k=1,kmaxx
2964  a(k+1)=a(k)+nseq(k+1)
2965  14 continue
2966 c----- Determine optimal row number for convergence.
2967  do 15 kopt=2,kmaxx-1
2968  if (a(kopt+1).gt.a(kopt)*alf(kopt-1,kopt)) goto 1
2969  15 continue
2970  1 kmax=kopt
2971  endif
2972  h=htry
2973 c----- Save the starting values.
2974  do 16 i=1,nv
2975  ysav(i)=y(i)
2976  16 continue
2977 c----- Evaluate Jacobian.
2978  call jac(x,y,dydx,dfdx,dfdy)
2979  if (kpr.eq.1) return
2980  if (h.ne.hnext.or.x.ne.xnew) then
2981  first=.true.
2982  kopt=kmax
2983  endif
2984  reduct=.false.
2985  2 do 18 k=1,kmax
2986  xnew=x+h
2987  if (xnew.eq.x) then
2988  icon=16000
2989  return
2990  endif
2991  call simpr (ysav,dydx,dfdx,dfdy,nmax,nv,x,h,nseq(k),yseq,
2992  & derivs)
2993  xest=(h/nseq(k))**2
2994  call pzextr (k,xest,yseq,y,yerr,nv)
2995  if (k.ne.1) then
2996  errmax=tiny
2997  do 17 i=1,nv
2998  errmax=max(errmax,abs(yerr(i)/yscal(i)))
2999  17 continue
3000  errmax=errmax/eps
3001  km=k-1
3002  err(km)=(errmax/safe1)**(1./(2*km+1))
3003  endif
3004  if (k.ne.1.and.(k.ge.kopt-1.or.first)) then
3005  if (errmax.lt.1.) goto 4
3006  if (k.eq.kmax.or.k.eq.kopt+1) then
3007  red=safe2/err(km)
3008  goto 3
3009  else if (k.eq.kopt) then
3010  if (alf(kopt-1,kopt).lt.err(km)) then
3011  red=1./err(km)
3012  goto 3
3013  endif
3014  else if (kopt.eq.kmax) then
3015  if (alf(km,kmax-1).lt.err(km)) then
3016  red=alf(km,kmax-1)*safe2/err(km)
3017  goto 3
3018  endif
3019  else if (alf(km,kopt).lt.err(km)) then
3020  red=alf(km,kopt-1)/err(km)
3021  goto 3
3022  endif
3023  endif
3024  18 continue
3025 c----- Reduce stepsize by at least REDMIN and at most REDMAX.
3026  3 red=min(red,redmin)
3027  red=max(red,redmax)
3028  h=h*red
3029  reduct=.true.
3030 c----- Try again!
3031  goto 2
3032 c----- Successful step taken.
3033  4 x=xnew
3034  hdid=h
3035  first=.false.
3036 c----- Compute optimalrow for convergence and corresponding stepsize.
3037  wrkmin=1.e35
3038  do 19 kk=1,km
3039  fact=max(err(kk),scalmx)
3040  work=fact*a(kk+1)
3041  if (work.lt.wrkmin) then
3042  scale=fact
3043  wrkmin=work
3044  kopt=kk+1
3045  endif
3046  19 continue
3047  hnext=h/scale
3048 c --- Check for possible order increase, but not if stepsize was just reduced.
3049  if (kopt.ge.k.and.kopt.ne.kmax.and..not.reduct) then
3050  fact=max(scale/alf(kopt-1,kopt),scalmx)
3051  if (a(kopt+1)*fact.le.wrkmin) then
3052  hnext=h/fact
3053  kopt=kopt+1
3054  endif
3055  endif
3056  icon=10
3057  return
3058  end
3059 *******************************************************************
3060 ***** simpr (Semi-implicit Extrapolation Method by NR )
3061 *******************************************************************
3062  subroutine simpr (y,dydx,dfdx,dfdy,nmax,n,xs,htot,nstep,yout,
3063  & derivs)
3064  implicit real*8 (a-h,o-z)
3065  common /lscal/ kpr
3066  integer n,nmax,nstep,nmaxx,kpr
3067 c real htot,xs,dfdx(n),dfdy(nmax,nmax),dydx(n),y(n),yout(n)
3068  dimension dfdx(n),dfdy(nmax,nmax),dydx(n),y(n),yout(n)
3069  external derivs
3070 c --- Maximum expected value of n.
3071  parameter(nmaxx=50)
3072  integer i,j,nn,indx(nmaxx)
3073 c real d,h,x,a(nmaxx,nmaxx),del(nmaxx),ytemp(nmaxx)
3074  dimension a(nmaxx,nmaxx),del(nmaxx),ytemp(nmaxx)
3075 c --- Stepsize
3076  h=htot/nstep
3077 c --- Set up the matrix 1-hf'.
3078  do 12 i=1,n
3079  do 11 j=1,n
3080  a(i,j)=-h*dfdy(i,j)
3081  11 continue
3082  a(i,i)=a(i,i)+1.
3083  12 continue
3084 c --- LU decommposition of the matrix.
3085  call ludcmp (a,n,nmaxx,indx,d)
3086 c --- Set up right hand side for first step.
3087 c --- Use yout for temporary storage.
3088  do 13 i=1,n
3089  yout(i)=h*(dydx(i)+h*dfdx(i))
3090  13 continue
3091  call lubksb (a,n,nmaxx,indx,yout)
3092 c --- First step.
3093  do 14 i=1,n
3094  del(i)=yout(i)
3095  ytemp(i)=y(i)+del(i)
3096  14 continue
3097  x=xs+h
3098 c --- Use yout for temporary storage of derivatives.
3099  call derivs(x,ytemp,yout)
3100  if (kpr.eq.1) return
3101 c --- General step.
3102  do 17 nn=2,nstep
3103 c --- set up right hand side for general step.
3104  do 15 i=1,n
3105  yout(i)=h*yout(i)-del(i)
3106  15 continue
3107  call lubksb (a,n,nmaxx,indx,yout)
3108  do 16 i=1,n
3109  del(i)=del(i)+2.*yout(i)
3110  ytemp(i)=ytemp(i)+del(i)
3111  16 continue
3112  x=x+h
3113  call derivs (x,ytemp,yout)
3114  if (kpr.eq.1) return
3115  17 continue
3116 c --- Set up right hand side for last step.
3117  do 18 i=1,n
3118  yout(i)=h*yout(i)-del(i)
3119  18 continue
3120  call lubksb(a,n,nmaxx,indx,yout)
3121 c --- take last step.
3122  do 19 i=1,n
3123  yout(i)=ytemp(i)+yout(i)
3124  19 continue
3125  kpr=0
3126  return
3127  end
3128 ***************************************************************
3129 ***** pzextr (The polynomial extrapolation routine by NR)
3130 ***************************************************************
3131  subroutine pzextr (iest,xest,yest,yz,dy,nv)
3132  implicit real*8 (a-h,o-z)
3133  common /lscal/ kpr
3134  integer iest,nv,imax,nmax,kpr
3135 c real xest,dy(nv),yest(nv),yz(nv)
3136  dimension dy(nv),yest(nv),yz(nv)
3137  parameter(imax=13,nmax=50)
3138  integer j,k1
3139 c real delta,f1,f2,q,d(nmax),qcol(nmax,imax),x(imax)
3140  dimension d(nmax),qcol(nmax,imax),x(imax)
3141  save qcol,x
3142 c --- Save current independent variable.
3143  x(iest)=xest
3144  do 11 j=1,nv
3145  dy(j)=yest(j)
3146  yz(j)=yest(j)
3147  11 continue
3148 c --- Store first estimate in first column.
3149  if (iest.eq.1) then
3150  do 12 j=1,nv
3151  qcol(j,1)=yest(j)
3152  12 continue
3153  else
3154  do 13 j=1,nv
3155  d(j)=yest(j)
3156  13 continue
3157  do 15 k1=1,iest-1
3158  delta=1./(x(iest-k1)-xest)
3159  f1=xest*delta
3160  f2=x(iest-k1)*delta
3161 c --- propagate tableau 1 diagonal more.
3162  do 14 j=1,nv
3163  q=qcol(j,k1)
3164  qcol(j,k1)=dy(j)
3165  delta=d(j)-q
3166  dy(j)=f1*delta
3167  d(j)=f2*delta
3168  yz(j)=yz(j)+dy(j)
3169  14 continue
3170  15 continue
3171  do 16 j=1,nv
3172  qcol(j,iest)=dy(j)
3173  16 continue
3174  endif
3175  return
3176  end
3177 ********************************************************************
3178 ***** ludcmp
3179 ********************************************************************
3180  subroutine ludcmp (a,n,np,indx,d)
3181  implicit real*8 (a-h,o-z)
3182  common /lscal/ kpr
3183  integer n,np,indx(n),nmax,kpr
3184 c real d,a(np,np),tiny
3185  dimension a(np,np)
3186  parameter(nmax=500,tiny=1.0e-20)
3187  integer i,imax,j,k
3188 c real aamax,dum,sum,vv(nmax)
3189  dimension vv(nmax)
3190  d=1.
3191  do 12 i=1,n
3192  aamax=0.
3193  do 11 j=1,n
3194  if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
3195  11 continue
3196 c if (aamax.eq.0.) pause "singular matrix in ludcmp"
3197  if (aamax.eq.0.) stop "singular matrix in ludcmp"
3198  vv(i)=1./aamax
3199  12 continue
3200  do 19 j=1,n
3201  do 14 i=1,j-1
3202  sum=a(i,j)
3203  do 13 k=1,i-1
3204  sum=sum-a(i,k)*a(k,j)
3205  13 continue
3206  a(i,j)=sum
3207  14 continue
3208  aamax=0.
3209  do 16 i=j,n
3210  sum=a(i,j)
3211  do 15 k=1,j-1
3212  sum=sum-a(i,k)*a(k,j)
3213  15 continue
3214  a(i,j)=sum
3215  dum=vv(i)*abs(sum)
3216  if (dum.ge.aamax) then
3217  imax=i
3218  aamax=dum
3219  endif
3220  16 continue
3221  if (j.ne.imax) then
3222  do 17 k=1,n
3223  dum=a(imax,k)
3224  a(imax,k)=a(j,k)
3225  a(j,k)=dum
3226  17 continue
3227  d=-d
3228  vv(imax)=vv(j)
3229  endif
3230  indx(j)=imax
3231  if (a(j,j).eq.0.) a(j,j)=tiny
3232  if (j.ne.n) then
3233  dum=1./a(j,j)
3234  do 18 i=j+1,n
3235  a(i,j)=a(i,j)*dum
3236  18 continue
3237  endif
3238  19 continue
3239  return
3240  end
3241 ***********************************************************************
3242 ***** lubksb
3243 ***********************************************************************
3244  subroutine lubksb (a,n,np,indx,b)
3245  implicit real*8 (a-h,o-z)
3246  common /lscal/ kpr
3247  integer n,np,indx(n),kpr
3248 c real a(np,np),b(n)
3249  dimension a(np,np),b(n)
3250  integer i,ii,j,ll
3251  ii=0
3252  do 12 i=1,n
3253  ll=indx(i)
3254  sum=b(ll)
3255  b(ll)=b(i)
3256  if (ii.ne.0) then
3257  do 11 j=ii,i-1
3258  sum=sum-a(i,j)*b(j)
3259  11 continue
3260  else if (sum.ne.0.) then
3261  ii=i
3262  endif
3263  b(i)=sum
3264  12 continue
3265  do 14 i=n,1,-1
3266  sum=b(i)
3267  do 13 j=i+1,n
3268  sum=sum-a(i,j)*b(j)
3269  13 continue
3270  b(i)=sum/a(i,i)
3271  14 continue
3272  return
3273  end
kes
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic xltk common lscal kpr common shoki kes(jmax) common/coula/nlmax
ki
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 ki
Definition: para.h:11
xkappa
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp xkappa
Definition: para.h:11
intmdl
subroutine intmdl
Definition: accdisk.f:199
pxmdot
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb pxmdot
Definition: para.h:11
pxlin
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb pxlin
Definition: para.h:11
intgrt
subroutine intgrt(R, ETA, DL, IC)
Definition: accdisk.f:803
rdfl7
subroutine rdfl7
Definition: accdisk.f:2877
lw1d
subroutine lw1d
Definition: accdisk.f:1615
pxmd
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb pxmd
Definition: para.h:11
file8
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 file8
Definition: para.h:11
file9
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic xltk common lscal kpr common shoki file9
Definition: para.h:34
temp
double precision function temp(DE, PR)
Definition: accdisk.f:1347
parameter
IMPLICIT REAL *A O Z parameter(jmax=500) common/BCLCK/ITIME common/BCNST/PI
hini
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error hini
Definition: para.h:11
tauo
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic tauo(jmax) common/lumi/xlumi
relerr
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error relerr
Definition: para.h:11
file
subroutine file
Definition: accdisk.f:2833
deini
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic xltk common lscal kpr common shoki deini(jmax)
cfl
subroutine cfl
Definition: accdisk.f:2406
sadmk2
subroutine sadmk2(R, ETA, DL, DE, VX, VP, TE)
Definition: accdisk.f:467
sadmn
subroutine sadmn(R, ETA, DL)
Definition: accdisk.f:626
fun
subroutine fun(XI, YI, YP)
Definition: accdisk.f:984
jac
subroutine jac(xi, yi, yp, dfdr, pdj)
Definition: accdisk.f:1151
lubksb
subroutine lubksb(a, n, np, indx, b)
Definition: accdisk.f:3245
xltn
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic xltn
Definition: para.h:22
print
subroutine print
Definition: accdisk.f:2615
prtrb3
subroutine prtrb3
Definition: accdisk.f:2542
vxini
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic xltk common lscal kpr common shoki vxini(jmax)
sadmk
subroutine sadmk(R, ETA, DL)
Definition: accdisk.f:266
rdprm
subroutine rdprm
Definition: accdisk.f:91
fnjn
double precision function fnjn(NPL)
Definition: accdisk.f:1478
kcool
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich kcool
Definition: para.h:11
file7
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 file7
Definition: para.h:11
simpr
subroutine simpr(y, dydx, dfdx, dfdy, nmax, n, xs, htot, nstep, yout, derivs)
Definition: accdisk.f:3064
fnin
double precision function fnin(NPL)
Definition: accdisk.f:1462
ondo
double precision function ondo(DE, PR, taup, tau)
Definition: accdisk.f:1420
vpini
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic xltk common lscal kpr common shoki vpini(jmax)
teini
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic xltk common lscal kpr common shoki teini(jmax) common/cauti/krs(jmax)
tauabso
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic tauabso(jmax)
pzextr
subroutine pzextr(iest, xest, yest, yz, dy, nv)
Definition: accdisk.f:3132
xmd
IMPLICIT REAL *A O Z RG B5 common BPARM xmd
Definition: para.h:4
at
program at
Definition: accdisk.f:25
prtrb2
subroutine prtrb2
Definition: accdisk.f:2520
sadmn2
subroutine sadmn2(R, ETA, DL, DE, VX, VP, TE)
Definition: accdisk.f:725
cnvrbs
subroutine cnvrbs
Definition: accdisk.f:1495
nl
IMPLICIT REAL *A O Z RG B5 common BPARM XNPL1 common BMESH JX common BSTEP lns common BWRT1 IRSTRT common BPHY0 kf common licur file10 common curve IMORE common swich ktcond c common error abserr common bradi mesh common artvp QU common radia kr common prtrb xmdi common optic xltk common lscal kpr common shoki nl(jmax) common/energ/es1(jmax) common/fluxr0/drepsi(jmax)
stifbs
subroutine stifbs(y, dydx, nv, x, htry, eps, yscal, hdid, hnext, derivs, icon)
Definition: accdisk.f:2922
prtrb4
subroutine prtrb4
Definition: accdisk.f:2572
prtrb1
subroutine prtrb1
Definition: accdisk.f:2495
ludcmp
subroutine ludcmp(a, n, np, indx, d)
Definition: accdisk.f:3181