************************************************************************
*        RENUMMAT - data preparation program for sire and animal models
* This version for maternal effects
*
*  This program:
* 1. Recodes single or multiple fields into consecutive numbers,
* 2. recodes the animals for the animal model:
*    a) validates animals and parents based on their year of birth
*       (if available)
*    b) removes most nocontributing animals,
*    c) can assign unknown parent groups for unknown parents and
*       noncontributing animals
*    d) assign consecutive numbers in the following order: animals with
*       records, animals with pedigrees, other animals, unknown parent groups,
*    e) provides a sample parameter file for program JAA (and related programs).
* 
* New addition: parameter file can contain an entry for the permanent
*               environmental effect
*
* Names of data files and parameters to the program are entered interactively.
* The data files should be in free format, where fields are separated by
* space(s). Animals can have up to 16 characters, other fields should
* be numeric with 9 decimal digits at the most. Upon return, the program
* returns the following files:
*  1. renumbered data file
*  2. renumber information for nonanimal factors
*  3. (if animal effect) pedigree file with the following fields:
*     a) animal number (from 1)
*     b) parent 1 number or unknown parent group number for parent 1
*     c) parent 2 number or unknown parent group number for parent 2
*     d) 3 minus number of known parents
*     e) known or estimated year of birth (0 if not provided)
*     f) number of known parents (parents might be eliminated if not
*        contributing
*     g) number of records (before elimination due to other effects)
*     h) number of progeny (incl. noncontributing animals)
*     i) original animal id
*  4. message file renum.msg
*  5. sample parameter file for program JAA, renum.jaa
*
*
*  Memory requirements
*  The program needs about 50 bytes of memory per animal. The size of memory
*  allocated in bytes is determined by the parameter maxmem. Its value equal
*  to 100000 is sufficient for approximately 2000 animals.
*
*  Disclaimer
*  This program was developed for my own research. It is believed to be
*  reasonably free of serious errors, but few attempts have been made to make
*  it resistant to errors in parameters and input files. Please use it at
*  your own risk.
*
*  last correction:, 06/06/99
*  Copyright 1987-1999 I. Misztal, University of Illinois
************************************************************************
      program renum
      integer maxx,m1,m2,m3,m4,maxmem
 
c memory assigned to the largest variables in the program; in bytes
      parameter (maxmem=250000)

      parameter(maxx=maxmem/12,m1=20,m2=40,m3=9,m4=maxmem/m3/6)
      integer x(maxx),xani(m3,m4),table(m1,m2),pass(m2),kc(m2),
     + cutlev(m1),i,nel,nfact,nelmax,rows,posmax,anitab(4,2),
     + groups(100),ngroups,genint(3),isrpt,ismat
      character name(5)*30,c(m2)*30
 
      print 10,int(m4*.9),maxx/5
10    format(' RENUMMAT - renumbering tool'/' has capacity to renumber'
     +       ,' up to approximately',i8, ' animals'/' and ',i8,
     +       ' other levels')
      print*,'input data file?'
      call readuntil('#',kc,c,m2,nel,1)
      name(1)=c(1)
      print*,'output data file?'
      call readuntil('#',kc,c,m2,nel,1)
      name(2)=c(1)
      print*,'printout data name?'
      call readuntil('#',kc,c,m2,nel,1)
      name(3)=c(1)
      open(1,file=name(1))
      open(2,file=name(2))
      open(3,file=name(3))
      open(10,file='renum.jaa')
      open(11,file='renum.msg')
 
 
      nelmax=0
      nfact=0       
      posmax=0
      print*, 'give positions of colums to be renumbered into single'
     *,' factor, 0 ends'
100   call readuntil('#',kc,c,m2,nel,-1)
      if (nel.ne.0 .and. kc(1).ne.0) then
          if (nel.gt.m2+2) then
              print*,'too many items, increase parameter m2'
              stop
          endif
          nfact=nfact+1
          if (nfact.gt.m1) then
              print*,'too many renumberings, increase parameter m1'
              stop
          endif
          nelmax=max(nelmax,nel)
          table(nfact,1)=nel
          do 110 i=1,nel
             posmax=max(posmax,kc(i))
110          table(nfact,i+1)=kc(i)
          goto 100
      endif
      
115   print*,'give minimum number of observations per level for each'
      print*,'    factor, 0 if all observations retained'
      call readuntil('#',cutlev,c,m2,nel,0)
      if (nfact.ne.nel) then
         print*,nfact,' numbers expected, missing numbers assigned 0'
         do i=nel+1,nfact
            cutlev(i)=0
         enddo
      endif
            
      print*,'give positions of columns to pass unchanged'
      call readuntil('#',kc,c,m2,nel,0)
      pass(1)=nel
      do 120 i=1,nel
         posmax=max(posmax,kc(i))
120      pass(i+1)=kc(i)
 
