anchor.diff:             enddo
           enddo
        enddo
  +     !$OMP END PARALLEL DO
   
        ! Initial velocities = 0
   
  
positive.code:      if (adv_h3 > domhi(3)) then
          do k = domhi(3)+1, adv_h3
             do j = adv_l2, adv_h2
                do i = adv_l1, adv_h1
                   adv(i,j,k,URHO) = dens_ambient
                   adv(i,j,k,UTEMP) = temp_ambient
                   adv(i,j,k,UFS:UFS-1+nspec) = dens_ambient*xn_ambient(:)
  
                   vx = adv(i,j,domhi(3),UMX)/adv(i,j,domhi(3),URHO)
                   vy = adv(i,j,domhi(3),UMY)/adv(i,j,domhi(3),URHO)
                   vz = adv(i,j,domhi(3),UMZ)/adv(i,j,domhi(3),URHO)
                   adv(i,j,k,UMX) = dens_ambient*vx
                   adv(i,j,k,UMY) = dens_ambient*vy
                   adv(i,j,k,UMZ) = dens_ambient*vz
  
negative.code:                                       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))
  
  