forrest@0: module realizer_mod forrest@0: forrest@0: ! Realizer module written for the MidPtFM2d program (now called realizer) by forrest@0: ! Forrest Hoffman on Fri May 4 16:56:34 EDT 2007 to eliminate duplicated forrest@0: ! parameters and common blocks. Parts of the code were converted to forrest@0: ! Fortran-90 style and intrinsics. Compiled and tested with gfortran forrest@0: ! version 4.1.1 under Fedora Core 6 Linux. forrest@0: forrest@0: implicit none forrest@0: save forrest@0: forrest@0: ! Parameters forrest@0: integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real forrest@0: integer, parameter :: r4 = selected_real_kind( 6) ! 4 byte real forrest@0: integer, parameter :: nr = kind(1.0) ! native real forrest@0: integer, parameter :: i8 = selected_int_kind(13) ! 8 byte integer forrest@0: integer, parameter :: i4 = selected_int_kind( 6) ! 4 byte integer forrest@0: integer, parameter :: ni = kind(1) ! native integer forrest@0: integer(i4), parameter :: maxN=515 forrest@0: integer(i4), parameter :: maxsp=10 forrest@0: real(r8), parameter :: rad_per_deg=(3.14159265358979323846/180.0) forrest@0: real(r8), parameter :: deg_per_rad=(180.0/3.14159265358979323846) forrest@0: forrest@0: ! Global Variables forrest@0: integer(i4) :: numcells forrest@0: integer(i4) :: itiedsplit forrest@0: integer(i4) :: irowcols forrest@0: integer(i4) :: mtsplit forrest@0: integer(i4) :: mtresolved forrest@0: integer(i4) :: isupply forrest@0: integer(i4) :: ihabs forrest@0: integer(i4) :: maxgray forrest@0: real(r4) :: northeast forrest@0: real(r4) :: southeast forrest@0: real(r4) :: southwest forrest@0: real(r4) :: northwest forrest@0: forrest@0: ! Private Member Functions forrest@0: public dump_map forrest@0: public mpfm2d forrest@0: private StackSort ! unused? forrest@0: private sort ! unused? forrest@0: private SwapDown forrest@0: public Interchange forrest@0: public HeapSort forrest@0: private SD2 forrest@0: private I2 forrest@0: public HS2 forrest@0: private Gauss forrest@0: private cosd forrest@0: private sind forrest@0: private atand forrest@0: private atan2d forrest@0: forrest@0: contains forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine dump_map(P) forrest@0: implicit none forrest@0: integer(i4), intent(in) :: P(maxN*maxN,0:maxsp+3) forrest@0: ! local variables forrest@0: integer(i4) :: i, j ! loop indices forrest@0: integer(i4) :: temp(maxsp+3) forrest@0: character (len=40) :: fmt forrest@0: ! forrest@0: write(fmt,'("(i5, 1x, ",I0,"(i5, 1x), 1x, 3(i5, 1x))")') ihabs+1 forrest@0: do i=1,numcells forrest@0: do j = 0,ihabs+3 forrest@0: temp(j) = P(i,j) forrest@0: end do forrest@0: print fmt, i, (temp(j), j=0,ihabs+3) forrest@0: end do forrest@0: ! forrest@0: return forrest@0: end subroutine dump_map forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: real(r4) function f3(delta, x0, x1, x2, sigma) forrest@0: implicit none forrest@0: real(r4), intent(in) :: delta forrest@0: real(r4), intent(in) :: x0, x1, x2 forrest@0: real(r4), intent(in) :: sigma forrest@0: forrest@0: f3 = (x0+x1+x2) / 3.0 + delta * Gauss(0.,sigma) forrest@0: forrest@0: return forrest@0: end function f3 forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: real(r4) function f4(delta, x0, x1, x2, x3, sigma) forrest@0: implicit none forrest@0: real(r4), intent(in) :: delta forrest@0: real(r4), intent(in) :: x0, x1, x2, x3 forrest@0: real(r4), intent(in) :: sigma forrest@0: forrest@0: f4 = (x0+x1+x2+x3) / 4.0 + delta * Gauss(0.,sigma) forrest@0: forrest@0: return forrest@0: end function f4 forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine mpfm2d(X, maxlevel, sigma, H, addition, wrap, gradient, iz) forrest@0: implicit none forrest@0: real(r4), intent(inout) :: X(0:maxN,0:maxN) forrest@0: integer(i4), intent(in) :: maxlevel forrest@0: real(r4), intent(in) :: sigma forrest@0: real(r4), intent(in) :: H(maxN,maxsp) forrest@0: logical, intent(in) :: addition forrest@0: logical, intent(in) :: wrap forrest@0: logical, intent(in) :: gradient forrest@0: integer(i4), intent(in) :: iz forrest@0: ! local variables forrest@0: integer(i4) :: N forrest@0: real(r4) :: delta forrest@0: integer(i4) :: DD forrest@0: integer(i4) :: d forrest@0: integer(i4) :: is, ix, iy ! loop indices forrest@0: ! forrest@0: print *, "Generating fractal" forrest@0: N = 2**maxlevel forrest@0: delta = sigma forrest@0: X(0,0) = delta * Gauss(0.,sigma) forrest@0: X(0,N) = delta * Gauss(0.,sigma) forrest@0: X(N,0) = delta * Gauss(0.,sigma) forrest@0: X(N,N) = delta * Gauss(0.,sigma) forrest@0: ! forrest@0: if (gradient) then forrest@0: X(0,0) = southwest forrest@0: X(0,N) = southeast forrest@0: X(N,0) = northwest forrest@0: X(N,N) = northeast forrest@0: end if forrest@0: ! forrest@0: DD = N forrest@0: d = N / 2 forrest@0: ! forrest@0: do is = 1, maxlevel forrest@0: ! going from grid type I to type II forrest@0: delta = delta * 0.5**(0.5 * H(is,iz)) forrest@0: do ix = d, N-d, DD forrest@0: do iy = d, N-d, DD forrest@0: X(ix,iy) = f4(delta, X(ix+d,iy+d), X(ix+d,iy-d), & forrest@0: X(ix-d,iy+d), X(ix-d,iy-d), sigma) forrest@0: end do forrest@0: end do forrest@0: ! displace other points also if needed forrest@0: if (addition) then forrest@0: do ix = 0, N, DD forrest@0: do iy = 0, N, DD forrest@0: X(ix,iy) = X(ix,iy) + delta * Gauss(0.,sigma) forrest@0: end do forrest@0: end do forrest@0: end if forrest@0: ! going from grid type II to type I forrest@0: delta = delta * 0.5**(0.5*H(is,iz)) forrest@0: ! interpolate and offset boundary grid points forrest@0: do ix = d, N-d, DD forrest@0: X(ix,0) = f3(delta, X(ix+d,0), X(ix-d,0), X(ix,d), sigma) forrest@0: X(ix,N) = f3(delta, X(ix+d,N), X(ix-d,N), X(ix,N-d), sigma) forrest@0: X(0,ix) = f3(delta, X(0,ix+d), X(0,ix-d), X(d,ix), sigma) forrest@0: X(N,ix) = f3(delta, X(N,ix+d), X(N,ix-d), X(N-d,ix), sigma) forrest@0: if (wrap) then forrest@0: X(ix,N) = X(ix,0) forrest@0: X(N,ix) = X(0,ix) forrest@0: end if forrest@0: end do forrest@0: ! interpolate and offset inter grid points forrest@0: do ix = d, N-d, DD forrest@0: do iy = DD, N-d, DD forrest@0: X(ix,iy) = f4(delta, X(ix,iy+d), X(ix,iy-d), & forrest@0: X(ix+d,iy), X(ix-d,iy), sigma) forrest@0: end do forrest@0: end do forrest@0: do ix = DD, N-d, DD forrest@0: do iy = d, N-d, DD forrest@0: X(ix,iy) = f4(delta, X(ix,iy+d), X(ix,iy-d), & forrest@0: X(ix+d,iy), X(ix-d,iy), sigma) forrest@0: end do forrest@0: end do forrest@0: ! displace other points also if needed forrest@0: if (addition) then forrest@0: do ix = 0, N, DD forrest@0: do iy = 0, N, DD forrest@0: X(ix,iy) = X(ix,iy) + delta * Gauss(0.,sigma) forrest@0: end do forrest@0: end do forrest@0: do ix = d, N-d, DD forrest@0: do iy = d, N-d, DD forrest@0: X(ix,iy) = X(ix,iy) + delta * Gauss(0.,sigma) forrest@0: end do forrest@0: end do forrest@0: end if forrest@0: ! forrest@0: DD = DD / 2 forrest@0: d = d / 2 forrest@0: ! forrest@0: end do ! for all levels forrest@0: return forrest@0: end subroutine mpfm2d forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: real(r4) function Gauss(xm, xs) forrest@0: implicit none forrest@0: real(r4), intent(in) :: xm forrest@0: real(r4), intent(in) :: xs forrest@0: ! local variables forrest@0: real(r4) :: total forrest@0: real(r4) :: p forrest@0: integer(i4) :: i forrest@0: ! forrest@0: ! create a normal number by summing 12 uniform random numbers forrest@0: total = 0.0 forrest@0: do i = 1,12 forrest@0: call random_number(p) forrest@0: total = total + p forrest@0: end do forrest@0: Gauss = (total - 6.0) * xs + xm forrest@0: ! forrest@0: return forrest@0: end function Gauss forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine StackSort(n,arr) forrest@0: implicit none forrest@0: integer(i4), intent(in) :: n forrest@0: real(r4), intent(inout) :: arr(n) forrest@0: ! local variables forrest@0: integer(i4) :: i, j ! loop indices forrest@0: real(r4) :: a forrest@0: ! forrest@0: ! pick out each element in turn forrest@0: do j = 2,n forrest@0: a = arr(j) forrest@0: ! look for the place to insert it forrest@0: do i = j-1,1,-1 forrest@0: if (arr(i) <= a) go to 10 forrest@0: arr(i+1) = arr(i) forrest@0: end do forrest@0: i = 0 forrest@0: ! insert the number here forrest@0: 10 arr(i+1) = a forrest@0: end do forrest@0: return forrest@0: end subroutine StackSort forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: ! sort routine from Numerical Recipes forrest@0: subroutine sort (n, ra) forrest@0: implicit none forrest@0: integer(i4), intent(in) :: n forrest@0: real(r4), intent(inout) :: ra(n) forrest@0: ! local variables forrest@0: integer(i4) :: i, j, l forrest@0: integer(i4) :: ir forrest@0: real(r4) :: rra forrest@0: forrest@0: l = n/2 + 1 forrest@0: ir = n forrest@0: ! forrest@0: 10 continue forrest@0: if (l .gt. 1) then forrest@0: l = l-1 forrest@0: rra = ra(l) forrest@0: else forrest@0: rra = ra(ir) forrest@0: ra(ir) = ra(1) forrest@0: ir = ir-1 forrest@0: if (ir .eq. 1) then forrest@0: ra(1) = rra forrest@0: return forrest@0: end if forrest@0: end if forrest@0: i = l forrest@0: j = l+l forrest@0: ! forrest@0: 20 if (j .le. ir) then forrest@0: if (j .lt. ir) then forrest@0: if (ra(j) .lt. ra(j+1)) j = j+1 forrest@0: end if forrest@0: if (rra .lt. ra(j)) then forrest@0: ra(i) = ra(j) forrest@0: i = j forrest@0: j = j+j forrest@0: else forrest@0: j = ir+1 forrest@0: end if forrest@0: goto 20 forrest@0: end if forrest@0: ! forrest@0: ra(i) = rra forrest@0: goto 10 forrest@0: ! forrest@0: end subroutine sort forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine SwapDown(Heap, r, isize, istart, keyfield) forrest@0: implicit none forrest@0: integer(i4), intent(inout) :: Heap(maxN*maxN,0:maxsp+3) forrest@0: integer(i4), intent(inout) :: r forrest@0: integer(i4), intent(in) :: isize forrest@0: integer(i4), intent(in) :: istart forrest@0: integer(i4), intent(in) :: keyfield forrest@0: ! local variables forrest@0: integer(i4) :: child forrest@0: integer(i4) :: ipos forrest@0: logical :: done forrest@0: ! forrest@0: ! print *, 'entering SwapDown' forrest@0: ! print *, 'params are ', r, isize, istart, keyfield forrest@0: ipos = istart - 1 forrest@0: done = .false. forrest@0: child = 2 * r forrest@0: do while (.not. done .and. child .le. isize) forrest@0: ! find the largest child forrest@0: if (child .lt. isize .and. Heap(child+ipos,keyfield) .lt. & forrest@0: Heap(child+ipos+1,keyfield)) child = child + 1 forrest@0: ! interchange node and largest child if necessary forrest@0: ! and move down to the next subtree forrest@0: if (Heap(r+ipos,keyfield) .lt. Heap(child+ipos,keyfield)) then forrest@0: Call Interchange(Heap,r+ipos,child+ipos) forrest@0: r = child forrest@0: child = 2 * child forrest@0: else forrest@0: done = .true. forrest@0: end if forrest@0: ! forrest@0: end do forrest@0: ! forrest@0: return forrest@0: end subroutine SwapDown forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine Interchange(X,m,n) forrest@0: implicit none forrest@0: integer(i4), intent(inout) :: X(maxN*maxN,0:maxsp+3) forrest@0: integer(i4), intent(in) :: m forrest@0: integer(i4), intent(in) :: n forrest@0: ! local variables forrest@0: integer(i4) :: i forrest@0: integer(i4) :: temp forrest@0: ! forrest@0: do i = 0, ihabs+3 forrest@0: temp = X(m,i) forrest@0: X(m,i) = X(n,i) forrest@0: X(n,i) = temp forrest@0: end do forrest@0: ! forrest@0: return forrest@0: end subroutine Interchange forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine HeapSort(X,isize,istart,keyfield) forrest@0: implicit none forrest@0: integer(i4), intent(inout) :: X(maxN*maxN,0:maxsp+3) forrest@0: integer(i4), intent(in) :: isize forrest@0: integer(i4), intent(in) :: istart forrest@0: integer(i4), intent(in) :: keyfield forrest@0: ! local variables forrest@0: integer(i4) :: i ! loop indices forrest@0: integer(i4) :: r forrest@0: integer(i4) :: tempr forrest@0: ! forrest@0: ! print *, "isize in HeapSort is ", isize forrest@0: ! print *, "istart in HeapSort is ", istart forrest@0: ! print *, "keyfield in HeapSort is ", keyfield forrest@0: ! heapify forrest@0: do r = isize/2, 1, -1 forrest@0: tempr = r forrest@0: !Call SwapDown(X,r,isize,istart,keyfield) forrest@0: Call SwapDown(X,tempr,isize,istart,keyfield) forrest@0: !r = tempr forrest@0: end do forrest@0: ! forrest@0: ! print *, "P after heapify is " forrest@0: ! do i=istart,isize forrest@0: ! print *, i, X(i,keyfield) forrest@0: ! end do forrest@0: ! forrest@0: ! repeatedly put largest item in root at end of list, forrest@0: ! prune it from the tree forrest@0: ! and apply SwapDown to the rest of the tree forrest@0: ! forrest@0: do i = isize, 2, -1 forrest@0: Call Interchange(X,istart,(i+istart-1)) forrest@0: ! r = isize forrest@0: r = 1 forrest@0: Call SwapDown(X,r,i-1,istart,keyfield) forrest@0: ! forrest@0: ! print *, "P after Swapdown is " forrest@0: ! do j=istart,isize forrest@0: ! print *, j, X(j,keyfield) forrest@0: ! end do forrest@0: end do forrest@0: ! forrest@0: return forrest@0: end subroutine HeapSort forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine XPMMap(N,II) forrest@0: ! forrest@0: ! output routine for xpm X-window display tool forrest@0: ! forrest@0: ! forrest@0: implicit none forrest@0: integer(i4), intent(in) :: N forrest@0: integer(i4), intent(in) :: II(0:maxN,0:maxN) forrest@0: ! local variables forrest@0: integer(i4) :: itemp(0:maxN) forrest@0: integer(i4) :: i, j ! loop indices forrest@0: character (len=40) :: fmt forrest@0: ! forrest@0: ! forrest@0: ! output map forrest@0: print *, "writing final map to xpm format file landscape.xpm" forrest@0: ! forrest@0: open (20,file='landscape.xpm',status='unknown') forrest@0: ! forrest@0: write (20, 100) forrest@0: 100 format ('/* XPM */',/'static char *realizer[]', & forrest@0: ' = {',/'/* cols rows ncolors chars_per_pixel', & forrest@0: ' */') forrest@0: write (20, 102) N+1, N+1 forrest@0: 102 format ('"',i6,1x,i6,1x,'11 1 0 0",',/'/* colors */') forrest@0: write (20, 104) forrest@0: 104 format ('"0 c white m white', & forrest@0: '",',/'"1 s Decid c BurlyWood m black', & forrest@0: '",',/'"2 s Mixed c DarkGoldenrod3 m gray40', & forrest@0: '",',/'"3 s Everg c MistyRose3 m gray50', & forrest@0: '",',/'"4 s Trans c LightGoldenrod4 m gray60', & forrest@0: '",',/'"5 s Urban c DarkGoldenrod1 m gray70', & forrest@0: '",',/'"6 s Water c gold m gray80', & forrest@0: '",',/'"7 s Maint c LightGoldenrod3 m white', & forrest@0: '",',/'"8 s Cropl c SaddleBrown m lightgray', & forrest@0: '",',/'"9 s Lawng c RosyBrown m gray', & forrest@0: '",',/'"* c chocolate m black",',/'/* pixels */') forrest@0: ! forrest@0: write(fmt,'("(",I0,"(i1)$)")') N+1 forrest@0: do i= 0, N forrest@0: do j = 0, N forrest@0: ! itemp(j) = II(i,j) forrest@0: itemp(j) = II(j,i) forrest@0: end do forrest@0: write (20,'(1A$)') '"' forrest@0: ! write (20,108) (itemp(j), j=0,N) forrest@0: write (20,fmt) (itemp(j), j=0,N) forrest@0: write (20,'(2A)') '",' forrest@0: !108 format (' "',(i1),'",') forrest@0: end do forrest@0: ! forrest@0: ! forrest@0: write (20,112) forrest@0: 112 format ('};') forrest@0: ! forrest@0: return forrest@0: end subroutine XPMMap forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine InputMap(mapfile,P,N,k,invert) forrest@0: ! forrest@0: ! input map from a data file forrest@0: ! forrest@0: implicit none forrest@0: character (len=60), intent(in) :: mapfile forrest@0: integer(i4), intent(inout) :: P(maxN*maxN,0:maxsp+3) forrest@0: integer(i4), intent(in) :: N forrest@0: integer(i4), intent(in) :: k forrest@0: logical, intent(in) :: invert forrest@0: ! local variables forrest@0: integer(i4) :: i, j, iv ! loop indices forrest@0: real(r4) :: X(0:maxN,0:maxN) forrest@0: real(r4) :: xmin forrest@0: real(r4) :: xmax forrest@0: real(r4) :: xrange forrest@0: real(r4) :: temp forrest@0: ! forrest@0: open (14, file=mapfile, status='old') forrest@0: ! forrest@0: xmin = 1e30 forrest@0: xmax = -1e30 forrest@0: do j = 0, N forrest@0: read (14,*) ( X(i,j), i = 0, N) forrest@0: ! read (14,118) X(i,j) forrest@0: !118 format (<(N+1)*(N+1)>i5, 1x) forrest@0: do i = 0, N forrest@0: if (X(i,j) .lt. xmin) xmin = X(i,j) forrest@0: if (X(i,j) .gt. xmax) xmax = X(i,j) forrest@0: end do forrest@0: end do forrest@0: ! forrest@0: ! invert if necessary and scale between 0 and 32767 forrest@0: xrange = xmax - xmin forrest@0: do j = 0, N forrest@0: do i = 0, N forrest@0: if (invert) then forrest@0: temp = ((-X(i,j) + xmax) / xrange) * 32767 forrest@0: else forrest@0: temp = ((X(i,j) - xmin) / xrange) * 32767 forrest@0: end if forrest@0: ! forrest@0: iv = (i*(N+1))+j+1 forrest@0: P(iv,k) = nint(temp) forrest@0: P(iv,ihabs+1) = i forrest@0: P(iv,ihabs+2) = j forrest@0: end do forrest@0: end do forrest@0: ! forrest@0: close(14) forrest@0: rewind(14) forrest@0: ! forrest@0: return forrest@0: end subroutine InputMap forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine SD2(Heap,r,isize,istart) forrest@0: implicit none forrest@0: integer(i4), intent(inout) :: Heap(maxN*maxN) forrest@0: integer(i4), intent(inout) :: r forrest@0: integer(i4), intent(in) :: isize forrest@0: integer(i4), intent(in) :: istart forrest@0: ! local variables forrest@0: integer(i4) :: ipos forrest@0: integer(i4) :: child forrest@0: logical :: done forrest@0: ! forrest@0: ! print *, 'entering SwapDown' forrest@0: ! print *, 'params are ', r, isize, istart forrest@0: ipos = istart - 1 forrest@0: done = .false. forrest@0: child = 2 * r forrest@0: do while (.not. done .and. child .le. isize) forrest@0: ! find the largest child forrest@0: if (child .lt. isize .and. Heap(child+ipos) .lt. & forrest@0: Heap(child+ipos+1)) child = child + 1 forrest@0: ! interchange node and largest child if necessary forrest@0: ! and move down to the next subtree forrest@0: if (Heap(r+ipos) .lt. Heap(child+ipos)) then forrest@0: Call I2(Heap,r+ipos,child+ipos) forrest@0: r = child forrest@0: child = 2 * child forrest@0: else forrest@0: done = .true. forrest@0: end if forrest@0: ! forrest@0: end do forrest@0: ! forrest@0: return forrest@0: end subroutine SD2 forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine I2(X,m,n) forrest@0: implicit none forrest@0: integer(i4), intent(inout) :: X(maxN*maxN) forrest@0: integer(i4), intent(in) :: m forrest@0: integer(i4), intent(in) :: n forrest@0: ! local variables forrest@0: integer(i4) :: temp forrest@0: ! forrest@0: temp = X(m) forrest@0: X(m) = X(n) forrest@0: X(n) = temp forrest@0: ! forrest@0: return forrest@0: end subroutine I2 forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine HS2(X,isize,istart) forrest@0: implicit none forrest@0: integer(i4), intent(inout) :: X(maxN*maxN) forrest@0: integer(i4), intent(in) :: isize forrest@0: integer(i4), intent(in) :: istart forrest@0: ! local variables forrest@0: integer(i4) :: i, j ! loop indices forrest@0: integer(i4) :: r forrest@0: integer(i4) :: tempr forrest@0: ! forrest@0: ! print *, "isize in HeapSort is ", isize forrest@0: ! print *, "istart in HeapSort is ", istart forrest@0: ! heapify forrest@0: do r = isize/2, 1, -1 forrest@0: tempr = r forrest@0: Call SD2(X,tempr,isize,istart) forrest@0: !Call SD2(X,r,isize,istart) forrest@0: !r = tempr forrest@0: end do forrest@0: ! forrest@0: ! print*, "P after heapify is " forrest@0: ! do i=istart,isize forrest@0: ! print *, i, X(i) forrest@0: ! end do forrest@0: ! forrest@0: ! repeatedly put largest item in root at end of list, forrest@0: ! prune it from the tree forrest@0: ! and apply SwapDown to the rest of the tree forrest@0: ! forrest@0: do i = isize, 2, -1 forrest@0: Call I2(X,istart,(i+istart-1)) forrest@0: ! r = isize forrest@0: r = 1 forrest@0: Call SD2(X,r,i-1,istart) forrest@0: ! forrest@0: ! print*, "P after Swapdown is " forrest@0: ! do j=istart,isize forrest@0: ! print *, j, X(j) forrest@0: ! end do forrest@0: forrest@0: end do forrest@0: ! forrest@0: return forrest@0: end subroutine HS2 forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine XPMTies(P) forrest@0: ! forrest@0: ! output routine for xpm X-window display tool forrest@0: ! forrest@0: ! write xpm format map with grayshades forrest@0: ! forrest@0: implicit none forrest@0: integer(i4), intent(in) :: P(maxN*maxN,0:maxsp+3) forrest@0: ! local variables forrest@0: integer(i4) :: itemp(0:maxN,0:maxN) forrest@0: integer(i4) :: i, j, m ! loop indices forrest@0: real(r4) :: scalelog forrest@0: character (len=40) :: fmt forrest@0: ! forrest@0: open (40,file='ties.xpm',status='unknown') forrest@0: ! forrest@0: print 99 forrest@0: 99 format ('XPM map file ties.xpm being written') forrest@0: write (40, 100) forrest@0: 100 format ('/* XPM */',/'static char *tiesmap[]', & forrest@0: ' = {',/'/* cols rows ncolors chars_per_pixel', & forrest@0: ' */') forrest@0: write (40, 102) irowcols+1, irowcols+1 forrest@0: 102 format (' "',i6,1x,i6,1x,'108 3 0 0",',/'/* colors */') forrest@0: write (40, 104) forrest@0: 104 format ('" 0 s 0.0prob c white g4 white g white m white', & forrest@0: '",',/'" 1 s 0.01prob c gray99 g4 gray99 g gray99 m gray99', & forrest@0: '",',/'" 2 s 0.02prob c gray98 g4 gray98 g gray98 m gray98', & forrest@0: '",',/'" 3 s 0.03prob c gray97 g4 gray97 g gray97 m gray97', & forrest@0: '",',/'" 4 s 0.04prob c gray96 g4 gray96 g gray96 m gray96', & forrest@0: '",',/'" 5 s 0.05prob c gray95 g4 gray95 g gray95 m gray95', & forrest@0: '",',/'" 6 s 0.06prob c gray94 g4 gray94 g gray94 m gray94', & forrest@0: '",',/'" 7 s 0.07prob c gray93 g4 gray93 g gray93 m gray93', & forrest@0: '",',/'" 8 s 0.08prob c gray92 g4 gray92 g gray92 m gray92', & forrest@0: '",',/'" 9 s 0.09prob c gray91 g4 gray91 g gray91 m gray91', & forrest@0: '",',/'" 10 s 0.10prob c gray90 g4 gray90 g gray90 m gray90', & forrest@0: '",',/'" 11 s 0.11prob c gray89 g4 gray89 g gray89 m gray89', & forrest@0: '",',/'" 12 s 0.12prob c gray88 g4 gray88 g gray88 m gray88', & forrest@0: '",',/'" 13 s 0.13prob c gray87 g4 gray87 g gray87 m gray87', & forrest@0: '",',/'" 14 s 0.14prob c gray86 g4 gray86 g gray86 m gray86', & forrest@0: '",',/'" 15 s 0.15prob c gray85 g4 gray85 g gray85 m gray85', & forrest@0: '",',/'" 16 s 0.16prob c gray84 g4 gray84 g gray84 m gray84', & forrest@0: '",',/'" 17 s 0.17prob c gray83 g4 gray83 g gray83 m gray83', & forrest@0: '",',/'" 18 s 0.18prob c gray82 g4 gray82 g gray82 m gray82', & forrest@0: '",',/'" 19 s 0.19prob c gray81 g4 gray81 g gray81 m gray81', & forrest@0: '",',/'" 20 s 0.20prob c gray80 g4 gray80 g gray80 m gray80', & forrest@0: '",',/'" 21 s 0.21prob c gray79 g4 gray79 g gray79 m gray79', & forrest@0: '",',/'" 22 s 0.22prob c gray78 g4 gray78 g gray78 m gray78', & forrest@0: '",',/'" 23 s 0.23prob c gray77 g4 gray77 g gray77 m gray77', & forrest@0: '",',/'" 24 s 0.24prob c gray76 g4 gray76 g gray76 m gray76', & forrest@0: '",',/'" 25 s 0.25prob c gray75 g4 gray75 g gray75 m gray75', & forrest@0: '",',/'" 26 s 0.26prob c gray74 g4 gray74 g gray74 m gray74', & forrest@0: '",',/'" 27 s 0.27prob c gray73 g4 gray73 g gray73 m gray73', & forrest@0: '",',/'" 28 s 0.28prob c gray72 g4 gray72 g gray72 m gray72', & forrest@0: '",',/'" 29 s 0.29prob c gray71 g4 gray71 g gray71 m gray71', & forrest@0: '",',/'" 30 s 0.30prob c gray70 g4 gray70 g gray70 m gray70', & forrest@0: '",',/'" 31 s 0.31prob c gray69 g4 gray69 g gray69 m gray69', & forrest@0: '",',/'" 32 s 0.32prob c gray68 g4 gray68 g gray68 m gray68', & forrest@0: '",',/'" 33 s 0.33prob c gray67 g4 gray67 g gray67 m gray67', & forrest@0: '",',/'" 34 s 0.34prob c gray66 g4 gray66 g gray66 m gray66', & forrest@0: '",',/'" 35 s 0.35prob c gray65 g4 gray65 g gray65 m gray65', & forrest@0: '",',/'" 36 s 0.36prob c gray64 g4 gray64 g gray64 m gray64', & forrest@0: '",',/'" 37 s 0.37prob c gray63 g4 gray63 g gray63 m gray63', & forrest@0: '",',/'" 38 s 0.38prob c gray62 g4 gray62 g gray62 m gray62', & forrest@0: '",',/'" 39 s 0.39prob c gray61 g4 gray61 g gray61 m gray61', & forrest@0: '",',/'" 40 s 0.40prob c gray60 g4 gray60 g gray60 m gray60', & forrest@0: '",',/'" 41 s 0.41prob c gray59 g4 gray59 g gray59 m gray59', & forrest@0: '",',/'" 42 s 0.42prob c gray58 g4 gray58 g gray58 m gray58', & forrest@0: '",',/'" 43 s 0.43prob c gray57 g4 gray57 g gray57 m gray57', & forrest@0: '",',/'" 44 s 0.44prob c gray56 g4 gray56 g gray56 m gray56', & forrest@0: '",',/'" 45 s 0.45prob c gray55 g4 gray55 g gray55 m gray55', & forrest@0: '",',/'" 46 s 0.46prob c gray54 g4 gray54 g gray54 m gray54', & forrest@0: '",',/'" 47 s 0.47prob c gray53 g4 gray53 g gray53 m gray53', & forrest@0: '",',/'" 48 s 0.48prob c gray52 g4 gray52 g gray52 m gray52', & forrest@0: '",',/'" 49 s 0.49prob c gray51 g4 gray51 g gray51 m gray51', & forrest@0: '",',/'" 50 s 0.50prob c gray50 g4 gray50 g gray50 m gray50', & forrest@0: '",',/'" 51 s 0.51prob c gray49 g4 gray49 g gray49 m gray49', & forrest@0: '",',/'" 52 s 0.52prob c gray48 g4 gray48 g gray48 m gray48', & forrest@0: '",',/'" 53 s 0.53prob c gray47 g4 gray47 g gray47 m gray47', & forrest@0: '",',/'" 54 s 0.54prob c gray46 g4 gray46 g gray46 m gray46', & forrest@0: '",',/'" 55 s 0.55prob c gray45 g4 gray45 g gray45 m gray45', & forrest@0: '",',/'" 56 s 0.56prob c gray44 g4 gray44 g gray44 m gray44', & forrest@0: '",',/'" 57 s 0.57prob c gray43 g4 gray43 g gray43 m gray43', & forrest@0: '",',/'" 58 s 0.58prob c gray42 g4 gray42 g gray42 m gray42', & forrest@0: '",',/'" 59 s 0.59prob c gray41 g4 gray41 g gray41 m gray41', & forrest@0: '",',/'" 60 s 0.60prob c gray40 g4 gray40 g gray40 m gray40', & forrest@0: '",',/'" 61 s 0.61prob c gray39 g4 gray39 g gray39 m gray39', & forrest@0: '",',/'" 62 s 0.62prob c gray38 g4 gray38 g gray38 m gray38', & forrest@0: '",',/'" 63 s 0.63prob c gray37 g4 gray37 g gray37 m gray37', & forrest@0: '",',/'" 64 s 0.64prob c gray36 g4 gray36 g gray36 m gray36', & forrest@0: '",',/'" 65 s 0.65prob c gray35 g4 gray35 g gray35 m gray35', & forrest@0: '",',/'" 66 s 0.66prob c gray34 g4 gray34 g gray34 m gray34', & forrest@0: '",',/'" 67 s 0.67prob c gray33 g4 gray33 g gray33 m gray33', & forrest@0: '",',/'" 68 s 0.68prob c gray32 g4 gray32 g gray32 m gray32', & forrest@0: '",',/'" 69 s 0.69prob c gray31 g4 gray31 g gray31 m gray31', & forrest@0: '",',/'" 70 s 0.70prob c gray30 g4 gray30 g gray30 m gray30', & forrest@0: '",',/'" 71 s 0.71prob c gray29 g4 gray29 g gray29 m gray29', & forrest@0: '",',/'" 72 s 0.72prob c gray28 g4 gray28 g gray28 m gray28', & forrest@0: '",',/'" 73 s 0.73prob c gray27 g4 gray27 g gray27 m gray27', & forrest@0: '",',/'" 74 s 0.74prob c gray26 g4 gray26 g gray26 m gray26', & forrest@0: '",',/'" 75 s 0.75prob c gray25 g4 gray25 g gray25 m gray25', & forrest@0: '",',/'" 76 s 0.76prob c gray24 g4 gray24 g gray24 m gray24', & forrest@0: '",',/'" 77 s 0.77prob c gray23 g4 gray23 g gray23 m gray23', & forrest@0: '",',/'" 78 s 0.78prob c gray22 g4 gray22 g gray22 m gray22', & forrest@0: '",',/'" 79 s 0.79prob c gray21 g4 gray21 g gray21 m gray21', & forrest@0: '",',/'" 80 s 0.80prob c gray20 g4 gray20 g gray20 m gray20', & forrest@0: '",',/'" 81 s 0.81prob c gray19 g4 gray19 g gray19 m gray19', & forrest@0: '",',/'" 82 s 0.82prob c gray18 g4 gray18 g gray18 m gray18', & forrest@0: '",',/'" 83 s 0.83prob c gray17 g4 gray17 g gray17 m gray17', & forrest@0: '",',/'" 84 s 0.84prob c gray16 g4 gray16 g gray16 m gray16', & forrest@0: '",',/'" 85 s 0.85prob c gray15 g4 gray15 g gray15 m gray15', & forrest@0: '",',/'" 86 s 0.86prob c gray14 g4 gray14 g gray14 m gray14', & forrest@0: '",',/'" 87 s 0.87prob c gray13 g4 gray13 g gray13 m gray13', & forrest@0: '",',/'" 88 s 0.88prob c gray12 g4 gray12 g gray12 m gray12', & forrest@0: '",',/'" 89 s 0.89prob c gray11 g4 gray11 g gray11 m gray11', & forrest@0: '",',/'" 90 s 0.90prob c gray10 g4 gray10 g gray10 m gray10', & forrest@0: '",',/'" 91 s 0.91prob c gray9 g4 gray9 g gray9 m gray9', & forrest@0: '",',/'" 92 s 0.92prob c gray8 g4 gray8 g gray8 m gray8', & forrest@0: '",',/'" 93 s 0.93prob c gray7 g4 gray7 g gray7 m gray7', & forrest@0: '",',/'" 94 s 0.94prob c gray6 g4 gray6 g gray6 m gray6', & forrest@0: '",',/'" 95 s 0.95prob c gray5 g4 gray5 g gray5 m gray5', & forrest@0: '",',/'" 96 s 0.96prob c gray4 g4 gray4 g gray4 m gray4', & forrest@0: '",',/'" 97 s 0.97prob c gray3 g4 gray3 g gray3 m gray3', & forrest@0: '",',/'" 98 s 0.98prob c gray2 g4 gray2 g gray2 m gray2', & forrest@0: '",',/'" 99 s 0.99prob c gray1 g4 gray1 g gray1 m gray1', & forrest@0: '",',/'"100 s 1.00prob c black g4 black g black m black', & forrest@0: '",',/'"101 c white m white s NotFlammable', & forrest@0: '",',/'"102 c GreenYellow m gray70 s LP0', & forrest@0: '",',/'"103 c LawnGreen m gray60 s LP1', & forrest@0: '",',/'"104 c LimeGreen m gray50 s LP2', & forrest@0: '",',/'"105 c ForestGreen m gray40 s LP3', & forrest@0: '",',/'"106 c yellow m gray80 s NF', & forrest@0: '",',/'"107 c blue m black s water ",', & forrest@0: /'/* pixels */') forrest@0: ! forrest@0: ! calculate a log scale for gray levels forrest@0: ! use this scale log for uniform scaling across realizations forrest@0: ! scalelog = 10**(log10(100.)/32767.) forrest@0: ! scalelog for simple summed priority forrest@0: ! scale max gray to ihabs-isupply-1 full-on categories forrest@0: ! (= ihabs-isupply-1 x 32767) forrest@0: !c scalelog = 10**(log10(100.)/(32767. * (ihabs-isupply-1))) forrest@0: ! use this scalelog for maximum priority resolution w/in one map forrest@0: scalelog = 10**(log10(100.) / real(maxgray,kind(scalelog))) forrest@0: ! print *, scalelog forrest@0: ! forrest@0: ! start at bottom of list (highest priority ties) forrest@0: do m = numcells, itiedsplit+1, -1 ! over all tied forrest@0: ! itemp(P(m,ihabs+1),P(m,ihabs+2)) = forrest@0: ! & (real(P(m,ihabs+3),kind(itemp))/32767)*100 forrest@0: ! forrest@0: ! print *, scalelog**P(m,ihabs+3) forrest@0: ! print *, P(m,ihabs+1),P(m,ihabs+2),P(m,ihabs+3) forrest@0: itemp(P(m,ihabs+1),P(m,ihabs+2)) = scalelog**P(m,ihabs+3) forrest@0: ! itemp(P(m,ihabs+1),P(m,ihabs+2)) = 100 forrest@0: ! all priorities greater than color index 100 set black forrest@0: if (itemp(P(m,ihabs+1),P(m,ihabs+2)) .gt. 100) & forrest@0: itemp(P(m,ihabs+1),P(m,ihabs+2)) = 100 forrest@0: end do forrest@0: ! forrest@0: ! color used blanks yellow forrest@0: do m = mtsplit,mtresolved+1,-1 forrest@0: itemp(P(m,ihabs+1),P(m,ihabs+2)) = 106 forrest@0: end do forrest@0: ! forrest@0: write (fmt,'("(",I0,"(i3)$)")') irowcols+1 forrest@0: do i = 0, irowcols forrest@0: write (40,'(1A$)') '"' forrest@0: !write (40,108) (itemp(j,i), j=0, irowcols) forrest@0: write (40,fmt) (itemp(j,i), j=0, irowcols) forrest@0: write (40,'(2A)') '",' forrest@0: !108 format ('"',(i3),'",') forrest@0: end do forrest@0: ! forrest@0: write (40,112) forrest@0: 112 format ('};') forrest@0: close (40) forrest@0: ! forrest@0: return forrest@0: end subroutine XPMTies forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: subroutine XPMRelief(P,iz,iflood) forrest@0: ! forrest@0: ! output routine for xpm X-window display tool forrest@0: ! forrest@0: ! write xpm format aspect map of painted fractal forrest@0: ! spatial probability surfaces with grayshades forrest@0: ! forrest@0: implicit none forrest@0: integer(i4), intent(in) :: P(maxN*maxN,0:maxsp+3) forrest@0: integer(i4), intent(in) :: iz forrest@0: integer(i4), intent(in) :: iflood forrest@0: ! local variables forrest@0: integer(i4) :: itemp(0:maxN,0:maxN) forrest@0: integer(i4) :: idispvec(0:maxN) forrest@0: integer(i4) :: i, j, k, l, m ! loop indices forrest@0: real(r4) :: alt ! altitude forrest@0: real(r4) :: az ! azimuth forrest@0: real(r4) :: x, y forrest@0: real(r4) :: slope ! slope forrest@0: real(r4) :: aspect ! aspect forrest@0: real(r4) :: c forrest@0: character (len=15) :: outfile(0:maxsp+1) forrest@0: character (len=40) :: fmt forrest@0: ! forrest@0: ! set filenames for the first 10 categories forrest@0: outfile(0) = 'probmap0.xpm ' forrest@0: outfile(1) = 'probmap1.xpm ' forrest@0: outfile(2) = 'probmap2.xpm ' forrest@0: outfile(3) = 'probmap3.xpm ' forrest@0: outfile(4) = 'probmap4.xpm ' forrest@0: outfile(5) = 'probmap5.xpm ' forrest@0: outfile(6) = 'probmap6.xpm ' forrest@0: outfile(7) = 'probmap7.xpm ' forrest@0: outfile(8) = 'probmap8.xpm ' forrest@0: outfile(9) = 'probmap9.xpm ' forrest@0: outfile(10) = 'probmap10.xpm ' forrest@0: ! forrest@0: ! forrest@0: print 99, outfile(iz) forrest@0: 99 format ('XPM map file ', a15 'being written') forrest@0: ! forrest@0: ! open the right file for xpm output forrest@0: open (50,file=outfile(iz),status='unknown') forrest@0: ! forrest@0: write (50, 100) forrest@0: 100 format ('/* XPM */',/'static char *fracasp[]', & forrest@0: ' = {',/'/* cols rows ncolors chars_per_pixel', ' */') forrest@0: write (50, 102) irowcols+1, irowcols+1 forrest@0: 102 format ('"',i6,1x,i6,1x,'108 3 0 0",',/'/* colors */') forrest@0: write (50, 104) forrest@0: 104 format ('" 0 s 0.0prob c white g4 white g white m white', & forrest@0: '",',/'" 1 s 0.01prob c gray99 g4 gray99 g gray99 m gray99', & forrest@0: '",',/'" 2 s 0.02prob c gray98 g4 gray98 g gray98 m gray98', & forrest@0: '",',/'" 3 s 0.03prob c gray97 g4 gray97 g gray97 m gray97', & forrest@0: '",',/'" 4 s 0.04prob c gray96 g4 gray96 g gray96 m gray96', & forrest@0: '",',/'" 5 s 0.05prob c gray95 g4 gray95 g gray95 m gray95', & forrest@0: '",',/'" 6 s 0.06prob c gray94 g4 gray94 g gray94 m gray94', & forrest@0: '",',/'" 7 s 0.07prob c gray93 g4 gray93 g gray93 m gray93', & forrest@0: '",',/'" 8 s 0.08prob c gray92 g4 gray92 g gray92 m gray92', & forrest@0: '",',/'" 9 s 0.09prob c gray91 g4 gray91 g gray91 m gray91', & forrest@0: '",',/'" 10 s 0.10prob c gray90 g4 gray90 g gray90 m gray90', & forrest@0: '",',/'" 11 s 0.11prob c gray89 g4 gray89 g gray89 m gray89', & forrest@0: '",',/'" 12 s 0.12prob c gray88 g4 gray88 g gray88 m gray88', & forrest@0: '",',/'" 13 s 0.13prob c gray87 g4 gray87 g gray87 m gray87', & forrest@0: '",',/'" 14 s 0.14prob c gray86 g4 gray86 g gray86 m gray86', & forrest@0: '",',/'" 15 s 0.15prob c gray85 g4 gray85 g gray85 m gray85', & forrest@0: '",',/'" 16 s 0.16prob c gray84 g4 gray84 g gray84 m gray84', & forrest@0: '",',/'" 17 s 0.17prob c gray83 g4 gray83 g gray83 m gray83', & forrest@0: '",',/'" 18 s 0.18prob c gray82 g4 gray82 g gray82 m gray82', & forrest@0: '",',/'" 19 s 0.19prob c gray81 g4 gray81 g gray81 m gray81', & forrest@0: '",',/'" 20 s 0.20prob c gray80 g4 gray80 g gray80 m gray80', & forrest@0: '",',/'" 21 s 0.21prob c gray79 g4 gray79 g gray79 m gray79', & forrest@0: '",',/'" 22 s 0.22prob c gray78 g4 gray78 g gray78 m gray78', & forrest@0: '",',/'" 23 s 0.23prob c gray77 g4 gray77 g gray77 m gray77', & forrest@0: '",',/'" 24 s 0.24prob c gray76 g4 gray76 g gray76 m gray76', & forrest@0: '",',/'" 25 s 0.25prob c gray75 g4 gray75 g gray75 m gray75', & forrest@0: '",',/'" 26 s 0.26prob c gray74 g4 gray74 g gray74 m gray74', & forrest@0: '",',/'" 27 s 0.27prob c gray73 g4 gray73 g gray73 m gray73', & forrest@0: '",',/'" 28 s 0.28prob c gray72 g4 gray72 g gray72 m gray72', & forrest@0: '",',/'" 29 s 0.29prob c gray71 g4 gray71 g gray71 m gray71', & forrest@0: '",',/'" 30 s 0.30prob c gray70 g4 gray70 g gray70 m gray70', & forrest@0: '",',/'" 31 s 0.31prob c gray69 g4 gray69 g gray69 m gray69', & forrest@0: '",',/'" 32 s 0.32prob c gray68 g4 gray68 g gray68 m gray68', & forrest@0: '",',/'" 33 s 0.33prob c gray67 g4 gray67 g gray67 m gray67', & forrest@0: '",',/'" 34 s 0.34prob c gray66 g4 gray66 g gray66 m gray66', & forrest@0: '",',/'" 35 s 0.35prob c gray65 g4 gray65 g gray65 m gray65', & forrest@0: '",',/'" 36 s 0.36prob c gray64 g4 gray64 g gray64 m gray64', & forrest@0: '",',/'" 37 s 0.37prob c gray63 g4 gray63 g gray63 m gray63', & forrest@0: '",',/'" 38 s 0.38prob c gray62 g4 gray62 g gray62 m gray62', & forrest@0: '",',/'" 39 s 0.39prob c gray61 g4 gray61 g gray61 m gray61', & forrest@0: '",',/'" 40 s 0.40prob c gray60 g4 gray60 g gray60 m gray60', & forrest@0: '",',/'" 41 s 0.41prob c gray59 g4 gray59 g gray59 m gray59', & forrest@0: '",',/'" 42 s 0.42prob c gray58 g4 gray58 g gray58 m gray58', & forrest@0: '",',/'" 43 s 0.43prob c gray57 g4 gray57 g gray57 m gray57', & forrest@0: '",',/'" 44 s 0.44prob c gray56 g4 gray56 g gray56 m gray56', & forrest@0: '",',/'" 45 s 0.45prob c gray55 g4 gray55 g gray55 m gray55', & forrest@0: '",',/'" 46 s 0.46prob c gray54 g4 gray54 g gray54 m gray54', & forrest@0: '",',/'" 47 s 0.47prob c gray53 g4 gray53 g gray53 m gray53', & forrest@0: '",',/'" 48 s 0.48prob c gray52 g4 gray52 g gray52 m gray52', & forrest@0: '",',/'" 49 s 0.49prob c gray51 g4 gray51 g gray51 m gray51', & forrest@0: '",',/'" 50 s 0.50prob c gray50 g4 gray50 g gray50 m gray50', & forrest@0: '",',/'" 51 s 0.51prob c gray49 g4 gray49 g gray49 m gray49', & forrest@0: '",',/'" 52 s 0.52prob c gray48 g4 gray48 g gray48 m gray48', & forrest@0: '",',/'" 53 s 0.53prob c gray47 g4 gray47 g gray47 m gray47', & forrest@0: '",',/'" 54 s 0.54prob c gray46 g4 gray46 g gray46 m gray46', & forrest@0: '",',/'" 55 s 0.55prob c gray45 g4 gray45 g gray45 m gray45', & forrest@0: '",',/'" 56 s 0.56prob c gray44 g4 gray44 g gray44 m gray44', & forrest@0: '",',/'" 57 s 0.57prob c gray43 g4 gray43 g gray43 m gray43', & forrest@0: '",',/'" 58 s 0.58prob c gray42 g4 gray42 g gray42 m gray42', & forrest@0: '",',/'" 59 s 0.59prob c gray41 g4 gray41 g gray41 m gray41', & forrest@0: '",',/'" 60 s 0.60prob c gray40 g4 gray40 g gray40 m gray40', & forrest@0: '",',/'" 61 s 0.61prob c gray39 g4 gray39 g gray39 m gray39', & forrest@0: '",',/'" 62 s 0.62prob c gray38 g4 gray38 g gray38 m gray38', & forrest@0: '",',/'" 63 s 0.63prob c gray37 g4 gray37 g gray37 m gray37', & forrest@0: '",',/'" 64 s 0.64prob c gray36 g4 gray36 g gray36 m gray36', & forrest@0: '",',/'" 65 s 0.65prob c gray35 g4 gray35 g gray35 m gray35', & forrest@0: '",',/'" 66 s 0.66prob c gray34 g4 gray34 g gray34 m gray34', & forrest@0: '",',/'" 67 s 0.67prob c gray33 g4 gray33 g gray33 m gray33', & forrest@0: '",',/'" 68 s 0.68prob c gray32 g4 gray32 g gray32 m gray32', & forrest@0: '",',/'" 69 s 0.69prob c gray31 g4 gray31 g gray31 m gray31', & forrest@0: '",',/'" 70 s 0.70prob c gray30 g4 gray30 g gray30 m gray30', & forrest@0: '",',/'" 71 s 0.71prob c gray29 g4 gray29 g gray29 m gray29', & forrest@0: '",',/'" 72 s 0.72prob c gray28 g4 gray28 g gray28 m gray28', & forrest@0: '",',/'" 73 s 0.73prob c gray27 g4 gray27 g gray27 m gray27', & forrest@0: '",',/'" 74 s 0.74prob c gray26 g4 gray26 g gray26 m gray26', & forrest@0: '",',/'" 75 s 0.75prob c gray25 g4 gray25 g gray25 m gray25', & forrest@0: '",',/'" 76 s 0.76prob c gray24 g4 gray24 g gray24 m gray24', & forrest@0: '",',/'" 77 s 0.77prob c gray23 g4 gray23 g gray23 m gray23', & forrest@0: '",',/'" 78 s 0.78prob c gray22 g4 gray22 g gray22 m gray22', & forrest@0: '",',/'" 79 s 0.79prob c gray21 g4 gray21 g gray21 m gray21', & forrest@0: '",',/'" 80 s 0.80prob c gray20 g4 gray20 g gray20 m gray20', & forrest@0: '",',/'" 81 s 0.81prob c gray19 g4 gray19 g gray19 m gray19', & forrest@0: '",',/'" 82 s 0.82prob c gray18 g4 gray18 g gray18 m gray18', & forrest@0: '",',/'" 83 s 0.83prob c gray17 g4 gray17 g gray17 m gray17', & forrest@0: '",',/'" 84 s 0.84prob c gray16 g4 gray16 g gray16 m gray16', & forrest@0: '",',/'" 85 s 0.85prob c gray15 g4 gray15 g gray15 m gray15', & forrest@0: '",',/'" 86 s 0.86prob c gray14 g4 gray14 g gray14 m gray14', & forrest@0: '",',/'" 87 s 0.87prob c gray13 g4 gray13 g gray13 m gray13', & forrest@0: '",',/'" 88 s 0.88prob c gray12 g4 gray12 g gray12 m gray12', & forrest@0: '",',/'" 89 s 0.89prob c gray11 g4 gray11 g gray11 m gray11', & forrest@0: '",',/'" 90 s 0.90prob c gray10 g4 gray10 g gray10 m gray10', & forrest@0: '",',/'" 91 s 0.91prob c gray9 g4 gray9 g gray9 m gray9', & forrest@0: '",',/'" 92 s 0.92prob c gray8 g4 gray8 g gray8 m gray8', & forrest@0: '",',/'" 93 s 0.93prob c gray7 g4 gray7 g gray7 m gray7', & forrest@0: '",',/'" 94 s 0.94prob c gray6 g4 gray6 g gray6 m gray6', & forrest@0: '",',/'" 95 s 0.95prob c gray5 g4 gray5 g gray5 m gray5', & forrest@0: '",',/'" 96 s 0.96prob c gray4 g4 gray4 g gray4 m gray4', & forrest@0: '",',/'" 97 s 0.97prob c gray3 g4 gray3 g gray3 m gray3', & forrest@0: '",',/'" 98 s 0.98prob c gray2 g4 gray2 g gray2 m gray2', & forrest@0: '",',/'" 99 s 0.99prob c gray1 g4 gray1 g gray1 m gray1', & forrest@0: '",',/'"100 s 1.00prob c black g4 black g black m black', & forrest@0: '",',/'"101 c white m white s NotFlammable', & forrest@0: '",',/'"102 c GreenYellow m gray70 s LP0', & forrest@0: '",',/'"103 c LawnGreen m gray60 s LP1', & forrest@0: '",',/'"104 c LimeGreen m gray50 s LP2', & forrest@0: '",',/'"105 c ForestGreen m gray40 s LP3', & forrest@0: '",',/'"106 c yellow m gray80 s NF', & forrest@0: '",',/'"107 c blue m black s water ",', & forrest@0: /'/* pixels */') forrest@0: ! forrest@0: ! forrest@0: ! load fractal probability terrain into itemp array forrest@0: do m = 1, numcells ! over all cells forrest@0: itemp(P(m,ihabs+1),P(m,ihabs+2)) = P(m,iz) forrest@0: end do forrest@0: ! forrest@0: ! forrest@0: ! equations for slope and aspect are from: forrest@0: ! Horn, B.K.P., 1981. Hill Shading and the Reflectance Map. forrest@0: ! Proceedings of the I.E.E.E., 69(1):14-47. forrest@0: ! forrest@0: ! calculate and display aspect map forrest@0: ! set sun altitude above horizon in degrees (0 to 90) forrest@0: alt = 40. forrest@0: ! set sun azimuth in degrees west of north (-1 to 360) forrest@0: az = 105. forrest@0: k = 0 forrest@0: do j = 0,irowcols forrest@0: do i = 0,irowcols forrest@0: ! forrest@0: x = (itemp(i-1,j-1)+2.*itemp(i,j-1)+itemp(i+1,j-1) & forrest@0: -itemp(i-1,j+1)-2.*itemp(i,j+1)-itemp(i+1,j+1)) & forrest@0: /(8.0*30.0) forrest@0: ! forrest@0: y = (itemp(i-1,j-1)+2.*itemp(i-1,j)+itemp(i-1,j+1) & forrest@0: -itemp(i+1,j-1)-2.*itemp(i+1,j)-itemp(i+1,j+1)) & forrest@0: /(8.0*30.0) forrest@0: ! forrest@0: slope = 90. - atand(sqrt(x*x + y*y)) forrest@0: if (x .eq. 0. .and. y .eq. 0.) then forrest@0: aspect = 360. forrest@0: else forrest@0: aspect = atan2d(x,y) forrest@0: end if forrest@0: ! forrest@0: c = sind(alt)*sind(slope) + cosd(alt)*cosd(slope)*cosd(az-aspect) forrest@0: ! forrest@0: if (c .le. 0.) then forrest@0: idispvec(k)=0 forrest@0: else forrest@0: idispvec(k)=nint(c*100.) forrest@0: end if forrest@0: ! forrest@0: ! flood with water if not painted forrest@0: if (itemp(i,j) .le. iflood) idispvec(k) = 107 forrest@0: ! forrest@0: ! print *, i,j,k,idispvec(k) forrest@0: k = k + 1 forrest@0: end do forrest@0: write (fmt,'("(",I0,"(i3)$)")') irowcols+1 forrest@0: write (50,'(2A$)') ' "' forrest@0: !write (50,108) (idispvec(l), l=0, irowcols) forrest@0: write (50,fmt) (idispvec(l), l=0, irowcols) forrest@0: write (50,'(2A)') '",' forrest@0: !108 format (' "',(i3),'",') forrest@0: k = 0 forrest@0: end do forrest@0: ! forrest@0: write (50,112) forrest@0: 112 format ('};') forrest@0: close (50) forrest@0: ! forrest@0: return forrest@0: end subroutine XPMRelief forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: !********************************************************** forrest@0: ! The following were added by Forrest Hoffman forrest@0: ! Sat Apr 7 23:20:48 EDT 2007 forrest@0: !********************************************************** forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: real function cosd(angle) forrest@0: implicit none forrest@0: real, intent(in) :: angle forrest@0: forrest@0: cosd = cos(angle * rad_per_deg) forrest@0: forrest@0: return forrest@0: end function cosd forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: real function sind(angle) forrest@0: implicit none forrest@0: real, intent(in) :: angle forrest@0: forrest@0: sind = sin(angle * rad_per_deg) forrest@0: forrest@0: return forrest@0: end function sind forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: real function atand(x) forrest@0: implicit none forrest@0: real, intent(in) :: x forrest@0: forrest@0: atand = atan(x) * deg_per_rad forrest@0: forrest@0: return forrest@0: end function atand forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: real function atan2d(x, y) forrest@0: implicit none forrest@0: real, intent(in) :: x forrest@0: real, intent(in) :: y forrest@0: forrest@0: atan2d = atan2(x, y) * deg_per_rad forrest@0: forrest@0: return forrest@0: end function atan2d forrest@0: ! forrest@0: !------------------------------------------------------------------------------ forrest@0: ! forrest@0: forrest@0: end module realizer_mod