c animal renumbering part
      print*,'give position of animal in data file, 0 if none'
      call readuntil('#',kc,c,m2,nel,1)
      anitab(1,1)=kc(1)
      if (anitab(1,1).gt.0) then
         print*,'name of the pedigree file'
         call readuntil('#',kc,c,m2,nel,1)
         name(4)=c(1)
         print*,'name of the output pedigree file?'
         call readuntil('#',kc,c,m2,nel,1)
         name(5)=c(1)
         open(8,file=name(4))
         open(9,file=name(5))
         print*,'positions of animal, sire, dam and year of birth'
         print*,'in the pedigree file, year of birth = 0 if missing'
         call readuntil('#',kc,c,m2,nel,4)
         do i=1,4
            anitab(i,2)=kc(i)
         enddo
         if (anitab(4,2).ne.0) then
           print*,'give years to separate unknown parent groups, or 0'
           call readuntil('#',groups,c,m2,ngroups,-1)
           if (groups(1).eq.0) ngroups=0
           print*,
     +       'give minumum, average, and maximum generation interval'
           call readuntil('#',kc,c,m2,nel,3)
           do i=1,3
              genint(i)=kc(i)
           enddo
          else
           ngroups=0
           genint(3)=0
         endif
      endif                        
      print*, 'Is it a repeatability model? (1=yes/0=no)'
      call readuntil('#',kc,c,m2,nel,1)
      isrpt=kc(1)
      if (isrpt.gt.1 .or. isrpt.lt.0) isrpt=0
 
      print*, 'Is it a maternal model? (1=yes/0=no)'
      call readuntil('#',kc,c,m2,nel,1)
      ismat=kc(1)
      if (isrpt.gt.1 .or. isrpt.lt.0) ismat=0

      print*, 'data files read to 200 column only'
 
      nelmax=nelmax+2
      rows=maxx/nelmax
      call ren(x,rows,nelmax,table,m1,m2,nfact,pass,cutlev,
     +         groups,ngroups,m3,m4,xani,name,anitab,genint,isrpt,ismat)
      end
 
      subroutine ren(x,rows,cols,table,m1,m2,nfact,pass,
     *cutlev,groups,ngroups,m3,m4,xani,name,anitab,genint,isrpt,ismat)
      integer m1,m2,rows,cols,nfact,cutlev(m1),m3,m4,
     +        xani(m3,m4),ngroups,groups(ngroups+1),nbuf,notachar,
     +        anitab(4,2),msglev,genint(3)
      parameter(nbuf=100)
      integer x(cols,rows),table(m1,m2),pass(m2),buf(nbuf),buf1(nbuf),n,
     +        bufb(4,nbuf)
      integer dat(100),scol(40),i,j,k,k1,l,m,hashves,nwritt
      integer ifactor,lastfactor,ix,ix1,ixani,nfactor(20),
     +        nrecord,nped,nparent,nelim,nantot,
     +        nanim,nfact1,nrepeat,kpar1,kpar2,upg,grstat(nbuf+1),
     +        known,iy1,isani,isrpt,ismat
      character bufa(nbuf)*16,name(5)*(*)

c fields for xani
      integer nyob,nrec,nprog,npar,ncons,nmat
      parameter (nyob=4,nrec=nyob+1,nprog=nrec+1,npar=nprog+1,
     +           ncons=npar+1,nmat=ncons+1)
      character a*150
      logical rd1,goodid
      data msglev/0/,ixani/0/,nrecord/0/,nped/0/,nmatunk/0/,nelim/0/,
     +     nrepeat/0/

 
      call assignchar(notachar,'%%%%')
 
c initialize
      call zero(cols*rows,x)
      call zero(m3*m4,xani)
      call zero(ngroups+1,grstat)
      n=0
      ix=0
      if (anitab(1,1).ne.0) then
         isani=1
        else
         isani=0
      endif
 
      if (nbuf.le.ngroups .and. isani.eq.1) then
         print*,'nbuf must be greater than ngroups:',nbuf,ngroups
         stop
      endif

      nanim=0
20    if (rd1(1,buf,bufa,nbuf,m)) then
          call convert(bufa,bufb,m)
          nrecord=nrecord+1
c        hash regular factors
          do 30 i=1,nfact
            dat(1)=i
            do 40 j=1,table(i,1)
40            dat(j+1)=buf(table(i,j+1))
            do 50 j=table(i,1)+1,cols-1
50            dat(j+1)=0                 
            k= hashves(dat,cols-1,x,rows,cols,1,ix)
30          x(cols,k)=x(cols,k)+1
c           hash animals  from the production file
            if (anitab(1,1).gt.0) then
               k=hashves(bufb(1,anitab(1,1)),3,xani,m4,m3,1,ixani)
               if (xani(nrec,k) .eq. 0) then
