diff --git a/src/boundaries.F90 b/src/boundaries.F90
index c8de844..ee6aee4 100644
--- a/src/boundaries.F90
+++ b/src/boundaries.F90
@@ -247,7 +247,7 @@ module boundaries
     use coordinates   , only : ib, jb, kb, ie, je, ke
     use coordinates   , only : ibl, jbl, kbl, iel, jel, kel
     use coordinates   , only : ibu, jbu, kbu, ieu, jeu, keu
-    use variables     , only : nt
+    use equations     , only : nv
 
 ! local variables are not implicit by default
 !
@@ -273,7 +273,7 @@ module boundaries
 
 ! iterate over all variables
 !
-      do n = 1, nt
+      do n = 1, nv
 
 ! edges
 !
@@ -450,7 +450,7 @@ module boundaries
     use mpitools      , only : periodic
 #ifdef MPI
     use mpitools      , only : nproc, nprocs, npmax
-    use variables     , only : nt
+    use equations     , only : nv
 #endif /* MPI */
 
 ! local variables are not implicit by default
@@ -677,11 +677,11 @@ module boundaries
 !
               select case(idir)
               case(1)
-                allocate(rbuf(nblocks,nt,nd,jm,km))
+                allocate(rbuf(nblocks,nv,nd,jm,km))
               case(2)
-                allocate(rbuf(nblocks,nt,im,nd,km))
+                allocate(rbuf(nblocks,nv,im,nd,km))
               case(3)
-                allocate(rbuf(nblocks,nt,im,jm,nd))
+                allocate(rbuf(nblocks,nv,im,jm,nd))
               end select
 
 ! if isend == nproc we are sending data
@@ -839,7 +839,7 @@ module boundaries
     use mpitools      , only : periodic
 #ifdef MPI
     use mpitools      , only : nproc, nprocs, npmax
-    use variables     , only : nt
+    use equations     , only : nv
 #endif /* MPI */
 
 ! local variables are not implicit by default
@@ -1068,11 +1068,11 @@ module boundaries
 !
           select case(idir)
           case(1)
-            allocate(rbuf(nblocks,nt,nh+2,jm,km))
+            allocate(rbuf(nblocks,nv,nh+2,jm,km))
           case(2)
-            allocate(rbuf(nblocks,nt,im,nh+2,km))
+            allocate(rbuf(nblocks,nv,im,nh+2,km))
           case(3)
-            allocate(rbuf(nblocks,nt,im,jm,nh+2))
+            allocate(rbuf(nblocks,nv,im,jm,nh+2))
           end select
 
 ! if isend == nproc we are sending data
@@ -1239,7 +1239,7 @@ module boundaries
     use mpitools      , only : periodic
 #ifdef MPI
     use mpitools      , only : nproc, nprocs, npmax
-    use variables     , only : nt
+    use equations     , only : nv
 #endif /* MPI */
 
 ! local variables are not implicit by default
@@ -1448,12 +1448,12 @@ module boundaries
 !
           select case(idir)
           case(1)
-            allocate(rbuf(nblocks,nt,ng,jm,km))
+            allocate(rbuf(nblocks,nv,ng,jm,km))
           case(2)
-            allocate(rbuf(nblocks,nt,im,ng,km))
+            allocate(rbuf(nblocks,nv,im,ng,km))
 #if NDIMS == 3
           case(3)
-            allocate(rbuf(nblocks,nt,im,jm,ng))
+            allocate(rbuf(nblocks,nv,im,jm,ng))
 #endif /* NDIMS == 3 */
           end select
 
@@ -1601,7 +1601,7 @@ module boundaries
     use coordinates, only : im, jm, km
     use mpitools , only : send_real_array, receive_real_array
     use mpitools , only : nprocs, nproc
-    use variables, only : nqt
+    use equations, only : nv
 #endif /* MPI */
 
     implicit none
@@ -1800,12 +1800,12 @@ module boundaries
 !
             select case(idir)
             case(1)
-              allocate(rbuf(nblocks,nqt,jm,km))
+              allocate(rbuf(nblocks,nv,jm,km))
             case(2)
-              allocate(rbuf(nblocks,nqt,im,km))
+              allocate(rbuf(nblocks,nv,im,km))
 #if NDIMS == 3
             case(3)
-              allocate(rbuf(nblocks,nqt,im,jm))
+              allocate(rbuf(nblocks,nv,im,jm))
 #endif /* NDIMS == 3 */
             end select
 
@@ -2194,7 +2194,7 @@ module boundaries
 
     use blocks   , only : block_data
     use coordinates, only : ng, im, jm, km, ibl, ieu, jbl, jeu, kbl, keu
-    use variables, only : nqt
+    use equations, only : nv
 
     implicit none
 
@@ -2211,26 +2211,26 @@ module boundaries
     case(1)
 
       if (iside .eq. 1) then
-        pdata%u(1:nqt,  1:ibl,1:jm,1:km) = u(1:nqt,1:ng,1:jm,1:km)
+        pdata%u(1:nv,  1:ibl,1:jm,1:km) = u(1:nv,1:ng,1:jm,1:km)
       else
-        pdata%u(1:nqt,ieu:im ,1:jm,1:km) = u(1:nqt,1:ng,1:jm,1:km)
+        pdata%u(1:nv,ieu:im ,1:jm,1:km) = u(1:nv,1:ng,1:jm,1:km)
       end if
 
     case(2)
 
       if (iside .eq. 1) then
-        pdata%u(1:nqt,1:im,  1:jbl,1:km) = u(1:nqt,1:im,1:ng,1:km)
+        pdata%u(1:nv,1:im,  1:jbl,1:km) = u(1:nv,1:im,1:ng,1:km)
       else
-        pdata%u(1:nqt,1:im,jeu:jm ,1:km) = u(1:nqt,1:im,1:ng,1:km)
+        pdata%u(1:nv,1:im,jeu:jm ,1:km) = u(1:nv,1:im,1:ng,1:km)
       end if
 
 #if NDIMS == 3
     case(3)
 
       if (iside .eq. 1) then
-        pdata%u(1:nqt,1:im,1:jm,  1:kbl) = u(1:nqt,1:im,1:jm,1:ng)
+        pdata%u(1:nv,1:im,1:jm,  1:kbl) = u(1:nv,1:im,1:jm,1:ng)
       else
-        pdata%u(1:nqt,1:im,1:jm,keu:km ) = u(1:nqt,1:im,1:jm,1:ng)
+        pdata%u(1:nv,1:im,1:jm,keu:km ) = u(1:nv,1:im,1:jm,1:ng)
       end if
 #endif /* NDIMS == 3 */
 
@@ -2253,7 +2253,7 @@ module boundaries
     use coordinates, only : ng, im, ih, ib, ie, ieu           &
                         , nd, jm, jh, jb, je, jeu           &
                         , nh, km, kh, kb, ke, keu
-    use variables, only : nqt
+    use equations, only : nv
 
     implicit none
 
@@ -2434,7 +2434,7 @@ module boundaries
                              , nd, jm, jh, jb, je, jeu                         &
                              , nh, km, kh, kb, ke, keu
     use interpolations, only : minmod3
-    use variables     , only : nqt
+    use equations     , only : nv
 
     implicit none
 
@@ -2707,7 +2707,7 @@ module boundaries
           it = 2 * (i - il) + is
           ip = it + 1
 
-          do q = 1, nqt
+          do q = 1, nv
 
             dul = u(q,i  ,j,k) - u(q,i-1,j,k)
             dur = u(q,i+1,j,k) - u(q,i  ,j,k)
@@ -2760,14 +2760,9 @@ module boundaries
     use blocks       , only : block_data
     use coordinates  , only : ng, im, jm, km, ib, ibl, ie, ieu, jb    &
                             , jbl, je, jeu, kb, kbl, ke, keu
+    use equations    , only : idn, imx, imy, imz, ibx, iby, ibz, ibp
     use error        , only : print_warning
