! Program conflict.f90              Calls chktrio to check parentage conflicts
!                                   Counts conflicts by animal and by chip
!                                   Corrects errors and fills missing using ped
!                                   Merges genotypes for clones
!
!                                   Paul VanRaden,  April 2009
!                                   USDA Animal Improvement Programs Lab
!
!            OPTIONS:
!                     Xchrom = chromosome number for X-specific (excluding PAR)
!                     Ychrom = chromosome number for Y
!                     Mchrom = chromosome number for Mitochondria
!                     genout = 0 does not output genotypes.txt
!                            = 1 outputs genotypes.txt
!                     chksnp = 0 does not output SNP or animal error counts
!                            = 1 outputs error files but doubles run time
! COMMAND LINE ARGUMENT:
!                     threads= number of parallel processors to use
  USE IFPORT
  use omp_lib
  implicit none
!include "omp_lib.h"
  integer(kind=1), allocatable :: genotype(:,:), genovec(:)
  integer,         allocatable :: ped(:,:), progeny(:), nmiss(:,:), chiplist(:,:)
  integer,         allocatable :: errloc(:,:), xhet(:,:), check(:,:), females(:), genmale(:,:)
  integer,         allocatable :: animkey(:), nprog(:), pntloc(:)
  integer,         allocatable :: eqnlist(:), ids(:), amiss(:), maxloc(:)
  integer,         allocatable :: genfreq(:,:), erra(:), errs(:), errd(:)
  integer,         allocatable :: errprg(:), errprog(:), errpnt(:)
  integer,         allocatable :: chip(:), loclist(:), nchip(:), chrovec(:)
  integer,         allocatable :: impute0(:), impute2(:), clone(:)
  real*8,          allocatable :: freq(:,:)
  real*8,          allocatable :: aprob(:), sprob(:), dprob(:), progsum(:)
  real*8,          allocatable :: errpct(:), misspct(:), probpnt(:)
  character*1,     allocatable :: sexv(:), chrome_type(:)
  character*3,     allocatable :: oldyoung(:)
  character*1 sex, last
  character*50 markername
  integer*8   location, nfill, mfill, pntmiss, progmiss
  integer genanml, gensire, gendam
  integer erranml, errsire, errdam
  integer fillanml, fillsire, filldam
  real    probanml, probsire, probdam
  integer chips, chipn, genout, chksnp, chkanml
  integer maxped,maxan,maxan2,loci,nloci,aloci,Xloci,Yloci,Mloci
  integer chrome,chromes,idnum,maxn,nhash,Xchrom,Ychrom,Mchrom
  integer n16,called,callm,totprog,maxprog,nclone/0/
  integer missing/5/
  integer loc,locus,nped,errtot,misstot,seg
  integer i,j,k,nn,pnt,pnt1,pnt2,maxpnt,biryr
  real    errate,HW_homo,HW_het,HW_hetm,MAF
  integer :: domain_get_max_threads, get_max_threads
  integer :: command_argument_count, length, status, threads, nthr
  character*2  cthrds
  logical exist
  print *,'Output from conflict.f90'
  print *
  inquire  (file='conflict.options',exist=exist)
  if (.not. exist) then
    print *,'file conflict.options is missing'
    call exit(9)
    stop
    end if
  inquire  (file='chromosome.data',exist=exist)
  if (.not. exist) then
    print *,'file chromosome.data is missing'
    call exit(14)
    stop
    end if
! ---------------------------------------------------------  Input data files
                     open ( 9,file='conflict.options')
                     open (10,file='pedigree.file')
                     open (14,file='chromosome.data')
                     open (16,file='genotypes.raw')
                     open (17,file='clone.keys')
!                                                                Read options
  read (9,*)     Xchrom, Ychrom, Mchrom, genout, chksnp
  print '(5i8)', Xchrom, Ychrom, Mchrom, genout, chksnp
  print *,    '  Xchrom, Ychrom, Mchrom, genout, chksnp'
