diff --git a/.github/workflows/amrclaw_tests.yml b/.github/workflows/amrclaw_tests.yml
new file mode 100644
index 0000000..c1e402a
--- /dev/null
+++ b/.github/workflows/amrclaw_tests.yml
@@ -0,0 +1,66 @@
+name: Test Riemann solvers with AMRClaw tests
+
+on:
+  push:
+    branches: [ "implicit_none" ]
+  pull_request:
+    branches: [ "implicit_none" ]
+
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+env:
+  CLAW: ${{ github.workspace }}
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+    - name: Set up Python 3.10
+      uses: actions/setup-python@v5
+      with:
+        python-version: "3.10"
+
+    - name: Install dependencies
+      run: |
+        sudo apt-get update
+        sudo apt-get install gfortran
+
+        python -m pip install --upgrade pip
+        pip install 'numpy<2.0'
+        pip install matplotlib #Some imports require matplotlib
+        pip install scipy #To not skip tests
+        pip install nose
+        pip install flake8 meson-python ninja pytest 
+
+    - name: Checkout Clawpack
+      uses: actions/checkout@v4.1.5
+      with:
+        repository: clawpack/clawpack
+        submodules: true
+
+    - name: Checkout AMRClaw branch
+      uses: actions/checkout@v4.1.5
+      with:
+        repository: clawpack/amrclaw  #
+        path: amrclaw
+        ref: master
+    
+    - name: Checkout implicit_none from this repo
+      uses: actions/checkout@v4.1.5
+      with:
+        repository: ${{ github.repository }}
+        path: ${{ env.CLAW }}/riemann
+        ref: implicit_none
+
+    - name: Install Clawpack
+      run: |
+        cd ${CLAW}
+        pip install --no-build-isolation --editable .
+
+    - name: Test with pytest
+      run: |
+        cd ${CLAW}/riemann
+        bash collect_pytests_amrclaw.sh
diff --git a/.github/workflows/classic_tests.yml b/.github/workflows/classic_tests.yml
new file mode 100644
index 0000000..1c3c75d
--- /dev/null
+++ b/.github/workflows/classic_tests.yml
@@ -0,0 +1,66 @@
+name: Test Riemann solvers with Classic Clawpack tests
+
+on:
+  push:
+    branches: [ "implicit_none" ]
+  pull_request:
+    branches: [ "implicit_none" ]
+
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+env:
+  CLAW: ${{ github.workspace }}
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+    - name: Set up Python 3.10
+      uses: actions/setup-python@v5
+      with:
+        python-version: "3.10"
+
+    - name: Install dependencies
+      run: |
+        sudo apt-get update
+        sudo apt-get install gfortran
+
+        python -m pip install --upgrade pip
+        pip install 'numpy<2.0'
+        pip install matplotlib #Some imports require matplotlib
+        pip install scipy #To not skip tests
+        pip install nose
+        pip install flake8 meson-python ninja pytest 
+
+    - name: Checkout Clawpack
+      uses: actions/checkout@v4.1.5
+      with:
+        repository: clawpack/clawpack
+        submodules: true
+
+    - name: Checkout Classic branch
+      uses: actions/checkout@v4.1.5
+      with:
+        repository: clawpack/classic  #
+        path: classic
+        ref: master
+    
+    - name: Checkout implicit_none from this repo
+      uses: actions/checkout@v4.1.5
+      with:
+        repository: ${{ github.repository }}
+        path: ${{ env.CLAW }}/riemann
+        ref: implicit_none
+
+    - name: Install Clawpack
+      run: |
+        cd ${CLAW}
+        pip install --no-build-isolation --editable .
+
+    - name: Test with pytest
+      run: |
+        cd ${CLAW}/classic
+        pytest tests
diff --git a/.github/workflows/pyclaw_tests.yml b/.github/workflows/pyclaw_tests.yml
new file mode 100644
index 0000000..40f0a01
--- /dev/null
+++ b/.github/workflows/pyclaw_tests.yml
@@ -0,0 +1,65 @@
+name: Test Riemann solvers with PyClaw tests
+
+on:
+  push:
+    branches: [ "implicit_none" ]
+  pull_request:
+    branches: [ "implicit_none" ]
+
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+env:
+  CLAW: ${{ github.workspace }}
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+    - name: Set up Python 3.10
+      uses: actions/setup-python@v5
+      with:
+        python-version: "3.10"
+
+    - name: Install dependencies
+      run: |
+        sudo apt-get update
+        sudo apt-get install gfortran
+
+        python -m pip install --upgrade pip
+        pip install 'numpy<2.0'
+        pip install matplotlib #Some imports require matplotlib
+        pip install scipy #To not skip tests
+        pip install flake8 meson-python ninja pytest
+
+    - name: Checkout Clawpack
+      uses: actions/checkout@v4.1.5
+      with:
+        repository: clawpack/clawpack
+        submodules: true
+
+    - name: Checkout PyClaw branch
+      uses: actions/checkout@v4.1.5
+      with:
+        repository: clawpack/pyclaw  #
+        path: pyclaw
+        ref: master
+    
+    - name: Checkout implicit_none from this repo
+      uses: actions/checkout@v4.1.5
+      with:
+        repository: ${{ github.repository }}
+        path: ${{ env.CLAW }}/riemann
+        ref: implicit_none
+
+    - name: Install Clawpack
+      run: |
+        cd ${CLAW}
+        pip install --no-build-isolation --editable .
+
+    - name: Test with pytest
+      run: |
+        cd ${CLAW}/pyclaw
+        pytest --ignore=development --ignore=examples/shallow_sphere
diff --git a/collect_pytests_amrclaw.sh b/collect_pytests_amrclaw.sh
new file mode 100644
index 0000000..a3e1f2c
--- /dev/null
+++ b/collect_pytests_amrclaw.sh
@@ -0,0 +1,19 @@
+#Collect Pytest tests for AMRClaw
+cd ${CLAW}/amrclaw/tests
+pathlist=()
+for dir in $(ls -d */); do
+    cd $dir
+    #Check if there are any test files
+    #Capture the error message and check if it is empty
+    if [ -z "$(ls *test*.py 2>/dev/null )" ]; then
+        cd ..
+        continue
+    fi
+    for file in $(ls *test*.py); do
+        pathlist+=($(pwd)/$file)
+    done
+    cd ..
+done
+
+#Run the tests
+pytest ${pathlist[@]}
\ No newline at end of file
diff --git a/src/rp1_acoustics.f90 b/src/rp1_acoustics.f90
index cf58cb4..a082090 100644
--- a/src/rp1_acoustics.f90
+++ b/src/rp1_acoustics.f90
@@ -29,23 +29,23 @@ subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 ! From the basic clawpack routines, this routine is called with ql = qr
 
 
-    implicit double precision (a-h,o-z)
+    implicit none
 
-    dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves,1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension apdq(meqn, 1-mbc:maxm+mbc)
-    dimension amdq(meqn, 1-mbc:maxm+mbc)
+    integer, intent(in) :: maxm, meqn, mwaves, mbc, mx, maux
+    real(kind=8), intent(in) :: ql(meqn, 1-mbc:maxm+mbc),qr(meqn, 1-mbc:maxm+mbc)
 
-!     local arrays
-!     ------------
-    dimension delta(2)
+    real(kind=8), intent(out)  :: wave(meqn, mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: apdq(meqn, 1-mbc:maxm+mbc), amdq(meqn, 1-mbc:maxm+mbc)
 
-!     # density, bulk modulus, and sound speed, and impedence of medium:
-!     # (should be set in setprob.f)
-    common /cparam/ rho,bulk,cc,zz
+    real(kind=8) :: delta(2)
+    real(kind=8) :: rho, bulk, cc, zz
+    real(kind=8) :: a1, a2, auxl, auxr
+    integer :: i, m
 
+    ! density, bulk modulus, and sound speed, and impedence of medium:
+    ! (should be set in setprob.f)
+    common /cparam/ rho,bulk,cc,zz
 
 !     # split the jump in q at each interface into waves
 
diff --git a/src/rp1_acoustics_adjoint.f90 b/src/rp1_acoustics_adjoint.f90
index a5b5343..a23c326 100644
--- a/src/rp1_acoustics_adjoint.f90
+++ b/src/rp1_acoustics_adjoint.f90
@@ -25,20 +25,26 @@ subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq)
 
 !     # aux arrays not used in this solver.
 
-    implicit double precision (a-h,o-z)
+    implicit none
 
-    dimension fwave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension apdq(meqn, 1-mbc:maxm+mbc)
-    dimension amdq(meqn, 1-mbc:maxm+mbc)
-    dimension auxl(maux, 1-mbc:maxm+mbc)
-    dimension auxr(maux, 1-mbc:maxm+mbc)
+    !Input
+    integer, intent(in) :: maxm, meqn, mwaves, maux, mbc, mx
+    double precision, intent(in) :: ql(meqn, 1-mbc:maxm+mbc)
+    double precision, intent(in) :: qr(meqn, 1-mbc:maxm+mbc)
+    double precision, intent(in) :: auxl(maux, 1-mbc:maxm+mbc)
+    double precision, intent(in) :: auxr(maux, 1-mbc:maxm+mbc)
+
+    double precision, intent(out) :: fwave(meqn, mwaves, 1-mbc:maxm+mbc)
+    double precision, intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    double precision, intent(out) :: amdq(meqn, 1-mbc:maxm+mbc)
+    double precision, intent(out) :: apdq(meqn, 1-mbc:maxm+mbc)
+
+    !local
+    integer :: i, m
+    double precision :: rhoi, bulki, beta1, beta2
+    double precision :: rho, bulk, cc, zz
+    double precision, dimension(2) :: delta
 