!                 new animal
                  xani(nrec,k)=xani(nrec,k)+1
                  nanim=nanim+1
                  xani(ncons,k)=nanim
               endif
            endif
          goto 20
      endif

      if (nrecord.eq.0) then
         print '('' file '',a,''  empty'')',name(1)
         stop
      endif
 
      if (msglev.gt.0) then
          print*,' after reading data file'
          call printhash(xani,nyob-1,m3,m4)
      endif
 
 
      if (isani.eq.1) then
         nparent=nanim
c hash animals from the pedigree file
80    if (rd1(8,buf,bufa,nbuf,m)) then
          call convert(bufa,bufb,m)
          nped=nped+1
            if (anitab(1,2).gt.0) then
c          animal
               k1=hashves(bufb(1,anitab(1,2)),3,xani,m4,m3,1,ixani)
               if (anitab(4,2).gt.0) xani(nyob,k1)=buf(anitab(4,2))
c          first parent
               if (goodid(bufa(anitab(2,2)))) then
                  k=hashves(bufb(1,anitab(2,2)),3,xani,m4,m3,1,ixani)
                  xani(nprog,k)=xani(nprog,k)+1
                  xani(npar,k1)=xani(npar,k1)+1
c            progeny's year of birth is the lowest of its parents - genint(2)
c                 assigned only if it does not have its own pedigree entry
                  if (xani(nyob,k).eq.0) then
                      xani(nyob,k)=buf(anitab(4,2))-genint(2)
                    else
                    xani(nyob,k)=min(xani(nyob,k),
     +                                    buf(anitab(4,2))-genint(2))
                  endif
               endif
c           second parent  (assumed dam for maternal effect)
               if (goodid(bufa(anitab(3,2)))) then
                  k=hashves(bufb(1,anitab(3,2)),3,xani,m4,m3,1,ixani)
                  xani(nprog,k)=xani(nprog,k)+1
                  xani(npar,k1)=xani(npar,k1)+1
                  if (xani(nyob,k).eq.0) then
                      xani(nyob,k)=buf(anitab(4,2))
                    else
                       xani(nyob,k)=min(xani(nyob,k),
     +                                    buf(anitab(4,2))-genint(2))
                  endif
                  if (ismat.eq.1) then
                     if (xani(ncons,k).eq.0.and.xani(nrec,k1).gt.0) then
                        nparent=nparent+1
                        xani(ncons,k)=nparent
                     endif
                     xani(nmat,k1)=xani(ncons,k)
                   endif
               else
                  if (ismat.eq.1 .and. xani(nrec,k1).gt.0) then
                     if (xani(nmat,k1).eq.0) then
                        nmatunk=nmatunk+1
                        xani(nmat,k1)=-nmatunk
                      endif
                  endif
               endif

            endif
          goto 80
      endif

      if (nped.eq.0 .and. isani.eq.1) then
         print '('' WARNING: file '',a,'' empty'')',name(4)
      endif
 
      if (msglev.gt.0) then
          print*,' after reading pedigree file'
          call printhash(xani,nyob-1,m3,m4)
      endif
 
c remove obviously noncontributing animals that do not have records and:
c   1. have no progeny, or
c   2. have 1 progeny but no parents
C In maternal models, removed animals cannot be dams of animals with records
      do i=1,m4
       if (xani(1,i).ne.0 .and. xani(nrec,i).eq.0 .and. 
     +            xani(ncons,i).eq.0 .and. xani(nmat,i).eq.0) then
          if (xani(nprog,i).eq.0 .or.
     +        xani(nprog,i).eq.1 .and. xani(npar,i).eq.0) then
              xani(1,i)=notachar
              nelim=nelim+1
          endif
        endif
      enddo
      nantot=ixani-nelim
 
      if (msglev.gt.0) then
          print*,' after removing redundant animals'
          call printhash(xani,nyob-1,m3,m4)
      endif

      endif
 
 
c cluster all nonzero x elements to speed-up sorting
100   call clusters(x,rows,cols,ix)
 
c quicksort
      do 105 i=1,cols
105      scol(i)=i
      call sortqms(x,rows,cols,ix,cols,scol,cols)
 
c print the renumber table. In occurences column place the value of renumbered
c levels
      ifactor=0
      write(3,200)nrecord
200   format(' records',i12/
     * //' factor    renumbered level      #      original levels',
     * /1x,60('-'))
      do 210 i=1,ix
        if (ifactor.ne.x(1,i)) then
            ifactor=x(1,i)
            lastfactor=table(ifactor,1)+1
            nfactor(ifactor)=0
            write(3,'(i2,2h :,10i3)')ifactor,(table(ifactor,j),
     *       j=2,lastfactor)
        endif
        if (x(cols,i).ge.cutlev(ifactor)) then
           nfactor(ifactor)=nfactor(ifactor)+1
           write(3,'(10x,10i12)')nfactor(ifactor),x(cols,i),
     *            (x(j,i),j=2,lastfactor)
           x(cols,i)=nfactor(ifactor)
          else
           write(3,'(22x,10i12)')x(cols,i),(x(j,i),j=2,lastfactor)
           x(cols,i)=0
        endif    
                
