* remove OMP pragmas from most routines, because they HEADmaster

  slow down pointwise calls

* make a copy of the press routine, call it OMP, add interface
  entry; use this routine for GF-wide calls.

* uglyness ensues



git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinEOS/EOS_Omni/trunk@84 8e189c6b-2ab8-4400-aa02-70a9cfce18b9

3 files changed, 176 insertions, 96 deletions

diff --git a/interface.ccl b/interface.ccl
index 57ef222..f652f69 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -22,6 +22,20 @@ void FUNCTION EOS_Omni_press(CCTK_INT IN eoskey,         \
 
 PROVIDES FUNCTION EOS_Omni_press WITH EOS_Omni_EOS_Press LANGUAGE Fortran
 
+void FUNCTION EOS_Omni_pressOMP(CCTK_INT IN eoskey,         \
+			     CCTK_INT IN havetemp,       \
+			     CCTK_REAL IN rf_precision,  \
+                             CCTK_INT IN npoints,        \
+			     CCTK_REAL IN ARRAY rho,     \
+			     CCTK_REAL INOUT ARRAY eps,  \
+			     CCTK_REAL INOUT ARRAY temp, \
+			     CCTK_REAL IN ARRAY ye,      \
+			     CCTK_REAL OUT ARRAY press,  \
+			     CCTK_INT OUT ARRAY keyerr,  \
+                             CCTK_INT OUT anyerr)
+
+PROVIDES FUNCTION EOS_Omni_pressOMP WITH EOS_Omni_EOS_PressOMP LANGUAGE Fortran
+
 
 void FUNCTION EOS_Omni_DPressByDEps(CCTK_INT IN eoskey,         \
 			     CCTK_INT IN havetemp,              \
diff --git a/src/EOS_Omni_MultiVarCalls.F90 b/src/EOS_Omni_MultiVarCalls.F90
index 5db911d..5d4a67d 100644
--- a/src/EOS_Omni_MultiVarCalls.F90
+++ b/src/EOS_Omni_MultiVarCalls.F90
@@ -47,8 +47,6 @@ subroutine EOS_Omni_EOS_short(eoskey,keytemp,rf_precision,npoints,&
   anyerr    = 0
   keyerr(:) = 0
 
-  !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xent,xcs2,&
-  !$OMP                    xdedt,xdpderho,xdpdrhoe,xmunu)
   do i=1,npoints
 
      xrho = rho(i) * inv_rho_gf
@@ -61,9 +59,7 @@ subroutine EOS_Omni_EOS_short(eoskey,keytemp,rf_precision,npoints,&
           keytemp,keyerr(i),rf_precision)
      
      if(keyerr(i).ne.0) then
-        !$OMP CRITICAL
         anyerr = 1
-        !$OMP END CRITICAL
      endif
      
      if(keytemp.eq.1) then
@@ -84,7 +80,6 @@ subroutine EOS_Omni_EOS_short(eoskey,keytemp,rf_precision,npoints,&
      munu(i)      = xmunu
      
   enddo
-  !$OMP END PARALLEL DO
 
 end subroutine EOS_Omni_EOS_short
 
@@ -123,45 +118,35 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
   case (1)
      ! polytropic EOS                                                                                      
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO
         do i=1,npoints
            dpdrhoe(i) = poly_k *  &
                 poly_gamma *        &
                 rho(i) ** (poly_gamma - 1.d0)
            dpderho(i) = 0.0d0
         enddo
-        !$OMP END PARALLEL DO
   case (2)
      ! gamma-law EOS                                                                                       
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO
         do i=1,npoints
            dpdrhoe(i) = (gl_gamma-1.0d0) * &
                 eps(i)
            dpderho(i) = (gl_gamma - 1.0d0) * &
                 rho(i)
         enddo
-        !$OMP END PARALLEL DO
   case (3)
      ! hybrid EOS                                                                                          
-     !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
-     !$OMP                     hybrid_dp_poly,hybrid_dp_th1,hybrid_dp_th2)
         do i=1,npoints
            if(rho(i).gt.hybrid_rho_nuc) then
               hybrid_local_gamma = hybrid_gamma2
@@ -187,9 +172,7 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
            dpdrhoe(i) = hybrid_dp_poly + hybrid_dp_th1 + hybrid_dp_th2
            dpderho(i) = (hybrid_gamma_th - 1.0d0) * rho(i)
         enddo
-        !$OMP END PARALLEL DO
     case (4)
-       !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xdpderho,xdpdrhoe)
         do i=1,npoints
            xrho = rho(i) * inv_rho_gf
            xtemp = temp(i)
@@ -200,9 +183,7 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
                 keytemp,keyerr(i),rf_precision)
 
            if(keyerr(i).ne.0) then
