Calculate sqrt's and divisions once, store the result and reused later.

These are optimisations that are not (in some cases) done by the
compiler except at the highest optimization levels. Leads to a 8.6%
speed improvement on my laptop. May not make any difference on some
architectures and compilers.


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

6 files changed, 109 insertions, 86 deletions

diff --git a/src/Utils.F90 b/src/Utils.F90
index 935a607..446b39d 100644
--- a/src/Utils.F90
+++ b/src/Utils.F90
@@ -94,14 +94,15 @@ subroutine UpperMetric(uxx, uxy, uxz, uyy, uyz, uzz, &
   implicit none
   
   CCTK_REAL :: uxx, uxy, uxz, uyy, uyz, uzz, det, &
-       gxx, gxy, gxz, gyy, gyz, gzz
+       gxx, gxy, gxz, gyy, gyz, gzz, invdet
   
-  uxx=(-gyz**2 + gyy*gzz)/det
-  uxy=(gxz*gyz - gxy*gzz)/det
-  uyy=(-gxz**2 + gxx*gzz)/det
-  uxz=(-gxz*gyy + gxy*gyz)/det
-  uyz=(gxy*gxz - gxx*gyz)/det
-  uzz=(-gxy**2 + gxx*gyy)/det
+  invdet = 1.d0 / det
+  uxx=(-gyz**2 + gyy*gzz)*invdet
+  uxy=(gxz*gyz - gxy*gzz)*invdet
+  uyy=(-gxz**2 + gxx*gzz)*invdet
+  uxz=(-gxz*gyy + gxy*gyz)*invdet
+  uyz=(gxy*gxz - gxx*gyz)*invdet
+  uzz=(-gxy**2 + gxx*gyy)*invdet
   
   return
   
diff --git a/src/Whisky_Con2Prim.F90 b/src/Whisky_Con2Prim.F90
index e68c36c..6361b70 100644
--- a/src/Whisky_Con2Prim.F90
+++ b/src/Whisky_Con2Prim.F90
@@ -873,7 +873,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
   CCTK_REAL s2, f, df, vlowx, vlowy, vlowz
   CCTK_INT count, handle, whisky_reflevel, whisky_C2P_failed
   CCTK_REAL udens, usx, usy, usz, rhoold, rhonew, epsold, epsnew, &
-       enthalpy, denthalpy
+       enthalpy, denthalpy, sqrtdet, invsqrtdet, invfac
   character(len=200) warnline
 
 #ifdef _EOS_BASE_INC_
@@ -884,10 +884,12 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
 
 !!$  Undensitize the variables 
 
-  udens = dens /sqrt(det)
-  usx = sx /sqrt(det)
-  usy = sy /sqrt(det)
-  usz = sz /sqrt(det)
+  sqrtdet = sqrt(det)
+  invsqrtdet = 1.d0/sqrtdet
+  udens = dens*invsqrtdet
+  usx = sx*invsqrtdet
+  usy = sy*invsqrtdet
+  usz = sz*invsqrtdet
   s2 = usx*usx*uxx + usy*usy*uyy + usz*usz*uzz + 2.d0*usx*usy*uxy + &
        2.d0*usx*usz*uxz + 2.d0*usy*usz*uyz
 
@@ -941,7 +943,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
       ! for safety, let's set the point to atmosphere
       rhonew = whisky_rho_min
       udens = rhonew
-      dens = sqrt(det) * rhonew
+      dens = sqrtdet * rhonew
       press = EOS_Pressure(handle, rhonew, 1.d0)
       epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
       sx = 0.d0
@@ -975,7 +977,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
     w_lorentz = sqrt(1.d0 + s2 / ( (udens*enthalpy)**2 ))
 
     press = EOS_Pressure(handle, rhonew, epsilon)
-    tau = sqrt(det) * ((rho * enthalpy) * w_lorentz**2 - press) - dens
+    tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens
     
     f = rhonew * w_lorentz - udens
 
@@ -987,7 +989,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
     rhonew = whisky_rho_min
     udens = rhonew
     ! before 2009/10/07 dens was not reset and was used with its negative value below; this has since been changed, but the change produces significant changes in the results
-    dens = sqrt(det) * rhonew
+    dens = sqrtdet * rhonew
     press = EOS_Pressure(handle, rhonew, 1.d0)
     epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
     ! w_lorentz=1, so the expression for utau reduces to:
@@ -1019,7 +1021,7 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
 
   press = EOS_Pressure(handle, rhonew, epsilon)
   epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
-  tau = sqrt(det) * ((rho * enthalpy) * w_lorentz**2 - press) - dens
+  tau = sqrtdet * ((rho * enthalpy) * w_lorentz**2 - press) - dens
 
   if (epsilon .lt. 0.0d0) then
     whisky_C2P_failed = 1
@@ -1036,11 +1038,11 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
     !$OMP END CRITICAL
     
     rho = whisky_rho_min
-    dens = sqrt(det) * rho
+    dens = sqrtdet * rho
     press = EOS_Pressure(handle, rhonew, 1.d0)
     epsilon = EOS_SpecificIntEnergy(handle, rhonew, press)
     ! w_lorentz=1, so the expression for tau reduces to:
-    tau  = sqrt(det) * (rho+rho*epsilon) - dens
+    tau  = sqrtdet * (rho+rho*epsilon) - dens
     usx = 0.d0
     usy = 0.d0
     usz = 0.d0
@@ -1048,9 +1050,10 @@ subroutine Con2Prim_ptPolytype(handle, dens, sx, sy, sz, tau, rho, &
 
   end if
 
-  vlowx = usx / ( (rho + rho*epsilon + press) * w_lorentz**2)
-  vlowy = usy / ( (rho + rho*epsilon + press) * w_lorentz**2)
-  vlowz = usz / ( (rho + rho*epsilon + press) * w_lorentz**2)
+  invfac = 1.d0 / ( (rho + rho*epsilon + press) * w_lorentz**2)
+  vlowx = usx * invfac
+  vlowy = usy * invfac
+  vlowz = usz * invfac
   velx = uxx * vlowx + uxy * vlowy + uxz * vlowz
   vely = uxy * vlowx + uyy * vlowy + uyz * vlowz
   velz = uxz * vlowx + uyz * vlowy + uzz * vlowz
diff --git a/src/Whisky_Eigenproblem_Marquina.F90 b/src/Whisky_Eigenproblem_Marquina.F90
index a40d8ee..34857e4 100644
--- a/src/Whisky_Eigenproblem_Marquina.F90
+++ b/src/Whisky_Eigenproblem_Marquina.F90
@@ -86,6 +86,9 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   CCTK_REAL axpr,axmr,vxpr,vxmr,xsir,dltr
   CCTK_REAL cxx,cxy,cxz,gam,vxa,vxb
   CCTK_INT handle
+  CCTK_REAL invrhol, betainvalp, invfacl
+  CCTK_REAL invrhor, invfacr
+  CCTK_REAL invulpfac, invulmfac, invurpfac, invurmfac
 
 !!$  Warning, warning. Nasty hack follows
 
@@ -104,10 +107,11 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   pressl = EOS_Pressure(handle,rhol,epsl)
   dpdrhol = EOS_DPressByDRho(handle,rhol,epsl)
   dpdepsl = EOS_DPressByDEps(handle,rhol,epsl)
-  cs2l = (dpdrhol + pressl * dpdepsl / (rhol**2))/ &
-       (1.0d0 + epsl + pressl/rhol)
+  invrhol = one / rhol
+  cs2l = (dpdrhol + pressl * dpdepsl * invrhol**2)/ &
+       (1.0d0 + epsl + pressl*invrhol)
 
-  enthalpyl = one + epsl + pressl / rhol
+  enthalpyl = one + epsl + pressl * invrhol
 
   vlowxl = gxx*velxl + gxy*velyl + gxz*velzl
   vlowyl = gxy*velxl + gyy*velyl + gyz*velzl
@@ -127,13 +131,15 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
 
 !!$  EIGENVALUES
 
-  lam1l = velxl - beta/alp
+  betainvalp = beta / alp
+  lam1l = velxl - betainvalp
   lam2l = lam1l
   lam3l = lam1l
+  invfacl = one/(one-v2l*cs2l)
   lamp_nobetal = (velxl*(one-cs2l) + sqrt(cs2l*(one-v2l)* &
-       (u*(one-v2l*cs2l) - velxl**2*(one-cs2l))))/(one-v2l*cs2l)
+       (u*(one-v2l*cs2l) - velxl**2*(one-cs2l))))*invfacl
   lamm_nobetal = (velxl*(one-cs2l) - sqrt(cs2l*(one-v2l)* &
-       (u*(one-v2l*cs2l) - velxl**2*(one-cs2l))))/(one-v2l*cs2l)
+       (u*(one-v2l*cs2l) - velxl**2*(one-cs2l))))*invfacl
 
 !!$ BEGIN: Check for NaN value (2nd check)
   if (lamp_nobetal .ne. lamp_nobetal) then