! ......................................................... Open output files
    if (chksnp == 1) open (20,file='parent.errors')
    if (chksnp == 1) open (21,file='marker.errors')
    if (chksnp == 1) open (22,file='marker.freq')
    if (chksnp == 1) open (51,file='hetx')
    if (genout == 1) open (26,file='genotypes.txt')

! -----------------------------------------------  Default number of threads
  print *
  nthr=omp_get_max_threads()
!                                                      Set number of threads
  if (command_argument_count() > 0) then
    call get_command_argument(1,cthrds,length,status)
    if (cthrds(2:2) == ' ') then
      cthrds(2:2) = cthrds(1:1)
      cthrds(1:1) = '0'
      end if
    read (cthrds,'(i)') threads
    threads = min(nthr, threads)
    print *,'Maximum number of processors set to       ',threads
  else
    threads = nthr
    print *,'Number of processors available and used =',nthr
    print *,'Can set threads by command line argument: conflict 5'
    print *,'Or can be controlled using: export OMP_NUM_THREADS=5'
    end if
  call omp_set_num_threads(threads)
!
  if (threads > 5) print *,'Speed may decrease with more than 5 processors'
  print *
! ----------------------------  Get loci, chromes, chips from chromosome.data
  loci = 0
  chromes = 0
  read(14,*,end=5)
  do
    read(14,*,end=5) markername,chrome,nloci,loc,location,chips
    loci = loci + 1
    chromes = max (chrome, chromes)
    end do
 5 continue
  print *,loci,' loci on ',chromes,' chromosomes' &
         ,' read from chromosome.data'
  rewind 14
  allocate (chiplist(loci,chips), errpct(0:chips), misspct(0:chips))
  allocate (loclist(chips), nchip(0:chips), maxloc(0:chips))
  allocate (chrovec(loci), chrome_type(chromes))
  allocate (impute0(loci), impute2(loci))
  nchip = 0
  loci = 0
  aloci = 0
  Xloci = 0
  Yloci = 0
  Mloci = 0
  maxloc = 0
! ---------------------------------------------------   Read chromosome.data
  read(14,*,end=10)
  do
    loci = loci + 1
    read(14,*,end=10) markername,chrome,nloci,loc,location,chips &
      ,(loclist(chipn),chipn=1,chips)
    do chipn = 1,chips
      if (loclist(chipn) > 0) chiplist(loclist(chipn),chipn) = loci
      maxloc(chipn) = max(maxloc(chipn),loclist(chipn))
      end do
    if (chrome > chromes) then
      print *,chrome,' exceeds # chromosomes ',chromes
      stop
      end if
    chrovec(loci) = chrome
    if (chrome < Xchrom)  aloci = aloci + 1
    if (chrome == Xchrom) Xloci = Xloci + 1
    if (chrome == Ychrom) Yloci = Yloci + 1
    if (chrome == Mchrom) Mloci = Mloci + 1
    if (chrome < Xchrom)  chrome_type(chrome) = 'A'
    if (chrome == Xchrom) chrome_type(chrome) = 'X'
    if (chrome == Ychrom) chrome_type(chrome) = 'Y'
    if (chrome == Mchrom) chrome_type(chrome) = 'M'
    end do
 10 loci = loci - 1
  print *,aloci,' autosomal and pseudo-autosomal loci'
  print *,Xloci,' X-specific loci on chromosome   ',Xchrom
  print *,Yloci,' Y-specific loci on chromosome   ',Ychrom
  print *,Mloci,' Mitochondrial loci on chromosome',Mchrom
  print *
  print *,'Missing genotype code = ',missing
!                                          Get maxan from genotypes.txt
  maxan = 0
  do
    read(16,*,end=12) idnum,chipn,loc
    maxan = maxan + 1
    end do
 12 continue
  rewind 16
  print *
  print *,maxan,' animals in genotype file'
!                                         Get maxped from pedigree.file
  maxped = 0
  do
    read(10,*,end=13) 
    maxped = maxped + 1
    end do
 13 continue
  print *
  print *,maxped,' animals in pedigree file'
  rewind 10