-    use variables    , only : nvr, nfl, idn, imx, imy, imz
-#ifdef MHD
-    use variables    , only : ibx, iby, ibz
-#ifdef GLM
-    use variables    , only : iph
-#endif /* GLM */
-#endif /* MHD */
+    use equations    , only : nv
 
     implicit none
 
@@ -2824,7 +2819,7 @@ module boundaries
           pdata%u(  iby,i,:,:) = pdata%u(  iby,ib,:,:)
           pdata%u(  ibz,i,:,:) = pdata%u(  ibz,ib,:,:)
 #ifdef GLM
-          pdata%u(  iph,i,:,:) = pdata%u(  iph,ib,:,:)
+          pdata%u(  ibp,i,:,:) = pdata%u(  ibp,ib,:,:)
 #endif /* GLM */
         end do
 
@@ -2869,7 +2864,7 @@ module boundaries
         do i = 1, ng
           pdata%u(  :,i,:,:) = pdata%u(:,ib,:,:)
 #if defined MHD && defined GLM
-          pdata%u(iph,i,:,:) = 0.0d0
+          pdata%u(ibp,i,:,:) = 0.0d0
 #endif /* MHD & GLM */
         end do
 
@@ -2901,7 +2896,7 @@ module boundaries
           pdata%u(  iby,i,:,:) = pdata%u(  iby,ie,:,:)
           pdata%u(  ibz,i,:,:) = pdata%u(  ibz,ie,:,:)
 #ifdef GLM
-          pdata%u(  iph,i,:,:) = pdata%u(  iph,ie,:,:)
+          pdata%u(  ibp,i,:,:) = pdata%u(  ibp,ie,:,:)
 #endif /* GLM */
         end do
 
@@ -2946,7 +2941,7 @@ module boundaries
         do i = ieu, im
           pdata%u(  :,i,:,:) = pdata%u(:,ie,:,:)
 #if defined MHD && defined GLM
-          pdata%u(iph,i,:,:) = 0.0d0
+          pdata%u(ibp,i,:,:) = 0.0d0
 #endif /* MHD & GLM */
         end do
 
@@ -2978,7 +2973,7 @@ module boundaries
           pdata%u(  ibx,:,j,:) = pdata%u(  ibx,:,jb,:)
           pdata%u(  ibz,:,j,:) = pdata%u(  ibz,:,jb,:)
 #ifdef GLM
-          pdata%u(  iph,:,j,:) = pdata%u(  iph,:,jb,:)
+          pdata%u(  ibp,:,j,:) = pdata%u(  ibp,:,jb,:)
 #endif /* GLM */
         end do
 
@@ -3029,7 +3024,7 @@ module boundaries
           pdata%u(iby,:,j,:) =  0.0d0
           pdata%u(ibz,:,j,:) =  0.0d0
 #ifdef GLM
-          pdata%u(iph,:,j,:) =  0.0d0
+          pdata%u(ibp,:,j,:) =  0.0d0
 #endif /* GLM */
 #endif /* MHD */
         end do
@@ -3039,7 +3034,7 @@ module boundaries
         do j = 1, ng
           pdata%u(  :,:,j,:) = pdata%u(:,:,jb,:)
 #if defined MHD && defined GLM
-          pdata%u(iph,:,j,:) = 0.0d0
+          pdata%u(ibp,:,j,:) = 0.0d0
 #endif /* MHD & GLM */
         end do
 
@@ -3071,7 +3066,7 @@ module boundaries
           pdata%u(  ibx,:,j,:) = pdata%u(  ibx,:,je,:)
           pdata%u(  ibz,:,j,:) = pdata%u(  ibz,:,je,:)
 #ifdef GLM
-          pdata%u(  iph,:,j,:) = pdata%u(  iph,:,je,:)
+          pdata%u(  ibp,:,j,:) = pdata%u(  ibp,:,je,:)
 #endif /* GLM */
         end do
 
@@ -3122,7 +3117,7 @@ module boundaries
           pdata%u(iby,:,j,:) =  0.0d0
           pdata%u(ibz,:,j,:) =  0.0d0
 #ifdef GLM
-          pdata%u(iph,:,j,:) =  0.0d0
+          pdata%u(ibp,:,j,:) =  0.0d0
 #endif /* GLM */
 #endif /* MHD */
         end do
@@ -3132,7 +3127,7 @@ module boundaries
         do j = jeu, jm
           pdata%u(  :,:,j,:) = pdata%u(:,:,je,:)
 #if defined MHD && defined GLM
-          pdata%u(iph,:,j,:) = 0.0d0
+          pdata%u(ibp,:,j,:) = 0.0d0
 #endif /* MHD & GLM */
         end do
 
@@ -3165,7 +3160,7 @@ module boundaries
           pdata%u(  ibx,:,:,k) = pdata%u(  ibx,:,:,kb)
           pdata%u(  iby,:,:,k) = pdata%u(  iby,:,:,kb)
 #ifdef GLM
-          pdata%u(  iph,:,:,k) = pdata%u(  iph,:,:,kb)
+          pdata%u(  ibp,:,:,k) = pdata%u(  ibp,:,:,kb)
 #endif /* GLM */
         end do
 
@@ -3202,7 +3197,7 @@ module boundaries
         do k = 1, ng
           pdata%u(  :,:,:,k) = pdata%u(:,:,:,kb)
 #if defined MHD && defined GLM
-          pdata%u(iph,:,:,k) = 0.0d0
+          pdata%u(ibp,:,:,k) = 0.0d0
 #endif /* MHD & GLM */
         end do
 
@@ -3234,7 +3229,7 @@ module boundaries
           pdata%u(  ibx,:,:,k) = pdata%u(  ibx,:,:,ke)
           pdata%u(  iby,:,:,k) = pdata%u(  iby,:,:,ke)
 #ifdef GLM
-          pdata%u(  iph,:,:,k) = pdata%u(  iph,:,:,ke)
+          pdata%u(  ibp,:,:,k) = pdata%u(  ibp,:,:,ke)
 #endif /* GLM */
         end do
 
@@ -3271,7 +3266,7 @@ module boundaries
         do k = keu, km
           pdata%u(  :,:,:,k) = pdata%u(:,:,:,ke)
 #if defined MHD && defined GLM
-          pdata%u(iph,:,:,k) = 0.0d0
+          pdata%u(ibp,:,:,k) = 0.0d0
 #endif /* MHD & GLM */
         end do
 
diff --git a/src/equations.F90 b/src/equations.F90
index 5ac8d24..f81cc8d 100644
--- a/src/equations.F90
+++ b/src/equations.F90
@@ -115,6 +115,9 @@ module equations
   public :: gamma
   public :: csnd
   public :: cmax, cmax2
+  public :: nv
+  public :: idn, ivx, ivy, ivz, imx, imy, imz
+  public :: ibx, iby, ibz, ibp, ipr, ien
 
 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 !
diff --git a/src/evolution.F90 b/src/evolution.F90
index ecd816c..c239461 100644
--- a/src/evolution.F90
+++ b/src/evolution.F90
@@ -187,7 +187,7 @@ module evolution
     use blocks        , only : block_data, list_data
     use coordinates   , only : adx, ady, adz
     use coordinates   , only : im, jm, km
-    use variables     , only : nfl
+    use equations     , only : nv
 
 ! local variables are not implicit by default
 !
@@ -199,8 +199,8 @@ module evolution
 
 ! local arrays
 !
-    real, dimension(3)            :: dh
-    real, dimension(nfl,im,jm,km) :: du
+    real, dimension(3)           :: dh
+    real, dimension(nv,im,jm,km) :: du
 !
 !-------------------------------------------------------------------------------
 !
@@ -225,7 +225,7 @@ module evolution
 
 ! update the solution for the fluid variables
 !
-      pblock%u1(1:nfl,:,:,:) = pblock%u0(1:nfl,:,:,:) + du(1:nfl,:,:,:)
+      pblock%u1(1:nv,:,:,:) = pblock%u0(1:nv,:,:,:) + du(1:nv,:,:,:)
 
 ! update the conservative variable pointer
 !
