MERGE divergence cleaning feature into trunk.

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

1 files changed, 41 insertions, 41 deletions

diff --git a/src/GRHydro_Con2PrimM.F90 b/src/GRHydro_Con2PrimM.F90
index 582eec5..e062976 100644
--- a/src/GRHydro_Con2PrimM.F90
+++ b/src/GRHydro_Con2PrimM.F90
@@ -129,8 +129,7 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
            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(det)*GRHydro_rho_min*(1.d0+GRHydro_atmo_tolerance) ) then
               
               b2=gxx(i,j,k)*Bvec(i,j,k,1)**2+gyy(i,j,k)*Bvec(i,j,k,2)**2+gzz(i,j,k)*Bvec(i,j,k,3)**2+ &
@@ -163,11 +162,11 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
            if(evolve_temper.eq.0) then
               call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam, dens(i,j,k), &
                    scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), tau(i,j,k), &
-                   rho(i,j,k),&
-                   vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
+                   Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),&
+                   vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),&
+                   Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(i,j,k),&
                    gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
                    uxx,uxy,uxz,uyy,uyz,uzz,det, &
-                   Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
                    epsnegative,GRHydro_C2P_failed(i,j,k))
            else
               stop "Please implement finite T con2prim routine in MHD part!"
@@ -239,13 +238,13 @@ subroutine Conservative2PrimitiveM(CCTK_ARGUMENTS)
 
               call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), &
                    scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, &
-                   rho(i,j,k),&
-                   vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
+                   Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),&
+                   vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),&
+                   Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(i,j,k),&
                    gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
                    uxx,uxy,uxz,uyy,uyz,uzz,det, &
-                   Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
                    epsnegative,GRHydro_C2P_failed(i,j,k))
-              
+             
 #endif          
               
            end if
@@ -294,7 +293,7 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
        uxxr, uxyr, uxzr, uyyr, uyzr, uzzr, pmin, epsmin
   CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,&
        gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr
-  CCTK_REAL :: b2minus, b2plus, local_gam, local_pgam,local_K,sc
+  CCTK_REAL :: b2minus, b2plus, local_gam, local_pgam,local_K,scminus,scplus
   CCTK_INT :: epsnegative
   character(len=100) warnline
  
@@ -387,11 +386,11 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
            
            call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam,densminus(i,j,k), &
                 sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), tauminus(i,j,k), &
-                rhominus(i,j,k),&
-                velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),&
+                Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), rhominus(i,j,k),&
+                velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),&
+                Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),&
                 gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
                 uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
-                Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,&
                 epsnegative,GRHydro_C2P_failed(i,j,k))
            
            if (epsnegative .ne. 0) then
@@ -407,16 +406,17 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
               call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                    xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
               local_pgam=log(xpress/local_K)/log(xrho)
-              sc = local_K*densminus(i,j,k)
+              scminus = local_K*densminus(i,j,k)
 
               call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densminus(i,j,k), &
-                   sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), sc, &
-                   rhominus(i,j,k),&
-                   velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),&
+                   sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus, &
+                   Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), rhominus(i,j,k),&
+                   velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),&
+                   Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),&
                    gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
                    uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
-                   Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,&
                    epsnegative,GRHydro_C2P_failed(i,j,k))
+          
            end if
            
            if (epsminus(i,j,k) .lt. 0.0d0) then
@@ -443,13 +443,13 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
            epsnegative = 0
            call GRHydro_Con2PrimM_pt(GRHydro_eos_handle, local_gam, densplus(i,j,k), &
                 sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), tauplus(i,j,k),&
-                rhoplus(i,j,k),&
-                velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),&
+                Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(i,j,k), rhoplus(i,j,k),&
+                velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+                Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),&
                 gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
                 uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
-                Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,&
                 epsnegative,GRHydro_C2P_failed(i,j,k))
-           
+ 
            if (epsnegative .ne. 0) then
               !$OMP CRITICAL
               call CCTK_WARN(GRHydro_NaN_verbose+2, 'Specific internal energy just went below 0, trying polytype!!')
@@ -463,15 +463,15 @@ subroutine Conservative2PrimitiveBoundsM(CCTK_ARGUMENTS)
               call EOS_Omni_press(GRHydro_polytrope_handle,keytemp,GRHydro_eos_rf_prec,n,&
                    xrho,1.0d0,xtemp,xye,xpress,keyerr,anyerr)
               local_pgam=log(xpress/local_K)/log(xrho)
