anchor.diff:            ! here the sign() function makes sure we only add the right-moving
             ! waves
             beta_xp(m) = 0.25_dp_t*dtdx*(eval(3) - eval(m))* &
  -               (sign(1.0_dp_t, eval(m)) + 1.0_dp_t)*sum
  +               (sign(ONE, eval(m)) + ONE)*sum
   
             ! here the sign() function makes sure we only add the left-moving
             ! waves
             beta_xm(m) = 0.25_dp_t*dtdx*(eval(1) - eval(m))* &
  -               (1.0_dp_t - sign(1.0_dp_t, eval(m)))*sum
  +               (ONE - sign(ONE, eval(m)))*sum
   
          enddo
   
          ! finally, sum up all the jumps 
          
          ! density
  -       sum_xm = 0.0_dp_t
  -       sum_xp = 0.0_dp_t
  +       sum_xm = ZERO
  +       sum_xp = ZERO
          do n = 1, nwaves
             sum_xm = sum_xm + beta_xm(n)*rvec(n,1)
             sum_xp = sum_xp + beta_xp(n)*rvec(n,1)
  
positive.code:     ! apply the source terms
      !-------------------------------------------------------------------------
      do i = U%grid%lo, U%grid%hi+1
         Q_l%data(i,iqxvel) = Q_l%data(i,iqxvel) + 0.5_dp_t*dt*grav
         Q_r%data(i,iqxvel) = Q_r%data(i,iqxvel) + 0.5_dp_t*dt*grav
      enddo
  
      ! special fixes at the boundary -- gravity must be reflected
  
negative.code:        ! These expressions are the V_{L,R} in Colella (1990) at the 
         ! bottom of page 191.  They are also in Saltzman (1994) as
         ! V_ref on page 161.
         r_xp = r + HALF*(ONE - dtdx*max(eval(3), ZERO))*ldr
         r_xm = r - HALF*(ONE + dtdx*min(eval(1), ZERO))*ldr
  
         u_xp = ux + HALF*(ONE - dtdx*max(eval(3), ZERO))*ldu
         u_xm = ux - HALF*(ONE + dtdx*min(eval(1), ZERO))*ldu
  
         p_xp = p + HALF*(ONE - dtdx*max(eval(3), ZERO))*ldp
         p_xm = p - HALF*(ONE + dtdx*min(eval(1), ZERO))*ldp
  
         !                                                   ^
         ! Now compute the interface states.   These are the V expressions
  