@@ -266,8 +266,8 @@ module evolution
 
 ! update the solution for the fluid variables
 !
-      pblock%u0(1:nfl,:,:,:) = 0.5d0 * (pblock%u0(1:nfl,:,:,:)                 &
-                                   + pblock%u1(1:nfl,:,:,:) + du(1:nfl,:,:,:))
+      pblock%u0(1:nv,:,:,:) = 0.5d0 * (pblock%u0(1:nv,:,:,:)                 &
+                                   + pblock%u1(1:nv,:,:,:) + du(1:nv,:,:,:))
 
 ! update the conservative variable pointer
 !
diff --git a/src/forcing.F90 b/src/forcing.F90
index 57e8291..421ff05 100644
--- a/src/forcing.F90
+++ b/src/forcing.F90
@@ -510,7 +510,7 @@ module forcing
     use constants, only : dpi
     use coordinates, only : im, jm, km, ib, ie, jb, je, kb, ke
     use coordinates, only : ax, ay, az, advol
-    use variables, only : idn, imx, imy, imz
+    use equations  , only : idn, imx, imy, imz
 
     implicit none
 
diff --git a/src/integrals.F90 b/src/integrals.F90
index 7f45392..20ff51d 100644
--- a/src/integrals.F90
+++ b/src/integrals.F90
@@ -146,12 +146,9 @@ module integrals
 #ifdef MPI
     use mpitools , only : reduce_sum_real_array
 #endif /* MPI */
-    use variables, only : idn, imx, imy, imz
-#ifdef ADI
-    use variables, only : ien
-#endif /* ADI */
+    use equations, only : idn, imx, imy, imz, ien
 #ifdef MHD
-    use variables, only : ibx, iby, ibz
+    use equations, only : ibx, iby, ibz
 #endif /* MHD */
 
     implicit none
diff --git a/src/interpolations.F90 b/src/interpolations.F90
index 64c1bd5..667cef9 100644
--- a/src/interpolations.F90
+++ b/src/interpolations.F90
@@ -581,7 +581,6 @@ module interpolations
   subroutine divergence(b, db, dx, dy, dz)
 
     use coordinates, only : im, jm, km
-    use variables, only : nvr
 
     implicit none
 
diff --git a/src/io.F90 b/src/io.F90
index a88d644..654bc1a 100644
--- a/src/io.F90
+++ b/src/io.F90
@@ -1268,7 +1268,7 @@ module io
     use hdf5     , only : hid_t, hsize_t
     use hdf5     , only : h5gopen_f, h5gclose_f, h5aget_num_attrs_f            &
                         , h5aopen_idx_f, h5aclose_f, h5aget_name_f
-    use variables, only : nqt
+    use equations, only : nv
 
 ! declare variables
 !
@@ -1296,7 +1296,7 @@ module io
 ! initiate the output vector
 !
     dm(:) = 0
-    dm(2) = nqt
+    dm(2) = nv
 
 ! open the global attributes group
 !
@@ -2467,7 +2467,7 @@ module io
     use error    , only : print_error
     use hdf5     , only : hid_t, hsize_t
     use hdf5     , only : h5gcreate_f, h5gclose_f
-    use variables, only : nqt
+    use equations, only : nv
 
 ! declare variables
 !
@@ -2515,13 +2515,13 @@ module io
         am(1) = get_dblocks()
         cm(1) = get_dblocks()
         dm(1) = get_dblocks()
-        dm(2) = nqt
+        dm(2) = nv
         dm(3) = im
         dm(4) = jm
         dm(5) = km
         qm(1) = get_dblocks()
         qm(2) = NDIMS
-        qm(3) = nqt
+        qm(3) = nv
         qm(4) = im
         qm(5) = jm
         qm(6) = km
@@ -2901,15 +2901,15 @@ module io
     use error        , only : print_error
     use hdf5         , only : hid_t, hsize_t
     use hdf5         , only : h5gcreate_f, h5gclose_f
-    use variables    , only : nvr, nqt
-    use variables    , only : idn, imx, imy, imz, ivx, ivy, ivz
+    use equations    , only : nv
+    use equations    , only : idn, imx, imy, imz, ivx, ivy, ivz
 #ifdef ADI
-    use variables    , only : ien, ipr
+    use equations    , only : ien, ipr
 #endif /* ADI */
 #ifdef MHD
-    use variables    , only : ibx, iby, ibz
+    use equations    , only : ibx, iby, ibz
 #ifdef GLM
-    use variables    , only : iph
+    use equations    , only : ibp
 #endif /* GLM */
 #endif /* MHD */
 
@@ -3015,7 +3015,7 @@ module io
           magy(l,1:im,1:jm,1:km) = pdata%u(iby,1:im,1:jm,1:km)
           magz(l,1:im,1:jm,1:km) = pdata%u(ibz,1:im,1:jm,1:km)
 #ifdef GLM
-          bpsi(l,1:im,1:jm,1:km) = pdata%q(iph,1:im,1:jm,1:km)
+          bpsi(l,1:im,1:jm,1:km) = pdata%q(ibp,1:im,1:jm,1:km)
 #endif /* GLM */
 #endif /* MHD */
 
@@ -3123,15 +3123,15 @@ module io
 #ifdef DEBUG
     use refinement   , only : check_refinement_criterion
 #endif /* DEBUG */
-    use variables    , only : nqt
-    use variables    , only : idn, ivx, ivy, ivz
+    use equations    , only : nv
+    use equations    , only : idn, ivx, ivy, ivz
 #ifdef ADI
-    use variables    , only : ipr
+    use equations    , only : ipr
 #endif /* ADI */
 #ifdef MHD
-    use variables    , only : ibx, iby, ibz
+    use equations    , only : ibx, iby, ibz
 #ifdef GLM
-    use variables    , only : iph
+    use equations    , only : ibp
 #endif /* GLM */
 #endif /* MHD */
 
@@ -3237,7 +3237,7 @@ module io
         magy(l,1:in,1:jn,1:kn) = real(pdata%q(iby,ib:ie,jb:je,kb:ke),kind=4)
         magz(l,1:in,1:jn,1:kn) = real(pdata%q(ibz,ib:ie,jb:je,kb:ke),kind=4)
 #ifdef GLM
-        bpsi(l,1:in,1:jn,1:kn) = real(pdata%q(iph,ib:ie,jb:je,kb:ke),kind=4)
+        bpsi(l,1:in,1:jn,1:kn) = real(pdata%q(ibp,ib:ie,jb:je,kb:ke),kind=4)
 #endif /* GLM */
 #endif /* MHD */
 #ifdef DEBUG
diff --git a/src/makefile b/src/makefile
index 119559e..d7cd855 100644
--- a/src/makefile
+++ b/src/makefile
@@ -162,11 +162,11 @@ sources = blocks.F90 boundaries.F90 constants.F90 coordinates.F90 domains.F90  \
           driver.F90 equations.F90 error.F90 evolution.F90 forcing.F90         \
           integrals.F90 interpolations.F90 io.F90 mesh.F90 mpitools.F90        \
           parameters.F90 problems.F90 random.F90 refinement.F90 schemes.F90    \
-          timers.F90 variables.F90
+          timers.F90
 objects = blocks.o boundaries.o constants.o coordinates.o domains.o driver.o   \
           equations.o error.o evolution.o forcing.o integrals.o                \
           interpolations.o io.o mesh.o mpitools.o parameters.o problems.o      \
-          random.o refinement.o schemes.o timers.o variables.o
+          random.o refinement.o schemes.o timers.o
 files   = $(sources) makefile make.default license.txt hosts
 
 $(name).x: $(objects)
@@ -201,40 +201,35 @@ clean-all: clean-bak clean-data clean-exec clean-logs clean-modules            \
 #-------------------------------------------------------------------------------
 
 blocks.o         : blocks.F90 error.o
