!:EntropyFix:!
            ! Compute flux differences amdq and apdq.
            ! First compute amdq as sum of s*wave for left going waves.
            ! Incorporate entropy fix by adding a modified fraction of wave
            ! if s should change sign.
            DO j=2-mbc,mx+mbc
              ! 1-wave. Check for fully supersonic case
              s0=u(j-1)-c(j-1)
              IF ( s0>=zero .AND. speed(j,1) > zero ) THEN
                ! Everything is right-going
                Amdq(j,1:NrVars) = zero
                CYCLE
              END IF
              ! u-c to right of 1-wave
              rho1 = q1D(j-1,1) + wave(j,1,1)
              rhou1 = q1D(j-1,mu) + wave(j,mu,1)
              rhov1 = q1D(j-1,mv) + wave(j,mv,1)
              en1 = q1D(j-1,4) + wave(j,4,1)
              p1 = gamma1*(en1  - 0.5*(rhou1**2+rhov1**2)/rho1)
              IF ((p1<zero) .OR. (rho1<zero)) THEN
                PRINT *,'Negative pressure/density between 1-wave and 2,3-waves'
                CYCLE
                STOP
              END IF
              c1 = SQRT(gamma*p1/rho1)
              s1 = rhou1/rho1 - c1
              IF ( s0<zero .AND. s1>zero ) THEN
                ! Transonic rarefaction in the 1-wave
                sfract = s0 * (s1-speed(j,1)) / (s1-s0)
              ELSE IF (speed(j,1) < zero) THEN
                ! Left-going 1-wave
                sfract=speed(j,1)
              ELSE
                ! Right-going 1-wave
                sfract=zero  ! Never should reach this instruction
              END IF
              Amdq(j,1:NrVars) = sfract*wave(j,1:NrVars,1)

              ! Check 2,3-wave (contact+shear discontinuities)
              IF (speed(j,2) >= zero) CYCLE ! 2- and 3- wave are right-going
              Amdq(j,1:NrVars) = Amdq(j,1:NrVars) + &
                                 uR(j)*(wave(j,1:NrVars,2) + wave(j,1:NrVars,3))
              ! Check 4-wave. Compute u+c to left of 4-wave
              rho2 = q1D(j,1) - wave(j,1,4)
              rhou2 = q1D(j,mu) - wave(j,mu,4)
              rhov2 = q1D(j,mv) - wave(j,mv,4)
              en2 = q1D(j,4) - wave(j,4,4)
              p2 = gamma1*(en2 - 0.5d0*(rhou2**2+rhov2**2)/rho2)
              IF ((p2<=zero) .OR. (rho2<=zero)) THEN
                PRINT *,'Negative pressure/density between 2,3-waves and 4-wave'
                CYCLE
                STOP
              END IF
              c2 = SQRT(gamma*p2/rho2)
              s2 = rhou2/rho2 + c2; s3=u(j)+c(j)
              IF ( s2 < zero .AND. s3 > zero ) THEN
                ! Transonic rarefaction in the 4-wave
                sfract = s2 * (s3-speed(j,4)) / (s3-s2)
              ELSE IF (speed(j,4) < zero) THEN
                ! Left-going 4-wave
                sfract = speed(j,4)
              ELSE
                ! Right-going 4-wave
                CYCLE
              END IF
              Amdq(j,1:NrVars) = Amdq(j,1:NrVars) + sfract*wave(j,1:NrVars,4)
            END DO