! -----------------------------------------------------------  allocate
  maxan2 = 2 * maxan
  allocate (genotype(loci,0:maxan), probpnt(loci))
  allocate (ped(2,maxped), ids(maxan2), eqnlist(0:maxan2))
  allocate (animkey(maxped), sexv(maxped), nmiss(loci,0:chips))
  allocate (errloc(loci,0:chips), xhet(loci,0:chips), freq(loci,0:chips), errprg(loci))
  allocate (erra(maxan), errs(maxan), errd(maxan), errpnt(loci))
  allocate (aprob(maxan), sprob(maxan), dprob(maxan), amiss(maxan))
  allocate (errprog(0:maxan), progsum(0:maxan), nprog(0:maxan))
  allocate (chip(maxped), genovec(loci), check(loci,0:chips))
  allocate (females(0:chips), genmale(0:2,loci), genfreq(0:2,loci))
  allocate (pntloc(0:maxan), progeny(2*maxan), clone(maxan))
  nn = nhash(0,-1,ids,maxan2,maxn)
! -----------------------------------------------------  Read genotypes
  eqnlist = 0
  n16 = 0
  genfreq = 0
  genmale = 0
  genotype = missing
  print *
  print *,'First 5 animals in genotype file'
  print *,'      IDnum        chip        loci'
!                                                    Read genotypes.txt
  do
    read(16,'(3i10,1x,2000000i1)',end=20) idnum,chipn,loc &
      ,(genovec(locus),locus=1,loc)
    n16 = n16 + 1
    if (n16 <= 5) print *,idnum,chipn,loc
    animkey(n16) = idnum
!   if (idnum == 2914728) chkanml = n16
    chip(n16) = chipn
    nchip(chipn) = nchip(chipn) + 1
    nchip(0) = nchip(0) + 1
    nn = nhash(idnum,1,ids,maxan2,maxn)
    eqnlist(nn) = n16
!                                           Pad low density to max loci
    do locus = 1,loc
      genotype(chiplist(locus,chipn),n16) = genovec(locus)
      end do
    end do
 20 continue
  print *,n16,' records read from genotypes.txt'
  print *
  do chipn = 1, chips
    if (nchip(chipn) > 0) print *,nchip(chipn),'records for chip ',chipn,maxloc(chipn)
    end do
! -------------------------------------------------  Read pedigree data
  nped = 0
  do
    read(10,*,end=40) sex,idnum,pnt1,pnt2,biryr
    nped = nped + 1
    if (nped > maxped) then
      print *,'Increase maxped = ',maxped
      stop
      end if
!                                                       Hash individual
    nn = nhash(idnum,0,ids,maxan2,maxn)
    if (nn > 0) then
      sexv(eqnlist(nn)) = sex
!                                                          Hash parents
      ped(1,eqnlist(nn)) = eqnlist(nhash(pnt1,0,ids,maxan2,maxn))
      ped(2,eqnlist(nn)) = eqnlist(nhash(pnt2,0,ids,maxan2,maxn))
      end if
    end do
 40 print *
  print *,nped,' records read from pedigree file'
! -------------------------------------------------  Count genotyped progeny
  nprog = 0
  maxprog = 0
  do i = 1,maxan
    do pnt = 1,2
!                                                               Parent known
      if (ped(pnt,i) > 0) then
        nprog(ped(pnt,i)) = nprog(ped(pnt,i)) + 1
        if (nprog(ped(pnt,i)) > maxprog) then
          maxprog = nprog(ped(pnt,i))
          maxpnt = animkey(ped(pnt,i))
          end if
        end if
      end do ! pnt
    end do ! i
  print '(i9,a,i9)',maxprog &
    ,' maximum genotyped progeny for parent id = ',maxpnt
  print *
!                                               Set up parent location vector
  totprog = 1
  do i = 1,maxan
    if (nprog(i) > 0) then
      pntloc(i) = totprog
      totprog = totprog + nprog(i)
      end if
    end do ! i
  nprog = 0