-              !$OMP CRITICAL
               anyerr = 1
-              !$OMP END CRITICAL
            endif
 
            if(keytemp.eq.1) then
@@ -257,45 +238,34 @@ subroutine EOS_Omni_EOS_dpderho_dpdrhoe(eoskey,keytemp,rf_precision,npoints,&
   case (1)
      ! polytropic EOS                                                                                      
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO
         do i=1,npoints
            depsdpress(i) = 1.0d0/(poly_gamma - 1.0d0)/rho(i)
            depsdrho(i) = depsdpress(i) * poly_k *  &
                 poly_gamma *         &
                 rho(i) ** (poly_gamma - 1.d0)
         enddo
-        !$OMP END PARALLEL DO
   case (2)
      ! gamma-law EOS                                                                                       
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO
         do i=1,npoints
            depsdpress(i) = 1.0/( (gl_gamma - 1.0d0) * &
                 rho(i))
            depsdrho(i) = -eps(i)/rho(i)
        enddo
-       !$OMP END PARALLEL DO
   case (3)
      ! hybrid EOS                                                                                          
-     !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
-     !$OMP                    hybrid_dp_poly,hybrid_dp_th1,&
-     !$OMP                    hybrid_dp_th2,xdpdrhoe,xdpderho)
         do i=1,npoints
            if(rho(i).gt.hybrid_rho_nuc) then
               hybrid_local_gamma = hybrid_gamma2
diff --git a/src/EOS_Omni_SingleVarCalls.F90 b/src/EOS_Omni_SingleVarCalls.F90
index fead6cb..be80a1e 100644
--- a/src/EOS_Omni_SingleVarCalls.F90
+++ b/src/EOS_Omni_SingleVarCalls.F90
@@ -47,6 +47,156 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
      case (1)
         ! polytropic EOS
         if(keytemp.eq.1) then
