anchor.diff:       double precision :: pres_zone, temp_zone
        double precision :: dist_P, dist_S
   
  +     type (eos_t) :: eos_state
  +
        integer :: i,j,k,ii,jj,kk,n
   
        !$OMP PARALLEL DO PRIVATE(i, j, k, xl, yl, zl, xx, yy, zz) &
  
positive.code:               temp_zone = temp_zone/(nsub*nsub*nsub)
                state(i,j,k,UFS:UFS-1+nspec) = state(i,j,k,UFS:UFS-1+nspec)/(nsub*nsub*nsub)
  
                ! thermodynamics
                if (interp_temp) then
                   state(i,j,k,UTEMP) = temp_zone
                   call eos_given_RTX(state(i,j,k,UEINT),pres_zone,state(i,j,k,URHO),state(i,j,k,UTEMP), &
                                        state(i,j,k,UFS:UFS-1+nspec))
                else
  
                   call eos_e_given_RPX(state(i,j,k,UEINT),state(i,j,k,UTEMP),state(i,j,k,URHO), &
                                        pres_zone,state(i,j,k,UFS:UFS-1+nspec))
  
                endif
  
                state(i,j,k,UEDEN) = state(i,j,k,URHO) * state(i,j,k,UEINT)
                state(i,j,k,UEINT) = state(i,j,k,URHO) * state(i,j,k,UEINT)
  
negative.code:      double precision :: a
       double precision :: length
  
       type (eos_t) :: eos_state
  
       ! Temporary storage variables in case we need to switch the primary and secondary.
  
       double precision, allocatable :: temp_model_state(:,:), temp_model_r(:)
  