@@ -142,8 +148,8 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   endif
 !!$ END: Check for NaN value (2nd check)
 
-  lampl = lamp_nobetal - beta/alp
-  lamml = lamm_nobetal - beta/alp
+  lampl = lamp_nobetal - betainvalp
+  lamml = lamm_nobetal - betainvalp
 
 !!$  RIGHT
 
@@ -152,10 +158,11 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   pressr = EOS_Pressure(handle,rhor,epsr)
   dpdrhor = EOS_DPressByDRho(handle,rhor,epsr)
   dpdepsr = EOS_DPressByDEps(handle,rhor,epsr)
-  cs2r = (dpdrhor + pressr * dpdepsr / (rhor**2))/ &
-       (1.0d0 + epsr + pressr/rhor)
+  invrhor = one / rhor
+  cs2r = (dpdrhor + pressr * dpdepsr * invrhor**2)/ &
+       (1.0d0 + epsr + pressr*invrhor)
 
-  enthalpyr = one + epsr + pressr / rhor
+  enthalpyr = one + epsr + pressr * invrhor
 
   vlowxr = gxx*velxr + gxy*velyr + gxz*velzr
   vlowyr = gxy*velxr + gyy*velyr + gyz*velzr
@@ -175,13 +182,14 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
 
 !!$  EIGENVALUES
 
