/*@@
   @file      num_dt.F
   @date      October 1999
   @author    Miguel Alcubierre
   @desc
     Calculate time derivatives of position, momentum and proper time.
   @enddesc
   @version   $Id$
@@*/

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"


/*@@
  @routine    num_dt
  @date       August 1997
  @author     Miguel Alcubierre
  @desc
              This subroutine calculates the time derivatives
              of position, momentum and proper time.

              On input, {XS,YS,ZS} and {PXS,PYS,PZS} are
              the position and momentum respectively.

              On output, {XS,YS,ZS} and {PXS,PYS,PZS} are
              the time derivatives of the position and momentum
              respectively, and S is the time derivative of the
              proper time.

              I need many auxilliary arrays because I have to
              interpolate over the whole grid at the same time.
              The obvious way to eliminate all these arrays is
              to interpolate on a point by point basis.  However,
              this is not a good idea because each call to the
              interpolation routine has A LARGE overhead since the
              routine is MPI aware and must set up communications
              between all processors.  This means that calling it
              on a point by point basis causes the code to run very
              slow on multiprocessors (about 10 times slower on 2
              processors than on just 1; believe me, I tried it).

  @enddesc
  @calls      getpoints
  @calledby   timegeodesic_num_integrate
@@*/

      subroutine num_dt(CCTK_ARGUMENTS,XS,YS,ZS,PXS,PYS,PZS,S,lgeos)

      implicit none

      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_FUNCTIONS

      integer nx,ny,nz
      integer i,l,m,p,q
      integer npoints
      CCTK_INT lgeos
      integer ierror,alp_index
      integer gxx_index,gyy_index,gzz_index
      integer gxy_index,gxz_index,gyz_index
      integer SHIFT_ACTIVE

      integer vindex,interp_handle,coord_system_handle,param_table_handle
      character(30) options_string
      CCTK_POINTER, dimension(3) :: interp_coords
      CCTK_INT, dimension(7) :: in_array_indices
      CCTK_POINTER, dimension(7) :: out_arrays
      CCTK_INT, dimension(7) :: out_array_type_codes

      CCTK_REAL dx,dy,dz,idx,idy,idz,hdx,hdy,hdz
      CCTK_REAL zero,half,one,two
      CCTK_REAL aux

      CCTK_REAL, dimension (3,3)   :: GU,DL,DU
      CCTK_REAL, dimension (lgeos) :: A,S
      CCTK_REAL, dimension (lgeos) :: G1,G2,G3,G4,G5,G6
      CCTK_REAL, dimension (lgeos) :: T1,T2,T3,T4,T5,T6
      CCTK_REAL, dimension (lgeos) :: U1,U2,U3,U4,U5,U6
      CCTK_REAL, dimension (lgeos) :: V1,V2,V3,V4,V5,V6
      CCTK_REAL, dimension (lgeos) :: XS,YS,ZS,PXS,PYS,PZS

!     Description of variables:
!
!     i,l,m,p,q   Counters.
!
!     XS          Array that holds x position of the geodesics
!                 on input to the subroutine, and its time derivative
!                 on output.
!
!     YS          Array that holds y position of the geodesics
!                 on input to the subroutine, and its time derivative
!                 on output.
!
!     ZS          Array that holds z position of the geodesics
!                 on input to the subroutine, and its time derivative
!                 on output.
!
!     PXS         Array that holds x momentum of the geodesics
!                 on input to the subroutine, and its time derivative
!                 on output.
!
!     PYS         Array that holds y momentum of the geodesics
!                 on input to the subroutine, and its time derivative
!                 on output.
!
!     PZS         Array that holds z momentum of the geodesics
!                 on input to the subroutine, and its time derivative
!                 on output.
!
!     S           Array that contains the time derivatives of the
!                 proper time.
!
!     A           Array with interpolated values of the lapse.
!
!     G#          Arrays that contain the interpolated inverse metric.
!
!     GU,DL,DU    Small auxilliary arrays.
!
!     T#,U#,V#    Large auxilliary arrays.