210     continue
c store and rehash the  table            
      open(20,form='unformatted')
      do 300 i=1,ix
300     write(20)(x(j,i),j=1,cols)
      do 310 i=1,rows
      do 310 j=1,cols
310     x(j,i)=0
      rewind 20
      ix1=0
      do 320 i=1,ix
        read(20)(dat(j),j=1,cols)
        k= hashves(dat,cols-1,x,rows,cols,1,ix1)
320     x(cols,k)=dat(cols)
                               
 
      write(11,330)name(1),name(2)
330   format(' input data file: ',a/' output data file: ',a)
      if (isani.gt.0) then
         write(11,340)name(4),name(5)
      endif
340   format(' input pedigree file: ',a/' output pedigree file: ',a)
 
 
c read the data again and renumber levels
      rewind 1
      nwritt=0
      nfact1=nfact
      if (anitab(1,1).ne.0) nfact1=nfact+1
 
400   if (rd1(1,buf,bufa,nbuf,m)) then
          call convert(bufa,bufb,m)
          do 410 i=1,nfact
            dat(1)=i
            do 405 j=1,table(i,1)
405           dat(j+1)=buf(table(i,j+1))
            do 407 j=table(i,1)+1,cols-1
407           dat(j+1)=0                 
            l= hashves(dat,cols-1,x,rows,cols,0,ix1)
            if (x(cols,l).eq.0) goto 400
410         buf1(i)=x(cols,l)
          do 420 j=1,pass(1)
420         buf1(nfact1+j)=buf(pass(j+1))
c         animal effect
          if (anitab(1,1).ne.0) then
             k=hashves(bufb(1,anitab(1,1)),3,xani,m4,m3,0,iy1)
             if (k.eq.0) then
                   print*,'cannot find animal in pedigree file:'
                   print*, bufa(anitab(1,1))
                   stop
             endif
             l=xani(ncons,k)
             if (l.eq.0) then
                 print*,'This should not happen'
                 stop
             endif
             buf1(nfact+1)=l
             if (ismat.eq.1) then
                 l=xani(nmat,k)
                 if (l.lt.0) then
                     l=abs(l)+nantot
                    else if (l .eq.0) then
                     print*,'Pedigree for animal with record missing: ',
     +                      bufa(anitab(1,1))

                 endif
                 buf1(nfact1+pass(1)+ismat)=l
             endif
           endif
          call pack(a,m,buf1,nfact1+pass(1)+ismat)
          write(2,'(150a1)')(a(j:j),j=1,m)
          nwritt=nwritt+1
          goto 400   
      endif
 
      j=0
      do i=1,nfact
         j=j+nfactor(i)
      enddo
      write(11,450)nrecord,ix,j
450   format(' data records:',t40,i8/' original levels:',t40,i8/
     +       ' levels after cutting:',t40,i8)
      write(11,'(/1x,''Renumbered factors'')')
      write(11,'(1x,''factor'',t20,''levels'')')
      do i=1,nfact
         write(11,'(1x,i5,t20,i8)')i,nfactor(i)
      enddo
      write(11,'(//)')
      if (isani.eq.0) goto 525
 
 
      if (msglev.gt.0) then
          print*,' after writing data output file'
          call printhash(xani,nyob-1,m3,m4)
      endif
 
 
c assign consecutive numbers for animals with pedigree records
      rewind 8
460    if (rd1(8,buf,bufa,nbuf,m)) then
          call convert(bufa,bufb,m)
            k=hashves(bufb(1,anitab(1,2)),3,xani,m4,m3,0,ixani)
            if (k.ne.0) then
               if (xani(ncons,k).eq.0) then
                  nparent=nparent+1
                  xani(ncons,k)=nparent
               endif
            endif
          goto 460
      endif
 
c write relationship file
      rewind 1
      rewind 8
500    if (rd1(8,buf,bufa,nbuf,m)) then
            call convert(bufa,bufb,m)
            k=hashves(bufb(1,anitab(1,2)),3,xani,m4,m3,0,ixani)
            if (k.eq.0) goto 500
c           don't do repeated assignments
            if (xani(ncons,k).lt.0) then
               nrepeat=nrepeat+1
              else
               kpar1=hashves(bufb(1,anitab(2,2)),3,xani,m4,m3,0,ixani)
               kpar2=hashves(bufb(1,anitab(3,2)),3,xani,m4,m3,0,ixani)
               known=1
               if (xani(ncons,k).ne.0) then
                  buf(1)=xani(ncons,k)
                 else
                  nparent=nparent+1
                  buf(1)=nparent
                  xani(ncons,k)=nparent
               endif
               xani(ncons,k)=-xani(ncons,k)
 
               if (kpar1.ne.0) then
                  if (xani(ncons,kpar1).ne.0) then
                     buf(2)=abs(xani(ncons,kpar1))
                    else
                     nparent=nparent+1
                     buf(2)=nparent
                     xani(ncons,kpar1)=nparent
