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:        ! total energy
         call eos(eos_input_p, Q_l%data(i,iqpres), e, Q_l%data(i,iqdens))
         U_l%data(i,iuener) = Q_l%data(i,iqdens)*e + &
              0.5_dp_t*Q_l%data(i,iqdens)*Q_l%data(i,iqxvel)**2
  
         call eos(eos_input_p, Q_r%data(i,iqpres), e, Q_r%data(i,iqdens))
         U_r%data(i,iuener) = Q_r%data(i,iqdens)*e + &
              0.5_dp_t*Q_r%data(i,iqdens)*Q_r%data(i,iqxvel)**2
  
      enddo
      
  
negative.code: 
  
         ! velocity
         sum_xm = ZERO
         sum_xp = ZERO
         do n = 1, nwaves
            sum_xm = sum_xm + beta_xm(n)*rvec(n,2)
            sum_xp = sum_xp + beta_xp(n)*rvec(n,2)
  