anchor.diff:  
       dtdx = dt/U%grid%dx
   
  -    ! We need the left and right eigenvectors and the eigenvalues for
  -    ! the system
  -    !
  -    ! The eigenvalues are u - c, u, u + c
  -    !
  -    ! The matrix A in the primitive variable system is
  -    !
  -    !       / u   r   0   \  
  -    !   A = | 0   u   1/r |  
  -    !       \ 0  rc^2 u   /
  -    !
  -    ! With the rows corresponding to rho, u, and p
  -    !
  -    ! The right eigenvectors are
  -    !
  -    !       /  1  \        / 1 \        /  1  \
  -    ! r1 =  |-c/r |   r2 = | 0 |   r3 = | c/r |
  -    !       \ c^2 /        \ 0 /        \ c^2 /
  -    !
  -    !
  -    ! The left eigenvectors are    
  -    !
  -    !   l1 =     ( 0,  -r/(2c),  1/(2c^2) )
  -    !   l2 =     ( 1,     0,     -1/c^2,  )
  -    !   l3 =     ( 0,   r/(2c),  1/(2c^2) )
  -    !
  -    ! The fluxes are going to be defined on the left edge of the
  -    ! computational zone.
  -    !
  -    !          |             |             |             | 
  -    !          |             |             |             |  
  -    !         -+------+------+------+------+------+------+--   
  -    !          |     i-1     |      i      |     i+1     |   
  -    !                       * *           * 
  -    !                   q_l,i  q_r,i   q_l,i+1 
  -    !           
  -    ! q_l,i+1 are computed using the information in zone i
  -    !
  -    !
       ! The basic idea here is that we do a characteristic
       ! decomposition.  The jump in primitive variables (Q) can be
       ! transformed to a jump in characteristic variables using the left
  
positive.code:     type(gridedgevar_t) :: Q_l, Q_r
      type(gridvar_t) :: ldelta, temp
  
      real (kind=dp_t) :: rvec(nprim,nprim), lvec(nprim,nprim), eval(nprim)
      real (kind=dp_t) :: dQ(nprim)
  
      real (kind=dp_t) :: r, ux, p, cs
      real (kind=dp_t) :: ldr, ldu, ldp
      real (kind=dp_t) :: r_xm, r_xp, u_xm, u_xp, p_xm, p_xp
      real (kind=dp_t) :: beta_xm(nprim), beta_xp(nprim)
      real (kind=dp_t) :: sum_xm, sum_xp
      real (kind=dp_t) :: sum
  
  
negative.code: module interface_states_plm_module
  
    use datatypes_module
    use grid_module
    use variables_module
    use params_module
    use eos_module
    use eigen_module
  
    implicit none
  
  