!                                                     Store progeny locations
  do i = 1,maxan
    do pnt = 1,2
      if (ped(pnt,i) > 0) then
        progeny(pntloc(ped(pnt,i)) + nprog(ped(pnt,i))) = i
        nprog(ped(pnt,i)) = nprog(ped(pnt,i)) + 1
        end if
      end do ! pnt
    end do ! i
! -------------------------------------------------------  Read list of clones
  do
    read (17,*,end=17) idnum, i, last
    nclone = nclone + 1
    nn = nhash(idnum,0,ids,maxan2,maxn)
    if (nn > 0) then
      clone(i) = eqnlist(nn)
!                                                          Process clone group
      if (last == 'Y') then
        do j = 2, i
          do loc = 1, loci
!                                                                 Fill missing
            if (genotype(loc,clone(1)) == 5) then
              genotype(loc,clone(1)) = genotype(loc,clone(j))
            else
!                                                           Suppress conflicts
              if (genotype(loc,clone(j)) /= 5 .and. &
                  genotype(loc,clone(j)) /= genotype(loc,clone(1))) &
                genotype(loc,clone(1)) = 5
              end if
            end do
            if (maxloc(chip(clone(j))) > maxloc(chip(clone(1)))) &
              chip(clone(1)) = chip(clone(j))
          end do
!                                                                Copy to group
        do j = 2, i
          do loc = 1, loci
            genotype(loc,clone(j)) = genotype(loc,clone(1))
            end do
          chip(clone(j)) = chip(clone(1))
          end do
        end if
    else
      print *,idnum,'animal in clone.keys not in genotype file'
      stop
      end if
    end do
 17 continue
 print *,nclone, 'clones read from clone.keys'
!                                                       Initialize trio arrays
 call trio_setup()
! -----------------------------------------------  Count errors by SNP, animal
 if (chksnp == 1) then
  females = 0
  erra = 0
  errs = 0
  errd = 0
  errtot = 0
  errloc = 0
  xhet = 0
  check = 0
  aprob = 1.
  sprob = 1.
  dprob = 1.
  errprog = 0
  progsum = 0.d0
  nprog = 0
  amiss = 0
  nmiss = 0
  misstot = 0
  freq = 0
  do locus = 1,loci
    nmiss(locus,0) = maxan
    end do
!                                                   Loop across animals
  do i = 1,maxan
    sex = sexv(i)
    chipn = chip(i)
    if (sex == 'F') females(chipn) = females(chipn) + 1
!                                                   Loop across loci
    do loc = 1,maxloc(chipn)
      locus = chiplist(loc,chipn)
      chrome = chrovec(locus)
!                                                   Genotypes of animal, sire, dam
      genanml = genotype(locus,i)
      gensire = genotype(locus,ped(1,i))
      gendam  = genotype(locus,ped(2,i))
!                                                   Check each locus
      call chktrio (genanml, gensire, gendam, sex, chrome_type(chrome)  &
                   ,erranml, errsire, errdam   &
                   ,probanml, probsire, probdam &
                   ,fillanml, fillsire, filldam)
		   
!                                                   Count x het errors
      if (chrome==Xchrom .and. sex=='M' .and. erranml==1 .and. errdam==0) then
        erranml=0
	xhet(locus,0) = xhet(locus,0) + 1
        xhet(locus,chip(i)) = xhet(locus,chip(i)) + 1
      end if
!                                                   Count errors
      erra(i) = erra(i) + erranml
      errs(i) = errs(i) + errsire
      errd(i) = errd(i) + errdam
      errloc(locus,0) = errloc(locus,0) + erranml
      errloc(locus,chip(i)) = errloc(locus,chip(i)) + erranml
      if (probanml > 0.) then
        check(locus,0) = check(locus,0) + 1
        if (chrome /= Ychrom .or. sex == 'M') &
          check(locus,chip(i)) = check(locus,chip(i)) + 1
        end if