-!     local arrays
-!     ------------
-    dimension delta(2)
 
 !     # density, bulk modulus, and sound speed, and impedence of medium:
 !     # (should be set in setprob.f)
@@ -47,7 +53,7 @@ subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq)
 !     # split the jump in f(q) at each interface into waves
 
 !     # find b1 and b2, the coefficients of the 2 eigenvectors:
-    do 20 i = 2-mbc, mx+mbc
+    do i = 2-mbc, mx+mbc
 !       # material properties
         rhoi = zz/cc
         bulki = rhoi*cc**2
@@ -70,7 +76,7 @@ subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq)
         fwave(2,2,i) = -beta2*zz
         s(2,i) = cc
     
-    20 END DO
+    END DO
 
 
 !     # compute the leftgoing and rightgoing flux differences:
diff --git a/src/rp1_acoustics_variable.f90 b/src/rp1_acoustics_variable.f90
index 7684fa5..5ad7b05 100644
--- a/src/rp1_acoustics_variable.f90
+++ b/src/rp1_acoustics_variable.f90
@@ -36,20 +36,21 @@ subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 ! From the basic clawpack routines, this routine is called with ql = qr
 
 
-    implicit double precision (a-h,o-z)
-
-    dimension auxl(maux, 1-mbc:maxm+mbc)
-    dimension auxr(maux, 1-mbc:maxm+mbc)
-    dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension apdq(meqn, 1-mbc:maxm+mbc)
-    dimension amdq(meqn, 1-mbc:maxm+mbc)
-
-!    local arrays
-!    ------------
-    dimension delta(2)
+    implicit none
+
+    integer, intent(in) :: maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn, 1-mbc:mx+mbc)
+    real(kind=8), intent(in) :: qr(meqn, 1-mbc:mx+mbc)
+    real(kind=8), intent(in) :: auxl(maux, 1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: auxr(maux, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: apdq(meqn, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: amdq(meqn, 1-mbc:maxm+mbc)
+
+    real(kind=8) :: delta(2)
+    real(kind=8) :: zi, zim, a1, a2
+    integer :: i, m
 
 
     ! split the jump in q at each interface into waves
diff --git a/src/rp1_acoustics_variable_adjoint.f90 b/src/rp1_acoustics_variable_adjoint.f90
index a9ff6dd..bf008aa 100644
--- a/src/rp1_acoustics_variable_adjoint.f90
+++ b/src/rp1_acoustics_variable_adjoint.f90
@@ -8,26 +8,29 @@ subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq)
 !     # auxl(1,i) should contain the impedance Z in cell i
 !     # auxl(2,i) should contain the sound speed c in cell i
 
-    implicit double precision (a-h,o-z)
+    implicit none
 
-    dimension auxl(maux, 1-mbc:maxm+mbc)
-    dimension auxr(maux, 1-mbc:maxm+mbc)
-    dimension fwave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension apdq(meqn, 1-mbc:maxm+mbc)
-    dimension amdq(meqn, 1-mbc:maxm+mbc)
+    !Input
+    integer, intent(in) :: maxm, meqn, mwaves, maux, mbc, mx
+    double precision, intent(in) :: ql(meqn, 1-mbc:maxm+mbc)
+    double precision, intent(in) :: qr(meqn, 1-mbc:maxm+mbc)
+    double precision, intent(in) :: auxl(maux, 1-mbc:maxm+mbc)
+    double precision, intent(in) :: auxr(maux, 1-mbc:maxm+mbc)
 
-!     local arrays
-!     ------------
-    dimension delta(2)
+    double precision, intent(out) :: fwave(meqn, mwaves, 1-mbc:maxm+mbc)
+    double precision, intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    double precision, intent(out) :: amdq(meqn, 1-mbc:maxm+mbc)
+    double precision, intent(out) :: apdq(meqn, 1-mbc:maxm+mbc)
 
+    !Local
+    integer :: i, m
+    double precision :: zi, zim, ci, cim, rhoi, rhoim, bulki, bulkim, beta1, beta2
+    double precision :: delta(2)
 
 !   # split the jump in f(q) at each interface into f-waves
 
 !   # find a1 and a2, the coefficients of the 2 eigenvectors:
-    do 20 i = 2-mbc, mx+mbc
+    do i = 2-mbc, mx+mbc
 !       # material properties
         zi = auxl(1,i)
         zim = auxr(1,i-1)
@@ -55,17 +58,18 @@ subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq)
         fwave(2,2,i) = -beta2*zi
         s(2,i) = ci
     
-    20 END DO
+    END DO
 
 !     # compute the leftgoing and rightgoing fluctuations:
 !     # Note s(1,i) < 0   and   s(2,i) > 0.
 
 !   # f-wave spitting, so do not multiply by wave speeds:
-    do 220 m=1,meqn
-        do 220 i = 2-mbc, mx+mbc
+    do m=1,meqn
+        do i = 2-mbc, mx+mbc
             amdq(m,i) = fwave(m,1,i)
             apdq(m,i) = fwave(m,2,i)
-    220 END DO
+        end do
+    END DO
 
     return
     end subroutine rp1
diff --git a/src/rp1_advection.f90 b/src/rp1_advection.f90
index b8b8a8c..de97b9f 100644
--- a/src/rp1_advection.f90
+++ b/src/rp1_advection.f90
@@ -23,17 +23,24 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 !                                    and right state ql(i,:)
 ! From the basic clawpack routine step1, rp is called with ql = qr = q.
 
+    implicit none
 
-    implicit double precision (a-h,o-z)
-    double precision ::   ql( meqn,1-mbc:maxmx+mbc)
-    double precision ::   qr( meqn,1-mbc:maxmx+mbc)
-    double precision ::   amdq( meqn,1-mbc:maxmx+mbc)
-    double precision ::   apdq( meqn,1-mbc:maxmx+mbc)
-    double precision ::    s(mwaves,1-mbc:maxmx+mbc)
-    double precision :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+    integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxmx+mbc)
+
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxmx+mbc)
+
+    real(kind=8) :: u
+    integer :: i
     common /cparam/ u
 
-    do 30 i=2-mbc,mx+mbc
+    do i = 2-mbc,mx+mbc
     
     !        # Compute the wave and speed
     
@@ -41,7 +48,7 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
         s(1,i) = u
         amdq(1,i) = dmin1(u, 0.d0) * wave(1,1,i)
         apdq(1,i) = dmax1(u, 0.d0) * wave(1,1,i)
-    30 END DO
+    end do
 
     return
     end subroutine rp1
diff --git a/src/rp1_advection_adjoint.f90 b/src/rp1_advection_adjoint.f90
index 2e4b73d..ba1fb88 100644
--- a/src/rp1_advection_adjoint.f90
+++ b/src/rp1_advection_adjoint.f90
@@ -23,17 +23,25 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 !                                    and right state ql(i,:)
 ! From the basic clawpack routine step1, rp is called with ql = qr = q.
 
-
-    implicit double precision (a-h,o-z)
-    double precision ::   ql( meqn,1-mbc:maxmx+mbc)
-    double precision ::   qr( meqn,1-mbc:maxmx+mbc)
-    double precision ::   amdq( meqn,1-mbc:maxmx+mbc)
-    double precision ::   apdq( meqn,1-mbc:maxmx+mbc)
-    double precision ::    s(mwaves,1-mbc:maxmx+mbc)
-    double precision :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+    implicit none
+    !Input
+    integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
+    double precision, intent(in) ::   ql( meqn,1-mbc:maxmx+mbc)
+    double precision, intent(in) ::   qr( meqn,1-mbc:maxmx+mbc)
+    double precision, intent(in) :: auxl(maux,1-mbc:maxmx+mbc)
+    double precision, intent(in) :: auxr(maux,1-mbc:maxmx+mbc)
+    
+    double precision, intent(out) ::   amdq( meqn,1-mbc:maxmx+mbc)
+    double precision, intent(out) ::   apdq( meqn,1-mbc:maxmx+mbc)
+    double precision, intent(out) ::    s(mwaves,1-mbc:maxmx+mbc)
+    double precision, intent(out) :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+
+    !Local
+    double precision :: u
+    integer :: i
     common /cparam/ u
 
-    do 30 i=2-mbc,mx+mbc
+    do i=2-mbc,mx+mbc
     
     !        # Compute the wave and speed
 !       # (speed negative because we are going backward in time)
@@ -42,7 +50,7 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
         s(1,i) = -u
         amdq(1,i) = dmin1(s(1,i), 0.d0) * wave(1,1,i)
         apdq(1,i) = dmax1(s(1,i), 0.d0) * wave(1,1,i)
-    30 END DO
+    END DO
 
     return
     end subroutine rp1
diff --git a/src/rp1_advection_color.f90 b/src/rp1_advection_color.f90
index eb45e50..05e2de4 100644
--- a/src/rp1_advection_color.f90
+++ b/src/rp1_advection_color.f90
@@ -26,15 +26,21 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 ! From the basic clawpack routine step1, rp is called with ql = qr = q.
 
 
-    implicit double precision (a-h,o-z)
-    dimension   ql(meqn,1-mbc:maxmx+mbc)
-    dimension   qr(meqn,1-mbc:maxmx+mbc)
-    dimension  auxl(maux,1-mbc:maxmx+mbc)
-    dimension  auxr(maux,1-mbc:maxmx+mbc)
-    dimension    s(mwaves,1-mbc:maxmx+mbc)
-    dimension wave(meqn, mwaves,1-mbc:maxmx+mbc)
-    dimension amdq(meqn,1-mbc:maxmx+mbc)
-    dimension apdq(meqn,1-mbc:maxmx+mbc)
+    implicit none
+
+    integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxmx+mbc)
+
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxmx+mbc)
+
+    real(kind=8) :: u
+    integer :: i
     common /comrp/ u
 
 
