====== Auxiliary programs for WOMBAT ======

This page collects a small number of simple, stand-alone programs to be used in conjunction with WOMBAT.
Change the output statements as you please. \\ NB: You need a FORTRAN 95 compiler to be able to generate `executables' from the source code given here; e.g. `gfortran' is free and available for most common operating systems.

===== Reading binary output files =====

==== NRM inverse ====

<file fortran read_nrminv.f90>
! ====================
  program read_nrminv
! ====================                                             

! auxiliary program to WOMBAT to examine contents of binary file(s) 
! "nrmivers_.bin"    ->  inverse of numerator relationship matrix

  implicit none
  integer           :: i, ii, nan, nnrm, jj
  real(8)           :: detl, xx
  character(len=14) :: nrmfile ='nrminvers1.bin' ! set no. 1, 2,...  as required

! open file & read header
  open(44,file = nrmfile, status ='old', form = 'unformatted', action ='read')

  read(44) nan
  read(44) nnrm, detl
  read(44) xx

! read elements of lower triangle, row-wise & write to screen
  do i = 1, nnrm
     read(44) ii, jj, xx
     write(*,'(2i10,f12.4)') ii, jj, xx
  end do

  close(44)
  stop 'end of READ_NRMINV'
  end program read_nrminv                                        ! @C  K.M.

</file>

==== Inbreeding coefficients ====

<file fortran read_inbreed.f90>
! ====================
  program read_inbreed
! ====================                                             

! auxiliary program to WOMBAT to examine contents of binary file
! "inbreed.bin"    ->  inbreeding coefficients (in %)

  implicit none
  integer                              :: m, n, ii, i, m1, n1
  integer, dimension(:), allocatable   :: nannew
  integer, dimension(:,:), allocatable :: idwec
  real(8), dimension(:,:), allocatable :: inbreed
  character(len=11)                    :: inbrdfile = 'inbreed.bin'
  character(len=30)                    :: pedfile 
  logical                              :: lex, ids

! does the file exist?
  inquire(file = inbrdfile, exist = lex) 
  if( .not. lex ) then
      write(*,*)'Binary file "',inbrdfile,'" does not seem to exist'; stop
  end if
! open file & read header
  open(1,file = inbrdfile, status ='old', form = 'unformatted', action ='read')
  read(1) m, n
! ... allocate arrays
  allocate( nannew(m), stat = ii); if(ii /= 0 ) stop 'alloc nannew'
  rewind(1); read(1) m, n, nannew, pedfile
  allocate(inbreed(n,m), stat=ii); if(ii/=0) stop 'alloc inbreed'
! ... read inbreeding coefficients
  read(1) inbreed
! ... look for list of IDs
  read(1,iostat = ii ) m1, n1
  if( ii /= 0 ) then
      ids = .false.; go to 10
  end if
  allocate( idwec(m1,n1), stat = ii ); if( ii /= 0 ) stop 'alloc idwec'
  read(1,iostat=ii) idwec; if( ii == 0 ) ids = .true.
10 close(1)

! now ready to write out information - modify to your requirements
  write(*,'(5a)')'Information from WOMBAT binary file "',inbrdfile,'"'
  write(*,'(5a)')'  ... inbreeding for pedigree file "',trim(pedfile),'"'
  write(*,'(a,i8)')'  ... no. of genetic effects       =',m
  write(*,'(a,8i8)')'  ... no. of levels                =',nannew

! ... list of inbreeding coefficients
  if( ids ) then
     do i = 1, n
        write(*,'(i8,i12,f9.3)') i, idwec(i,1), 100.d0*inbreed(i,1)
     end do
  else                 ! no. ids (older version of wombat)
     do i = 1, n
        write(*,'(i8,f9.3)') i, 100.d0*inbreed(i,1)
     end do
  end if

  stop 'end of READ_INBREED'
  end program read_inbreed              ! @(C)  K.M. 2009
</file>