!                                                   Multiply prob(error)
      aprob(i) = aprob(i) * (1. - probanml)
      sprob(i) = sprob(i) * (1. - probsire)
      dprob(i) = dprob(i) * (1. - probdam)
!                                                   Count errors for parents
      if (ped(1,i) > 0) then
        errprog(ped(1,i)) = errprog(ped(1,i)) + errsire
        progsum(ped(1,i)) = progsum(ped(1,i)) + probsire
        end if
      if (ped(2,i) > 0) then
        errprog(ped(2,i)) = errprog(ped(2,i)) + errdam
        progsum(ped(2,i)) = progsum(ped(2,i)) + probdam
        end if
!                                                   Count missing
      if (genanml == 5) then
        if (chrome /= Ychrom .or. sex == 'M') &
          nmiss(locus,chip(i)) = nmiss(locus,chip(i)) + 1
        amiss(i) = amiss(i) + 1
      else
        freq(locus,0) = freq(locus,0) + genanml - 1.
        freq(locus,chip(i)) = freq(locus,chip(i)) + genanml - 1.
        nmiss(locus,0) = nmiss(locus,0) - 1
!                                                   Count BB, AB, AA (exclude male X)
        if (chrome < Xchrom .or. sex == 'F') then
          genfreq(genanml,locus) = genfreq(genanml,locus) + 1
        else
!                                                  Separate count for male X
          genmale(genanml,locus) = genmale(genanml,locus) + 1
          end if
        end if
      end do ! locus
! ------------------------------------------------  Report problem animals
    if (sprob(i) == 0) errs(i) = -1
    if (dprob(i) == 0) errd(i) = -1

    if (erra(i) + errs(i) + errd(i) > .02 * (loci - amiss(i))) then
      print *
      print '(i9,3i7,a)',animkey(i),erra(i),errs(i),errd(i) &
           ,' genotype had many anml, sire, or dam conflicts'
      end if
!                                                  -1 indicates parent not genotyped
    if (ped(1,i) <= 0) errs(i) = -1
    if (ped(2,i) <= 0) errd(i) = -1
!                                                   Number of genotyped progeny
    nprog(ped(1,i)) = nprog(ped(1,i)) + 1
    nprog(ped(2,i)) = nprog(ped(2,i)) + 1
    errtot = errtot + erra(i)
    misstot = misstot + amiss(i)
    end do ! i
  print *
  print *,'       chip     females       males        total'
  do chipn = 1, chips
    if (nchip(chipn) > 0) print *,chipn, females(chipn) &
       ,nchip(chipn) - females(chipn), nchip(chipn)
    end do
!
  write(20,'(a)') 'animal    errors errsire errdam nmiss  errate  nprogeny'
  write(22,'(a)') '    locus       BB       AB       AA'
!
!                                                   Output animal.errors
  do i = 1,maxan
    errate = 0.
    if (nprog(i) > 0) errate = float(errprog(i))/nprog(i)

    write(20,'(i9,3i7,i8,f10.3,i7)') animkey(i) &
      ,erra(i),errs(i),errd(i),amiss(i),errate,nprog(i)
    end do ! i
!                                                  Output animals / chip
    write(21,'(a)') '            observations total and on each chip'
    write(21,'(a,999i8)') 'locus     stat',(nchip(j),j=0,chips)
    write(21,*)
!                                                   Output marker.errors

    write(51,'(a,999i8)') 'locus     ', ((nchip(j) - females(j)),j=0,chips)
    
  do locus = 1,loci
    chrome = chrovec(locus)
    do j = 0,chips
      if (nchip(j) - nmiss(locus,j) > 0) then
        freq(locus,j) = .5 + .5 * freq(locus,j) / (nchip(j) - nmiss(locus,j))
        errpct(j) = 100. * errloc(locus,j) / check(locus,j)
        misspct(j) = 100. * nmiss(locus,j) / nchip(j)
      else
        freq(locus,j) = 0.
        errpct(j) = 0
        misspct(j) = 0
        end if
      end do ! j