-boundaries.o     : boundaries.F90 blocks.o coordinates.o error.o               \
-                   interpolations.o mpitools.o variables.o
+boundaries.o     : boundaries.F90 blocks.o coordinates.o equations.o error.o   \
+                   interpolations.o mpitools.o
 constants.o      : constants.F90
 coordinates.o    : coordinates.F90 parameters.o
 driver.o         : driver.F90 blocks.o coordinates.o equations.o evolution.o   \
                    forcing.o integrals.o interpolations.o io.o mesh.o          \
                    mpitools.o parameters.o problems.o random.o refinement.o
-equations.o      : equations.F90 coordinates.o parameters.o variables.o
+equations.o      : equations.F90 coordinates.o parameters.o
 error.o          : error.F90
 evolution.o      : evolution.F90 blocks.o boundaries.o coordinates.o           \
-                   equations.o mesh.o mpitools.o parameters.o schemes.o        \
-                   variables.o
+                   equations.o mesh.o mpitools.o parameters.o schemes.o
 domains.o        : domains.F90 blocks.o boundaries.o coordinates.o parameters.o
-forcing.o        : forcing.F90 constants.o coordinates.o mpitools.o            \
-                   parameters.o random.o variables.o
-integrals.o      : integrals.F90 blocks.o coordinates.o evolution.o mpitools.o \
-                   variables.o
-interpolations.o : interpolations.F90 blocks.o coordinates.o parameters.o      \
-                   variables.o
-io.o             : io.F90 blocks.o coordinates.o error.o evolution.o           \
-                   mpitools.o random.o refinement.o variables.o
-mesh.o           : mesh.F90 blocks.o coordinates.o domains.o error.o           \
-                   interpolations.o mpitools.o problems.o refinement.o         \
-                   variables.o
+forcing.o        : forcing.F90 constants.o coordinates.o equations.o           \
+                   mpitools.o parameters.o random.o
+integrals.o      : integrals.F90 blocks.o coordinates.o equations.o \
+                   evolution.o mpitools.o
+interpolations.o : interpolations.F90 blocks.o coordinates.o parameters.o
+io.o             : io.F90 blocks.o coordinates.o equations.o error.o           \
+                   evolution.o mpitools.o random.o refinement.o
+mesh.o           : mesh.F90 blocks.o coordinates.o domains.o equations.o       \
+                   error.o interpolations.o mpitools.o problems.o refinement.o
 mpitools.o       : mpitools.F90 timers.o
 parameters.o     : parameters.F90 mpitools.o
 problems.o       : problems.F90 blocks.o coordinates.o equations.o error.o     \
-                   parameters.o variables.o
-refinement.o     : refinement.F90 blocks.o coordinates.o parameters.o          \
-                   variables.o
-schemes.o        : schemes.F90 coordinates.o equations.o interpolations.o      \
-                   variables.o
+                   parameters.o
+refinement.o     : refinement.F90 blocks.o coordinates.o equations.o           \
+                   parameters.o
+schemes.o        : schemes.F90 coordinates.o equations.o interpolations.o
 random.o         : random.F90 mpitools.o parameters.o
 timers.o         : timers.F90
-variables.o      : variables.F90
 
 #-------------------------------------------------------------------------------
diff --git a/src/mesh.F90 b/src/mesh.F90
index cdb7284..66259f1 100644
--- a/src/mesh.F90
+++ b/src/mesh.F90
@@ -55,7 +55,7 @@ module mesh
     use coordinates, only : xmin, xmax, ymin, ymax, zmin, zmax
     use coordinates, only : toplev, im, jm, km
     use mpitools   , only : master, nprocs
-    use variables  , only : nqt, nvr
+    use equations  , only : nv
 
 ! local variables are not implicit by default
 !
@@ -78,7 +78,7 @@ module mesh
 !
 ! set data block dimensions
 !
-    call datablock_set_dims(nqt, nvr, im, jm, km)
+    call datablock_set_dims(nv, nv, im, jm, km)
 
 ! print general information about resolutions
 !
@@ -439,7 +439,7 @@ module mesh
     use mpitools , only : master, nprocs, nproc
 #endif /* MPI */
     use refinement , only : check_refinement_criterion
-    use variables, only : nqt
+    use equations  , only : nv
 
     implicit none
 
@@ -464,7 +464,7 @@ module mesh
 
 ! local buffer for data block exchange
 !
-    real(kind=8)   , dimension(nqt,im,jm,km)   :: rbuf
+    real(kind=8)   , dimension(nv,im,jm,km)     :: rbuf
 #endif /* MPI */
 
 ! local pointers
@@ -912,7 +912,7 @@ module mesh
     use coordinates, only : im, jm, km
     use mpitools , only : send_real_array, receive_real_array
     use mpitools , only : nprocs, nproc
-    use variables, only : nqt
+    use equations, only : nv
 
     implicit none
 
@@ -931,7 +931,7 @@ module mesh
 
 ! local buffer for data block exchange
 !
-    real(kind=8)   , dimension(nqt,im,jm,km)   :: rbuf
+    real(kind=8)   , dimension(nv,im,jm,km)   :: rbuf
 
 ! local pointers
 !
@@ -1045,7 +1045,7 @@ module mesh
     use coordinates   , only : ng, nh, in, jn, kn, im, jm, km
     use coordinates   , only : ib, ie, jb, je, kb, ke
     use interpolations, only : minmod3
-    use variables     , only : nqt
+    use equations     , only : nv
 
     implicit none
 
@@ -1088,7 +1088,7 @@ module mesh
 
 ! allocate array for the expanded arrays
 !
-    allocate(u(nqt, dm(1), dm(2), dm(3)))
+    allocate(u(nv, dm(1), dm(2), dm(3)))
 
 ! prepare indices
 !
@@ -1119,7 +1119,7 @@ module mesh
           ic = 2 * (i - il) + 1
           ip = ic + 1
 
-          do p = 1, nqt
+          do p = 1, nv
 
             dul = pdata%u(p,i  ,j,k) - pdata%u(p,i-1,j,k)
             dur = pdata%u(p,i+1,j,k) - pdata%u(p,i  ,j,k)
@@ -1190,10 +1190,10 @@ module mesh
 ! copy data to the current child
 !
 #if NDIMS == 2
-      pchild%data%u(1:nqt,1:im,1:jm,1:km) = u(1:nqt,il:iu,jl:ju, 1:km)
+      pchild%data%u(1:nv,1:im,1:jm,1:km) = u(1:nv,il:iu,jl:ju, 1:km)
 #endif /* NDIMS == 2 */
 #if NDIMS == 3
-      pchild%data%u(1:nqt,1:im,1:jm,1:km) = u(1:nqt,il:iu,jl:ju,kl:ku)
+      pchild%data%u(1:nv,1:im,1:jm,1:km) = u(1:nv,il:iu,jl:ju,kl:ku)
 #endif /* NDIMS == 3 */
 
     end do
@@ -1217,11 +1217,11 @@ module mesh
     use blocks     , only : block_meta, block_data, nchild
     use coordinates, only : ng, nh, in, jn, kn, im, jm, km
     use coordinates, only : ih, jh, kh, ib, jb, kb, ie, je, ke
-    use variables  , only : nfl
+    use equations  , only : nv
 #ifdef MHD
-    use variables  , only : ibx, iby, ibz
+    use equations  , only : ibx, iby, ibz
 #ifdef GLM
-    use variables  , only : iph
+    use equations  , only : ibp
 #endif /* GLM */
 #endif /* MHD */
 
@@ -1310,11 +1310,11 @@ module mesh
 ! the parent block
 !
 #if NDIMS == 2