diff --git a/src/rp1_burgers.f90 b/src/rp1_burgers.f90
index 3e35f9e..387069f 100644
--- a/src/rp1_burgers.f90
+++ b/src/rp1_burgers.f90
@@ -12,16 +12,19 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 ! Conserved quantities:
 !       1 q
     
-    implicit double precision (a-h,o-z)
-
-    integer :: maxmx, meqn, mwaves, mbc, mx
-        
-    double precision :: ql(meqn,1-mbc:maxmx+mbc)
-    double precision :: qr(meqn,1-mbc:maxmx+mbc)
-    double precision :: s(mwaves, 1-mbc:maxmx+mbc)
-    double precision :: wave(meqn, mwaves, 1-mbc:maxmx+mbc)
-    double precision :: amdq(meqn, 1-mbc:maxmx+mbc)
-    double precision :: apdq(meqn, 1-mbc:maxmx+mbc)
+    implicit none
+
+    integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxmx+mbc)
+
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxmx+mbc)
+
     integer :: i
     logical :: efix
  
diff --git a/src/rp1_cubic.f90 b/src/rp1_cubic.f90
index 720f5f6..a6e23c0 100644
--- a/src/rp1_cubic.f90
+++ b/src/rp1_cubic.f90
@@ -23,14 +23,18 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 !                                    and right state ql(:, i)
 ! From the basic clawpack routine step1, rp is called with ql = qr = q.
 
-    implicit none
+     implicit none
 
     integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
-    double precision, dimension(meqn, 1-mbc:maxmx+mbc),         intent(in)  :: ql, qr
-    double precision, dimension(maux, 1-mbc:maxmx+mbc),         intent(in)  :: auxl, auxr
-    double precision, dimension(meqn, mwaves, 1-mbc:maxmx+mbc), intent(out) :: wave
-    double precision, dimension(meqn, 1-mbc:maxmx+mbc),         intent(out) :: amdq, apdq
-    double precision, dimension(mwaves, 1-mbc:maxmx+mbc),       intent(out) :: s
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxmx+mbc)
+
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxmx+mbc)
 
     integer :: i
 
diff --git a/src/rp1_euler_hlle.f90 b/src/rp1_euler_hlle.f90
index bd948d3..774adf1 100644
--- a/src/rp1_euler_hlle.f90
+++ b/src/rp1_euler_hlle.f90
@@ -15,17 +15,21 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 
     implicit none
 
-    integer, intent(in) :: maxmx, meqn, mwaves, mbc, mx, maux
-    double precision, dimension(meqn,1-mbc:maxmx+mbc), intent(in) :: ql, qr
-    double precision, dimension(maux,1-mbc:maxmx+mbc), intent(in) :: auxl, auxr
-    double precision, dimension(meqn, mwaves, 1-mbc:maxmx+mbc), intent(out) :: wave
-    double precision, dimension(meqn, 1-mbc:maxmx+mbc), intent(out) :: amdq, apdq
-    double precision, dimension(mwaves, 1-mbc:maxmx+mbc), intent(out) :: s
-
-    double precision :: u_l, u_r, c_l, c_r, u_hat, c_hat
-    double precision :: p_l, p_r, H_l, H_r, rhsqrt_l, rhsqrt_r
-    double precision :: rhsq2, H_hat, gamma, gamma1, E_l, E_r
-    double precision :: rho_m, rhou_m, E_m, s1, s2
+    integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxmx+mbc)
+
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxmx+mbc)
+
+    real(kind=8) :: u_l, u_r, c_l, c_r, u_hat, c_hat
+    real(kind=8) :: p_l, p_r, H_l, H_r, rhsqrt_l, rhsqrt_r
+    real(kind=8) :: rhsq2, H_hat, gamma, gamma1, E_l, E_r
+    real(kind=8) :: rho_m, rhou_m, E_m, s1, s2
     integer :: m, i, mw
 
     common /cparam/  gamma
diff --git a/src/rp1_euler_with_efix.f90 b/src/rp1_euler_with_efix.f90
index 571b62c..828e348 100644
--- a/src/rp1_euler_with_efix.f90
+++ b/src/rp1_euler_with_efix.f90
@@ -24,21 +24,28 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 !                                    and right state ql(:,i)
 ! From the basic clawpack routine step1, rp is called with ql = qr = q.
 
-
-    implicit double precision (a-h,o-z)
-    dimension   ql(meqn,1-mbc:maxmx+mbc)
-    dimension   qr(meqn,1-mbc:maxmx+mbc)
-    dimension    s(mwaves,1-mbc:maxmx+mbc)
-    dimension wave(meqn, mwaves,1-mbc:maxmx+mbc)
-    dimension amdq(meqn,1-mbc:maxmx+mbc)
-    dimension apdq(meqn,1-mbc:maxmx+mbc)
-
-!     # local storage
-!     ---------------
-    dimension delta(3)
-    dimension u(1-mbc:maxmx+mbc),enth(1-mbc:maxmx+mbc)
-    dimension a(1-mbc:maxmx+mbc)
+    implicit none
+
+    integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxmx+mbc)
+
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxmx+mbc)
+
+    real(kind=8) delta(3)
+    real(kind=8) u(1-mbc:maxmx+mbc),enth(1-mbc:maxmx+mbc)
+    real(kind=8) a(1-mbc:maxmx+mbc)
     logical :: efix
+    real(kind=8) :: gamma, gamma1, rhsqrtl, rhsqrtr, pl, pr, rhsq2, a1, a2, a3
+    real(kind=8) :: rhoim1, pim1, cim1, s0, rho1, rhou1, en1, p1, c1, s1, sfract
+    real(kind=8) :: rhoi, pi, ci, s3, rho2, rhou2, en2, p2, c2, s2, df
+    integer :: i, m, mw
+
     common /cparam/  gamma
 
     data efix /.true./    !# use entropy fix for transonic rarefactions
diff --git a/src/rp1_nonlinear_elasticity_fwave.f90 b/src/rp1_nonlinear_elasticity_fwave.f90
index 87351a7..22ced6d 100644
--- a/src/rp1_nonlinear_elasticity_fwave.f90
+++ b/src/rp1_nonlinear_elasticity_fwave.f90
@@ -1,5 +1,5 @@
 ! =====================================================
-subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq)
+subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq)
 ! =====================================================
 
 ! This version uses interleaved arrays.
@@ -44,18 +44,23 @@ subroutine rp1(maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,apdq)
 ! From the basic clawpack routines, this routine is called with ql = qr
 
 
-    implicit double precision (a-h,o-z)
+    implicit none
 
-    dimension auxl(maux,1-mbc:maxm+mbc)
-    dimension auxr(maux,1-mbc:maxm+mbc)
-    dimension fwave(meqn,mwaves,1-mbc:maxm+mbc)
-    dimension    s(mwaves,1-mbc:maxm+mbc)
-    dimension   ql(meqn,1-mbc:maxm+mbc)
-    dimension   qr(meqn,1-mbc:maxm+mbc)
-    dimension apdq(meqn,1-mbc:maxm+mbc)
-    dimension amdq(meqn,1-mbc:maxm+mbc)
+    integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxmx+mbc)
 
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: fwave(meqn, mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxmx+mbc)
 
+    integer :: i, m
+    real(kind=8) :: rhoi, rhoim, epsi, epsim, urhoi, urhoim
+    real(kind=8) :: bulki, bulkim, ci, cim, zi, zim, du, dsig, b1, b2
+    real(kind=8) :: sigma, sigmap
 !     # split the jump in q at each interface into waves
 
     do i = 2-mbc, mx+mbc
@@ -124,7 +129,10 @@ end subroutine rp1
 !     --------------------------------------------
     double precision function sigma(eps,i,coefA,coefB)
 !     --------------------------------------------
-    implicit double precision (a-h,o-z)
+    implicit none
+
+    real(kind=8) :: eps, coefA, coefB
+    integer :: i
 
 !     # stress-strain relation in ith cell
 
@@ -156,7 +164,10 @@ double precision function sigma(eps,i,coefA,coefB)
 !     --------------------------------------------
     double precision function sigmap(eps,i,coefA,coefB)
 !     --------------------------------------------
-    implicit double precision (a-h,o-z)
+    implicit none
+
+    real(kind=8) :: eps, coefA, coefB
+    integer :: i
 
 !     # derivative of stress-strain relation in ith cell
 
@@ -182,4 +193,4 @@ double precision function sigmap(eps,i,coefA,coefB)
 
 
     return
-    END function
+    END function
\ No newline at end of file
diff --git a/src/rp1_reactive_euler_with_efix.f90 b/src/rp1_reactive_euler_with_efix.f90
index a6c718f..721a345 100644
--- a/src/rp1_reactive_euler_with_efix.f90
+++ b/src/rp1_reactive_euler_with_efix.f90
@@ -26,20 +26,27 @@ subroutine rp1(maxmx,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apdq)
 !                                    and right state ql(:,i)
 ! From the basic clawpack routine step1, rp is called with ql = qr = q.
 
