c  **********************************************************************
c  *									*
c  *  File: mesh.for							*
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  *  20.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 angle_step
  
      real*8 anglez(dims),r2(dims), z2(dims), z(dims)
      real*8 anglexy(nmax)
      real*8 r(dims,nmax),angle(dims,nmax) 
      real*8 xx(dims,nmax),yy(dims,nmax),zz(dims,nmax)
      integer*4 mpTimerElapsed

c compile directives for CM5
cmf$  layout xx(:news,:news), yy(:news,:news),zz(:news,:news)  
cmf$  layout r(:news,:news),angle(:news,:news)
cmf$  layout anglez(:news),r2(:news), z2(:news), z(:news),anglexy(:news)

c compile directives for DECmpp-12000 machine 
CMPF  ONDPU xx, yy, zz, r,angle,anglez,r2, z2, z,anglexy

      character*20 fname
      integer i, j

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

c     --- CM Timer
c      call CM_timer_clear(1)
c      call CM_timer_start(1)

c     --- Maspar Timer
      call mpTimerStart()

      r2 = points*points
      angle_step = pi/float(dims-1)

      forall (i=1:dims) anglez(i) = float(i-1)*angle_step
      forall (i=1:nmax) anglexy(i) = float(i-1)*angle_step
      angle = spread(anglexy, dim=1, ncopies=dims)

      z = points*cos(anglez)
      z2 = z*z
      r = spread(sqrt(abs(r2 - z2)), dim=2, ncopies=nmax)

      xx = r*cos(angle)
      yy = r*sin(angle)
      zz = spread(z, dim=2, ncopies=nmax)

c     --- CM Timer
c      call CM_timer_stop(1)
c      call CM_timer_print(1)

      print *, 'Time = ',mpTimerElapsed()/1000.0,' 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

c  100  format(<3*nmax>f15.6)  ! CMfortran/MPfortran doesn't allow this format
 100  format(10000(3f15.6))

      close(11)
 
      print *,'total points:',nmax*dims
      stop
      end