-  lam1r = velxr - beta/alp
+  lam1r = velxr - betainvalp
   lam2r = lam1r
   lam3r = lam1r
+  invfacr = one/(one-v2r*cs2r)
   lamp_nobetar = (velxr*(one-cs2r) + sqrt(cs2r*(one-v2r)* &
-       (u*(one-v2r*cs2r) - velxr**2*(one-cs2r))))/(one-v2r*cs2r)
+       (u*(one-v2r*cs2r) - velxr**2*(one-cs2r))))*invfacr
   lamm_nobetar = (velxr*(one-cs2r) - sqrt(cs2r*(one-v2r)* &
-       (u*(one-v2r*cs2r) - velxr**2*(one-cs2r))))/(one-v2r*cs2r)
+       (u*(one-v2r*cs2r) - velxr**2*(one-cs2r))))*invfacr
 
 !!$ BEGIN: Check for NaN value (4th check)
   if (lamp_nobetar .ne. lamp_nobetar) then
@@ -190,8 +198,8 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
   endif
 !!$ END: Check for NaN value (4th check)
 
-  lampr = lamp_nobetar - beta/alp
-  lammr = lamm_nobetar - beta/alp
+  lampr = lamp_nobetar - betainvalp
+  lammr = lamm_nobetar - betainvalp
 
 !!$  FINAL
 
@@ -217,10 +225,12 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
 
 !!$  Compute some auxiliary quantities
 
