c
define(IQX,59)
define(IQY,61)
define(NBDY,9)
define(NVAR,5)
define(MAXCOMMENTS,6)

      program test2d


      ifelse(ifdouble,1,implicit real*8 (a-h,o-z), )
c     ifelse(ifdouble,1,8,4)
      parameter(nxzones=IQX)
      parameter(nyzones=IQY)
      parameter(nbdy=NBDY)
      parameter(nvar=NVAR)
      parameter(maxcomments=MAXCOMMENTS)
c
c   These are 1d arrays which contain information defining the grid,
c      that is, the locations of the left hand zone interfaces.  One 
c      extra zone is requires on the left to fully specify the grid.

      dimension xintrl(1-nbdy:nxzones+nbdy+1)
      dimension yintrl(1-nbdy:nyzones+nbdy+1)
     
c
      character*20 filename
c
c   This holdes some commentary for inclusion in the compressed
c      dump file.
      character*80 comments(maxcomments)
c
c   This array contains zone averaged information of the state variables
c     in interlieaved order, making it a 3d array for two spatial dimensions
c     They are:
c        ddd(1,*,*) = density
c        ddd(2,*,*) = pressure
c        ddd(3,*,*) = x-velocity
c        ddd(4,*,*) = y-velocity
c        ddd(5,*,*) = diagnostic (an arbitrary quantity)
      dimension ddd(nvar,1-nbdy:nxzones+nbdy,
     ^                   1-nbdy:nyzones+nbdy)

c
c   Initialize the disk file name
      filename = 'q42a0000'
c
c   Initialize the commentary
      comments(1)='  This is the result from the program cdumps2dexample'
      comments(2)='  Its purpose is to demonstrate the PPMLIB compressed'
      comments(3)='  dumps routines. The arrays are REAL*ifelse(ifdouble,1,8,4)'
      ncomments = 3

c
c   Initialize the grid
      call init_grid2d(xintrl, yintrl, nxzones, nyzones, nbdy)
c
c   Initialize the zone average grid
      call init_ddd(ddd,nvar,nxzones,nyzones,nbdy)

c
c   Open the file, write the comments including the date of file
c     creation.
      iunit = 66
      call ppm98_binopen(filename,iunit)
      do i=1,ncomments
         call ppm98_write_comment(comments(i),iunit)
      enddo
      call ppm98_timestamp(iunit)
c
c   Now record the grid information in the compressed dump routine.  The
c     three axies will be assigned the names 'X', 'Y', and 'Z'.  These
c     should be capitol letters, as their derivatives will be references
c     with the equivalent lower case letters in a3d.
      call ifelse(ifdouble,1,d_,)ppm98_mesho2
     ^                           ('X', 'Y', nbdy, nxzones, nbdy,
     ^                            nbdy, nyzones, nbdy,
     ^                            xintrl, yintrl, iunit )

c
c   This illustrates adding a constant value to the compressed dumps file
      pi =  atan2(0.,-1.)
      call ppm98_var0d('ratio of circumference of circle to diameter', 
     ^               'PI',pi,iunit)

c
c   Now add the individual zone average quantities to the compressed dump
      call ifelse(ifdouble,1,d_,) ppm98_var2d
     ^                      ('Density', 'RHO', nbdy, nxzones, nbdy, 
     ^                                         nbdy, nyzones, nbdy,
     ^                                    nvar, 1, ddd, 0, iunit)
      call ifelse(ifdouble,1,d_,) ppm98_var2d
     ^                          ('Pressure', 'P', nbdy, nxzones, nbdy, 
     ^                                        nbdy, nyzones, nbdy,
     ^                                    nvar, 2, ddd, 0, iunit)
      call ifelse(ifdouble,1,d_,) ppm98_var2d
     ^                          ('X velocity', 'VX', nbdy, nxzones, nbdy, 
     ^                                           nbdy, nyzones, nbdy,
     ^                                    nvar, 3, ddd, 0, iunit)
      call ifelse(ifdouble,1,d_,) ppm98_var2d
     ^                          ('Y velocity', 'VY', nbdy, nxzones, nbdy, 
     ^                                           nbdy, nyzones, nbdy,
     ^                                    nvar, 4, ddd, 0, iunit)
      call ifelse(ifdouble,1,d_,) ppm98_var2d
     ^                          ('2D chess board diagnostic', 'DIAG',
     ^                                           nbdy, nxzones, nbdy,
     ^                                           nbdy, nyzones, nbdy,
     ^                                    nvar, 5, ddd, 0, iunit)

c   Close the file and exit
      call ppm98_binclose(iunit)


      return
      end

      subroutine init_ddd(ddd,nvar,nx,ny,nbdy)
      ifelse(ifdouble,1,implicit real*8 (a-h,o-z), )
      dimension ddd(nvar, 1-nbdy:nx+nbdy,
     ^                    1-nbdy:ny+nbdy)
c
c    ddd(1,*,*) = density
c    ddd(2,*,*) = pressure
c    ddd(3,*,*) = x-velocity
c    ddd(4,*,*) = y-velocity
c    ddd(5,*,*) = a diagnostic
c
c   The domain will be divided into quadrants for density and pressure
c
c   The velocities shall consist of a sine wave
      twopi =  2.0*atan2(0.,-1.)
      xfac = twopi/float(nx)
      yfac = twopi/float(ny)
c
c   The diagnostic shall consist of alternating cubes of 0.,1. in
c     a chess board pattern. 
      lcube = 16
c
      ycube = 0.0
      do iy=1-nbdy,ny+nbdy
	 if (iy.le.ny/2) then
	    y = 1.0
	 else
	    y = 2.0
	 endif
	 yvel = sin(float(iy)*yfac)
	 if (mod(iy,lcube) .eq. 0) ycube = 1.0 - ycube
	 xcube = ycube
         do ix=1-nbdy,nx+nbdy
	    if (ix.le.nx/2) then
	       x = 3.0
	    else
	       x = 4.0
	    endif
	    xvel = sin(float(ix)*xfac)
	    if (mod(ix,lcube) .eq. 0) xcube = 1.0 - xcube
	       ddd(1,ix,iy) = x*y
	       ddd(2,ix,iy) = x*y
	       ddd(3,ix,iy) = xvel
	       ddd(4,ix,iy) = yvel
	       ddd(5,ix,iy) = xcube
         enddo
      enddo
 
       return
      end
      subroutine init_grid2d(x1intrl,x2intrl,n1zones,n2zones,nbdy)
      ifelse(ifdouble,1,implicit real*8 (a-h,o-z), )
      dimension x1intrl(1-nbdy:n1zones+nbdy+1)
      dimension x2intrl(1-nbdy:n2zones+nbdy+1)

c   Set up a uniform 2d grid where the non boundary portion 
c      ranges from 0 to 1.

      dx = 1.0 / float(n1zones)
      do i=1-nbdy,n1zones+nbdy+1
        x1intrl(i) = float(i-1)*dx
      enddo

      dx = 1.0 / float(n2zones)
      do i=1-nbdy,n2zones+nbdy+1
        x2intrl(i) = float(i-1)*dx
      enddo

      return
      end