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