-  axpl = (u - velxl*velxl)/(u - velxl*lamp_nobetal)
-  axml = (u - velxl*velxl)/(u - velxl*lamm_nobetal)
-  vxpl = (velxl - lamp_nobetal)/(u - velxl * lamp_nobetal)
-  vxml = (velxl - lamm_nobetal)/(u - velxl * lamm_nobetal)
+  invulpfac = one / (u - velxl*lamp_nobetal)
+  invulmfac = one / (u - velxl*lamm_nobetal)
+  axpl = (u - velxl*velxl)*invulpfac
+  axml = (u - velxl*velxl)*invulmfac
+  vxpl = (velxl - lamp_nobetal)*invulpfac
+  vxml = (velxl - lamm_nobetal)*invulmfac
 
 !!$  Calculate associated right eigenvectors
 
@@ -270,10 +280,12 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
 
 !!$  Compute some auxiliary quantities
 
-  axpr = (u - velxr*velxr)/(u - velxr*lamp_nobetar)
-  axmr = (u - velxr*velxr)/(u - velxr*lamm_nobetar)
-  vxpr = (velxr - lamp_nobetar)/(u - velxr * lamp_nobetar)
-  vxmr = (velxr - lamm_nobetar)/(u - velxr * lamm_nobetar)
+  invurpfac = one / (u - velxr*lamp_nobetar)
+  invurmfac = one / (u - velxr*lamm_nobetar)
+  axpr = (u - velxr*velxr)*invurpfac
+  axmr = (u - velxr*velxr)*invurmfac
+  vxpr = (velxr - lamp_nobetar)*invurpfac
+  vxmr = (velxr - lamm_nobetar)*invurmfac
 
 !!$  Calculate associated right eigenvectors
 
@@ -560,9 +572,9 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
 !!$  Get lower left triangle to be all zeros
     do i=1,5
 !!$  First, make diagonal element 1
-      tmp1 = paug(i,i)
+      tmp1 = one / paug(i,i)
       do j=i,6
-        paug(i,j) = paug(i,j)/tmp1
+        paug(i,j) = paug(i,j)*tmp1
       enddo
 !!$  Now, get rid of everything below that diagonal
       do j=i+1,5
