c  **********************************************************************
c  *									*
c  *  File: mesh.f     							*
c  *									*
c  *  Description:							*
c  *       example program  to generate a mesh geometry data in AVS for *
c  *		a sphere                                                *
c  *									*
c  *       Usage:							*
c  *       mesh								*
c  *	   (type in output file name <fname.mesh>)                      *
c  *									*
c  *  History:								*
c  *									*
c  *  17.06.93      Gang Cheng      Created 				*
c  *									*
c  *									*
c  *  Copyright 1993, Northeast Parallel Architectures Center		*
c  *									*
c  **********************************************************************

      program mesh_to_sphere

      implicit none

      integer dims,nmax
      real*8 points,step,pi
      parameter (step=0.025)
      parameter (points = 1.0, pi=3.141592)
      parameter (dims = INT(2*points/step+1))
      parameter (nmax=2*(dims-1)+1)

      real*8 x,y,z,z2,r2,r,angle_step,angle,anglez
      real*8 xx(dims,nmax),yy(dims,nmax),zz(dims,nmax)
      real timediff,timearray(2),dtime

      character*20 fname
      integer i, j

      print *, 'Output data file name(.mesh):' 
      read *,fname

c     --- Timer
      timediff = dtime(timearray)

      r2 = points*points
      angle_step = pi/float(dims-1.0)
      anglez = 0.0
      do i = 1,dims
         z = points*cos(anglez)
         z2 = z*z
         r = sqrt(abs(r2 - z2))
         angle = 0.0
         do j = 1, nmax
            x = r*cos(angle)
            y = r*sin(angle)
            xx(i,j) = x
            yy(i,j) = y
            zz(i,j) = z
            angle = angle + angle_step
         enddo
         anglez = anglez + angle_step 
      enddo
         
c     --- Timer
      timediff = dtime(timearray)
      print *, 'Calculation time = ',timediff,' seconds.'

      open(unit=11,file=fname,status='UNKNOWN')
      write(11,*) 'vertex ',dims,nmax
      
      do i=1, dims
            write(11,100)(xx(i,j),yy(i,j),zz(i,j), j=1,nmax)
      enddo

 100  format(<3*nmax>f15.6)
      close(11)
 
      print *,'total points:',nmax*dims
      stop
      end