!     **************************
!     ***   DEFINE NUMBERS   ***
!     **************************

      zero = 0.0d0
      half = 0.5d0
      one  = 1.0d0
      two  = 2.0d0

!     Grid parameters.

      nx = cctk_lsh(1)
      ny = cctk_lsh(2)
      nz = cctk_lsh(3)

      dx = cctk_delta_space(1)
      dy = cctk_delta_space(2)
      dz = cctk_delta_space(3)

      hdx = half*dx
      hdy = half*dy
      hdz = half*dz

      idx = one/dx
      idy = one/dy
      idz = one/dz

c     This was added in order to fit into the squeme of Einstein2
      SHIFT_ACTIVE = 1

!     **************************************
!     ***   INTERPOLATION PRELIMINARIES  ***
!     **************************************

      npoints = lgeos

!     parameter, local interpolator, and coordinate system handle.
      param_table_handle = -1
      interp_handle = -1
      coord_system_handle = -1

      options_string = "order = " // char(ichar('0') + interpolation_order)
      call Util_TableCreateFromString (param_table_handle, options_string)
      if (param_table_handle .lt. 0) then
        call CCTK_WARN(0,"Cannot create parameter table for interpolator")
      endif

      call CCTK_InterpHandle (interp_handle, "uniform cartesian")
      if (interp_handle .lt. 0) then
        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
      endif

      call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
      if (coord_system_handle .lt. 0) then
        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
      endif

!     fill in the input/output arrays for the interpolator
      interp_coords(1) = CCTK_PointerTo(XS)
      interp_coords(2) = CCTK_PointerTo(YS)
      interp_coords(3) = CCTK_PointerTo(ZS)

!     ************************************************
!     ***   INTERPOLATE LAPSE AND SPATIAL METRIC   ***
!     ************************************************

!     Interpolate the lapse and store it in the array A.
!     Interpolate the components of the 3-metric and store
!     them in the arrays {G1,G2,G3,G4,G5,G6}.
      call CCTK_VarIndex (vindex, "admbase::alp")
      in_array_indices(1) = vindex
      call CCTK_VarIndex (vindex, "admbase::gxx")
      in_array_indices(2) = vindex
      call CCTK_VarIndex (vindex, "admbase::gyy")
      in_array_indices(3) = vindex
      call CCTK_VarIndex (vindex, "admbase::gzz")
      in_array_indices(4) = vindex
      call CCTK_VarIndex (vindex, "admbase::gxy")
      in_array_indices(5) = vindex
      call CCTK_VarIndex (vindex, "admbase::gxz")
      in_array_indices(6) = vindex
      call CCTK_VarIndex (vindex, "admbase::gyz")
      in_array_indices(7) = vindex

      out_arrays(1) = CCTK_PointerTo(A)
      out_arrays(2) = CCTK_PointerTo(G1)
      out_arrays(3) = CCTK_PointerTo(G2)
      out_arrays(4) = CCTK_PointerTo(G3)
      out_arrays(5) = CCTK_PointerTo(G4)
      out_arrays(6) = CCTK_PointerTo(G5)
      out_arrays(7) = CCTK_PointerTo(G6)

      out_array_type_codes = CCTK_VARIABLE_REAL

      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            7, in_array_indices,
     .                            7, out_array_type_codes, out_arrays)
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif


!     *******************************
!     ***   FIND INVERSE METRIC   ***
!     *******************************

!     Find the inverse of the determinant of the 3-metric
!     and store it in the array T1.

      T1 = one/(G1*G2*G3 + two*G4*G5*G6
     &   - (G1*G6**2 + G2*G5**2 + G3*G4**2))

!     Find inverse metric components and store them in the
!     arrays {U1,U2,U3,U4,U5,U6}.

      U1 = T1*(G2*G3 - G6**2)
      U2 = T1*(G1*G3 - G5**2)
      U3 = T1*(G1*G2 - G4**2)

      U4 = T1*(G5*G6 - G3*G4)
      U5 = T1*(G4*G6 - G2*G5)
      U6 = T1*(G4*G5 - G1*G6)

