module realizer_mod

! Realizer module written for the MidPtFM2d program (now called realizer) by

! Forrest Hoffman on Fri May  4 16:56:34 EDT 2007 to eliminate duplicated

! parameters and common blocks.  Parts of the code were converted to

! Fortran-90 style and intrinsics.  Compiled and tested with gfortran

! version 4.1.1 under Fedora Core 6 Linux.

   implicit none

   save

   ! Parameters

   integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real

   integer, parameter :: r4 = selected_real_kind( 6) ! 4 byte real

   integer, parameter :: nr = kind(1.0)              ! native real

   integer, parameter :: i8 = selected_int_kind(13)  ! 8 byte integer

   integer, parameter :: i4 = selected_int_kind( 6)  ! 4 byte integer

   integer, parameter :: ni = kind(1)                ! native integer

   integer(i4), parameter :: maxN=515

   integer(i4), parameter :: maxsp=10

   real(r8),    parameter :: rad_per_deg=(3.14159265358979323846/180.0)

   real(r8),    parameter :: deg_per_rad=(180.0/3.14159265358979323846)

   ! Global Variables

   integer(i4) :: numcells

   integer(i4) :: itiedsplit

   integer(i4) :: irowcols

   integer(i4) :: mtsplit

   integer(i4) :: mtresolved

   integer(i4) :: isupply

   integer(i4) :: ihabs

   integer(i4) :: maxgray

   real(r4)    :: northeast

   real(r4)    :: southeast

   real(r4)    :: southwest

   real(r4)    :: northwest

   ! Private Member Functions

   public  dump_map

   public  mpfm2d

   private StackSort ! unused?

   private sort      ! unused?

   private SwapDown

   public  Interchange

   public  HeapSort

   private SD2

   private I2

   public  HS2

   private Gauss

   private cosd

   private sind

   private atand

   private atan2d

contains

!

!------------------------------------------------------------------------------

!

   subroutine dump_map(P)

      implicit none

      integer(i4), intent(in) :: P(maxN*maxN,0:maxsp+3)

      ! local variables

      integer(i4) :: i, j ! loop indices

      integer(i4) :: temp(maxsp+3)

      character (len=40) :: fmt

!

      write(fmt,'("(i5, 1x, ",I0,"(i5, 1x), 1x, 3(i5, 1x))")') ihabs+1

      do i=1,numcells

         do j = 0,ihabs+3

           temp(j) = P(i,j)

         end do

         print fmt, i, (temp(j), j=0,ihabs+3)

      end do

!

      return

   end subroutine dump_map

!

!------------------------------------------------------------------------------

!

   real(r4) function f3(delta, x0, x1, x2, sigma)

      implicit none

      real(r4), intent(in) :: delta

      real(r4), intent(in) :: x0, x1, x2

      real(r4), intent(in) :: sigma

      f3 = (x0+x1+x2) / 3.0 + delta * Gauss(0.,sigma)

      return

   end function f3

!

!------------------------------------------------------------------------------

!

   real(r4) function f4(delta, x0, x1, x2, x3, sigma)

      implicit none

      real(r4), intent(in) :: delta

      real(r4), intent(in) :: x0, x1, x2, x3

      real(r4), intent(in) :: sigma

      f4 = (x0+x1+x2+x3) / 4.0 + delta * Gauss(0.,sigma)

      return

   end function f4

!

!------------------------------------------------------------------------------

!

   subroutine mpfm2d(X, maxlevel, sigma, H, addition, wrap, gradient, iz)

      implicit none

      real(r4), intent(inout) :: X(0:maxN,0:maxN)

      integer(i4), intent(in) :: maxlevel

      real(r4), intent(in)    :: sigma

      real(r4), intent(in)    :: H(maxN,maxsp)

      logical, intent(in)     :: addition

      logical, intent(in)     :: wrap

      logical, intent(in)     :: gradient

      integer(i4), intent(in) :: iz

      ! local variables

      integer(i4) :: N

      real(r4)    :: delta

      integer(i4) :: DD

      integer(i4) :: d

      integer(i4) :: is, ix, iy ! loop indices

!

       print *, "Generating fractal"

       N = 2**maxlevel

       delta = sigma

       X(0,0) = delta * Gauss(0.,sigma)

       X(0,N) = delta * Gauss(0.,sigma)

       X(N,0) = delta * Gauss(0.,sigma)

       X(N,N) = delta * Gauss(0.,sigma)

!

       if (gradient) then

         X(0,0) = southwest

         X(0,N) = southeast

         X(N,0) = northwest

         X(N,N) = northeast

       end if

!

       DD = N

       d = N / 2

!

       do is = 1, maxlevel

!  going from grid type I to type II

         delta = delta * 0.5**(0.5 * H(is,iz))

         do ix = d, N-d, DD

           do iy = d, N-d, DD

             X(ix,iy) = f4(delta, X(ix+d,iy+d), X(ix+d,iy-d), &

                           X(ix-d,iy+d), X(ix-d,iy-d), sigma)

           end do

         end do

!  displace other points also if needed

         if (addition) then

           do ix = 0, N, DD

             do iy = 0, N, DD

               X(ix,iy) = X(ix,iy) + delta * Gauss(0.,sigma)

             end do

           end do

         end if

!  going from grid type II to type I

         delta = delta * 0.5**(0.5*H(is,iz))

!  interpolate and offset boundary grid points

         do ix = d, N-d, DD

           X(ix,0) = f3(delta, X(ix+d,0), X(ix-d,0), X(ix,d), sigma)

           X(ix,N) = f3(delta, X(ix+d,N), X(ix-d,N), X(ix,N-d), sigma)

           X(0,ix) = f3(delta, X(0,ix+d), X(0,ix-d), X(d,ix), sigma)

           X(N,ix) = f3(delta, X(N,ix+d), X(N,ix-d), X(N-d,ix), sigma)

           if (wrap) then

             X(ix,N) = X(ix,0)

             X(N,ix) = X(0,ix)

           end if

         end do

