C.hr gchirv
C@
*....*...1.........2.........3.........4.........5.........6.........7.*
      real*4 function gchirv(a,iseed)
      implicit none
      save

*     algorithm to generate random variables with the chi distribution
*     with a degrees of freedom,  for a greater than or equal to one

*     j f monahan (may, 1986) dept of statistics, ncsu, raleigh, nc usa
*      modified 1/24/90 by a r gallant to set seed differently

*     If you call with x=gchirv(alpha,iseed) then x*x/2 is gamma with 
*     shape alpha/2 and unit scale.  

*     John F. Monahan, An Algorithm for Generating Chi Random Variables, 
*     ACM TOMS vol. 13 (1987) pp. 168-172.  

      real*4 a
      integer*4 iseed

      real*4 alpha,alphm1,beta,u,v,vmin,vdif,z,zz,rnum,w,s,emhlf,vmaxu
      real*4 rsqrt2,emhlf4,eqtrt2,c
      real*4 ran

      intrinsic SQRT, ALOG

      data emhlf/ .6065307 /, vmaxu/ .8577639 /, rsqrt2/ .7071068 /
      data emhlf4/ .1516327 /, eqtrt2/ 2.568051 /, c/ 1.036961 /

      data alpha/ 0. /
*                     is this alpha the same as the last one?
      if( a .eq. alpha ) go to 1
*                     do a little setup
      alpha = a
      alphm1 = alpha - 1.
      beta = SQRT( alphm1 )
*                     get dimensions of box
      vmaxu = emhlf * ( rsqrt2 + beta )/( .5 + beta )
      vmin =  -beta
      if( beta .gt. 0.483643 ) vmin = emhlf4/alpha - emhlf
      vdif = vmaxu - vmin
*                     start ( and restart ) algorithm here
  1   continue
      u = ran(iseed)
      v = vdif*ran(iseed) + vmin
      z = v / u
*                     do some quick reject checks first
      if( z .le. - beta ) go to 1
      zz = z * z
      rnum = 2.5 - zz
      if( z .lt. 0. ) rnum = rnum + zz * z / ( 3. * (z + beta ) )
*                     do quick inner check
      if( u .lt. rnum/eqtrt2 ) go to 9
      if( zz .ge. c / u + 1.4 ) go to 1
*                     above was knuth's normal outer check
      w = 2. * ALOG( u )
*                     now the real check
      s = - ( zz / 2. + z * beta )
      if( beta .gt. 0. ) s = alphm1*ALOG(1.+z/beta) + s
      if( w .gt. s ) go to 1
*                     success here -- transform and go
  9   gchirv = z + beta
      return
      end