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
