diff options
Diffstat (limited to 'src/GRHydro_Prim2Con.F90')
| -rw-r--r-- | src/GRHydro_Prim2Con.F90 | 224 |
1 files changed, 112 insertions, 112 deletions
diff --git a/src/GRHydro_Prim2Con.F90 b/src/GRHydro_Prim2Con.F90 index fbe118b..2bd2d29 100644 --- a/src/GRHydro_Prim2Con.F90 +++ b/src/GRHydro_Prim2Con.F90 @@ -13,9 +13,9 @@ #include "cctk_Functions.h" #include "GRHydro_Macros.h" -#define velx(i,j,k) vel(i,j,k,1) -#define vely(i,j,k) vel(i,j,k,2) -#define velz(i,j,k) vel(i,j,k,3) +#define velx(i,j,k) lvel(i,j,k,1) +#define vely(i,j,k) lvel(i,j,k,2) +#define velz(i,j,k) lvel(i,j,k,3) #define sx(i,j,k) scon(i,j,k,1) #define sy(i,j,k) scon(i,j,k,2) #define sz(i,j,k) scon(i,j,k,3) @@ -44,44 +44,44 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& + gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr CCTK_REAL :: xtemp(1) if(evolve_temper.ne.1) then !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr,& - !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + !$OMP gaal,gabl,gacl,gbbl,gbcl,gccl, & + !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr) 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 - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl) + avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr) - call prim2con(GRHydro_eos_handle, gxxl,gxyl,gxzl,gyyl,& - gyzl,gzzl, & + call prim2con(GRHydro_eos_handle, gaal,gabl,gacl,gbbl,& + gbcl,gccl, & avg_detl,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 prim2con(GRHydro_eos_handle, gxxr,gxyr,gxzr,& - gyyr,gyzr,gzzr, & + call prim2con(GRHydro_eos_handle, gaar,gabr,gacr,& + gbbr,gbcr,gccr, & avg_detr, 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),& @@ -94,27 +94,27 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) !$OMP END PARALLEL DO else !$OMP PARALLEL DO PRIVATE(i, j, k, avg_detl, avg_detr, xtemp,& - !$OMP gxxl,gxyl,gxzl,gyyl,gyzl,gzzl, & - !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr) + !$OMP gaal,gabl,gacl,gbbl,gbcl,gccl, & + !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr) 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 - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) + gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl) + avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr) ! we do not have plus/minus vars for temperature since ! eps is reconstructed. Hence, we do not update the @@ -124,8 +124,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) xtemp(1) = 0.5d0*(temperature(i,j,k) + & temperature(i-xoffset,j-yoffset,k-zoffset)) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel,& - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxl,gxyl,gxzl,gyyl,& - gyzl,gzzl, & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaal,gabl,gacl,gbbl,& + gbcl,gccl, & avg_detl,densminus(i,j,k),sxminus(i,j,k),& syminus(i,j,k),szminus(i,j,k),tauminus(i,j,k),& rhominus(i,j,k), & @@ -141,8 +141,8 @@ subroutine primitive2conservative(CCTK_ARGUMENTS) xtemp(1) = 0.5d0*(temperature(i,j,k) + & temperature(i+xoffset,j+yoffset,k+zoffset)) call prim2con_hot(GRHydro_eos_handle, GRHydro_reflevel, & - i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gxxr,gxyr,gxzr,& - gyyr,gyzr,gzzr, & + i,j,k,x(i,j,k),y(i,j,k),z(i,j,k),gaar,gabr,gacr,& + gbbr,gbcr,gccr, & avg_detr, 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),& @@ -353,12 +353,12 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) - call prim2con(GRHydro_eos_handle,gxx(i,j,k),& - gxy(i,j,k),gxz(i,j,k),& - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + call prim2con(GRHydro_eos_handle,gaa(i,j,k),& + gab(i,j,k),gac(i,j,k),& + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k),& det, 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)) @@ -372,13 +372,13 @@ subroutine Primitive2ConservativeCells(CCTK_ARGUMENTS) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k), \ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k), \ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) call prim2con_hot(GRHydro_eos_handle,GRHydro_reflevel,i,j,k,& - x(i,j,k),y(i,j,k),z(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),& + x(i,j,k),y(i,j,k),z(i,j,k),gaa(i,j,k),& + gab(i,j,k),gac(i,j,k),& + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k),& det, 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), & @@ -429,40 +429,40 @@ subroutine Prim2ConservativePolytype(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& + gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr - !$OMP PARALLEL DO PRIVATE(i, j, k, gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - !$OMP gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr) + !$OMP PARALLEL DO PRIVATE(i, j, k, gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& + !$OMP gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr) 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 - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) - - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) - - call prim2conpolytype(GRHydro_eos_handle, gxxl,gxyl,gxzl,& - gyyl,gyzl,gzzl, & + gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl) + avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr) + + call prim2conpolytype(GRHydro_eos_handle, gaal,gabl,gacl,& + gbbl,gbcl,gccl, & avg_detl,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, gxxr,gxyr,gxzr,& - gyyr,gyzr,gzzr, & + call prim2conpolytype(GRHydro_eos_handle, gaar,gabr,gacr,& + gbbr,gbcr,gccr, & avg_detr, 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),& @@ -586,11 +586,11 @@ subroutine Primitive2ConservativePolyCells(CCTK_ARGUMENTS) do j = 1,cctk_lsh(2) do i = 1,cctk_lsh(1) - det = SPATIAL_DETERMINANT(gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),\ - gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)) + det = SPATIAL_DETERMINANT(gaa(i,j,k),gab(i,j,k),gac(i,j,k),\ + gbb(i,j,k),gbc(i,j,k),gcc(i,j,k)) - call prim2conpolytype(GRHydro_eos_handle,gxx(i,j,k),gxy(i,j,k),& - gxz(i,j,k),gyy(i,j,k),gyz(i,j,k),gzz(i,j,k),& + call prim2conpolytype(GRHydro_eos_handle,gaa(i,j,k),gab(i,j,k),& + gac(i,j,k),gbb(i,j,k),gbc(i,j,k),gcc(i,j,k),& det, 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)) @@ -637,46 +637,46 @@ subroutine Prim2ConservativeTracer(CCTK_ARGUMENTS) DECLARE_CCTK_PARAMETERS integer :: i, j, k - CCTK_REAL :: gxxl,gxyl,gxzl,gyyl,gyzl,gzzl,avg_detl,& - gxxr,gxyr,gxzr,gyyr,gyzr,gzzr,avg_detr + CCTK_REAL :: gaal,gabl,gacl,gbbl,gbcl,gccl,avg_detl,& + gaar,gabr,gacr,gbbr,gbcr,gccr,avg_detr 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 - gxxl = 0.5d0 * (gxx(i,j,k) + gxx(i-xoffset,j-yoffset,k-zoffset)) - gxyl = 0.5d0 * (gxy(i,j,k) + gxy(i-xoffset,j-yoffset,k-zoffset)) - gxzl = 0.5d0 * (gxz(i,j,k) + gxz(i-xoffset,j-yoffset,k-zoffset)) - gyyl = 0.5d0 * (gyy(i,j,k) + gyy(i-xoffset,j-yoffset,k-zoffset)) - gyzl = 0.5d0 * (gyz(i,j,k) + gyz(i-xoffset,j-yoffset,k-zoffset)) - gzzl = 0.5d0 * (gzz(i,j,k) + gzz(i-xoffset,j-yoffset,k-zoffset)) - gxxr = 0.5d0 * (gxx(i,j,k) + gxx(i+xoffset,j+yoffset,k+zoffset)) - gxyr = 0.5d0 * (gxy(i,j,k) + gxy(i+xoffset,j+yoffset,k+zoffset)) - gxzr = 0.5d0 * (gxz(i,j,k) + gxz(i+xoffset,j+yoffset,k+zoffset)) - gyyr = 0.5d0 * (gyy(i,j,k) + gyy(i+xoffset,j+yoffset,k+zoffset)) - gyzr = 0.5d0 * (gyz(i,j,k) + gyz(i+xoffset,j+yoffset,k+zoffset)) - gzzr = 0.5d0 * (gzz(i,j,k) + gzz(i+xoffset,j+yoffset,k+zoffset)) - - avg_detl = SPATIAL_DETERMINANT(gxxl,gxyl,gxzl,gyyl, gyzl,gzzl) - avg_detr = SPATIAL_DETERMINANT(gxxr,gxyr,gxzr,gyyr, gyzr,gzzr) + gaal = 0.5d0 * (gaa(i,j,k) + gaa(i-xoffset,j-yoffset,k-zoffset)) + gabl = 0.5d0 * (gab(i,j,k) + gab(i-xoffset,j-yoffset,k-zoffset)) + gacl = 0.5d0 * (gac(i,j,k) + gac(i-xoffset,j-yoffset,k-zoffset)) + gbbl = 0.5d0 * (gbb(i,j,k) + gbb(i-xoffset,j-yoffset,k-zoffset)) + gbcl = 0.5d0 * (gbc(i,j,k) + gbc(i-xoffset,j-yoffset,k-zoffset)) + gccl = 0.5d0 * (gcc(i,j,k) + gcc(i-xoffset,j-yoffset,k-zoffset)) + gaar = 0.5d0 * (gaa(i,j,k) + gaa(i+xoffset,j+yoffset,k+zoffset)) + gabr = 0.5d0 * (gab(i,j,k) + gab(i+xoffset,j+yoffset,k+zoffset)) + gacr = 0.5d0 * (gac(i,j,k) + gac(i+xoffset,j+yoffset,k+zoffset)) + gbbr = 0.5d0 * (gbb(i,j,k) + gbb(i+xoffset,j+yoffset,k+zoffset)) + gbcr = 0.5d0 * (gbc(i,j,k) + gbc(i+xoffset,j+yoffset,k+zoffset)) + gccr = 0.5d0 * (gcc(i,j,k) + gcc(i+xoffset,j+yoffset,k+zoffset)) + + avg_detl = SPATIAL_DETERMINANT(gaal,gabl,gacl,gbbl, gbcl,gccl) + avg_detr = SPATIAL_DETERMINANT(gaar,gabr,gacr,gbbr, gbcr,gccr) cons_tracerplus(i,j,k,:) = tracerplus(i,j,k,:) * & sqrt(avg_detr) * rhoplus(i,j,k) / & sqrt(1.d0 - & - (gxxr * velxplus(i,j,k)**2 + & - gyyr * velyplus(i,j,k)**2 + & - gzzr * velzplus(i,j,k)**2 + & - 2.d0 * (gxyr * velxplus(i,j,k) * velyplus(i,j,k) + & - gxzr * velxplus(i,j,k) * velzplus(i,j,k) + & - gyzr * velyplus(i,j,k) * velzplus(i,j,k) ) ) ) + (gaar * velxplus(i,j,k)**2 + & + gbbr * velyplus(i,j,k)**2 + & + gccr * velzplus(i,j,k)**2 + & + 2.d0 * (gabr * velxplus(i,j,k) * velyplus(i,j,k) + & + gacr * velxplus(i,j,k) * velzplus(i,j,k) + & + gbcr * velyplus(i,j,k) * velzplus(i,j,k) ) ) ) cons_tracerminus(i,j,k,:) = tracerminus(i,j,k,:) * & sqrt(avg_detl) * rhominus(i,j,k) / & sqrt(1.d0 - & - (gxxl * velxminus(i,j,k)**2 + & - gyyl * velyminus(i,j,k)**2 + & - gzzl * velzminus(i,j,k)**2 + & - 2.d0 * (gxyl * velxminus(i,j,k) * velyminus(i,j,k) + & - gxzl * velxminus(i,j,k) * velzminus(i,j,k) + & - gyzl * velyminus(i,j,k) * velzminus(i,j,k) ) ) ) + (gaal * velxminus(i,j,k)**2 + & + gbbl * velyminus(i,j,k)**2 + & + gccl * velzminus(i,j,k)**2 + & + 2.d0 * (gabl * velxminus(i,j,k) * velyminus(i,j,k) + & + gacl * velxminus(i,j,k) * velzminus(i,j,k) + & + gbcl * velyminus(i,j,k) * velzminus(i,j,k) ) ) ) end do end do |