c
c     Compute inbreeding coefficients
c    add phantom parent group numbers
c   create coded pedigree file too
c
c  Updated: May 28, 2018
c
      Parameter (nam=1000000, nyr=200, ngrp=2000)

      Character*12 aid(nam),cid
      Integer sir(nam),dam(nam),year(nam),brd(nam),
     x kanm,ksir,kdam,inum(nyr),lbc,pgs(20,20),
     x pgd(20,20),npg(nam,3),sirk,damk
      Real*8 fi(nam),fd,fs,b,bii,yobf(nyr,2)
      Real*8 fa
c
c  Open files
c
      itest = 10
      mgrp = 0

      if(itest.gt.0)print *,'Opening files'
      open(9,file='PEDA3S.d',form='formatted',status='old')
      open(12,file='PEDA4.d',form='formatted',status='unknown')
      open(13,file='CPEDA4.d',form='formatted',status='unknown')
      open(17,file='ped4.out',form='formatted',status='unknown')
c
      mam=0
      madd=ngrp
      sir = 0
      dam = 0
      pgs = 0
      pgd = 0
      npg = 0
      year = 0
      aid = '               '
 8    read(9,1001,end=9)kanm,ksir,kdam,kgen,cid
 1001 format(1x,3i10,1x,i5,1x,a12)
      mam = mam + 1
      kanm=kanm+madd
       if(ksir.eq.0)then
          mgrp = (kgen-1)*40+(kbr-1)*2 + 1
          if(mgrp.gt.ngrp)print *,'ERROR IN ped04.f, must fix'
          npg(mgrp,1)=kgen
           npg(mgrp,2)=kbr
           npg(mgrp,3)=1
           pgs(kgen,kbr)=mgrp
         sir(kanm) = mgrp
       else
         sir(kanm) = ksir+madd
       endif

       if(kdam.eq.0)then
          mgrp = (kgen-1)*40 + (kbr-1)*2 + 2
       if(mgrp.gt.ngrp)print *,'ERROR IN ped04.f, must fix'
          npg(mgrp,1)=kgen
          npg(mgrp,2)=kbr
          npg(mgrp,3)=2
          pgd(kgen,kbr)=mgrp
        dam(kanm) = mgrp
       else
        dam(kanm) = kdam+madd
       endif
      year(kanm) = kgen
      brd(kanm) = kbr
      aid(kanm) = cid

      go to 8
c
c  Initialize arrays and subroutines, if necessary
c
 9    inum = 0
      call fxini(nam)
      yobf = 0.0
      fi = 0.0
      write(17,1901)madd,ngrp
 1901 format(1x,i8,'  parent groups formed',i8)
        iss=0
         bii=0.d0
         fs=0.d0
         cid='    '
      im = 0
      do 19 ik=1,ngrp
     
        call fxadd(iss,iss,fa)
          bii=1.d0
c        if(ik.gt.mgrp)go to 19
c        kgen = npg(ik,1)
c        kbr = npg(ik,2)
c        isx = npg(ik,3)
        fi(ik) = fa

        write(12,3938)ik,iss,iss,bii,fa,ik
 3938 format(1x,3i10,1x,f20.12,2x,f20.12,2x,' PG ',i8)
c                 31 32   52   54    74  76   91
       icode=0
       write(13,3966)ik,icode,iss,iss
 19   continue  
      
      if(itest.gt.0)print *,'initializing inbreeding routine'
c  inbreeding routine, initialization
c  Go through pedigrees
c
      madd = ngrp 
      mtall = mam + ngrp
      do 10 k=(ngrp+1),mtall
       kgen = year(k)
       kbr = brd(k)
       kanm = k 
       ksir = sir(kanm) 
       kdam = dam(kanm)
       damk = dam(kanm)
         fs=-1.0
         fd=-1.0
          jsir=0
          jdam=0
         if(ksir.gt.ngrp)then
              fs=fi(ksir)
              jsir=ksir
            endif
         if(kdam.gt.ngrp)then
             fd=fi(kdam)
             jdam=kdam
           endif
         bii = 1.d0
         if(jsir.gt.0.or.jdam.gt.0)then
         b = 0.5 - 0.25*(fs+fd)
         bii = 1.d0/b
         endif
        call fxadd(jsir,jdam,fa)
        fi(kanm) = fa
           
      write(12,3908)kanm,ksir,kdam,bii,fi(kanm),aid(kanm)
 3908 format(1x,3i10,1x,f20.12,2x,f20.12,2x,a12)
c                 31 32   52   54   74   76  88
c  Ped4.d has length 88 -r89  -ka77a12
c
c  write coded pedigree file - to be sorted 
c
      icode=0
      write(13,3966)kanm,icode,ksir,kdam
 3966 format(1x,i10,i3,1x,2i10)
c                11 14 15  35
c      record length is 35 -r36 -ka2a13
        mped = mped + 1
      icode=1
      if(ksir.gt.0)then
         write(13,3966)ksir,icode,kanm,kdam
         mped = mped + 1
         endif
      icode=2
      if(kdam.gt.0)then
         write(13,3966)kdam,icode,kanm,ksir
         mped = mped + 1
        endif 
c
c  sum for averages (By generation groups
c
      jy = year(kanm)
      if(jy.lt.1)jy=1
      if(jy.gt.nyr)jy=200
      yobf(jy,1) = yobf(jy,1) + 1.0
      yobf(jy,2) = yobf(jy,2) + fa
      if(fa.gt.0.0)inum(jy)=inum(jy)+1

 10   continue
      mam = mam + ngrp
      go to 9991
 9001 write(17,1790)
 1790 format(' Number of records in is greater than nam')
 9002  write(17,1795) kanm, ksir, kdam
 1795  format('  Animal ID greater than nam',3i8)
 9003  write(17,1801)mam,kanm
 1801  format(' Input file not sorted correctly?',2i10)
       go to 9995
 9991 write(17,1799)mam
 1799 format(1x,i10,' Number of Animals and PG'/'  Average inbreeding',
     x ' coefficient by generation number')
c
      do 30 jy=1,200
         d = yobf(jy,1)
         if(d.lt.0.5)go to 30
         fd = yobf(jy,2)/d
         kyr=jy
         write(17,7003)kyr,d,fd,inum(jy)
 7003    format(3x,i6,f10.0,d15.6,i10)
 30   continue   

      write(17,7088)mped
 7088 format(1x,i10,' Coded pedigree records - to be sorted')
c
c
      write(17,7092)ngrp
 7092 format(1x,i10,' Phantom Groups')
      close(12)
      close(13)
 9995 write(17,7008)
 7008 format(' %%%%%%%%%%%%%%%% End of ped04.f %%%%%%%%%%%%%%%%')
      call fxfree
      stop
      end