! genosim.f90:         Simulates genotypes.raw for genomic predictions
!
!                      Fortran program written by Paul VanRaden  2006
!                      Animal Improvement Programs Lab, USDA
!                      Paul.VanRaden@ars.usda.gov
!
! OPTION DEFINITIONS:         (reads from genosim.options file)
!                Morgans   =   length of each chromosome
!                LDparam   = 0 generates founding alleles in equilibrium
!                          = x from > 0 to < 1 generates more linkage
!                genout    = 0 does not output true genotype file
!                          = 1 outputs true genotypes for animals with chip > -1
!                hapout    = 0 does not output haplotype file
!                          = 1 outputs haplotypes for only genotyped animals
!                          = 2 outputs haplotypes for all animals
!                          = 3 outputs haplotypes for SNPs and QTLs of all
!                chromes   =   number of chromosome pairs
!                Xchrome   =   X chromosome number (usually last, or 0 if none)
!                notread   =   fraction of missing genotypes generated
!                errate    =   fraction of incorrect genotypes generated
!                mutate    =   fraction of mutated haplotypes / generation
  USE IFPORT
  implicit none
  real*8, dimension(:,:),  allocatable :: BV
  real*8, dimension(:),    allocatable :: qtlvar
  real,   dimension(:,:),  allocatable :: effect
  real,   dimension(:),    allocatable :: freq,fmiss,ageSNP
  real,   dimension(:),    allocatable :: corr,MAF
  integer*1,dimension(:,:,:),allocatable :: inherit
  integer,dimension(:,:),  allocatable :: ped,chro,cmark,inheritT
  integer,dimension(:,:),  allocatable :: chiplist
  integer,dimension(:),    allocatable :: loclist,chipcount,chipsize
  integer,dimension(:),    allocatable :: locate,location
  integer,dimension(:),    allocatable :: eqnlist,animkey,ids,eqns
  integer,dimension(:),    allocatable :: marker,biryr,ancestor
  integer,dimension(:),    allocatable :: listeqn,chromnum,nloci,ploci
  integer,dimension(:),    allocatable :: genotype,idvec,chip,chipvec
  integer*1,dimension(:),  allocatable :: directn
  character*1,dimension(:),allocatable :: sexv,tmp
  character*1 sex
  character*8 birth
  character*2 chrnum
  character*50 markername
  integer*8    nmiss
  integer maxped,maxgen,loci,qtls,maxqtl,maxped2
  integer seednum,chromes,chrome,Xchrome,idnum,maxn,seed1(1)
  integer nhash,pnteqn,nanc
  integer traits,hapout,genout/1/
  integer chips,chipn,chip0
  integer prevmin,maxloci/0/,maxloc,marks
  integer which,mloci/0/,loc,position,n20,n26,trt,bigtrt/1/
  integer cross(0:99999),maxcross/0/,missing/5/
  integer nped,total,totcross,mutations/0/
  integer prev,next,ncross,len,recomb,ngen
  integer i,j,nn,pnt,pnt1,pnt2,l,locus,anim
  integer storage
  parameter (storage=10000000)
  real    bigqtl,thisqtl,notread,errate,mutate
  real    random_normal,corr1/1./
  real    var,uniform,corruni,ageuni,Morgans,LDparam
  real*8  clength
! -----------------------------------------------  Read option file
  print *,'Output from genosim.f90'

  open ( 9,file='genosim.options')
  read(9,*)              Morgans,LDparam,hapout
  print *,           '   Morgans,LDparam,hapout'
  print '(f8.3,f8.5,i7)',Morgans,LDparam,hapout
  print *
  read(9,*)
  read(9,*)
  read(9,*)           chromes,Xchrome,notread,errate,mutate
  print *,       '    chromes,Xchrome,notread,errate,mutate'
  print '(2i7,2f10.5,f15.10)',chromes,Xchrome,notread,errate,mutate
!
! -----------------------------------------------  Input data files
! sort pedfile by ascending birth year !
                    open (10,file='pedigree.file')
                    open (14,file='chromosome.data')
                    open (15,file='true.effects')
                    open (25,file='simseed')
                    open (26,file='genotype.data0')