c                   verify generation intervals, if given
                     call checkgen(k,kpar1,xani,nyob,genint,m3,m4)
                   endif
                  else
                   buf(2)=nantot+nmatunk+
     +                     upg(xani(nyob,k),groups,ngroups,grstat)
                   known=known+1
                endif
 
                if (kpar2.ne.0) then
                  if (xani(ncons,kpar2).ne.0) then
                     buf(3)=abs(xani(ncons,kpar2))
                    else
                     nparent=nparent+1
                     buf(3)=nparent
                     xani(ncons,kpar2)=nparent
                     call checkgen(k,kpar2,xani,nyob,genint,m3,m4)
                   endif
                  else
                   if (xani(nmat,k).eq.0) then
                     buf(3)=nantot+nmatunk+
     +                    upg(xani(nyob,k),groups,ngroups,grstat)
                     known=known+1
                    else
                     buf(3)=nantot+abs(xani(nmat,k))
                     k1=nantot+nmatunk+upg(xani(nyob,k)-
     +                             genint(2),groups,ngroups,grstat)
c                    line doubled to have correct count of groups
                     k1=nantot+nmatunk+upg(xani(nyob,k)-
     +                             genint(2),groups,ngroups,grstat)
                     write(9,'(3i8,i2)')nantot+abs(xani(nmat,k)),k1,k1,3
                   endif
                endif
                buf(4)=known
                buf(5)=xani(nyob,k)
                buf(6)=xani(npar,k)
                buf(7)=xani(nrec,k)
                buf(8)=xani(nprog,k)
                call pack(a,m,buf,8)
                write(9,'(1x,a,1x,3a4)')a(1:m),(xani(i,k),i=1,3)
              endif
          goto 500
      endif
 
      if (msglev.gt.0) then
          print*,' after writing pedigree output file'
          call printhash(xani,nyob-1,m3,m4)
      endif
 
c create relationships for remaining animals
      do i=1,m4
         if (xani(1,i).ne.0 .and. xani(ncons,i).ge.0 .and.
     +                                xani(1,i).ne.notachar) then
              if (xani(ncons,i).ne.0) then
                  buf(1)=xani(ncons,i)
                else
                  nparent=nparent+1
                  buf(1)=nparent
                  xani(ncons,i)=nparent
              endif
              buf(2)=nantot+nmatunk+
     +                         upg(xani(nyob,i),groups,ngroups,grstat)
c             line doubled for correct count of UPG
              buf(2)=nantot+nmatunk+
     +                         upg(xani(nyob,i),groups,ngroups,grstat)
              buf(3)=buf(2)
              buf(4)=3
              buf(5)=xani(nyob,i)
              buf(6)=xani(npar,i)
              buf(7)=xani(nrec,i)
              buf(8)=xani(nprog,i)
              call pack(a,m,buf,8)
              write(9,'(1x,a,1x,3a4)')a(1:m),(xani(j,i),j=1,3)
              if (buf(7).gt.0 .and. buf(5).le.0
     +                      .and. anitab(4,2).ne.0) then
                  write(11,470)(xani(j,i),j=1,3)
470               format(' missing pedigree or missing year of birth'
     +                   ,' in animal:',3a4)
             endif
          endif
       enddo
 
      if (msglev.gt.0) then
          print*,' after writing extra pedigree entries'
          call printhash(xani,nyob-1,m3,m4)
      endif
 
      write(11,505)nped,ixani,nelim,nrepeat,1,nanim,nanim+1,
     +        nparent,nparent+1,nparent+nmatunk, nparent+nmatunk+1,
     +        nparent+nmatunk+ngroups+1
505   format(///' Animal factor summary'//
     +       ' pedigree records:',t40,i8/
     +       ' animals found:',t40,i8/' animals eliminated:',t40,i8/
     +       ' repeated pedigrees:',t40,i8//
     +       ' Consecutive number assignment'//
     +       ' animals with records',t40,i8,'  -',i8/
     +       ' parents with pedigrees or ties',t40,i8,'  -',i8/
     +       ' parents for maternal effect',t40,i8,'  -',i8/
     +       ' unknown parent groups',t40,i8,'  -',i8/)
 
      if (anitab(4,2).ne.0 .and. ngroups.gt.0) then
         write(11,510)
510      format(//' Unknown parent group allocation'//
     +     t5,' group no.',t16,'consecutive no',t44,'#',t52,' years')
         write(11,520)1,nantot+nmatunk+1,grstat(1),0,groups(1)
         do i=2,ngroups
            write(11,520)i,nantot+nmatunk+i,grstat(i),groups(i-1)+1,
     +                                 groups(i)
         enddo
         write(11,520)ngroups+1,nantot+nmatunk+ngroups+1,
     +                grstat(ngroups+1),groups(ngroups)+1
      endif
520   format(3i15,5x,i4,1h-,i4)
 
