Dumping rasters

The following very simply modified version of main.f will accomplish the task:

PROGRAM main

  USE parameters; USE diffusion; USE fourier; USE pessl_hpf
  IMPLICIT NONE
  INTEGER                   :: n, ix, jx, iy, jy, k, i, j, stat, it, ij
  INTEGER                   :: num_errors
  REAL(long)                :: t1, t2, delx, dely
  INTEGER                   :: j1, j2
  REAL(long), ALLOCATABLE   :: dif_x(:), dif_y(:)
  REAL(long), ALLOCATABLE   :: lambda(:)
  REAL(long), ALLOCATABLE   :: gap(:)
  REAL(long), ALLOCATABLE   :: f(:,:), b(:,:), a(:), a_rep(:)
  REAL(long), ALLOCATABLE   :: x(:,:), ab(:), abt(:,:)
  REAL(long), ALLOCATABLE   :: xsines(:,:), ysines(:,:), temp(:,:)
  INTEGER                   :: biga_index, num_eigvalues, num_vectors
!HPF$  DISTRIBUTE (*,BLOCK) :: f, b, x, abt, temp
!HPF$  DISTRIBUTE (BLOCK)   :: a, ab

! Raster related variables and constants

  INTERFACE
     EXTRINSIC(hpf_serial) FUNCTION d8aimg(file_name, raster, nx, ny, mode)
       INTEGER :: d8aimg
       INTEGER :: nx, ny, mode
       INTEGER(kind=1) :: raster(nx * ny)
       CHARACTER(len=*) :: file_name
     END FUNCTION d8aimg
     EXTRINSIC(hpf_serial) FUNCTION d8pimg(file_name, raster, nx, ny, mode)
       INTEGER :: d8pimg
       INTEGER :: nx, ny, mode
       INTEGER(kind=1) :: raster(nx * ny)
       CHARACTER(len=*) :: file_name       
     END FUNCTION d8pimg
  END INTERFACE

  INTEGER(kind=1), ALLOCATABLE :: raster(:)
  INTEGER, PARAMETER :: COMP_RLE = 11
  REAL(long) :: temp_min, temp_max

  CALL init_diffusion

  ALLOCATE( dif_x(dif_numx*dif_numy) ); ALLOCATE( dif_y(dif_numy*dif_numx) )

  delx = PI/ ( dif_numx + 1.d0); dely = PI/ ( dif_numy + 1.d0)
  DO i = 1, dif_numx
     DO j = 1, dif_numy
        ij = dif_numx*(j-1)  + i
        dif_x(ij) = delx* i; dif_y(ij) = dely * j
     ENDDO
  ENDDO

  num_errors = 0; n = dif_nx * dif_ny

  ALLOCATE( lambda(n), stat=stat); IF( stat .NE. 0) num_errors = num_errors + 1
  ALLOCATE( a(n), stat=stat); IF( stat .NE. 0) num_errors = num_errors + 1
  ALLOCATE( ab(n), stat=stat); IF( stat .NE. 0) num_errors = num_errors + 1
  ALLOCATE( a_rep(n), stat=stat); IF( stat .NE. 0) num_errors = num_errors + 1
  ALLOCATE( xsines(dif_npts, dif_nx) , stat=stat ); IF( stat .NE. 0) &
       num_errors = num_errors + 1
  ALLOCATE( ysines(dif_npts, dif_ny) , stat=stat ); IF( stat .NE. 0) &
       num_errors = num_errors + 1
  ALLOCATE( temp(dif_npts, dif_ntemps) , stat=stat ); IF( stat .NE. 0) &
       num_errors = num_errors + 1
  ALLOCATE( raster(dif_npts), stat=stat ); IF( stat .NE. 0) &
       num_errors = num_errors + 1
  ALLOCATE( f(n,n) , stat=stat ); IF( stat .NE. 0) num_errors = num_errors + 1
  ALLOCATE( b(n,n) , stat=stat ); IF( stat .NE. 0) num_errors = num_errors + 1
  ALLOCATE( abt(n,dif_ntemps) , stat=stat ); IF( stat .NE. 0) &
       num_errors = num_errors + 1
  ALLOCATE( x(n,dif_ntemps) , stat=stat ); IF( stat .NE. 0) &
       num_errors = num_errors + 1

  IF( num_errors .GT. 0 ) THEN
     WRITE(*,*) 'Error in allocating arrays in main'; STOP
  ENDIF

  CALL expand_temp_profile(a)
  a_rep = a

  CALL get_diffusion_matrix(f)

  DO i = 1, dif_nx
     DO j = 1, dif_npts
        t1 = i*dif_x(j); xsines(j,i) = SQRT(2.d0/pi) * dsin(t1)
     ENDDO
  ENDDO

  DO i = 1, dif_ny
     DO j = 1, dif_npts
        t2 = i*dif_y(j); ysines(j,i) = SQRT(2.d0/pi) * dsin( t2)
     ENDDO
  ENDDO

  IF( num_errors .GT. 0 ) THEN
     WRITE(*,*) 'Error in allocating arrays for syevx in main'; STOP
  ENDIF

  CALL syevx(f,lambda,'u', z=b)

  ab(i) = 0.d0
  DO i = 1, n
     DO j = 1, n
        ab(i) = b(j,i)*a_rep(j) + ab(i)
     ENDDO
  ENDDO

  DO it = 1, dif_ntemps     
     i = 0
     DO i = 1, n
        abt(i,it) = EXP( -lambda(i) *(it-1)*dif_delta_t) * ab(i)
     ENDDO
  ENDDO

  CALL gemm(1.d0,b,abt,0.d0,x,'N','N')

  DO it = 1, dif_ntemps
     DO k = 1, dif_npts
        temp(k, it) = 0.d0
     ENDDO
  ENDDO

  DO it = 1, dif_ntemps
     DO i = 1, n  ! NOTE: there is a bug here in /usr/lpp/pessl/example/hpf
        iy = (i-1)/dif_nx+1; ix = MOD(i-1,dif_nx)+1
        DO k = 1, dif_npts
           temp(k,it) = temp(k,it) + xsines(k,ix) * ysines(k,iy) * x(i,it)
        ENDDO
     ENDDO
  ENDDO

  temp_min = MINVAL(temp)
  temp_max = MAXVAL(temp)

  DO it = 1, dif_ntemps
     raster = INT((temp(:,it) - temp_min)/(temp_max - temp_min) * 230 + 15, 1)
     IF (it .EQ. 1) THEN
        stat=d8pimg('movie.hdf', raster, dif_numx, dif_numy, COMP_RLE)
     ELSE
        stat=d8aimg('movie.hdf', raster, dif_numx, dif_numy, COMP_RLE)
     END IF
  END DO

  DEALLOCATE(xsines); DEALLOCATE(ysines); DEALLOCATE(lambda)
  DEALLOCATE(temp); DEALLOCATE(a); DEALLOCATE(abt); DEALLOCATE(a)
  DEALLOCATE(f); DEALLOCATE(dif_x); DEALLOCATE(dif_y); DEALLOCATE(raster)

  STOP

