diff options
Diffstat (limited to 'src/CalcK.F90')
| -rw-r--r-- | src/CalcK.F90 | 15 |
1 files changed, 11 insertions, 4 deletions
diff --git a/src/CalcK.F90 b/src/CalcK.F90 index f77d500..2dacc50 100644 --- a/src/CalcK.F90 +++ b/src/CalcK.F90 @@ -22,6 +22,7 @@ subroutine CalcK (CCTK_ARGUMENTS) CCTK_REAL :: dt, dx(3) + CCTK_REAL :: ps, dps(3) CCTK_REAL :: gama(3,3), gama_dot(3,3), dgama(3,3,3) CCTK_REAL :: alfa CCTK_REAL :: beta(3), dbeta(3,3) @@ -53,6 +54,8 @@ subroutine CalcK (CCTK_ARGUMENTS) do j = imin(2), imax(2) do i = imin(1), imax(1) + ps = psi(i,j,k) + gama(1,1) = gxx(i,j,k) gama(1,2) = gxy(i,j,k) gama(1,3) = gxz(i,j,k) @@ -82,6 +85,10 @@ subroutine CalcK (CCTK_ARGUMENTS) gama_dot(3,1) = gama_dot(1,3) gama_dot(3,2) = gama_dot(2,3) + dps(1) = (psi(i+1,j,k) - psi(i-1,j,k)) / (2*dx(1)) + dps(2) = (psi(i,j+1,k) - psi(i,j-1,k)) / (2*dx(2)) + dps(3) = (psi(i,j,k+1) - psi(i,j,k-1)) / (2*dx(3)) + dgama(1,1,1) = (gxx(i+1,j,k) - gxx(i-1,j,k)) / (2*dx(1)) dgama(1,2,1) = (gxy(i+1,j,k) - gxy(i-1,j,k)) / (2*dx(1)) dgama(1,3,1) = (gxz(i+1,j,k) - gxz(i-1,j,k)) / (2*dx(1)) @@ -127,11 +134,11 @@ subroutine CalcK (CCTK_ARGUMENTS) do a=1,3 do b=1,3 - kk(a,b) = - gama_dot(a,b) + kk(a,b) = - ps**4 * gama_dot(a,b) do c=1,3 - kk(a,b) = kk(a,b) + gama(k,j) * dbeta(k,i) & - & + gama(i,k) * dbeta(k,j) & - & + beta(k) * dgama(i,j,k) + kk(a,b) = kk(a,b) + ps**4 * gama(k,j) * dbeta(k,i) & + & + ps**4 * gama(i,k) * dbeta(k,j) & + & + beta(k) * (ps**4 * dgama(i,j,k) + 4*ps**3 * dps(k) * gama(i,j)) end do kk(a,b) = kk(a,b) / (2*alfa) end do |