!  interpolate and offset inter grid points

         do ix = d, N-d, DD

           do iy = DD, N-d, DD

              X(ix,iy) = f4(delta, X(ix,iy+d), X(ix,iy-d), &

                            X(ix+d,iy), X(ix-d,iy), sigma)

           end do

         end do

         do ix = DD, N-d, DD

           do iy = d, N-d, DD

             X(ix,iy) = f4(delta, X(ix,iy+d), X(ix,iy-d), &

                           X(ix+d,iy), X(ix-d,iy), sigma)

           end do

         end do

!  displace other points also if needed

         if (addition) then

           do ix = 0, N, DD

             do iy = 0, N, DD

               X(ix,iy) = X(ix,iy) + delta * Gauss(0.,sigma)

             end do

           end do

           do ix = d, N-d, DD

             do iy = d, N-d, DD

               X(ix,iy) = X(ix,iy) + delta * Gauss(0.,sigma)

             end do

           end do

         end if

!

         DD = DD / 2

         d = d / 2

!

      end do ! for all levels

      return

   end subroutine mpfm2d

!

!------------------------------------------------------------------------------

!

   real(r4) function Gauss(xm, xs)

      implicit none

      real(r4), intent(in) :: xm

      real(r4), intent(in) :: xs

      ! local variables

      real(r4)    :: total

      real(r4)    :: p

      integer(i4) :: i

!

! create a normal number by summing 12 uniform random numbers

      total = 0.0

      do i = 1,12

          call random_number(p)

          total = total + p

      end do

      Gauss = (total - 6.0) * xs + xm

!

      return

   end function Gauss

!

!------------------------------------------------------------------------------

!

   subroutine StackSort(n,arr)

      implicit none

      integer(i4), intent(in) :: n

      real(r4), intent(inout) :: arr(n)

      ! local variables

      integer(i4) :: i, j ! loop indices

      real(r4)    :: a

!

!  pick out each element in turn

      do j = 2,n

         a = arr(j)

!  look for the place to insert it

         do i = j-1,1,-1

           if (arr(i) <= a) go to 10

           arr(i+1) = arr(i)

         end do

         i = 0

!  insert the number here

10       arr(i+1) = a

      end do

      return

   end subroutine StackSort

!

!------------------------------------------------------------------------------

!

!  sort routine from Numerical Recipes

   subroutine sort (n, ra)

      implicit none

      integer(i4), intent(in) :: n

      real(r4), intent(inout) :: ra(n)

      ! local variables

      integer(i4) :: i, j, l

      integer(i4) :: ir

      real(r4)    :: rra

       l = n/2 + 1

       ir = n

!

10     continue

          if (l .gt. 1) then

             l = l-1

             rra = ra(l)

          else

             rra = ra(ir)

             ra(ir) = ra(1)

             ir = ir-1

             if (ir .eq. 1) then

                ra(1) = rra

                return

             end if

          end if

          i = l

          j = l+l

!

20     if (j .le. ir) then

           if (j .lt. ir) then

              if (ra(j) .lt. ra(j+1)) j = j+1

           end if

           if (rra .lt. ra(j)) then

               ra(i) = ra(j)

               i = j

               j = j+j

           else

               j = ir+1

           end if

           goto 20

       end if

!

       ra(i) = rra

       goto 10

!

   end subroutine sort

!

!------------------------------------------------------------------------------

!

   subroutine SwapDown(Heap, r, isize, istart, keyfield)

      implicit none

      integer(i4), intent(inout) :: Heap(maxN*maxN,0:maxsp+3)

      integer(i4), intent(inout) :: r

      integer(i4), intent(in)    :: isize

      integer(i4), intent(in)    :: istart

      integer(i4), intent(in)    :: keyfield

      ! local variables

      integer(i4) :: child

      integer(i4) :: ipos

      logical :: done

!

!      print *, 'entering SwapDown'

!      print *, 'params are ', r, isize, istart, keyfield

      ipos = istart - 1

      done = .false.

      child = 2 * r

      do while (.not. done .and. child .le. isize)

! find the largest child

         if (child .lt. isize .and. Heap(child+ipos,keyfield) .lt. &

            Heap(child+ipos+1,keyfield)) child = child + 1

! interchange node and largest child if necessary

! and move down to the next subtree

         if (Heap(r+ipos,keyfield) .lt. Heap(child+ipos,keyfield)) then

            Call Interchange(Heap,r+ipos,child+ipos)

            r = child

            child = 2 * child

         else

            done = .true.

         end if

!

      end do

!

      return

   end subroutine SwapDown

!

!------------------------------------------------------------------------------

!

   subroutine Interchange(X,m,n)

      implicit none

      integer(i4), intent(inout) :: X(maxN*maxN,0:maxsp+3)

      integer(i4), intent(in) :: m

      integer(i4), intent(in) :: n

      ! local variables

      integer(i4) :: i

      integer(i4) :: temp

!

      do i = 0, ihabs+3

         temp = X(m,i)

         X(m,i) = X(n,i)

         X(n,i) = temp

      end do

!

      return

   end subroutine Interchange

!

!------------------------------------------------------------------------------

!

   subroutine HeapSort(X,isize,istart,keyfield)

      implicit none

      integer(i4), intent(inout) :: X(maxN*maxN,0:maxsp+3)

      integer(i4), intent(in) :: isize

      integer(i4), intent(in) :: istart

      integer(i4), intent(in) :: keyfield

      ! local variables

      integer(i4) :: i ! loop indices

      integer(i4) :: r

      integer(i4) :: tempr

!

!      print *, "isize in HeapSort is ", isize

!      print *, "istart in HeapSort is ", istart

!      print *, "keyfield in HeapSort is ", keyfield