+           do i=1,npoints
+              eps(i) =  poly_k * &
+                   rho(i)**(poly_gamma) / &
+                   (poly_gamma - 1.0d0) / rho(i)
+           enddo
+        endif
+        do i=1,npoints
+           press(i) = poly_k * &
+                rho(i)**poly_gamma
+        enddo
+     case (2)
+        ! gamma-law EOS
+        if(keytemp.eq.1) then
+           do i=1,npoints
+              eps(i) = gl_k * &
+                   rho(i)**(gl_gamma) / &
+                   (gl_gamma - 1.0d0) / rho(i)
+           enddo
+        endif
+        do i=1,npoints
+           press(i) = (gl_gamma - 1.0d0) * rho(i) * eps(i)
+        enddo
+
+     case (3)
+        ! hybrid EOS
+        do i=1,npoints
+           if(rho(i).gt.hybrid_rho_nuc) then
+              hybrid_local_gamma = hybrid_gamma2
+              hybrid_local_k = hybrid_k2
+           else
+              hybrid_local_gamma = hybrid_gamma1
+              hybrid_local_k = hybrid_k1
+           endif
+           hybrid_p_poly = hybrid_local_k * &
+                     rho(i)**hybrid_local_gamma
+           hybrid_p_th = - hybrid_local_k * (hybrid_gamma_th - 1.d0) /      &
+                (hybrid_local_gamma - 1.0d0) * rho(i)**hybrid_local_gamma + &
+                (hybrid_gamma_th - 1.0d0) * rho(i) * eps(i) -                              &
+                (hybrid_gamma_th - 1.d0) * (hybrid_local_gamma - hybrid_gamma1) /          &
+                (hybrid_gamma1 - 1.d0) / (hybrid_gamma2 - 1.d0) *                          &
+                hybrid_k1 *                     &
+                hybrid_rho_nuc**(hybrid_gamma1 - 1.d0) * rho(i)
+           hybrid_p_th = max(zero, hybrid_p_th)
+           press(i) = hybrid_p_poly + hybrid_p_th
+        enddo
+
+     case (4)
+        ! nuc eos
+        do i=1,npoints
+           xrho = rho(i) * inv_rho_gf
+           xtemp = temp(i)
+           xye = ye(i)
+           xenr = eps(i) * inv_eps_gf
+           call nuc_eos_press_eps(xrho,xtemp,xye,xenr,xprs,& 
+                keytemp,keyerr(i),rf_precision)
+           
+           if(keyerr(i).ne.0) then
+              anyerr = 1
+           endif
+
+           if(keytemp.eq.1) then
+              eps(i) = xenr * eps_gf
+           else
+              temp(i) = xtemp
+           endif
+
+           press(i) = xprs * press_gf
+        enddo
+
+     case (5)
+        ! cold tabular EOS with gamma law
+        do i=1,npoints
+           if(rho(i).lt.coldeos_rhomin) then
+              press(i) = coldeos_low_kappa * rho(i)**coldeos_low_gamma
+              eps(i) = press(i) / (coldeos_low_gamma - 1.0) / rho(i) 
+              cycle
+           else if(rho(i).gt.coldeos_rhomax) then
+              keyerr(i) = 103
+              anyerr = 1
+           else
+              xrho = log10(rho(i))
+              ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
+           endif
+           ir = max(2, min(ir,coldeos_nrho))
+
+           gamma = (coldeos_gamma(ir) - coldeos_gamma(ir-1)) / &
+                (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                (xrho - coldeos_logrho(ir-1)) + coldeos_gamma(ir-1)
+
+           eps_cold = (coldeos_eps(ir) - coldeos_eps(ir-1)) / &
+                    (coldeos_logrho(ir) - coldeos_logrho(ir-1)) * &
+                    (xrho - coldeos_logrho(ir-1)) + coldeos_eps(ir-1)
+
+           if(keytemp.eq.0) then
+              eps_th = max(0.0d0,eps(i) - eps_cold)
+           else if (keytemp.eq.1) then
+              eps_th = 0.0d0
+              eps(i) = eps_cold
+           endif
+
+           press_cold = coldeos_kappa * rho(i)**gamma
+           press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th
+           press(i) = press_cold + press_th
+        enddo
+
+     case DEFAULT
+        write(warnstring,*) "eoskey ",eoskey," not implemented!"
+        call CCTK_WARN(0,warnstring)
+     end select
+
+end subroutine EOS_Omni_EOS_Press
+
+subroutine EOS_Omni_EOS_PressOMP(eoskey,keytemp,rf_precision,npoints,&
+                              rho,eps,temp,ye,press,keyerr,anyerr)
+
+  use EOS_Omni_Module
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+
+  CCTK_INT, intent(in)     :: eoskey,keytemp,npoints
+  CCTK_INT, intent(out)    :: keyerr(npoints)
+  CCTK_INT, intent(out)    :: anyerr
+  CCTK_REAL, intent(in)    :: rf_precision
+  CCTK_REAL, intent(in)    :: rho(npoints),ye(npoints)
+  CCTK_REAL, intent(inout) :: eps(npoints), temp(npoints)
+  CCTK_REAL, intent(out)   :: press(npoints)
+
+  ! local vars
+  integer          :: i
+  character(256)   :: warnstring
+  real*8           :: hybrid_local_gamma, hybrid_local_k, &
+                      hybrid_p_poly, hybrid_p_th
+  real*8,parameter :: zero = 0.0d0
+  ! temporary vars for nuc_eos
+  real*8           :: xrho,xye,xtemp,xenr,xent
+  real*8           :: xprs,xmunu,xcs2
+  real*8           :: xdedt,xdpderho,xdpdrhoe
+  ! temporary vars for coldeos + gamma law
+  integer :: ir
+  real*8 :: press_cold, press_th
+  real*8 :: eps_cold, eps_th
+  real*8 :: gamma
+
+  anyerr    = 0
+  keyerr(:) = 0
+
+  select case (eoskey)
+     case (1)
+        ! polytropic EOS
+        if(keytemp.eq.1) then
            !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) =  poly_k * &
@@ -177,8 +327,7 @@ subroutine EOS_Omni_EOS_Press(eoskey,keytemp,rf_precision,npoints,&
         call CCTK_WARN(0,warnstring)
      end select
 
-end subroutine EOS_Omni_EOS_Press
-
+end subroutine EOS_Omni_EOS_PressOMP
 
 
 subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
@@ -219,43 +368,32 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
      case (1)
         ! polytropic EOS
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO
         do i=1,npoints
            dpdrhoe(i) = poly_k *  &
                 poly_gamma *         & 
                 rho(i) ** (poly_gamma - 1.d0)
         enddo
-        !$OMP END PARALLEL DO
      case (2)
         ! gamma-law EOS
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO
         do i=1,npoints
            dpdrhoe(i) = (gl_gamma-1.0d0) * &
                 eps(i)
         enddo
-        !$OMP END PARALLEL DO
      case (3)
         ! hybrid EOS
-        !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
-        !$OMP                     hybrid_dp_poly, hybrid_dp_th1, &
-        !$OMP                     hybrid_dp_th2)
         do i=1,npoints
            if(rho(i).gt.hybrid_rho_nuc) then
               hybrid_local_gamma = hybrid_gamma2