-      pparent%u(1:nfl,is:it,js:jt,1) =                                         &
-                                   0.25 * (pchild%u(1:nfl,il:iu:2,jl:ju:2,1)   &
-                                         + pchild%u(1:nfl,ip:iu:2,jl:ju:2,1)   &
-                                         + pchild%u(1:nfl,il:iu:2,jp:ju:2,1)   &
-                                         + pchild%u(1:nfl,ip:iu:2,jp:ju:2,1))
+      pparent%u(1:nv,is:it,js:jt,1) =                                         &
+                                   0.25 * (pchild%u(1:nv,il:iu:2,jl:ju:2,1)   &
+                                         + pchild%u(1:nv,ip:iu:2,jl:ju:2,1)   &
+                                         + pchild%u(1:nv,il:iu:2,jp:ju:2,1)   &
+                                         + pchild%u(1:nv,ip:iu:2,jp:ju:2,1))
 
 #ifdef MHD
       pparent%u(ibx:ibz,is:it,js:jt,1) =                                       &
@@ -1323,24 +1323,24 @@ module mesh
                                        + pchild%u(ibx:ibz,il:iu:2,jp:ju:2,1)   &
                                        + pchild%u(ibx:ibz,ip:iu:2,jp:ju:2,1))
 #ifdef GLM
-      pparent%u(iph    ,is:it,js:jt,1) =                                       &
-                                 0.25 * (pchild%u(iph    ,il:iu:2,jl:ju:2,1)   &
-                                       + pchild%u(iph    ,ip:iu:2,jl:ju:2,1)   &
-                                       + pchild%u(iph    ,il:iu:2,jp:ju:2,1)   &
-                                       + pchild%u(iph    ,ip:iu:2,jp:ju:2,1))
+      pparent%u(ibp    ,is:it,js:jt,1) =                                       &
+                                 0.25 * (pchild%u(ibp    ,il:iu:2,jl:ju:2,1)   &
+                                       + pchild%u(ibp    ,ip:iu:2,jl:ju:2,1)   &
+                                       + pchild%u(ibp    ,il:iu:2,jp:ju:2,1)   &
+                                       + pchild%u(ibp    ,ip:iu:2,jp:ju:2,1))
 #endif /* GLM */
 #endif /* MHD */
 #endif /* NDIMS == 2 */
 #if NDIMS == 3
-      pparent%u(1:nfl,is:it,js:jt,ks:kt) =                                     &
-                             0.125 * (pchild%u(1:nfl,il:iu:2,jl:ju:2,kl:ku:2)  &
-                                    + pchild%u(1:nfl,ip:iu:2,jl:ju:2,kl:ku:2)  &
-                                    + pchild%u(1:nfl,il:iu:2,jp:ju:2,kl:ku:2)  &
-                                    + pchild%u(1:nfl,ip:iu:2,jp:ju:2,kl:ku:2)  &
-                                    + pchild%u(1:nfl,il:iu:2,jl:ju:2,kp:ku:2)  &
-                                    + pchild%u(1:nfl,ip:iu:2,jl:ju:2,kp:ku:2)  &
-                                    + pchild%u(1:nfl,il:iu:2,jp:ju:2,kp:ku:2)  &
-                                    + pchild%u(1:nfl,ip:iu:2,jp:ju:2,kp:ku:2))
+      pparent%u(1:nv,is:it,js:jt,ks:kt) =                                     &
+                             0.125 * (pchild%u(1:nv,il:iu:2,jl:ju:2,kl:ku:2)  &
+                                    + pchild%u(1:nv,ip:iu:2,jl:ju:2,kl:ku:2)  &
+                                    + pchild%u(1:nv,il:iu:2,jp:ju:2,kl:ku:2)  &
+                                    + pchild%u(1:nv,ip:iu:2,jp:ju:2,kl:ku:2)  &
+                                    + pchild%u(1:nv,il:iu:2,jl:ju:2,kp:ku:2)  &
+                                    + pchild%u(1:nv,ip:iu:2,jl:ju:2,kp:ku:2)  &
+                                    + pchild%u(1:nv,il:iu:2,jp:ju:2,kp:ku:2)  &
+                                    + pchild%u(1:nv,ip:iu:2,jp:ju:2,kp:ku:2))
 #ifdef MHD
       pparent%u(ibx:ibz,is:it,js:jt,ks:kt) =                                   &
                            0.125 * (pchild%u(ibx:ibz,il:iu:2,jl:ju:2,kl:ku:2)  &
@@ -1352,15 +1352,15 @@ module mesh
                                   + pchild%u(ibx:ibz,il:iu:2,jp:ju:2,kp:ku:2)  &
                                   + pchild%u(ibx:ibz,ip:iu:2,jp:ju:2,kp:ku:2))
 #ifdef GLM
-      pparent%u(iph    ,is:it,js:jt,ks:kt) =                                   &
-                           0.125 * (pchild%u(iph    ,il:iu:2,jl:ju:2,kl:ku:2)  &
-                                  + pchild%u(iph    ,ip:iu:2,jl:ju:2,kl:ku:2)  &
-                                  + pchild%u(iph    ,il:iu:2,jp:ju:2,kl:ku:2)  &
-                                  + pchild%u(iph    ,ip:iu:2,jp:ju:2,kl:ku:2)  &
-                                  + pchild%u(iph    ,il:iu:2,jl:ju:2,kp:ku:2)  &
-                                  + pchild%u(iph    ,ip:iu:2,jl:ju:2,kp:ku:2)  &
-                                  + pchild%u(iph    ,il:iu:2,jp:ju:2,kp:ku:2)  &
-                                  + pchild%u(iph    ,ip:iu:2,jp:ju:2,kp:ku:2))
+      pparent%u(ibp    ,is:it,js:jt,ks:kt) =                                   &
+                           0.125 * (pchild%u(ibp    ,il:iu:2,jl:ju:2,kl:ku:2)  &
+                                  + pchild%u(ibp    ,ip:iu:2,jl:ju:2,kl:ku:2)  &
+                                  + pchild%u(ibp    ,il:iu:2,jp:ju:2,kl:ku:2)  &
+                                  + pchild%u(ibp    ,ip:iu:2,jp:ju:2,kl:ku:2)  &
+                                  + pchild%u(ibp    ,il:iu:2,jl:ju:2,kp:ku:2)  &
+                                  + pchild%u(ibp    ,ip:iu:2,jl:ju:2,kp:ku:2)  &
+                                  + pchild%u(ibp    ,il:iu:2,jp:ju:2,kp:ku:2)  &
+                                  + pchild%u(ibp    ,ip:iu:2,jp:ju:2,kp:ku:2))
 #endif /* GLM */
 #endif /* MHD */
 #endif /* NDIMS == 3 */
diff --git a/src/problems.F90 b/src/problems.F90
index 3b1889a..fa0f84a 100644
--- a/src/problems.F90
+++ b/src/problems.F90
@@ -163,18 +163,9 @@ module problems
     use coordinates, only : ax, ay, az, adx, ady, adz
     use equations  , only : prim2cons
     use equations  , only : gamma
+    use equations  , only : idn, ivx, ivy, ivz, ipr, ibx, iby, ibz, ibp
     use parameters , only : get_parameter_real
-    use variables  , only : nt
-    use variables  , only : idn, ivx, ivy, ivz
-#ifdef ADI
-    use variables  , only : ipr
-#endif /* ADI */
-#ifdef MHD
-    use variables  , only : ibx, iby, ibz
-#ifdef GLM
-    use variables  , only : iph
-#endif /* GLM */
-#endif /* MHD */
+    use equations  , only : nv
 
 ! local variables are not implicit by default
 !
@@ -207,7 +198,7 @@ module problems
 
 ! local arrays
 !
-    real, dimension(nt,im) :: q, u
+    real, dimension(nv,im) :: q, u
     real, dimension(im)    :: x
     real, dimension(jm)    :: y
     real, dimension(km)    :: z
@@ -287,9 +278,9 @@ module problems
 #ifdef MHD
 ! set the uniform magnetic field
 !
-    q(ibx,:) = bext(inx)
-    q(iby,:) = bext(iny)
-    q(ibz,:) = bext(inz)
+    q(ibx,:) = bext(1)
+    q(iby,:) = bext(2)
+    q(ibz,:) = bext(3)
 #endif /* MHD */
 
 ! iterate over all positions in the YZ plane
@@ -426,15 +417,15 @@ module problems
 
 ! convert the primitive variables to conservative ones
 !