!     Copy the inverse metric into the arrays {G1,G2,G3,G4,G5,G6}.

      G1 = U1
      G2 = U2
      G3 = U3

      G4 = U4
      G5 = U5
      G6 = U6


!     ******************************************
!     ***   FIND DERIVATIVE OF PROPER TIME   ***
!     ******************************************

!     Find the time derivative of the proper time
!     and store it in the array S.

      S = A/dsqrt(one + G1*PXS**2 + G2*PYS**2 + G3*PZS**2
     &  + two*(G4*PXS*PYS + G5*PXS*PZS + G6*PYS*PZS))


!     *********************************************
!     ***   FIND TIME DERIVATIVES OF POSITION   ***
!     *********************************************

!     Find shift independent contributions to time derivatives
!     of the positions and store them in {T1,T2,T3}.

      T1 = S*(G1*PXS + G4*PYS + G5*PZS)
      T2 = S*(G4*PXS + G2*PYS + G6*PZS)
      T3 = S*(G5*PXS + G6*PYS + G3*PZS)

!     Add contributions from the shift (only if the shift is active).

      if (shift_state == SHIFT_ACTIVE) then

         call CCTK_VarIndex (vindex, "admbase::betax")
         in_array_indices(1) = vindex
         call CCTK_VarIndex (vindex, "admbase::betay")
         in_array_indices(2) = vindex
         call CCTK_VarIndex (vindex, "admbase::betaz")
         in_array_indices(3) = vindex

         out_arrays(1) = CCTK_PointerTo(U1)
         out_arrays(2) = CCTK_PointerTo(U2)
         out_arrays(3) = CCTK_PointerTo(U3)

         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               3, in_array_indices,
     .                               3, out_array_type_codes, out_arrays)
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif


         T1 = T1 - U1
         T2 = T2 - U2
         T3 = T3 - U3

      end if


!     *********************************************
!     ***   FIND TIME DERIVATIVES OF MOMENTUM   ***
!     *********************************************

!     Find contributions from derivatives of the lapse to
!     the time derivatives of the momentum and store them
!     in arrays {T4,T5,T6}.  For this I need to interpolate
!     the lapse in displaced points and then take the difference.

      call CCTK_VarIndex (vindex, "admbase::alp")
      in_array_indices(1) = vindex

      out_arrays(1) = CCTK_PointerTo(U1)
      out_arrays(2) = CCTK_PointerTo(V1)
      out_arrays(3) = CCTK_PointerTo(U2)
      out_arrays(4) = CCTK_PointerTo(V2)
      out_arrays(5) = CCTK_PointerTo(U3)
      out_arrays(6) = CCTK_PointerTo(V3)

      XS = XS + hdx
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            1, in_array_indices,
     .                            1, out_array_type_codes, out_arrays(1))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      XS = XS - hdx

      XS = XS - hdx
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            1, in_array_indices,
     .                            1, out_array_type_codes, out_arrays(2))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      XS = XS + hdx
      T4 = - A/S*(U1-V1)*idx

      YS = YS + hdy
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            1, in_array_indices,
     .                            1, out_array_type_codes, out_arrays(3))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      YS = YS - hdy

      YS = YS - hdy
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            1, in_array_indices,
     .                            1, out_array_type_codes, out_arrays(4))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      YS = YS + hdy
      T5 = - A/S*(U2-V2)*idy

      ZS = ZS + hdz
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            1, in_array_indices,
     .                            1, out_array_type_codes, out_arrays(5))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      ZS = ZS - hdz

      ZS = ZS - hdz
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            1, in_array_indices,
     .                            1, out_array_type_codes, out_arrays(6))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      ZS = ZS + hdz
      T6 = - A/S*(U3-V3)*idz