! ...................................................  Output files
   if (hapout > 0)  open (17,file='haplotype.data')
                    open (16,file='genomic.bv')
                    open (20,file='genotypes.raw')
                    open (21,file='genotypes.true')
!                  (unit  30+ later opened separately by chromosome)
! ------------------------------------------------ Declare matrices
  read (25,*)
  read (25,*) seednum, loci, qtls, traits
  print *
  print *,traits,' traits will be generated (from markersim.options)'
  maxloc = 3 * loci / chromes
!                                     Get maxped from pedigree.file
  maxped = 0
  do
    read(10,*,end=1) sex,idnum,pnt1,pnt2,birth
    maxped = maxped + 1
    end do
 1 rewind (10)
  maxped2 = 2 * maxped
  allocate (sexv(maxped),chip(0:maxped),chipvec(0:maxped))
  allocate (locate(chromes*storage),location(loci),eqns(0:maxped2))
  allocate (ped(2,maxped),ids(maxped2),biryr(maxped))
  allocate (eqnlist(maxped),animkey(maxped))
  allocate (freq(loci),fmiss(loci),marker(loci),ageSNP(loci))
  allocate (listeqn(0:maxped),corr(loci),MAF(loci))
  allocate (directn(0:loci),ancestor(maxped))
  allocate (genotype(loci),tmp(loci),idvec(maxped))
  allocate (chromnum(0:loci),chro(maxloc,chromes))
  allocate (cmark(maxloc,chromes))
  allocate (nloci(chromes),ploci(chromes))
  allocate (BV(traits,maxped),effect(traits,loci))
  allocate (qtlvar(traits))
!                                                          Set seed
  seed1(1) = seednum
  if (seednum > 0) call random_seed (put = seed1)
!                            rand function is for Absoft compiler
! if (seednum > 0) uniform = rand(seednum)
! -------------------------------------------   Read chromosome.data (twice)
!                                                 Get loci and chips
  read (14,*,end=2)
  do
    read(14,*,end=2) markername,chrome,loc,mloci,position,chips
!   if (mloci < 10) print *,markername,chrome,loc,mloci,position,chips
    end do
 2 continue
  rewind 14
  allocate (chiplist(mloci,-1:chips),chipsize(0:chips))
  allocate (loclist(0:chips),chipcount(0:chips))
!                                                    Store marker list
  chiplist = 0
  read (14,*,end=10)
  do
    read(14,*,end=10) markername,chrome,loc,mloci,position,chips &
      ,(chiplist(mloci,chipn),chipn=1,chips)
!   if (mloci < 10) print '(a,5i12)', markername,chrome,loc,mloci,position,chips
    if (chrome > chromes) then
      print *,chrome,' exceeds # chromosomes ',chromes
      stop
      end if
    end do
 10 print *
  print *,mloci,' loci on ',chromes,' chromosomes' &
         ,' read from chromosome.data'
  chipsize = 0
  do loc = 1, mloci
    do chipn = 1, chips
      chipsize(chipn) = max(chiplist(loc,chipn),chipsize(chipn))
      end do
    end do
! ----------------------------- Read allele frequencies and effects
  mloci = loci
  if (qtls < loci) mloci = loci - qtls
  do locus = 1,loci
    read(15,*,end=995) freq(locus),position &
      ,(effect(trt,locus),trt=1,traits)
    end do ! locus
  print *
  print *,loci,' allele effects read from true.effects'
  print *
  locate = 0
! ------------------------------  Generate marker and QTL locations
  do locus = 1,loci
!                                       Random: find empty location
    do i = 1,1000
      call random_number (uniform)
      next = drand(0)*chromes*storage + 1
!                                                 Absoft statements
!     uniform = rand(0)
!     next = uniform*chromes*storage + 1
!
      if (locate(next) == 0) then
        locate(next) = locus
        go to 3
        end if
      end do ! i
    print *,'QTL location table too full!'
    stop
 3  continue
    end do ! locus
