* remove support for "General" EOS interface (EOSG_*)

* test suites pass with Intel 11 (no optimizaton) on bethe 
  and gcc 4.4.5




git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEvolve/GRHydro/trunk@225 c83d129a-5a75-4d5a-9c4d-ed3a5855bf45

1 files changed, 0 insertions, 754 deletions

diff --git a/src/GRHydro_Con2Prim.F90 b/src/GRHydro_Con2Prim.F90
index ae83741..4d0f015 100644
--- a/src/GRHydro_Con2Prim.F90
+++ b/src/GRHydro_Con2Prim.F90
@@ -1724,760 +1724,6 @@ subroutine Con2Prim_ptBoundsTracer(cons_tracer, tracer, rho, one_over_w_lorentz,
 end subroutine Con2Prim_ptBoundsTracer
 
 
- /*@@
-   @routine    Conservative2PrimitiveGeneral
-   @date       Sat Jan 26 01:08:33 2002
-   @author     Ian Hawke
-   @desc 
-   Wrapper routine that converts from conserved to primitive variables
-   at every grid cell centre. Converted to the EOSGeneral format
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
-@@*/
-
-subroutine Conservative2PrimitiveGeneral(CCTK_ARGUMENTS)
-  
-  USE GRHydro_Scalars
-
-  implicit none
-  
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_FUNCTIONS
-  
-  CCTK_INT :: i, j, k, itracer, nx, ny, nz, ierr
-  character(len=100) warnline
-
-  integer :: maxerrloc(3)
-  integer :: ii,jj,kk
-  CCTK_REAL :: maxerr
-
-
-  CCTK_INT :: count
-  CCTK_REAL :: s2, df, v2, w, vlowx, vlowy, vlowz 
-  CCTK_REAL :: temp1, temp2, temp3
-  CCTK_REAL :: udens, usx, usy, usz, utau
-  CCTK_REAL :: drhodpress, depsdpress
-  logical :: loop_condition
-
-  CCTK_REAL, dimension(:,:,:), allocatable :: f
-  logical, dimension(:,:,:), allocatable :: point_loop_condition
-
-  CCTK_INT :: type_bits, atmosphere
-  CCTK_INT :: type2_bits
-
-  CCTK_REAL :: local_min_tracer
-
-  CCTK_INT, dimension(3) :: loc, loc2
-
-  allocate(f(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
-           point_loop_condition(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3)),&
-           STAT=ierr)
-  if (ierr .ne. 0) then
-    call CCTK_WARN(0, "Allocation problems with f")
-  end if
-  
-  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
-  call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
-  type2_bits = -1
-
-  nx = cctk_lsh(1)
-  ny = cctk_lsh(2)
-  nz = cctk_lsh(3)
-  
-  count = 0
-  press_old = max(eosgeneral_pmin, press)
-
-  if (use_min_tracer .ne. 0) then
-    local_min_tracer = min_tracer
-  else
-    local_min_tracer = 0d0
-  end if
-
-  loop_condition = .true.
-
-!!$  Set up rho and epsilon
-
-  do k = 1, nz
-    do j = 1, ny
-      do i = 1, nx
-
-        if (evolve_tracer .ne. 0) then
-          do itracer=1,number_of_tracers
-            call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), &
-                 tracer(i,j,k,itracer), dens(i,j,k))
-            if (use_min_tracer .ne. 0) then
-              if (tracer(i,j,k,itracer) .le. local_min_tracer) then
-                tracer(i,j,k,itracer) = local_min_tracer
-              end if
-            end if
-          enddo
-        endif
-
-        if(evolve_Y_e.ne.0) then
-           Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
-        endif
-
-        if ( (dens(i,j,k).le.sqrt(GRHydro_Det(i,j,k)) * &
-             GRHydro_rho_min*(1.0d0+GRHydro_atmo_tolerance)) &
-             .or.(tau(i,j,k) .le. 0d0)) then
-
-           rho(i,j,k) = GRHydro_rho_min
-           dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k)
-           scon(i,j,k,:) = 0.d0
-           w_lorentz(i,j,k) = 1.d0
-           press(i,j,k) = 0.1d0 * eosgeneral_pmin
-!           eps(i,j,k) = press(i,j,k) / rho(i,j,k) ! Note that this should be improved
-           eps(i,j,k) = 1.0e-12
-           tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) * eps(i,j,k)
-
-           SpaceMask_SetStateBitsF90(space_mask, i, j, k, type_bits, atmosphere)
-           atmosphere_mask(i,j,k) = 1
-
-        end if
-        
-        udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
-        usx = scon(i,j,k,1) / sqrt(GRHydro_det(i,j,k))
-        usy = scon(i,j,k,2) / sqrt(GRHydro_det(i,j,k))
-        usz = scon(i,j,k,3) / sqrt(GRHydro_det(i,j,k))
-        utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k))
-        s2 = usx*usx*GRHydro_uxx(i,j,k) + &
-             usy*usy*GRHydro_uyy(i,j,k) + &
-             usz*usz*GRHydro_uzz(i,j,k) + &
-             2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
-             2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
-             2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
-!!$ Check that the variables have a chance of being physical
-  
-        if( (utau + press_old(i,j,k) + udens)**2 - s2 .le. 0.0d0) then
-          if(GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
-            !$OMP CRITICAL
-            call CCTK_WARN(1,'Con2PrimGeneral: variables unphysical!')
-            write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,g15.6)') 'undensitized energy: ',utau
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,g15.6)') 'undensitized density: ',udens
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,g15.6)') 'undensitized s2: ',s2
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,g15.6)') 'pressure guess: ',press_old(i,j,k)
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,i4)') 'refinement level: ',GRHydro_reflevel
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',&
-                 x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k)
-            call CCTK_WARN(1,warnline)
-            call CCTK_WARN(GRHydro_c2p_warnlevel, "Unphysical variables")
-            !$OMP END CRITICAL
-            exit
-          else
-            !$OMP CRITICAL
-            write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',GRHydro_reflevel
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',&
-                 x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k)
-            call CCTK_WARN(1,warnline)
-            !$OMP END CRITICAL
-            exit
-          endif
-        endif
-
-        temp1 = (utau + udens + press_old(i,j,k))
-        temp2 = temp1**2 - s2
-        temp3 = sqrt(temp2)
-
-        rho(i,j,k) = udens * temp3 / temp1
-        w_lorentz(i,j,k) = temp1 / temp3
-
-!!$ Eps not previously set?
-        eps(i,j,k) = &
-             (temp3 - press_old(i,j,k) * w_lorentz(i,j,k) - udens) / &
-             udens
-
-      end do
-    end do
-  end do
-  
-  ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
-  where (press .le. eosgeneral_pmin)
-    press = eosgeneral_pmin
-    rho = GRHydro_rho_min
-    dens = sqrt(GRHydro_det) * rho
-    scon(:,:,:,1) = 0.d0
-    scon(:,:,:,2) = 0.d0
-    scon(:,:,:,3) = 0.d0
-    w_lorentz = 1d0
-!    eps = press / rho
-    eps = 1.0e-12
-    tau = sqrt(GRHydro_det) * rho * eps
-  end where
-  f = press_old - press
-  press_new = press_old
-
-  do while (loop_condition)
-
-    count = count + 1
-
-    if (count > GRHydro_countmax) then
-       if(GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
-          !$OMP CRITICAL
-          call CCTK_WARN(1, 'Con2PrimGeneral: error: did not converge in ')
-          write(warnline,'(a20,i12,a10)') '              ',&
-               GRHydro_countmax,' steps'
-          call CCTK_WARN(1,warnline)
-          write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
-          call CCTK_WARN(1,warnline)
-          maxerr =  maxval(abs(press_new - press_old) & 
-               / abs(press_new))
-          maxerrloc = maxloc(abs(press_new - press_old) & 
-               / abs(press_new)) 
-          ii = maxerrloc(1)
-          jj = maxerrloc(2)
-          kk = maxerrloc(3)
-          write(warnline,'(1P10E15.6)') maxerr
-          call CCTK_WARN(1,warnline)
-          write(warnline,'(1P10E15.6)') press_new(ii,jj,kk), press_old(ii,jj,kk)
-          call CCTK_WARN(1,warnline)
-          write(warnline,'(3i4)') ii,jj,kk
-          call CCTK_WARN(1,warnline)
-          write(warnline,'(1P10E15.6)') x(ii,jj,kk), &
-               y(ii,jj,kk), &
-               z(ii,jj,kk)
-          call CCTK_WARN(1,warnline)
-          call CCTK_WARN(0,"Abort!")
-          !$OMP END CRITICAL
-          exit 
-       else
-         !$OMP CRITICAL
-         write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',&
-              GRHydro_reflevel
-         call CCTK_WARN(1,warnline)
-         !$OMP END CRITICAL
-         exit
-       endif
-    end if
-
-    do k = 1, nz 
-      do j = 1, ny 
-        do i = 1, nx
-
-          !do not compute if in atmosphere or in an excised region
-          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
-             GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) &
-               then
-             press_new(i,j,k) = press(i,j,k)
-             press_old(i,j,k) = press(i,j,k)
-!             cycle
-          endif
-
-
-          ! All right, is this point already converged?!?
-          loop_condition = abs(press_new(i,j,k)-press_old(i,j,k)) / abs(press_new(i,j,k)) > &
-               GRHydro_perc_ptol
-
-          loop_condition = loop_condition .and. &
-               abs(press_new(i,j,k)-press_old(i,j,k)) > GRHydro_del_ptol
-
-          loop_condition = loop_condition .or. (count < GRHydro_countmin) &
-               .or. (count < 2)  ! we have to at least go once through c2p.
-
-          point_loop_condition(i,j,k) = loop_condition
-
-          if(.not.loop_condition) cycle
-
-          udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
-          usx   = scon(i,j,k,1)   / sqrt(GRHydro_det(i,j,k))
-          usy   = scon(i,j,k,2)   / sqrt(GRHydro_det(i,j,k))
-          usz   = scon(i,j,k,3)   / sqrt(GRHydro_det(i,j,k))
-          utau  = tau(i,j,k)  / sqrt(GRHydro_det(i,j,k))
-          s2 = usx*usx*GRHydro_uxx(i,j,k) + &
-               usy*usy*GRHydro_uyy(i,j,k) + &
-               usz*usz*GRHydro_uzz(i,j,k) + &
-               2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
-               2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
-               2.d0*usy*usz*GRHydro_uyz(i,j,k)
-          
-          temp1 = (utau + udens + press_new(i,j,k))
-          temp2 = temp1**2 - s2
-          temp3 = sqrt(temp2)
-          
-          drhodpress = udens * s2 / (temp3 * temp1**2)
-          depsdpress = press_new(i,j,k) * s2 / (rho(i,j,k) * temp1 * temp2)
-
-          df = 1.d0 - eos_dpdeps_temp(i,j,k) * depsdpress - &
-                      eos_dpdrho_temp(i,j,k) * drhodpress
-         
-          press_old(i,j,k) = press_new(i,j,k) 
-                  
-          press_new(i,j,k) = press_old(i,j,k) - f(i,j,k) / df
-          if (press_new(i,j,k) .le. eosgeneral_pmin) then
-            press_new(i,j,k) = eosgeneral_pmin
-            rho(i,j,k) = GRHydro_rho_min
-            dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k)
-            scon(i,j,k,:) = 0.d0
-            w_lorentz(i,j,k) = 1d0
-            eps(i,j,k) = press(i,j,k) / rho(i,j,k)
-            tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k) * eps(i,j,k)
-            udens = rho(i,j,k)
-            utau = rho(i,j,k) * eps(i,j,k)
-            s2 = 0.d0
-          end if
-         
-          temp1 = (utau + udens + press_new(i,j,k))
-          temp2 = temp1**2 - s2
-          temp3 = sqrt(temp2)
-          
-          rho(i,j,k) = udens * temp3 / temp1
-          w_lorentz(i,j,k) = temp1 / temp3
-          eps(i,j,k) = &
-               (temp3 - press_new(i,j,k) * w_lorentz(i,j,k) - udens) / &
-               udens         
-
-          point_loop_condition(i,j,k) = &
-               (abs(press_new(i,j,k) - press_old(i,j,k)) / abs(press_new(i,j,k)) > &
-                      GRHydro_perc_ptol)
-          point_loop_condition(i,j,k) = point_loop_condition(i,j,k) .and. &
-               (abs(press_new(i,j,k) - press_old(i,j,k))  > &
-                      GRHydro_del_ptol)          
-          
-         end do
-      end do
-    end do
-    
-    ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
-    press = max(press, eosgeneral_pmin)
-
-    f = press_new - press
-
-    loop_condition = any(point_loop_condition)
-
-    loop_condition = loop_condition .or. &
-         (count < GRHydro_countmin) 
-
-   end do
-
-  deallocate(point_loop_condition)
-
-
-  do count = 1, GRHydro_polish
-
-    do k = 1, nz 
-      do j = 1, ny 
-        do i = 1, nx
-
-          !do not compute if in atmosphere or in an excised region
-          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
-              GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
- 
-          f(i,j,k) = press_new(i,j,k) - press(i,j,k)
-
-          udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
-          usx   = scon(i,j,k,1)   / sqrt(GRHydro_det(i,j,k))
-          usy   = scon(i,j,k,2)   / sqrt(GRHydro_det(i,j,k))
-          usz   = scon(i,j,k,3)   / sqrt(GRHydro_det(i,j,k))
-          utau  = tau(i,j,k)  / sqrt(GRHydro_det(i,j,k))
-          s2 = usx*usx*GRHydro_uxx(i,j,k) + &
-               usy*usy*GRHydro_uyy(i,j,k) + &
-               usz*usz*GRHydro_uzz(i,j,k) + &
-               2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
-               2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
-               2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
-          temp1 = (utau + udens + press_new(i,j,k))
-          temp2 = temp1**2 - s2
-          temp3 = sqrt(temp2)
-
-          drhodpress = udens * s2 / (temp3 * temp1**2)
-          depsdpress = press_new(i,j,k) * s2 / (rho(i,j,k) * temp1 * temp2)
-
-          df = 1.d0 - eos_dpdeps_temp(i,j,k) * depsdpress - &
-                      eos_dpdrho_temp(i,j,k) * drhodpress
-
-          press_new(i,j,k) = max(press_new(i,j,k) - f(i,j,k) / df, eosgeneral_pmin)
-
-          temp1 = (utau + udens + press_new(i,j,k))
-          temp2 = temp1**2 - s2
-          temp3 = sqrt(temp2)
-
-          rho(i,j,k) = udens * temp3 / temp1
-          w_lorentz(i,j,k) = temp1 / temp3
-          eps(i,j,k) = &
-               (temp3 - press_new(i,j,k) * w_lorentz(i,j,k) - udens) / &
-               udens
-
-        end do
-      end do
-    end do
-    
-    ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
-    press = max(press, eosgeneral_pmin)
-  
-  end do
-  
-  
-  do k = 1, nz 
-    do j = 1, ny 
-      do i = 1, nx
-        
-        !do not compute if in atmosphere or in an excised region
-        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
-            GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
-
-!!$         Note that in the following we enforce eps > 0
-!!$         This makes the later check redundant
-       
-        if ( (rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) .or. &
-             (eps(i,j,k) .lt. 0d0) ) then
-            
-          rho(i,j,k) = GRHydro_rho_min
-          udens = rho(i,j,k)
-          utau = tau(i,j,k) / sqrt(GRHydro_det(i,j,k))
-          dens(i,j,k) = sqrt(GRHydro_det(i,j,k)) * rho(i,j,k)
-          eps(i,j,k) = utau / udens
-          scon(i,j,k,:) = 0.d0
-          w_lorentz(i,j,k) = 1.d0
-          usx = 0.d0
-          usy = 0.d0
-          usz = 0.d0
-
-        else
-
-          udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
-          usx   = scon(i,j,k,1)   / sqrt(GRHydro_det(i,j,k))
-          usy   = scon(i,j,k,2)   / sqrt(GRHydro_det(i,j,k))
-          usz   = scon(i,j,k,3)   / sqrt(GRHydro_det(i,j,k))
-          utau  = tau(i,j,k)  / sqrt(GRHydro_det(i,j,k))
-
-        end if
-
-        vlowx = usx / &
-             ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
-               w_lorentz(i,j,k)**2)
-        vlowy = usy / &
-             ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
-               w_lorentz(i,j,k)**2)
-        vlowz = usz / &
-             ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
-               w_lorentz(i,j,k)**2)
-        vel(i,j,k,1) = GRHydro_uxx(i,j,k) * vlowx + &
-                       GRHydro_uxy(i,j,k) * vlowy + &
-                       GRHydro_uxz(i,j,k) * vlowz
-        vel(i,j,k,2) = GRHydro_uxy(i,j,k) * vlowx + &
-                       GRHydro_uyy(i,j,k) * vlowy + &
-                       GRHydro_uyz(i,j,k) * vlowz
-        vel(i,j,k,3) = GRHydro_uxz(i,j,k) * vlowx + &
-                       GRHydro_uyz(i,j,k) * vlowy + &
-                       GRHydro_uzz(i,j,k) * vlowz
-  
-!!$  If eps negative then complain vociferously
-
-        if (eps(i,j,k) .le. 0.d0) then
-          if(GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then          
-            !$OMP CRITICAL
-            call CCTK_WARN(1,'Con2PrimGeneral: stopping the code.')
-            call CCTK_WARN(1, '   specific internal energy just went below 0! ')
-            write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a20,3i6)') 'ijk location: ',&
-                 i,j,k
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a20,1P10E15.6)') 'xyz location: ',&
-                 x(i,j,k),y(i,j,k),z(i,j,k)
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,i5)') 'refinement level: ',GRHydro_reflevel
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a20,1P1E15.6)') 'radius: ',r(i,j,k)
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a20,1P1E15.6)') 'velocities: ',&
-                 vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3)
-            call CCTK_WARN(1,warnline)
-            call CCTK_WARN(GRHydro_c2p_warnlevel, \
-            "Specific internal energy negative")
-            !$OMP END CRITICAL
-            exit
-          else
-            !$OMP CRITICAL
-            write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but \
-                        I was told to ignore level: ',GRHydro_reflevel
-            call CCTK_WARN(1,warnline)
-            write(warnline,'(a25,4g15.6)') 'coordinates: x,y,z,r:',&
-                 x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k)
-            call CCTK_WARN(1,warnline)
-            !$OMP END CRITICAL
-            exit
-          endif
-        endif
-
-      end do
-    end do
-  end do
-
-  deallocate(f)
-
-end subroutine Conservative2PrimitiveGeneral
-
-
- /*@@
-   @routine    Conservative2PrimitivePolytypeGeneral
-   @date       Sat Jan 26 01:08:33 2002
-   @author     Ian Hawke
-   @desc 
-   Wrapper routine that converts from conserved to primitive variables
-   at every grid cell centre. Converted to the EOSGeneral format
-   @enddesc 
-   @calls     
-   @calledby   
-   @history 
-   @endhistory 
-
-@@*/
-
-subroutine Con2PrimPolytypeGeneral(CCTK_ARGUMENTS)
-  
-  USE GRHydro_Scalars
-
-  implicit none
-  
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_FUNCTIONS
-  
-  CCTK_INT :: i, j, k, itracer, nx, ny, nz
-  character(len=100) warnline
-
-  CCTK_INT :: count, ierr
-  CCTK_REAL :: s2, f, df, v2, w, vlowx, vlowy, vlowz 
-  CCTK_REAL :: temp1, temp2, temp3
-  CCTK_REAL :: udens, usx, usy, usz
-  CCTK_REAL :: enthalpy, denthalpy
-  logical :: loop_condition
-  
-  CCTK_INT :: type_bits
-  CCTK_INT :: atmosphere
-  CCTK_INT :: type2_bits
-
-  CCTK_REAL :: local_min_tracer
-                                                                               
-  call SpaceMask_GetTypeBits(type_bits, "Hydro_Atmosphere")
-  call SpaceMask_GetStateBits(atmosphere, "Hydro_Atmosphere", "in_atmosphere")
-  type2_bits = -1
-
-  nx = cctk_lsh(1)
-  ny = cctk_lsh(2)
-  nz = cctk_lsh(3)
-  
-  if (use_min_tracer .ne. 0) then
-    local_min_tracer = min_tracer
-  else
-    local_min_tracer = 0d0
-  end if
-  
-!!$  In what follows press_temp is really rho_temp
-
-  count = 0
-  press_new = max(GRHydro_rho_min, rho)
-
-  loop_condition = .true.
-  
-  ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
-
-!!$  Set up W and epsilon
-
-  do k = 1, nz 
-    do j = 1, ny 
-      do i = 1, nx       
-
-        if (evolve_tracer .ne. 0) then
-          do itracer=1,number_of_tracers
-            call Con2Prim_ptTracer(cons_tracer(i,j,k,itracer), & 
-                 tracer(i,j,k,itracer), dens(i,j,k))
-
-            if (use_min_tracer .ne. 0) then
-              if (tracer(i,j,k,itracer) .le. local_min_tracer) then
-                tracer(i,j,k,itracer) = local_min_tracer
-              end if
-            end if
-            
-          enddo
-        end if
-
-        if(evolve_Y_e.ne.0) then
-           Y_e(i,j,k) = Y_e_con(i,j,k) / dens(i,j,k)
-        endif
-
-        udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
-        usx   = scon(i,j,k,1)   / sqrt(GRHydro_det(i,j,k))
-        usy   = scon(i,j,k,2)   / sqrt(GRHydro_det(i,j,k))
-        usz   = scon(i,j,k,3)   / sqrt(GRHydro_det(i,j,k))
-        s2 = usx*usx*GRHydro_uxx(i,j,k) + &
-             usy*usy*GRHydro_uyy(i,j,k) + &
-             usz*usz*GRHydro_uzz(i,j,k) + &
-             2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
-             2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
-             2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
-        enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k)
-        w_lorentz(i,j,k) = sqrt( 1.d0 + s2 / ( (udens*enthalpy)**2 ) )
-  
-      end do
-    end do
-  end do
-
-  do while (loop_condition)
-
-    count = count + 1
-
-    if (count > GRHydro_countmax) then 
-      if (GRHydro_reflevel.ge.GRHydro_c2p_warn_from_reflevel) then
-        !$OMP CRITICAL
-        call CCTK_WARN(1, 'Con2PrimPolytypeGeneral: error: did not converge in ')
-        write(warnline,'(a20,g12.7,a10)') '              ',&
-             GRHydro_countmax,' steps'
-        call CCTK_WARN(1,warnline)
-        write(warnline,'(a28,i2)') 'on carpet reflevel: ',GRHydro_reflevel
-        call CCTK_WARN(1,warnline)
-        call CCTK_WARN(GRHydro_c2p_warnlevel, "Did not converge")
-        !$OMP END CRITICAL
-        exit
-      else 
-        !$OMP CRITICAL
-        write(warnline,'(a60,i2)') 'Con2Prim: eps negative, but I was told to ignore level: ',GRHydro_reflevel
-        call CCTK_WARN(1,warnline)
-        !$OMP END CRITICAL
-        exit
-      endif
-    end if
-
-    do k = 1, nz 
-      do j = 1, ny 
-        do i = 1, nx
-
-          !do not compute if in atmosphere or in an excised region
-          if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
-              GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
-
-          enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k)
-
-          tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * &
-               ( press_new(i,j,k) * enthalpy * w_lorentz(i,j,k)**2 - &
-                 press(i,j,k) ) - dens(i,j,k)
-
-          udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
-          usx   = scon(i,j,k,1)   / sqrt(GRHydro_det(i,j,k))
-          usy   = scon(i,j,k,2)   / sqrt(GRHydro_det(i,j,k))
-          usz   = scon(i,j,k,3)   / sqrt(GRHydro_det(i,j,k))
-          s2 = usx*usx*GRHydro_uxx(i,j,k) + &
-               usy*usy*GRHydro_uyy(i,j,k) + &
-               usz*usz*GRHydro_uzz(i,j,k) + &
-               2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
-               2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
-               2.d0*usy*usz*GRHydro_uyz(i,j,k)
-
-          f = press_new(i,j,k) * w_lorentz(i,j,k) - udens
-
-          denthalpy = eos_dpdrho_temp(i,j,k) / press_new(i,j,k)
-
-          df = w_lorentz(i,j,k) - press_new(i,j,k) * s2 * denthalpy / &
-               ( w_lorentz(i,j,k) * (udens**2) * (enthalpy**3) )
-
-          rho(i,j,k) = press_new(i,j,k)
-          press_new(i,j,k) = max(press_new(i,j,k) - f / df, GRHydro_rho_min)
-
-          enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / press_new(i,j,k)
-          w_lorentz(i,j,k) = sqrt( 1.d0 + s2 / ( (udens*enthalpy)**2 ) )
-
-        end do
-      end do
-    end do
-    
-    ierr = EOS_SetGFs(cctkGH, EOS_Con2PrimCall)
-
-    loop_condition = (maxval(abs(rho - press_new) / abs(press_new)) > &
-                      GRHydro_perc_ptol)
-    loop_condition = loop_condition .and. &
-         (maxval(abs(rho - press_new)) > GRHydro_del_ptol)
-    loop_condition = loop_condition .or. &
-         (count < GRHydro_countmin) 
-  
-  end do  
-  
-  do k = 1, nz 
-    do j = 1, ny 
-      do i = 1, nx
-
-        !do not compute if in atmosphere or in an excised region
-        if (SpaceMask_CheckStateBitsF90(space_mask, i, j, k, type_bits, atmosphere) .or. &
-            GRHydro_enable_internal_excision /= 0 .and. (hydro_excision_mask(i,j,k) .ne. 0)) cycle
- 
-        udens = dens(i,j,k) / sqrt(GRHydro_det(i,j,k))
-
-        if (rho(i,j,k) .le. GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then        
-
-          rho(i,j,k) = GRHydro_rho_min
-          scon(i,j,k,:) = 0.d0
-          w_lorentz(i,j,k) = 1.d0
-          usx = 0.d0
-          usy = 0.d0
-          usz = 0.d0
-          s2 = 0.d0
-          
-        else
-          
-          usx   = scon(i,j,k,1)   / sqrt(GRHydro_det(i,j,k))
-          usy   = scon(i,j,k,2)   / sqrt(GRHydro_det(i,j,k))
-          usz   = scon(i,j,k,3)   / sqrt(GRHydro_det(i,j,k))
-          s2 = usx*usx*GRHydro_uxx(i,j,k) + &
-               usy*usy*GRHydro_uyy(i,j,k) + &
-               usz*usz*GRHydro_uzz(i,j,k) + &
-               2.d0*usx*usy*GRHydro_uxy(i,j,k) + &
-               2.d0*usx*usz*GRHydro_uxz(i,j,k) + &
-               2.d0*usy*usz*GRHydro_uyz(i,j,k)
-          
-        end if
-        
-        enthalpy = 1.d0 + eps(i,j,k) + press(i,j,k) / rho(i,j,k)
-        w_lorentz(i,j,k) = sqrt( 1.d0 + s2 / ( (udens*enthalpy)**2 ) )
-        tau(i,j,k) = sqrt(GRHydro_det(i,j,k)) * &
-             ( (rho(i,j,k) * enthalpy) * w_lorentz(i,j,k)**2 - &
-                press(i,j,k) ) - dens(i,j,k)
-        vlowx = usx / &
-             ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
-               w_lorentz(i,j,k)**2)
-        vlowy = usy / &
-             ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
-               w_lorentz(i,j,k)**2)
-        vlowz = usz / &
-             ( (rho(i,j,k) + rho(i,j,k) * eps(i,j,k) + press(i,j,k)) * &
-               w_lorentz(i,j,k)**2)
-        vel(i,j,k,1) = GRHydro_uxx(i,j,k) * vlowx + &
-                       GRHydro_uxy(i,j,k) * vlowy + &
-                       GRHydro_uxz(i,j,k) * vlowz
-        vel(i,j,k,2) = GRHydro_uxy(i,j,k) * vlowx + &
-                       GRHydro_uyy(i,j,k) * vlowy + &
-                       GRHydro_uyz(i,j,k) * vlowz
-        vel(i,j,k,3) = GRHydro_uxz(i,j,k) * vlowx + &
-                       GRHydro_uyz(i,j,k) * vlowy + &
-                       GRHydro_uzz(i,j,k) * vlowz
-  
-!!$  If eps negative then complain vociferously
-
-        if (eps(i,j,k) .lt. 0.0d0) then
-          call CCTK_WARN(1, "c2p failed: eps went negative")
-          GRHydro_C2P_failed(i,j,k) = 1
-        endif
-
-      end do
-    end do
-  end do
-  
-end subroutine Con2PrimPolytypeGeneral
-
 ! subroutines to manage the C2P failure mask
 subroutine reset_GRHydro_C2P_failed(CCTK_ARGUMENTS)