!     Add contributions from derivatives of the shift
!     (only if the shift is active).

      if (shift_state == SHIFT_ACTIVE) then

         call CCTK_VarIndex (vindex, "admbase::betax")
         in_array_indices(1) = vindex
         call CCTK_VarIndex (vindex, "admbase::betay")
         in_array_indices(2) = vindex
         call CCTK_VarIndex (vindex, "admbase::betaz")
         in_array_indices(3) = vindex

         out_arrays(1) = CCTK_PointerTo(U1)
         out_arrays(2) = CCTK_PointerTo(U2)
         out_arrays(3) = CCTK_PointerTo(U3)
         out_arrays(4) = CCTK_PointerTo(V1)
         out_arrays(5) = CCTK_PointerTo(V2)
         out_arrays(6) = CCTK_PointerTo(V3)

         XS = XS + hdx
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               3, in_array_indices,
     .                               3, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif
         XS = XS - hdx

         XS = XS - hdx
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               3, in_array_indices,
     .                               3, out_array_type_codes, out_arrays(4))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif
         XS = XS + hdx

         U1 = (U1-V1)*idx
         U2 = (U2-V2)*idx
         U3 = (U3-V3)*idx

         T4 = T4 + (PXS*U1 + PYS*U2 + PZS*U3)

         YS = YS + hdy
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               3, in_array_indices,
     .                               3, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif
         YS = YS - hdy

         YS = YS - hdy
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               3, in_array_indices,
     .                               3, out_array_type_codes, out_arrays(4))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif
         YS = YS + hdy

         U1 = (U1-V1)*idy
         U2 = (U2-V2)*idy
         U3 = (U3-V3)*idy

         T5 = T5 + (PXS*U1 + PYS*U2 + PZS*U3)

         ZS = ZS + hdz
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               3, in_array_indices,
     .                               3, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif
         ZS = ZS - hdz

         ZS = ZS - hdz
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                               param_table_handle, coord_system_handle,
     .                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                               3, in_array_indices,
     .                               3, out_array_type_codes, out_arrays(4))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif
         ZS = ZS + hdz

         U1 = (U1-V1)*idz
         U2 = (U2-V2)*idz
         U3 = (U3-V3)*idz

         T6 = T6 + (PXS*U1 + PYS*U2 + PZS*U3)

      end if

!     Add contributions from derivatives of the metric.
!     This is somewhat complicated and done in stages to
!     avoid introducing even more large arrays.

!     Interpolate x derivatives of the metric and store
!     them in {U1,U2,U3,U4,U5,U6}.

      call CCTK_VarIndex (vindex, "admbase::gxx")
      in_array_indices(1) = vindex
      call CCTK_VarIndex (vindex, "admbase::gyy")
      in_array_indices(2) = vindex
      call CCTK_VarIndex (vindex, "admbase::gzz")
      in_array_indices(3) = vindex
      call CCTK_VarIndex (vindex, "admbase::gxy")
      in_array_indices(4) = vindex
      call CCTK_VarIndex (vindex, "admbase::gxz")
      in_array_indices(5) = vindex
      call CCTK_VarIndex (vindex, "admbase::gyz")
      in_array_indices(6) = vindex

      out_arrays(1) = CCTK_PointerTo(U1)
      out_arrays(2) = CCTK_PointerTo(U2)
      out_arrays(3) = CCTK_PointerTo(U3)
      out_arrays(4) = CCTK_PointerTo(U4)
      out_arrays(5) = CCTK_PointerTo(U5)
      out_arrays(6) = CCTK_PointerTo(U6)

      XS = XS + hdx
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            6, in_array_indices,
     .                            6, out_array_type_codes, out_arrays)
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      XS = XS - hdx

      out_arrays(1) = CCTK_PointerTo(V1)
      out_arrays(2) = CCTK_PointerTo(V2)
      out_arrays(3) = CCTK_PointerTo(V3)
      out_arrays(4) = CCTK_PointerTo(V4)
      out_arrays(5) = CCTK_PointerTo(V5)
      out_arrays(6) = CCTK_PointerTo(V6)

      XS = XS - hdx
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            6, in_array_indices,
     .                            6, out_array_type_codes, out_arrays)
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      XS = XS + hdx

      U1 = (U1-V1)*idx
      U2 = (U2-V2)*idx
      U3 = (U3-V3)*idx
      U4 = (U4-V4)*idx
      U5 = (U5-V5)*idx
      U6 = (U6-V6)*idx