-
-    implicit double precision (a-h,o-z)
-    dimension   ql(meqn,1-mbc:maxmx+mbc)
-    dimension   qr(meqn,1-mbc:maxmx+mbc)
-    dimension    s(mwaves,1-mbc:maxmx+mbc)
-    dimension wave(meqn, mwaves,1-mbc:maxmx+mbc)
-    dimension amdq(meqn,1-mbc:maxmx+mbc)
-    dimension apdq(meqn,1-mbc:maxmx+mbc)
-
-!     # local storage
-!     ---------------
-    dimension delta(meqn)
-    dimension u(1-mbc:maxmx+mbc),enth(1-mbc:maxmx+mbc)
-    dimension a(1-mbc:maxmx+mbc)
+    implicit none
+
+    integer, intent(in) :: maxmx, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxmx+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxmx+mbc)
+
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxmx+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxmx+mbc)
+
+    real(kind=8) :: delta(meqn), gamma, gamma1, qheat
+    real(kind=8) :: u(1-mbc:maxmx+mbc),enth(1-mbc:maxmx+mbc)
+    real(kind=8) :: a(1-mbc:maxmx+mbc)
+    real(kind=8) :: a1, a2, a3, df, s0, s1, s2, s3, sfract
+    real(kind=8) :: rhoim1, pim1, cim1, rho1, rhou1, en1, qlam1, p1, c1
+    real(kind=8) :: rhoi, pi, ci, rho2, rhou2, en2, qlam2, p2, c2
+    real(kind=8) :: rhsqrtl, rhsqrtr, rhsq2, pl, pr
+    integer :: i, m, mw
     logical :: efix
     common /cparam/  gamma, qheat
 
diff --git a/src/rpn2_acoustics.f90 b/src/rpn2_acoustics.f90
index afd28b2..7875b53 100644
--- a/src/rpn2_acoustics.f90
+++ b/src/rpn2_acoustics.f90
@@ -31,21 +31,21 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apd
 !                                    and right state ql(i,:)
 ! From the basic clawpack routines, this routine is called with ql = qr
 
+    implicit none
 
-    implicit double precision (a-h,o-z)
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxm+mbc)
 
-    dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension apdq(meqn, 1-mbc:maxm+mbc)
-    dimension amdq(meqn, 1-mbc:maxm+mbc)
-    dimension auxl(maux, 1-mbc:maxm+mbc)
-    dimension auxr(maux, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxm+mbc)
 
-!     local arrays
-!     ------------
-    dimension delta(3)
+    real(kind=8) :: delta(3), a1, a2, rho, bulk, cc, zz
+    integer :: i, m, mu, mv
 
 !     # density, bulk modulus, and sound speed, and impedence of medium:
 !     # (should be set in setprob.f)
diff --git a/src/rpn2_acoustics_adjoint.f90 b/src/rpn2_acoustics_adjoint.f90
index c21ec65..fc45840 100644
--- a/src/rpn2_acoustics_adjoint.f90
+++ b/src/rpn2_acoustics_adjoint.f90
@@ -25,20 +25,24 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,ap
 
 !     # aux arrays not used in this solver.
 
-    implicit double precision (a-h,o-z)
-
-    dimension fwave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension apdq(meqn, 1-mbc:maxm+mbc)
-    dimension amdq(meqn, 1-mbc:maxm+mbc)
-    dimension auxl(maux, 1-mbc:maxm+mbc)
-    dimension auxr(maux, 1-mbc:maxm+mbc)
-
-!     local arrays
-!     ------------
-    dimension delta(3)
+    implicit none
+    !Input
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(in) :: ql, qr
+    real(kind=8), dimension(maux, 1-mbc:maxm+mbc), intent(in) :: auxl, auxr
+
+    !Output
+    real(kind=8), intent(out) :: fwave(meqn, mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(out) :: apdq, amdq
+
+    !Local
+    integer :: i, m, mu, mv
+    real(kind=8) :: rho, bulk, cc, zz
+    real(kind=8) :: rhoi, bulki, beta1, beta2
+    real(kind=8) :: delta(3)
+    
+
 
 !     # density, bulk modulus, and sound speed, and impedence of medium:
 !     # (should be set in setprob.f)
diff --git a/src/rpn2_advection.f90 b/src/rpn2_advection.f90
index 34da3ee..8e9f531 100644
--- a/src/rpn2_advection.f90
+++ b/src/rpn2_advection.f90
@@ -31,17 +31,22 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apd
 ! From the basic clawpack routines, this routine is called with ql = qr
 ! maux=0 and aux arrays are unused in this example.
 
+    implicit none
 
-    implicit double precision (a-h,o-z)
-    common /cparam/ u,v
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxm+mbc)
 
-    dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension  apdq(meqn, 1-mbc:maxm+mbc)
-    dimension  amdq(meqn, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxm+mbc)
 
+    real(kind=8) :: u, v
+    integer :: i
+    common /cparam/ u,v
 
 !$$$      do 30 i = 2-mbc, mx+mbc-1
     do 30 i = 2-mbc, mx+mbc
diff --git a/src/rpn2_burgers.f90 b/src/rpn2_burgers.f90
index f3deffc..167af11 100644
--- a/src/rpn2_burgers.f90
+++ b/src/rpn2_burgers.f90
@@ -25,23 +25,28 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apd
 !                                    and right state ql(i,:)
 ! From the basic clawpack routines, this routine is called with ql = qr
 