END PROGRAM main

The program is easiest to understand by looking at the differences between the original version of main.f and the new version with   diff -c.

Here is the first difference:

*** main.f	Mon Jan  4 17:43:59 1999
--- main.f.ORIG	Mon Jan  4 17:30:15 1999
***************
*** 16,42 ****
  !HPF$  DISTRIBUTE (*,BLOCK) :: f, b, x, abt, temp
  !HPF$  DISTRIBUTE (BLOCK)   :: a, ab
  
- ! Raster related variables and constants
- 
-   INTERFACE
-      EXTRINSIC(hpf_serial) FUNCTION d8aimg(file_name, raster, nx, ny, mode)
-        INTEGER :: d8aimg
-        INTEGER :: nx, ny, mode
-        INTEGER(kind=1) :: raster(nx * ny)
-        CHARACTER(len=*) :: file_name
-      END FUNCTION d8aimg
-      EXTRINSIC(hpf_serial) FUNCTION d8pimg(file_name, raster, nx, ny, mode)
-        INTEGER :: d8pimg
-        INTEGER :: nx, ny, mode
-        INTEGER(kind=1) :: raster(nx * ny)
-        CHARACTER(len=*) :: file_name       
-      END FUNCTION d8pimg
-   END INTERFACE
- 
-   INTEGER(kind=1), ALLOCATABLE :: raster(:)
-   INTEGER, PARAMETER :: comp_rle = 11
-   REAL(long) :: temp_min, temp_max
- 
    CALL init_diffusion
  
    ALLOCATE( dif_x(dif_numx*dif_numy) ); ALLOCATE( dif_y(dif_numy*dif_numx) )

The hyphens  at the beginning of the lines indicate text insertion.

What we have inserted is the definition of the interface, which you have already seen in one of the hints. We define a new alocatable array called raster too. COMP_RLE is a constant defined in one of the HDF header files, which specifies that RLE compression should be used when writing the images. Finally we have variables for maximum and minimum values of temperature.

Observe that raster is defined both in the interface and in the program as a one-dimensional variable. I have done it like that, because temp is a one dimensional variable in (x,y) space too. That way mapping between the two is much easier.

The second difference is as follows:

--- 16,21 ----
***************
*** 61,68 ****
         num_errors = num_errors + 1
    ALLOCATE( temp(dif_npts, dif_ntemps) , stat=stat ); IF( stat .NE. 0) &
         num_errors = num_errors + 1