! heapify

      do r = isize/2, 1, -1

         tempr = r

         !Call SwapDown(X,r,isize,istart,keyfield)

         Call SwapDown(X,tempr,isize,istart,keyfield)

         !r = tempr

      end do

!

!      print *, "P after heapify is "

!      do i=istart,isize

!         print *, i, X(i,keyfield)

!      end do

!

! repeatedly put largest item in root at end of list,

! prune it from the tree

! and apply SwapDown to the rest of the tree

!

      do i = isize, 2, -1

         Call Interchange(X,istart,(i+istart-1))

!         r = isize

         r = 1

         Call SwapDown(X,r,i-1,istart,keyfield)

!

!         print *, "P after Swapdown is  "

!         do j=istart,isize

!           print *, j, X(j,keyfield)

!         end do

      end do

!

      return

   end subroutine HeapSort

!

!------------------------------------------------------------------------------

!

   subroutine XPMMap(N,II)

!

!  output routine for xpm X-window display tool

!

!

      implicit none

      integer(i4), intent(in) :: N

      integer(i4), intent(in) :: II(0:maxN,0:maxN)

      ! local variables

      integer(i4) :: itemp(0:maxN)

      integer(i4) :: i, j ! loop indices

      character (len=40) :: fmt

!

!

!  output map

      print *, "writing final map to xpm format file landscape.xpm"

!

      open (20,file='landscape.xpm',status='unknown')

!

      write (20, 100)

100   format ('/* XPM */',/'static char *realizer[]', &

              ' = {',/'/* cols rows ncolors chars_per_pixel', &

              ' */')

      write (20, 102) N+1, N+1

102   format ('"',i6,1x,i6,1x,'11 1 0 0",',/'/* colors */')

      write (20, 104)

104   format ('"0  c white    m white', &

              '",',/'"1  s Decid  c BurlyWood    m black', &

              '",',/'"2  s Mixed  c DarkGoldenrod3       m gray40', &

              '",',/'"3  s Everg  c MistyRose3         m gray50', &

              '",',/'"4  s Trans  c LightGoldenrod4         m gray60', &

              '",',/'"5  s Urban  c DarkGoldenrod1       m gray70', &

              '",',/'"6  s Water  c gold    m gray80', &

              '",',/'"7  s Maint  c LightGoldenrod3       m white', &

              '",',/'"8  s Cropl  c SaddleBrown    m lightgray', &

              '",',/'"9  s Lawng  c RosyBrown    m gray', &

              '",',/'"*  c chocolate     m black",',/'/* pixels */')

!

      write(fmt,'("(",I0,"(i1)$)")') N+1

      do i= 0, N

         do j = 0, N

!            itemp(j) = II(i,j)

            itemp(j) = II(j,i)

         end do

         write (20,'(1A$)') '"'

!         write (20,108) (itemp(j), j=0,N)

         write (20,fmt) (itemp(j), j=0,N)

         write (20,'(2A)') '",'

!108      format (' "',<N+1>(i1),'",')

      end do

!

!

      write (20,112)

112   format ('};')

!

      return

   end subroutine XPMMap

!

!------------------------------------------------------------------------------

!

   subroutine InputMap(mapfile,P,N,k,invert)

!

!   input map from a data file

!

      implicit none

      character (len=60), intent(in) :: mapfile

      integer(i4), intent(inout)     :: P(maxN*maxN,0:maxsp+3)

      integer(i4), intent(in)        :: N

      integer(i4), intent(in)        :: k

      logical, intent(in)            :: invert

      ! local variables

      integer(i4) :: i, j, iv ! loop indices

      real(r4)    :: X(0:maxN,0:maxN)

      real(r4)    :: xmin

      real(r4)    :: xmax

      real(r4)    :: xrange

      real(r4)    :: temp

!

      open (14, file=mapfile, status='old')

!

      xmin = 1e30

      xmax = -1e30

      do j = 0, N

         read (14,*) ( X(i,j), i = 0, N)

!         read (14,118) X(i,j)

!118      format (<(N+1)*(N+1)>i5, 1x)

         do i = 0, N

            if (X(i,j) .lt. xmin) xmin = X(i,j)

            if (X(i,j) .gt. xmax) xmax = X(i,j)

         end do

      end do

!

!      invert if necessary and scale between 0 and 32767

       xrange = xmax - xmin

       do j = 0, N

          do i = 0, N

             if (invert) then

                temp = ((-X(i,j) + xmax) / xrange) * 32767

              else

                temp = ((X(i,j) - xmin) / xrange) * 32767

              end if

!

              iv = (i*(N+1))+j+1

              P(iv,k) = nint(temp)

              P(iv,ihabs+1) = i

              P(iv,ihabs+2) = j

          end do

      end do

!

      close(14)

      rewind(14)

!

      return

   end subroutine InputMap

!

!------------------------------------------------------------------------------

!

   subroutine SD2(Heap,r,isize,istart)

      implicit none

      integer(i4), intent(inout) :: Heap(maxN*maxN)

      integer(i4), intent(inout) :: r

      integer(i4), intent(in)    :: isize

      integer(i4), intent(in)    :: istart

      ! local variables

      integer(i4) :: ipos

      integer(i4) :: child

      logical     :: done

!

!      print *, 'entering SwapDown'

!      print *, 'params are ', r, isize, istart

      ipos = istart - 1

      done = .false.

      child = 2 * r

      do while (.not. done .and. child .le. isize)

!  find the largest child

         if (child .lt. isize .and. Heap(child+ipos) .lt. &

             Heap(child+ipos+1)) child = child + 1

!  interchange node and largest child if necessary

!   and move down to the next subtree

         if (Heap(r+ipos) .lt. Heap(child+ipos)) then

            Call I2(Heap,r+ipos,child+ipos)

            r = child

            child = 2 * child

         else

            done = .true.

         end if

!

      end do

!

      return

   end subroutine SD2

!

!------------------------------------------------------------------------------

