diff options
Diffstat (limited to 'src/GRHydro_Prim2Con.F90')
| -rw-r--r-- | src/GRHydro_Prim2Con.F90 | 133 |
1 files changed, 61 insertions, 72 deletions
diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index 680ecfc..20a6c9d 100644 --- a/src/GRHydro_Prim2Con.F90 +++ b/src/GRHydro_Prim2Con.F90 @@ -51,8 +51,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) integer :: i, j, k integer :: keytemp - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr CCTK_REAL :: xtemp(1) ! save memory when MP is not used @@ -79,7 +79,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) end if if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr,& + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_sdetl, avg_sdetr,& !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 @@ -99,12 +99,12 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2con(GRHydro_eos_handle, g11l,g12l,g13l,g22l,& g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,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),& @@ -112,7 +112,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) call prim2con(GRHydro_eos_handle, g11r,g12r,g13r,& g22r,g23r,g33r, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, 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),& @@ -125,7 +125,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) else if(reconstruct_temper.ne.0) then keytemp = 1 - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_sdetl, avg_sdetr, xtemp,& !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 @@ -145,15 +145,15 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2con_hot(GRHydro_eos_handle, keytemp, GRHydro_reflevel,& cctk_iteration,i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),& r(i,j,k),& g11l,g12l,g13l,g22l,& g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,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),& @@ -165,7 +165,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) r(i,j,k),& g11r,g12r,g13r,& g22r,g23r,g33r,& - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, 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),& @@ -178,7 +178,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) !$OMP END PARALLEL DO else keytemp = 0 - !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& + !$OMP PARALLEL DO PRIVATE(i, j, k, avg_sdetl, avg_sdetr, xtemp,& !$OMP g11l,g12l,g13l,g22l,g23l,g33l, & !$OMP g11r,g12r,g13r,g22r,g23r,g33r) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 @@ -198,8 +198,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) xtemp(1) = 0.5d0*(temperature(i,j,k) + & temperature(i-xoffset,j-yoffset,k-zoffset)) @@ -208,7 +208,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) r(i,j,k),& g11l,g12l,g13l,g22l,& g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,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),& @@ -222,7 +222,7 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) r(i,j,k),& g11r,g12r,g13r,& g22r,g23r,g33r,& - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, 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),& @@ -254,14 +254,14 @@ end subroutine primitive2conservative @@*/ -subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & +subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, sdetg, ddens, & dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdetg CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,& deps, dpress, w, vlowx, vlowy, vlowz CCTK_INT :: handle @@ -284,17 +284,17 @@ subroutine prim2con(handle, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & 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 + ddens = sdetg * drho * w + dsx = sdetg * (drho*(1+deps)+dpress)*w*w * vlowx + dsy = sdetg * (drho*(1+deps)+dpress)*w*w * vlowy + dsz = sdetg * (drho*(1+deps)+dpress)*w*w * vlowz + dtau = sdetg * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens end subroutine prim2con subroutine prim2con_hot(handle, keytemp, GRHydro_reflevel, cctk_iteration, ii, jj, kk, & - x, y, z, r, gxx, gxy, gxz, gyy, gyz, gzz, det, ddens, & + x, y, z, r, gxx, gxy, gxz, gyy, gyz, gzz, sdetg, ddens, & dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w, & temp,ye) @@ -303,12 +303,12 @@ subroutine prim2con_hot(handle, keytemp, GRHydro_reflevel, cctk_iteration, ii, j DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdetg CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho(1), dvelx, dvely, dvelz,& deps(1), dpress(1), w, vlowx, vlowy, vlowz CCTK_REAL :: temp(1),ye(1), x, y, z, r CCTK_INT :: handle, GRHydro_reflevel, cctk_iteration, ii, jj, kk - CCTK_REAL :: sdet, h + CCTK_REAL :: h character(len=512) warnline ! begin EOS Omni vars @@ -439,13 +439,12 @@ subroutine prim2con_hot(handle, keytemp, GRHydro_reflevel, cctk_iteration, ii, j vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz - sdet = sqrt(det) h = drho(1)*(1.0d0+deps(1))+dpress(1) - ddens = sdet * drho(1) * w - dsx = sdet * h*w*w * vlowx - dsy = sdet * h*w*w * vlowy - dsz = sdet * h*w*w * vlowz - dtau = sdet * (h*w*w - dpress(1)) - ddens + ddens = sdetg * drho(1) * w + dsx = sdetg * h*w*w * vlowx + dsy = sdetg * h*w*w * vlowy + dsz = sdetg * h*w*w * vlowz + dtau = sdetg * (h*w*w - dpress(1)) - ddens end subroutine prim2con_hot @@ -483,8 +482,8 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) CCTK_INT :: i, j, k CCTK_INT :: keytemp - CCTK_REAL :: det + character(len=512) :: warnline ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates CCTK_REAL, DIMENSION(:,:,:), POINTER :: g11, g12, g13, g22, g23, g33 @@ -509,41 +508,36 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) end if if(evolve_temper.ne.1) then - !$OMP PARALLEL DO PRIVATE(i, j, k, det) + !$OMP PARALLEL DO PRIVATE(i, j, k) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k), \ - g22(i,j,k),g23(i,j,k),g33(i,j,k)) - call prim2con(GRHydro_eos_handle,g11(i,j,k),& g12(i,j,k),g13(i,j,k),& g22(i,j,k),g23(i,j,k),g33(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdetg(i,j,k), dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) + end do end do end do !$OMP END PARALLEL DO else keytemp = 0 - !$OMP PARALLEL DO PRIVATE(i, j, k, det) + !$OMP PARALLEL DO PRIVATE(i, j, k) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k), \ - g22(i,j,k),g23(i,j,k),g33(i,j,k)) - call prim2con_hot(GRHydro_eos_handle,keytemp,& GRHydro_reflevel,cctk_iteration,& i,j,k,& x(i,j,k),y(i,j,k),z(i,j,k),r(i,j,k),g11(i,j,k),& g12(i,j,k),g13(i,j,k),& g22(i,j,k),g23(i,j,k),g33(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdetg(i,j,k), dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k), & temperature(i,j,k),y_e(i,j,k)) @@ -599,8 +593,8 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -625,8 +619,8 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) vup => vel end if - !$OMP PARALLEL DO PRIVATE(i, j, k, g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - !$OMP g11r,g12r,g13r,g22r,g23r,g33r,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + !$OMP g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr) do k = GRHydro_stencil,cctk_lsh(3)-GRHydro_stencil+1 do j = GRHydro_stencil,cctk_lsh(2)-GRHydro_stencil+1 do i = GRHydro_stencil,cctk_lsh(1)-GRHydro_stencil+1 @@ -644,19 +638,19 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) call prim2conpolytype(GRHydro_eos_handle, g11l,g12l,g13l,& g22l,g23l,g33l, & - avg_detl,densminus(i,j,k),sxminus(i,j,k),& + avg_sdetl,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)) call prim2conpolytype(GRHydro_eos_handle, g11r,g12r,g13r,& g22r,g23r,g33r, & - avg_detr, densplus(i,j,k),sxplus(i,j,k),& + avg_sdetr, 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),& @@ -689,14 +683,14 @@ end subroutine Prim2ConservativePolytype @@*/ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, & - gzz, det, ddens, & + gzz, sdetg, ddens, & dsx, dsy, dsz, dtau , drho, dvelx, dvely, dvelz, deps, dpress, w) implicit none DECLARE_CCTK_PARAMETERS - CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, det + CCTK_REAL :: gxx, gxy, gxz, gyy, gyz, gzz, sdetg CCTK_REAL :: ddens, dsx, dsy, dsz, dtau, drho, dvelx, dvely, dvelz,& deps, dpress, w_tmp, w, vlowx, vlowy, vlowz, sqrtdet CCTK_INT :: handle @@ -737,12 +731,11 @@ subroutine prim2conpolytype(handle, gxx, gxy, gxz, gyy, gyz, & vlowy = gxy*dvelx + gyy*dvely + gyz*dvelz vlowz = gxz*dvelx + gyz*dvely + gzz*dvelz - 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 + ddens = sdetg * drho * w + dsx = sdetg * (drho*(1+deps)+dpress)*w*w * vlowx + dsy = sdetg * (drho*(1+deps)+dpress)*w*w * vlowy + dsz = sdetg * (drho*(1+deps)+dpress)*w*w * vlowz + dtau = sdetg * ((drho*(1+deps)+dpress)*w*w - dpress) - ddens end subroutine prim2conpolytype @@ -778,7 +771,6 @@ subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS CCTK_INT :: i, j, k - CCTK_REAL :: det ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -803,17 +795,14 @@ subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS) vup => vel end if - !$OMP PARALLEL DO PRIVATE(i, j, k, det) + !$OMP PARALLEL DO PRIVATE(i, j, k) do k = 1,cctk_lsh(3) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - det = SPATIAL_DETERMINANT(g11(i,j,k),g12(i,j,k),g13(i,j,k),\ - g22(i,j,k),g23(i,j,k),g33(i,j,k)) - call prim2conpolytype(GRHydro_eos_handle,g11(i,j,k),g12(i,j,k),& g13(i,j,k),g22(i,j,k),g23(i,j,k),g33(i,j,k),& - det, dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& + sdetg(i,j,k), dens(i,j,k),sx(i,j,k),sy(i,j,k),sz(i,j,k),& tau(i,j,k),rho(i,j,k),velx(i,j,k),vely(i,j,k),velz(i,j,k),& eps(i,j,k),press(i,j,k),w_lorentz(i,j,k)) @@ -864,8 +853,8 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_detl,& - g11r,g12r,g13r,g22r,g23r,g33r,avg_detr + CCTK_REAL :: g11l,g12l,g13l,g22l,g23l,g33l,avg_sdetl,& + g11r,g12r,g13r,g22r,g23r,g33r,avg_sdetr ! save memory when MP is not used CCTK_INT :: GRHydro_UseGeneralCoordinates @@ -904,10 +893,10 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) g23r = 0.5d0 * (g23(i,j,k) + g23(i+xoffset,j+yoffset,k+zoffset)) g33r = 0.5d0 * (g33(i,j,k) + g33(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l) - avg_detr = SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r) + avg_sdetl = sqrt(SPATIAL_DETERMINANT(g11l,g12l,g13l,g22l, g23l,g33l)) + avg_sdetr = sqrt(SPATIAL_DETERMINANT(g11r,g12r,g13r,g22r, g23r,g33r)) cons_tracerplus(i,j,k,:) = tracerplus(i,j,k,:) * & - sqrt(avg_detr) * rhoplus(i,j,k) / & + avg_sdetr * rhoplus(i,j,k) / & sqrt(1.d0 - & (g11r * velxplus(i,j,k)**2 + & g22r * velyplus(i,j,k)**2 + & @@ -916,7 +905,7 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) g13r * velxplus(i,j,k) * velzplus(i,j,k) + & g23r * velyplus(i,j,k) * velzplus(i,j,k) ) ) ) cons_tracerminus(i,j,k,:) = tracerminus(i,j,k,:) * & - sqrt(avg_detl) * rhominus(i,j,k) / & + avg_sdetl * rhominus(i,j,k) / & sqrt(1.d0 - & (g11l * velxminus(i,j,k)**2 + & g22l * velyminus(i,j,k)**2 + & |