diff -r 000000000000 -r 5fda18b64dcb src/realizer_mod.f90
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/realizer_mod.f90	Wed Sep 26 17:50:53 2007 -0400
@@ -0,0 +1,1091 @@
+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