!

   subroutine I2(X,m,n)

      implicit none

      integer(i4), intent(inout) :: X(maxN*maxN)

      integer(i4), intent(in)    :: m

      integer(i4), intent(in)    :: n

      ! local variables

      integer(i4) :: temp

!

      temp = X(m)

      X(m) = X(n)

      X(n) = temp

!

      return

   end subroutine I2

!

!------------------------------------------------------------------------------

!

   subroutine HS2(X,isize,istart)

      implicit none

      integer(i4), intent(inout) :: X(maxN*maxN)

      integer(i4), intent(in)    :: isize

      integer(i4), intent(in)    :: istart

      ! local variables

      integer(i4) :: i, j ! loop indices

      integer(i4) :: r

      integer(i4) :: tempr

!

!      print *, "isize in HeapSort is ", isize

!      print *, "istart in HeapSort is ", istart

!     heapify

      do r = isize/2, 1, -1

         tempr = r

         Call SD2(X,tempr,isize,istart)

         !Call SD2(X,r,isize,istart)

         !r = tempr

      end do

!

!      print*, "P after heapify is "

!      do i=istart,isize

!         print *, i, X(i)

!      end do

!

!   repeatedly put largest item in root at end of list,

!   prune it from the tree

!   and apply SwapDown to the rest of the tree

!

      do i = isize, 2, -1

         Call I2(X,istart,(i+istart-1))

!         r = isize

         r = 1

         Call SD2(X,r,i-1,istart)

!

!         print*, "P after Swapdown is  "

!         do j=istart,isize

!           print *, j, X(j)

!         end do

      end do

!

      return

   end subroutine HS2

!

!------------------------------------------------------------------------------

!

   subroutine XPMTies(P)

!

!  output routine for xpm X-window display tool

!

!  write xpm format map with grayshades

!

      implicit none

      integer(i4), intent(in) :: P(maxN*maxN,0:maxsp+3)

      ! local variables

      integer(i4) :: itemp(0:maxN,0:maxN)

      integer(i4) :: i, j, m ! loop indices

      real(r4)    :: scalelog

      character (len=40) :: fmt

!

      open (40,file='ties.xpm',status='unknown')

!

      print 99

99    format ('XPM map file ties.xpm being written')

      write (40, 100)

100   format ('/* XPM */',/'static char *tiesmap[]', &

              ' = {',/'/* cols rows ncolors chars_per_pixel', &

              ' */')

      write (40, 102) irowcols+1, irowcols+1

102   format (' "',i6,1x,i6,1x,'108 3 0 0",',/'/* colors */')

      write (40, 104)