create JAA parameter file

      write(10,'(a)')name(2),name(5),'jaaout'
      write(10,'(2i4,f6.2,i4)')nfact+isani+isrpt+ismat,nfact+1+isani,
     +                         isani*.8,10
      do i=1,nfact
         write(10,'(i4,i8,f8.2)')i,nfactor(i),0.0
      enddo
      if (isani.eq.1) then
         if (isrpt.eq.1)write(10,'(i4,i8,f8.2)')nfact+1,nanim,2.0
         write(10,'(i4,i8,f8.2)')nfact+1,nparent+nmatunk+ngroups+1,-2.0
      endif
      if (ismat.eq.1) then
         write(10,'(i4,i8,f8.2)')nfact1+pass(1)+1,
     +                            nparent+nmatunk+ngroups+1,-2.0
      endif

      print*,'nantot,nparent',nantot,nparent
c Create Blupf90 parameter file
      close(10)
      open(10,file='renum.b90')
      write(10,'(''# Generated by RENUMMAT; must be edited'')')
      write(10,'(a/a)')'DATAFILE',name(2)
      write(10,'(a/i2)')'NUMBER_OF_TRAITS',1
      write(10,'(a/i2)')'NUMBER_OF_EFFECTS',nfact+isani+isrpt+ismat
      write(10,'(a/i2)')'OBSERVATION(S)',nfact+1+isani
      write(10,'(a/)')'WEIGHT(S)'
      write(10,'(2a)')'EFFECTS: POSITIONS_IN_DATAFILE NUMBER_OF_LEVELS',
     +                ' TYPE_OF_EFFECT  [EFFECT NESTED]'
      do i=1,nfact
         write(10,'(i4,i8,'' cross'')')i,nfactor(i)
      enddo
      if (isani.eq.1) then
         if (isrpt.eq.1)write(10,'(i4,i8,'' cross'')')nfact+1,nanim
         write(10,'(i4,i8,'' cross'')')nfact+1,nparent+nmatunk+ngroups+1
      endif
      if (ismat.eq.1) then
         write(10,'(i4,i8,'' cross'')')nfact1+pass(1)+1,
     +                            nparent+nmatunk+ngroups+1
      endif
      write(10,'(a/f5.1)')'RANDOM_RESIDUAL VALUES',1.0

      if (isani.eq.1) then
         if (isrpt.eq.1)then
              write(10,'(a/i2)')'RANDOM_GROUP',nfact+1
              write(10,'(a/a/a/)')'RANDOM_TYPE','diagonal','FILE'
              write(10,'(a/f5.1)')'(CO)VARIANCES',1.0
         endif
         j=nfact+1+isrpt
         write(10,'(a/2i2)')'RANDOM_GROUP',(i,i=j,j+ismat)
         write(10,'(a/a/a/a)')'RANDOM_TYPE','add_an_upg','FILE',name(5)
         write(10,'(a)')'(CO)VARIANCES'
         do i=1,1+ismat
            write(10,'(2f5.1)')(1.0,j=1,1+ismat)
         enddo
      endif

525   print 530
      write(11,530)
530   format(//' Example parameter file for program JAA was stored in'
     +       ,' file renum.jaa'/' Example parameter file for program',
     +        ' blupf90 was stored in file renum.b90'/
     +        ' Both files must be edited before use')

      print 540
540   format(' Summaries and messages were stored in file renum.msg')
 
      end
 
      subroutine nums(a,x,xc,m,n)
c separates line a into items delimited by blanks. character elements 
c kept in a character vector xc, decoded numeric values in vector x.
c version corrected for negative numbers 4/14/93
      integer m,x(m),decode,d,i,j,e,b,lena,n
      character a*(*),xc(m)*(*),buf*1,sp*1,sp12*20
      data sp/' '/,sp12/'                    '/
      d=0
      n=0
      lena=len(a)
      do 5 i=1,m
5       xc(m)=sp12
      j=1
10    buf=a(j:j)
      if (buf.eq.sp) then
        if (d.ne.0) then
          e=j-1
          n=n+1
           xc(n)=a(b:e)
          if (e-b.lt.10)  x(n)=decode(xc(n),e-b+1)
          d=0
        endif
        else
        if (d.eq.0) then
          b=j
          d=1
        endif
      endif
      j=j+1
      if (j.le.lena) goto 10
      return
      end
 
      integer function decode(buf,b)
      integer i,j,dd,b0,b,sign
      character*1 d0,buf*(*),minus
      data d0/'0'/,minus/'-'/
      decode=0
      b0=ichar(d0)
      if (buf(1:1).eq.minus) then
           sign=-1
           j=2
         else
           sign=1
           j=1
      endif
      do 10 i=j,b
      dd=ichar(buf(i:i))-b0
      if (dd.ge.0 .and. dd.le.9) then
        decode=10*decode+dd
      else
        goto 20
      endif
10    continue
20    decode=decode*sign
      return
      end
 
      subroutine pack(a,m,x,n)