!                                                     Hardy-Weinberg test
    called = genfreq(0,locus) + genfreq(1,locus) + genfreq(2,locus)
    callm  = genmale(0,locus) + genmale(1,locus) + genmale(2,locus)
    MAF = min(freq(locus,0), 1.-freq(locus,0))
    if (freq(locus,0) > .01 .and. freq(locus,0) < .99  &
       .and. called > 0) then
      HW_het  = 1. * genfreq(1,locus) / called  &
              / (2. * freq(locus,0) * (1. - freq(locus,0)))
      HW_homo = 1. * min(genfreq(0,locus),genfreq(2,locus)) &
              / (called * MAF**2)
      HW_hetm = 1. * genmale(1,locus) / callm  &
              / (2. * freq(locus,0) * (1. - freq(locus,0)))
!                                                       Test X-specific vs PAR
      if (chrome == Xchrom) then
        if (HW_hetm > .3 .and. HW_hetm < 1.5) &
          print *,locus,' locus on X should be on PAR', freq(locus,0)
        end if
      if (chrome == Xchrom - 1) then
        if (HW_hetm < .1) &
          print *,locus,' locus on PAR should be on X', freq(locus,0)
        end if
    else
      HW_het  = 1.
      HW_homo = 1.
      end if
      
    if (chrome==Xchrom) write(51,'(i9,1x,999i8)')   locus,(xhet(locus,j),j=0,chips)
      
    write(21,'(i9,1x,a,999i8)')   locus,'error',(errloc(locus,j),j=0,chips)
    write(21,'(i9,1x,a,999i8)')   locus,'nmiss',(nmiss(locus,j),j=0,chips)
    write(21,'(i9,1x,a,999i8)')   locus,'check',(check(locus,j),j=0,chips)
    write(21,'(i9,1x,a,999f8.4)') locus,' freq',(freq(locus,j),j=0,chips)
    write(22,'(4i9)')     locus,(genfreq(i,locus),i=0,2)
    end do ! locus
  print *
  print '(a,f9.6,a)','Overall missing rate = ',100. * misstot / loci / maxan,' %'
  print '(a,f9.6,a)','Overall error   rate = ',100. * errtot / loci / maxan,' %'
  print *
 end if ! chksnp
!                                             End chksnp section, begin genout
 if (genout == 1) then
  pntmiss = 0
  progmiss = 0
  nfill = 0
  mfill = 0
! ------------------------------------------  Resolve conflicts across animals
  do i = 1,maxan
    sex = sexv(i)
!                                                             Loop across loci
    do loc = 1,maxloc(chip(i))
      locus = chiplist(loc,chip(i))
      chrome = chrovec(locus)
!                                               Genotypes of animal, sire, dam
      genanml = genotype(locus,i)
      gensire = genotype(locus,ped(1,i))
      gendam  = genotype(locus,ped(2,i))
!                                                      Check each parent locus
      call chktrio (genanml, gensire, gendam, sex, chrome_type(chrome)  &
                   ,erranml, errsire, errdam   &
                   ,probanml, probsire, probdam &
                   ,fillanml, fillsire, filldam)
      errpnt(locus) = erranml
      probpnt(locus) = probanml
!                                                        Fill missing genotpes
      if (genanml == 5 .and. fillanml /= 5) then
        genotype(locus,i) = fillanml
        nfill = nfill + 1
        end if
!                                                          Zero progeny counts
      errprg(locus) = 0
      impute0(locus) = 0
      impute2(locus) = 0
      end do ! loc

!$OMP PARALLEL DO REDUCTION(+: errprg,progmiss,impute0,impute2) &
!$OMP             PRIVATE(k,sex,loc,locus,chrome,genanml,gensire,gendam, &
!$OMP             erranml,errsire,errdam,probanml,probsire,probdam, &
!$OMP             fillanml,fillsire,filldam)
!                                                    Process progeny of animal
    do j = pntloc(i), pntloc(i) + nprog(i) - 1