!                                            Condense QTL locations
  next = 0
  do i = 1,chromes*storage
    if (locate(i) > 0) then
      next = next + 1
      location(next) = i
      end if
    end do ! i
  if (next .ne. loci) print *,'loci mismatch ',next,loci
  qtlvar = 0.d0
  bigqtl = 0.
  print *,'  marker, chrome, position,  freq   (first 10 loci)'
  prevmin = 0
  ploci = 0
  chromnum(0) = 0
!                                  Renumber loci within chromosomes
  do locus = 1,loci
    chromnum(locus) = (location(locus) - 1)/storage + 1
    if (chromnum(locus) .ne. chromnum(locus - 1)) loc = 0
    loc = loc + 1
    nloci(chromnum(locus)) = loc
    maxloci = max(maxloci,loc)
    if (maxloci > maxloc) then
      print *,maxloc,' maxloc too small, increase !'
      stop
      end if
    chro(loc,chromnum(locus)) = locus
!                                         Genetic variance from QTL
    if (locate(location(locus)) <= qtls) then
      var = 2. * freq(locus) * (1. - freq(locus))
!     effect(locus) = random_normal() * sqrt(1./(qtls * var))
      do trt = 1,traits
        thisqtl = effect(trt,locus)**2 * var
        qtlvar(trt) = qtlvar(trt) + thisqtl
        if (thisqtl > bigqtl) then
          bigqtl = thisqtl
          bigtrt = trt
          maxqtl = locus
          end if
        end do ! trt
      end if
    if (locus <= 10 .or. locus == loci) &
      print '(3i9,f9.5,2x,99f10.5)',locus,chromnum(locus) &
            ,location(locus),freq(locus)
    end do ! locus
  print *
  if (qtls < loci) then
    print *,qtls,' of the ',loci,' loci had effects' 
    print *,'     QTL,  chrome,  position,   freq,' &
           ,'     effects on traits'
    next = 0
    do locus = 1,loci
      if (locate(location(locus)) <= qtls) then
        next = next + 1
        if (next <= 10) print '(3i9,f9.5,2x,99f10.5)',locus &
          ,chromnum(locus),location(locus),freq(locus) &
          ,(effect(trt,locus),trt=1,traits)
      else
        marker(locus - next) = locus
!                                           Markers on each chrome
        ploci(chromnum(locus)) = ploci(chromnum(locus)) + 1
        marks = ploci(chromnum(locus))
        cmark(marks,chromnum(locus)) = locus - next
        end if
      end do ! locus
  else
    print *,qtls,'QTL and marker locations assumed to be identical'
    do locus = 1,loci
      marker(locus) = locus
      end do ! locus
    end if
  print *
  print *,'Largest QTL var ',bigqtl,' was '
  print *,100*bigqtl/qtlvar(bigtrt) &
         ,'% of genetic variance = ',qtlvar(bigtrt) &
         ,'  for trait ',bigtrt
  print *
  do i = 1,chromes
    print *,'chromosome ',i,' contains ',nloci(i),' loci'
    end do
  print *
  print *,'locus   freq  minorfreq   corr(n-1)  corr(1)    ageSNP'
!                                          Set linkage correlations
  do locus = 1,loci
    corr(locus) = LDparam
!                                       Random direction of linkage
    call random_number (uniform)
    if (uniform > .5) corr(locus) = -1. * corr(locus)
    MAF(locus) = .5 - abs(freq(locus) - .5)
!                                                   Age of mutation
    if (locus > 1) then
      if (abs(freq(locus) - freq(locus-1)) > .02) then
        call random_number (ageuni)
        ageSNP(locus) = ageuni
      else
        ageSNP(locus) = ageSNP(locus-1)
        end if
      end if
    if (locus <= 10) print '(i5,5f10.4)' &
       ,locus,freq(locus),MAF(locus),corr(locus),corr1,ageSNP(locus)
    corr1 = corr1 * corr(locus)
    end do ! locus
  eqns(0) = 0
  nn = nhash(0,-1,ids,maxped2,maxn)
! ----------------------------------------------  Read pedigree data
  eqnlist = 0
  listeqn = 0
  nped = 0
  do
    read(10,*,end=5) sex,idnum,pnt1,pnt2,birth
    nped = nped + 1
    read(birth,'(i4)') biryr(nped)
    sexv(nped) = sex
    ped(1,nped) = pnt1
    ped(2,nped) = pnt2
    if(pnt1 <= 0) ped(1,nped) = 0
    if(pnt2 <= 0) ped(2,nped) = 0
    idvec(nped) = idnum