-    
-        implicit double precision (a-h,o-z)
-        
-        dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
-        dimension    s(mwaves, 1-mbc:maxm+mbc)
-        dimension   ql(meqn, 1-mbc:maxm+mbc)
-        dimension   qr(meqn, 1-mbc:maxm+mbc)
-        dimension  apdq(meqn, 1-mbc:maxm+mbc)
-        dimension  amdq(meqn, 1-mbc:maxm+mbc)
-        logical efix
+    implicit none
+
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), intent(in) :: ql(meqn,1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: qr(meqn,1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: auxl(maux,1-mbc:maxm+mbc)
+    real(kind=8), intent(in) :: auxr(maux,1-mbc:maxm+mbc)
+
+    real(kind=8), intent(out) :: s(mwaves,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: wave(meqn, mwaves,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: amdq(meqn,1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: apdq(meqn,1-mbc:maxm+mbc)
+
+    integer :: i
+    logical efix
         
 !     # x- and y- Riemann problems are identical, so it doesn't matter if
 !     # ixy=1 or 2.
 
         efix = .true.
 
-        do 10 i = 2-mbc, mx+mbc
+        do i = 2-mbc, mx+mbc
 !        # wave is jump in q, speed comes from R-H condition:
             wave(1,1,i) = ql(1,i) - qr(1,i-1)
             s(1,i) = 0.5d0*(qr(1,i-1) + ql(1,i))
@@ -58,7 +63,7 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apd
                     apdq(1,i) =   0.5d0 * ql(1,i)**2
                 end if
             end if
-        10 continue
+        end do
         
     return
     end subroutine rpn2
diff --git a/src/rpn2_euler_4wave.f90 b/src/rpn2_euler_4wave.f90
index a553536..d8fd8ee 100644
--- a/src/rpn2_euler_4wave.f90
+++ b/src/rpn2_euler_4wave.f90
@@ -34,27 +34,45 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apd
 ! From the basic clawpack routines, this routine is called with ql = qr
 !
 !
-      implicit double precision (a-h,o-z)
-!
-      dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
-      dimension    s(mwaves, 1-mbc:maxm+mbc)
-      dimension   ql(meqn,1-mbc:maxm+mbc)
-      dimension   qr(meqn,1-mbc:maxm+mbc)
-      dimension  apdq(meqn, 1-mbc:maxm+mbc)
-      dimension  amdq(meqn, 1-mbc:maxm+mbc)
-!
+    implicit none
+    !Input
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    double precision, intent(in) :: ql(meqn, 1-mbc:maxm+mbc)
+    double precision, intent(in) :: qr(meqn, 1-mbc:maxm+mbc)
+    double precision, intent(in) :: auxl(maux, 1-mbc:maxm+mbc)
+    double precision, intent(in) :: auxr(maux, 1-mbc:maxm+mbc)
+
+    double precision, intent(out) :: wave(meqn, mwaves, 1-mbc:maxm+mbc)
+    double precision, intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    double precision, intent(out) :: amdq(meqn, 1-mbc:maxm+mbc)
+    double precision, intent(out) :: apdq(meqn, 1-mbc:maxm+mbc)
+
+    !Local
+    integer :: i, m, mu, mv, mw
+    integer, parameter :: maxm2 = 1800
+    double precision :: delta(4)
+    logical :: efix
+    double precision :: gamma, gamma1
+    double precision :: u2v2(-6:maxm2+7), u(-6:maxm2+7), v(-6:maxm2+7)
+    double precision :: enth(-6:maxm2+7), a(-6:maxm2+7), g1a2(-6:maxm2+7)
+    double precision :: euv(-6:maxm2+7)
+    double precision :: rhsqrtl, rhsqrtr, pl, pr, rhsq2, a2, a3, a1, a4
+    double precision :: s0, s1, s2, s3, sfract
+    double precision :: rhoim1, pim1, cim1, rho1, rhou1, rhov1, en1, p1, c1
+    double precision :: rhoi, pi, ci, rho2, rhou2, rhov2, en2, p2, c2
+    double precision :: df
+
+
 !     local arrays -- common block comroe is passed to rpt2eu
 !     ------------
-      parameter (maxm2 = 1800)  !# assumes at most 200x200 grid with mbc=2
-      dimension delta(4)
-      logical efix
+      ! parameter (maxm2 = 1800)  !# assumes at most 200x200 grid with mbc=2
+      ! dimension delta(4)
+      ! logical efix
       common /cparam/  gamma
 !
       data efix /.true./    !# use entropy fix for transonic rarefactions
 
-      double precision u2v2(-6:maxm2+7), &
-            u(-6:maxm2+7),v(-6:maxm2+7),enth(-6:maxm2+7),a(-6:maxm2+7), &
-            g1a2(-6:maxm2+7),euv(-6:maxm2+7) 
+
 !
       if (mbc > 7 .OR. maxm2 < maxm) then
          write(6,*) 'need to increase maxm2 or 7 in rpn'
diff --git a/src/rpn2_psystem.f90 b/src/rpn2_psystem.f90
index 66783c2..e0c5fde 100644
--- a/src/rpn2_psystem.f90
+++ b/src/rpn2_psystem.f90
@@ -42,15 +42,43 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,ap
 !                                    and right state ql(i,:)
 ! From the basic clawpack routines, this routine is called with ql = qr
 
-    implicit double precision (a-h,o-z)
-    dimension fwave(meqn,mwaves,1-mbc:maxm+mbc)
-    dimension    s(mwaves,1-mbc:maxm+mbc)
-    dimension   ql(meqn,1-mbc:maxm+mbc)
-    dimension   qr(meqn,1-mbc:maxm+mbc)
-    dimension apdq(meqn,1-mbc:maxm+mbc)
-    dimension amdq(meqn,1-mbc:maxm+mbc)
-    dimension auxl(maux,1-mbc:maxm+mbc)
-    dimension auxr(maux,1-mbc:maxm+mbc)
+    implicit none
+
+    ! Input
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(in) :: ql, qr
+    real(kind=8), dimension(maux, 1-mbc:maxm+mbc), intent(in) :: auxl, auxr
+
+    !Output
+    real(kind=8), intent(out) :: fwave(meqn, mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(out) :: apdq, amdq
+
+    !Local
+    integer :: i
+    real(kind=8) :: epsi, urhoi, vrhoi, epsim, urhoim, vrhoim
+    real(kind=8) :: sigmai, sigmaim, sigmapi, sigmapim
+    real(kind=8) :: pim, pi, Eim, Ei
+    integer :: linearity_mati, linearity_matim
+    real(kind=8) :: r11, r13
+    real(kind=8) :: dF1, dF2, dF3
+    real(kind=8) :: beta1, beta3
+
+    interface
+        function sigma(eps,E,linearity_mat)
+            implicit none
+            integer, intent(in) :: linearity_mat
+            real(kind=8), intent(in) :: eps, E
+            real(kind=8) :: sigma
+        end function
+        function sigmap(eps,E,linearity_mat)
+            implicit none
+            integer, intent(in) :: linearity_mat
+            real(kind=8), intent(in) :: eps, E
+            real(kind=8) :: sigmap
+        end function
+    end interface
+
           
     do 10 i=2-mbc,mx+mbc
     ! material properties
@@ -90,11 +118,11 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,ap
             beta3=(-dF1+r11*dF2)/(r11-r13)
         ! compute f-waves
             fwave(1,1,i)=beta1*r11
-            fwave(2,1,i)=beta1*1
-            fwave(3,1,i)=beta1*0
+            fwave(2,1,i)=beta1*1.0d0
+            fwave(3,1,i)=beta1*0.0d0
             fwave(1,2,i)=beta3*r13
-            fwave(2,2,i)=beta3*1
-            fwave(3,2,i)=beta3*0
+            fwave(2,2,i)=beta3*1.0d0
+            fwave(3,2,i)=beta3*0.0d0
         else                   !y direction
         ! compute jump in flux
             dF1=-(vrhoi/pi-vrhoim/pim)
@@ -104,11 +132,11 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,ap
             beta3=(-dF1+r11*dF3)/(r11-r13)
         ! compute f-waves
             fwave(1,1,i)=beta1*r11
-            fwave(2,1,i)=beta1*0
-            fwave(3,1,i)=beta1*1
+            fwave(2,1,i)=beta1*0.0d0
+            fwave(3,1,i)=beta1*1.0d0
             fwave(1,2,i)=beta3*r13
-            fwave(2,2,i)=beta3*0
-            fwave(3,2,i)=beta3*1
+            fwave(2,2,i)=beta3*0.0d0
+            fwave(3,2,i)=beta3*1.0d0
         endif
     ! computation of the fluctuations
         amdq(1,i)=fwave(1,1,i)
@@ -124,33 +152,41 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,ap
     end subroutine rpn2
 
     double precision function sigma(eps,E,linearity_mat)
-!     Returns the flux sigma for a given
-!     eps, E and depending the linearity of the material
-    implicit double precision (a-h,o-z)
-    beta=5
-    select case (linearity_mat)
-    case (1)
-    sigma=E*eps
-    case (2)
-    sigma=dexp(E*eps)-1
-    case (3)
-    sigma=0.1*E*eps+beta*eps**3*E**3
-    end select
-    return
+    !Returns the flux sigma for a given
+    !eps, E and depending the linearity of the material
+        implicit none
+        !Input
+        integer, intent(in) :: linearity_mat
+        real(kind=8), intent(in) :: eps, E
+        real(kind=8) :: beta
+        beta=5.d0
+        select case (linearity_mat)
+        case (1)
+        sigma=E*eps
+        case (2)
+        sigma=dexp(E*eps)-1
+        case (3)
+        sigma=0.1*E*eps+beta*eps**3*E**3
+        end select
+        return
     END function
 
     double precision function sigmap(eps,E,linearity_mat)
-    implicit double precision (a-h,o-z)
-! Returns the derivative of sigma wrt eps for a given
-!     eps, E and depending the linearity of the material
-    beta=5
-    select case (linearity_mat)
-    case (1)
-    sigmap=E
-    case (2)
-    sigmap=E*dexp(E*eps)
-    case (3)
-    sigmap=0.1*E+3*beta*eps**2*E**3
-    end select
-    return
+    ! Returns the derivative of sigma wrt eps for a given
+    !     eps, E and depending the linearity of the material
+        implicit none
+        !Input
+        integer, intent(in) :: linearity_mat
+        real(kind=8), intent(in) :: eps, E
+        real(kind=8) :: beta
+        beta=5.d0
+        select case (linearity_mat)
+        case (1)
+        sigmap=E
+        case (2)
+        sigmap=E*dexp(E*eps)
+        case (3)
+        sigmap=0.1*E+3*beta*eps**2*E**3
+        end select
+        return
     END function
diff --git a/src/rpn2_shallow_roe_with_efix.f90 b/src/rpn2_shallow_roe_with_efix.f90
index 2f1c39f..dc6572f 100644
--- a/src/rpn2_shallow_roe_with_efix.f90
+++ b/src/rpn2_shallow_roe_with_efix.f90
@@ -36,29 +36,33 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apd
 ! This Riemann solver differs from the original clawpack Riemann solver
 ! for the interleaved indices
 
-    implicit double precision (a-h,o-z)
-
-    dimension   ql(meqn,           1-mbc:maxm+mbc)
-    dimension   qr(meqn,           1-mbc:maxm+mbc)
-    dimension    s(mwaves,         1-mbc:maxm+mbc)
-    dimension wave(meqn,   mwaves, 1-mbc:maxm+mbc)
-    dimension  apdq(meqn,          1-mbc:maxm+mbc)
-    dimension  amdq(meqn,          1-mbc:maxm+mbc)
+    implicit none
+    !Input
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(in) :: ql, qr
+    real(kind=8), dimension(maux, 1-mbc:maxm+mbc), intent(in) :: auxl, auxr
+
+    !Output
+    real(kind=8), intent(out) :: wave(meqn, mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(out) :: apdq, amdq
+
+    !Local
+    integer :: i, m, mu, mv, mw
+    integer, parameter :: maxm2 = 1800
+    real(kind=8) :: grav, hsqrtl, hsqrtr, hsq2, him1, s0, h1, hu1, s1, sfract
+    real(kind=8) :: hi, s03, h3, hu3, s3, df
+    real(kind=8) :: delta(3), a1, a2, a3
+    logical :: efix
 
 !   # Gravity constant set in the shallow1D.py or setprob.f90 file
     common /cparam/ grav
 
 !   # Roe averages quantities of each interface
-    parameter (maxm2 = 1800)
-    double precision u(-6:maxm2+7),v(-6:maxm2+7),a(-6:maxm2+7), &
+    double precision :: u(-6:maxm2+7),v(-6:maxm2+7),a(-6:maxm2+7), &
                      h(-6:maxm2+7)
 
 
-!   local arrays
-!   ------------
-    dimension delta(3)
-    logical :: efix
-
     data efix /.true./    !# Use entropy fix for transonic rarefactions
 
 !   # Set mu to point to  the component of the system that corresponds
diff --git a/src/rpn2_shallow_sphere.f90 b/src/rpn2_shallow_sphere.f90
index 10a430c..eb2945f 100644
--- a/src/rpn2_shallow_sphere.f90
+++ b/src/rpn2_shallow_sphere.f90
@@ -54,26 +54,34 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apd
 ! From the basic clawpack routines, this routine is called with ql = qr
 
 
-    implicit double precision (a-h,o-z)
-
-    dimension wave(meqn, mwaves,1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension  apdq(meqn, 1-mbc:maxm+mbc)
-    dimension  amdq(meqn, 1-mbc:maxm+mbc)
-    dimension auxl(maux, 1-mbc:maxm+mbc)
-    dimension auxr(maux, 1-mbc:maxm+mbc)
-
-!     local arrays
-!     ------------
-    dimension delta(3)
+    implicit none
+    !Input
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(in) :: ql, qr
+    real(kind=8), dimension(maux, 1-mbc:maxm+mbc), intent(in) :: auxl, auxr
+
+    !Output
+    real(kind=8), intent(out) :: wave(meqn, mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(out) :: apdq, amdq
+
+    !Local
+    integer :: i, mw, ioff
+    real(kind=8) :: enx, eny, enz, etx, ety, etz
+    real(kind=8) :: gamma, hunl, hunr, hutl, hutr, hsqr, hsql, hsq
+    real(kind=8) :: hl, hr, delta(3), a1, a2, a3, amn, apn
+    real(kind=8) :: df, erx, ery, erz, sfract
+    real(kind=8) :: h1, hu1, s1, h3, hu3, s3, s0, s03
+    real(kind=8) :: him1, h3m1
+    real(kind=8) :: g, sw !Not sure about sw, it is not used anywhere
+    real(kind=8) :: dtcom, dxcom, dycom, tcom
+    integer :: icom, jcom, m
+    real(kind=8) :: dy, hi
+    real(kind=8), dimension(1-mbc:maxm+mbc) :: u, v, a, h
+    integer, parameter :: maxm2 = 1800
     logical :: efix
 
-    parameter (maxm2 = 1800)
     common /sw/  g
-    dimension u(1-mbc:maxm+mbc),v(1-mbc:maxm+mbc),a(1-mbc:maxm+mbc), &
-    h(1-mbc:maxm+mbc)
     common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
 
     data efix /.true./    !# use entropy fix for transonic rarefactions
diff --git a/src/rpn2_vc_acoustics.f90 b/src/rpn2_vc_acoustics.f90
index 84611b2..99ba0ac 100644
--- a/src/rpn2_vc_acoustics.f90
+++ b/src/rpn2_vc_acoustics.f90
@@ -41,20 +41,24 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apd
 ! From the basic clawpack routines, this routine is called with ql = qr
 
 
-    implicit double precision (a-h,o-z)
-
-    dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension apdq(meqn, 1-mbc:maxm+mbc)
-    dimension amdq(meqn, 1-mbc:maxm+mbc)
-    dimension auxl(maux, 1-mbc:maxm+mbc)
-    dimension auxr(maux, 1-mbc:maxm+mbc)
-
-!     local arrays
-!     ------------
-    dimension delta(3)
+    implicit none
+
+    ! Input
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(in) :: ql, qr
+    real(kind=8), dimension(maux, 1-mbc:maxm+mbc), intent(in) :: auxl, auxr
+
+    !Output
+    real(kind=8), intent(out) :: wave(meqn, mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(out) :: apdq, amdq
+
+    !Local
+    integer :: i, m, mu, mv
+    real(kind=8) :: delta(3)
+    real(kind=8) :: zi, zim, a1, a2
+
+
 
 !     # set mu to point to  the component of the system that corresponds
 !     # to velocity in the direction of this slice, mv to the orthogonal
diff --git a/src/rpn2_vc_acoustics_adjoint.f90 b/src/rpn2_vc_acoustics_adjoint.f90
index 643f420..6547f25 100644
--- a/src/rpn2_vc_acoustics_adjoint.f90
+++ b/src/rpn2_vc_acoustics_adjoint.f90
@@ -39,20 +39,24 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,fwave,s,amdq,ap
 !     #                                    and right state ql(i,:)
 !     # From the basic clawpack routines, this routine is called with ql = qr
 
-    implicit double precision (a-h,o-z)
-
-    dimension fwave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension apdq(meqn, 1-mbc:maxm+mbc)
-    dimension amdq(meqn, 1-mbc:maxm+mbc)
-    dimension auxl(maux, 1-mbc:maxm+mbc)
-    dimension auxr(maux, 1-mbc:maxm+mbc)
+    implicit none
+    !Input
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(in) :: ql, qr
+    real(kind=8), dimension(maux, 1-mbc:maxm+mbc), intent(in) :: auxl, auxr
+
+    !Output
+    real(kind=8), intent(out) :: fwave(meqn, mwaves, 1-mbc:maxm+mbc) 
+    real(kind=8), intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(out) :: apdq, amdq
+
+    !Local
+    integer :: i, m, mu, mv
+    real(kind=8) :: zi, zim, ci, cim, rhoi, rhoim, bulki, bulkim, beta1, beta2
+    real(kind=8) :: delta(3)
 
 !     local arrays
 !     ------------
-    dimension delta(3)
 
 !     # set mu to point to  the component of the system that corresponds
 !     # to velocity in the direction of this slice, mv to the orthogonal
diff --git a/src/rpn2_vc_advection.f90 b/src/rpn2_vc_advection.f90
index 18eccae..28d6063 100644
--- a/src/rpn2_vc_advection.f90
+++ b/src/rpn2_vc_advection.f90
@@ -33,16 +33,20 @@ subroutine rpn2(ixy,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,auxl,auxr,wave,s,amdq,apd
 !                                    and right state ql(i,:)
 ! From the basic clawpack routines, this routine is called with ql = qr
 
-    implicit real*8(a-h,o-z)
-
-    dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(mwaves, 1-mbc:maxm+mbc)
-    dimension   ql(meqn, 1-mbc:maxm+mbc)
-    dimension   qr(meqn, 1-mbc:maxm+mbc)
-    dimension apdq(meqn,1-mbc:maxm+mbc)
-    dimension amdq(meqn,1-mbc:maxm+mbc)
-    dimension auxl(maux,1-mbc:maxm+mbc)
-    dimension auxr(maux,1-mbc:maxm+mbc)
+    implicit none
+
+    !Input
+    integer, intent(in) :: ixy, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(in) :: ql, qr
+    real(kind=8), dimension(maux, 1-mbc:maxm+mbc), intent(in) :: auxl, auxr
+
+    !Output
+    real(kind=8), intent(out) :: wave(meqn, mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), intent(out) :: s(mwaves, 1-mbc:maxm+mbc)
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(out) :: apdq, amdq
+
+    !Local
+    integer :: i
 
 
 !     # Set wave, speed, and flux differences:
diff --git a/src/rpt2_acoustics.f90 b/src/rpt2_acoustics.f90
index 50d161e..0ff603a 100644
--- a/src/rpt2_acoustics.f90
+++ b/src/rpt2_acoustics.f90
@@ -1,17 +1,26 @@
 !   =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 !   =====================================================
-    implicit double precision (a-h,o-z)
 
 !     # Riemann solver in the transverse direction for the acoustics equations.
 
 !     # Split asdq into down-going flux bmasdq and up-going flux bpasdq.
 
-    dimension     ql(meqn, 1-mbc:maxm+mbc)
-    dimension     qr(meqn, 1-mbc:maxm+mbc)
-    dimension   asdq(meqn, 1-mbc:maxm+mbc)
-    dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-    dimension bpasdq(meqn, 1-mbc:maxm+mbc)
+    implicit none
+
+    !Input
+    integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+    
+    !Local
+    integer :: i,mu,mv
+    double precision :: a1,a2
+    double precision :: rho,bulk,cc,zz
 
 !     # density, bulk modulus, and sound speed, and impedence of medium:
 !     # (should be set in setprob.f)
diff --git a/src/rpt2_acoustics_adjoint.f90 b/src/rpt2_acoustics_adjoint.f90
index d3a2751..937bbe2 100644
--- a/src/rpt2_acoustics_adjoint.f90
+++ b/src/rpt2_acoustics_adjoint.f90
@@ -1,17 +1,25 @@
 !   =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 !   =====================================================
-    implicit double precision (a-h,o-z)
 
 !     # Riemann solver in the transverse direction for the acoustics equations.
 
 !     # Split asdq into down-going flux bmasdq and up-going flux bpasdq.
-
-    dimension     ql(meqn, 1-mbc:maxm+mbc)
-    dimension     qr(meqn, 1-mbc:maxm+mbc)
-    dimension   asdq(meqn, 1-mbc:maxm+mbc)
-    dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-    dimension bpasdq(meqn, 1-mbc:maxm+mbc)
+    implicit none
+    !Input
+    integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+
+    !Local
+    integer :: i,mu,mv
+    double precision :: a1,a2
+    double precision :: rho,bulk,cc,zz
+    
 
 !     # density, bulk modulus, and sound speed, and impedence of medium:
 !     # (should be set in setprob.f)
diff --git a/src/rpt2_advection.f90 b/src/rpt2_advection.f90
index d255ba7..d435e21 100644
--- a/src/rpt2_advection.f90
+++ b/src/rpt2_advection.f90
@@ -1,20 +1,26 @@
 ! =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 ! =====================================================
-    implicit double precision (a-h,o-z)
-
 !     # Riemann solver in the transverse direction for the scalar equation
 
 !     # Split asdq (= A^* \Delta q, where * = + or -)
 !     # into down-going flux difference bmasdq (= B^- A^* \Delta q)
 !     #    and up-going flux difference bpasdq (= B^+ A^* \Delta q)
 
-
-    dimension     ql(meqn,1-mbc:maxm+mbc)
-    dimension     qr(meqn,1-mbc:maxm+mbc)
-    dimension   asdq(meqn,1-mbc:maxm+mbc)
-    dimension bmasdq(meqn,1-mbc:maxm+mbc)
-    dimension bpasdq(meqn,1-mbc:maxm+mbc)
+    implicit none
+    !Input
+    integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+    !Local
+    integer :: i
+    double precision :: stran
+    double precision :: u,v
+    double precision :: stranm,stranp
     common /cparam/ u,v
 
 !     # transverse wave speeds have been computed in rpn2
diff --git a/src/rpt2_burgers.f90 b/src/rpt2_burgers.f90
index 230ed4d..f3c658e 100644
--- a/src/rpt2_burgers.f90
+++ b/src/rpt2_burgers.f90
@@ -2,19 +2,24 @@
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 !   =====================================================
 
-        implicit double precision (a-h,o-z)
-
 !     # Riemann solver in the transverse direction for 2D Burgers' equation
 !
 !     # Split asdq into eigenvectors of Roe matrix B.
 !     # For the scalar equation, this simply amounts to computing the
 !     # transverse wave speed from the opposite Riemann problem.
+        implicit none
+        !Input
+        integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+        double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+        double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+        double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+        !Output
+        double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
 
-        dimension    ql(meqn, 1-mbc:maxm+mbc)
-        dimension    qr(meqn, 1-mbc:maxm+mbc)
-        dimension   asdq(meqn, 1-mbc:maxm+mbc)
-        dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-        dimension bpasdq(meqn, 1-mbc:maxm+mbc)
+        !Local
+        integer :: i
+        double precision :: sb
         
 !     # x- and y- Riemann problems are identical, so it doesn't matter if
 !     # ixy=1 or 2.
diff --git a/src/rpt2_dummy.f90 b/src/rpt2_dummy.f90
index 3afdbb5..234df68 100644
--- a/src/rpt2_dummy.f90
+++ b/src/rpt2_dummy.f90
@@ -1,8 +1,15 @@
 ! =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 ! =====================================================
-    implicit double precision (a-h,o-z)
-
+    implicit none
+    !Input
+    integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+    
 !     # Dummy transverse Riemann solver, for use in dimensionally-split algorithm.
 
     write(*,*) 'Error: Dummy transverse Riemann solver called!'
diff --git a/src/rpt2_euler.f90 b/src/rpt2_euler.f90
index 7bd3a60..6ea891e 100644
--- a/src/rpt2_euler.f90
+++ b/src/rpt2_euler.f90
@@ -1,7 +1,6 @@
 ! =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 ! =====================================================
-      implicit double precision (a-h,o-z)
 !
 !     # Riemann solver in the transverse direction for the Euler equations.
 !     # Split asdq (= A^* \Delta q, where * = + or -)
@@ -11,19 +10,26 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,b
 !     # Uses Roe averages and other quantities which were 
 !     # computed in rpn2eu and stored in the common block comroe.
 !
-      dimension     ql(meqn, 1-mbc:maxm+mbc)
-      dimension     qr(meqn, 1-mbc:maxm+mbc)
-      dimension   asdq(meqn, 1-mbc:maxm+mbc)
-      dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-      dimension bpasdq(meqn, 1-mbc:maxm+mbc)
+      implicit none
+      !Input
+      integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+      double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+      double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+      double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+      !Output
+      double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+
+      !Local
+      integer :: i,m,mu,mv,mw
+      double precision :: waveb(4,3),sb(3)
+      integer, parameter :: maxm2 = 1800 !# assumes at most 200x200 grid with mbc=2
+      double precision, dimension(-6:maxm2+7) :: u2v2,u,v,enth,a,g1a2,euv
+      double precision :: rhsqrtl,rhsqrtr,pl,pr,rhsq2
+      double precision :: gamma,gamma1,a1,a2,a3,a4
 !
       common /cparam/  gamma
-      dimension waveb(4,3),sb(3)
-      parameter (maxm2 = 1800)  !# assumes at most 200x200 grid with mbc=2
-!
-      double precision u2v2(-6:maxm2+7), &
-            u(-6:maxm2+7),v(-6:maxm2+7),enth(-6:maxm2+7),a(-6:maxm2+7), &
-            g1a2(-6:maxm2+7),euv(-6:maxm2+7) 
+
 
       if (mbc > 7 .OR. maxm2 < maxm) then
          write(6,*) 'need to increase maxm2 or 7 in rpn'
diff --git a/src/rpt2_psystem.f90 b/src/rpt2_psystem.f90
index 6b59f2f..4e0dd48 100644
--- a/src/rpt2_psystem.f90
+++ b/src/rpt2_psystem.f90
@@ -14,17 +14,33 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,b
 !     # into down-going flux difference bmasdq (= B^- A^* \Delta q)
 !     #    and up-going flux difference bpasdq (= B^+ A^* \Delta q)
 
-
-    implicit double precision (a-h,o-z)
-    dimension     ql(meqn,1-mbc:maxm+mbc)
-    dimension     qr(meqn,1-mbc:maxm+mbc)
-    dimension   asdq(meqn,1-mbc:maxm+mbc)
-    dimension bmasdq(meqn,1-mbc:maxm+mbc)
-    dimension bpasdq(meqn,1-mbc:maxm+mbc)
-    dimension   aux1(maux,1-mbc:maxm+mbc)
-    dimension   aux2(maux,1-mbc:maxm+mbc)
-    dimension   aux3(maux,1-mbc:maxm+mbc)
-    dimension s (2,1-mbc:maxm+mbc)
+    implicit none
+    !Input 
+    integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+
+    !Local
+    integer :: i,ix
+    double precision :: pjm,pj,pjp,Ejm,Ej,Ejp
+    double precision :: epsjm,epsj,epsjp
+    double precision :: sigmapjm,sigmapj,sigmapjp
+    double precision :: r11,r13, gamma1,gamma3
+    integer :: linearity_matjm,linearity_matj,linearity_matjp
+    double precision :: s(2,1-mbc:maxm+mbc)
+
+    interface
+        function sigmap(eps,E,linearity_mat)
+            implicit none
+            integer, intent(in) :: linearity_mat
+            real(kind=8), intent(in) :: eps, E
+            real(kind=8) :: sigmap
+        end function
+    end interface
 
     do 10 i = 2-mbc, mx+mbc
         if (imp == 1) then !left going fluctuation
diff --git a/src/rpt2_shallow_roe_with_efix.f90 b/src/rpt2_shallow_roe_with_efix.f90
index c13585d..311cc10 100644
--- a/src/rpt2_shallow_roe_with_efix.f90
+++ b/src/rpt2_shallow_roe_with_efix.f90
@@ -1,7 +1,6 @@
 ! =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 ! =====================================================
-    implicit double precision (a-h,o-z)
 
 !     # Riemann solver in the transverse direction for the shallow water
 !     equations .
@@ -9,21 +8,28 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,b
 !     # into down-going flux difference bmasdq (= B^- A^* \Delta q)
 !     #    and up-going flux difference bpasdq (= B^+ A^* \Delta q)
 
-    dimension     ql(meqn, 1-mbc:maxm+mbc)
-    dimension     qr(meqn, 1-mbc:maxm+mbc)
-    dimension   asdq(meqn, 1-mbc:maxm+mbc)
-    dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-    dimension bpasdq(meqn, 1-mbc:maxm+mbc)
+    implicit none
+    !Input
+    integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
 
-    dimension waveb(3,3),sb(3)
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+
+    !Local
+    integer :: i,m,mw,mu,mv
+    double precision :: hsqrtl,hsqrtr,hsq2,a1,a2,a3
+    integer, parameter :: maxm2 = 1800
+    !   # Roe averages quantities of each interface
+    double precision, dimension(-6:maxm2+7) :: u,v,a,h
+    double precision :: waveb(3,3),sb(3)
+    double precision :: grav
 
 !   # grav must be set elsewhere
     common /cparam/ grav
 
-!   # Roe averages quantities of each interface
-    parameter (maxm2 = 1800)
-    double precision u(-6:maxm2+7),v(-6:maxm2+7),a(-6:maxm2+7), &
-                     h(-6:maxm2+7)
 
     if (ixy == 1) then
         mu = 2
diff --git a/src/rpt2_shallow_sphere.f90 b/src/rpt2_shallow_sphere.f90
index 51d64ab..f606271 100644
--- a/src/rpt2_shallow_sphere.f90
+++ b/src/rpt2_shallow_sphere.f90
@@ -1,7 +1,6 @@
 ! =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 ! =====================================================
-    implicit double precision (a-h,o-z)
 
 !     # Riemann solver in the transverse direction for the shallow water
 !     # equations  on a quadrilateral grid.
@@ -10,27 +9,32 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,b
 !     # into down-going flux difference bmasdq (= B^- A^* \Delta q)
 !     #    and up-going flux difference bpasdq (= B^+ A^* \Delta q)
 
-    dimension     ql(meqn, 1-mbc:maxm+mbc)
-    dimension     qr(meqn, 1-mbc:maxm+mbc)
-    dimension   asdq(meqn, 1-mbc:maxm+mbc)
-    dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-    dimension bpasdq(meqn, 1-mbc:maxm+mbc)
-    dimension   aux1(maux, 1-mbc:maxm+mbc)
-    dimension   aux2(maux, 1-mbc:maxm+mbc)
-    dimension   aux3(maux, 1-mbc:maxm+mbc)
-
-!      parameter (maxm2 = 1002)  !# assumes at most 1000x1000 grid with mbc=2
-    dimension u(1-mbc:maxm+mbc),v(1-mbc:maxm+mbc),a(1-mbc:maxm+mbc), &
-    h(1-mbc:maxm+mbc)
-    dimension wave(meqn, mwaves, 1-mbc:maxm+mbc)
-    dimension    s(3, 1-mbc:maxm+mbc)
-    dimension enx(1-mbc:maxm+mbc), eny(1-mbc:maxm+mbc), &
-    enz(1-mbc:maxm+mbc)
-    dimension etx(1-mbc:maxm+mbc), ety(1-mbc:maxm+mbc), &
-    etz(1-mbc:maxm+mbc)
-    dimension gamma(1-mbc:maxm+mbc)
-         
-    dimension delta(4)
+    implicit none
+    !Input
+    integer, intent(in) :: ixy, imp, maxm, meqn, mwaves, maux, mbc, mx
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(in) :: ql, qr
+    real(kind=8), dimension(maux, 1-mbc:maxm+mbc), intent(in) :: aux1, aux2, aux3
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(in) :: asdq
+
+    !Output
+    real(kind=8), dimension(meqn, 1-mbc:maxm+mbc), intent(out) :: bmasdq, bpasdq
+
+    !Local
+    integer :: i, i1, ioff, ix1, ixm1
+    real(kind=8) :: dx
+    real(kind=8) :: a1, a2, a3
+    real(kind=8), dimension(1-mbc:maxm+mbc) :: h, u, v, a
+    real(kind=8) :: wave(meqn, mwaves, 1-mbc:maxm+mbc)
+    real(kind=8) :: s(3, 1-mbc:maxm+mbc)
+    real(kind=8), dimension(1-mbc:maxm+mbc) :: enx, eny, enz, etx, ety, etz
+    real(kind=8) :: gamma(1-mbc:maxm+mbc)
+    real(kind=8) :: delta(4)
+    real(kind=8) :: sw, g
+    real(kind=8) :: dtcom, dxcom, dycom, tcom
+    real(kind=8) :: erx, ery, erz, bn
+    integer :: icom, jcom, m, mw
+
+
     common /sw/  g
     common /comxyt/ dtcom,dxcom,dycom,tcom,icom,jcom
 
diff --git a/src/rpt2_vc_acoustics.f90 b/src/rpt2_vc_acoustics.f90
index e40e832..00e7177 100644
--- a/src/rpt2_vc_acoustics.f90
+++ b/src/rpt2_vc_acoustics.f90
@@ -1,7 +1,6 @@
 ! =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 ! =====================================================
-    implicit double precision (a-h,o-z)
 
 !     # Riemann solver in the transverse direction for the acoustics equations
 !     # with varying material properties rho and kappa
@@ -16,14 +15,22 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,b
 
 !     # imp=1  means  asdq=amdq,    imp=2 means asdq=apdq
 
-    dimension    ql(meqn, 1-mbc:maxm+mbc)
-    dimension    qr(meqn, 1-mbc:maxm+mbc)
-    dimension    asdq(meqn, 1-mbc:maxm+mbc)
-    dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-    dimension bpasdq(meqn, 1-mbc:maxm+mbc)
-    dimension   aux1(maux, 1-mbc:maxm+mbc)
-    dimension   aux2(maux, 1-mbc:maxm+mbc)
-    dimension   aux3(maux, 1-mbc:maxm+mbc)
+    implicit none
+
+    !Input
+    integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+
+    !Local
+    integer :: i,i1,mu,mv
+    double precision :: a1,a2
+    double precision :: cm,c,cp,zm,zz,zp
+
 
 
 
diff --git a/src/rpt2_vc_acoustics_adjoint.f90 b/src/rpt2_vc_acoustics_adjoint.f90
index 9bb0cec..4fece62 100644
--- a/src/rpt2_vc_acoustics_adjoint.f90
+++ b/src/rpt2_vc_acoustics_adjoint.f90
@@ -1,7 +1,6 @@
 !   =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 !   =====================================================
-    implicit double precision (a-h,o-z)
 
 !     # Riemann solver in the transverse direction for the acoustics equations.
 
@@ -13,14 +12,20 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,b
 
 !     # Split asdq into down-going flux bmasdq and up-going flux bpasdq.
 
-    dimension     ql(meqn, 1-mbc:maxm+mbc)
-    dimension     qr(meqn, 1-mbc:maxm+mbc)
-    dimension   asdq(meqn, 1-mbc:maxm+mbc)
-    dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-    dimension bpasdq(meqn, 1-mbc:maxm+mbc)
-    dimension   aux1(maux, 1-mbc:maxm+mbc)
-    dimension   aux2(maux, 1-mbc:maxm+mbc)
-    dimension   aux3(maux, 1-mbc:maxm+mbc)
+    implicit none
+    !Input
+    integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+
+    !Local
+    integer :: i,i1,mu,mv
+    double precision :: a1,a2
+    double precision :: cm,c,cp,zm,zz,zp
 
 
     if (ixy == 1) then
diff --git a/src/rpt2_vc_advection.f90 b/src/rpt2_vc_advection.f90
index 7b47996..5062c85 100644
--- a/src/rpt2_vc_advection.f90
+++ b/src/rpt2_vc_advection.f90
@@ -1,18 +1,21 @@
 ! =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 ! =====================================================
-    implicit double precision(a-h,o-z)
 
 !     # Riemann solver in the transverse direction for the advection equation.
 
-    dimension     ql(meqn, 1-mbc:maxm+mbc)
-    dimension     qr(meqn, 1-mbc:maxm+mbc)
-    dimension   asdq(meqn, 1-mbc:maxm+mbc)
-    dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-    dimension bpasdq(meqn, 1-mbc:maxm+mbc)
-    dimension   aux1(maux, 1-mbc:maxm+mbc)
-    dimension   aux2(maux, 1-mbc:maxm+mbc)
-    dimension   aux3(maux, 1-mbc:maxm+mbc)
+    implicit none
+    !Input
+    integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+
+    !Local
+    integer :: i,i1,kv
 
     kv = 3-ixy  !#  = 1 if ixy=2  or  = 2 if ixy=1
     do 10 i=2-mbc,mx+mbc
diff --git a/src/rpt2_vc_elasticity.f90 b/src/rpt2_vc_elasticity.f90
index 6dcf48d..e9f4a76 100644
--- a/src/rpt2_vc_elasticity.f90
+++ b/src/rpt2_vc_elasticity.f90
@@ -1,7 +1,6 @@
 ! =====================================================
 subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,bmasdq,bpasdq)
 ! =====================================================
-      implicit double precision (a-h,o-z)
 !
 !     # Riemann solver in the transverse direction for the elastic equations
 !     # with varying material properties 
@@ -32,14 +31,24 @@ subroutine rpt2(ixy,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3,asdq,b
 !
 !     # imp=1  means  asdq=amdq,    imp=2 means asdq=apdq
 !
-      dimension     ql(meqn, 1-mbc:maxm+mbc)
-      dimension     qr(meqn, 1-mbc:maxm+mbc)
-      dimension   asdq(meqn, 1-mbc:maxm+mbc)
-      dimension bmasdq(meqn, 1-mbc:maxm+mbc)
-      dimension bpasdq(meqn, 1-mbc:maxm+mbc)
-      dimension   aux1(maux, 1-mbc:maxm+mbc)
-      dimension   aux2(maux, 1-mbc:maxm+mbc)
-      dimension   aux3(maux, 1-mbc:maxm+mbc)
+      implicit none
+      !Input
+      integer, intent(in)  :: ixy,imp,maxm,meqn,mwaves,maux,mbc,mx
+      double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+      double precision, dimension(maux,1-mbc:maxm+mbc), intent(in) :: aux1,aux2,aux3
+      double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+      !Output
+      double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+
+      !Local
+      integer :: i,i1,ksig11,ksig22,ku,kv
+      double precision :: a1,a2,a3,a4
+      double precision :: alamm,alam,alamp,amum,amu,amup,bulkm,bulk,bulkp
+      double precision :: cpm,cp,cpp,csm,cs,csp
+      double precision :: dsig11,dsig22,dsig12,du,dv
+      double precision :: det
+
 !
 !
 !
diff --git a/src/rpt3_vc_acoustics.f90 b/src/rpt3_vc_acoustics.f90
index 2623543..572c131 100644
--- a/src/rpt3_vc_acoustics.f90
+++ b/src/rpt3_vc_acoustics.f90
@@ -57,15 +57,21 @@ subroutine rpt3(ixyz,icoor,imp,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux2,aux3
 !     #                z-like direction
 
 
-    implicit real*8(a-h,o-z)
-    dimension     ql(meqn,1-mbc:maxm+mbc)
-    dimension     qr(meqn,1-mbc:maxm+mbc)
-    dimension   asdq(meqn,1-mbc:maxm+mbc)
-    dimension bmasdq(meqn,1-mbc:maxm+mbc)
-    dimension bpasdq(meqn,1-mbc:maxm+mbc)
-    dimension   aux1(maux,1-mbc:maxm+mbc,3)
-    dimension   aux2(maux,1-mbc:maxm+mbc,3)
-    dimension   aux3(maux,1-mbc:maxm+mbc,3)
+    ! implicit real*8(a-h,o-z)
+    implicit none
+    !Input
+    integer, intent(in)  :: ixyz,icoor,imp,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc,3), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: asdq
+
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: bmasdq,bpasdq
+
+    !Local
+    integer :: i,i1,iuvw
+    double precision :: zm,zz,zp,cm,c,cp,a1,a2
+
 
 
 !     # set iuvw = 2,3,4, depending on which component of q represents
diff --git a/src/rptt3_vc_acoustics.f90 b/src/rptt3_vc_acoustics.f90
index e70dce2..863c7ca 100644
--- a/src/rptt3_vc_acoustics.f90
+++ b/src/rptt3_vc_acoustics.f90
@@ -65,15 +65,19 @@ subroutine rptt3(ixyz,icoor,imp,impt,maxm,meqn,mwaves,maux,mbc,mx,ql,qr,aux1,aux
 !     #                z-like direction
 
 
-    implicit real*8(a-h,o-z)
-    dimension       ql(meqn,1-mbc:maxm+mbc)
-    dimension       qr(meqn,1-mbc:maxm+mbc)
-    dimension   bsasdq(meqn,1-mbc:maxm+mbc)
-    dimension cmbsasdq(meqn,1-mbc:maxm+mbc)
-    dimension cpbsasdq(meqn,1-mbc:maxm+mbc)
-    dimension     aux1(maux,1-mbc:maxm+mbc,3)
-    dimension     aux2(maux,1-mbc:maxm+mbc,3)
-    dimension     aux3(maux,1-mbc:maxm+mbc,3)
+    implicit none
+    !Input
+    integer, intent(in)  :: ixyz,icoor,imp,impt,maxm,meqn,mwaves,maux,mbc,mx
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: ql,qr
+    double precision, dimension(maux,1-mbc:maxm+mbc,3), intent(in) :: aux1,aux2,aux3
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(in) :: bsasdq
+    
+    !Output
+    double precision, dimension(meqn,1-mbc:maxm+mbc), intent(out) :: cmbsasdq,cpbsasdq
+    
+    !Local
+    integer :: i,i1,iuvw
+    double precision :: zm,zz,zp,cm,c,cp,a1,a2
 
 
 !     # set iuvw = 2,3,4, depending on which component of q represents