diff --git a/src/simulation/m_data_output.fpp b/src/simulation/m_data_output.fpp
index cb41a187c..4f48958f6 100644
--- a/src/simulation/m_data_output.fpp
+++ b/src/simulation/m_data_output.fpp
@@ -1049,6 +1049,8 @@ contains
         real(kind(0d0)) :: E_e
         real(kind(0d0)), dimension(6) :: tau_e
         real(kind(0d0)) :: G
+        real(kind(0d0)) :: dyn_p
+
 
         integer :: i, j, k, l, s, q !< Generic loop iterator
 
@@ -1128,22 +1130,20 @@ contains
                         vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k, l)/rho
                     end do
 
+                    dyn_p = 0.5d0*rho*dot_product(vel, vel)
+
                     if (hypoelasticity) then
                         call s_compute_pressure( &
                             q_cons_vf(1)%sf(j - 2, k, l), &
                             q_cons_vf(alf_idx)%sf(j - 2, k, l), &
-                            0.5d0*(q_cons_vf(2)%sf(j - 2, k, l)**2.d0)/ &
-                            q_cons_vf(1)%sf(j - 2, k, l), &
-                            pi_inf, gamma, rho, qv, pres, &
+                            dyn_p, pi_inf, gamma, rho, qv, pres, &
                             q_cons_vf(stress_idx%beg)%sf(j - 2, k, l), &
                             q_cons_vf(mom_idx%beg)%sf(j - 2, k, l), G)
                     else
                         call s_compute_pressure( &
                             q_cons_vf(1)%sf(j - 2, k, l), &
                             q_cons_vf(alf_idx)%sf(j - 2, k, l), &
-                            0.5d0*(q_cons_vf(2)%sf(j - 2, k, l)**2.d0)/ &
-                            q_cons_vf(1)%sf(j - 2, k, l), &
-                            pi_inf, gamma, rho, qv, pres)
+                            dyn_p, pi_inf, gamma, rho, qv, pres)
                     end if
 
                     if (model_eqns == 4) then
@@ -1230,19 +1230,19 @@ contains
                             vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l)/rho
                         end do
 
+                        dyn_p = 0.5d0*rho*dot_product(vel, vel)
+
                         if (hypoelasticity) then
                             call s_compute_pressure( &
                                 q_cons_vf(1)%sf(j - 2, k - 2, l), &
                                 q_cons_vf(alf_idx)%sf(j - 2, k - 2, l), &
-                                0.5d0*(q_cons_vf(2)%sf(j - 2, k - 2, l)**2.d0)/ &
-                                q_cons_vf(1)%sf(j - 2, k - 2, l), &
-                                pi_inf, gamma, rho, qv, pres, &
+                                dyn_p, pi_inf, gamma, rho, qv, pres, &
                                 q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l), &
                                 q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l), G)
                         else
                             call s_compute_pressure(q_cons_vf(E_idx)%sf(j - 2, k - 2, l), &
                                                     q_cons_vf(alf_idx)%sf(j - 2, k - 2, l), &
-                                                    0.5d0*rho*dot_product(vel, vel), pi_inf, gamma, rho, qv, pres)
+                                                    dyn_p, pi_inf, gamma, rho, qv, pres)
                         end if
 
                         if (model_eqns == 4) then
@@ -1313,19 +1313,19 @@ contains
                                 vel(s) = q_cons_vf(cont_idx%end + s)%sf(j - 2, k - 2, l - 2)/rho
                             end do
 
+                            dyn_p = 0.5d0*rho*dot_product(vel, vel)
+
                             if (hypoelasticity) then
                                 call s_compute_pressure( &
                                     q_cons_vf(1)%sf(j - 2, k - 2, l - 2), &
                                     q_cons_vf(alf_idx)%sf(j - 2, k - 2, l - 2), &
-                                    0.5d0*(q_cons_vf(2)%sf(j - 2, k - 2, l - 2)**2.d0)/ &
-                                    q_cons_vf(1)%sf(j - 2, k - 2, l - 2), &
-                                    pi_inf, gamma, rho, qv, pres, &
+                                    dyn_p, pi_inf, gamma, rho, qv, pres, &
                                     q_cons_vf(stress_idx%beg)%sf(j - 2, k - 2, l - 2), &
                                     q_cons_vf(mom_idx%beg)%sf(j - 2, k - 2, l - 2), G)
                             else
                                 call s_compute_pressure(q_cons_vf(E_idx)%sf(j - 2, k - 2, l - 2), &
                                                         q_cons_vf(alf_idx)%sf(j - 2, k - 2, l - 2), &
-                                                        0.5d0*rho*dot_product(vel, vel), pi_inf, gamma, rho, qv, pres)
+                                                        dyn_p, pi_inf, gamma, rho, qv, pres)
                             end if
 
                             ! Compute mixture sound speed