!     Calculate the x derivatives of the inverse metric
!     and store them again in {U1,U2,U3,U4,U5,U6}.

      do i=1,lgeos

         GU(1,1) = G1(i)
         GU(1,2) = G4(i)
         GU(1,3) = G5(i)
         GU(2,1) = G4(i)
         GU(2,2) = G2(i)
         GU(2,3) = G6(i)
         GU(3,1) = G5(i)
         GU(3,2) = G6(i)
         GU(3,3) = G3(i)

         DL(1,1) = U1(i)
         DL(1,2) = U4(i)
         DL(1,3) = U5(i)
         DL(2,1) = U4(i)
         DL(2,2) = U2(i)
         DL(2,3) = U6(i)
         DL(3,1) = U5(i)
         DL(3,2) = U6(i)
         DL(3,3) = U3(i)

         do m=1,3
            do l=1,3

               aux = zero

               do q=1,3
                  do p=1,3
                     aux = aux - GU(l,p)*GU(m,q)*DL(p,q)
                  end do
               end do

               DU(l,m) = aux

            end do
         end do

         U1(i) = DU(1,1)
         U2(i) = DU(2,2)
         U3(i) = DU(3,3)
         U4(i) = DU(1,2)
         U5(i) = DU(1,3)
         U6(i) = DU(2,3)

      end do

!     Add contributions from derivatives of the inverse metric
!     to the time derivative of the x component of the momentum.

      T4 = T4 - half*S*(PXS**2*U1 + PYS**2*U2 + PZS**2*U3
     &   + two*(PXS*PYS*U4 + PXS*PZS*U5 + PYS*PZS*U6))

!     Interpolate y derivatives of the metric and store
!     them in {U1,U2,U3,U4,U5,U6}.

      out_arrays(1) = CCTK_PointerTo(U1)
      out_arrays(2) = CCTK_PointerTo(U2)
      out_arrays(3) = CCTK_PointerTo(U3)
      out_arrays(4) = CCTK_PointerTo(U4)
      out_arrays(5) = CCTK_PointerTo(U5)
      out_arrays(6) = CCTK_PointerTo(U6)

      YS = YS + hdy
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            6, in_array_indices,
     .                            6, out_array_type_codes, out_arrays)
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      YS = YS - hdy

      out_arrays(1) = CCTK_PointerTo(V1)
      out_arrays(2) = CCTK_PointerTo(V2)
      out_arrays(3) = CCTK_PointerTo(V3)
      out_arrays(4) = CCTK_PointerTo(V4)
      out_arrays(5) = CCTK_PointerTo(V5)
      out_arrays(6) = CCTK_PointerTo(V6)

      YS = YS - hdy
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            6, in_array_indices,
     .                            6, out_array_type_codes, out_arrays)
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      YS = YS + hdy

      U1 = (U1-V1)*idy
      U2 = (U2-V2)*idy
      U3 = (U3-V3)*idy
      U4 = (U4-V4)*idy
      U5 = (U5-V5)*idy
      U6 = (U6-V6)*idy

!     Calculate the y derivatives of the inverse metric
!     and store them again in {U1,U2,U3,U4,U5,U6}.

      do i=1,lgeos

         GU(1,1) = G1(i)
         GU(1,2) = G4(i)
         GU(1,3) = G5(i)
         GU(2,1) = G4(i)
         GU(2,2) = G2(i)
         GU(2,3) = G6(i)
         GU(3,1) = G5(i)
         GU(3,2) = G6(i)
         GU(3,3) = G3(i)

         DL(1,1) = U1(i)
         DL(1,2) = U4(i)
         DL(1,3) = U5(i)
         DL(2,1) = U4(i)
         DL(2,2) = U2(i)
         DL(2,3) = U6(i)
         DL(3,1) = U5(i)
         DL(3,2) = U6(i)
         DL(3,3) = U3(i)

         do m=1,3
            do l=1,3

               aux = zero

               do q=1,3
                  do p=1,3
                     aux = aux - GU(l,p)*GU(m,q)*DL(p,q)
                  end do
               end do

               DU(l,m) = aux

            end do
         end do

         U1(i) = DU(1,1)
         U2(i) = DU(2,2)
         U3(i) = DU(3,3)
         U4(i) = DU(1,2)
         U5(i) = DU(1,3)
         U6(i) = DU(2,3)

      end do

