anchor.diff:      ! 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
  +       Q_l%data(i,iqxvel) = Q_l%data(i,iqxvel) + HALF*dt*grav
  +       Q_r%data(i,iqxvel) = Q_r%data(i,iqxvel) + HALF*dt*grav
       enddo
   
       ! special fixes at the boundary -- gravity must be reflected
  
positive.code:           ! 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
  
            ! 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
  
         enddo
  
         ! finally, sum up all the jumps 
         
         ! density
         sum_xm = 0.0_dp_t
         sum_xp = 0.0_dp_t
         do n = 1, nwaves
            sum_xm = sum_xm + beta_xm(n)*rvec(n,1)
            sum_xp = sum_xp + beta_xp(n)*rvec(n,1)
  
negative.code: 
  
         ! pressure
         sum_xm = ZERO
         sum_xp = ZERO
         do n = 1, nwaves
            sum_xm = sum_xm + beta_xm(n)*rvec(n,3)
            sum_xp = sum_xp + beta_xp(n)*rvec(n,3)
  