@@ -700,9 +712,9 @@ subroutine eigenproblem_marquina(handle,rhor,velxr,velyr,&
 
     do i=1,5
 !!$  First, make diagonal element 1
-      tmp1 = paug(i,i)
+      tmp1 = one / paug(i,i)
       do j=i,6
-        paug(i,j) = paug(i,j)/tmp1
+        paug(i,j) = paug(i,j)*tmp1
       enddo
 !!$  Now, get rid of everything below that diagonal
       do j=i+1,5
diff --git a/src/Whisky_Flux.F90 b/src/Whisky_Flux.F90
index bd59379..62866d2 100644
--- a/src/Whisky_Flux.F90
+++ b/src/Whisky_Flux.F90
@@ -83,20 +83,23 @@ subroutine num_x_flux(dens,sx,sy,sz,tau,&
      densf,sxf,syf,szf,tauf,vel,press,det,alp,beta)
 
   implicit none
-  
+
   CCTK_REAL :: dens,sx,sy,sz,tau,vel,det
   CCTK_REAL :: densf,sxf,syf,szf,tauf
   CCTK_REAL :: alp,beta,press 
-  
-  densf = dens * ( vel  - beta / alp )
-  
-  sxf = sx * ( vel  - beta/alp) + sqrt(det)*press
-               
-  syf = sy * ( vel  - beta/alp )
-              
-  szf = sz * ( vel - beta/alp )
-               
-  tauf = tau * ( vel - beta/alp ) + sqrt(det)*press*vel
+  CCTK_REAL :: velmbetainvalp,sqrtdetpress
+
+  velmbetainvalp = vel - beta / alp
+  sqrtdetpress = sqrt(det)*press
+  densf = dens * ( velmbetainvalp )
+
+  sxf = sx * velmbetainvalp + sqrtdetpress
+
+  syf = sy * velmbetainvalp
+
+  szf = sz * velmbetainvalp
+
+  tauf = tau * velmbetainvalp + sqrtdetpress*vel
  
 end subroutine num_x_flux
 
diff --git a/src/Whisky_Prim2Con.F90 b/src/Whisky_Prim2Con.F90
index 74d5fcf..3669516 100644
--- a/src/Whisky_Prim2Con.F90
+++ b/src/Whisky_Prim2Con.F90
@@ -326,7 +326,7 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, &
   
   CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det
   CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,&
-       deps, dpress, w, vlowx, vlowy, vlowz   
+       deps, dpress, w, vlowx, vlowy, vlowz, sqrtdet   
   CCTK_INT :: handle
 
 #ifdef _EOS_BASE_INC_
@@ -345,11 +345,12 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, &
   vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz
   vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz
 
-  ddens = sqrt(det) * drho * w 
-  dsx = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowx
-  dsy = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowy
-  dsz = sqrt(det) * (drho*(1+deps)+dpress)*w*w * vlowz
-  dtau = sqrt(det) * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens 
+  sqrtdet = sqrt(det)
+  ddens = sqrtdet * drho * w 
+  dsx = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowx
+  dsy = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowy
+  dsz = sqrtdet * (drho*(1+deps)+dpress)*w*w * vlowz
+  dtau = sqrtdet * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens 
 
 end subroutine prim2conpolytype
 
diff --git a/src/Whisky_Source.F90 b/src/Whisky_Source.F90
index f96f317..e6400cb 100644
--- a/src/Whisky_Source.F90
+++ b/src/Whisky_Source.F90
@@ -71,6 +71,7 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
   CCTK_REAL :: sumTK
   CCTK_REAL :: halfshiftdgx, halfshiftdgy, halfshiftdgz
   CCTK_REAL :: halfTdgx, halfTdgy, halfTdgz
+  CCTK_REAL :: invalp, invalp2
 
   logical, allocatable, dimension (:,:,:) :: force_spatial_second_order
 
@@ -140,7 +141,7 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
   !$OMP dz_gxx, dz_gxy, dz_gxz, dz_gyy, dz_gyz, dz_gzz,&
   !$OMP dx_psi4,dy_psi4,dz_psi4,shiftshiftk,shiftkx,shiftky,shiftkz,&
   !$OMP sumTK,halfshiftdgx,halfshiftdgy,halfshiftdgz,&
-  !$OMP halfTdgx,halfTdgy,halfTdgz)
+  !$OMP halfTdgx,halfTdgy,halfTdgz,invalp,invalp2)
   do k=1 + whisky_stencil,nz - whisky_stencil
     do j=1 + whisky_stencil,ny - whisky_stencil
       do i=1 + whisky_stencil,nx - whisky_stencil
@@ -243,9 +244,11 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
           
         endif
 
-        velxshift = velx(i,j,k) - shiftx/alp(i,j,k)
-        velyshift = vely(i,j,k) - shifty/alp(i,j,k)
-        velzshift = velz(i,j,k) - shiftz/alp(i,j,k)
+        invalp = 1.0d0 / alp(i,j,k)
+        invalp2 = invalp**2
+        velxshift = velx(i,j,k) - shiftx*invalp
+        velyshift = vely(i,j,k) - shifty*invalp
+        velzshift = velz(i,j,k) - shiftz*invalp
         vlowx = velx(i,j,k)*localgxx + vely(i,j,k)*localgxy +&
              velz(i,j,k)*localgxz
         vlowy = velx(i,j,k)*localgxy + vely(i,j,k)*localgyy +&
@@ -255,25 +258,25 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
 
 !!$        For a change, these are T^{ij}
 
-        t00 = (rhoenthalpyW2 - press(i,j,k))/(alp(i,j,k)**2)
+        t00 = (rhoenthalpyW2 - press(i,j,k))*invalp2
         t0x = rhoenthalpyW2*velxshift/alp(i,j,k) +&
-             press(i,j,k)*shiftx/(alp(i,j,k)**2)
+             press(i,j,k)*shiftx*invalp2
         t0y = rhoenthalpyW2*velyshift/alp(i,j,k) +&