!                                                    Hash individual
    nn = nhash(idnum,1,ids,maxped2,maxn)
    if (maxn .ne. nped) then
      print *,'Duplicate IDs in pedigree.file !'
      stop
      end if
    eqns(nn) = nped
    end do
 5 print *
   print *,nped,' records read from pedigree file'
  ancestor = 0
  nanc = 0
!                                                       Hash parents
  do i = 1,nped
    do pnt = 1,2
      pnteqn = eqns(nhash(ped(pnt,i),0,ids,maxped2,maxn))
      if (pnteqn >= i) then
        print *,'progeny ',idvec(i),' must precede parent ' &
         ,ped(pnt,i),' in pedigree.file'
        ped(pnt,i) = 0
      else
        ped(pnt,i) = pnteqn
        end if
!                                               Fill missing parents
      if (ped(pnt,i) == 0) then
        do j = i, 1, -1
          if ((sexv(j) == 'M' .and. pnt == 1) .or. &
              (sexv(j) == 'F' .and. pnt == 2)) then
            if (ancestor(j) == 0 .and. biryr(j) <= biryr(i)-5) then
              ped(pnt,i) = j
              ancestor(j) = 1
              nanc = nanc + 1
              exit
              end if
            end if
          end do ! j
        end if
      end do ! pnt
    end do ! i
  print *,nanc,' missing parents set to known 5 years earlier'
!
! -----------------------  Assign chip densities from genotype.data0
  chip  = 1
  chipvec  = 1
  chip(0) = -1
  chips = 1
  n26 = 0
  ngen = 0
!                                                     2 header lines
  read (26,*)
  read (26,*)
  do
    n26 = n26 + 1
    read (26,*,end=6) idnum, chip0
!                                           Reset default chip to -1
    if (n26 == 1) chip = -1
    if (n26 == 1) chipvec = -1
    chips = max(chip0,chips)
    chipcount(chip0) = chipcount(chip0) + 1
    nn = nhash(idnum,0,ids,maxped2,maxn)
    if (nn == 0) then
      print *,idnum,' ID in genotype.data0 ' &
             ,'not found in pedigree.file for n26 ',n26
      cycle
      end if
    if (chip0 > -1) then
      ngen = ngen + 1
      eqnlist(ngen) = eqns(nn)
      animkey(ngen) = idnum
      listeqn(eqns(nn)) = ngen
      chip(listeqn(eqns(nn))) = chip0
      end if
    chipvec(eqns(nn)) = chip0
    end do ! read 26
 6  n26 = n26 - 1
  print *,n26,' records read from genotype.data0'
  print *
  print *,'  animals  chip  markers'
  do i = 0,chips
    if (chipcount(i) > 0) print '(i9,i6,i9)', &
      chipcount(i),i,chipsize(i)
    end do
  print *
  n20 = 0
! ----------------------------------------------- Simulate genotypes
  clength = 1. / Morgans
  print *,'Simulate genotypes'
  allocate (inherit(maxloci,2,maxped))
  allocate (inheritT(loci,2))
  BV = 0.
  totcross = 0
  total = 0
  cross(0) = 0
!                               Process by chromosome to save memory
  do chrome = 1,chromes
!                    Reset seed to get same genotypes with added ped
!   seed1(1) = seednum + chrome
!   call random_seed (put = seed1)
!                                                    For each animal
  do i = 1,nped
!                                         Maternal, paternal gametes
    do pnt = 1,2
!                                         Marker and QTL inheritance
      if (ped(pnt,i) > 0) then
        ncross = 0
        recomb = 0
!                                         Chromosome segment lengths
        call random_number (uniform)
        if (uniform == 0.) uniform = 1.
        len = storage * log(1./uniform) * clength + 1
        next = chrome*storage
        prev = (chrome-1)*storage