-   ALLOCATE( raster(dif_npts), stat=stat ); IF( stat .NE. 0) &
-        num_errors = num_errors + 1
    ALLOCATE( f(n,n) , stat=stat ); IF( stat .NE. 0) num_errors = num_errors + 1
    ALLOCATE( b(n,n) , stat=stat ); IF( stat .NE. 0) num_errors = num_errors + 1
    ALLOCATE( abt(n,dif_ntemps) , stat=stat ); IF( stat .NE. 0) &

Here, again, we have just one insertion: we ALLOCATE space for raster.

The third and last difference is a little more elaborate:

--- 40,45 ----
***************
*** 128,148 ****
       ENDDO
    ENDDO
  
!   temp_min = MINVAL(temp)
!   temp_max = MAXVAL(temp)
! 
!   DO it = 1, dif_ntemps
!      raster = INT((temp(:,it) - temp_min)/(temp_max - temp_min) * 230 + 15, 1)
!      IF (it .EQ. 1) THEN
!         stat=d8pimg('movie.hdf', raster, dif_numx, dif_numy, COMP_RLE)
!      ELSE
!         stat=d8aimg('movie.hdf', raster, dif_numx, dif_numy, COMP_RLE)
!      END IF
!   END DO
  
    DEALLOCATE(xsines); DEALLOCATE(ysines); DEALLOCATE(lambda)
    DEALLOCATE(temp); DEALLOCATE(a); DEALLOCATE(abt); DEALLOCATE(a)
!   DEALLOCATE(f); DEALLOCATE(dif_x); DEALLOCATE(dif_y); DEALLOCATE(raster)
  
    STOP
  
--- 105,126 ----
       ENDDO
    ENDDO
  
!   WRITE(*,*) '   point #      X         Y'
!   DO i = 1, dif_npts
!      WRITE(*,'(2x, i6, 2x, 2f11.4)') i, dif_x(i), dif_y(i)
!   ENDDO
!   WRITE(*,*)
!   DO k = 1, dif_npts, 6
!      WRITE(*,*); WRITE(*,'(30X,''Points'')')
!      WRITE(*,'(''   TIME   '',6(5x,''#'', i4))') (i, i=k, k+5)
!      DO i = 1, dif_ntemps
!         WRITE(*,'(7f10.5)') i*dif_delta_t, (temp(j,i),j=k,MIN(k+5,dif_npts))
!      ENDDO
!   ENDDO
  
    DEALLOCATE(xsines); DEALLOCATE(ysines); DEALLOCATE(lambda)
    DEALLOCATE(temp); DEALLOCATE(a); DEALLOCATE(abt); DEALLOCATE(a)
!   DEALLOCATE(f); DEALLOCATE(dif_x); DEALLOCATE(dif_y)
  
    STOP

Exclamation marks   in this listing indicate a replacement. So program diff noticed that we have replaced the two loops that did writing on standard output with a new loop that writes the images.

Before entering that loop we find the minimum and maximum values of temperature reached throughout the whole evolution of the system. We need those in order to scale all images uniformly, so that they would fit into the colour map and that the map itself would be the same for all images. Within the loop itself we cast scaled temperature slices onto raster and then write the rasters on a file movie.hdf. Observe that even though raster is a one-dimensional variable, we give its x and y dimensions here.

In Fortran-77 all arrays are one-dimensional anyway, and they can be indexed either with a collapsed index or with multidimensional indexes. You cannot make the same assumption in HPF though, especially with distributed arrays, which can be laid out as is required for the most efficient transfer of data between nodes. HPF vendors are no longer obliged to store distributed Fortran arrays in a column-major mode. But raster is a replicated array, and d8pimg and d8aimg are F77 subroutines, so this should still work.

The last difference is that we DEALLOCATE raster when it's no longer needed.

You can see (e.g., with NCSA Collage) rasters generated by this program on the following three files in

/afs/ovpit.indiana.edu/common/www/htdocs/gustav

1.

movie-1-1.hdf which corresponds to the following namelist

&input ly_ratio =1.d0, delta_t =0.1, numx =100, numy =200,
nx =7, ny =7, numt = 20, init_f = 1, diff_f = 1, /

2.

movie-4-3.hdf which corresponds to the following namelist

&input ly_ratio =1.d0, delta_t =0.1, numx =100, numy =200,
nx =7, ny =7, numt = 20, init_f = 4, diff_f = 3, /

3.

movie-0.25-4-3.hdf which corresponds to the following namelist:

&input ly_ratio =0.25d0, delta_t =0.05, numx =100, numy =200,
nx =7, ny =7, numt = 40, init_f = 4, diff_f = 3, /

This is a more interesting example, because here we have changed the geometry of the beam, so that l_r = l_x/l_y = 1/2. Recall that ly_ratio corresponds to l_r². Also, observe that the generated image has the same geometry: numx/numy = 1/2.