subroutine prin(prefix)
C
C Outputs ascii array if ndim = 1, else if ndim > 1 then
C write out SDS data file containing all variables (plus time).
C
      include 'global.h'
      include 'sweep.h'
      include 'zone.h'
      character*1 char, coord
      character*4 tmp
      character*5 prefix
      character*10 filename
      integer  dssdims, dssdist, dssdisc, dssdast
      integer dsadata, dspdata, isdn, rank, shape(3)
      character*8 datalabel, dataunit, datafmt, dimlabels(3)
      character*8 dimunits(3), dimfmts(3)
C
C------------------------------------------------------------------------------
C
C Create filename from integer nfile (in global.h) and prefix such that filename
C looks like prefx.1000 where 1000 is the value of nfile
C
      write(tmp,910) nfile
 910  format(i4)
      do 920 i = 1, 4
         if ((tmp(i:i)) .eq. ' ') then
            tmp(i:i) = '0'
         endif
 920  continue
      write(filename,930) prefix,tmp
 930  format(a5,'.',a4)
      nfile = nfile + 1
c
C coordinate geometry
        if(ngeomy.lt.3) then
           coord = 'c'
        else
           coord = 's'
        endif
C
c Size of dimensions
	shape(1) = imax
	shape(2) = jmax
	shape(3) = kmax
c
c Set parameters for dimension arrays
	dimlabels(1) = 'X'
	dimlabels(2) = 'Y'
	dimlabels(3) = 'Z'
	dimunits(1)  = 'cm'
	dimunits(2)  = 'cm'
	dimunits(3)  = 'cm'
	dimfmts(1)   = 'E12.5'
	dimfmts(2)   = 'E12.5'
	dimfmts(3)   = 'E12.5'
C
      if (ndim .eq. 1) then
C
C Keep 1D output simple, just write out in ascii...
C
         j = 1
         k = 1
         open(unit=3,file=filename,form='formatted')
         do 10 i = 1, imax
            write(3, 1003) zxa(i), zrho(i,j,k)
     &                 ,zp(i,j,k), zux (i,j,k)
10       continue
         close(3)
C
      else if (ndim .eq. 2) then
C
c Treat higher D cases by writing an SDS file 
C
c Dimension of data array
	rank = 2
        k = 1
C
c  Load density in first HDF data file, including X and Y scales
	datalabel = 'Density'
	dataunit  = 'g*cm^-3'
	datafmt   = '1pe11.3'
	isdn = dssdims(rank, shape)
	isdn = dssdast(datalabel, dataunit, datafmt, coord )
	isdn = dssdist(1, dimlabels(1), dimunits(1), dimfmts(1))
	isdn = dssdist(2, dimlabels(2), dimunits(2), dimfmts(2))
	isdn = dssdisc(1, shape(1), zxa)
	isdn = dssdisc(2, shape(2), zya)
	isdn = dspdata(filename, 2, shape, zrho)
c
c  Load pressure in second HDF data file
c	
	datalabel = 'Pressure'
	dataunit  = 'g/cm/s^2'
	isdn = dssdast(datalabel, dataunit, datafmt, coord )
	isdn = dsadata(filename, 2, shape, zp)
c
c  Load Y velocity in third HDF data file
c
	datalabel = 'Yvelocty'
	dataunit  = 'cm *s^-1'
	isdn = dssdast(datalabel, dataunit, datafmt, coord )
	isdn = dsadata(filename, 2, shape, zuy)
c
c  Load X velocity in fourth HDF data file
c
	datalabel = 'Xvelocty'
	dataunit  = 'cm *s^-1'
	isdn = dssdast(datalabel, dataunit, datafmt, coord )
	isdn = dsadata(filename, 2, shape, zux)
c
c  Load time slice as a fifth data file of dimension (1)
c
	datalabel = 'time slice'
	dataunit  = 'second'
	isdn = dssdims(1, 1)
	isdn = dssdast(datalabel, dataunit, datafmt, coord)
	isdn = dsadata(filename, 1, 1, time)
C
      else 
C
c Dimension of data array
	rank = 3
C
c  Load density in first HDF data file, including X and Y scales
	datalabel = 'Density'
	dataunit  = 'g*cm^-3'
	datafmt   = '1pe11.3'
	isdn = dssdims(rank, shape)
	isdn = dssdast(datalabel, dataunit, datafmt, coord )
	isdn = dssdist(1, dimlabels(1), dimunits(1), dimfmts(1))
	isdn = dssdist(2, dimlabels(2), dimunits(2), dimfmts(2))
	isdn = dssdist(3, dimlabels(3), dimunits(3), dimfmts(3))
	isdn = dssdisc(1, shape(1), zxa)
	isdn = dssdisc(2, shape(2), zya)
	isdn = dssdisc(3, shape(3), zza)
	isdn = dspdata(filename, 3, shape, zrho)
c
c  Load pressure in second HDF data file
c	
	datalabel = 'Pressure'
	dataunit  = 'g/cm/s^2'
	isdn = dssdast(datalabel, dataunit, datafmt, coord )
	isdn = dsadata(filename, 3, shape, zp)
c
c  Load X velocity in third HDF data file
c
	datalabel = 'Xvelocty'
	dataunit  = 'cm/s    '
	isdn = dssdast(datalabel, dataunit, datafmt, coord )
	isdn = dsadata(filename, 3, shape, zux)
c
c  Load Y velocity in fourth HDF data file
c
	datalabel = 'Yvelocty'
	isdn = dssdast(datalabel, dataunit, datafmt, coord )
	isdn = dsadata(filename, 3, shape, zuy)
C
c  Load Z velocity in fourth HDF data file
c
	datalabel = 'Zvelocty'
	isdn = dssdast(datalabel, dataunit, datafmt, coord )
	isdn = dsadata(filename, 3, shape, zuz)

c  Load time slice as a fifth data file of dimension (1)
c
	datalabel = 'time slice'
	dataunit  = 'second'
	isdn = dssdims(1,1)
	isdn = dssdast(datalabel, dataunit, datafmt, coord)
	isdn = dsadata(filename, 1, 1, time)
C
      endif
C
      write(8,6000) filename, time, ncycle
C
 6000 format('Wrote ',a10,' to disk at time =',1pe9.2,' (ncycle =',
     &        i5,')')
 1003 format(' ',4e13.5)
C
      return
      end