c packs an x(n) vector of integers into a*80 character variable. Variables
c there are separated by spaces
      character*80 a
      integer n,x(n),y(10),xi,iz,i,iy,j,m
      iz=ichar('0')
      m=0
      do 10 i=1,n
        xi=x(i)
        if (xi.lt.0) then
          xi=-xi
          m=m+1
          a(m:m)='-'
        endif
        iy=0
20      iy=iy+1
        y(iy)=mod(xi,10)
        xi=xi/10
        if (xi.ne.0) goto 20
        do 30 j=iy,1,-1
          m=m+1
30        a(m:m)=char(y(j)+iz)
        if (i.ne.n) then
          m=m+1
          a(m:m)=' '
        endif
10      continue                    
      return
      end
 
      subroutine sortqms(a,n,m,ix,iy,col,icol)
c sorts matrix of integers by icol columns given by col, in ascending
c order. Matrix a is n x m, only first ix rows and iy columns are sorted
c !! same as sortqm except rows and columns switched ''
      integer i,j,s,l,r,j1,j2,k,a1,a2,iy,icol,ix,m,n
      integer a(m,n),col(icol),x(20)
      integer stack(50,2) 
      s=1
      stack(1,1)=1
      stack(1,2)=ix
10    l=stack(s,1)
      r=stack(s,2)
      s=s-1
20    i=l
      j=r
      j1=(l+r)/2
      do 25 k=1,icol
25    x(k)=a(col(k),j1)
30    do 35 k=1,icol
      a1=a(col(k),i)
      a2=x(k)
      if (a1.lt.a2) then
      i=i+1
      goto 30
      else if(a1.gt.a2) then 
      goto 40
      endif
35    continue      
40    do 45 k=1,icol
      a1=a(col(k),j)
      a2=x(k)
      if (a1.gt.a2) then
      j=j-1
      goto 40
      else if(a1.lt.a2) then
      goto 47
      endif
45    continue
47    if (i.le.j) then
      do 50 j2=1,iy
      a1=a(j2,i)
      a(j2,i)=a(j2,j)
50    a(j2,j)=a1
      i=i+1
      j=j-1
      endif
      if (i.le.j) goto 30
      if (i.lt.r) then
      s=s+1
      stack(s,1)=i
      stack(s,2)=r
      endif
      r=j
      if (l.lt.r) goto 20
      if (s.ne.0) goto 10
      return
      end
 
      integer function hashves(dat,r,ind,m,n,mode,nr)
c returns address of  element of ind()  containing vector dat. 
c Searching goes with a HASH technique. When mode=1 and dat was not
c there before, dat is written  into ind("hash",1..r). When mod=0
c and dat is not found, thi function returns 0.
c rank of ind is m x n,
c dat contains r numbers.
c !! same as hashvec except that rows and columns switched !!
      integer r,mode,nr,m,n
      integer ind(n,m),dat(r)
      integer*4 iaddr,prime(10),ie,iaddress,k,j,izer,ieq,i1,i
      data ie/17/,prime/433,53,761,197,1039,31,887,389,277,97/
c        
5     iaddr=0
      do 7 i=1,r
7     iaddr=iaddr+prime(i)*dat(i)
c      iaddress=mod(iabs(iaddr),m)+1
      iaddress=iabs(mod(iaddr,m))+1
 
      do 10 k=1,400
         j=ind(1,iaddress)
         izer=0
         ieq=0
         do 20 i=1,r
            i1=ind(i,iaddress)
            if (i1.ne.dat(i)) ieq=1
 20         if (ind(i,iaddress).ne.0) izer=1
 
 
         if (izer.eq.0 .or. ieq.eq.0) then
             if (izer.eq.0 .and. mode.eq.1) then
                 do 30 i=1,r
30                    ind(i,iaddress)=dat(i)
                 nr=nr+1
             endif
             if (mode.eq.0 .and.izer.eq.0) then
                 hashves=0
               else
                 hashves=iaddress
             endif
             return
         endif
         iaddress=mod(iaddress+ie-1,m)+1
10    continue
 
      write(*,60)
60    format(' hash matrix too small')
      stop
      end                              
      
      subroutine clusters(x,m1,m2,ix)
c clusters ix rows of x containing nonzero last element
c at the beginning of the matrix
c !! same as cluster except rows and columns switched !!
 
      integer m1,m2,i,ix,x(m2,m1),j1,k1,k2
      j1=1
      k1=1
      do 120 i=1,ix
130   if (x(m2,j1).ne.0 .and. j1.le.m1) then
        j1=j1+1
        goto 130
      endif
140   if (x(m2,k1).eq.0 .and. k1.le.m1) then
        k1=k1+1
        goto 140
      endif
      if (k1.gt.j1) then
        do 150 k2=1,m2
          x(k2,j1)=x(k2,k1)
150       x(k2,k1)=0
      else
        k1=k1+1      
      endif
120   continue
      return
      end
 
 
      logical function rd(un,buf,n,free,ff)