@@ -281,10 +419,8 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
            dpdrhoe(i) = hybrid_dp_poly + max(0.0d0,hybrid_dp_th1 + hybrid_dp_th2)
 
         enddo
-        !$OMP END PARALLEL DO
+
      case (4)
-        !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs,xent,xcs2,xdedt,&
-        !$OMP                     xdpderho,xdpdrhoe,xmunu,h)
         do i=1,npoints
            xrho = rho(i) * inv_rho_gf
            xtemp = temp(i)
@@ -295,9 +431,7 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
                 keytemp,keyerr(i),rf_precision)
            
            if(keyerr(i).ne.0) then
-              !$OMP CRITICAL
               anyerr = 1
-              !$OMP END CRITICAL
            endif
 
            if(keytemp.eq.1) then
@@ -308,13 +442,11 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
 
            dpdrhoe(i) = xdpdrhoe * press_gf * inv_rho_gf
         enddo
-        !$OMP END PARALLEL DO
+        
      case (5)
         ! with the cold eos we have to assume P = P(rho), so
         ! by definition dPdrho is at constant internal energy
         ! and entropy (the latter, because T=0)
-        !$OMP PARALLEL DO PRIVATE(xrho,ir,gamma,eps_cold,eps_th, &
-        !$OMP                     cs2,press_cold,press_th)
         do i=1,npoints
            if(rho(i).lt.coldeos_rhomin) then
               dpdrhoe(i) = coldeos_low_kappa * coldeos_low_gamma * &
@@ -322,9 +454,7 @@ subroutine EOS_Omni_EOS_DPressByDRho(eoskey,keytemp,rf_precision,npoints,&
               cycle
            else if(rho(i).gt.coldeos_rhomax) then
               keyerr(i) = 103
-              !$OMP CRITICAL
               anyerr = 1
-              !$OMP END CRITICAL
            else
               xrho = log10(rho(i))
               ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
@@ -401,49 +531,38 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
      case (1)
         ! polytropic EOS
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO
         do i=1,npoints
            dpdepsrho(i) = 0.0d0
         enddo
-        !$OMP END PARALLEL DO
      case (2)
         ! gamma-law EOS
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO
         do i=1,npoints
            dpdepsrho(i) = (gl_gamma - 1.0d0) * &
                 rho(i)
         enddo
-        !$OMP END PARALLEL DO
      case (3)
         ! hybrid EOS
-        !$OMP PARALLEL DO
         do i=1,npoints
            dpdepsrho(i) = (hybrid_gamma_th - 1.0d0) * rho(i)
         enddo
-        !$OMP END PARALLEL DO
+
      case (4)
         ! nuc_eos
-        !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs,&
-        !$OMP                     xent,xcs2,xdedt,xdpderho,xdpdrhoe,&
-        !$OMP                     xmunu)
         do i=1,npoints
+
            xrho = rho(i) * inv_rho_gf
            xtemp = temp(i)
            xye = ye(i)
@@ -453,9 +572,7 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
                 keytemp,keyerr(i),rf_precision)
            
            if(keyerr(i).ne.0) then
-              !$OMP CRITICAL
               anyerr = 1
-              !$OMP END CRITICAL
            endif
 
            if(keytemp.eq.1) then