-        call prim2cons(im, q(1:nt,1:im), u(1:nt,1:im))
+        call prim2cons(im, q(1:nv,1:im), u(1:nv,1:im))
 
 ! copy the conserved variables to the current block
 !
-        pdata%u(1:nt,1:im,j,k) = u(1:nt,1:im)
+        pdata%u(1:nv,1:im,j,k) = u(1:nv,1:im)
 
 ! copy the primitive variables to the current block
 !
-        pdata%q(1:nt,1:im,j,k) = q(1:nt,1:im)
+        pdata%q(1:nv,1:im,j,k) = q(1:nv,1:im)
 
       end do
     end do
diff --git a/src/refinement.F90 b/src/refinement.F90
index 35ecb8e..592c164 100644
--- a/src/refinement.F90
+++ b/src/refinement.F90
@@ -112,13 +112,7 @@ module refinement
 !
     use blocks     , only : block_data
     use coordinates, only : ib, jb, kb, ie, je, ke
-    use variables  , only : idn
-#ifdef ADI
-    use variables  , only : ipr
-#endif /* ADI */
-#ifdef MHD
-    use variables  , only : ibx, iby, ibz
-#endif /* MHD */
+    use equations  , only : idn, ipr, ibx, iby, ibz
 
 ! local variables are not implicit by default
 !
diff --git a/src/schemes.F90 b/src/schemes.F90
index 80eadde..99f081b 100644
--- a/src/schemes.F90
+++ b/src/schemes.F90
@@ -75,17 +75,11 @@ module schemes
 ! include external variables
 !
     use coordinates   , only : im, jm, km
-    use variables     , only : idn, ivx, ivy, ivz, imx, imy, imz
-#ifdef ADI
-    use variables     , only : ipr, ien
-#endif /* ADI */
+    use equations     , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
 #ifdef MHD
-    use variables     , only : ibx, iby, ibz
-#ifdef GLM
-    use variables     , only : iph
-#endif /* GLM */
+    use equations     , only : ibx, iby, ibz, ibp
 #endif /* MHD */
-    use variables     , only : nt
+    use equations     , only : nv
 
 ! local variables are not implicit by default
 !
@@ -95,8 +89,8 @@ module schemes
 !
     integer                     , intent(in)    :: idir
     real                        , intent(in)    :: dx
-    real, dimension(nt,im,jm,km), intent(in)    :: q
-    real, dimension(nt,im,jm,km), intent(inout) :: f
+    real, dimension(nv,im,jm,km), intent(in)    :: q
+    real, dimension(nv,im,jm,km), intent(inout) :: f
 
 ! local variables
 !
@@ -104,10 +98,10 @@ module schemes
 
 ! local temporary arrays
 !
-    real, dimension(nt,im)       :: qx, fx
-    real, dimension(nt,jm)       :: qy, fy
+    real, dimension(nv,im)       :: qx, fx
+    real, dimension(nv,jm)       :: qy, fy
 #if NDIMS == 3
-    real, dimension(nt,km)       :: qz, fz
+    real, dimension(nv,km)       :: qz, fz
 #endif /* NDIMS == 3 */
 !
 !-------------------------------------------------------------------------------
@@ -140,13 +134,13 @@ module schemes
           qx(iby,1:im) = q(iby,1:im,j,k)
           qx(ibz,1:im) = q(ibz,1:im,j,k)
 #ifdef GLM
-          qx(iph,1:im) = q(iph,1:im,j,k)
+          qx(ibp,1:im) = q(ibp,1:im,j,k)
 #endif /* GLM */
 #endif /* MHD */
 
 ! call one dimensional Riemann solver in order to obtain numerical fluxes
 !
-          call riemann(im, dx, qx(1:nt,1:im), fx(1:nt,1:im))
+          call riemann(im, dx, qx(1:nv,1:im), fx(1:nv,1:im))
 
 ! update the array of fluxes
 !
@@ -162,7 +156,7 @@ module schemes
           f(iby,1:im,j,k) = fx(iby,1:im)
           f(ibz,1:im,j,k) = fx(ibz,1:im)
 #ifdef GLM
-          f(iph,1:im,j,k) = fx(iph,1:im)
+          f(ibp,1:im,j,k) = fx(ibp,1:im)
 #endif /* GLM */
 #endif /* MHD */
 
@@ -190,13 +184,13 @@ module schemes
           qy(iby,1:jm) = q(ibz,i,1:jm,k)
           qy(ibz,1:jm) = q(ibx,i,1:jm,k)
 #ifdef GLM
-          qy(iph,1:jm) = q(iph,i,1:jm,k)
+          qy(ibp,1:jm) = q(ibp,i,1:jm,k)
 #endif /* GLM */
 #endif /* MHD */
 
 ! call one dimensional Riemann solver in order to obtain numerical fluxes
 !
-          call riemann(jm, dx, qy(1:nt,1:jm), fy(1:nt,1:jm))
+          call riemann(jm, dx, qy(1:nv,1:jm), fy(1:nv,1:jm))
 
 ! update the array of fluxes
 !
@@ -212,7 +206,7 @@ module schemes
           f(iby,i,1:jm,k) = fy(ibx,1:jm)
           f(ibz,i,1:jm,k) = fy(iby,1:jm)
 #ifdef GLM
-          f(iph,i,1:jm,k) = fy(iph,1:jm)
+          f(ibp,i,1:jm,k) = fy(ibp,1:jm)
 #endif /* GLM */
 #endif /* MHD */
 
@@ -241,13 +235,13 @@ module schemes
           qz(iby,1:km) = q(ibx,i,j,1:km)
           qz(ibz,1:km) = q(iby,i,j,1:km)
 #ifdef GLM
-          qz(iph,1:km) = q(iph,i,j,1:km)
+          qz(ibp,1:km) = q(ibp,i,j,1:km)
 #endif /* GLM */
 #endif /* MHD */
 
 ! call one dimensional Riemann solver in order to obtain numerical fluxes
 !
-          call riemann(km, dx, qz(1:nt,1:km), fz(1:nt,1:km))
+          call riemann(km, dx, qz(1:nv,1:km), fz(1:nv,1:km))
 
 ! update the array of fluxes
 !
@@ -263,7 +257,7 @@ module schemes
           f(iby,i,j,1:km) = fz(ibz,1:km)
           f(ibz,i,j,1:km) = fz(ibx,1:km)
 #ifdef GLM
-          f(iph,i,j,k) = fz(iph,k)
+          f(ibp,i,j,k) = fz(ibp,k)
 #endif /* GLM */
 #endif /* MHD */
 
@@ -298,7 +292,7 @@ module schemes
 ! include external variables
 !
     use coordinates   , only : im, jm, km
-    use variables     , only : nfl
+    use equations     , only : nv
 
 ! local variables are not implicit by default
 !
@@ -306,9 +300,9 @@ module schemes
 
 ! subroutine arguments
 !
-    real, dimension(3)                 , intent(in)    :: dh
-    real, dimension(NDIMS,nfl,im,jm,km), intent(in)    :: f
-    real, dimension(      nfl,im,jm,km), intent(inout) :: du
+    real, dimension(3)                , intent(in)    :: dh
+    real, dimension(NDIMS,nv,im,jm,km), intent(in)    :: f
+    real, dimension(      nv,im,jm,km), intent(inout) :: du
 
 ! local variables
 !
@@ -391,11 +385,8 @@ module schemes
 
 ! include external variables
 !
-    use variables     , only : nt
-    use variables     , only : ivx, ivy, ivz
-#ifdef ADI
-    use variables     , only : ien
-#endif /* ADI */
+    use equations     , only : nv
+    use equations     , only : ivx, ivy, ivz, ien
 
 ! local variables are not implicit by default
 !
@@ -405,8 +396,8 @@ module schemes
 !
     integer              , intent(in)  :: n
     real                 , intent(in)  :: h
