c These parameter statements can only be used for a 2-D, 2-eq problem,
c spatially discretized using quadrilateral bilinear elements
c
c nelx is number of element in latitudinal direction
c nely is number of element in radial direction
c    (storage space is smallest if nelx<=nely)
c ityp is element type (=23 for bilinear)
c   nd is number of node per element (=4 for bilinear)
c  ngp is number of gauss point per dimension, for FEM element integration
c ntmx is maximal number of time steps
c---------------------------------------------------------------------
      integer    nelx,nely,ngp,nd,ityp,ntmx
      parameter(nelx=64, nely=96, ngp=2, nd=4, ityp=23, ntmx=12000 )
c---------------------------------------------------------------------
c from these are calculated:
c
c nel = total number of elements
c nnd = total number of mesh nodes
c nm  = total number of FEM coefficients
c nbc = maximal possible number of Dirichlet boundary conditions
c ng  = total number of gauss points per element
c mhb = half bandwidth of system matrices
c mfb = full bandwidth of system matrices 
c---------------------------------------------------------------------
      integer    nel,nm,ng,nnd,nbc,mhb,mfb,nelxp1,nelyp1
      parameter( nelxp1=nelx+1, nelyp1=nely+1 )
      parameter( nel=nelx*nely, nnd=(nelx+1)*(nely+1), ng=ngp*ngp )
      parameter( nm=2*nnd, nbc=4*(nelx+nely)+nely*nelx/10 )
      parameter( mhb=2*nelx+7, mfb=2*mhb-1 )
c---------------------------------------------------------------------
c Other constants and parameters defined here:
c 
c 1. time stepping (See Burnett, Sec. X.Y)
c
c theta = time-stepping parameter (0.5<theta<=1.0)
c nsteps= number of sub-timesteps
c---------------------------------------------------------------------
      real       theta
      integer    nsteps
      parameter( nsteps=1000,theta=2./3. )