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