-    real, dimension(nt,n), intent(in)  :: q
-    real, dimension(nt,n), intent(out) :: f
+    real, dimension(nv,n), intent(in)  :: q
+    real, dimension(nv,n), intent(out) :: f
 
 ! local variables
 !
@@ -415,14 +406,14 @@ module schemes
 
 ! local arrays to store the states
 !
-    real, dimension(nt,n) :: ql, qr, ul, ur, fl, fr
+    real, dimension(nv,n) :: ql, qr, ul, ur, fl, fr
     real, dimension(n)    :: cl, cr
 !
 !-------------------------------------------------------------------------------
 !
 ! reconstruct the left and right states of primitive variables
 !
-    do p = 1, nt
+    do p = 1, nv
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:))
     end do
 
@@ -525,11 +516,8 @@ module schemes
 
 ! include external variables
 !
-    use variables     , only : nt
-    use variables     , only : idn, ivx, ivy, ivz, imx, imy, imz
-#ifdef ADI
-    use variables     , only : ipr, ien
-#endif /* ADI */
+    use equations     , only : nv
+    use equations     , only : idn, ivx, ivy, ivz, imx, imy, imz, ipr, ien
 
 ! local variables are not implicit by default
 !
@@ -539,8 +527,8 @@ module schemes
 !
     integer              , intent(in)  :: n
     real                 , intent(in)  :: h
-    real, dimension(nt,n), intent(in)  :: q
-    real, dimension(nt,n), intent(out) :: f
+    real, dimension(nv,n), intent(in)  :: q
+    real, dimension(nv,n), intent(out) :: f
 
 ! local variables
 !
@@ -551,15 +539,15 @@ module schemes
 
 ! local arrays to store the states
 !
-    real, dimension(nt,n) :: ql, qr, ul, ur, fl, fr
+    real, dimension(nv,n) :: ql, qr, ul, ur, fl, fr
     real, dimension(n)    :: cl, cr
-    real, dimension(nt)   :: q1l, q1r, u1l, u1r
+    real, dimension(nv)   :: q1l, q1r, u1l, u1r
 !
 !-------------------------------------------------------------------------------
 !
 ! reconstruct the left and right states of primitive variables
 !
-    do p = 1, nt
+    do p = 1, nv
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:))
     end do
 
@@ -732,11 +720,8 @@ module schemes
 ! include external variables
 !
     use equations     , only : gamma
-    use variables     , only : nt
-    use variables     , only : idn, ivx, ivy, ivz
-#ifdef ADI
-    use variables     , only : ipr, ien
-#endif /* ADI */
+    use equations     , only : nv
+    use equations     , only : idn, ivx, ivy, ivz, ipr, ien
 
 ! local variables are not implicit by default
 !
@@ -746,8 +731,8 @@ module schemes
 !
     integer              , intent(in)  :: n
     real                 , intent(in)  :: h
-    real, dimension(nt,n), intent(in)  :: q
-    real, dimension(nt,n), intent(out) :: f
+    real, dimension(nv,n), intent(in)  :: q
+    real, dimension(nv,n), intent(out) :: f
 
 ! local variables
 !
@@ -757,16 +742,16 @@ module schemes
 
 ! local arrays to store the states
 !
-    real, dimension(nt,n)  :: ql, qr, ul, ur, fl, fr
+    real, dimension(nv,n)  :: ql, qr, ul, ur, fl, fr
     real, dimension(n)     :: cl, cr
-    real, dimension(nt,nt) :: li, ri
-    real, dimension(nt)    :: qi, ci, et, du
+    real, dimension(nv,nv) :: li, ri
+    real, dimension(nv)    :: qi, ci, et, du
 !
 !-------------------------------------------------------------------------------
 !
 ! reconstruct the left and right states of primitive variables
 !
-    do p = 1, nt
+    do p = 1, nv
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:))
     end do
 
@@ -829,7 +814,7 @@ module schemes
 ! calculate vector (Ur - Ul).L
 !
       et(:) = 0.0d0
-      do p = 1, nt
+      do p = 1, nv
         et(:) = et(:) + du(p) * li(p,:)
       end do
 
@@ -839,7 +824,7 @@ module schemes
 
 ! correct the flux by adding the term abs(lambda) R.(Ur - Ul).L
 !
-      do p = 1, nt
+      do p = 1, nv
         f(:,i) = f(:,i) - 0.5d0 * abs(ci(p)) * et(p) * ri(p,:)
       end do
 
@@ -873,8 +858,8 @@ module schemes
 ! include external variables
 !
     use equations     , only : gammam1
-    use variables     , only : nt
-    use variables     , only : idn, ivx, ivy, ivz, ien
+    use equations     , only : nv
+    use equations     , only : idn, ivx, ivy, ivz, ien
 
 ! local variables are not implicit by default
 !
@@ -882,9 +867,9 @@ module schemes
 
 ! subroutine arguments
 !
-    real, dimension(nt)   , intent(in)    :: q
-    real, dimension(nt)   , intent(inout) :: c
-    real, dimension(nt,nt), intent(inout) :: l, r
+    real, dimension(nv)   , intent(in)    :: q
+    real, dimension(nv)   , intent(inout) :: c
+    real, dimension(nv,nv), intent(inout) :: l, r
 
 ! local variables
 !
@@ -976,8 +961,8 @@ module schemes
 ! include external variables
 !
     use equations     , only : csnd
-    use variables     , only : nt
-    use variables     , only : idn, ivx, ivy, ivz, ien
+    use equations     , only : nv
+    use equations     , only : idn, ivx, ivy, ivz, ien
 
 ! local variables are not implicit by default
 !
@@ -985,9 +970,9 @@ module schemes
 
 ! subroutine arguments
 !
-    real, dimension(nt)   , intent(in)    :: q
-    real, dimension(nt)   , intent(inout) :: c
-    real, dimension(nt,nt), intent(inout) :: l, r
+    real, dimension(nv)   , intent(in)    :: q
+    real, dimension(nv)   , intent(inout) :: c
+    real, dimension(nv,nv), intent(inout) :: l, r
 
 ! local variables
 !
@@ -1080,14 +1065,14 @@ module schemes
 
 ! include external variables
 !
-    use variables     , only : nt
-    use variables     , only : ivx, ivy, ivz
+    use equations     , only : nv
+    use equations     , only : ivx, ivy, ivz
 #ifdef ADI
-    use variables     , only : ien
+    use equations     , only : ien
 #endif /* ADI */
-    use variables     , only : ibx, iby, ibz
+    use equations     , only : ibx, iby, ibz
 #ifdef GLM
-    use variables     , only : iph
+    use equations     , only : ibp
 #endif /* GLM */
 
 ! local variables are not implicit by default
@@ -1098,8 +1083,8 @@ module schemes
 !
     integer              , intent(in)  :: n
     real                 , intent(in)  :: h
-    real, dimension(nt,n), intent(in)  :: q
-    real, dimension(nt,n), intent(out) :: f
+    real, dimension(nv,n), intent(in)  :: q
+    real, dimension(nv,n), intent(out) :: f
 
 ! local variables
 !
@@ -1108,14 +1093,14 @@ module schemes
 
 ! local arrays to store the states
 !
-    real, dimension(nt,n) :: ql, qr, ul, ur, fl, fr
+    real, dimension(nv,n) :: ql, qr, ul, ur, fl, fr
     real, dimension(n)    :: cl, cr
 !
 !-------------------------------------------------------------------------------
 !
 ! reconstruct the left and right states of primitive variables
 !
-    do p = 1, nt
+    do p = 1, nv
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:))
     end do
 
@@ -1128,16 +1113,16 @@ module schemes
 #ifdef GLM
 ! reconstruct the left and right states of the scalar potential
 !
-    call reconstruct(n, h, q(iph,:), ql(iph,:), qr(iph,:))
+    call reconstruct(n, h, q(ibp,:), ql(ibp,:), qr(ibp,:))
 
 ! obtain the state values for Bx and Psi for the GLM-MHD equations
 !