*
c reads buf(n) integers from unit un using either free format (free=true)
c or format ff (free=false). the function is set false if end of file was
c found, it is set true otherwise
*
      logical free
      integer un,n,buf(n),i
      character ff*(*)
      if (free) then
           read(un,*,end=100)(buf(i),i=1,n)
         else
           read(un,ff,end=100)(buf(i),i=1,n)
      endif
      rd=.true.
      return
100   rd=.false.           
      return
      end
 
      logical function rd1(un,buf,bufa,m,n)
c from unit un reads words delimited by spaces into bufa, decodes them
c into integers in buf, and returns the number of words in n. buf and bufa
c have a dimension of m
c !! Input line must have less than 200 characters !!
      integer m,n,un,buf(m)
      character bufa(m)*(*),a*200
      read(un,'(a)',end=100)a
      call nums(a,buf,bufa,m,n)
      rd1=.true.
      return
100   rd1=.false.
      end
 
      subroutine zero(n,x)
      integer n,x(n),i
      do i=1,n
         x(i)=0
      enddo
      end
 
      subroutine printhash(x,l,m,n)
c print columns of hash matrix, only with nonzero first row. First
c l integers are printed as characters
      integer*4 l,m,n,x(m,n),y(20),i,j
      character c*80
      equivalence(c,y)
 
      do i=1,n
         if (x(1,i).ne.0) then
            do j=1,l
               y(j)=x(j,i)
            enddo
            print '(1x,a,10i6)',c(1:4*l),(x(j,i),j=l+1,m)
         endif
      enddo
      end
 
      logical function goodid(id)
c     return false if id contains 0s only
      character id*(*),zero*16
      data zero/'0000000000000000'/
      if (id(1:1).eq.'0') then
         if (id.eq.'0') then
            goodid=.false.
         elseif (id.le.zero) then
            goodid=.false.
         endif
        else
         goodid=.true.
       endif
       end
 
       integer function upg(yob,groups,ngroups,grstat)
c assigns unknown parent groups based on yob (year of birth) and ngroups
c years stored in table groups
       integer yob,ngroups,groups(ngroups),grstat(ngroups+1)
 
       if (ngroups.eq.0 .or.
     +           ngroups.eq.1 .and. groups(1).eq.0) then
           upg=1
         else
           upg=1
           do while (upg.le.ngroups)
              if (yob.le.groups(upg)) goto 10
c           print*,'yob,groups(upg),ngroups',yob,groups(upg),ngroups
              upg=upg+1
           enddo
10         continue
        endif
 
        grstat(upg)=grstat(upg)+1
c        print*,'upg,yob',upg,yob
        end
 
 
      subroutine checkgen(k,kpar,xani,nyob,genint,m3,m4)
c prints messages if year year of births of animal and its parent are
c outside ranges given in generation intervals. Disabled if year
c of birth=0 or genint(3)=0
      integer k,kpar,m3,m4,xani(m3,m4),nyob,genint(3),i,j,l
 
      if (genint(3).ne.0) then
         i=xani(nyob,k)
         j=xani(nyob,kpar)
         if (i.gt.0 .and. j.gt.0) then
            if (i-j .gt. genint(3)) write(11,10)(xani(l,kpar),l=1,3),
     +                                       (xani(l,k),l=1,3),i-j
         endif
      endif
 
10    format(' parent ',3a4,' older than ',3a4,' by ',i4,' years')
      end
 
      subroutine readuntil(char,buf,bufa,n,m,eln)
c reads a line from input, and deletes everything from the first
c occurence of char. If line is empty, reads another line.
c if eln is positive, new line is read if the number of
c items separated by spaces is not equal eln. If eln is negative, the
c number of items should be -eln or greater
      integer n,m,buf(n),i,eln
      character bufa(n)*(*),a*80,char
 
10    read(*,'(a)',end=100)a
      i=index(a,char)
      if (i.eq.1) goto 10
      if (i.eq.0) i=80
      call nums(a(1:i),buf,bufa,n,m)
      if (eln.gt.0 .and. eln.ne.m) then
         print '(i4,'' components expected, repeat line'')',eln
         goto 10
      endif
      if (eln.le.0 .and. -eln.gt.m) then
         print '('' at least'',i4,'' components expected''
     +                ,'', repeat line'')',-eln
         goto 10
      endif
      return
100   print*,' end of input'
      end
 
 
      subroutine assignchar(a,b)
c assigns a  character*4 variable b to a.
      character b*4
      integer*4  a
      read(b,'(a4)') a
      end

      subroutine convert(bufa,bufb,n)
c places alphanumeric entries of bufa(n) into integer columns of bufb(3,n)
c written for compatibility with some compilers, including WATCOM and DEC
      integer n,bufb(4,n),b(4),i,j
      character bufa(n)*(*), a*16
      equivalence (a,b)
c
      do i=1,n
         a=bufa(i)
         do j=1,4
            bufb(j,i)=b(j)
         enddo
      enddo
      end