1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
|
#include "cctk.h"
c ========================================================================
SUBROUTINE met_rad_der(origin,Dx,Dy,Dz,x,y,z,g,dg)
c ------------------------------------------------------------------------
c
c Calculate isotropic radial derivative of unphysical metric,
c component or conformal factor at each point on the 3D grid,
c using derivatives in the Cartesian directions.
c
c ------------------------------------------------------------------------
IMPLICIT NONE
c Input variables
CCTK_REAL,INTENT(IN) ::
& Dx,Dy,Dz,origin(3)
CCTK_REAL,DIMENSION(:),INTENT(IN) ::
& x,y,z
CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) ::
& g
c Output variables
CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) ::
& dg
c Local variables, here only
INTEGER ::
& i,j,k
CCTK_REAL,PARAMETER ::
& zero = 0.0D0,
& half = 0.5D0
CCTK_REAL ::
& rad
c ------------------------------------------------------------------------
DO k = 2, SIZE(z)-1
DO j = 2, SIZE(y)-1
DO i = 2, SIZE(x)-1
rad = SQRT((x(i)-origin(1))**2
& +(y(j)-origin(2))**2
& +(z(k)-origin(3))**2)
IF (rad.NE.0) THEN
dg(i,j,k) = half/rad*
& ((x(i)-origin(1))/Dx*(g(i+1,j,k)-g(i-1,j,k))
& +(y(j)-origin(2))/Dy*(g(i,j+1,k)-g(i,j-1,k))
& +(z(k)-origin(3))/Dz*(g(i,j,k+1)-g(i,j,k-1)))
ELSE
dg(i,j,k) = zero
ENDIF
ENDDO
ENDDO
ENDDO
c This is needed when the grid is an octant, but it does not hurt
c if it is not
DO k = 2, size(z)-1
DO j = 2, size(y)-1
dg(1,j,k) = dg(2,j,k)
ENDDO
ENDDO
DO k = 2, size(z)-1
DO i = 1, size(x)-1
dg(i,1,k) = dg(i,2,k)
ENDDO
ENDDO
DO j = 1, size(y)-1
DO i = 1, size(x)-1
dg(i,j,1) = dg(i,j,2)
ENDDO
ENDDO
END SUBROUTINE met_rad_der
|
