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)

 !

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

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

                 /(8.0*30.0)

 !

             slope = 90. - atand(sqrt(x*x + y*y))

             if (x .eq. 0. .and. y .eq. 0.) then

                aspect = 360.

             else

                aspect = atan2d(x,y)

             end if

 !

             c = sind(alt)*sind(slope) + cosd(alt)*cosd(slope)*cosd(az-aspect)

 !

             if (c .le. 0.) then

                idispvec(k)=0

             else

                idispvec(k)=nint(c*100.)

             end if

 !

 !           flood with water if not painted

             if (itemp(i,j) .le. iflood) idispvec(k) = 107

 !

 !            print *, i,j,k,idispvec(k)

             k = k + 1

          end do

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

          write (50,'(2A$)') ' "'

          !write (50,108) (idispvec(l), l=0, irowcols)

          write (50,fmt) (idispvec(l), l=0, irowcols)

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

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

          k = 0

       end do

 !

       write (50,112)

 112   format ('};')

       close (50)

 !

       return

    end subroutine XPMRelief

 !

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

 !

 !**********************************************************

 ! The following were added by Forrest Hoffman

 ! Sat Apr  7 23:20:48 EDT 2007

 !**********************************************************

 !

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

 !

    real function cosd(angle)

       implicit none

       real, intent(in) :: angle

       cosd = cos(angle * rad_per_deg)

       return

    end function cosd

 !

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

 !

    real function sind(angle)

       implicit none

       real, intent(in) :: angle

       sind = sin(angle * rad_per_deg)

       return

    end function sind

 !

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

 !

    real function atand(x)

       implicit none

       real, intent(in) :: x

       atand = atan(x) * deg_per_rad

       return

    end function atand

 !

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

 !

    real function atan2d(x, y)

       implicit none

       real, intent(in) :: x

       real, intent(in) :: y

       atan2d = atan2(x, y) * deg_per_rad

       return

    end function atan2d

 !

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

 !

 end module realizer_mod