       Program makeq
C***********************************************************************
C
C    This program will generate an inviscid vortex in a cartesian grid.
C
c    Dimension parameter : maxpt >= max (number of points in input grid,
c    -------------------                 number of points in output grid)
c
c    Internal memory required = 6*maxpt
c    ------------------------
c
C***********************************************************************
      parameter (maxpt=300001)
      PARAMETER ( NQ=6 )
      dimension x(maxpt),y(maxpt),z(maxpt)
      dimension q(maxpt,nq)
c
c    Read input grid
      call initio(maxpt,jmax,kmax,lmax)
      call readgr(jmax,kmax,lmax,x,y,z)
c
c    Menu for generating vortex file
      call dynm(jmax,kmax,lmax,x, z, y, q,qn,nq)
c
c    Write out final final q-file
      call writq(jmax,kmax,lmax,q,qn,nq)
c
      stop
      end
c***********************************************************************
      subroutine initio(maxpt,jmax,kmax,lmax)
 5    continue
c
c    Open files and read dimensions
       open(unit=2,file='mkgrid.fmt.x',status='old',form='formatted')
       read(2,*) jmax,kmax,lmax
c
       write(6,*)'Input Grid Dimensions:  ', jmax,kmax,lmax
       write(*,*)

c    Check dimensions
      npt = jmax*kmax*lmax
      if (npt.gt.maxpt) then
       write(*,*)'Total number of points in input grid =',npt
       write(*,*)'Maximum number in dimension statement =',maxpt
       write(*,*)'Max dimension exceeded. Please recompile code with ',
     >           'larger maxpt'
       stop
      endif
c
 100  continue
      return
      end
C
C***********************************************************************
C
      subroutine dynm(jmax,kmax,lmax,x,y,z,q,qn,nq)
      PARAMETER ( PI=3.141592653589793D0 )
c
      dimension x(jmax,kmax,lmax),y(jmax,kmax,lmax),z(jmax,kmax,lmax)
      dimension q(jmax,kmax,lmax,nq)
c
      character*80 file
C
C      Default values
C
      DT     = 0.1
      TIME   = 0.
      FSMACH = 0.5

c     initial x location (x00) and last x location (xm1)

      X00    = 5.
      XM1    = 5.-DT*FSMACH
      Z00    = 5.
      Y00    = 0.
      STRNTH = 5.
      GAMINF = 1.4
      GM1    = GAMINF - 1.
C
      RINF   = 1.
      CINF   = 1.
      PINF   = RINF*CINF**2/GAMINF
      VMINF  = FSMACH*CINF
      UINF   = VMINF
      VINF   = 0.
      WINF   = 0.
      AFAC   = STRNTH/(2.*PI)
      BFAC   = -0.5*GM1/GAMINF*AFAC**2
      TINF   = PINF/RINF
      SINF   = PINF/RINF**GAMINF
c
c--- Perturbation due to vortex ---
C    present time
c
      DO 11 L = 1,LMAX
      DO 11 K = 1,KMAX
      DO 11 J = 1,JMAX

       XBAR      = X(J,K,L) - X00
C      YBAR      = Y(J,K,L) - Y00
       ZBAR      = Z(J,K,L) - Z00
       RSQ       = XBAR**2 + ZBAR**2

       EE     = EXP(0.5*(1.-RSQ))
       U      = UINF - (AFAC*EE)*ZBAR
       V      = VINF
       W      = WINF + (AFAC*EE)*XBAR
       V2     = 0.5*(U**2 + V**2 + W**2)
C
       TPR    = BFAC*EE**2
       T      = TINF + TPR
       P      = (T**GAMINF/SINF)**(1./GM1)
       RHO    = P/T

       q(j,k,l,1) = RHO
       q(j,k,l,2) = RHO*U
       q(j,k,l,3) = RHO*W
       q(j,k,l,4) = RHO*V
       q(j,k,l,5) = P/GM1 + RHO*V2
       q(j,k,l,6) = GAMINF

11    CONTINUE
C
      return
      end
c***********************************************************************
      subroutine readgr(jmax,kmax,lmax,x,y,z)
c
      dimension x(jmax,kmax,lmax),y(jmax,kmax,lmax),z(jmax,kmax,lmax)
c
c    Read input grid
        read(2,*) (((x(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax),
     >          (((y(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax),
     >          (((z(j,k,l),j=1,jmax),k=1,kmax),l=1,lmax)
c
      return
      end
c***********************************************************************
      subroutine writq(jmax,kmax,lmax,q,qn,nq)
c
      dimension q(jmax,kmax,lmax,nq)
      character*40 name
      parameter ( zero=0. )
c
      fsmach  = 0.5
      name    = 'vortex'

      open(unit=4,file='voric.fmt.q',status='unknown',form='formatted')
      write(4,*) jmax,kmax,lmax
      WRITE(4,*) fsmach,zero,zero,zero,zero
      write(4,*)((((q(j,k,l,n),j=1,jmax),k=1,kmax),l=1,lmax),n=1,5)


      close( 4 )
c
 100  continue
c
      return
      end