!                                      Crossovers on each chromosome
        do
          if (prev+len < next) then
            ncross = ncross + 1
            if (ncross > 99999) then
              print *,'Crossover limit reached ' &
                     ,ncross,prev,len,next
              stop
              end if
            cross(ncross) = prev+len
            prev = prev+len
            call random_number (uniform)
            if (uniform == 0.) uniform = 1.
            len = storage * log(1./uniform) * clength + 1
          else
            exit
            end if
          end do ! chomosome segments
!       print *,'Number of crossovers, recombinations = ' &
!              ,ncross - recomb,recomb
        totcross = totcross + ncross - recomb
        total = total + 1
        maxcross = max(maxcross,ncross - recomb)
        next = 1
        which = 1
        corr1 = 1.
        call random_number (uniform)
        if (uniform < .5) which = 2
!                                         Alleles from known parents
        do loc=1,nloci(chrome)
          locus = chro(loc,chrome)
!                                            Are crossovers between?
 30       if (next > ncross) go to 35
          if (cross(next) <= location(locus)) then
!                                           Toggle which from 1 to 2
            which = 3 - which
            next = next + 1
            if (cross(next) > location(locus)) go to 35
            go to 30
            end if
 35       continue
!                                             Transfer parent allele
          inherit(loc,pnt,i) = inherit(loc,which,ped(pnt,i))
          call random_number (uniform)
!                                                Introduce mutations
          if (uniform < mutate) then
            inherit(loc,pnt,i) = 3 - inherit(loc,pnt,i)
            mutations = mutations + 1
            end if
          end do ! loc
      else
!                                             Unknown parent alleles
        call random_number (uniform)
        call random_number (corruni)
!                                                 Starting direction
        if (uniform < .5) then
          directn(0) = 1
        else
          directn(0) = 2
          end if
        do loc=1,nloci(chrome)
          locus = chro(loc,chrome)
!                                      Underlying (.5 freq) alleles
          call random_number (uniform)
!                                             Blocks of linked loci
          if (uniform < abs(corr(locus))) then
            if (corr(locus) > 0.) then
              directn(loc) = directn(loc - 1)
            else
              directn(loc) = 3 - directn(loc - 1)
              end if
          else
!                                       Break points between blocks
            call random_number (uniform)
            if (uniform < .5) then
              directn(loc) = 1
            else
              directn(loc) = 2
              end if
!                                  Correlated alleles within blocks
            call random_number (corruni)
            end if
!         inherit(loc,pnt,i) = directn(loc)
!                                       Account for age of mutation
!                                            for SNPs with high MAF
          if (uniform > (ageSNP(locus) * 2.*MAF(locus))**2) then
            inherit(loc,pnt,i) = directn(loc)
          else
            inherit(loc,pnt,i) = 3 - directn(loc)
            end if
!                                              Set allele frequency
          if (corruni > 2.*MAF(locus)) then
            if (freq(locus) < .5) then
              inherit(loc,pnt,i) = 1
            else
              inherit(loc,pnt,i) = 2
              end if
            end if
          end do ! loc
!       print '(i5,2x,99999i1)',pnt,(directn(loc),loc=1,80)
        end if
      end do ! pnt
!                                                      X chromosome
    if (chrome == Xchrome .and. sexv(i) == 'M') then
      do loc=1,nloci(chrome)
        inherit(loc,1,i) = inherit(loc,2,i)
        end do ! loc
      end if
!                                                    Breeding values
    do pnt = 1,2
      do loc=1,nloci(chrome)
        locus = chro(loc,chrome)
!                                                    Add QTL effects !!! (it subtracts the effects instead)
        do trt = 1,traits
          BV(trt,i) = BV(trt,i) + (1. - inherit(loc,pnt,i) &
             + freq(locus))*effect(trt,locus)
          end do ! trt
        end do ! loc
      end do ! pnt
!                                        Remove QTLs from haplotypes
    do loc = 1,nloci(chrome)
      locus = chro(loc,chrome)
      do pnt = 1,2
        inheritT(locus,pnt) = inherit(loc,pnt,i)
        end do
      end do
