00001
00002 c Generate a King model.
00003 c
00004 c Parameters: -n number
00005 c -w depth
00006 c -i random seed
00007 c
00008 c Note: Uses functions qromb, rkode, and random numbers
00009 c from "Numerical Recipes"
00010
00011 program kingmodel
00012
00013 implicit real*8 (a-h, o-z)
00014 parameter (NMAX = 50000)
00015
00016 real*8 mass(NMAX), pos(NMAX, 3), vel(NMAX, 3)
00017 real*8 w0
00018 real*8 r, random
00019 integer n, iargc
00020 integer iseed
00021 character*80 arg
00022
00023 c Settable parameters: w0 = dimensionless depth
00024 c n = number of particles
00025
00026 iout = 1
00027
00028 c Initialize random number generator using iseed.
00029 iseed = 42
00030
00031 n = 0
00032 do i=1,iargc()
00033 call getarg(i, arg)
00034 if (arg(1:2) .eq. '-w') then
00035 call getarg(i+1, arg)
00036 read(arg, '(i10)', iostat=io) iw
00037 if (io.eq.0) then
00038 w0 = iw
00039 else
00040 read(arg, '(f8.3)', iostat=io) w0
00041 if (io.ne.0) then
00042 write(6,'(a)')'Error reading w0.'
00043 stop
00044 end if
00045 end if
00046 call setking(w0)
00047 else if (arg(1:2) .eq. '-n') then
00048 call getarg(i+1, arg)
00049 read(arg, '(i10)', iostat=io) n
00050 if (io.ne.0) then
00051 write(6,'(a)')'Error reading N.'
00052 stop
00053 end if
00054 else if (arg(1:2) .eq. '-i') then
00055 call getarg(i+1, arg)
00056 read(arg, '(i10)', iostat=io) iseed
00057 if (io.ne.0) then
00058 write(6,'(a)')'Error reading iseed.'
00059 stop
00060 end if
00061 r = random(iseed)
00062 end if
00063 end do
00064
00065 if (n .le. 0) then
00066 write(6,'(a)')'Number of stars (N) not set: use "-n #".'
00067 stop
00068 end if
00069
00070 call king(mass, pos, vel, n, NMAX, iout)
00071
00072 do i=1,n
00073 print *, i, mass(i), (pos(i,k),k=1,3), (vel(i,k),k=1,3)
00074 end do
00075
00076 end
00077
00078
00079 block data init_king
00080
00081 implicit real*8 (a-h,o-z)
00082
00083 common/kpars1/w0,v0,e0,pi,twopi
00084 common/kpars2/nprof
00085 common/therm/g(0:1000)
00086
00087 data w0/0./
00088 data g(0)/-1./
00089
00090 end
00091
00092
00093 subroutine setking(w)
00094
00095 implicit real*8 (a-h,o-z)
00096
00097 common/kpars1/w0,v0,e0,pi,twopi
00098 common/kpars2/nprof
00099
00100 w0=w
00101
00102 end
00103
00104
00105 subroutine king(body,x,xdot,n,nmx,iout)
00106
00107 c Set up a King model, with total mass = 1, core radius = 1..
00108
00109 implicit real*8 (a-h,o-z)
00110
00111 dimension body(nmx),xdot(nmx,3)
00112 real*8 x(nmx,3)
00113
00114 parameter (nm=2500)
00115 common/profile/rr(0:nm),d(0:nm),v2(0:nm),psi(0:nm),
00116 & zm(0:nm),indx(0:100)
00117 common/kpars1/w0,v0,e0,pi,twopi
00118 common/kpars2/nprof
00119
00120 if (w0.eq.0.) then
00121 write(6,'(a)')'Dimensionless depth of'
00122 & ' cluster potential (w0) not set: use "-w #".'
00123 stop
00124 end if
00125
00126 c Compute the cluster density/velocity/potential profile.
00127
00128 call poisson(rr,d,v2,psi,zm,nm,w0,nprof,v20,iout)
00129
00130 c do i=0,nprof
00131 c write(0,*), i, rr(i), d(i)
00132 c end do
00133
00134 c Index the mass distribution.
00135
00136 do 10 i=0,nprof
00137 zm(i)=zm(i)/zm(nprof)
00138 10 continue
00139 indx(0)=0
00140 indx(100)=nprof
00141
00142 c Determine core and half-mass radii.
00143
00144 dz=.01
00145 z=dz
00146 iz=1
00147 do 20 j=1,nprof-1
00148 if (rr(j).lt.1.)jcore=j
00149 if (zm(j).lt..5)jhalf=j
00150 if (zm(j).gt.z) then
00151 indx(iz)=j-1
00152 z=z+dz
00153 iz=iz+1
00154 end if
00155 20 continue
00156 rhalf=rr(jhalf)+(rr(jhalf+1)-rr(jhalf))*(.5-zm(jhalf))
00157 & /(zm(jhalf+1)-zm(jhalf))
00158 zmcore=zm(jcore)+(zm(jcore+1)-zm(jcore))*(1.-rr(jcore))
00159 & /(rr(jcore+1)-rr(jcore))
00160
00161 if (iout.ne.0)write(6,30)w0,nprof,v20,rr(nprof),rhalf,zmcore
00162 30 format(/' King model: W0 =',f6.2,', NPROF = ',i4/
00163 & ' V20/E0 =',f7.3,', Rt/Rc =',f7.2,
00164 & ', Rh/Rc =',f6.2,', Mc/M =',f7.4/)
00165
00166 zmass=0.
00167 body(1)=1.
00168
00169 do 50 i=1,n
00170 body(i)=body(1)
00171 zmass=zmass+body(i)
00172 50 continue
00173 d0=zmass/zm(nprof)
00174 pi=4.*atan(1.)
00175 twopi=2.*pi
00176 four3pi=4.*pi/3.
00177 e0=four3pi*d0/v20
00178 v0=.004*sqrt(2.*e0)
00179
00180 c Assign positions and velocities. Note that it may actually
00181 c be preferable to do this in layers instead.
00182
00183 c Note also that rcore is taken to be unity.
00184
00185 do 100 i=1,n
00186 call posvel(xx,yy,zz,uu,vv,ww)
00187 x(i,1)=xx
00188 x(i,2)=yy
00189 x(i,3)=zz
00190 xdot(i,1)=uu
00191 xdot(i,2)=vv
00192 xdot(i,3)=ww
00193 100 continue
00194
00195 end
00196
00197
00198 subroutine posvel(x,y,z,vx,vy,vz)
00199
00200 implicit real*8 (a-h,o-z)
00201
00202 parameter (nm=2500)
00203 common/profile/xx(0:nm),d(0:nm),v2(0:nm),psi(0:nm),
00204 & zm(0:nm),indx(0:100)
00205
00206 common/kpars1/w0,v0,e0,pi,twopi
00207 common/kpars2/nprof
00208 common/therm/g(0:1000)
00209
00210 dimension v33(0:1000)
00211 real*8 random,r
00212 data v33(1000)/0./
00213
00214 if (v33(1000).eq.0.) then
00215 do 5 i=0,1000
00216 v33(i)=(.004*i)**3/3.
00217 5 continue
00218 end if
00219
00220 c Choose radius randomly from the mass distribution.
00221
00222 r=random(0)
00223 i=100*r
00224 do 10 i1=indx(i),indx(i+1)+1
00225 if (zm(i1).gt.r) go to 15
00226 10 continue
00227
00228 write(6,*)' **** ERROR in king.'
00229 stop
00230
00231 15 rfac=(r-zm(i1-1))/(zm(i1)-zm(i1-1))
00232 r=xx(i1-1)+rfac*(xx(i1)-xx(i1-1))
00233
00234 c Angular position random.
00235
00236 cth=2.*random(0)-1.
00237 sth=sqrt(1.-cth**2)
00238 ph=twopi*random(0)
00239 cph=cos(ph)
00240 sph=sin(ph)
00241 x=r*sth*cph
00242 y=r*sth*sph
00243 z=r*cth
00244
00245 c Choose speed randomly from the distribution at this radius.
00246
00247 p=-(psi(i1-1)+rfac*(psi(i1)-psi(i1-1)))
00248 pfac=exp(p)
00249 il=0
00250 rl=0.
00251 iu=250.*sqrt(p)
00252 ru=pfac*g(iu)-v33(iu)
00253 r=random(0)*ru
00254 100 if (iu-il.gt.1) then
00255 im=(il+iu)/2
00256 rm=pfac*g(im)-v33(im)
00257 if (rm.gt.r) then
00258 iu=im
00259 ru=rm
00260 else
00261 il=im
00262 rl=rm
00263 end if
00264 go to 100
00265 end if
00266 v=v0*(il+(r-rl)/(ru-rl))
00267
00268 c Direction is random.
00269
00270 cth=2.*random(0)-1.
00271 sth=sqrt(1.-cth**2)
00272 ph=twopi*random(0)
00273 cph=cos(ph)
00274 sph=sin(ph)
00275 vx=v*sth*cph
00276 vy=v*sth*sph
00277 vz=v*cth
00278
00279 end
00280
00281
00282 subroutine poisson(x,d,v2,psi,zm,nmax,w0,nprof,v20,iout)
00283
00284 c Self-contained 1-D (spherical) Poisson's equation solver.
00285 c Currently knows only about King models.
00286 c
00287 c Input: nmax is the maximum number of points allowed
00288 c w0 is the dimensionless central potential
00289 c iout allows messages if nonzero
00290 c
00291 c Output: x is scaled radius
00292 c d is scaled density
00293 c v2 is scaled velocity dispersion
00294 c psi is potential
00295 c zm is cumulative mass (scaling from x, d scalings)
00296 c nprof is the actual number of points generated
00297 c v20 is the central velocity dispersion
00298
00299 implicit real*8 (a-h,o-z)
00300
00301 external rhs
00302 parameter (RLIN=.25, NLIN=100, RMAX=1.e4, TOL=1.e-6)
00303
00304 dimension x(0:nmax),d(0:nmax),v2(0:nmax),psi(0:nmax),
00305 & zm(0:nmax),y(2)
00306
00307 common/center/di
00308
00309 parameter (four3=1.333333333333333, fourpi=12.5663706)
00310
00311 psi0=-abs(w0)
00312
00313 c Initialize at center of cluster.
00314
00315 xn=0.
00316 y(1)=0.
00317 y(2)=psi0
00318 x(0)=0.
00319 psi(0)=psi0
00320 v2(0)=1.
00321 call densvel(psi0,di,v2(0))
00322 d(0)=di
00323 di=1./di
00324 zm(0)=0.
00325 fac=10.**(log10(RMAX/RLIN)/(NMAX-NLIN))
00326
00327 c Equation to be solved is
00328 c
00329 c (x psi)'' = 9 x d
00330 c
00331 c by defining y(1) = (x psi)
00332 c y(2) = y(1)'
00333 c
00334 c Cover the first RLIN core radii linearly with NLIN points;
00335 c the remaining coverage is logarithmic, out to RMAX core radii,
00336 c if necessary. Stop when d <= 0.
00337
00338 do 1000 i=1,nmax
00339 xo=xn
00340 if (i.le.NLIN) then
00341 xn=(RLIN*i)/NLIN
00342 else
00343 xn=fac*xo
00344 end if
00345 call rkode(y,2,xo,xn,TOL,.1*(xn-xo),.0001*(xn-xo),
00346 & nok,nbad,rhs,ier)
00347 c
00348 c N.B. Remember that y(1) is x*psi
00349 c
00350 if (ier.ne.0) then
00351 if (iout.ne.0)write(6,100)ier
00352 100 format(1x,'ERROR #',i3,' IN RKODE')
00353 stop
00354 end if
00355 x(i)=xn
00356 psi(i)=y(1)/xn
00357 v2(i)=1.
00358 call densvel(psi(i),d(i),v2(i))
00359 if (d(i).lt.0.) then
00360 x(i)=x(i-1)+(x(i)-x(i-1))/(1.-d(i)/d(i-1))
00361 d(i)=0.
00362 v2(i)=0.
00363 end if
00364 zm(i)=x(i)*y(2)-y(1)
00365 if (d(i).le.0.) go to 2000
00366 1000 continue
00367 i=nmax
00368 2000 nprof=i
00369 v20=v2(0)
00370 do 2500 i=nprof,0,-1
00371 d(i)=d(i)/d(0)
00372 v2(i)=v2(i)/v2(0)
00373 zm(i)=fourpi/9.*zm(i)
00374 2500 continue
00375
00376 99999 return
00377 end
00378
00379 c
00380 subroutine rhs(x,y,ypr)
00381
00382 c Define RHS of ODE, for use by rkode.
00383
00384 implicit real*8 (a-h,o-z)
00385
00386 dimension y(2),ypr(2)
00387 common/center/di
00388 data zero/0./
00389
00390 ypr(1)=y(2)
00391
00392 if (x.le.0.) then
00393 d=1.
00394 else
00395 call densvel(y(1)/x,d,zero)
00396 d=d*di
00397 end if
00398 ypr(2)=9.*x*d
00399
00400 end
00401
00402
00403 subroutine densvel(psi,d,v2)
00404
00405 c Return scaled density d and velocity dispersion v2,
00406 c given scaled potential psi.
00407
00408 implicit real*8 (a-h,o-z)
00409
00410 d=0.
00411 if (psi.ge.0.) return
00412
00413 dw=sqrt(-psi)
00414 w=dw
00415 g4=v2
00416 call gg(w,g2,g4)
00417 ep=exp(-psi)
00418 g2e=g2*ep
00419 d=g2e+dw*psi/3.
00420 if (v2.gt.0.and.d.gt.0.) v2=2.*(g4*ep-.2*dw*psi**2)/d
00421
00422 end
00423
00424
00425 subroutine gg(w,g2,g4)
00426
00427 implicit real*8 (a-h,o-z)
00428
00429 common/therm/g(0:1000)
00430
00431 c Array g contains the integral (0-->x) of y**2 exp(-y**2) dy
00432
00433 external gaus2
00434 data eps/1.e-6/
00435
00436 if (g(0).lt.0.) then
00437 g(0)=0.
00438 xo=0.
00439 dx=.004
00440 do 100 i=1,1000
00441 x=xo+dx
00442 g(i)=g(i-1)+qromb(gaus2,xo,x,eps,ntr)
00443 xo=x
00444 100 continue
00445 end if
00446
00447 ww=250.*w
00448 iw=ww
00449 if (iw.ge.1000) then
00450 g2=g(1000)
00451 else
00452 g2=g(iw)+(ww-iw)*(g(iw+1)-g(iw))
00453 end if
00454 if (g4.lt.0.) return
00455 w2=w*w
00456 ew=exp(-w2)
00457 wew=w*ew
00458 g4=1.5*g2-.5*w2*w*ew
00459
00460 end
00461
00462
00463 function gaus2(x)
00464
00465 implicit real*8 (a-h,o-z)
00466
00467 gaus2=x*x*exp(-x*x)
00468
00469 end
00470 c
00471 c=======================================================================
00472
00473 c Routines from "Numerical Recipes".
00474
00475 c qromb - Romberg function integrator
00476 c rkode - Runge-Kutta ODE integrator
00477
00478 function qromb(f,a,b,eps,ntrap)
00479
00480 c Integrate the function f from a to b using Romberg integration
00481 c with (relative) error tolerance eps.
00482 c |ntrap| on return is number of trapezoidal subdivisions.
00483 c ntrap < 0 indicates an error.
00484
00485 c See "Numerical Recipes," Chapter 4.3
00486
00487
00488 implicit real*8 (a-h,o-z)
00489
00490 parameter (one=1.,quarter=.25,zero=0.)
00491 cd parameter (one=1.d0,quarter=.25d0,zero=0.d0)
00492 parameter (jmax=20,jmaxp=jmax+1,k=5,km=k-1)
00493 dimension s(jmaxp),h(jmaxp)
00494 external f
00495
00496 ntrap=0
00497 h(1)=one
00498 do 10 j=1,jmax
00499 call trapzd(f,a,b,s(j),j)
00500 ntrap=j
00501 if (j.ge.k) then
00502 call polint(h(j-km),s(j-km),k,zero,qromb,dss,ier)
00503 if (ier.ne.0) then
00504 ntrap=-ntrap
00505 qromb=0.
00506 return
00507 end if
00508 if (abs(dss).lt.eps*abs(qromb)) return
00509 end if
00510 s(j+1)=s(j)
00511 h(j+1)=quarter*h(j)
00512 10 continue
00513 ntrap=-1-ntrap
00514
00515 end
00516
00517
00518 subroutine trapzd(f,a,b,s,n)
00519
00520 implicit real*8 (a-h,o-z)
00521
00522 external f
00523 parameter (half=.5)
00524 save it
00525
00526 if (n.eq.1) then
00527 s=half*(b-a)*(f(a)+f(b))
00528 it=1
00529 else
00530 tnm=it
00531 del=(b-a)/tnm
00532 x=a+half*del
00533 sum=0.
00534 do 10 j=1,it
00535 sum=sum+f(x)
00536 x=x+del
00537 10 continue
00538 s=half*(s+del*sum)
00539 it=2*it
00540 end if
00541
00542 end
00543
00544
00545 subroutine polint(xa,ya,n,x,y,dy,ier)
00546
00547 implicit real*8 (a-h,o-z)
00548
00549 parameter (nmax=10)
00550 dimension xa(n),ya(n),c(nmax),d(nmax)
00551
00552 ns=1
00553 dif=abs(x-xa(1))
00554 do 10 i=1,n
00555 dift=abs(x-xa(i))
00556 if (dift.lt.dif) then
00557 ns=i
00558 dif=dift
00559 end if
00560 c(i)=ya(i)
00561 d(i)=ya(i)
00562 10 continue
00563 y=ya(ns)
00564 ns=ns-1
00565 do 20 m=1,n-1
00566 do 15 i=1,n-m
00567 ho=xa(i)-x
00568 hp=xa(i+m)-x
00569 w=c(i+1)-d(i)
00570 den=ho-hp
00571 if (den.eq.0.) then
00572 ier=1
00573 return
00574 end if
00575 den=w/den
00576 d(i)=hp*den
00577 c(i)=ho*den
00578 15 continue
00579 if (2*ns.lt.n-m) then
00580 dy=c(ns+1)
00581 else
00582 dy=d(ns)
00583 ns=ns-1
00584 end if
00585 y=y+dy
00586 20 continue
00587 ier=0
00588
00589 end
00590
00591 c
00592 c-----------------------------------------------------------------------
00593
00594 subroutine rkode(ystart,nvar,x1,x2,eps,h1,hmin,nok,nbad,
00595 & derivs,ier)
00596
00597 c FIFTH-order Runge-Kutta integrator with adaptive step size.
00598 c See "Numerical Recipes," Chapter 15.2
00599 c -- Note that their common block /path/ is not used.
00600
00601 implicit real*8 (a-h,o-z)
00602
00603 parameter (maxstp=10000,nmax=10,two=1.,zero=0.,tiny=1.e-30)
00604 dimension ystart(nvar),yscal(nmax),y(nmax),dydx(nmax)
00605 external derivs
00606
00607 ier=0
00608 x=x1
00609 h=sign(h1,x2-x1)
00610 nok=0
00611 nbad=0
00612 kount=0
00613
00614 do 10 i=1,nvar
00615 y(i)=ystart(i)
00616 10 continue
00617
00618 do 50 nstp=1,maxstp
00619 call derivs(x,y,dydx)
00620 do 20 i=1,nvar
00621 yscal(i)=abs(y(i))+abs(h*dydx(i))+tiny
00622 20 continue
00623 if ((x+h-x2)*(x+h-x1).gt.zero)h=x2-x
00624 if (x+h.eq.x) go to 39
00625 call rk(y,dydx,nvar,x,h,eps,yscal,hdid,hnext,derivs,ier)
00626 if (ier.ne.0) go to 9999
00627 if (hdid.eq.h) then
00628 nok=nok+1
00629 else
00630 nbad=nbad+1
00631 end if
00632 if ((x-x2)*(x2-x1).lt.zero) go to 45
00633 39 do 40 i=1,nvar
00634 ystart(i)=y(i)
00635 40 continue
00636 return
00637 45 if (abs(hnext).lt.hmin) go to 99999
00638 h=hnext
00639 50 continue
00640
00641 9999 ier=ier+1
00642 99999 ier=ier+1
00643
00644 end
00645
00646
00647 subroutine rk(y,dydx,n,x,htry,eps,yscal,
00648 & hdid,hnext,derivs,ier)
00649
00650 implicit real*8 (a-h,o-z)
00651
00652 parameter (nmax=10,pgrow=-.2,pshrink=-.25,fcor=1./15.,
00653 & one=1.,safety=.9,errcon=6.e-4)
00654 external derivs
00655 dimension y(n),dydx(n),yscal(n),ytemp(nmax),ysav(nmax),
00656 & dysav(nmax)
00657
00658 ier=0
00659 xsav=x
00660 do 10 i=1,n
00661 ysav(i)=y(i)
00662 dysav(i)=dydx(i)
00663 10 continue
00664 h=htry
00665 1 hh=.5*h
00666 call rk4(ysav,dysav,n,xsav,hh,ytemp,derivs)
00667 x=xsav+hh
00668 call derivs(x,ytemp,dydx)
00669 call rk4(ytemp,dydx,n,x,hh,y,derivs)
00670 x=xsav+h
00671 if (x.eq.xsav) go to 99999
00672 call rk4(ysav,dysav,n,xsav,h,ytemp,derivs)
00673 errmax=0.
00674 do 20 i=1,n
00675 ytemp(i)=y(i)-ytemp(i)
00676 errmax=max(errmax,abs(ytemp(i)/yscal(i)))
00677 20 continue
00678 errmax=errmax/eps
00679 if (errmax.gt.one) then
00680 h=safety*h*(errmax**pshrink)
00681 go to 1
00682 else
00683 hdid=h
00684 if (errmax.gt.errcon) then
00685 hnext=safety*h*(errmax**pgrow)
00686 else
00687 hnext=4.*h
00688 end if
00689 end if
00690 do 30 i=1,n
00691 y(i)=y(i)+ytemp(i)*fcor
00692 30 continue
00693 return
00694
00695 99999 ier=1
00696
00697 end
00698
00699
00700 subroutine rk4(y,dydx,n,x,h,yout,derivs)
00701
00702 implicit real*8 (a-h,o-z)
00703
00704 parameter (nmax=10)
00705 dimension y(n),dydx(n),yout(n),yt(nmax),dyt(nmax),dym(nmax)
00706 external derivs
00707
00708 hh=.5*h
00709 h6=h/6.
00710 xh=x+hh
00711 do 10 i=1,n
00712 yt(i)=y(i)+hh*dydx(i)
00713 10 continue
00714 call derivs(xh,yt,dyt)
00715 do 20 i=1,n
00716 yt(i)=y(i)+hh*dyt(i)
00717 20 continue
00718 call derivs(xh,yt,dym)
00719 do 30 i=1,n
00720 yt(i)=y(i)+h*dym(i)
00721 dym(i)=dyt(i)+dym(i)
00722 30 continue
00723 call derivs(x+h,yt,dyt)
00724 do 40 i=1,n
00725 yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.*dym(i))
00726 40 continue
00727
00728 end
00729
00730 c
00731 c=======================================================================
00732
00733 c Don't need to define random or anything below this line if a
00734 c built-in exists...
00735
00736 function random(iseed)
00737
00738 c Return a random number in the range [0,1), using the congruential
00739 c random number generator CONGRNO below.
00740
00741 c Convention: iseed = 0 ==> return next number in the sqruence.
00742 c iseed > 0 ==> initialize the random generator.
00743
00744 implicit real*8 (a-h, o-z)
00745 real*8 random
00746 real*8 congrno
00747
00748 save irset
00749 data irset/0/
00750
00751 if (irset.eq.0.or.iseed.gt.0) then
00752
00753 c No system built-in random number generator is known:
00754
00755 dum = congrno(iseed)
00756
00757 irset = 1
00758 end if
00759
00760 c Return the next random number in the sequence.
00761
00762 random = congrno(0)
00763
00764 end
00765
00766
00767 function congrno(iseed)
00768
00769 c Random number generator (Press et al. p. 195).
00770 c ---------------------------------------------
00771
00772 implicit real*8 (a-h,o-z)
00773 real*8 congrno
00774
00775 parameter (M=714025, IA=1366, IC=150889, RM=1.d0/M)
00776
00777 integer ir(97)
00778 save iy,iff,ir,idum
00779 data iff /0/
00780
00781 if (iseed.lt.0.or.iff.eq.0) then
00782 idum = -abs(iseed)
00783 iff = 1
00784 idum = mod(IC-idum,M)
00785 do 11 j = 1,97
00786 idum = mod(IA*idum+IC,M)
00787 ir(j) = idum
00788 11 continue
00789 idum = mod(IA*idum+IC,M)
00790 iy = idum
00791 end if
00792
00793 j = 1 + (97*iy)/m
00794 if (j.gt.97.or.j.lt.1) write (6,12) j, idum
00795 12 format (/' RAN2: j, idum = ',2i12)
00796
00797 iy = ir(j)
00798 congrno = iy*RM
00799 idum = mod(IA*idum+IC,M)
00800 ir(j) = idum
00801
00802 end
