1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/realizer_mod.f90 Wed Sep 26 17:50:53 2007 -0400
1.3 @@ -0,0 +1,1091 @@
1.4 +module realizer_mod
1.5 +
1.6 +! Realizer module written for the MidPtFM2d program (now called realizer) by
1.7 +! Forrest Hoffman on Fri May 4 16:56:34 EDT 2007 to eliminate duplicated
1.8 +! parameters and common blocks. Parts of the code were converted to
1.9 +! Fortran-90 style and intrinsics. Compiled and tested with gfortran
1.10 +! version 4.1.1 under Fedora Core 6 Linux.
1.11 +
1.12 + implicit none
1.13 + save
1.14 +
1.15 + ! Parameters
1.16 + integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real
1.17 + integer, parameter :: r4 = selected_real_kind( 6) ! 4 byte real
1.18 + integer, parameter :: nr = kind(1.0) ! native real
1.19 + integer, parameter :: i8 = selected_int_kind(13) ! 8 byte integer
1.20 + integer, parameter :: i4 = selected_int_kind( 6) ! 4 byte integer
1.21 + integer, parameter :: ni = kind(1) ! native integer
1.22 + integer(i4), parameter :: maxN=515
1.23 + integer(i4), parameter :: maxsp=10
1.24 + real(r8), parameter :: rad_per_deg=(3.14159265358979323846/180.0)
1.25 + real(r8), parameter :: deg_per_rad=(180.0/3.14159265358979323846)
1.26 +
1.27 + ! Global Variables
1.28 + integer(i4) :: numcells
1.29 + integer(i4) :: itiedsplit
1.30 + integer(i4) :: irowcols
1.31 + integer(i4) :: mtsplit
1.32 + integer(i4) :: mtresolved
1.33 + integer(i4) :: isupply
1.34 + integer(i4) :: ihabs
1.35 + integer(i4) :: maxgray
1.36 + real(r4) :: northeast
1.37 + real(r4) :: southeast
1.38 + real(r4) :: southwest
1.39 + real(r4) :: northwest
1.40 +
1.41 + ! Private Member Functions
1.42 + public dump_map
1.43 + public mpfm2d
1.44 + private StackSort ! unused?
1.45 + private sort ! unused?
1.46 + private SwapDown
1.47 + public Interchange
1.48 + public HeapSort
1.49 + private SD2
1.50 + private I2
1.51 + public HS2
1.52 + private Gauss
1.53 + private cosd
1.54 + private sind
1.55 + private atand
1.56 + private atan2d
1.57 +
1.58 +contains
1.59 +!
1.60 +!------------------------------------------------------------------------------
1.61 +!
1.62 + subroutine dump_map(P)
1.63 + implicit none
1.64 + integer(i4), intent(in) :: P(maxN*maxN,0:maxsp+3)
1.65 + ! local variables
1.66 + integer(i4) :: i, j ! loop indices
1.67 + integer(i4) :: temp(maxsp+3)
1.68 + character (len=40) :: fmt
1.69 +!
1.70 + write(fmt,'("(i5, 1x, ",I0,"(i5, 1x), 1x, 3(i5, 1x))")') ihabs+1
1.71 + do i=1,numcells
1.72 + do j = 0,ihabs+3
1.73 + temp(j) = P(i,j)
1.74 + end do
1.75 + print fmt, i, (temp(j), j=0,ihabs+3)
1.76 + end do
1.77 +!
1.78 + return
1.79 + end subroutine dump_map
1.80 +!
1.81 +!------------------------------------------------------------------------------
1.82 +!
1.83 + real(r4) function f3(delta, x0, x1, x2, sigma)
1.84 + implicit none
1.85 + real(r4), intent(in) :: delta
1.86 + real(r4), intent(in) :: x0, x1, x2
1.87 + real(r4), intent(in) :: sigma
1.88 +
1.89 + f3 = (x0+x1+x2) / 3.0 + delta * Gauss(0.,sigma)
1.90 +
1.91 + return
1.92 + end function f3
1.93 +!
1.94 +!------------------------------------------------------------------------------
1.95 +!
1.96 + real(r4) function f4(delta, x0, x1, x2, x3, sigma)
1.97 + implicit none
1.98 + real(r4), intent(in) :: delta
1.99 + real(r4), intent(in) :: x0, x1, x2, x3
1.100 + real(r4), intent(in) :: sigma
1.101 +
1.102 + f4 = (x0+x1+x2+x3) / 4.0 + delta * Gauss(0.,sigma)
1.103 +
1.104 + return
1.105 + end function f4
1.106 +!
1.107 +!------------------------------------------------------------------------------
1.108 +!
1.109 + subroutine mpfm2d(X, maxlevel, sigma, H, addition, wrap, gradient, iz)
1.110 + implicit none
1.111 + real(r4), intent(inout) :: X(0:maxN,0:maxN)
1.112 + integer(i4), intent(in) :: maxlevel
1.113 + real(r4), intent(in) :: sigma
1.114 + real(r4), intent(in) :: H(maxN,maxsp)
1.115 + logical, intent(in) :: addition
1.116 + logical, intent(in) :: wrap
1.117 + logical, intent(in) :: gradient
1.118 + integer(i4), intent(in) :: iz
1.119 + ! local variables
1.120 + integer(i4) :: N
1.121 + real(r4) :: delta
1.122 + integer(i4) :: DD
1.123 + integer(i4) :: d
1.124 + integer(i4) :: is, ix, iy ! loop indices
1.125 +!
1.126 + print *, "Generating fractal"
1.127 + N = 2**maxlevel
1.128 + delta = sigma
1.129 + X(0,0) = delta * Gauss(0.,sigma)
1.130 + X(0,N) = delta * Gauss(0.,sigma)
1.131 + X(N,0) = delta * Gauss(0.,sigma)
1.132 + X(N,N) = delta * Gauss(0.,sigma)
1.133 +!
1.134 + if (gradient) then
1.135 + X(0,0) = southwest
1.136 + X(0,N) = southeast
1.137 + X(N,0) = northwest
1.138 + X(N,N) = northeast
1.139 + end if
1.140 +!
1.141 + DD = N
1.142 + d = N / 2
1.143 +!
1.144 + do is = 1, maxlevel
1.145 +! going from grid type I to type II
1.146 + delta = delta * 0.5**(0.5 * H(is,iz))
1.147 + do ix = d, N-d, DD
1.148 + do iy = d, N-d, DD
1.149 + X(ix,iy) = f4(delta, X(ix+d,iy+d), X(ix+d,iy-d), &
1.150 + X(ix-d,iy+d), X(ix-d,iy-d), sigma)
1.151 + end do
1.152 + end do
1.153 +! displace other points also if needed
1.154 + if (addition) then
1.155 + do ix = 0, N, DD
1.156 + do iy = 0, N, DD
1.157 + X(ix,iy) = X(ix,iy) + delta * Gauss(0.,sigma)
1.158 + end do
1.159 + end do
1.160 + end if
1.161 +! going from grid type II to type I
1.162 + delta = delta * 0.5**(0.5*H(is,iz))
1.163 +! interpolate and offset boundary grid points
1.164 + do ix = d, N-d, DD
1.165 + X(ix,0) = f3(delta, X(ix+d,0), X(ix-d,0), X(ix,d), sigma)
1.166 + X(ix,N) = f3(delta, X(ix+d,N), X(ix-d,N), X(ix,N-d), sigma)
1.167 + X(0,ix) = f3(delta, X(0,ix+d), X(0,ix-d), X(d,ix), sigma)
1.168 + X(N,ix) = f3(delta, X(N,ix+d), X(N,ix-d), X(N-d,ix), sigma)
1.169 + if (wrap) then
1.170 + X(ix,N) = X(ix,0)
1.171 + X(N,ix) = X(0,ix)
1.172 + end if
1.173 + end do
1.174 +! interpolate and offset inter grid points
1.175 + do ix = d, N-d, DD
1.176 + do iy = DD, N-d, DD
1.177 + X(ix,iy) = f4(delta, X(ix,iy+d), X(ix,iy-d), &
1.178 + X(ix+d,iy), X(ix-d,iy), sigma)
1.179 + end do
1.180 + end do
1.181 + do ix = DD, N-d, DD
1.182 + do iy = d, N-d, DD
1.183 + X(ix,iy) = f4(delta, X(ix,iy+d), X(ix,iy-d), &
1.184 + X(ix+d,iy), X(ix-d,iy), sigma)
1.185 + end do
1.186 + end do
1.187 +! displace other points also if needed
1.188 + if (addition) then
1.189 + do ix = 0, N, DD
1.190 + do iy = 0, N, DD
1.191 + X(ix,iy) = X(ix,iy) + delta * Gauss(0.,sigma)
1.192 + end do
1.193 + end do
1.194 + do ix = d, N-d, DD
1.195 + do iy = d, N-d, DD
1.196 + X(ix,iy) = X(ix,iy) + delta * Gauss(0.,sigma)
1.197 + end do
1.198 + end do
1.199 + end if
1.200 +!
1.201 + DD = DD / 2
1.202 + d = d / 2
1.203 +!
1.204 + end do ! for all levels
1.205 + return
1.206 + end subroutine mpfm2d
1.207 +!
1.208 +!------------------------------------------------------------------------------
1.209 +!
1.210 + real(r4) function Gauss(xm, xs)
1.211 + implicit none
1.212 + real(r4), intent(in) :: xm
1.213 + real(r4), intent(in) :: xs
1.214 + ! local variables
1.215 + real(r4) :: total
1.216 + real(r4) :: p
1.217 + integer(i4) :: i
1.218 +!
1.219 +! create a normal number by summing 12 uniform random numbers
1.220 + total = 0.0
1.221 + do i = 1,12
1.222 + call random_number(p)
1.223 + total = total + p
1.224 + end do
1.225 + Gauss = (total - 6.0) * xs + xm
1.226 +!
1.227 + return
1.228 + end function Gauss
1.229 +!
1.230 +!------------------------------------------------------------------------------
1.231 +!
1.232 + subroutine StackSort(n,arr)
1.233 + implicit none
1.234 + integer(i4), intent(in) :: n
1.235 + real(r4), intent(inout) :: arr(n)
1.236 + ! local variables
1.237 + integer(i4) :: i, j ! loop indices
1.238 + real(r4) :: a
1.239 +!
1.240 +! pick out each element in turn
1.241 + do j = 2,n
1.242 + a = arr(j)
1.243 +! look for the place to insert it
1.244 + do i = j-1,1,-1
1.245 + if (arr(i) <= a) go to 10
1.246 + arr(i+1) = arr(i)
1.247 + end do
1.248 + i = 0
1.249 +! insert the number here
1.250 +10 arr(i+1) = a
1.251 + end do
1.252 + return
1.253 + end subroutine StackSort
1.254 +!
1.255 +!------------------------------------------------------------------------------
1.256 +!
1.257 +! sort routine from Numerical Recipes
1.258 + subroutine sort (n, ra)
1.259 + implicit none
1.260 + integer(i4), intent(in) :: n
1.261 + real(r4), intent(inout) :: ra(n)
1.262 + ! local variables
1.263 + integer(i4) :: i, j, l
1.264 + integer(i4) :: ir
1.265 + real(r4) :: rra
1.266 +
1.267 + l = n/2 + 1
1.268 + ir = n
1.269 +!
1.270 +10 continue
1.271 + if (l .gt. 1) then
1.272 + l = l-1
1.273 + rra = ra(l)
1.274 + else
1.275 + rra = ra(ir)
1.276 + ra(ir) = ra(1)
1.277 + ir = ir-1
1.278 + if (ir .eq. 1) then
1.279 + ra(1) = rra
1.280 + return
1.281 + end if
1.282 + end if
1.283 + i = l
1.284 + j = l+l
1.285 +!
1.286 +20 if (j .le. ir) then
1.287 + if (j .lt. ir) then
1.288 + if (ra(j) .lt. ra(j+1)) j = j+1
1.289 + end if
1.290 + if (rra .lt. ra(j)) then
1.291 + ra(i) = ra(j)
1.292 + i = j
1.293 + j = j+j
1.294 + else
1.295 + j = ir+1
1.296 + end if
1.297 + goto 20
1.298 + end if
1.299 +!
1.300 + ra(i) = rra
1.301 + goto 10
1.302 +!
1.303 + end subroutine sort
1.304 +!
1.305 +!------------------------------------------------------------------------------
1.306 +!
1.307 + subroutine SwapDown(Heap, r, isize, istart, keyfield)
1.308 + implicit none
1.309 + integer(i4), intent(inout) :: Heap(maxN*maxN,0:maxsp+3)
1.310 + integer(i4), intent(inout) :: r
1.311 + integer(i4), intent(in) :: isize
1.312 + integer(i4), intent(in) :: istart
1.313 + integer(i4), intent(in) :: keyfield
1.314 + ! local variables
1.315 + integer(i4) :: child
1.316 + integer(i4) :: ipos
1.317 + logical :: done
1.318 +!
1.319 +! print *, 'entering SwapDown'
1.320 +! print *, 'params are ', r, isize, istart, keyfield
1.321 + ipos = istart - 1
1.322 + done = .false.
1.323 + child = 2 * r
1.324 + do while (.not. done .and. child .le. isize)
1.325 +! find the largest child
1.326 + if (child .lt. isize .and. Heap(child+ipos,keyfield) .lt. &
1.327 + Heap(child+ipos+1,keyfield)) child = child + 1
1.328 +! interchange node and largest child if necessary
1.329 +! and move down to the next subtree
1.330 + if (Heap(r+ipos,keyfield) .lt. Heap(child+ipos,keyfield)) then
1.331 + Call Interchange(Heap,r+ipos,child+ipos)
1.332 + r = child
1.333 + child = 2 * child
1.334 + else
1.335 + done = .true.
1.336 + end if
1.337 +!
1.338 + end do
1.339 +!
1.340 + return
1.341 + end subroutine SwapDown
1.342 +!
1.343 +!------------------------------------------------------------------------------
1.344 +!
1.345 + subroutine Interchange(X,m,n)
1.346 + implicit none
1.347 + integer(i4), intent(inout) :: X(maxN*maxN,0:maxsp+3)
1.348 + integer(i4), intent(in) :: m
1.349 + integer(i4), intent(in) :: n
1.350 + ! local variables
1.351 + integer(i4) :: i
1.352 + integer(i4) :: temp
1.353 +!
1.354 + do i = 0, ihabs+3
1.355 + temp = X(m,i)
1.356 + X(m,i) = X(n,i)
1.357 + X(n,i) = temp
1.358 + end do
1.359 +!
1.360 + return
1.361 + end subroutine Interchange
1.362 +!
1.363 +!------------------------------------------------------------------------------
1.364 +!
1.365 + subroutine HeapSort(X,isize,istart,keyfield)
1.366 + implicit none
1.367 + integer(i4), intent(inout) :: X(maxN*maxN,0:maxsp+3)
1.368 + integer(i4), intent(in) :: isize
1.369 + integer(i4), intent(in) :: istart
1.370 + integer(i4), intent(in) :: keyfield
1.371 + ! local variables
1.372 + integer(i4) :: i ! loop indices
1.373 + integer(i4) :: r
1.374 + integer(i4) :: tempr
1.375 +!
1.376 +! print *, "isize in HeapSort is ", isize
1.377 +! print *, "istart in HeapSort is ", istart
1.378 +! print *, "keyfield in HeapSort is ", keyfield
1.379 +! heapify
1.380 + do r = isize/2, 1, -1
1.381 + tempr = r
1.382 + !Call SwapDown(X,r,isize,istart,keyfield)
1.383 + Call SwapDown(X,tempr,isize,istart,keyfield)
1.384 + !r = tempr
1.385 + end do
1.386 +!
1.387 +! print *, "P after heapify is "
1.388 +! do i=istart,isize
1.389 +! print *, i, X(i,keyfield)
1.390 +! end do
1.391 +!
1.392 +! repeatedly put largest item in root at end of list,
1.393 +! prune it from the tree
1.394 +! and apply SwapDown to the rest of the tree
1.395 +!
1.396 + do i = isize, 2, -1
1.397 + Call Interchange(X,istart,(i+istart-1))
1.398 +! r = isize
1.399 + r = 1
1.400 + Call SwapDown(X,r,i-1,istart,keyfield)
1.401 +!
1.402 +! print *, "P after Swapdown is "
1.403 +! do j=istart,isize
1.404 +! print *, j, X(j,keyfield)
1.405 +! end do
1.406 + end do
1.407 +!
1.408 + return
1.409 + end subroutine HeapSort
1.410 +!
1.411 +!------------------------------------------------------------------------------
1.412 +!
1.413 + subroutine XPMMap(N,II)
1.414 +!
1.415 +! output routine for xpm X-window display tool
1.416 +!
1.417 +!
1.418 + implicit none
1.419 + integer(i4), intent(in) :: N
1.420 + integer(i4), intent(in) :: II(0:maxN,0:maxN)
1.421 + ! local variables
1.422 + integer(i4) :: itemp(0:maxN)
1.423 + integer(i4) :: i, j ! loop indices
1.424 + character (len=40) :: fmt
1.425 +!
1.426 +!
1.427 +! output map
1.428 + print *, "writing final map to xpm format file landscape.xpm"
1.429 +!
1.430 + open (20,file='landscape.xpm',status='unknown')
1.431 +!
1.432 + write (20, 100)
1.433 +100 format ('/* XPM */',/'static char *realizer[]', &
1.434 + ' = {',/'/* cols rows ncolors chars_per_pixel', &
1.435 + ' */')
1.436 + write (20, 102) N+1, N+1
1.437 +102 format ('"',i6,1x,i6,1x,'11 1 0 0",',/'/* colors */')
1.438 + write (20, 104)
1.439 +104 format ('"0 c white m white', &
1.440 + '",',/'"1 s Decid c BurlyWood m black', &
1.441 + '",',/'"2 s Mixed c DarkGoldenrod3 m gray40', &
1.442 + '",',/'"3 s Everg c MistyRose3 m gray50', &
1.443 + '",',/'"4 s Trans c LightGoldenrod4 m gray60', &
1.444 + '",',/'"5 s Urban c DarkGoldenrod1 m gray70', &
1.445 + '",',/'"6 s Water c gold m gray80', &
1.446 + '",',/'"7 s Maint c LightGoldenrod3 m white', &
1.447 + '",',/'"8 s Cropl c SaddleBrown m lightgray', &
1.448 + '",',/'"9 s Lawng c RosyBrown m gray', &
1.449 + '",',/'"* c chocolate m black",',/'/* pixels */')
1.450 +!
1.451 + write(fmt,'("(",I0,"(i1)$)")') N+1
1.452 + do i= 0, N
1.453 + do j = 0, N
1.454 +! itemp(j) = II(i,j)
1.455 + itemp(j) = II(j,i)
1.456 + end do
1.457 + write (20,'(1A$)') '"'
1.458 +! write (20,108) (itemp(j), j=0,N)
1.459 + write (20,fmt) (itemp(j), j=0,N)
1.460 + write (20,'(2A)') '",'
1.461 +!108 format (' "',<N+1>(i1),'",')
1.462 + end do
1.463 +!
1.464 +!
1.465 + write (20,112)
1.466 +112 format ('};')
1.467 +!
1.468 + return
1.469 + end subroutine XPMMap
1.470 +!
1.471 +!------------------------------------------------------------------------------
1.472 +!
1.473 + subroutine InputMap(mapfile,P,N,k,invert)
1.474 +!
1.475 +! input map from a data file
1.476 +!
1.477 + implicit none
1.478 + character (len=60), intent(in) :: mapfile
1.479 + integer(i4), intent(inout) :: P(maxN*maxN,0:maxsp+3)
1.480 + integer(i4), intent(in) :: N
1.481 + integer(i4), intent(in) :: k
1.482 + logical, intent(in) :: invert
1.483 + ! local variables
1.484 + integer(i4) :: i, j, iv ! loop indices
1.485 + real(r4) :: X(0:maxN,0:maxN)
1.486 + real(r4) :: xmin
1.487 + real(r4) :: xmax
1.488 + real(r4) :: xrange
1.489 + real(r4) :: temp
1.490 +!
1.491 + open (14, file=mapfile, status='old')
1.492 +!
1.493 + xmin = 1e30
1.494 + xmax = -1e30
1.495 + do j = 0, N
1.496 + read (14,*) ( X(i,j), i = 0, N)
1.497 +! read (14,118) X(i,j)
1.498 +!118 format (<(N+1)*(N+1)>i5, 1x)
1.499 + do i = 0, N
1.500 + if (X(i,j) .lt. xmin) xmin = X(i,j)
1.501 + if (X(i,j) .gt. xmax) xmax = X(i,j)
1.502 + end do
1.503 + end do
1.504 +!
1.505 +! invert if necessary and scale between 0 and 32767
1.506 + xrange = xmax - xmin
1.507 + do j = 0, N
1.508 + do i = 0, N
1.509 + if (invert) then
1.510 + temp = ((-X(i,j) + xmax) / xrange) * 32767
1.511 + else
1.512 + temp = ((X(i,j) - xmin) / xrange) * 32767
1.513 + end if
1.514 +!
1.515 + iv = (i*(N+1))+j+1
1.516 + P(iv,k) = nint(temp)
1.517 + P(iv,ihabs+1) = i
1.518 + P(iv,ihabs+2) = j
1.519 + end do
1.520 + end do
1.521 +!
1.522 + close(14)
1.523 + rewind(14)
1.524 +!
1.525 + return
1.526 + end subroutine InputMap
1.527 +!
1.528 +!------------------------------------------------------------------------------
1.529 +!
1.530 + subroutine SD2(Heap,r,isize,istart)
1.531 + implicit none
1.532 + integer(i4), intent(inout) :: Heap(maxN*maxN)
1.533 + integer(i4), intent(inout) :: r
1.534 + integer(i4), intent(in) :: isize
1.535 + integer(i4), intent(in) :: istart
1.536 + ! local variables
1.537 + integer(i4) :: ipos
1.538 + integer(i4) :: child
1.539 + logical :: done
1.540 +!
1.541 +! print *, 'entering SwapDown'
1.542 +! print *, 'params are ', r, isize, istart
1.543 + ipos = istart - 1
1.544 + done = .false.
1.545 + child = 2 * r
1.546 + do while (.not. done .and. child .le. isize)
1.547 +! find the largest child
1.548 + if (child .lt. isize .and. Heap(child+ipos) .lt. &
1.549 + Heap(child+ipos+1)) child = child + 1
1.550 +! interchange node and largest child if necessary
1.551 +! and move down to the next subtree
1.552 + if (Heap(r+ipos) .lt. Heap(child+ipos)) then
1.553 + Call I2(Heap,r+ipos,child+ipos)
1.554 + r = child
1.555 + child = 2 * child
1.556 + else
1.557 + done = .true.
1.558 + end if
1.559 +!
1.560 + end do
1.561 +!
1.562 + return
1.563 + end subroutine SD2
1.564 +!
1.565 +!------------------------------------------------------------------------------
1.566 +!
1.567 + subroutine I2(X,m,n)
1.568 + implicit none
1.569 + integer(i4), intent(inout) :: X(maxN*maxN)
1.570 + integer(i4), intent(in) :: m
1.571 + integer(i4), intent(in) :: n
1.572 + ! local variables
1.573 + integer(i4) :: temp
1.574 +!
1.575 + temp = X(m)
1.576 + X(m) = X(n)
1.577 + X(n) = temp
1.578 +!
1.579 + return
1.580 + end subroutine I2
1.581 +!
1.582 +!------------------------------------------------------------------------------
1.583 +!
1.584 + subroutine HS2(X,isize,istart)
1.585 + implicit none
1.586 + integer(i4), intent(inout) :: X(maxN*maxN)
1.587 + integer(i4), intent(in) :: isize
1.588 + integer(i4), intent(in) :: istart
1.589 + ! local variables
1.590 + integer(i4) :: i, j ! loop indices
1.591 + integer(i4) :: r
1.592 + integer(i4) :: tempr
1.593 +!
1.594 +! print *, "isize in HeapSort is ", isize
1.595 +! print *, "istart in HeapSort is ", istart
1.596 +! heapify
1.597 + do r = isize/2, 1, -1
1.598 + tempr = r
1.599 + Call SD2(X,tempr,isize,istart)
1.600 + !Call SD2(X,r,isize,istart)
1.601 + !r = tempr
1.602 + end do
1.603 +!
1.604 +! print*, "P after heapify is "
1.605 +! do i=istart,isize
1.606 +! print *, i, X(i)
1.607 +! end do
1.608 +!
1.609 +! repeatedly put largest item in root at end of list,
1.610 +! prune it from the tree
1.611 +! and apply SwapDown to the rest of the tree
1.612 +!
1.613 + do i = isize, 2, -1
1.614 + Call I2(X,istart,(i+istart-1))
1.615 +! r = isize
1.616 + r = 1
1.617 + Call SD2(X,r,i-1,istart)
1.618 +!
1.619 +! print*, "P after Swapdown is "
1.620 +! do j=istart,isize
1.621 +! print *, j, X(j)
1.622 +! end do
1.623 +
1.624 + end do
1.625 +!
1.626 + return
1.627 + end subroutine HS2
1.628 +!
1.629 +!------------------------------------------------------------------------------
1.630 +!
1.631 + subroutine XPMTies(P)
1.632 +!
1.633 +! output routine for xpm X-window display tool
1.634 +!
1.635 +! write xpm format map with grayshades
1.636 +!
1.637 + implicit none
1.638 + integer(i4), intent(in) :: P(maxN*maxN,0:maxsp+3)
1.639 + ! local variables
1.640 + integer(i4) :: itemp(0:maxN,0:maxN)
1.641 + integer(i4) :: i, j, m ! loop indices
1.642 + real(r4) :: scalelog
1.643 + character (len=40) :: fmt
1.644 +!
1.645 + open (40,file='ties.xpm',status='unknown')
1.646 +!
1.647 + print 99
1.648 +99 format ('XPM map file ties.xpm being written')
1.649 + write (40, 100)
1.650 +100 format ('/* XPM */',/'static char *tiesmap[]', &
1.651 + ' = {',/'/* cols rows ncolors chars_per_pixel', &
1.652 + ' */')
1.653 + write (40, 102) irowcols+1, irowcols+1
1.654 +102 format (' "',i6,1x,i6,1x,'108 3 0 0",',/'/* colors */')
1.655 + write (40, 104)
1.656 +104 format ('" 0 s 0.0prob c white g4 white g white m white', &
1.657 + '",',/'" 1 s 0.01prob c gray99 g4 gray99 g gray99 m gray99', &
1.658 + '",',/'" 2 s 0.02prob c gray98 g4 gray98 g gray98 m gray98', &
1.659 + '",',/'" 3 s 0.03prob c gray97 g4 gray97 g gray97 m gray97', &
1.660 + '",',/'" 4 s 0.04prob c gray96 g4 gray96 g gray96 m gray96', &
1.661 + '",',/'" 5 s 0.05prob c gray95 g4 gray95 g gray95 m gray95', &
1.662 + '",',/'" 6 s 0.06prob c gray94 g4 gray94 g gray94 m gray94', &
1.663 + '",',/'" 7 s 0.07prob c gray93 g4 gray93 g gray93 m gray93', &
1.664 + '",',/'" 8 s 0.08prob c gray92 g4 gray92 g gray92 m gray92', &
1.665 + '",',/'" 9 s 0.09prob c gray91 g4 gray91 g gray91 m gray91', &
1.666 + '",',/'" 10 s 0.10prob c gray90 g4 gray90 g gray90 m gray90', &
1.667 + '",',/'" 11 s 0.11prob c gray89 g4 gray89 g gray89 m gray89', &
1.668 + '",',/'" 12 s 0.12prob c gray88 g4 gray88 g gray88 m gray88', &
1.669 + '",',/'" 13 s 0.13prob c gray87 g4 gray87 g gray87 m gray87', &
1.670 + '",',/'" 14 s 0.14prob c gray86 g4 gray86 g gray86 m gray86', &
1.671 + '",',/'" 15 s 0.15prob c gray85 g4 gray85 g gray85 m gray85', &
1.672 + '",',/'" 16 s 0.16prob c gray84 g4 gray84 g gray84 m gray84', &
1.673 + '",',/'" 17 s 0.17prob c gray83 g4 gray83 g gray83 m gray83', &
1.674 + '",',/'" 18 s 0.18prob c gray82 g4 gray82 g gray82 m gray82', &
1.675 + '",',/'" 19 s 0.19prob c gray81 g4 gray81 g gray81 m gray81', &
1.676 + '",',/'" 20 s 0.20prob c gray80 g4 gray80 g gray80 m gray80', &
1.677 + '",',/'" 21 s 0.21prob c gray79 g4 gray79 g gray79 m gray79', &
1.678 + '",',/'" 22 s 0.22prob c gray78 g4 gray78 g gray78 m gray78', &
1.679 + '",',/'" 23 s 0.23prob c gray77 g4 gray77 g gray77 m gray77', &
1.680 + '",',/'" 24 s 0.24prob c gray76 g4 gray76 g gray76 m gray76', &
1.681 + '",',/'" 25 s 0.25prob c gray75 g4 gray75 g gray75 m gray75', &
1.682 + '",',/'" 26 s 0.26prob c gray74 g4 gray74 g gray74 m gray74', &
1.683 + '",',/'" 27 s 0.27prob c gray73 g4 gray73 g gray73 m gray73', &
1.684 + '",',/'" 28 s 0.28prob c gray72 g4 gray72 g gray72 m gray72', &
1.685 + '",',/'" 29 s 0.29prob c gray71 g4 gray71 g gray71 m gray71', &
1.686 + '",',/'" 30 s 0.30prob c gray70 g4 gray70 g gray70 m gray70', &
1.687 + '",',/'" 31 s 0.31prob c gray69 g4 gray69 g gray69 m gray69', &
1.688 + '",',/'" 32 s 0.32prob c gray68 g4 gray68 g gray68 m gray68', &
1.689 + '",',/'" 33 s 0.33prob c gray67 g4 gray67 g gray67 m gray67', &
1.690 + '",',/'" 34 s 0.34prob c gray66 g4 gray66 g gray66 m gray66', &
1.691 + '",',/'" 35 s 0.35prob c gray65 g4 gray65 g gray65 m gray65', &
1.692 + '",',/'" 36 s 0.36prob c gray64 g4 gray64 g gray64 m gray64', &
1.693 + '",',/'" 37 s 0.37prob c gray63 g4 gray63 g gray63 m gray63', &
1.694 + '",',/'" 38 s 0.38prob c gray62 g4 gray62 g gray62 m gray62', &
1.695 + '",',/'" 39 s 0.39prob c gray61 g4 gray61 g gray61 m gray61', &
1.696 + '",',/'" 40 s 0.40prob c gray60 g4 gray60 g gray60 m gray60', &
1.697 + '",',/'" 41 s 0.41prob c gray59 g4 gray59 g gray59 m gray59', &
1.698 + '",',/'" 42 s 0.42prob c gray58 g4 gray58 g gray58 m gray58', &
1.699 + '",',/'" 43 s 0.43prob c gray57 g4 gray57 g gray57 m gray57', &
1.700 + '",',/'" 44 s 0.44prob c gray56 g4 gray56 g gray56 m gray56', &
1.701 + '",',/'" 45 s 0.45prob c gray55 g4 gray55 g gray55 m gray55', &
1.702 + '",',/'" 46 s 0.46prob c gray54 g4 gray54 g gray54 m gray54', &
1.703 + '",',/'" 47 s 0.47prob c gray53 g4 gray53 g gray53 m gray53', &
1.704 + '",',/'" 48 s 0.48prob c gray52 g4 gray52 g gray52 m gray52', &
1.705 + '",',/'" 49 s 0.49prob c gray51 g4 gray51 g gray51 m gray51', &
1.706 + '",',/'" 50 s 0.50prob c gray50 g4 gray50 g gray50 m gray50', &
1.707 + '",',/'" 51 s 0.51prob c gray49 g4 gray49 g gray49 m gray49', &
1.708 + '",',/'" 52 s 0.52prob c gray48 g4 gray48 g gray48 m gray48', &
1.709 + '",',/'" 53 s 0.53prob c gray47 g4 gray47 g gray47 m gray47', &
1.710 + '",',/'" 54 s 0.54prob c gray46 g4 gray46 g gray46 m gray46', &
1.711 + '",',/'" 55 s 0.55prob c gray45 g4 gray45 g gray45 m gray45', &
1.712 + '",',/'" 56 s 0.56prob c gray44 g4 gray44 g gray44 m gray44', &
1.713 + '",',/'" 57 s 0.57prob c gray43 g4 gray43 g gray43 m gray43', &
1.714 + '",',/'" 58 s 0.58prob c gray42 g4 gray42 g gray42 m gray42', &
1.715 + '",',/'" 59 s 0.59prob c gray41 g4 gray41 g gray41 m gray41', &
1.716 + '",',/'" 60 s 0.60prob c gray40 g4 gray40 g gray40 m gray40', &
1.717 + '",',/'" 61 s 0.61prob c gray39 g4 gray39 g gray39 m gray39', &
1.718 + '",',/'" 62 s 0.62prob c gray38 g4 gray38 g gray38 m gray38', &
1.719 + '",',/'" 63 s 0.63prob c gray37 g4 gray37 g gray37 m gray37', &
1.720 + '",',/'" 64 s 0.64prob c gray36 g4 gray36 g gray36 m gray36', &
1.721 + '",',/'" 65 s 0.65prob c gray35 g4 gray35 g gray35 m gray35', &
1.722 + '",',/'" 66 s 0.66prob c gray34 g4 gray34 g gray34 m gray34', &
1.723 + '",',/'" 67 s 0.67prob c gray33 g4 gray33 g gray33 m gray33', &
1.724 + '",',/'" 68 s 0.68prob c gray32 g4 gray32 g gray32 m gray32', &
1.725 + '",',/'" 69 s 0.69prob c gray31 g4 gray31 g gray31 m gray31', &
1.726 + '",',/'" 70 s 0.70prob c gray30 g4 gray30 g gray30 m gray30', &
1.727 + '",',/'" 71 s 0.71prob c gray29 g4 gray29 g gray29 m gray29', &
1.728 + '",',/'" 72 s 0.72prob c gray28 g4 gray28 g gray28 m gray28', &
1.729 + '",',/'" 73 s 0.73prob c gray27 g4 gray27 g gray27 m gray27', &
1.730 + '",',/'" 74 s 0.74prob c gray26 g4 gray26 g gray26 m gray26', &
1.731 + '",',/'" 75 s 0.75prob c gray25 g4 gray25 g gray25 m gray25', &
1.732 + '",',/'" 76 s 0.76prob c gray24 g4 gray24 g gray24 m gray24', &
1.733 + '",',/'" 77 s 0.77prob c gray23 g4 gray23 g gray23 m gray23', &
1.734 + '",',/'" 78 s 0.78prob c gray22 g4 gray22 g gray22 m gray22', &
1.735 + '",',/'" 79 s 0.79prob c gray21 g4 gray21 g gray21 m gray21', &
1.736 + '",',/'" 80 s 0.80prob c gray20 g4 gray20 g gray20 m gray20', &
1.737 + '",',/'" 81 s 0.81prob c gray19 g4 gray19 g gray19 m gray19', &
1.738 + '",',/'" 82 s 0.82prob c gray18 g4 gray18 g gray18 m gray18', &
1.739 + '",',/'" 83 s 0.83prob c gray17 g4 gray17 g gray17 m gray17', &
1.740 + '",',/'" 84 s 0.84prob c gray16 g4 gray16 g gray16 m gray16', &
1.741 + '",',/'" 85 s 0.85prob c gray15 g4 gray15 g gray15 m gray15', &
1.742 + '",',/'" 86 s 0.86prob c gray14 g4 gray14 g gray14 m gray14', &
1.743 + '",',/'" 87 s 0.87prob c gray13 g4 gray13 g gray13 m gray13', &
1.744 + '",',/'" 88 s 0.88prob c gray12 g4 gray12 g gray12 m gray12', &
1.745 + '",',/'" 89 s 0.89prob c gray11 g4 gray11 g gray11 m gray11', &
1.746 + '",',/'" 90 s 0.90prob c gray10 g4 gray10 g gray10 m gray10', &
1.747 + '",',/'" 91 s 0.91prob c gray9 g4 gray9 g gray9 m gray9', &
1.748 + '",',/'" 92 s 0.92prob c gray8 g4 gray8 g gray8 m gray8', &
1.749 + '",',/'" 93 s 0.93prob c gray7 g4 gray7 g gray7 m gray7', &
1.750 + '",',/'" 94 s 0.94prob c gray6 g4 gray6 g gray6 m gray6', &
1.751 + '",',/'" 95 s 0.95prob c gray5 g4 gray5 g gray5 m gray5', &
1.752 + '",',/'" 96 s 0.96prob c gray4 g4 gray4 g gray4 m gray4', &
1.753 + '",',/'" 97 s 0.97prob c gray3 g4 gray3 g gray3 m gray3', &
1.754 + '",',/'" 98 s 0.98prob c gray2 g4 gray2 g gray2 m gray2', &
1.755 + '",',/'" 99 s 0.99prob c gray1 g4 gray1 g gray1 m gray1', &
1.756 + '",',/'"100 s 1.00prob c black g4 black g black m black', &
1.757 + '",',/'"101 c white m white s NotFlammable', &
1.758 + '",',/'"102 c GreenYellow m gray70 s LP0', &
1.759 + '",',/'"103 c LawnGreen m gray60 s LP1', &
1.760 + '",',/'"104 c LimeGreen m gray50 s LP2', &
1.761 + '",',/'"105 c ForestGreen m gray40 s LP3', &
1.762 + '",',/'"106 c yellow m gray80 s NF', &
1.763 + '",',/'"107 c blue m black s water ",', &
1.764 + /'/* pixels */')
1.765 +!
1.766 +! calculate a log scale for gray levels
1.767 +! use this scale log for uniform scaling across realizations
1.768 +! scalelog = 10**(log10(100.)/32767.)
1.769 +! scalelog for simple summed priority
1.770 +! scale max gray to ihabs-isupply-1 full-on categories
1.771 +! (= ihabs-isupply-1 x 32767)
1.772 +!c scalelog = 10**(log10(100.)/(32767. * (ihabs-isupply-1)))
1.773 +! use this scalelog for maximum priority resolution w/in one map
1.774 + scalelog = 10**(log10(100.) / real(maxgray,kind(scalelog)))
1.775 +! print *, scalelog
1.776 +!
1.777 +! start at bottom of list (highest priority ties)
1.778 + do m = numcells, itiedsplit+1, -1 ! over all tied
1.779 +! itemp(P(m,ihabs+1),P(m,ihabs+2)) =
1.780 +! & (real(P(m,ihabs+3),kind(itemp))/32767)*100
1.781 +!
1.782 +! print *, scalelog**P(m,ihabs+3)
1.783 +! print *, P(m,ihabs+1),P(m,ihabs+2),P(m,ihabs+3)
1.784 + itemp(P(m,ihabs+1),P(m,ihabs+2)) = scalelog**P(m,ihabs+3)
1.785 +! itemp(P(m,ihabs+1),P(m,ihabs+2)) = 100
1.786 +! all priorities greater than color index 100 set black
1.787 + if (itemp(P(m,ihabs+1),P(m,ihabs+2)) .gt. 100) &
1.788 + itemp(P(m,ihabs+1),P(m,ihabs+2)) = 100
1.789 + end do
1.790 +!
1.791 +! color used blanks yellow
1.792 + do m = mtsplit,mtresolved+1,-1
1.793 + itemp(P(m,ihabs+1),P(m,ihabs+2)) = 106
1.794 + end do
1.795 +!
1.796 + write (fmt,'("(",I0,"(i3)$)")') irowcols+1
1.797 + do i = 0, irowcols
1.798 + write (40,'(1A$)') '"'
1.799 + !write (40,108) (itemp(j,i), j=0, irowcols)
1.800 + write (40,fmt) (itemp(j,i), j=0, irowcols)
1.801 + write (40,'(2A)') '",'
1.802 +!108 format ('"',<irowcols+1>(i3),'",')
1.803 + end do
1.804 +!
1.805 + write (40,112)
1.806 +112 format ('};')
1.807 + close (40)
1.808 +!
1.809 + return
1.810 + end subroutine XPMTies
1.811 +!
1.812 +!------------------------------------------------------------------------------
1.813 +!
1.814 + subroutine XPMRelief(P,iz,iflood)
1.815 +!
1.816 +! output routine for xpm X-window display tool
1.817 +!
1.818 +! write xpm format aspect map of painted fractal
1.819 +! spatial probability surfaces with grayshades
1.820 +!
1.821 + implicit none
1.822 + integer(i4), intent(in) :: P(maxN*maxN,0:maxsp+3)
1.823 + integer(i4), intent(in) :: iz
1.824 + integer(i4), intent(in) :: iflood
1.825 + ! local variables
1.826 + integer(i4) :: itemp(0:maxN,0:maxN)
1.827 + integer(i4) :: idispvec(0:maxN)
1.828 + integer(i4) :: i, j, k, l, m ! loop indices
1.829 + real(r4) :: alt ! altitude
1.830 + real(r4) :: az ! azimuth
1.831 + real(r4) :: x, y
1.832 + real(r4) :: slope ! slope
1.833 + real(r4) :: aspect ! aspect
1.834 + real(r4) :: c
1.835 + character (len=15) :: outfile(0:maxsp+1)
1.836 + character (len=40) :: fmt
1.837 +!
1.838 +! set filenames for the first 10 categories
1.839 + outfile(0) = 'probmap0.xpm '
1.840 + outfile(1) = 'probmap1.xpm '
1.841 + outfile(2) = 'probmap2.xpm '
1.842 + outfile(3) = 'probmap3.xpm '
1.843 + outfile(4) = 'probmap4.xpm '
1.844 + outfile(5) = 'probmap5.xpm '
1.845 + outfile(6) = 'probmap6.xpm '
1.846 + outfile(7) = 'probmap7.xpm '
1.847 + outfile(8) = 'probmap8.xpm '
1.848 + outfile(9) = 'probmap9.xpm '
1.849 + outfile(10) = 'probmap10.xpm '
1.850 +!
1.851 +!
1.852 + print 99, outfile(iz)
1.853 +99 format ('XPM map file ', a15 'being written')
1.854 +!
1.855 +! open the right file for xpm output
1.856 + open (50,file=outfile(iz),status='unknown')
1.857 +!
1.858 + write (50, 100)
1.859 +100 format ('/* XPM */',/'static char *fracasp[]', &
1.860 + ' = {',/'/* cols rows ncolors chars_per_pixel', ' */')
1.861 + write (50, 102) irowcols+1, irowcols+1
1.862 +102 format ('"',i6,1x,i6,1x,'108 3 0 0",',/'/* colors */')
1.863 + write (50, 104)
1.864 +104 format ('" 0 s 0.0prob c white g4 white g white m white', &
1.865 + '",',/'" 1 s 0.01prob c gray99 g4 gray99 g gray99 m gray99', &
1.866 + '",',/'" 2 s 0.02prob c gray98 g4 gray98 g gray98 m gray98', &
1.867 + '",',/'" 3 s 0.03prob c gray97 g4 gray97 g gray97 m gray97', &
1.868 + '",',/'" 4 s 0.04prob c gray96 g4 gray96 g gray96 m gray96', &
1.869 + '",',/'" 5 s 0.05prob c gray95 g4 gray95 g gray95 m gray95', &
1.870 + '",',/'" 6 s 0.06prob c gray94 g4 gray94 g gray94 m gray94', &
1.871 + '",',/'" 7 s 0.07prob c gray93 g4 gray93 g gray93 m gray93', &
1.872 + '",',/'" 8 s 0.08prob c gray92 g4 gray92 g gray92 m gray92', &
1.873 + '",',/'" 9 s 0.09prob c gray91 g4 gray91 g gray91 m gray91', &
1.874 + '",',/'" 10 s 0.10prob c gray90 g4 gray90 g gray90 m gray90', &
1.875 + '",',/'" 11 s 0.11prob c gray89 g4 gray89 g gray89 m gray89', &
1.876 + '",',/'" 12 s 0.12prob c gray88 g4 gray88 g gray88 m gray88', &
1.877 + '",',/'" 13 s 0.13prob c gray87 g4 gray87 g gray87 m gray87', &
1.878 + '",',/'" 14 s 0.14prob c gray86 g4 gray86 g gray86 m gray86', &
1.879 + '",',/'" 15 s 0.15prob c gray85 g4 gray85 g gray85 m gray85', &
1.880 + '",',/'" 16 s 0.16prob c gray84 g4 gray84 g gray84 m gray84', &
1.881 + '",',/'" 17 s 0.17prob c gray83 g4 gray83 g gray83 m gray83', &
1.882 + '",',/'" 18 s 0.18prob c gray82 g4 gray82 g gray82 m gray82', &
1.883 + '",',/'" 19 s 0.19prob c gray81 g4 gray81 g gray81 m gray81', &
1.884 + '",',/'" 20 s 0.20prob c gray80 g4 gray80 g gray80 m gray80', &
1.885 + '",',/'" 21 s 0.21prob c gray79 g4 gray79 g gray79 m gray79', &
1.886 + '",',/'" 22 s 0.22prob c gray78 g4 gray78 g gray78 m gray78', &
1.887 + '",',/'" 23 s 0.23prob c gray77 g4 gray77 g gray77 m gray77', &
1.888 + '",',/'" 24 s 0.24prob c gray76 g4 gray76 g gray76 m gray76', &
1.889 + '",',/'" 25 s 0.25prob c gray75 g4 gray75 g gray75 m gray75', &
1.890 + '",',/'" 26 s 0.26prob c gray74 g4 gray74 g gray74 m gray74', &
1.891 + '",',/'" 27 s 0.27prob c gray73 g4 gray73 g gray73 m gray73', &
1.892 + '",',/'" 28 s 0.28prob c gray72 g4 gray72 g gray72 m gray72', &
1.893 + '",',/'" 29 s 0.29prob c gray71 g4 gray71 g gray71 m gray71', &
1.894 + '",',/'" 30 s 0.30prob c gray70 g4 gray70 g gray70 m gray70', &
1.895 + '",',/'" 31 s 0.31prob c gray69 g4 gray69 g gray69 m gray69', &
1.896 + '",',/'" 32 s 0.32prob c gray68 g4 gray68 g gray68 m gray68', &
1.897 + '",',/'" 33 s 0.33prob c gray67 g4 gray67 g gray67 m gray67', &
1.898 + '",',/'" 34 s 0.34prob c gray66 g4 gray66 g gray66 m gray66', &
1.899 + '",',/'" 35 s 0.35prob c gray65 g4 gray65 g gray65 m gray65', &
1.900 + '",',/'" 36 s 0.36prob c gray64 g4 gray64 g gray64 m gray64', &
1.901 + '",',/'" 37 s 0.37prob c gray63 g4 gray63 g gray63 m gray63', &
1.902 + '",',/'" 38 s 0.38prob c gray62 g4 gray62 g gray62 m gray62', &
1.903 + '",',/'" 39 s 0.39prob c gray61 g4 gray61 g gray61 m gray61', &
1.904 + '",',/'" 40 s 0.40prob c gray60 g4 gray60 g gray60 m gray60', &
1.905 + '",',/'" 41 s 0.41prob c gray59 g4 gray59 g gray59 m gray59', &
1.906 + '",',/'" 42 s 0.42prob c gray58 g4 gray58 g gray58 m gray58', &
1.907 + '",',/'" 43 s 0.43prob c gray57 g4 gray57 g gray57 m gray57', &
1.908 + '",',/'" 44 s 0.44prob c gray56 g4 gray56 g gray56 m gray56', &
1.909 + '",',/'" 45 s 0.45prob c gray55 g4 gray55 g gray55 m gray55', &
1.910 + '",',/'" 46 s 0.46prob c gray54 g4 gray54 g gray54 m gray54', &
1.911 + '",',/'" 47 s 0.47prob c gray53 g4 gray53 g gray53 m gray53', &
1.912 + '",',/'" 48 s 0.48prob c gray52 g4 gray52 g gray52 m gray52', &
1.913 + '",',/'" 49 s 0.49prob c gray51 g4 gray51 g gray51 m gray51', &
1.914 + '",',/'" 50 s 0.50prob c gray50 g4 gray50 g gray50 m gray50', &
1.915 + '",',/'" 51 s 0.51prob c gray49 g4 gray49 g gray49 m gray49', &
1.916 + '",',/'" 52 s 0.52prob c gray48 g4 gray48 g gray48 m gray48', &
1.917 + '",',/'" 53 s 0.53prob c gray47 g4 gray47 g gray47 m gray47', &
1.918 + '",',/'" 54 s 0.54prob c gray46 g4 gray46 g gray46 m gray46', &
1.919 + '",',/'" 55 s 0.55prob c gray45 g4 gray45 g gray45 m gray45', &
1.920 + '",',/'" 56 s 0.56prob c gray44 g4 gray44 g gray44 m gray44', &
1.921 + '",',/'" 57 s 0.57prob c gray43 g4 gray43 g gray43 m gray43', &
1.922 + '",',/'" 58 s 0.58prob c gray42 g4 gray42 g gray42 m gray42', &
1.923 + '",',/'" 59 s 0.59prob c gray41 g4 gray41 g gray41 m gray41', &
1.924 + '",',/'" 60 s 0.60prob c gray40 g4 gray40 g gray40 m gray40', &
1.925 + '",',/'" 61 s 0.61prob c gray39 g4 gray39 g gray39 m gray39', &
1.926 + '",',/'" 62 s 0.62prob c gray38 g4 gray38 g gray38 m gray38', &
1.927 + '",',/'" 63 s 0.63prob c gray37 g4 gray37 g gray37 m gray37', &
1.928 + '",',/'" 64 s 0.64prob c gray36 g4 gray36 g gray36 m gray36', &
1.929 + '",',/'" 65 s 0.65prob c gray35 g4 gray35 g gray35 m gray35', &
1.930 + '",',/'" 66 s 0.66prob c gray34 g4 gray34 g gray34 m gray34', &
1.931 + '",',/'" 67 s 0.67prob c gray33 g4 gray33 g gray33 m gray33', &
1.932 + '",',/'" 68 s 0.68prob c gray32 g4 gray32 g gray32 m gray32', &
1.933 + '",',/'" 69 s 0.69prob c gray31 g4 gray31 g gray31 m gray31', &
1.934 + '",',/'" 70 s 0.70prob c gray30 g4 gray30 g gray30 m gray30', &
1.935 + '",',/'" 71 s 0.71prob c gray29 g4 gray29 g gray29 m gray29', &
1.936 + '",',/'" 72 s 0.72prob c gray28 g4 gray28 g gray28 m gray28', &
1.937 + '",',/'" 73 s 0.73prob c gray27 g4 gray27 g gray27 m gray27', &
1.938 + '",',/'" 74 s 0.74prob c gray26 g4 gray26 g gray26 m gray26', &
1.939 + '",',/'" 75 s 0.75prob c gray25 g4 gray25 g gray25 m gray25', &
1.940 + '",',/'" 76 s 0.76prob c gray24 g4 gray24 g gray24 m gray24', &
1.941 + '",',/'" 77 s 0.77prob c gray23 g4 gray23 g gray23 m gray23', &
1.942 + '",',/'" 78 s 0.78prob c gray22 g4 gray22 g gray22 m gray22', &
1.943 + '",',/'" 79 s 0.79prob c gray21 g4 gray21 g gray21 m gray21', &
1.944 + '",',/'" 80 s 0.80prob c gray20 g4 gray20 g gray20 m gray20', &
1.945 + '",',/'" 81 s 0.81prob c gray19 g4 gray19 g gray19 m gray19', &
1.946 + '",',/'" 82 s 0.82prob c gray18 g4 gray18 g gray18 m gray18', &
1.947 + '",',/'" 83 s 0.83prob c gray17 g4 gray17 g gray17 m gray17', &
1.948 + '",',/'" 84 s 0.84prob c gray16 g4 gray16 g gray16 m gray16', &
1.949 + '",',/'" 85 s 0.85prob c gray15 g4 gray15 g gray15 m gray15', &
1.950 + '",',/'" 86 s 0.86prob c gray14 g4 gray14 g gray14 m gray14', &
1.951 + '",',/'" 87 s 0.87prob c gray13 g4 gray13 g gray13 m gray13', &
1.952 + '",',/'" 88 s 0.88prob c gray12 g4 gray12 g gray12 m gray12', &
1.953 + '",',/'" 89 s 0.89prob c gray11 g4 gray11 g gray11 m gray11', &
1.954 + '",',/'" 90 s 0.90prob c gray10 g4 gray10 g gray10 m gray10', &
1.955 + '",',/'" 91 s 0.91prob c gray9 g4 gray9 g gray9 m gray9', &
1.956 + '",',/'" 92 s 0.92prob c gray8 g4 gray8 g gray8 m gray8', &
1.957 + '",',/'" 93 s 0.93prob c gray7 g4 gray7 g gray7 m gray7', &
1.958 + '",',/'" 94 s 0.94prob c gray6 g4 gray6 g gray6 m gray6', &
1.959 + '",',/'" 95 s 0.95prob c gray5 g4 gray5 g gray5 m gray5', &
1.960 + '",',/'" 96 s 0.96prob c gray4 g4 gray4 g gray4 m gray4', &
1.961 + '",',/'" 97 s 0.97prob c gray3 g4 gray3 g gray3 m gray3', &
1.962 + '",',/'" 98 s 0.98prob c gray2 g4 gray2 g gray2 m gray2', &
1.963 + '",',/'" 99 s 0.99prob c gray1 g4 gray1 g gray1 m gray1', &
1.964 + '",',/'"100 s 1.00prob c black g4 black g black m black', &
1.965 + '",',/'"101 c white m white s NotFlammable', &
1.966 + '",',/'"102 c GreenYellow m gray70 s LP0', &
1.967 + '",',/'"103 c LawnGreen m gray60 s LP1', &
1.968 + '",',/'"104 c LimeGreen m gray50 s LP2', &
1.969 + '",',/'"105 c ForestGreen m gray40 s LP3', &
1.970 + '",',/'"106 c yellow m gray80 s NF', &
1.971 + '",',/'"107 c blue m black s water ",', &
1.972 + /'/* pixels */')
1.973 +!
1.974 +!
1.975 +! load fractal probability terrain into itemp array
1.976 + do m = 1, numcells ! over all cells
1.977 + itemp(P(m,ihabs+1),P(m,ihabs+2)) = P(m,iz)
1.978 + end do
1.979 +!
1.980 +!
1.981 +! equations for slope and aspect are from:
1.982 +! Horn, B.K.P., 1981. Hill Shading and the Reflectance Map.
1.983 +! Proceedings of the I.E.E.E., 69(1):14-47.
1.984 +!
1.985 +! calculate and display aspect map
1.986 +! set sun altitude above horizon in degrees (0 to 90)
1.987 + alt = 40.
1.988 +! set sun azimuth in degrees west of north (-1 to 360)
1.989 + az = 105.
1.990 + k = 0
1.991 + do j = 0,irowcols
1.992 + do i = 0,irowcols
1.993 +!
1.994 + x = (itemp(i-1,j-1)+2.*itemp(i,j-1)+itemp(i+1,j-1) &
1.995 + -itemp(i-1,j+1)-2.*itemp(i,j+1)-itemp(i+1,j+1)) &
1.996 + /(8.0*30.0)
1.997 +!
1.998 + y = (itemp(i-1,j-1)+2.*itemp(i-1,j)+itemp(i-1,j+1) &
1.999 + -itemp(i+1,j-1)-2.*itemp(i+1,j)-itemp(i+1,j+1)) &
1.1000 + /(8.0*30.0)
1.1001 +!
1.1002 + slope = 90. - atand(sqrt(x*x + y*y))
1.1003 + if (x .eq. 0. .and. y .eq. 0.) then
1.1004 + aspect = 360.
1.1005 + else
1.1006 + aspect = atan2d(x,y)
1.1007 + end if
1.1008 +!
1.1009 + c = sind(alt)*sind(slope) + cosd(alt)*cosd(slope)*cosd(az-aspect)
1.1010 +!
1.1011 + if (c .le. 0.) then
1.1012 + idispvec(k)=0
1.1013 + else
1.1014 + idispvec(k)=nint(c*100.)
1.1015 + end if
1.1016 +!
1.1017 +! flood with water if not painted
1.1018 + if (itemp(i,j) .le. iflood) idispvec(k) = 107
1.1019 +!
1.1020 +! print *, i,j,k,idispvec(k)
1.1021 + k = k + 1
1.1022 + end do
1.1023 + write (fmt,'("(",I0,"(i3)$)")') irowcols+1
1.1024 + write (50,'(2A$)') ' "'
1.1025 + !write (50,108) (idispvec(l), l=0, irowcols)
1.1026 + write (50,fmt) (idispvec(l), l=0, irowcols)
1.1027 + write (50,'(2A)') '",'
1.1028 +!108 format (' "',<irowcols+1>(i3),'",')
1.1029 + k = 0
1.1030 + end do
1.1031 +!
1.1032 + write (50,112)
1.1033 +112 format ('};')
1.1034 + close (50)
1.1035 +!
1.1036 + return
1.1037 + end subroutine XPMRelief
1.1038 +!
1.1039 +!------------------------------------------------------------------------------
1.1040 +!
1.1041 +!**********************************************************
1.1042 +! The following were added by Forrest Hoffman
1.1043 +! Sat Apr 7 23:20:48 EDT 2007
1.1044 +!**********************************************************
1.1045 +!
1.1046 +!------------------------------------------------------------------------------
1.1047 +!
1.1048 + real function cosd(angle)
1.1049 + implicit none
1.1050 + real, intent(in) :: angle
1.1051 +
1.1052 + cosd = cos(angle * rad_per_deg)
1.1053 +
1.1054 + return
1.1055 + end function cosd
1.1056 +!
1.1057 +!------------------------------------------------------------------------------
1.1058 +!
1.1059 + real function sind(angle)
1.1060 + implicit none
1.1061 + real, intent(in) :: angle
1.1062 +
1.1063 + sind = sin(angle * rad_per_deg)
1.1064 +
1.1065 + return
1.1066 + end function sind
1.1067 +!
1.1068 +!------------------------------------------------------------------------------
1.1069 +!
1.1070 + real function atand(x)
1.1071 + implicit none
1.1072 + real, intent(in) :: x
1.1073 +
1.1074 + atand = atan(x) * deg_per_rad
1.1075 +
1.1076 + return
1.1077 + end function atand
1.1078 +!
1.1079 +!------------------------------------------------------------------------------
1.1080 +!
1.1081 + real function atan2d(x, y)
1.1082 + implicit none
1.1083 + real, intent(in) :: x
1.1084 + real, intent(in) :: y
1.1085 +
1.1086 + atan2d = atan2(x, y) * deg_per_rad
1.1087 +
1.1088 + return
1.1089 + end function atan2d
1.1090 +!
1.1091 +!------------------------------------------------------------------------------
1.1092 +!
1.1093 +
1.1094 +end module realizer_mod