c  Mahc_Evolve: Numerical Evolution of the General Relativistic Hydro Equations
c  Copyright (C) 2001 Mark Miller

/*@@
  @file      apply_del_hydro.F
  @date      December 1999
  @author    Mark Miller
  @desc
      Add the del_hydro variables to the evolved variables.
  @enddesc
@@*/
#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h" 

/*@@
  @routine    apply_del_hydro
  @date       December 1999
  @author     Mark Miller
  @desc
      Add the del_hydro variables to the evolved variables
 
      Since different primative solvers have different personalities,
      the application criterion is different for each case. 
  @enddesc
  @calls
  @calledby
  @history
  @endhistory
@@*/ 

      subroutine apply_del_hydro(CCTK_FARGUMENTS,frac)
      implicit none
      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS
      CCTK_REAL frac

c     LOCAL VARS
      integer i,j,k,nx,ny,nz,CCTK_Equals
      CCTK_REAL apply_densmin

      nx = cctk_lsh(1)
      ny = cctk_lsh(2)
      nz = cctk_lsh(3) 
      if(CCTK_EQUALS(mahc_eos,'polytrope')) then
c     do the updates
      apply_densmin = mahc_densmin/100.0d0
      do k=1+mahc_stencil_size,nz-mahc_stencil_size
      do j=1+mahc_stencil_size,ny-mahc_stencil_size
      do i=1+mahc_stencil_size,nx-mahc_stencil_size
         if((dens(i,j,k) + frac * del_dens(i,j,k).gt. apply_densmin)) then
            dens(i,j,k) = dens(i,j,k) + frac * del_dens(i,j,k)
            sx(i,j,k)   = sx(i,j,k)   + frac * del_sx(i,j,k)
            sy(i,j,k)   = sy(i,j,k)   + frac * del_sy(i,j,k)
            sz(i,j,k)   = sz(i,j,k)   + frac * del_sz(i,j,k)
            tau(i,j,k)  = tau(i,j,k)  + frac * del_tau(i,j,k)
         else
            dens(i,j,k) = apply_densmin
            sx(i,j,k) = 0.0d0
            sy(i,j,k) = 0.0d0
            sz(i,j,k) = 0.0d0
         endif
      enddo
      enddo
      enddo

      else

c     do the updates
      apply_densmin = mahc_densmin/100.0d0
      do k=1+mahc_stencil_size,nz-mahc_stencil_size
      do j=1+mahc_stencil_size,ny-mahc_stencil_size
      do i=1+mahc_stencil_size,nx-mahc_stencil_size
         if((tau(i,j,k)  + frac * del_tau(i,j,k) .gt. 0.0d0).and.
     .      (dens(i,j,k) + frac * del_dens(i,j,k).gt. 0.0d0))then
         if((dens(i,j,k) + frac * del_dens(i,j,k).gt. apply_densmin)) then
            dens(i,j,k) = dens(i,j,k) + frac * del_dens(i,j,k)
            sx(i,j,k)   = sx(i,j,k)   + frac * del_sx(i,j,k)
            sy(i,j,k)   = sy(i,j,k)   + frac * del_sy(i,j,k)
            sz(i,j,k)   = sz(i,j,k)   + frac * del_sz(i,j,k)
            tau(i,j,k)  = tau(i,j,k)  + frac * del_tau(i,j,k)
         else
            dens(i,j,k) = apply_densmin
            sx(i,j,k) = 0.0d0
            sy(i,j,k) = 0.0d0
            sz(i,j,k) = 0.0d0
            tau(i,j,k) = apply_densmin * eps(i,j,k)
         endif
         endif
      enddo
      enddo
      enddo

      endif


      return
      end