!                                                               Sex of progeny
      k = progeny(j)
      sex = sexv(progeny(j))
!                                                             Loop across loci
      do loc = 1, maxloc(chip(k))
        locus = chiplist(loc,chip(k))
        chrome = chrovec(locus)
!                                                Genotypes of progeny and mate
        genanml = genotype(locus,k)
        gensire = genotype(locus,ped(1,k))
        gendam  = genotype(locus,ped(2,k))
!                                                     Check each progeny locus
        call chktrio (genanml, gensire, gendam, sex, chrome_type(chrome)  &
                     ,erranml, errsire, errdam   &
                     ,probanml, probsire, probdam &
                     ,fillanml, fillsire, filldam)
!                                                   Suppress progeny conflicts
!                                            Errors from this parent, not mate
        if (sexv(i) == 'M') then
          if (errsire == 1 .and. errdam /= 1) then
            errprg(locus) = errprg(locus) + 1
            if (errpnt(locus) == 0 .and. probpnt(locus) > 0) then
              genotype(locus,k) = missing
              progmiss = progmiss + 1
              end if
            end if
          if (gendam /= missing .and. gendam /= 1) then
            if (fillsire == 3) impute0(locus) = impute0(locus) + 1
            if (fillsire == 4) impute2(locus) = impute2(locus) + 1
            end if
          end if
        if (sexv(i) == 'F') then
          if (errdam == 1 .and. errsire /= 1) then
            errprg(locus) = errprg(locus) + 1
            if (errpnt(locus) == 0 .and. probpnt(locus) > 0) then
              genotype(locus,k) = missing
              progmiss = progmiss + 1
              end if
            end if
          if (gensire /= missing .and. gensire /= 1) then
            if (filldam == 3) impute0(locus) = impute0(locus) + 1
            if (filldam == 4) impute2(locus) = impute2(locus) + 1
            end if
          end if
        end do ! loc
      end do ! j
!END DO PARALLEL
!                                                      Loop across parent loci
    do loc = 1, maxloc(chip(i))
      locus = chiplist(loc,chip(i))
      chrome = chrovec(locus)
!                                                     Resolve parent conflicts
      if (errpnt(locus) == 1 .or. errprg(locus) > .1 * nprog(i)) then
        genotype(locus,i) = missing
        pntmiss = pntmiss + 1
        end if
!                                                        Fill missing genotpes
!                                                                 heterozygous
!                                                          (Require 2 of each)
      if (impute0(locus) >= 2 .and. impute2(locus) >= 2) then
        genotype(locus,i) = 1
        mfill = mfill + 1
        end if
!                                                                   homozygous
      if (impute0(locus) >= 10 .and. impute2(locus)*10 < impute0(locus)) then
        genotype(locus,i) = 0
        mfill = mfill + 1
        end if
      if (impute2(locus) >= 10 .and. impute0(locus)*10 < impute2(locus)) then
        genotype(locus,i) = 2
        mfill = mfill + 1
        end if
      errpnt(locus) = 0
      probpnt(locus) = 0.
      end do ! loc
    end do ! i
  print *,pntmiss, 'parent conflicts set to missing'
  print *,progmiss, 'progeny conflicts set to missing'
  print *,nfill, 'missing genotypes imputed from parents'
  print *,mfill, 'missing genotypes imputed from progeny'
!                                                         Write genotypes.txt
  do i = 1,maxan
!                                                  Shrink loci to low density
    do locus = 1, maxloc(chip(i))
      genovec(locus) = genotype(chiplist(locus,chip(i)),i)
      end do
    write(26,'(3i10,1x,2000000i1)') animkey(i),chip(i),maxloc(chip(i)) &
      ,(genovec(locus),locus=1,maxloc(chip(i)))
    end do ! i
  print *,maxan,' records written to genotypes.txt'
 end if
  print *
  print *,'Program conflict.f90 completed normally'
  stop
  end
