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