104   format ('"  0 s 0.0prob c white g4 white g white m white', &

              '",',/'"  1 s 0.01prob c gray99 g4 gray99 g gray99 m gray99', &

              '",',/'"  2 s 0.02prob c gray98 g4 gray98 g gray98 m gray98', &

              '",',/'"  3 s 0.03prob c gray97 g4 gray97 g gray97 m gray97', &

              '",',/'"  4 s 0.04prob c gray96 g4 gray96 g gray96 m gray96', &

              '",',/'"  5 s 0.05prob c gray95 g4 gray95 g gray95 m gray95', &

              '",',/'"  6 s 0.06prob c gray94 g4 gray94 g gray94 m gray94', &

              '",',/'"  7 s 0.07prob c gray93 g4 gray93 g gray93 m gray93', &

              '",',/'"  8 s 0.08prob c gray92 g4 gray92 g gray92 m gray92', &

              '",',/'"  9 s 0.09prob c gray91 g4 gray91 g gray91 m gray91', &

              '",',/'" 10 s 0.10prob c gray90 g4 gray90 g gray90 m gray90', &

              '",',/'" 11 s 0.11prob c gray89 g4 gray89 g gray89 m gray89', &

              '",',/'" 12 s 0.12prob c gray88 g4 gray88 g gray88 m gray88', &

              '",',/'" 13 s 0.13prob c gray87 g4 gray87 g gray87 m gray87', &

              '",',/'" 14 s 0.14prob c gray86 g4 gray86 g gray86 m gray86', &

              '",',/'" 15 s 0.15prob c gray85 g4 gray85 g gray85 m gray85', &

              '",',/'" 16 s 0.16prob c gray84 g4 gray84 g gray84 m gray84', &

              '",',/'" 17 s 0.17prob c gray83 g4 gray83 g gray83 m gray83', &

              '",',/'" 18 s 0.18prob c gray82 g4 gray82 g gray82 m gray82', &

              '",',/'" 19 s 0.19prob c gray81 g4 gray81 g gray81 m gray81', &

              '",',/'" 20 s 0.20prob c gray80 g4 gray80 g gray80 m gray80', &

              '",',/'" 21 s 0.21prob c gray79 g4 gray79 g gray79 m gray79', &

              '",',/'" 22 s 0.22prob c gray78 g4 gray78 g gray78 m gray78', &

              '",',/'" 23 s 0.23prob c gray77 g4 gray77 g gray77 m gray77', &

              '",',/'" 24 s 0.24prob c gray76 g4 gray76 g gray76 m gray76', &

              '",',/'" 25 s 0.25prob c gray75 g4 gray75 g gray75 m gray75', &

              '",',/'" 26 s 0.26prob c gray74 g4 gray74 g gray74 m gray74', &

              '",',/'" 27 s 0.27prob c gray73 g4 gray73 g gray73 m gray73', &

              '",',/'" 28 s 0.28prob c gray72 g4 gray72 g gray72 m gray72', &

              '",',/'" 29 s 0.29prob c gray71 g4 gray71 g gray71 m gray71', &

              '",',/'" 30 s 0.30prob c gray70 g4 gray70 g gray70 m gray70', &

              '",',/'" 31 s 0.31prob c gray69 g4 gray69 g gray69 m gray69', &

              '",',/'" 32 s 0.32prob c gray68 g4 gray68 g gray68 m gray68', &

              '",',/'" 33 s 0.33prob c gray67 g4 gray67 g gray67 m gray67', &

              '",',/'" 34 s 0.34prob c gray66 g4 gray66 g gray66 m gray66', &

              '",',/'" 35 s 0.35prob c gray65 g4 gray65 g gray65 m gray65', &

              '",',/'" 36 s 0.36prob c gray64 g4 gray64 g gray64 m gray64', &

              '",',/'" 37 s 0.37prob c gray63 g4 gray63 g gray63 m gray63', &

              '",',/'" 38 s 0.38prob c gray62 g4 gray62 g gray62 m gray62', &

              '",',/'" 39 s 0.39prob c gray61 g4 gray61 g gray61 m gray61', &

              '",',/'" 40 s 0.40prob c gray60 g4 gray60 g gray60 m gray60', &

              '",',/'" 41 s 0.41prob c gray59 g4 gray59 g gray59 m gray59', &

              '",',/'" 42 s 0.42prob c gray58 g4 gray58 g gray58 m gray58', &

              '",',/'" 43 s 0.43prob c gray57 g4 gray57 g gray57 m gray57', &

              '",',/'" 44 s 0.44prob c gray56 g4 gray56 g gray56 m gray56', &

              '",',/'" 45 s 0.45prob c gray55 g4 gray55 g gray55 m gray55', &

              '",',/'" 46 s 0.46prob c gray54 g4 gray54 g gray54 m gray54', &

              '",',/'" 47 s 0.47prob c gray53 g4 gray53 g gray53 m gray53', &

              '",',/'" 48 s 0.48prob c gray52 g4 gray52 g gray52 m gray52', &

              '",',/'" 49 s 0.49prob c gray51 g4 gray51 g gray51 m gray51', &

              '",',/'" 50 s 0.50prob c gray50 g4 gray50 g gray50 m gray50', &

              '",',/'" 51 s 0.51prob c gray49 g4 gray49 g gray49 m gray49', &

              '",',/'" 52 s 0.52prob c gray48 g4 gray48 g gray48 m gray48', &

              '",',/'" 53 s 0.53prob c gray47 g4 gray47 g gray47 m gray47', &

              '",',/'" 54 s 0.54prob c gray46 g4 gray46 g gray46 m gray46', &

              '",',/'" 55 s 0.55prob c gray45 g4 gray45 g gray45 m gray45', &

              '",',/'" 56 s 0.56prob c gray44 g4 gray44 g gray44 m gray44', &

              '",',/'" 57 s 0.57prob c gray43 g4 gray43 g gray43 m gray43', &

              '",',/'" 58 s 0.58prob c gray42 g4 gray42 g gray42 m gray42', &

              '",',/'" 59 s 0.59prob c gray41 g4 gray41 g gray41 m gray41', &

              '",',/'" 60 s 0.60prob c gray40 g4 gray40 g gray40 m gray40', &

              '",',/'" 61 s 0.61prob c gray39 g4 gray39 g gray39 m gray39', &

              '",',/'" 62 s 0.62prob c gray38 g4 gray38 g gray38 m gray38', &

              '",',/'" 63 s 0.63prob c gray37 g4 gray37 g gray37 m gray37', &

              '",',/'" 64 s 0.64prob c gray36 g4 gray36 g gray36 m gray36', &

              '",',/'" 65 s 0.65prob c gray35 g4 gray35 g gray35 m gray35', &

              '",',/'" 66 s 0.66prob c gray34 g4 gray34 g gray34 m gray34', &

              '",',/'" 67 s 0.67prob c gray33 g4 gray33 g gray33 m gray33', &

              '",',/'" 68 s 0.68prob c gray32 g4 gray32 g gray32 m gray32', &

              '",',/'" 69 s 0.69prob c gray31 g4 gray31 g gray31 m gray31', &

              '",',/'" 70 s 0.70prob c gray30 g4 gray30 g gray30 m gray30', &

              '",',/'" 71 s 0.71prob c gray29 g4 gray29 g gray29 m gray29', &

              '",',/'" 72 s 0.72prob c gray28 g4 gray28 g gray28 m gray28', &

              '",',/'" 73 s 0.73prob c gray27 g4 gray27 g gray27 m gray27', &

              '",',/'" 74 s 0.74prob c gray26 g4 gray26 g gray26 m gray26', &

              '",',/'" 75 s 0.75prob c gray25 g4 gray25 g gray25 m gray25', &

              '",',/'" 76 s 0.76prob c gray24 g4 gray24 g gray24 m gray24', &

              '",',/'" 77 s 0.77prob c gray23 g4 gray23 g gray23 m gray23', &

              '",',/'" 78 s 0.78prob c gray22 g4 gray22 g gray22 m gray22', &

              '",',/'" 79 s 0.79prob c gray21 g4 gray21 g gray21 m gray21', &

              '",',/'" 80 s 0.80prob c gray20 g4 gray20 g gray20 m gray20', &

              '",',/'" 81 s 0.81prob c gray19 g4 gray19 g gray19 m gray19', &

              '",',/'" 82 s 0.82prob c gray18 g4 gray18 g gray18 m gray18', &

              '",',/'" 83 s 0.83prob c gray17 g4 gray17 g gray17 m gray17', &

              '",',/'" 84 s 0.84prob c gray16 g4 gray16 g gray16 m gray16', &

              '",',/'" 85 s 0.85prob c gray15 g4 gray15 g gray15 m gray15', &

              '",',/'" 86 s 0.86prob c gray14 g4 gray14 g gray14 m gray14', &

              '",',/'" 87 s 0.87prob c gray13 g4 gray13 g gray13 m gray13', &

              '",',/'" 88 s 0.88prob c gray12 g4 gray12 g gray12 m gray12', &

              '",',/'" 89 s 0.89prob c gray11 g4 gray11 g gray11 m gray11', &

              '",',/'" 90 s 0.90prob c gray10 g4 gray10 g gray10 m gray10', &

              '",',/'" 91 s 0.91prob c gray9 g4 gray9 g gray9 m gray9', &

              '",',/'" 92 s 0.92prob c gray8 g4 gray8 g gray8 m gray8', &

              '",',/'" 93 s 0.93prob c gray7 g4 gray7 g gray7 m gray7', &

              '",',/'" 94 s 0.94prob c gray6 g4 gray6 g gray6 m gray6', &

              '",',/'" 95 s 0.95prob c gray5 g4 gray5 g gray5 m gray5', &

              '",',/'" 96 s 0.96prob c gray4 g4 gray4 g gray4 m gray4', &

              '",',/'" 97 s 0.97prob c gray3 g4 gray3 g gray3 m gray3', &

              '",',/'" 98 s 0.98prob c gray2 g4 gray2 g gray2 m gray2', &

              '",',/'" 99 s 0.99prob c gray1 g4 gray1 g gray1 m gray1', &

              '",',/'"100 s 1.00prob c black g4 black g black m black', &

              '",',/'"101  c white    m white  s NotFlammable', &

              '",',/'"102  c GreenYellow       m gray70  s LP0', &

              '",',/'"103  c LawnGreen         m gray60 s LP1', &

              '",',/'"104  c LimeGreen         m gray50   s LP2', &

              '",',/'"105  c ForestGreen       m gray40   s LP3', &

              '",',/'"106  c yellow    m gray80  s NF', &

              '",',/'"107  c blue    m black  s water ",', &

              /'/* pixels */')

