1 files changed, 91 insertions, 0 deletions
diff --git a/src/met_rad_der.F b/src/met_rad_der.F
new file mode 100644
index 0000000..31821f0
--- /dev/null
+++ b/src/met_rad_der.F
@@ -0,0 +1,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
+
+
+
+
+
+