-             press(i,j,k)*shifty/(alp(i,j,k)**2)
+             press(i,j,k)*shifty*invalp2
         t0z = rhoenthalpyW2*velzshift/alp(i,j,k) +&
-             press(i,j,k)*shiftz/(alp(i,j,k)**2)
+             press(i,j,k)*shiftz*invalp2
         txx = rhoenthalpyW2*velxshift*velxshift +&
-             press(i,j,k)*(uxx - shiftx*shiftx/(alp(i,j,k)**2))
+             press(i,j,k)*(uxx - shiftx*shiftx*invalp2)
         txy = rhoenthalpyW2*velxshift*velyshift +&
-             press(i,j,k)*(uxy - shiftx*shifty/(alp(i,j,k)**2))
+             press(i,j,k)*(uxy - shiftx*shifty*invalp2)
         txz = rhoenthalpyW2*velxshift*velzshift +&
-             press(i,j,k)*(uxz - shiftx*shiftz/(alp(i,j,k)**2))
+             press(i,j,k)*(uxz - shiftx*shiftz*invalp2)
         tyy = rhoenthalpyW2*velyshift*velyshift +&
-             press(i,j,k)*(uyy - shifty*shifty/(alp(i,j,k)**2))
+             press(i,j,k)*(uyy - shifty*shifty*invalp2)
         tyz = rhoenthalpyW2*velyshift*velzshift +&
-             press(i,j,k)*(uyz - shifty*shiftz/(alp(i,j,k)**2))
+             press(i,j,k)*(uyz - shifty*shiftz*invalp2)
         tzz = rhoenthalpyW2*velzshift*velzshift +&
-             press(i,j,k)*(uzz - shiftz*shiftz/(alp(i,j,k)**2))
+             press(i,j,k)*(uzz - shiftz*shiftz*invalp2)
 
 !!$        Derivatives of the lapse, metric and shift
 
@@ -457,24 +460,24 @@ subroutine SourceTerms(CCTK_ARGUMENTS)
              t0y*(shiftx*dx_gxy + shifty*dx_gyy + shiftz*dx_gyz) +&
              t0z*(shiftx*dx_gxz + shifty*dx_gyz + shiftz*dx_gzz) +&
              halfTdgx + rhoenthalpyW2*&
-             (vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz)/&
-             alp(i,j,k)
+             (vlowx*dx_betax + vlowy*dx_betay + vlowz*dx_betaz)*&
+             invalp
         sy_source = t00*&
              (halfshiftdgy - alp(i,j,k)*dy_alp) +&
              t0x*(shiftx*dy_gxx + shifty*dy_gxy + shiftz*dy_gxz) +&
              t0y*(shiftx*dy_gxy + shifty*dy_gyy + shiftz*dy_gyz) +&
              t0z*(shiftx*dy_gxz + shifty*dy_gyz + shiftz*dy_gzz) +&
              halfTdgy + rhoenthalpyW2*&
-             (vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz)/&
-             alp(i,j,k)
+             (vlowx*dy_betax + vlowy*dy_betay + vlowz*dy_betaz)*&
+             invalp
         sz_source = t00*&
              (halfshiftdgz - alp(i,j,k)*dz_alp) +&
              t0x*(shiftx*dz_gxx + shifty*dz_gxy + shiftz*dz_gxz) +&
              t0y*(shiftx*dz_gxy + shifty*dz_gyy + shiftz*dz_gyz) +&
              t0z*(shiftx*dz_gxz + shifty*dz_gyz + shiftz*dz_gzz) +&
              halfTdgz + rhoenthalpyW2*&
-             (vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz)/&
-             alp(i,j,k)
+             (vlowx*dz_betax + vlowy*dz_betay + vlowz*dz_betaz)*&
+             invalp
 
         densrhs(i,j,k) = 0.d0
         srhs(i,j,k,1)  = alp(i,j,k)*sqrtdet*sx_source