!     Add contributions from derivatives of the inverse metric
!     to the time derivative of the y component of the momentum.

      T5 = T5 - half*S*(PXS**2*U1 + PYS**2*U2 + PZS**2*U3
     &   + two*(PXS*PYS*U4 + PXS*PZS*U5 + PYS*PZS*U6))

!     Interpolate z derivatives of the metric and store
!     them in {U1,U2,U3,U4,U5,U6}.

      out_arrays(1) = CCTK_PointerTo(U1)
      out_arrays(2) = CCTK_PointerTo(U2)
      out_arrays(3) = CCTK_PointerTo(U3)
      out_arrays(4) = CCTK_PointerTo(U4)
      out_arrays(5) = CCTK_PointerTo(U5)
      out_arrays(6) = CCTK_PointerTo(U6)

      ZS = ZS + hdz
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            6, in_array_indices,
     .                            6, out_array_type_codes, out_arrays)
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      ZS = ZS - hdz

      out_arrays(1) = CCTK_PointerTo(V1)
      out_arrays(2) = CCTK_PointerTo(V2)
      out_arrays(3) = CCTK_PointerTo(V3)
      out_arrays(4) = CCTK_PointerTo(V4)
      out_arrays(5) = CCTK_PointerTo(V5)
      out_arrays(6) = CCTK_PointerTo(V6)

      ZS = ZS - hdz
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     .                            param_table_handle, coord_system_handle,
     .                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     .                            6, in_array_indices,
     .                            6, out_array_type_codes, out_arrays)
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif
      ZS = ZS + hdz

      U1 = (U1-V1)*idz
      U2 = (U2-V2)*idz
      U3 = (U3-V3)*idz
      U4 = (U4-V4)*idz
      U5 = (U5-V5)*idz
      U6 = (U6-V6)*idz

!     Calculate the z derivatives of the inverse metric
!     and store them again in {U1,U2,U3,U4,U5,U6}.

      do i=1,lgeos

         GU(1,1) = G1(i)
         GU(1,2) = G4(i)
         GU(1,3) = G5(i)
         GU(2,1) = G4(i)
         GU(2,2) = G2(i)
         GU(2,3) = G6(i)
         GU(3,1) = G5(i)
         GU(3,2) = G6(i)
         GU(3,3) = G3(i)

         DL(1,1) = U1(i)
         DL(1,2) = U4(i)
         DL(1,3) = U5(i)
         DL(2,1) = U4(i)
         DL(2,2) = U2(i)
         DL(2,3) = U6(i)
         DL(3,1) = U5(i)
         DL(3,2) = U6(i)
         DL(3,3) = U3(i)

         do m=1,3
            do l=1,3

               aux = zero

               do q=1,3
                  do p=1,3
                     aux = aux - GU(l,p)*GU(m,q)*DL(p,q)
                  end do
               end do

               DU(l,m) = aux

            end do
         end do

         U1(i) = DU(1,1)
         U2(i) = DU(2,2)
         U3(i) = DU(3,3)
         U4(i) = DU(1,2)
         U5(i) = DU(1,3)
         U6(i) = DU(2,3)

      end do

!     Add contributions from derivatives of the inverse metric
!     to the time derivative of the z component of the momentum.

      T6 = T6 - half*S*(PXS**2*U1 + PYS**2*U2 + PZS**2*U3
     &   + two*(PXS*PYS*U4 + PXS*PZS*U5 + PYS*PZS*U6))


!     ****************************************************
!     ***   COPY TIME DERIVATIVES INTO OUTPUT ARRAYS   ***
!     ****************************************************

!     Copy time derivatives of position into arrays {XS,YS,ZS}.

      XS = T1
      YS = T2
      ZS = T3

!     Copy time derivatives of momentum into arrays {PXS,PYS,PZS}.

      PXS = T4
      PYS = T5
      PZS = T6


!     ***************
!     ***   END   ***
!     ***************

      return
      end