@@ -467,20 +584,17 @@ subroutine EOS_Omni_EOS_DPressByDEps(eoskey,keytemp,rf_precision,npoints,&
            dpdepsrho(i) = xdpderho * press_gf * inv_eps_gf 
 
         enddo
-        !$OMP END PARALLEL DO
+
      case (5)
         ! with the cold eos we have to assume P = P(rho), so
         ! only the gamma law has non-zero dPdeps
-        !$OMP PARALLEL DO PRIVATE(xrho,ir,eps_cold,eps_th)
         do i=1,npoints
            if(rho(i).lt.coldeos_rhomin) then
               dpdepsrho(i) = 0.0d0
               cycle
            else if(rho(i).gt.coldeos_rhomax) then
               keyerr(i) = 103
-              !$OMP CRITICAL
               anyerr = 1
-              !$OMP END CRITICAL
            else
               xrho = log10(rho(i))
               ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
@@ -547,34 +661,27 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
      case (1)
         ! polytropic EOS
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = poly_k * &
                    rho(i)**(poly_gamma) / &
                    (poly_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO
         do i=1,npoints
            xpress = poly_k * &
                 rho(i)**(poly_gamma)
            cs2(i) = poly_gamma * xpress / rho(i) / &
                 (1 + eps(i) + xpress/rho(i))
         enddo
-        !$OMP END PARALLEL DO
      case (2)
         ! gamma-law EOS
         if(keytemp.eq.1) then
-           !$OMP PARALLEL DO
            do i=1,npoints
               eps(i) = gl_k * &
                    rho(i)**(gl_gamma) / &
                    (gl_gamma - 1.0d0) / rho(i)
            enddo
-           !$OMP END PARALLEL DO
         endif
-        !$OMP PARALLEL DO PRIVATE(xpress,xdpdrhoe,xdpderho)
         do i=1,npoints
            xpress = (gl_gamma-1.0d0)*rho(i)*eps(i)
            xdpdrhoe = (gl_gamma-1.0d0)*eps(i)
@@ -582,12 +689,8 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
            cs2(i) = (xdpdrhoe + xpress * xdpderho / (rho(i)**2)) / &
                 (1.0d0 + eps(i) + xpress/rho(i))
         enddo
-        !$OMP END PARALLEL DO
      case(3)
         ! hybrid EOS
-        !$OMP PARALLEL DO PRIVATE(hybrid_local_gamma,hybrid_local_k,&
-        !$OMP                     hybrid_p_poly,hybrid_p_th,&
-        !$OMP                     xpress)
         do i=1,npoints
            if(rho(i).gt.hybrid_rho_nuc) then
               hybrid_local_gamma = hybrid_gamma2
@@ -611,12 +714,11 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
            cs2(i) = (hybrid_local_gamma * hybrid_p_poly + hybrid_gamma_th * hybrid_p_th) / &
                 rho(i) / (1.0d0 + eps(i) + xpress/rho(i))
         enddo
-        !$OMP END PARALLEL DO
      case(4)
         ! nuc_eos
-        !$OMP PARALLEL DO PRIVATE(xrho,xtemp,xye,xenr,xprs,xent,xcs2,&
-        !$OMP                     xdedt,xdpderho,xdpdrhoe,xmunu)
+
         do i=1,npoints
+           
            xrho = rho(i) * inv_rho_gf
            xtemp = temp(i)
            xye = ye(i)
@@ -626,9 +728,7 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
                 keytemp,keyerr(i),rf_precision)
            
            if(keyerr(i).ne.0) then
-              !$OMP CRITICAL
               anyerr = 1
-              !$OMP END CRITICAL
            endif
 
            if(keytemp.eq.1) then
@@ -639,15 +739,12 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
 
            cs2(i) = xcs2 * cliteinv2 / &
                 (1.0d0 + eps(i) + xprs * press_gf / rho(i))
+
         enddo
-        !$OMP END PARALLEL DO
      case (5)
         ! with the cold eos we have to assume P = P(rho), so
         ! by definition dPdrho is at constant internal energy
         ! and entropy (the latter, because T=0)
-        !$OMP PARALLEL DO PRIVATE(xprs,xrho,ir,gamma,cs2_cold,&
-        !$OMP                     eps_cold,eps_th,press_cold,press_th,&
-        !$OMP                     xdpdrhoe,xdpderho,cs2_th,h,h_cold)
         do i=1,npoints
            if(rho(i).lt.coldeos_rhomin) then
               xprs = coldeos_low_kappa * rho(i)**coldeos_low_gamma
@@ -658,10 +755,8 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
               cycle
            else if(rho(i).gt.coldeos_rhomax) then
               keyerr(i) = 103
-              !$OMP CRITICAL
               anyerr = 1
-              !$OMP END CRITICAL
-           else
+                 else
               xrho = log10(rho(i))
               ir = 2 + INT( (xrho - coldeos_logrho(1) - 1.0d-10) * coldeos_dlrhoi )
            endif
@@ -686,6 +781,7 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
               eps_th = 0.0d0
               eps(i) = eps_cold
            endif
+
            press_cold = coldeos_kappa * rho(i)**gamma
            press_th = coldeos_thfac*(coldeos_gammath - 1.0d0)*rho(i)*eps_th
 
@@ -698,7 +794,7 @@ subroutine EOS_Omni_EOS_cs2(eoskey,keytemp,rf_precision,npoints,&
            
            cs2(i) = (cs2_cold * h_cold + cs2_th) / h
         enddo
-        !$OMP END PARALLEL DO
+           
      case DEFAULT
         write(warnstring,*) "eoskey ",eoskey," not implemented!"
         call CCTK_WARN(0,warnstring)