diff -r 000000000000 -r 5fda18b64dcb src/realizer_pro.f90
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/realizer_pro.f90	Wed Sep 26 17:50:53 2007 -0400
@@ -0,0 +1,942 @@
+program realizer
+!       program MidPtFM2d
+!
+!   fractal generator adapted from algorithm MidPointFM2D in the 
+!     Science of Fractal Images, M.F. Barnsley, R. L. Devaney
+!     B. B. Mandelbrot, H.-O. Peitgen, D. Saupe, and R. F. Voss
+!     Springer-Verlag
+!
+!  this program modified to combine multiple fractal topographies
+!     then normalize them relative to each other,
+!     and perform multiple horizontal slices
+!     to obtain requested p values for each habitat
+!     in a new synthetic realization of the landscape
+!
+!     William W. Hargrove
+!     Paul M Schwartz
+!     5/7/96
+!     (423) 241-2748
+!
+      use realizer_mod
+      implicit none
+
+      integer(i4) :: P(maxN*maxN,0:maxsp+3)
+      real(r4)    :: X(0:maxN,0:maxN)
+      real(r4)    :: prob(0:maxsp)
+      real(r4)    :: ssum
+      real(r4)    :: ztemp
+      integer(i4) :: igreatest
+      integer(i4) :: hslice(0:maxsp)
+      real(r4)    :: elev(maxN*maxN)
+      real(r4)    :: H(maxN,maxsp)
+      integer(i4) :: II(0:maxN,0:maxN)
+      integer(i4) :: kcnt(0:maxsp)
+      integer(i4) :: ties(maxsp)
+      integer(i4) :: ifcnt
+      integer(i4) :: iend
+      integer(i4) :: istart
+      integer(i4) :: icnt
+      real(r4)    :: sigma
+      real(r4)    :: xmin, xmax, xmean, xss
+      real(r4)    :: range1
+      real(r4)    :: zmean
+      real(r4)    :: pfract
+      real(r4)    :: xlow
+      integer(i4) :: maxlevel
+      integer(i4) :: numpainted
+      integer(i4) :: numtied
+      integer(i4) :: numblank
+      integer(i4) :: ihighest
+      integer(i4) :: ibiggest
+      integer(i4) :: idummy
+      integer(i4) :: idtype
+      integer(i4) :: idumm
+      integer(i4) :: isum, jsum
+      integer(i4) :: ipixid
+      integer(i4) :: kount
+      integer(i4) :: i, j, k, l, m ! loop indices
+      integer(i4) :: iv, iz
+      integer(i4) :: N
+      integer :: useed                 ! user seed value
+      integer :: rsize                 ! length of array of seeds
+      integer :: ierr                  ! error
+      integer, allocatable :: iseed(:) ! array of random number seeds
+      integer(i4) :: idispvec(0:maxN*maxN)
+      character (len=1)  :: ans
+      character (len=28) :: maptitle
+      character (len=60) :: mapfile
+      character (len=10) :: mapcolor
+      character (len=40) :: fmt  ! For dynamic format strings in gfortran
+      logical :: addition, wrap, gradient, slice, invert
+      logical :: grassout, xpmout, probout, vizout
+      logical :: constraint(0:maxsp)
+!
+!  write execution script for next run
+      open (9, file='input.scr',status='unknown')
+!
+!
+! enter random number seed
+      print 196
+196   format (t10,'Do you want to enter a random number seed <Y>', &
+              /t12,' or use the system clock? [clock]')
+      read (5,97) ans
+97    format (a1)
+      if (ans .eq. 'Y' .or. ans .eq. 'y') then
+         write (9,197) ans
+197      format (a1,t20,"(enter seed or clock?)")
+         print 100
+100      format (t10,'Enter random number seed ')
+         read *, useed
+         write (9,101) useed
+101      format (i12,t20,"(random number seed)")
+      else
+21       useed = int(time8(),kind(useed))
+         if (useed .eq. 0) then
+            print *, "useed was zero."
+            go to 21
+         end if
+         print 298
+298      format (t15,'Random number seed from system clock.'/)
+         write (9,29) ans, useed
+29       format (a1,t20,"(random seed from clock = ",i12,")")
+      end if
+      ! Set the random number seed
+      call random_seed(size=rsize)
+      allocate(iseed(rsize),stat=ierr)
+      if (ierr /= 0) then
+         print *,'Error allocating iseed(rsize)'
+         stop
+      end if
+      iseed(1:rsize) = useed
+      call random_seed(put=iseed)
+!
+      print 1017, useed
+1017  format (t15,'Random number seed = ',i12/)
+!
+!
+      print 512
+512   format (t10,'Do you want to visualize', &
+              /t12, ' the synthetic landscape? <Y,N>')
+      read (5,522) ans
+522   format (a1)
+      write (9,57) ans
+57    format (a1,t20,"(visualize the landscape?)")
+      if (ans .eq. 'y' .or. ans .eq. 'Y') then
+         vizout = .true.
+      else
+         vizout = .false.
+      end if
+!
+1     print 225
+225   format (t10,'How many levels? (map is 2**levels on a side)')
+      read *, maxlevel
+!
+      N = 2**maxlevel
+      irowcols = N
+      numcells = (N+1)*(N+1)
+      if (N .gt. maxN) goto 1
+      write (9,102) maxlevel
+102   format (i3,t20,"(Number of levels)")
+!
+      print 114
+114   format (t10,'Highest data category (excluding no-data)', &
+              /t12,' in the realized synthetic landscape?')
+      read *, ihabs
+      write (9,199) ihabs
+199   format (i3,t20,"(highest data category)")
+!
+      print 232
+232   format (t10,'Will you supply external probability maps', &
+              /t12,' for any of these categories?')
+      read (5,111) ans
+111   format (a1)
+      write (9,161) ans
+161   format (a1,t20,"(supply external prob maps?)")
+!
+      do j = 0, ihabs
+         constraint(j) = .false.
+      end do
+!
+      isupply = 0
+      if (ans .eq. 'Y' .or. ans .eq. 'y') then
+!
+         print 134
+134      format (t10,'For how many categories in the synthetic', &
+                 ' landscape',/t12,' will you supply external probability', &
+                 '  maps?')
+         read *, isupply
+         write (9,189) isupply
+189      format (i3,t20,"(number of external prob maps)")
+!
+         do i = 1, isupply
+            print 216
+216         format (t10,'Constrain which map category?')
+            read *, k
+            write (9,799) k
+799         format (i3,t20,"(constrain which category?)")
+!           set constraint flag for this category
+            constraint(k) = .true.
+!
+            print 220
+220         format (t10,'Enter name of input map file (integer*4, ', &
+                    /t12,' delimited by single spaces): ')
+            read (5,230) mapfile
+230         format (a60)
+            print 2030, mapfile
+2030        format (t15,'Input map name: ',a60/)
+            write (9,231) mapfile
+231         format (a60)
+!
+            print 144
+144         format (t10, 'Invert probabilities (i.e., DEM)? <Y,N>')
+            read (5,321) ans
+321         format (a1)
+            write (9,327) ans
+327         format (a1,t20,"(invert probabilities?)")
+            if (ans .eq. 'y' .or. ans .eq. 'Y') then
+               invert = .true.
+            else
+               invert = .false.
+            end if
+!
+!           load constraint maps in proper column of P array
+!           invert probabilities if necessary, i.e., a DEM
+            Call InputMap(mapfile,P,N,k,invert)
+         end do
+!
+      end if
+!
+!
+      do iz = 0, ihabs
+         if (.not. constraint(iz)) then
+            do i = 1, maxlevel
+118            print 159, i, iz
+159            format (t10,'Input value of H for level ',i3,', category ',i3)
+               read *, H(i,iz)
+               print *, H(i,iz)
+               write (9,155) H(i,iz), i, iz
+155            format (f4.2,t20,"(value of H for level ",i2,&
+                       ", category", i3,")")
+               if (H(i,iz) .lt. 0.0 .or. H(i,iz) .gt. 1.0) goto 118
+            end do
+         end if
+664   end do
+!
+!
+      ssum = 0.0
+      do i = 1, ihabs
+115      print 116, i
+116      format (t10,'Probability for habitat ',i3)
+         read *, prob(i)
+         if (prob(i) .lt. 0.0) goto 115
+         if (prob(i) .gt. 1.0) goto 115
+         ssum = ssum + prob(i)
+         if (ssum .gt. 1.0000) then
+            print *, 'Total probability exceeds 1.0 - Re-enter'
+            ssum = ssum - prob(i)
+            goto 115
+         end if
+         write (9,151) prob(i), i
+151      format (f6.5,t20,"(probability of habitat type", i2,")")
+      end do
+!
+!     find greatest non-constraint p to make background
+      ztemp = 0._r4
+      igreatest = 0
+      do i = 0, ihabs
+         if (.not. constraint(i)) then
+            if (prob(i) .gt. ztemp) then
+               ztemp = prob(i)
+               igreatest = i
+               !print *, 'at i of ',i ,'new greatest is ', igreatest
+            end if
+         end if
+662   end do
+!
+!     zero category by difference
+      prob(0) = 1.0 - ssum
+      if (prob(0) .lt. 0.00001) prob(0) = 0.0
+      print *, 'prob(0) is  ', prob(0)
+!
+!      print *, 'igreatest is ', igreatest
+      if (prob(0) .gt. prob(igreatest)) igreatest = 0
+      print *, 'Cat ', igreatest,' will be calculated by absence'
+!
+!      The Fractal Realizer assumes that the user-requested category
+! with the greatest coverage on the landscape (the highest requested p)
+! is the background or matrix category.  In order to give the other
+! categories the greatest freedom for fractal tuning and to minimize
+! contention ties, this majority background color is computed by absence,
+! and all other categories are reclassed.  Furthermore, categories which
+! are requested but for which the requested p = 0.0 are not actually
+! computed.
+!
+!      calculate reclassed categories
+!      ihabs = 0
+!      do i = 0, ihabsmax
+!         if (i .eq. igreatest) then
+!            ireclass(0) = i  
+!            print *, 'Cat ', i,' reassigned to background'
+!            goto 215
+!         end if
+!         if (prob(i) .eq. 0.0) then
+!            print *, 'Cat ', i,' requested but not used - skipping'
+!            goto 215
+!         end if
+!         ihabs = ihabs + 1
+!         ireclass(ihabs) = i
+!         print *, 'cat ', i,' will now be cat ', ihabs
+!
+!215   end do
+!
+!
+!      do i = 0, ihabs
+!         print 2002, ireclass(i), i, prob(ireclass(i))
+!2002     format ('Original category ', i2, ' re-mapped to ', &
+!                 i2, ' with requested probability of ', f6.4)
+!      end do
+!
+      print 212
+212   format (t10,'Generate GRASS r.in.ascii file? <Y,N>')
+      read (5,222) ans
+222   format (a1)
+      write (9,27) ans
+27    format (a1,t20,"(generate GRASS input file?)")
+      if (ans .eq. 'y' .or. ans .eq. 'Y') then
+         grassout = .true.
+      else
+         grassout = .false.
+      end if
+!
+      print 122
+122   format (t10,'Generate xpm file of synthetic landscape? <Y,N>')
+      read (5,123) ans
+123   format (a1)
+      write (9,124) ans
+124   format (a1,t20,"(xpm file of realized landscape?)")
+      if (ans .eq. 'y' .or. ans .eq. 'Y') then
+         xpmout = .true.
+      else
+         xpmout = .false.
+      end if
+!
+      print 322
+322   format (t10,'Generate xpm file of fractal probability' & 
+              ' landscapes? <Y,N>')
+      read (5,323) ans
+323   format (a1)
+      write (9,324) ans
+324   format (a1,t20,"(xpm file of probability landscapes?)")
+      if (ans .eq. 'y' .or. ans .eq. 'Y') then
+         probout = .true.
+      else
+         probout = .false.
+      end if
+!
+! *****************************************************************
+!  fill P array with a separate fractal probability map
+!    for each desired habitat in the landscape
+      do iz = 0, ihabs
+         if (constraint(iz)) goto 22
+!
+!        don't bother if igreatest or no cells requested of this category
+         if (prob(iz) .eq. 0.0 .or. iz .eq. igreatest) then
+            do i = 1, numcells
+               P(i,iz) = 0
+            end do
+            goto 22
+         end if
+!
+!         print 102
+!102      format (t10,'Input the value for sigma')
+!         read *, sigma
+         sigma = 1.
+         print 103
+103      format (t10,'Input the threshold for detection (xlow) ')
+         read (5,119) xlow
+119      format (g12.6)
+         write (9,125) xlow
+125      format (f12.6,t20,"(threshold for detection (xlow))")
+         if (xlow .eq. 0) xlow = -1e30
+         print *, xlow
+         print 104
+104      format (t10,'Wrap? <Y,N>')
+         read (5,145) ans
+145      format (a1)
+         write (9,198) ans
+198      format (a1,t20,"(wrap?)")
+         if (ans .eq. 'y' .or. ans .eq. 'Y') then
+            wrap = .true.
+         else
+            wrap = .false.
+         end if
+!
+!  no gradient possible with wrap
+         if (.not. wrap) then
+            print 106
+106         format (t10,'Gradient? <Y,N>')
+            read (5,107) ans
+107         format (a1)
+            write (9,148) ans
+148         format (a1,t20,"(gradient?)")
+            if (ans .eq. 'y' .or. ans .eq. 'Y') then
+               gradient = .true.
+!
+               print 1850
+1850           format (t10,'Gradient: Initial values for four map corners')
+               print 1851
+1851           format (t10,'Northeast? ')
+               read *, northeast
+               print 1852, northeast
+1852           format (5x,g12.6)
+               write (9,135) northeast
+135            format (f12.6,t20,"(Gradient: northeast)")
+!
+               print 1853
+1853           format (t10,'Southeast? ')
+               read *, southeast
+               print 1852, southeast
+               write (9,136) southeast
+136            format (f12.6,t20,"(Gradient: southeast)")
+!
+               print 1854
+1854           format (t10,'Southwest? ')
+               read *, southwest
+               print 1852, southwest
+               write (9,137) southwest
+137            format (f12.6,t20,"(Gradient: southwest)")
+!
+               print 1855
+1855           format (t10,'Northwest? ')
+               read *, northwest
+               print 1852, northwest
+               write (9,138) northwest
+138            format (f12.6,t20,"(Gradient: northwest)")
+!
+            else
+               gradient = .false.
+            end if
+!
+         end if !wrap = .false.
+!
+!         print 111
+!111      format (t10,'Slice into discrete habitats? <Y,N>')
+!         read (5,110) ans
+!110      format (a1)
+!         if (ans .eq. 'y' .or. ans .eq. 'Y') then
+!            slice = .true.
+!         else
+!            slice = .false.
+!         end if
+!
+         slice = .false.
+!
+!         print 104
+! 104     format (t10,'Addition <Y,N>?')
+!         read (5,105) ans
+! 105     format (a1)
+!         if (ans .eq. 'y' .or. ans .eq. 'Y') then 
+!            addition = .true.
+!         else
+!            addition = .false.
+!         end if
+         addition = .false.
+!
+!  call 2d multi-fractal generator
+         print *,"Calling fractal generator for category ", iz
+         call mpfm2d(X, maxlevel, sigma, H, addition, wrap, gradient, iz)
+!
+         print *, "fractal generated"
+         xmin = 1e30
+         xmax = -1e30
+         xmean = 0.0
+         xss = 0.0
+         icnt = 0
+!
+!  obtain xmin and xmax for range normalization
+!   calculate initial raw statistics
+!   while you're at it
+         do i = 0,N
+            do j = 0,N
+               if (X(i,j) .lt. xmin) xmin = X(i,j)
+               if (X(i,j) .gt. xmax) xmax = X(i,j)
+!               xmean = xmean + X(i,j)
+!               xss = xss + X(i,j)*X(i,j)
+!              icnt = icnt+1
+            end do
+         end do
+!         xmean = xmean / icnt
+!         print 200, xmean, xss, xmin, xmax
+!200      format (/t10,'Hfract Summary'/ &
+!                  t10,'Mean =    ',g12.6 &
+!                 /t10,'Sum Sq. = ',g12.6 &
+!                 /t10,'Minimum = ',g12.6 &
+!                 /t10,'Maximum = ',g12.6)
+!
+!  rescale and normalize fractal map between 0 and 32767, inclusive
+!
+         range1 = xmax - xmin
+         do i = 0,N
+            do j = 0,N
+               iv = (i*(N+1))+j+1
+               P(iv,iz) = ((X(i,j) - xmin)/range1) * 32767
+               P(iv,ihabs+1) = i
+               P(iv,ihabs+2) = j
+!              write(10,*) P(iv,iz)
+            end do
+         end do
+         print*,'done normalizing now'
+!
+!c        write GRASS ascii map
+!         write(10,791), N
+!791      format ('north: ',i4)
+!         write(10,792)
+!792      format ('south: 0')
+!         write(10,793), N
+!793      format ('east: ',i4)
+!         write(10,794)
+!794      format ('west: 0')
+!         write(10,795), N+1
+!795      format ('rows: ',i4)
+!         write(10,796), N+1
+!796      format ('cols: ',i4)
+!         do i = 0,N
+!            write(10,*) (P((i*(N+1))+j+1,iz), j=0,N)
+!         end do
+!
+!  end of iz loop over all landscape categories in synthetic map
+22       continue
+      end do
+! **********************************************************************
+!
+!  ihabs fractal maps now installed in P matrix:
+!
+!                 ihabs cols + 3
+!		|---------------|
+!		|		|
+!		|		|
+!		|  empty cells  |
+!  mtresolved =>|		|
+!		|		|
+!   mtsplit   =>|---------------|
+!		|		|
+!		|		|
+!		|		|
+!		|		|
+!		|    singly-    |
+!		|    painted    |
+!		|		|
+!		|		|
+!		|		|
+!		|		|
+!		|		|
+!  itiedsplit =>|---------------|
+!		|		|
+!		|		|
+!		|     ties      |
+!           m =>|		|
+!		|		|
+!    numcells =>|_______________|
+!
+!  ihabs+1 column carries x-coord of cell
+!  ihabs+2 column carries y-coord of cell
+!  ihabs+3 column carries number of ties for this cell
+! 
+!
+!  initialize II array with background value for re-use as output map
+      II(:,:) = 32767
+!
+!      set paint thresholds
+      print *,"Calculating paint thresholds"
+      do i = 0,ihabs ! over all map categories
+!
+!c        don't bother if no cells requested of this cat
+         if (prob(i) .ne. 0.0) then
+!
+!            istart = 0
+!            print *, "numcells is ", numcells
+!            print *, "istart is ", istart
+!            print *, "i is ", i
+!            Call HeapSort(P,numcells,istart,i) ! sort P matrix by ith column
+!
+!
+!            next category to be painted goes to single vector for sorting
+            idispvec(0) = 0
+            do l = 1, numcells
+               idispvec(l) = P(l,i)
+            end do
+!
+!            call special vector version of Heapsort
+            istart = 1
+            Call HS2(idispvec,numcells+1,istart)
+!
+!c           diagnostic print
+!            do k=16636,16646
+!               print *, i, k, idispvec(k)
+!            end do
+!
+!            set paint threshold for this map category
+            hslice(i) = idispvec(numcells - (int((prob(i) * numcells)+.499)))
+!           for absolute constraint maps
+            if (hslice(i) .eq. 32767 .and. constraint(i)) hslice(i) = 32766
+            print *, "trying to assign ", int((prob(i) * numcells)+.499), &
+                     " cells for category", i
+            jsum = jsum + int((prob(i) * numcells)+.499)
+         end if
+33    end do ! over all map categories
+!
+      print*, 'trying to assign a total of', jsum,' cells of', numcells, &
+              ' total cells'
+!c
+!c     set paint threshold for this map category
+!      hslice(i) = P(numcells - (int((prob(i) * numcells)+.499)),i)
+!      print *, "trying to assign ", int((prob(i) * numcells)+.499), &
+!               " cells for category", i
+!      end do ! over all map categories
+!
+!      dump hslice array for checking
+      print *, 'horizontal slices at these elevations:'
+      do i = 0,ihabs
+         print *, i, hslice(i)
+      end do
+!
+!     write XPM file of fractal probability terrain aspect
+      if (probout) then
+         do iz = 0, ihabs
+            if (.not. (constraint(iz) .or. prob(iz) .eq. 0.0 .or. &
+                iz .eq. igreatest)) then
+               Call XPMRelief(P,iz,hslice(iz))
+            end if
+335      end do
+      end if
+!
+!
+      numpainted = 0
+      numtied = 0
+      numblank = 0
+!
+      print *,"Beginning tie resolution process"
+      do l = 1, numcells ! start at top of P matrix
+         idummy = 0
+         do j = 0, ihabs ! over all map categories
+!           probably don't need to do this if prob(j) .eq. 0.0
+            if (.not. (prob(j) .eq. 0.0 .or. j .eq. igreatest)) then
+!
+!              if this cell is painted this category
+               if ( P(l,j) .gt. hslice(j) ) then
+!               if ( P(l,j) .ge. hslice(j) ) then
+!
+!                 count the multiple assignments
+                  P(l,ihabs+3) = P(l,ihabs+3) + 1
+!
+!                 find the winner category for this cell
+!                 winner category has the highest probability surface
+                  if ( P(l,j) .gt. idummy ) then
+                     ihighest = j
+                     idummy = P(l,ihighest)
+                  end if
+!
+               end if ! this cell painted this category
+!
+            end if
+         end do ! over all map categories
+!         print *, "highest for cell ", l," was category", ihighest
+!
+!        keep count of each type of cell for displacement into
+!         sorted array later
+!        if singly-painted or currently tied, paint outmap with ihighest
+         if ( P(l,ihabs+3) .eq. 1) then
+!           paint this cell with the winning category
+!            print *, "painting uncontested cell", P(l,ihabs+1), &
+!                     P(l,ihabs+2)," with cat", ihighest
+!
+            numpainted = numpainted + 1
+            II( P(l,ihabs+1),P(l,ihabs+2) ) = ihighest
+         else if ( P(l,ihabs+3) .eq. 0) then
+            numblank = numblank + 1
+         else
+            numtied = numtied + 1
+!            print *, "painting tied cell", P(l,ihabs+1),P(l,ihabs+2), &
+!                     " with winning cat", ihighest
+            II( P(l,ihabs+1),P(l,ihabs+2) ) = ihighest
+         end if
+!
+      end do
+!
+!     at this point, II array has all winners and uncontested
+!     cells painted
+!
+!     sort P array into 3 sections by number of assignments
+!     ascending sort results in empty, singly-painted, and
+!     ties sections
+      print *,"Preparing to sort ties for resolution"
+!      
+!     sort by ties column
+      istart = 1
+      Call HeapSort(P,numcells,istart,ihabs+3)
+!      print *, "got back from HeapSort"
+!
+!     calculate pointers to section breaks
+      mtsplit = numblank
+      itiedsplit = numblank + numpainted
+!
+      mtresolved = mtsplit
+!
+      print *, "numblank is ", numblank
+      print *, "numtied is ", numtied
+      print *, "itiedsplit is ", itiedsplit
+      print *, "numpainted is ", numpainted
+      print *, "mtresolved and mtsplit are ", mtresolved
+
+      if (numtied .gt. numblank) print*, 'WARNING: More tied cells', &
+         ' than leftover blank cells to use for tie resolution!'
+!
+!     sum multiwayness of ties and compare to number of empty cells
+      do m = itiedsplit+1, numcells
+         isum = isum + ( P(m,ihabs+3) - 1 )
+      end do
+      print *, "sum of multiwayness of ties is ", isum
+      if (isum .gt. numblank) print *, 'Warning: cells are tied more', &
+         ' ways than there are blank cells to assign!'
+!
+!      call dump_map(P)
+!
+!  sort ties into priority order for resolution
+!   by calculating and sorting on summed elevation 
+!   for all tied categories
+!
+      do m = itiedsplit+1, numcells ! over all tied
+         kount = 0
+         idummy = 0
+         do k = 0, ihabs                ! over all categories
+!
+!           probably not necessary if prob(k) .eq. 0
+            if (.not. (prob(k) .eq. 0.0 .or. k .eq. igreatest)) then
+!           if this cell is painted this category
+               if ( P(m,k) .gt. hslice(k) ) then
+!               if ( P(m,k) .ge. hslice(k) ) then
+                  kount = kount + 1          ! number of ties for this cell
+                  ties(kount) = k            ! keep categories of ties
+!
+               end if
+            end if
+933      end do
+!
+         zmean = 0.0
+         do k = 1, kount
+            zmean = zmean + P(m,ties(k))
+         end do
+!         P(m,ihabs+3) = zmean/kount
+!        try weighting priority by simple sum, not average
+!        therefore, max sum = (ihabs-1) * 32767
+!        rescale as proportion of 32767 for integer*4 array cell
+         P(m,ihabs+3) = zmean/(ihabs-1)
+!
+      end do
+!
+!     sort ties by mean priority
+      Call HeapSort(P,numtied,itiedsplit+1,ihabs+3)
+!
+!     save min and max tie priorities for gray levels in XPMTies
+      maxgray = P(numcells,ihabs+3)
+      print*, "maxgray is ", maxgray
+!
+!     generate xpm picture of tied cells
+!      colored by priority gray shades
+!      Call XPMTies(P)
+!
+!
+!     start at bottom of list (highest priority ties)
+      do m = numcells, itiedsplit+1, -1 ! over all tied
+         kount = 0
+         idummy = 0
+         do k = 0, ihabs ! over all categories
+!
+!           probably not necessary if prob(k) .eq. 0
+            if (.not. (prob(k) .eq. 0.0 .or. k .eq. igreatest)) then
+!              if this cell is painted this category
+               if ( P(m,k) .gt. hslice(k) ) then
+!               if ( P(m,k) .ge. hslice(k) ) then
+                  kount = kount + 1 ! number of ties for this  cell
+                  ties(kount) = k ! keep categories of ties
+!                 save most likely category for this cell
+                  if ( P(m,k) .gt. idummy ) then
+                     ihighest = k
+                     idummy = P(m,ihighest)
+                  end if
+               end if
+            end if
+656      end do
+!
+!
+!        resolve all ties -
+!        over all tied categories that are not the winner
+!         (which has already been painted)
+!         we need an empty cell -
+!         search the empty cell list by this cat to find next
+!         most likely available empty cell
+!
+!        empty cell list will be searched once for each tie category
+!         except the winning category
+!
+!         print *, "deciding among ", kount, " ties for cell ", m
+         do k = 1, kount ! over all ties for this cell
+            if ( ties(k) .ne. ihighest ) then ! don't re-do the winner
+!               print *, ties(k), " was a tied category for cell ", m
+!
+               ibiggest = 0
+               idummy = 0
+               do i = 1,mtresolved
+                  if (P(i,ties(k)) .gt. idummy) then
+                     ibiggest = i
+                     idummy = P(i,ties(k))
+                  end if
+               end do
+!
+               Call Interchange(P,ibiggest,mtresolved)
+!
+!               istart = 1
+!               Call HeapSort(P,mtresolved,istart,ties(k))
+!
+!              now mtresolved points to highest
+!              next available empty cell of category ties(k)
+!              paint it!
+!               print *, "painting empty cell", mtresolved," at coords ", &
+!                  P(mtresolved,ihabs+1), P(mtresolved,ihabs+2)," with prob ", &
+!                  P(mtresolved,ties(k))," with cat", ties(k), &
+!                  " for tie at cell", m
+               II( P(mtresolved,ihabs+1),P(mtresolved,ihabs+2) ) = ties(k)
+!              burn one empty
+               mtresolved = mtresolved - 1
+            end if
+!
+992      end do ! over all ties for this cell
+!         print *, "resolving ties at cell ", m
+!
+      end do ! over all tied cells
+!
+!**********************************************************************
+!
+!
+!      Call dump_map(P)
+!
+!     generate xpm picture of tied cells
+!      colored by priority gray shades
+!      and used blank cells colored in yellow
+      Call XPMTies(P)
+!
+!
+!     set background category and
+!     perform summary of realized map
+      do i = 0,N
+         do j = 0,N
+            if (II(i,j) .eq. 32767) II(i,j) = igreatest
+            kcnt(II(i,j)) = kcnt(II(i,j)) + 1
+!           write(12,*) II(i,j)
+         end do
+      end do
+!
+!     if requested,
+!      write GRASS ascii map
+      if (grassout) then
+         write(10,91), N
+91       format ('north: ',i4)
+         write(10,92)
+92       format ('south: 0')
+         write(10,93), N
+93       format ('east: ',i4)
+         write(10,94)
+94       format ('west: 0')
+         write(10,95), N+1
+95       format ('rows: ',i4)
+         write(10,96), N+1
+96       format ('cols: ',i4)
+         write(fmt,'("(",I0,"(i5,1x))")') j
+         do i = 0,N
+!            write(10,47) (II(i,j), j=0,N)
+            write(10,fmt) (II(j,i), j=0,N)
+!            write(10,47) (II(j,i), j=0,N)
+!47          format (<j>(i5,1x))
+         end do
+      end if
+!
+!     write final map to xpm file
+      if (xpmout) Call XPMMap(N,II)
+!
+      print 299
+299   format (/t10,'Actual Map Summary:'/t12,'Habitat',5x, &
+              'P',3x,'Real Cells',3x,'Real P')
+!      pfract = 0.0
+      do i = 0, ihabs
+         pfract = real(kcnt(i),kind(pfract)) / real(numcells,kind(pfract))
+         print 300, i, prob(i), kcnt(i), pfract
+300      format (t10, i6, 4x, f6.5, 2x, i6, 2x, f6.5)
+      end do
+!
+!**********************************************************************
+!
+!     visualize the synthetic landscape
+      if (vizout) then
+!
+!        print final map to screen
+         do i=0,N
+            do j=0,N
+               idispvec(j) = II(i,j)
+            end do
+!            print 199, (idispvec(j), j=0,N)
+!199         format (<N+1>(i1, 1x))
+         end do
+!
+!        print current state of P array
+!         Call print(P)
+!
+!        visualize II array
+         mapcolor = 'bobs.cmap\0'
+         maptitle = 'Fractal Landscape Realizer \0'
+         idtype = 4
+         idumm = 0
+         ipixid = -1
+!
+!   build display vector
+!  make small maps more visible by repeating cells and lines
+!         irepeat = maxN / N
+!         print *, irepeat
+         k = 0
+         do i = 0,N
+!            kstart = k
+            do j=0,N
+!               do kk = 1, irepeat
+!                  idispvec(k)=II(i,j)
+                  idispvec(k)=II(j,i)
+                  k = k+1
+!               end do ! repeating cells in this row
+            end do ! all cells in this row
+!            kend = k-1
+!c           repeat lines
+!            if (irepeat .ge. 2) then
+!               do jj = 1,irepeat-1
+!                  do kk = kstart, kend
+!                     idispvec(k) = idispvec(kk)
+!                     k = k+1
+!                  end do ! dup as many times as necessary
+!               end do ! over all repeated lines
+!            end if
+!
+         end do
+!
+!c  adjust rows and columns by multiplier
+!         isize = (N+1) * irepeat
+!         jsize = (N+1) * irepeat
+!         print *, isize, jsize
+!         ierr=DrawPixmap(ipixid, isize, jsize, idtype, idispvec,
+!         ierr=DrawPixmap(ipixid,N+1,N+1,idtype,idispvec,idumm,mapcolor, &
+!                         maptitle)
+!         print *,'ierr is', ierr
+!         print *, 'done visualizing'
+!         call sleep(500)
+!
+      end if ! if visualizing landscape
+!
+      stop
+!      end program MidPtFM2d
+end program realizer