anchor.diff:  
   
   
  -    ! loop over interfaces and fill the states 
  -    do i = U%grid%lo, U%grid%hi+1
  -       U_l%data(i,:) = U%data(i-1,:)
  -       U_r%data(i,:) = U%data(i,:)
  -    enddo
   
     end subroutine make_interface_states_plm
   
  
positive.code:     !                   central difference term     
    
  
      !-------------------------------------------------------------------------
      ! convert to primitve variables
      !-------------------------------------------------------------------------
  
  
  
  
      !-------------------------------------------------------------------------
  
negative.code:     ! compute the monotonized central differences
      !-------------------------------------------------------------------------
  
      ! 4th order MC limiting.  See Colella (1985) Eq. 2.5 and 2.6 or 
      ! Saltzman 1994, page 156
      
      call build(temp, U%grid, 1)
      call build(ldelta, U%grid, nprim)
  
      do n = 1, nprim
  
         ! first do the normal MC limiting 
         do i = Q%grid%lo-3, Q%grid%hi+3
  
            test = (Q%data(i+1,n) - Q%data(i,n))*(Q%data(i,n) - Q%data(i-1,n))
  
            if (test > 0.0_dp_t) then
               temp%data(i,1) = min(0.5_dp_t*abs(Q%data(i+1,n) - Q%data(i-1,n)), &
                                    min(2.0_dp_t*abs(Q%data(i+1,n) - Q%data(i,n)), &
                                        2.0_dp_t*abs(Q%data(i,n) - Q%data(i-1,n)))) * &
                                sign(1.0_dp_t,Q%data(i+1,n) - Q%data(i-1,n))
  
            else
               temp%data(i,1) = 0.0_dp_t
            endif
  
         enddo
      
  
         ! now do the fourth order part
         do i = Q%grid%lo-2, Q%grid%hi+2
            
            test = (Q%data(i+1,n) - Q%data(i,n))*(Q%data(i,n) - Q%data(i-1,n))
  
            if (test > 0.0_dp_t) then
               ldelta%data(i,n) = &
                    min((2.0_dp_t/3.0_dp_t)*abs(Q%data(i+1,n) - Q%data(i-1,n) - &
                        0.25_dp_t*(temp%data(i+1,1) + temp%data(i-1,1))), &
                    min(2.0*abs(Q%data(i+1,n) - Q%data(i,n)), &
                        2.0*abs(Q%data(i,n) - Q%data(i-1,n)))) * &
                    sign(1.0_dp_t, Q%data(i+1,n) - Q%data(i-1,n))
            else
               ldelta%data(i,n) = 0.0_dp_t
            endif
  
         enddo
  
      enddo
  
  
      !-------------------------------------------------------------------------
      ! compute left and right primitive variable states
      !-------------------------------------------------------------------------
      call build(Q_l, U%grid, nprim)
      call build(Q_r, U%grid, nprim)
  
      ! We need the left and right eigenvectors and the eigenvalues for
      ! the system
      !
      ! The eigenvalues are u - c, u, u + c
      !
      ! The Jacobian matrix for the primitive variable formulation of the
      ! Euler equations 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     |   
      !                       * *           * 
      !                   V_l,i  V_r,i   V_l,i+1 
      !           
      ! V_l,i+1 are computed using the information in zone i,j.
      !
  
  
  
      ! clean-up
      call destroy(ldelta)
  
  
      !-------------------------------------------------------------------------
      ! transform the states into conserved variables
      !-------------------------------------------------------------------------
      do i = U%grid%lo, U%grid%hi+1
  
         ! density
         U_l%data(i,iudens) = Q_l%data(i,iqdens) 
         U_r%data(i,iudens) = Q_r%data(i,iqdens) 
  
         ! momentum
         U_l%data(i,iumomx) = Q_l%data(i,iqdens)*Q_l%data(i,iqxvel)
         U_r%data(i,iumomx) = Q_r%data(i,iqdens)*Q_r%data(i,iqxvel)
  
         ! 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
      
  
      ! clean-up
      call destroy(Q_l)
      call destroy(Q_r)
  
  
      !-------------------------------------------------------------------------
  