Fix compiler warnings.

Most of them could be fixed by renaming .F77 files to .F
Some had to be fixed by explicitly declaring some variables using
CCTK_DECLARE() (which also only works for .F, not for .F77)



git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Exact/trunk@287 e296648e-0e4f-0410-bd07-d597d9acff87

1 files changed, 157 insertions, 0 deletions

diff --git a/src/metrics/Kerr_KerrSchild.F b/src/metrics/Kerr_KerrSchild.F
new file mode 100644
index 0000000..36ea76f
--- /dev/null
+++ b/src/metrics/Kerr_KerrSchild.F
@@ -0,0 +1,157 @@
+C     Kerr-Schild form of boosted rotating black hole.
+C     Program g_ab = eta_ab + H l_a l_b, g^ab = eta^ab - H l^a l^b.
+C     Here eta_ab is Minkowski in Cartesian coordinates, H is a scalar,
+C     and l is a null vector.
+C
+C Author: unknown
+C Copyright/License: unknown
+C
+C $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Functions.h"
+
+      subroutine Exact__Kerr_KerrSchild(
+     $     x, y, z, t,
+     $     gdtt, gdtx, gdty, gdtz, 
+     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
+     $     gutt, gutx, guty, gutz, 
+     $     guxx, guyy, guzz, guxy, guyz, guzx,
+     $     psi, Tmunu_flag)
+
+      implicit none
+
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+c input arguments
+      CCTK_REAL x, y, z, t
+
+c output arguments
+      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
+     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
+     $       gutt, gutx, guty, gutz, 
+     $       guxx, guyy, guzz, guxy, guyz, guzx
+      CCTK_DECLARE(CCTK_REAL, psi,)
+      LOGICAL   Tmunu_flag
+
+c local variables
+      CCTK_REAL boostv, eps, m, a
+      integer   power
+
+c local variables
+      CCTK_REAL gamma, t0, z0, x0, y0, rho02, r02, r0, costheta0, 
+     $     lt0, lx0, ly0, lz0, hh, lt, lx, ly, lz
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+C This is a vacuum spacetime with no cosmological constant
+      Tmunu_flag = .false.
+
+C     Get parameters of the exact solution,
+C     and convert from parameter file spin parameter J/m^2
+C     to the J/m definition used in the code here.
+
+      boostv = Kerr_KerrSchild__boost_v
+      eps    = Kerr_KerrSchild__epsilon
+      power  = Kerr_KerrSchild__power
+      m      = Kerr_KerrSchild__mass
+      a      = m*Kerr_KerrSchild__spin
+
+C     Boost factor.
+
+      gamma = 1.d0 / sqrt(1.d0 - boostv**2)
+
+C     Lorentz transform t,x,y,z -> t0,x0,y0,z0. 
+C     t0 is never used, but is here for illustration, and we introduce
+C     x0 and y0 also only for clarity.
+C     Note that z0 = 0 means z = vt for the BH.
+
+      t0 = gamma * ((t - Kerr_KerrSchild__t) - boostv * (z - Kerr_KerrSchild__z))
+      z0 = gamma * ((z - Kerr_KerrSchild__z) - boostv * (t - Kerr_KerrSchild__t))
+      x0 = x - Kerr_KerrSchild__x
+      y0 = y - Kerr_KerrSchild__y
+
+C     Coordinate distance to center of black hole. Note it moves!
+
+      rho02 = x0**2 + y0**2 + z0**2
+
+C     Spherical auxiliary coordinate r and angle theta in BH rest frame.
+
+      r02 = 0.5d0 * (rho02 - a**2) 
+     $     + sqrt(0.25d0 * (rho02 - a**2)**2 + a**2 * z0**2)
+      r0 = sqrt(max(0.0d0,r02))
+      if (Kerr_KerrSchild__parabolic .eq. 0) then
+C        Use a power law to avoid the singularity
+         r0 = (r0**power + eps**power)**(1.0d0/power)
+      else
+         if (r0 .lt. eps) then
+            if (power .eq. 0) then
+               r0 = eps
+            else if (power .eq. 2) then
+               r0 = eps/2 + r0**2 * 1/(2*eps)
+            else if (power .eq. 4) then
+               r0 = 3*eps/8 + r0**2 * (3/(4*eps) - r0**2 * 1/(8*eps**3))
+            else if (power .eq. 6) then
+               r0 = 5*eps/16 + r0**2 * (15/(16*eps) + r0**2 * (-5/(16*eps**3) + r0**2 * 1/(16*eps**5)))
+            else if (power .eq. 8) then
+               r0 = 35*eps/128 + r0**2 * (35/(32*eps) + r0**2 * (-35/(64*eps**3) + r0**2 * (7/(32*eps**5) - r0**2 * 5/(128*eps**7))))
+            else
+               call CCTK_WARN (CCTK_WARN_ABORT, "Unsupported value of parameter Kerr_KerrSchild__power")
+            end if
+         end if
+      end if
+C     Another idea:
+C        r0 = r0 + eps * exp(-x/eps)
+
+      costheta0 = z0 / r0
+      
+C     Coefficient H. Note this transforms as a scalar, so does not carry
+C     the suffix 0.
+      hh = m * r0 / (r0**2 + a**2 * costheta0**2)
+      
+C     Components of l_a in rest frame. Note indices down.
+      lt0 = 1.d0
+      lx0 = (r0 * x0 + a * y0) / (r0**2 + a**2)
+      ly0 = (r0 * y0 - a * x0) / (r0**2 + a**2)
+      lz0 = z0 / r0
+
+C     Now boost it to coordinates x, y, z, t.
+C     This is the reverse Lorentz transformation, but applied
+C     to a one-form, so the sign of boostv is the same as the forward
+C     Lorentz transformation applied to the coordinates.
+
+      lt = gamma * (lt0 - boostv * lz0) 
+      lz = gamma * (lz0 - boostv * lt0) 
+      lx = lx0
+      ly = ly0
+      
+C     Down metric. g_ab = flat_ab + H l_a l_b
+
+      gdtt = - 1.d0 + 2.d0 * hh * lt * lt
+      gdtx = 2.d0 * hh * lt * lx
+      gdty = 2.d0 * hh * lt * ly
+      gdtz = 2.d0 * hh * lt * lz
+      gdxx = 1.d0 + 2.d0 * hh * lx * lx
+      gdyy = 1.d0 + 2.d0 * hh * ly * ly
+      gdzz = 1.d0 + 2.d0 * hh * lz * lz
+      gdxy = 2.d0 * hh * lx * ly
+      gdyz = 2.d0 * hh * ly * lz
+      gdzx = 2.d0 * hh * lz * lx
+
+C     Up metric. g^ab = flat^ab - H l^a l^b.
+C     Notice that g^ab = g_ab and l^i = l_i and l^0 = - l_0 in flat spacetime.
+      gutt = - 1.d0 - 2.d0 * hh * lt * lt 
+      gutx = 2.d0 * hh * lt * lx
+      guty = 2.d0 * hh * lt * ly
+      gutz = 2.d0 * hh * lt * lz
+      guxx = 1.d0 - 2.d0 * hh * lx * lx
+      guyy = 1.d0 - 2.d0 * hh * ly * ly
+      guzz = 1.d0 - 2.d0 * hh * lz * lz
+      guxy = - 2.d0 * hh * lx * ly
+      guyz = - 2.d0 * hh * ly * lz
+      guzx = - 2.d0 * hh * lz * lx
+
+      return
+      end