c-----
c A simple FORTRAN code to demonstrate the creation of an ADPAC mesh
c for a flat plate
c
c Christopher J. Miller		
c 12 Oct 1999
c-----
	parameter (imx = 81)
	parameter (jmx = 81)
	parameter (kmx =  1)
	dimension x(imx,jmx,kmx)
	dimension y(imx,jmx,kmx)
	dimension z(imx,jmx,kmx)
	dimension pi(imx),pj(jmx),pk(kmx)
c-----
c The mesh block coordinates, i.e. x(i,j,k), are such that increasing..
c	x and i are downstream
c	y and j are radially outboard
c	z and k are orthogonal to x & y in a LEFT HANDED sense.
c
c NOTE: ADPAC uses a left handed grid.  If you don't make it left handed
c ADPAC will find the cell volumes are negative and complain.
c-----
c Generate an algebraic grid: uniform in x and z, packed towards the
c inner and outer walls of the duct.  Adpac likes stretching ratios
c that are 1/1.2 <= ratio <= 1.2.  Outside that is ok but...
c-----
	call pack(pi,imx,1.1,1/1.1)
	call pack(pj,jmx,1.1,1/1.1)
	pk(1) = 0.0
	do k = 1,kmx
	    do j = 1,jmx
		do i = 1,imx
		    x(i,j,k) = pi(i)
		    y(i,j,k) = 0.05*pj(j)
		    z(i,j,k) = pk(k)
		enddo
	    enddo
	enddo
c-----
c Output the grid block:
c-----
	open (10,file='block1.unf',form='unformatted')
	write(10) imx,jmx,kmx
	print *,'block1.unf: ',imx,jmx,kmx
	write(10)
     >		(((x(i,j,k),i=1,imx),j=1,jmx),k=1,kmx),
     >		(((y(i,j,k),i=1,imx),j=1,jmx),k=1,kmx),
     >		(((z(i,j,k),i=1,imx),j=1,jmx),k=1,kmx)
	close(10)
c-----
	stop
	end
c-----------------------------------------------------------------------
        subroutine pack(pz,kmx,r1,r2)
c-----
c lay out a geometric packing
c cjm 960214
c-----
        dimension pz(kmx)
c
	pz(1) = 0.0
        pz(2) = 1.0
	if(kmx.lt.3) return
c-----
c lay out the ratios
c-----
	if(kmx.eq.3) then
	    ratio = sqrt(r1 * 1/r2)
	    pz(2) = ratio
	else
            do k = 2,kmx-1
		fz = float(k-2)/float(kmx-3)
		ratio = r1**(1-fz)/r2**fz
		pz(k) = ratio
	    enddo
	endif
c-----
c lay out the spacings
c-----
        savrat = pz(2)
	pz(1)  = 0.0
        pz(2)  = 1.0
	delta  = pz(2)-pz(1)
        do k = 3,kmx
	    ratio = savrat
	    delta = delta*ratio
	    savrat = pz(k)
	    pz(k)  = pz(k-1) + delta
        enddo
        do k = 2,kmx
	    pz(k) = pz(k)/pz(kmx)
        enddo
        return
        end