-              sc = local_K*densplus(i,j,k)
+              scplus = local_K*densplus(i,j,k)
 
               call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, densplus(i,j,k), &
-                   sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), sc,&
-                   rhoplus(i,j,k),&
-                   velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),&
+                   sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,&
+                   Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(i,j,k), rhoplus(i,j,k),&
+                   velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+                   Bvecxplus(i,j,k), Bvecyplus(i,j,k), Bveczplus(i,j,k),b2plus, w_lorentzplus(i,j,k),&
                    gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
                    uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
-                   Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,&
                    epsnegative,GRHydro_C2P_failed(i,j,k))
            end if
            
@@ -610,14 +610,14 @@ subroutine Conservative2PrimitivePolytypeM(CCTK_ARGUMENTS)
            sc = local_K*dens(i,j,k)
            
            call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam, dens(i,j,k), &
-                scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc,&
-                rho(i,j,k),&
-                vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),w_lorentz(i,j,k),&
-                gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
-                uxx,uxy,uxz,uyy,uyz,uzz,det, &
-                Bvec(i,j,k,1), Bvec(i,j,k,2),Bvec(i,j,k,3),b2,&
-                epsnegative,GRHydro_C2P_failed(i,j,k))
-
+                   scon(i,j,k,1),scon(i,j,k,2),scon(i,j,k,3), sc, &
+                   Bcons(i,j,k,1), Bcons(i,j,k,2), Bcons(i,j,k,3), rho(i,j,k),&
+                   vel(i,j,k,1),vel(i,j,k,2),vel(i,j,k,3),eps(i,j,k),press(i,j,k),&
+                   Bvec(i,j,k,1), Bvec(i,j,k,2), Bvec(i,j,k,3),b2, w_lorentz(i,j,k),&
+                   gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k), &
+                   uxx,uxy,uxz,uyy,uyz,uzz,det, &
+                   epsnegative,GRHydro_C2P_failed(i,j,k))
+            
         end do
      end do
   end do
@@ -723,20 +723,20 @@ subroutine Con2PrimBoundsPolytypeM(CCTK_ARGUMENTS)
            
            call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densminus(i,j,k), &
                 sxminus(i,j,k),syminus(i,j,k),szminus(i,j,k), scminus,&
-                rhominus(i,j,k),&
-                velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),w_lorentzminus(i,j,k),&
+                Bconsxminus(i,j,k), Bconsyminus(i,j,k), Bconszminus(i,j,k), rhominus(i,j,k),&
+                velxminus(i,j,k),velyminus(i,j,k),velzminus(i,j,k),epsminus(i,j,k),pressminus(i,j,k),&
+                Bvecxminus(i,j,k), Bvecyminus(i,j,k), Bveczminus(i,j,k),b2minus, w_lorentzminus(i,j,k),&
                 gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, &
                 uxxl,uxyl,uxzl,uyyl,uyzl,uzzl,avg_detl, &
-                Bvecxminus(i,j,k), Bvecyminus(i,j,k),Bveczminus(i,j,k),b2minus,&
                 epsnegative,GRHydro_C2P_failed(i,j,k))
 
            call GRHydro_Con2PrimM_Polytype_pt(GRHydro_eos_handle, local_pgam,densplus(i,j,k), &
                 sxplus(i,j,k),syplus(i,j,k),szplus(i,j,k), scplus,&
-                rhoplus(i,j,k),&
-                velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),w_lorentzplus(i,j,k),&
+                Bconsxplus(i,j,k), Bconsyplus(i,j,k), Bconszplus(i,j,k), rhoplus(i,j,k),&
+                velxplus(i,j,k),velyplus(i,j,k),velzplus(i,j,k),epsplus(i,j,k),pressplus(i,j,k),&
+                Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,w_lorentzplus(i,j,k),&
                 gxxr,gxyr,gxzr,gyyr,gyzr,gzzr, &
                 uxxr,uxyr,uxzr,uyyr,uyzr,uzzr,avg_detr, &
-                Bvecxplus(i,j,k), Bvecyplus(i,j,k),Bveczplus(i,j,k),b2plus,&
                 epsnegative,GRHydro_C2P_failed(i,j,k))
         end do
      end do