!

!      calculate a log scale for gray levels

!      use this scale log for uniform scaling across realizations

!       scalelog = 10**(log10(100.)/32767.)

!      scalelog for simple summed priority

!       scale max gray to ihabs-isupply-1 full-on categories

!        (= ihabs-isupply-1 x 32767)

!c       scalelog = 10**(log10(100.)/(32767. * (ihabs-isupply-1)))

!      use this scalelog for maximum priority resolution w/in one map

       scalelog = 10**(log10(100.) / real(maxgray,kind(scalelog)))

!       print *, scalelog

!

!     start at bottom of list (highest priority ties)

      do m = numcells, itiedsplit+1, -1 ! over all tied

!         itemp(P(m,ihabs+1),P(m,ihabs+2)) = 

!     &    (real(P(m,ihabs+3),kind(itemp))/32767)*100

!

!         print *, scalelog**P(m,ihabs+3)

!         print *, P(m,ihabs+1),P(m,ihabs+2),P(m,ihabs+3)

         itemp(P(m,ihabs+1),P(m,ihabs+2)) = scalelog**P(m,ihabs+3)

!         itemp(P(m,ihabs+1),P(m,ihabs+2)) = 100

!        all priorities greater than color index 100 set black

         if (itemp(P(m,ihabs+1),P(m,ihabs+2)) .gt. 100) &

            itemp(P(m,ihabs+1),P(m,ihabs+2)) = 100

      end do

!

!      color used blanks yellow

      do m = mtsplit,mtresolved+1,-1

         itemp(P(m,ihabs+1),P(m,ihabs+2)) = 106

      end do

!

      write (fmt,'("(",I0,"(i3)$)")') irowcols+1

      do i = 0, irowcols

         write (40,'(1A$)') '"'

         !write (40,108) (itemp(j,i), j=0, irowcols)

         write (40,fmt) (itemp(j,i), j=0, irowcols)

         write (40,'(2A)') '",'

!108      format ('"',<irowcols+1>(i3),'",')

      end do

!

      write (40,112)

112   format ('};')

      close (40)

!

      return

   end subroutine XPMTies

!

!------------------------------------------------------------------------------

!

   subroutine XPMRelief(P,iz,iflood)

!

!  output routine for xpm X-window display tool

!

!  write xpm format aspect map of painted fractal

!   spatial probability surfaces with grayshades

!

      implicit none

      integer(i4), intent(in) :: P(maxN*maxN,0:maxsp+3)

      integer(i4), intent(in) :: iz

      integer(i4), intent(in) :: iflood

      ! local variables

      integer(i4)        :: itemp(0:maxN,0:maxN)

      integer(i4)        :: idispvec(0:maxN)

      integer(i4)        :: i, j, k, l, m        ! loop indices

      real(r4)           :: alt                  ! altitude

      real(r4)           :: az                   ! azimuth

      real(r4)           :: x, y

      real(r4)           :: slope                ! slope

      real(r4)           :: aspect               ! aspect

      real(r4)           :: c

      character (len=15) :: outfile(0:maxsp+1)

      character (len=40) :: fmt

!

!     set filenames for the first 10 categories

      outfile(0)  = 'probmap0.xpm   '

      outfile(1)  = 'probmap1.xpm   '

      outfile(2)  = 'probmap2.xpm   '

      outfile(3)  = 'probmap3.xpm   '

      outfile(4)  = 'probmap4.xpm   '

      outfile(5)  = 'probmap5.xpm   '

      outfile(6)  = 'probmap6.xpm   '

      outfile(7)  = 'probmap7.xpm   '

      outfile(8)  = 'probmap8.xpm   '

      outfile(9)  = 'probmap9.xpm   '

      outfile(10) = 'probmap10.xpm  '

!

!

      print 99, outfile(iz)

99    format ('XPM map file ', a15 'being written')

!

!     open the right file for xpm output

      open (50,file=outfile(iz),status='unknown')

!

      write (50, 100)

100   format ('/* XPM */',/'static char *fracasp[]', &

              ' = {',/'/* cols rows ncolors chars_per_pixel', ' */')

      write (50, 102) irowcols+1, irowcols+1

102   format ('"',i6,1x,i6,1x,'108 3 0 0",',/'/* colors */')

      write (50, 104)