!                                Write all haplotypes (not including QTLs)
!
    nn = nhash(idvec(i),0,ids,maxped2,maxn)
    if ((hapout == 1 .and. chipvec(i) > -1) .or. &
         hapout == 2) write(17,'(3i10,1x,2000000i1)') idvec(i) &
      ,cmark(1,chrome),cmark(ploci(chrome),chrome) &
      ,(inheritT(marker(loc),1),inheritT(marker(loc),2) &
      ,loc=cmark(1,chrome),cmark(ploci(chrome),chrome))
!                                    Write all haplotypes (including QTLs)
!
    if (hapout == 3) write(17,'(3i10,1x,2000000i1)') idvec(i) &
      ,chro(1,chrome),chro(nloci(chrome),chrome) &
      ,(inherit(loc,1,i),inherit(loc,2,i) &
      ,loc=1,nloci(chrome))
    end do ! i
! ---------------------------------------------------- Finished chromosome
  write (chrnum,'(i2)') chrome
  if (chrnum(1:1) == ' ') chrnum(1:1) = '0'
  open (30 + chrome,file='genotypes.raw.'//chrnum)
  do i = 1,ngen
!                                              Store genotypes across loci
!                                                        Exclude ancestors
    do loc = 1,nloci(chrome)
      locus = chro(loc,chrome)
      do pnt = 1,2
        inheritT(locus,pnt) = inherit(loc,pnt,eqnlist(i))
        end do ! pnt
      end do ! loc
!                                        Convert haplotypes to genotype
    loc = 0
    do locus = cmark(1,chrome),cmark(ploci(chrome),chrome)
      loc = loc + 1
      genotype(loc) = inheritT(marker(locus),1) &
                    + inheritT(marker(locus),2) - 2
      end do ! locus
    write(30 + chrome,'(3i10,1x,100000000i1)') animkey(i) &
      ,chip(i),loc,(genotype(locus),locus=1,loc)
    end do ! i
  print *,'Finished chromosome ',chrome
  end do ! chrome
! -----------------------------------------------------------  Write BV
  do i = 1,nped
    write(16,'(i9,i4,99f10.3)') idvec(i) &
      ,traits,(bv(trt,i),trt=1,traits)
    end do ! i
  nmiss = 0
!                                  Reread separate chromosome genotypes
  do chrome = 1, chromes
    rewind (30 + chrome)
    end do ! chrome
  do i = 1,ngen
    l = 0
    do chrome = 1, chromes
      read(30 + chrome,'(3i10,1x,100000000i1)') anim &
        ,chip0,loc,(genotype(locus + l),locus=1,loc)
      l = l + loc
      end do
!                                                  Write true genotypes
    if (genout == 1) &
      write(21,'(3i10,1x,100000000i1)') animkey(i) &
        ,chip(i),mloci,(genotype(locus),locus=1,mloci)
!                                           Introduce missing genotypes
    do loc = 1, mloci
      call random_number (uniform)
      if (uniform < notread) genotype(loc) = missing
!                                           Introduce genotyping errors
      call random_number (uniform)
      if (uniform < errate) then
!                                            Homozygote to heterozygote
        if (genotype(loc) == 0 .or. genotype(loc) == 2) then
          genotype(loc) = 1
!                                                        Display errors
!         genotype(loc) = 6
        else
!                                            Heterozygote to homozygote
          if (uniform >  errate/2.) genotype(loc) = 0
          if (uniform <= errate/2.) genotype(loc) = 2
!                                                        Display errors
!         if (uniform >  errate/2.) genotype(loc) = 7
!         if (uniform <= errate/2.) genotype(loc) = 8
          end if
        end if
      end do ! loc
!                                                         Lower density chips
    chipn = chip(i)
    loc = 0
    do locus = 1,mloci
      if (chiplist(locus,chipn) > 0) then
        genotype(chiplist(locus,chipn)) = genotype(locus)
        if (genotype(locus) == missing) nmiss = nmiss + 1
        loc = loc + 1
      else
        if (chipn > 0) nmiss = nmiss + 1
        end if
      end do ! locus
! ------------------------------------------------------  Write genotype data
    if (loc > 0) then
      write(20,'(3i10,1x,100000000i1)') animkey(i) &
        ,chip(i),loc,(genotype(locus),locus=1,loc)
      n20 = n20 + 1
      end if
    end do ! i
  print *,n20,' genotype records written to unit 20'
!                                                    Delete chromosome files
  do chrome = 1, chromes
    close (30 + chrome, dispose='delete')
    end do ! chrome
  print *
  print '(f9.4,a)',100.*nmiss / n20 / mloci &
                  ,'% of genotypes were missing'
  print *
  print *,'Average crossovers per chromosome = ' &
         ,totcross/float(total)
  print *,'Maximum crossovers per chromosome = ',maxcross
  print *
  print *,'Total number of mutations introduced = ',mutations
  print *
  print *,'Program genosim completed normally'
    stop
 990 print *,'phenotype.data0 file is missing'
    stop
 995 print *,'Unexpected end of file, unit 15 (true.effects)'
    stop
 996 print *,'Unexpected end of file, unit 16 (genomic.bv)'
    stop
 999 print *,'Unexpected end of file, unit 10 (pedigree.file)'
  end
! -----------------------------------------------------------------------
! Function from NAG
! This file contains:
! 1. The module RANDOM for generating random deviates from a number of
!    continuous and discrete distributions.   The module also includes a
!    simple quicksort routine.
! 2. A program to demonstrate its use.   The program includes code for the
!    log-gamma function.
! The module uses the intrinsic RANDOM_NUMBER function supplied with the
! compiler.   It has only been tested using the Lahey LF90 compiler.

! The two parts are separated by a line of exclamation marks.

! Alan Miller

! The compilers own random number generator, SUBROUTINE RANDOM_NUMBER(r),
! is used to provide a source of uniformly distributed random numbers.

REAL FUNCTION random_normal()

! Adapted from the following Fortran 77 code
!      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
!      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
!      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.

!  The function random_normal() returns a normally distributed pseudo-random
!  number with zero mean and unit variance.

!  The algorithm uses the ratio of uniforms method of A.J. Kinderman
!  and J.F. Monahan augmented with quadratic bounding curves.

IMPLICIT NONE
REAL                      :: zero = 0.0, half = 0.5, one = 1.0, two = 2.0

!     Local variables
REAL     :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472,           &
            r1 = 0.27597, r2 = 0.27846, u, v, x, y, q, rand

!     Generate P = (u,v) uniform in rectangle enclosing acceptance region

DO
  u = rand(0)
  v = rand(0)
  call random_number (u)
  call random_number (v)
  v = 1.7156 * (v - half)

!     Evaluate the quadratic form
  x = u - s
  y = ABS(v) - t
  q = x**2 + y*(a*y - b*x)

!     Accept P if inside inner ellipse
  IF (q .LT. r1) EXIT
!     Reject P if outside outer ellipse
  IF (q .GT. r2) CYCLE
!     Reject P if outside acceptance region
  IF (v**2 .LT. -4.0*LOG(u)*u**2) EXIT
END DO

!     Return ratio of P's coordinates as the normal deviate
random_normal = v/u
RETURN

END FUNCTION random_normal
! --------------------------------------------------------------------
integer function nhash (id,mode,idtab,isize,maxn)
  integer mode,nprime/433/,ie/16/,i,j,isize,maxn &
    ,idtab(isize),id
!                                   instructions for hashing:
!                           mode=-1 initiates parameters, zeroes idtab
!                               = 0 looks up id only
!                               = 1 puts new id's into idtab
  if (mode == -1) then
    maxn = 0
    do i=1,isize
      idtab(i) = 0
    enddo
    nhash = -1
    return
  endif
  if (id == 0) then
    nhash = 0
    return
  endif
  nhash = abs(mod(nprime*id,isize)) + 1
  do j = 1,isize
    if (idtab(nhash) == 0) then
      if (mode == 1) then
        idtab(nhash) = id
        maxn = maxn + 1
        if (maxn .gt. .9*isize) then
          print *,'hash table size',isize &
           ,'  is too small, increase it'
          stop 261
        endif
      else
        nhash = 0
      endif
      return
    endif
    if (idtab(nhash) == id) return
    nhash = mod(nhash+ie,isize) + 1
  enddo
  print *,'Function NHASH dropped out of the loop'
  stop 41
end