-    cl(:) = 0.5d0 * ((qr(ibx,:) + ql(ibx,:)) - (qr(iph,:) - ql(iph,:)) / cmax)
-    cr(:) = 0.5d0 * ((qr(iph,:) + ql(iph,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax)
+    cl(:) = 0.5d0 * ((qr(ibx,:) + ql(ibx,:)) - (qr(ibp,:) - ql(ibp,:)) / cmax)
+    cr(:) = 0.5d0 * ((qr(ibp,:) + ql(ibp,:)) - (qr(ibx,:) - ql(ibx,:)) * cmax)
     ql(ibx,:) = cl(:)
     qr(ibx,:) = cl(:)
-    ql(iph,:) = cr(:)
-    qr(iph,:) = cr(:)
+    ql(ibp,:) = cr(:)
+    qr(ibp,:) = cr(:)
 #endif /* GLM */
 
 #ifdef FIX_POSITIVITY
@@ -1239,10 +1224,10 @@ module schemes
 !
     use equations     , only : prim2cons, fluxspeed
     use interpolations, only : reconstruct, fix_positivity
-    use variables     , only : nt
-    use variables     , only : idn, imx, imy, imz, ien
-    use variables     , only : ivx, ivy, ivz, ipr
-    use variables     , only : ibx, iby, ibz
+    use equations     , only : nv
+    use equations     , only : idn, imx, imy, imz, ien
+    use equations     , only : ivx, ivy, ivz, ipr
+    use equations     , only : ibx, iby, ibz
 
 ! local variables are not implicit by default
 !
@@ -1252,9 +1237,9 @@ module schemes
 !
     integer                 , intent(in)  :: n
     real                    , intent(in)  :: h
-    real   , dimension(nt,n), intent(in)  :: q
+    real   , dimension(nv,n), intent(in)  :: q
     real   , dimension(   n), intent(in)  :: b
-    real   , dimension(nt,n), intent(out) :: f
+    real   , dimension(nv,n), intent(out) :: f
     logical, dimension(n)   , intent(in)  :: s
 
 ! local variables
@@ -1271,9 +1256,9 @@ module schemes
 
 ! local arrays to store the states
 !
-    real, dimension(nt,n) :: ul, ur, ql, qr, fl, fr
+    real, dimension(nv,n) :: ul, ur, ql, qr, fl, fr
     real, dimension(n)    :: cl, cr
-    real, dimension(nt)   :: q1l, q1r, u1l, u1r
+    real, dimension(nv)   :: q1l, q1r, u1l, u1r
 !
 !-------------------------------------------------------------------------------
 !
@@ -1282,7 +1267,7 @@ module schemes
     do p = 1, ibx - 1
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:), s(:))
     end do
-    do p = ibx + 1, nt
+    do p = ibx + 1, nv
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:), s(:))
     end do
 
@@ -1809,10 +1794,10 @@ module schemes
     use equations     , only : gamma
     use equations     , only : prim2cons, fluxspeed
     use interpolations, only : reconstruct, fix_positivity
-    use variables     , only : nt
-    use variables     , only : idn, imx, imy, imz, ien
-    use variables     , only : ivx, ivy, ivz, ipr
-    use variables     , only : ibx, iby, ibz
+    use equations     , only : nv
+    use equations     , only : idn, imx, imy, imz, ien
+    use equations     , only : ivx, ivy, ivz, ipr
+    use equations     , only : ibx, iby, ibz
 
 ! local variables are not implicit by default
 !
@@ -1822,9 +1807,9 @@ module schemes
 !
     integer                 , intent(in)  :: n
     real                    , intent(in)  :: h
-    real   , dimension(nt,n), intent(in)  :: q
+    real   , dimension(nv,n), intent(in)  :: q
     real   , dimension(   n), intent(in)  :: b
-    real   , dimension(nt,n), intent(out) :: f
+    real   , dimension(nv,n), intent(out) :: f
     logical, dimension(n)   , intent(in)  :: s
 
 ! local variables
@@ -1841,9 +1826,9 @@ module schemes
 
 ! local arrays to store the states
 !
-    real, dimension(nt,n) :: ul, ur, ql, qr, fl, fr
+    real, dimension(nv,n) :: ul, ur, ql, qr, fl, fr
     real, dimension(n)    :: cl, cr
-    real, dimension(nt)   :: q1l, q1r, u1l, u1r, q2, u2
+    real, dimension(nv)   :: q1l, q1r, u1l, u1r, q2, u2
 !
 !-------------------------------------------------------------------------------
 !
@@ -1852,7 +1837,7 @@ module schemes
     do p = 1, ibx - 1
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:), s(:))
     end do
-    do p = ibx + 1, nt
+    do p = ibx + 1, nv
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:), s(:))
     end do
 
@@ -2289,10 +2274,10 @@ module schemes
     use equations     , only : csnd
     use equations     , only : prim2cons, fluxspeed
     use interpolations, only : reconstruct, fix_positivity
-    use variables     , only : nt
-    use variables     , only : idn, imx, imy, imz
-    use variables     , only : ivx, ivy, ivz
-    use variables     , only : ibx, iby, ibz
+    use equations     , only : nv
+    use equations     , only : idn, imx, imy, imz
+    use equations     , only : ivx, ivy, ivz
+    use equations     , only : ibx, iby, ibz
 
 ! local variables are not implicit by default
 !
@@ -2302,9 +2287,9 @@ module schemes
 !
     integer                 , intent(in)  :: n
     real                    , intent(in)  :: h
-    real   , dimension(nt,n), intent(in)  :: q
+    real   , dimension(nv,n), intent(in)  :: q
     real   , dimension(   n), intent(in)  :: b
-    real   , dimension(nt,n), intent(out) :: f
+    real   , dimension(nv,n), intent(out) :: f
     logical, dimension(n)   , intent(in)  :: s
 
 ! local variables
@@ -2318,9 +2303,9 @@ module schemes
 
 ! local arrays to store the states
 !
-    real, dimension(nt,n) :: ul, ur, ql, qr, fl, fr
+    real, dimension(nv,n) :: ul, ur, ql, qr, fl, fr
     real, dimension(n)    :: cl, cr
-    real, dimension(nt)   :: u1l, u1r, u2
+    real, dimension(nv)   :: u1l, u1r, u2
 !
 !-------------------------------------------------------------------------------
 !
@@ -2329,7 +2314,7 @@ module schemes
     do p = 1, ibx - 1
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:), s(:))
     end do
-    do p = ibx + 1, nt
+    do p = ibx + 1, nv
       call reconstruct(n, h, q(p,:), ql(p,:), qr(p,:), s(:))
     end do
 
diff --git a/src/variables.F90 b/src/variables.F90
index 81037b9..c92928d 100644
--- a/src/variables.F90
+++ b/src/variables.F90
@@ -28,33 +28,6 @@
 module variables
 
   implicit none
-
-! parameters
-!
-  integer(kind=4), parameter :: inx =  1, iny =  2, inz =  3
-  integer(kind=4), parameter :: idn =  1, imx =  2, imy =  3, imz =  4 &
-                                        , ivx =  2, ivy =  3, ivz =  4
-#ifdef ISO
-  integer(kind=4), parameter :: nfl =  4
-#endif /* ISO */
-#ifdef ADI
-  integer(kind=4), parameter :: ien =  5, ipr =  5
-  integer(kind=4), parameter :: nfl =  5
-#endif /* ADI */
-#ifdef HYDRO
-  integer(kind=4), parameter :: nqt = nfl
-#endif /* HYDRO */
-#ifdef MHD
-  integer(kind=4), parameter :: ibx = nfl + 1, iby = ibx + 1, ibz = iby + 1
-#ifdef GLM
-  integer(kind=4), parameter :: iph =  ibz + 1
-  integer(kind=4), parameter :: nqt = iph
-#else /* GLM */
-  integer(kind=4), parameter :: nqt = ibz
-#endif /* GLM */
-#endif /* MHD */
-  integer(kind=4), parameter :: nvr = nqt
-  integer(kind=4), parameter :: nt  = nqt
 !
 !===============================================================================
 !