104   format ('"  0 s 0.0prob c white g4 white g white m white', &

              '",',/'"  1 s 0.01prob c gray99 g4 gray99 g gray99 m gray99', &

              '",',/'"  2 s 0.02prob c gray98 g4 gray98 g gray98 m gray98', &

              '",',/'"  3 s 0.03prob c gray97 g4 gray97 g gray97 m gray97', &

              '",',/'"  4 s 0.04prob c gray96 g4 gray96 g gray96 m gray96', &

              '",',/'"  5 s 0.05prob c gray95 g4 gray95 g gray95 m gray95', &

              '",',/'"  6 s 0.06prob c gray94 g4 gray94 g gray94 m gray94', &

              '",',/'"  7 s 0.07prob c gray93 g4 gray93 g gray93 m gray93', &

              '",',/'"  8 s 0.08prob c gray92 g4 gray92 g gray92 m gray92', &

              '",',/'"  9 s 0.09prob c gray91 g4 gray91 g gray91 m gray91', &

              '",',/'" 10 s 0.10prob c gray90 g4 gray90 g gray90 m gray90', &

              '",',/'" 11 s 0.11prob c gray89 g4 gray89 g gray89 m gray89', &

              '",',/'" 12 s 0.12prob c gray88 g4 gray88 g gray88 m gray88', &

              '",',/'" 13 s 0.13prob c gray87 g4 gray87 g gray87 m gray87', &

              '",',/'" 14 s 0.14prob c gray86 g4 gray86 g gray86 m gray86', &

              '",',/'" 15 s 0.15prob c gray85 g4 gray85 g gray85 m gray85', &

              '",',/'" 16 s 0.16prob c gray84 g4 gray84 g gray84 m gray84', &

              '",',/'" 17 s 0.17prob c gray83 g4 gray83 g gray83 m gray83', &

              '",',/'" 18 s 0.18prob c gray82 g4 gray82 g gray82 m gray82', &

              '",',/'" 19 s 0.19prob c gray81 g4 gray81 g gray81 m gray81', &

              '",',/'" 20 s 0.20prob c gray80 g4 gray80 g gray80 m gray80', &

              '",',/'" 21 s 0.21prob c gray79 g4 gray79 g gray79 m gray79', &

              '",',/'" 22 s 0.22prob c gray78 g4 gray78 g gray78 m gray78', &

              '",',/'" 23 s 0.23prob c gray77 g4 gray77 g gray77 m gray77', &

              '",',/'" 24 s 0.24prob c gray76 g4 gray76 g gray76 m gray76', &

              '",',/'" 25 s 0.25prob c gray75 g4 gray75 g gray75 m gray75', &

              '",',/'" 26 s 0.26prob c gray74 g4 gray74 g gray74 m gray74', &

              '",',/'" 27 s 0.27prob c gray73 g4 gray73 g gray73 m gray73', &

              '",',/'" 28 s 0.28prob c gray72 g4 gray72 g gray72 m gray72', &

              '",',/'" 29 s 0.29prob c gray71 g4 gray71 g gray71 m gray71', &

              '",',/'" 30 s 0.30prob c gray70 g4 gray70 g gray70 m gray70', &

              '",',/'" 31 s 0.31prob c gray69 g4 gray69 g gray69 m gray69', &

              '",',/'" 32 s 0.32prob c gray68 g4 gray68 g gray68 m gray68', &

              '",',/'" 33 s 0.33prob c gray67 g4 gray67 g gray67 m gray67', &

              '",',/'" 34 s 0.34prob c gray66 g4 gray66 g gray66 m gray66', &

              '",',/'" 35 s 0.35prob c gray65 g4 gray65 g gray65 m gray65', &

              '",',/'" 36 s 0.36prob c gray64 g4 gray64 g gray64 m gray64', &

              '",',/'" 37 s 0.37prob c gray63 g4 gray63 g gray63 m gray63', &

              '",',/'" 38 s 0.38prob c gray62 g4 gray62 g gray62 m gray62', &

              '",',/'" 39 s 0.39prob c gray61 g4 gray61 g gray61 m gray61', &

              '",',/'" 40 s 0.40prob c gray60 g4 gray60 g gray60 m gray60', &

              '",',/'" 41 s 0.41prob c gray59 g4 gray59 g gray59 m gray59', &

              '",',/'" 42 s 0.42prob c gray58 g4 gray58 g gray58 m gray58', &

              '",',/'" 43 s 0.43prob c gray57 g4 gray57 g gray57 m gray57', &

              '",',/'" 44 s 0.44prob c gray56 g4 gray56 g gray56 m gray56', &

              '",',/'" 45 s 0.45prob c gray55 g4 gray55 g gray55 m gray55', &

              '",',/'" 46 s 0.46prob c gray54 g4 gray54 g gray54 m gray54', &

              '",',/'" 47 s 0.47prob c gray53 g4 gray53 g gray53 m gray53', &

              '",',/'" 48 s 0.48prob c gray52 g4 gray52 g gray52 m gray52', &

              '",',/'" 49 s 0.49prob c gray51 g4 gray51 g gray51 m gray51', &

              '",',/'" 50 s 0.50prob c gray50 g4 gray50 g gray50 m gray50', &

              '",',/'" 51 s 0.51prob c gray49 g4 gray49 g gray49 m gray49', &

              '",',/'" 52 s 0.52prob c gray48 g4 gray48 g gray48 m gray48', &

              '",',/'" 53 s 0.53prob c gray47 g4 gray47 g gray47 m gray47', &

              '",',/'" 54 s 0.54prob c gray46 g4 gray46 g gray46 m gray46', &

              '",',/'" 55 s 0.55prob c gray45 g4 gray45 g gray45 m gray45', &

              '",',/'" 56 s 0.56prob c gray44 g4 gray44 g gray44 m gray44', &

              '",',/'" 57 s 0.57prob c gray43 g4 gray43 g gray43 m gray43', &

              '",',/'" 58 s 0.58prob c gray42 g4 gray42 g gray42 m gray42', &

              '",',/'" 59 s 0.59prob c gray41 g4 gray41 g gray41 m gray41', &

              '",',/'" 60 s 0.60prob c gray40 g4 gray40 g gray40 m gray40', &

              '",',/'" 61 s 0.61prob c gray39 g4 gray39 g gray39 m gray39', &

              '",',/'" 62 s 0.62prob c gray38 g4 gray38 g gray38 m gray38', &

              '",',/'" 63 s 0.63prob c gray37 g4 gray37 g gray37 m gray37', &

              '",',/'" 64 s 0.64prob c gray36 g4 gray36 g gray36 m gray36', &

              '",',/'" 65 s 0.65prob c gray35 g4 gray35 g gray35 m gray35', &

              '",',/'" 66 s 0.66prob c gray34 g4 gray34 g gray34 m gray34', &

              '",',/'" 67 s 0.67prob c gray33 g4 gray33 g gray33 m gray33', &

              '",',/'" 68 s 0.68prob c gray32 g4 gray32 g gray32 m gray32', &

              '",',/'" 69 s 0.69prob c gray31 g4 gray31 g gray31 m gray31', &

              '",',/'" 70 s 0.70prob c gray30 g4 gray30 g gray30 m gray30', &

              '",',/'" 71 s 0.71prob c gray29 g4 gray29 g gray29 m gray29', &

              '",',/'" 72 s 0.72prob c gray28 g4 gray28 g gray28 m gray28', &

              '",',/'" 73 s 0.73prob c gray27 g4 gray27 g gray27 m gray27', &

              '",',/'" 74 s 0.74prob c gray26 g4 gray26 g gray26 m gray26', &

              '",',/'" 75 s 0.75prob c gray25 g4 gray25 g gray25 m gray25', &

              '",',/'" 76 s 0.76prob c gray24 g4 gray24 g gray24 m gray24', &

              '",',/'" 77 s 0.77prob c gray23 g4 gray23 g gray23 m gray23', &

              '",',/'" 78 s 0.78prob c gray22 g4 gray22 g gray22 m gray22', &

              '",',/'" 79 s 0.79prob c gray21 g4 gray21 g gray21 m gray21', &

              '",',/'" 80 s 0.80prob c gray20 g4 gray20 g gray20 m gray20', &

              '",',/'" 81 s 0.81prob c gray19 g4 gray19 g gray19 m gray19', &

              '",',/'" 82 s 0.82prob c gray18 g4 gray18 g gray18 m gray18', &

              '",',/'" 83 s 0.83prob c gray17 g4 gray17 g gray17 m gray17', &

              '",',/'" 84 s 0.84prob c gray16 g4 gray16 g gray16 m gray16', &

              '",',/'" 85 s 0.85prob c gray15 g4 gray15 g gray15 m gray15', &

              '",',/'" 86 s 0.86prob c gray14 g4 gray14 g gray14 m gray14', &

              '",',/'" 87 s 0.87prob c gray13 g4 gray13 g gray13 m gray13', &

              '",',/'" 88 s 0.88prob c gray12 g4 gray12 g gray12 m gray12', &

              '",',/'" 89 s 0.89prob c gray11 g4 gray11 g gray11 m gray11', &

              '",',/'" 90 s 0.90prob c gray10 g4 gray10 g gray10 m gray10', &

              '",',/'" 91 s 0.91prob c gray9 g4 gray9 g gray9 m gray9', &

              '",',/'" 92 s 0.92prob c gray8 g4 gray8 g gray8 m gray8', &

              '",',/'" 93 s 0.93prob c gray7 g4 gray7 g gray7 m gray7', &

              '",',/'" 94 s 0.94prob c gray6 g4 gray6 g gray6 m gray6', &

              '",',/'" 95 s 0.95prob c gray5 g4 gray5 g gray5 m gray5', &

              '",',/'" 96 s 0.96prob c gray4 g4 gray4 g gray4 m gray4', &

              '",',/'" 97 s 0.97prob c gray3 g4 gray3 g gray3 m gray3', &

              '",',/'" 98 s 0.98prob c gray2 g4 gray2 g gray2 m gray2', &

              '",',/'" 99 s 0.99prob c gray1 g4 gray1 g gray1 m gray1', &

              '",',/'"100 s 1.00prob c black g4 black g black m black', &

              '",',/'"101  c white    m white  s NotFlammable', &

              '",',/'"102  c GreenYellow       m gray70  s LP0', &

              '",',/'"103  c LawnGreen         m gray60 s LP1', &

              '",',/'"104  c LimeGreen         m gray50   s LP2', &

              '",',/'"105  c ForestGreen       m gray40   s LP3', &

              '",',/'"106  c yellow    m gray80  s NF', &

              '",',/'"107  c blue    m black  s water ",', &

              /'/* pixels */')

!

!

!      load fractal probability terrain into itemp array

      do m = 1, numcells       ! over all cells

         itemp(P(m,ihabs+1),P(m,ihabs+2)) = P(m,iz)

      end do

!

!

!  equations for slope and aspect are from:

!  Horn, B.K.P., 1981.  Hill Shading and the Reflectance Map.

!       Proceedings of the I.E.E.E., 69(1):14-47.

!

!  calculate and display aspect map

!  set sun altitude above horizon in degrees (0 to 90)

      alt = 40.

!  set sun azimuth in degrees west of north (-1 to 360)

      az = 105.

      k = 0

      do j = 0,irowcols

         do i = 0,irowcols

!

            x = (itemp(i-1,j-1)+2.*itemp(i,j-1)+itemp(i+1,j-1) &

                -itemp(i-1,j+1)-2.*itemp(i,j+1)-itemp(i+1,j+1)) &

                /(8.0*30.0)

!